Changeset e17904a in sasmodels


Ignore:
Timestamp:
Mar 18, 2016 9:39:00 AM (9 years ago)
Author:
wojciech
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
a629d8e
Parents:
fad5dc1 (diff), 08376e7 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' of https://github.com/SasView/sasmodels

Files:
4 edited

Legend:

Unmodified
Added
Removed
  • explore/J1.py

    ra5af4e1 raceac39  
    22Show numerical precision of $2 J_1(x)/x$. 
    33""" 
     4import sys; sys.path.insert(0, '..') 
    45 
    56import numpy as np 
    6 #from sympy.mpmath import mp 
    7 from mpmath import mp 
     7try: 
     8    from mpmath import mp 
     9except: 
     10    # CRUFT: mpmath split out into its own package 
     11    from sympy.mpmath import mp 
    812#import matplotlib; matplotlib.use('TkAgg') 
    913import pylab 
  • sasmodels/sasview_model.py

    r2622b3f r08376e7  
    2323from . import core 
    2424 
    25 def make_class(model_info, dtype='single', namestyle='name'): 
     25def standard_models(): 
     26    return [make_class(model_name) for model_name in core.list_models()] 
     27 
     28def make_class(model_name, namestyle='name'): 
    2629    """ 
    2730    Load the sasview model defined in *kernel_module*. 
    2831 
    29     Returns a class that can be used directly as a sasview model. 
     32    Returns a class that can be used directly as a sasview model.t 
    3033 
    3134    Defaults to using the new name for a model.  Setting 
     
    3336    compatible with SasView. 
    3437    """ 
    35     model = core.build_model(model_info, dtype=dtype) 
     38    model_info = core.load_model_info(model_name) 
    3639    def __init__(self, multfactor=1): 
    37         SasviewModel.__init__(self, model) 
    38     attrs = dict(__init__=__init__) 
    39     ConstructedModel = type(model.info[namestyle], (SasviewModel,), attrs) 
     40        SasviewModel.__init__(self) 
     41    attrs = dict(__init__=__init__, _model_info=model_info) 
     42    ConstructedModel = type(model_info[namestyle], (SasviewModel,), attrs) 
    4043    return ConstructedModel 
    4144 
     
    4447    Sasview wrapper for opencl/ctypes model. 
    4548    """ 
    46     def __init__(self, model): 
    47         """Initialization""" 
    48         self._model = model 
    49  
    50         self.name = model.info['name'] 
    51         self.oldname = model.info['oldname'] 
    52         self.description = model.info['description'] 
     49    def __init__(self): 
     50        self._kernel = None 
     51        model_info = self._model_info 
     52 
     53        self.name = model_info['name'] 
     54        self.oldname = model_info['oldname'] 
     55        self.description = model_info['description'] 
    5356        self.category = None 
    5457        self.multiplicity_info = None 
     
    6063        self.params = collections.OrderedDict() 
    6164        self.dispersion = dict() 
    62         partype = model.info['partype'] 
    63  
    64         for p in model.info['parameters']: 
     65        partype = model_info['partype'] 
     66 
     67        for p in model_info['parameters']: 
    6568            self.params[p.name] = p.default 
    6669            self.details[p.name] = [p.units] + p.limits 
     
    8386 
    8487        ## independent parameter name and unit [string] 
    85         self.input_name = model.info.get("input_name", "Q") 
    86         self.input_unit = model.info.get("input_unit", "A^{-1}") 
    87         self.output_name = model.info.get("output_name", "Intensity") 
    88         self.output_unit = model.info.get("output_unit", "cm^{-1}") 
     88        self.input_name = model_info.get("input_name", "Q") 
     89        self.input_unit = model_info.get("input_unit", "A^{-1}") 
     90        self.output_name = model_info.get("output_name", "Intensity") 
     91        self.output_unit = model_info.get("output_unit", "cm^{-1}") 
    8992 
    9093        ## _persistency_dict is used by sas.perspectives.fitting.basepage 
     
    9598        ## New fields introduced for opencl rewrite 
    9699        self.cutoff = 1e-5 
     100 
     101    def __get_state__(self): 
     102        state = self.__dict__.copy() 
     103        model_id = self._model_info['id'] 
     104        state.pop('_kernel') 
     105        # May need to reload model info on set state since it has pointers 
     106        # to python implementations of Iq, etc. 
     107        #state.pop('_model_info') 
     108        return state 
     109 
     110    def __set_state__(self, state): 
     111        self.__dict__ = state 
     112        self._kernel = None 
    97113 
    98114    def __str__(self): 
     
    187203        # TODO: fix test so that parameter order doesn't matter 
    188204        ret = ['%s.%s' % (d.lower(), p) 
    189                for d in self._model.info['partype']['pd-2d'] 
     205               for d in self._model_info['partype']['pd-2d'] 
    190206               for p in ('npts', 'nsigmas', 'width')] 
    191207        #print(ret) 
     
    261277            # Check whether we have a list of ndarrays [qx,qy] 
    262278            qx, qy = qdist 
    263             partype = self._model.info['partype'] 
     279            partype = self._model_info['partype'] 
    264280            if not partype['orientation'] and not partype['magnetic']: 
    265281                return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) 
     
    284300        to the card for each evaluation. 
    285301        """ 
     302        if self._kernel is None: 
     303            self._kernel = core.build_model(self._model_info) 
    286304        q_vectors = [np.asarray(q) for q in args] 
    287         fn = self._model(q_vectors) 
     305        fn = self._kernel(q_vectors) 
    288306        pars = [self.params[v] for v in fn.fixed_pars] 
    289307        pd_pars = [self._get_weights(p) for p in fn.pd_pars] 
     
    299317        :return: the value of the effective radius 
    300318        """ 
    301         ER = self._model.info.get('ER', None) 
     319        ER = self._model_info.get('ER', None) 
    302320        if ER is None: 
    303321            return 1.0 
     
    314332        :return: the value of the volf ratio 
    315333        """ 
    316         VR = self._model.info.get('VR', None) 
     334        VR = self._model_info.get('VR', None) 
    317335        if VR is None: 
    318336            return 1.0 
     
    358376        parameter set in the vector. 
    359377        """ 
    360         pars = self._model.info['partype']['volume'] 
     378        pars = self._model_info['partype']['volume'] 
    361379        return core.dispersion_mesh([self._get_weights(p) for p in pars]) 
    362380 
     
    368386        from . import weights 
    369387 
    370         relative = self._model.info['partype']['pd-rel'] 
    371         limits = self._model.info['limits'] 
     388        relative = self._model_info['partype']['pd-rel'] 
     389        limits = self._model_info['limits'] 
    372390        dis = self.dispersion[par] 
    373391        value, weight = weights.get_weights( 
     
    375393            self.params[par], limits[par], par in relative) 
    376394        return value, weight / np.sum(weight) 
     395 
  • sasmodels/models/lib/j1_cephes.c

    ra5af4e1 rfad5dc1  
    3939Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier 
    4040*/ 
    41 double J1(double ); 
     41 
     42constant double RP[8] = { 
     43    -8.99971225705559398224E8, 
     44    4.52228297998194034323E11, 
     45    -7.27494245221818276015E13, 
     46    3.68295732863852883286E15, 
     47    0.0, 
     48    0.0, 
     49    0.0, 
     50    0.0 }; 
     51 
     52constant double RQ[8] = { 
     53        6.20836478118054335476E2, 
     54        2.56987256757748830383E5, 
     55        8.35146791431949253037E7, 
     56        2.21511595479792499675E10, 
     57        4.74914122079991414898E12, 
     58        7.84369607876235854894E14, 
     59        8.95222336184627338078E16, 
     60        5.32278620332680085395E18 
     61    }; 
     62 
     63constant double PP[8] = { 
     64    7.62125616208173112003E-4, 
     65    7.31397056940917570436E-2, 
     66    1.12719608129684925192E0, 
     67    5.11207951146807644818E0, 
     68    8.42404590141772420927E0, 
     69    5.21451598682361504063E0, 
     70    1.00000000000000000254E0, 
     71    0.0} ; 
     72 
     73 
     74constant double PQ[8] = { 
     75    5.71323128072548699714E-4, 
     76    6.88455908754495404082E-2, 
     77    1.10514232634061696926E0, 
     78    5.07386386128601488557E0, 
     79    8.39985554327604159757E0, 
     80    5.20982848682361821619E0, 
     81    9.99999999999999997461E-1, 
     82    0.0 }; 
     83 
     84constant double QP[8] = { 
     85    5.10862594750176621635E-2, 
     86    4.98213872951233449420E0, 
     87    7.58238284132545283818E1, 
     88    3.66779609360150777800E2, 
     89    7.10856304998926107277E2, 
     90    5.97489612400613639965E2, 
     91    2.11688757100572135698E2, 
     92    2.52070205858023719784E1 }; 
     93 
     94constant double QQ[8] = { 
     95    7.42373277035675149943E1, 
     96    1.05644886038262816351E3, 
     97    4.98641058337653607651E3, 
     98    9.56231892404756170795E3, 
     99    7.99704160447350683650E3, 
     100    2.82619278517639096600E3, 
     101    3.36093607810698293419E2, 
     102    0.0 }; 
     103 
     104constant double JP[8] = { 
     105        -4.878788132172128E-009, 
     106        6.009061827883699E-007, 
     107        -4.541343896997497E-005, 
     108        1.937383947804541E-003, 
     109        -3.405537384615824E-002, 
     110        0.0, 
     111        0.0, 
     112        0.0 
     113    }; 
     114 
     115constant double MO1[8] = { 
     116        6.913942741265801E-002, 
     117        -2.284801500053359E-001, 
     118        3.138238455499697E-001, 
     119        -2.102302420403875E-001, 
     120        5.435364690523026E-003, 
     121        1.493389585089498E-001, 
     122        4.976029650847191E-006, 
     123        7.978845453073848E-001 
     124    }; 
     125 
     126constant double PH1[8] = { 
     127        -4.497014141919556E+001, 
     128        5.073465654089319E+001, 
     129        -2.485774108720340E+001, 
     130        7.222973196770240E+000, 
     131        -1.544842782180211E+000, 
     132        3.503787691653334E-001, 
     133        -1.637986776941202E-001, 
     134        3.749989509080821E-001 
     135    }; 
    42136 
    43137double J1(double x) { 
     
    60154        { 
    61155            z = x * x; 
    62             w = polevlRP( z, 3 ) / p1evlRQ( z, 8 ); 
     156            w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 ); 
    63157            w = w * x * (z - Z1) * (z - Z2); 
    64158            return( w ); 
     
    68162    z = w * w; 
    69163 
    70     p = polevlPP( z, 6)/polevlPQ( z, 6 ); 
    71     q = polevlQP( z, 7)/p1evlQQ( z, 7 ); 
     164    p = polevl( z, PP, 6)/polevl( z, PQ, 6 ); 
     165    q = polevl( z, QP, 7)/p1evl( z, QQ, 7 ); 
    72166 
    73167    xn = x - THPIO4; 
     
    95189        { 
    96190            z = xx * xx; 
    97             p = (z-Z1) * xx * polevlJP( z, 4 ); 
     191            p = (z-Z1) * xx * polevl( z, JP, 4 ); 
    98192            return( p ); 
    99193        } 
     
    102196    w = sqrt(q); 
    103197 
    104     p = w * polevlMO1( q, 7); 
     198    p = w * polevl( q, MO1, 7); 
    105199    w = q*q; 
    106     xn = q * polevlPH1( w, 7) - THPIO4F; 
     200    xn = q * polevl( w, PH1, 7) - THPIO4F; 
    107201    p = p * cos(xn + xx); 
    108202 
     
    110204#endif 
    111205} 
     206 
     207 
     208/*double J1c(double x) { 
     209    return 
     210}*/ 
  • sasmodels/models/lib/polevl.c

    r3b12dea rfad5dc1  
    5151*/ 
    5252 
    53 double polevlRP(double x, int N ) { 
    5453 
    55     double coef[8] = { 
    56     -8.99971225705559398224E8, 
    57     4.52228297998194034323E11, 
    58     -7.27494245221818276015E13, 
    59     3.68295732863852883286E15, 
    60     0.0, 
    61     0.0, 
    62     0.0, 
    63     0.0 }; 
    64  
    65     int i = 0; 
    66     double ans = coef[i]; 
    67  
    68     while (i < N) { 
    69         i++; 
    70         ans = ans * x + coef[i]; 
    71     } 
    72     return ans ; 
    73 } 
    74  
    75 double polevlRQ(double x, int N ) { 
    76  
    77     double coef[8] = { 
    78         6.20836478118054335476E2, 
    79         2.56987256757748830383E5, 
    80         8.35146791431949253037E7, 
    81         2.21511595479792499675E10, 
    82         4.74914122079991414898E12, 
    83         7.84369607876235854894E14, 
    84         8.95222336184627338078E16, 
    85         5.32278620332680085395E18 
    86     }; 
    87  
    88     int i = 0; 
    89     double ans = coef[i]; 
    90  
    91     while (i < N) { 
    92         i++; 
    93         ans = ans * x + coef[i]; 
    94     } 
    95     return ans ; 
    96 } 
    97  
    98 double polevlPP(double x, int N ) { 
    99  
    100     double coef[8] = { 
    101     7.62125616208173112003E-4, 
    102     7.31397056940917570436E-2, 
    103     1.12719608129684925192E0, 
    104     5.11207951146807644818E0, 
    105     8.42404590141772420927E0, 
    106     5.21451598682361504063E0, 
    107     1.00000000000000000254E0, 
    108     0.0} ; 
    109  
    110     int i = 0; 
    111     double ans = coef[i]; 
    112  
    113     while (i < N) { 
    114         i++; 
    115         ans = ans * x + coef[i]; 
    116     } 
    117     return ans ; 
    118 } 
    119  
    120 double polevlPQ(double x, int N ) { 
    121  
    122     double coef[8] = { 
    123     5.71323128072548699714E-4, 
    124     6.88455908754495404082E-2, 
    125     1.10514232634061696926E0, 
    126     5.07386386128601488557E0, 
    127     8.39985554327604159757E0, 
    128     5.20982848682361821619E0, 
    129     9.99999999999999997461E-1, 
    130     0.0 }; 
    131  
    132     int i = 0; 
    133     double ans = coef[i]; 
    134  
    135     while (i < N) { 
    136         i++; 
    137         ans = ans * x + coef[i]; 
    138     } 
    139     return ans ; 
    140 } 
    141  
    142 double polevlQP(double x, int N ) { 
    143  
    144     double coef[8] = { 
    145     5.10862594750176621635E-2, 
    146     4.98213872951233449420E0, 
    147     7.58238284132545283818E1, 
    148     3.66779609360150777800E2, 
    149     7.10856304998926107277E2, 
    150     5.97489612400613639965E2, 
    151     2.11688757100572135698E2, 
    152     2.52070205858023719784E1 }; 
    153  
    154     int i = 0; 
    155     double ans = coef[i]; 
    156  
    157     while (i < N) { 
    158         i++; 
    159         ans = ans * x + coef[i]; 
    160     } 
    161     return ans ; 
    162  
    163 } 
    164  
    165 double polevlQQ(double x, int N ) { 
    166  
    167     double coef[8] = { 
    168     7.42373277035675149943E1, 
    169     1.05644886038262816351E3, 
    170     4.98641058337653607651E3, 
    171     9.56231892404756170795E3, 
    172     7.99704160447350683650E3, 
    173     2.82619278517639096600E3, 
    174     3.36093607810698293419E2, 
    175     0.0 }; 
    176  
    177     int i = 0; 
    178     double ans = coef[i]; 
    179  
    180     while (i < N) { 
    181         i++; 
    182         ans = ans * x + coef[i]; 
    183     } 
    184     return ans ; 
    185  
    186 } 
    187  
    188 double polevlJP( double x, int N ) { 
    189     double coef[8] = { 
    190         -4.878788132172128E-009, 
    191         6.009061827883699E-007, 
    192         -4.541343896997497E-005, 
    193         1.937383947804541E-003, 
    194         -3.405537384615824E-002, 
    195         0.0, 
    196         0.0, 
    197         0.0 
    198     }; 
    199  
    200     int i = 0; 
    201     double ans = coef[i]; 
    202  
    203     while (i < N) { 
    204         i++; 
    205         ans = ans * x + coef[i]; 
    206     } 
    207     return ans ; 
    208  
    209 } 
    210  
    211 double polevlMO1( double x, int N ) { 
    212     double coef[8] = { 
    213         6.913942741265801E-002, 
    214         -2.284801500053359E-001, 
    215         3.138238455499697E-001, 
    216         -2.102302420403875E-001, 
    217         5.435364690523026E-003, 
    218         1.493389585089498E-001, 
    219         4.976029650847191E-006, 
    220         7.978845453073848E-001 
    221     }; 
    222  
    223     int i = 0; 
    224     double ans = coef[i]; 
    225  
    226     while (i < N) { 
    227         i++; 
    228         ans = ans * x + coef[i]; 
    229     } 
    230     return ans ; 
    231 } 
    232  
    233 double polevlPH1( double x, int N ) { 
    234  
    235     double coef[8] = { 
    236         -4.497014141919556E+001, 
    237         5.073465654089319E+001, 
    238         -2.485774108720340E+001, 
    239         7.222973196770240E+000, 
    240         -1.544842782180211E+000, 
    241         3.503787691653334E-001, 
    242         -1.637986776941202E-001, 
    243         3.749989509080821E-001 
    244     }; 
    245  
    246     int i = 0; 
    247     double ans = coef[i]; 
    248  
    249     while (i < N) { 
    250         i++; 
    251         ans = ans * x + coef[i]; 
    252     } 
    253     return ans ; 
    254 } 
    255  
    256 /*double polevl( double x, double coef[8], int N ) { 
     54double polevl( double x, constant double *coef, int N ) { 
    25755 
    25856    int i = 0; 
     
    26664    return ans ; 
    26765 
    268 }*/ 
     66} 
    26967 
    27068/*                                                      p1evl() */ 
     
    27472 */ 
    27573 
    276 double p1evlRP( double x, int N ) { 
    277     double coef[8] = { 
    278     -8.99971225705559398224E8, 
    279     4.52228297998194034323E11, 
    280     -7.27494245221818276015E13, 
    281     3.68295732863852883286E15, 
    282     0.0, 
    283     0.0, 
    284     0.0, 
    285     0.0 }; 
    286  
    287     int i=0; 
    288     double ans = x+coef[i]; 
    289  
    290     while (i < N-1) { 
    291         i++; 
    292         ans = ans*x + coef[i]; 
    293     } 
    294  
    295     return( ans ); 
    296  
    297 } 
    298  
    299 double p1evlRQ( double x, int N ) { 
    300 //1: RQ 
    301     double coef[8] = { 
    302     6.20836478118054335476E2, 
    303     2.56987256757748830383E5, 
    304     8.35146791431949253037E7, 
    305     2.21511595479792499675E10, 
    306     4.74914122079991414898E12, 
    307     7.84369607876235854894E14, 
    308     8.95222336184627338078E16, 
    309     5.32278620332680085395E18}; 
    310  
     74double p1evl( double x, constant double *coef, int N ) { 
    31175    int i=0; 
    31276    double ans = x+coef[i]; 
     
    31983    return( ans ); 
    32084} 
    321  
    322 double p1evlPP( double x, int N ) { 
    323 //3 : PP 
    324     double coef[8] = { 
    325     7.62125616208173112003E-4, 
    326     7.31397056940917570436E-2, 
    327     1.12719608129684925192E0, 
    328     5.11207951146807644818E0, 
    329     8.42404590141772420927E0, 
    330     5.21451598682361504063E0, 
    331     1.00000000000000000254E0, 
    332     0.0}; 
    333  
    334     int i=0; 
    335     double ans = x+coef[i]; 
    336  
    337     while (i < N-1) { 
    338         i++; 
    339         ans = ans*x + coef[i]; 
    340     } 
    341  
    342     return( ans ); 
    343 } 
    344  
    345 double p1evlPQ( double x, int N ) { 
    346 //4: PQ 
    347     double coef[8] = { 
    348     5.71323128072548699714E-4, 
    349     6.88455908754495404082E-2, 
    350     1.10514232634061696926E0, 
    351     5.07386386128601488557E0, 
    352     8.39985554327604159757E0, 
    353     5.20982848682361821619E0, 
    354     9.99999999999999997461E-1, 
    355     0.0}; 
    356  
    357     int i=0; 
    358     double ans = x+coef[i]; 
    359  
    360     while (i < N-1) { 
    361         i++; 
    362         ans = ans*x + coef[i]; 
    363     } 
    364  
    365     return( ans ); 
    366 } 
    367  
    368 double p1evlQP( double x, int N ) { 
    369 //5: QP 
    370     double coef[8] = { 
    371     5.10862594750176621635E-2, 
    372     4.98213872951233449420E0, 
    373     7.58238284132545283818E1, 
    374     3.66779609360150777800E2, 
    375     7.10856304998926107277E2, 
    376     5.97489612400613639965E2, 
    377     2.11688757100572135698E2, 
    378     2.52070205858023719784E1 }; 
    379  
    380     int i=0; 
    381     double ans = x+coef[i]; 
    382  
    383     while (i < N-1) { 
    384         i++; 
    385         ans = ans*x + coef[i]; 
    386     } 
    387  
    388     return( ans ); 
    389 } 
    390  
    391 double p1evlQQ( double x, int N ) { 
    392     double coef[8] = { 
    393     7.42373277035675149943E1, 
    394     1.05644886038262816351E3, 
    395     4.98641058337653607651E3, 
    396     9.56231892404756170795E3, 
    397     7.99704160447350683650E3, 
    398     2.82619278517639096600E3, 
    399     3.36093607810698293419E2, 
    400     0.0}; 
    401  
    402     int i=0; 
    403     double ans = x+coef[i]; 
    404  
    405     while (i < N-1) { 
    406         i++; 
    407         ans = ans*x + coef[i]; 
    408     } 
    409  
    410     return( ans ); 
    411  
    412 } 
    413  
    414 double p1evlJP( double x, int N ) { 
    415     double coef[8] = { 
    416         -4.878788132172128E-009, 
    417         6.009061827883699E-007, 
    418         -4.541343896997497E-005, 
    419         1.937383947804541E-003, 
    420         -3.405537384615824E-002, 
    421         0.0, 
    422         0.0, 
    423         0.0}; 
    424  
    425     int i=0; 
    426     double ans = x+coef[i]; 
    427  
    428     while (i < N-1) { 
    429         i++; 
    430         ans = ans*x + coef[i]; 
    431     } 
    432  
    433     return( ans ); 
    434 } 
    435  
    436 double p1evlMO1( double x, int N ) { 
    437     double coef[8] = { 
    438         6.913942741265801E-002, 
    439         -2.284801500053359E-001, 
    440         3.138238455499697E-001, 
    441         -2.102302420403875E-001, 
    442         5.435364690523026E-003, 
    443         1.493389585089498E-001, 
    444         4.976029650847191E-006, 
    445         7.978845453073848E-001}; 
    446  
    447     int i=0; 
    448     double ans = x+coef[i]; 
    449  
    450     while (i < N-1) { 
    451         i++; 
    452         ans = ans*x + coef[i]; 
    453     } 
    454  
    455     return( ans ); 
    456  
    457 } 
    458  
    459 double p1evlPH1( double x, int N ) { 
    460     double coef[8] = { 
    461         -4.497014141919556E+001, 
    462         5.073465654089319E+001, 
    463         -2.485774108720340E+001, 
    464         7.222973196770240E+000, 
    465         -1.544842782180211E+000, 
    466         3.503787691653334E-001, 
    467         -1.637986776941202E-001, 
    468         3.749989509080821E-001}; 
    469     int i=0; 
    470     double ans = x+coef[i]; 
    471  
    472     while (i < N-1) { 
    473         i++; 
    474         ans = ans*x + coef[i]; 
    475     } 
    476  
    477     return( ans ); 
    478 } 
    479  
    480 /*double p1evl( double x, double coef[8], int N ) { 
    481     int i=0; 
    482     double ans = x+coef[i]; 
    483  
    484     while (i < N-1) { 
    485         i++; 
    486         ans = ans*x + coef[i]; 
    487     } 
    488  
    489     return( ans ); 
    490  
    491 }*/ 
Note: See TracChangeset for help on using the changeset viewer.