Changeset e17904a in sasmodels
- Timestamp:
- Mar 18, 2016 11:39:00 AM (9 years ago)
- 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. - Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
explore/J1.py
ra5af4e1 raceac39 2 2 Show numerical precision of $2 J_1(x)/x$. 3 3 """ 4 import sys; sys.path.insert(0, '..') 4 5 5 6 import numpy as np 6 #from sympy.mpmath import mp 7 from mpmath import mp 7 try: 8 from mpmath import mp 9 except: 10 # CRUFT: mpmath split out into its own package 11 from sympy.mpmath import mp 8 12 #import matplotlib; matplotlib.use('TkAgg') 9 13 import pylab -
sasmodels/sasview_model.py
r2622b3f r08376e7 23 23 from . import core 24 24 25 def make_class(model_info, dtype='single', namestyle='name'): 25 def standard_models(): 26 return [make_class(model_name) for model_name in core.list_models()] 27 28 def make_class(model_name, namestyle='name'): 26 29 """ 27 30 Load the sasview model defined in *kernel_module*. 28 31 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 30 33 31 34 Defaults to using the new name for a model. Setting … … 33 36 compatible with SasView. 34 37 """ 35 model = core.build_model(model_info, dtype=dtype)38 model_info = core.load_model_info(model_name) 36 39 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) 40 43 return ConstructedModel 41 44 … … 44 47 Sasview wrapper for opencl/ctypes model. 45 48 """ 46 def __init__(self , model):47 """Initialization"""48 self._model = model49 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'] 53 56 self.category = None 54 57 self.multiplicity_info = None … … 60 63 self.params = collections.OrderedDict() 61 64 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']: 65 68 self.params[p.name] = p.default 66 69 self.details[p.name] = [p.units] + p.limits … … 83 86 84 87 ## 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}") 89 92 90 93 ## _persistency_dict is used by sas.perspectives.fitting.basepage … … 95 98 ## New fields introduced for opencl rewrite 96 99 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 97 113 98 114 def __str__(self): … … 187 203 # TODO: fix test so that parameter order doesn't matter 188 204 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'] 190 206 for p in ('npts', 'nsigmas', 'width')] 191 207 #print(ret) … … 261 277 # Check whether we have a list of ndarrays [qx,qy] 262 278 qx, qy = qdist 263 partype = self._model .info['partype']279 partype = self._model_info['partype'] 264 280 if not partype['orientation'] and not partype['magnetic']: 265 281 return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) … … 284 300 to the card for each evaluation. 285 301 """ 302 if self._kernel is None: 303 self._kernel = core.build_model(self._model_info) 286 304 q_vectors = [np.asarray(q) for q in args] 287 fn = self._ model(q_vectors)305 fn = self._kernel(q_vectors) 288 306 pars = [self.params[v] for v in fn.fixed_pars] 289 307 pd_pars = [self._get_weights(p) for p in fn.pd_pars] … … 299 317 :return: the value of the effective radius 300 318 """ 301 ER = self._model .info.get('ER', None)319 ER = self._model_info.get('ER', None) 302 320 if ER is None: 303 321 return 1.0 … … 314 332 :return: the value of the volf ratio 315 333 """ 316 VR = self._model .info.get('VR', None)334 VR = self._model_info.get('VR', None) 317 335 if VR is None: 318 336 return 1.0 … … 358 376 parameter set in the vector. 359 377 """ 360 pars = self._model .info['partype']['volume']378 pars = self._model_info['partype']['volume'] 361 379 return core.dispersion_mesh([self._get_weights(p) for p in pars]) 362 380 … … 368 386 from . import weights 369 387 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'] 372 390 dis = self.dispersion[par] 373 391 value, weight = weights.get_weights( … … 375 393 self.params[par], limits[par], par in relative) 376 394 return value, weight / np.sum(weight) 395 -
sasmodels/models/lib/j1_cephes.c
ra5af4e1 rfad5dc1 39 39 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier 40 40 */ 41 double J1(double ); 41 42 constant 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 52 constant 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 63 constant 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 74 constant 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 84 constant 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 94 constant 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 104 constant 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 115 constant 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 126 constant 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 }; 42 136 43 137 double J1(double x) { … … 60 154 { 61 155 z = x * x; 62 w = polevl RP( z, 3 ) / p1evlRQ( z, 8 );156 w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 ); 63 157 w = w * x * (z - Z1) * (z - Z2); 64 158 return( w ); … … 68 162 z = w * w; 69 163 70 p = polevl PP( z, 6)/polevlPQ( z, 6 );71 q = polevl QP( 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 ); 72 166 73 167 xn = x - THPIO4; … … 95 189 { 96 190 z = xx * xx; 97 p = (z-Z1) * xx * polevl JP( z, 4 );191 p = (z-Z1) * xx * polevl( z, JP, 4 ); 98 192 return( p ); 99 193 } … … 102 196 w = sqrt(q); 103 197 104 p = w * polevl MO1( q, 7);198 p = w * polevl( q, MO1, 7); 105 199 w = q*q; 106 xn = q * polevl PH1( w, 7) - THPIO4F;200 xn = q * polevl( w, PH1, 7) - THPIO4F; 107 201 p = p * cos(xn + xx); 108 202 … … 110 204 #endif 111 205 } 206 207 208 /*double J1c(double x) { 209 return 210 }*/ -
sasmodels/models/lib/polevl.c
r3b12dea rfad5dc1 51 51 */ 52 52 53 double polevlRP(double x, int N ) {54 53 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 ) { 54 double polevl( double x, constant double *coef, int N ) { 257 55 258 56 int i = 0; … … 266 64 return ans ; 267 65 268 } */66 } 269 67 270 68 /* p1evl() */ … … 274 72 */ 275 73 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 74 double p1evl( double x, constant double *coef, int N ) { 311 75 int i=0; 312 76 double ans = x+coef[i]; … … 319 83 return( ans ); 320 84 } 321 322 double p1evlPP( double x, int N ) {323 //3 : PP324 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: PQ347 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: QP370 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.