Changeset 33e91b1 in sasmodels
- Timestamp:
- Mar 3, 2015 4:28:50 PM (10 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:
- 53d0e24, 3a45c2c
- Parents:
- 6c8db9e
- Location:
- sasmodels
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/generate.py
ra5d0d00 r33e91b1 201 201 # Scale and background, which are parameters common to every form factor 202 202 COMMON_PARAMETERS = [ 203 [ "scale", "", 1, [0, np.inf], "", "Source intensity"],204 [ "background", "1/cm", 0, [0, np.inf], "", "Source background"],203 ["scale", "", 1, [0, np.inf], "", "Source intensity"], 204 ["background", "1/cm", 0, [0, np.inf], "", "Source background"], 205 205 ] 206 206 … … 233 233 # description and its parameter table. The remainder of the doc comes 234 234 # from the module docstring. 235 DOC_HEADER =""".. _%(name)s:235 DOC_HEADER = """.. _%(name)s: 236 236 237 237 %(label)s … … 259 259 ] 260 260 column_widths = [max(w, len(h)) 261 for w, h in zip(column_widths, PARTABLE_HEADERS)]261 for w, h in zip(column_widths, PARTABLE_HEADERS)] 262 262 263 263 sep = " ".join("="*w for w in column_widths) 264 264 lines = [ 265 265 sep, 266 " ".join("%-*s" %(w,h) for w,h in zip(column_widths, PARTABLE_HEADERS)),266 " ".join("%-*s" % (w, h) for w, h in zip(column_widths, PARTABLE_HEADERS)), 267 267 sep, 268 268 ] 269 269 for p in pars: 270 270 lines.append(" ".join([ 271 "%-*s" %(column_widths[0],p[0]),272 "%-*s" %(column_widths[1],p[-1]),273 "%-*s" %(column_widths[2],RST_UNITS[p[1]]),274 "%*g" %(column_widths[3],p[2]),271 "%-*s" % (column_widths[0], p[0]), 272 "%-*s" % (column_widths[1], p[-1]), 273 "%-*s" % (column_widths[2], RST_UNITS[p[1]]), 274 "%*g" % (column_widths[3], p[2]), 275 275 ])) 276 276 lines.append(sep) … … 287 287 if exists(target): 288 288 return target 289 raise ValueError("%r not found in %s" %(filename, search_path))289 raise ValueError("%r not found in %s" % (filename, search_path)) 290 290 291 291 def sources(info): … … 293 293 Return a list of the sources file paths for the module. 294 294 """ 295 search_path = [ 296 abspath(joinpath(dirname(__file__),'models'))]295 search_path = [dirname(info['filename']), 296 abspath(joinpath(dirname(__file__), 'models'))] 297 297 return [_search(search_path, f) for f in info['source']] 298 298 … … 333 333 334 334 for p in pars: 335 name, ptype = p[0],p[4]335 name, ptype = p[0], p[4] 336 336 if ptype == 'volume': 337 337 partype['pd-1d'].append(name) … … 346 346 partype['fixed-2d'].append(name) 347 347 else: 348 raise ValueError("unknown parameter type %r" %ptype)348 raise ValueError("unknown parameter type %r" % ptype) 349 349 partype[ptype].append(name) 350 350 … … 356 356 """ 357 357 spaces = " "*depth 358 sep = "\n" +spaces358 sep = "\n" + spaces 359 359 return spaces + sep.join(s.split("\n")) 360 360 … … 366 366 Returns loop opening and loop closing 367 367 """ 368 LOOP_OPEN ="""\368 LOOP_OPEN = """\ 369 369 for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 370 370 const double %(name)s = loops[2*(%(name)s_i%(offset)s)]; … … 376 376 loop_end = [] 377 377 for name in pd_pars: 378 subst = { 'name': name, 'offset': offset}379 loop_head.append(indent(LOOP_OPEN %subst, depth))378 subst = {'name': name, 'offset': offset} 379 loop_head.append(indent(LOOP_OPEN % subst, depth)) 380 380 loop_end.insert(0, (" "*depth) + "}") 381 offset += '+N' +name381 offset += '+N' + name 382 382 depth += 2 383 383 return "\n".join(loop_head), "\n".join(loop_end) 384 384 385 C_KERNEL_TEMPLATE =None385 C_KERNEL_TEMPLATE = None 386 386 def make_model(info): 387 387 """ … … 416 416 417 417 iq_parameters = [p[0] 418 for p in info['parameters'][2:] # skip scale, background419 if p[0] in set(fixed_1d+pd_1d)]418 for p in info['parameters'][2:] # skip scale, background 419 if p[0] in set(fixed_1d + pd_1d)] 420 420 iqxy_parameters = [p[0] 421 for p in info['parameters'][2:] # skip scale, background422 if p[0] in set(fixed_2d+pd_2d)]421 for p in info['parameters'][2:] # skip scale, background 422 if p[0] in set(fixed_2d + pd_2d)] 423 423 volume_parameters = [p[0] 424 for p in info['parameters']425 if p[4]=='volume']424 for p in info['parameters'] 425 if p[4] == 'volume'] 426 426 427 427 # Fill in defintions for volume parameters … … 430 430 ','.join(volume_parameters))) 431 431 defines.append(('VOLUME_WEIGHT_PRODUCT', 432 '*'.join(p +'_w' for p in volume_parameters)))432 '*'.join(p + '_w' for p in volume_parameters))) 433 433 434 434 # Generate form_volume function from body only 435 435 if info['form_volume'] is not None: 436 436 if volume_parameters: 437 vol_par_decl = ', '.join('double ' +p for p in volume_parameters)437 vol_par_decl = ', '.join('double ' + p for p in volume_parameters) 438 438 else: 439 439 vol_par_decl = 'void' … … 445 445 %(body)s 446 446 } 447 """ %{'body':info['form_volume']}447 """ % {'body':info['form_volume']} 448 448 source.append(fn) 449 449 450 450 # Fill in definitions for Iq parameters 451 defines.append(('IQ_KERNEL_NAME', info['name'] +'_Iq'))451 defines.append(('IQ_KERNEL_NAME', info['name'] + '_Iq')) 452 452 defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters))) 453 453 if fixed_1d: 454 454 defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS', 455 ', \\\n '.join('const double %s' %p for p in fixed_1d)))455 ', \\\n '.join('const double %s' % p for p in fixed_1d))) 456 456 if pd_1d: 457 457 defines.append(('IQ_WEIGHT_PRODUCT', 458 '*'.join(p +'_w' for p in pd_1d)))458 '*'.join(p + '_w' for p in pd_1d))) 459 459 defines.append(('IQ_DISPERSION_LENGTH_DECLARATIONS', 460 ', \\\n '.join('const int N%s' %p for p in pd_1d)))460 ', \\\n '.join('const int N%s' % p for p in pd_1d))) 461 461 defines.append(('IQ_DISPERSION_LENGTH_SUM', 462 '+'.join('N' +p for p in pd_1d)))462 '+'.join('N' + p for p in pd_1d))) 463 463 open_loops, close_loops = build_polydispersity_loops(pd_1d) 464 464 defines.append(('IQ_OPEN_LOOPS', 465 open_loops.replace('\n', ' \\\n')))465 open_loops.replace('\n', ' \\\n'))) 466 466 defines.append(('IQ_CLOSE_LOOPS', 467 close_loops.replace('\n', ' \\\n')))467 close_loops.replace('\n', ' \\\n'))) 468 468 if info['Iq'] is not None: 469 469 defines.append(('IQ_PARAMETER_DECLARATIONS', 470 ', '.join('double '+p for p in iq_parameters)))470 ', '.join('double ' + p for p in iq_parameters))) 471 471 fn = """\ 472 472 double Iq(double q, IQ_PARAMETER_DECLARATIONS); … … 474 474 %(body)s 475 475 } 476 """ %{'body':info['Iq']}476 """ % {'body':info['Iq']} 477 477 source.append(fn) 478 478 479 479 # Fill in definitions for Iqxy parameters 480 defines.append(('IQXY_KERNEL_NAME', info['name'] +'_Iqxy'))480 defines.append(('IQXY_KERNEL_NAME', info['name'] + '_Iqxy')) 481 481 defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters))) 482 482 if fixed_2d: 483 483 defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS', 484 ', \\\n '.join('const double %s' %p for p in fixed_2d)))484 ', \\\n '.join('const double %s' % p for p in fixed_2d))) 485 485 if pd_2d: 486 486 defines.append(('IQXY_WEIGHT_PRODUCT', 487 '*'.join(p +'_w' for p in pd_2d)))487 '*'.join(p + '_w' for p in pd_2d))) 488 488 defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS', 489 ', \\\n '.join('const int N%s' %p for p in pd_2d)))489 ', \\\n '.join('const int N%s' % p for p in pd_2d))) 490 490 defines.append(('IQXY_DISPERSION_LENGTH_SUM', 491 '+'.join('N' +p for p in pd_2d)))491 '+'.join('N' + p for p in pd_2d))) 492 492 open_loops, close_loops = build_polydispersity_loops(pd_2d) 493 493 defines.append(('IQXY_OPEN_LOOPS', 494 open_loops.replace('\n', ' \\\n')))494 open_loops.replace('\n', ' \\\n'))) 495 495 defines.append(('IQXY_CLOSE_LOOPS', 496 close_loops.replace('\n', ' \\\n')))496 close_loops.replace('\n', ' \\\n'))) 497 497 if info['Iqxy'] is not None: 498 498 defines.append(('IQXY_PARAMETER_DECLARATIONS', 499 ', '.join('double '+p for p in iqxy_parameters)))499 ', '.join('double ' + p for p in iqxy_parameters))) 500 500 fn = """\ 501 501 double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS); … … 503 503 %(body)s 504 504 } 505 """ %{'body':info['Iqxy']}505 """ % {'body':info['Iqxy']} 506 506 source.append(fn) 507 507 … … 513 513 514 514 #for d in defines: print d 515 DEFINES ='\n'.join('#define %s %s'%(k,v) for k,v in defines)516 SOURCES ='\n\n'.join(source)517 return C_KERNEL_TEMPLATE %{515 DEFINES = '\n'.join('#define %s %s' % (k, v) for k, v in defines) 516 SOURCES = '\n\n'.join(source) 517 return C_KERNEL_TEMPLATE % { 518 518 'DEFINES':DEFINES, 519 519 'SOURCES':SOURCES, … … 531 531 #print kernelfile 532 532 info = dict( 533 filename =abspath(kernel_module.__file__),534 name =kernel_module.name,535 title =kernel_module.title,536 description =kernel_module.description,537 parameters =COMMON_PARAMETERS + kernel_module.parameters,538 source =getattr(kernel_module, 'source', []),539 oldname =kernel_module.oldname,540 oldpars =kernel_module.oldpars,533 filename=abspath(kernel_module.__file__), 534 name=kernel_module.name, 535 title=kernel_module.title, 536 description=kernel_module.description, 537 parameters=COMMON_PARAMETERS + kernel_module.parameters, 538 source=getattr(kernel_module, 'source', []), 539 oldname=kernel_module.oldname, 540 oldpars=kernel_module.oldpars, 541 541 ) 542 542 # Fill in attributes which default to None 543 info.update((k, getattr(kernel_module, k, None))543 info.update((k, getattr(kernel_module, k, None)) 544 544 for k in ('ER', 'VR', 'form_volume', 'Iq', 'Iqxy')) 545 545 # Fill in the derived attributes 546 info['limits'] = dict((p[0], p[3]) for p in info['parameters'])546 info['limits'] = dict((p[0], p[3]) for p in info['parameters']) 547 547 info['partype'] = categorize_parameters(info['parameters']) 548 info['defaults'] = dict((p[0], p[2]) for p in info['parameters'])548 info['defaults'] = dict((p[0], p[2]) for p in info['parameters']) 549 549 550 550 # Assume if one part of the kernel is python then all parts are. … … 556 556 Return the documentation for the model. 557 557 """ 558 subst = dict(name=kernel_module.name.replace('_', '-'),558 subst = dict(name=kernel_module.name.replace('_', '-'), 559 559 label=" ".join(kernel_module.name.split('_')).capitalize(), 560 560 title=kernel_module.title, 561 561 parameters=make_partable(kernel_module.parameters), 562 562 docs=kernel_module.__doc__) 563 return DOC_HEADER %subst563 return DOC_HEADER % subst 564 564 565 565 … … 568 568 import datetime 569 569 from .models import cylinder 570 toc = lambda: (datetime.datetime.now()-tic).total_seconds()571 570 tic = datetime.datetime.now() 572 source, info = make(cylinder) 573 print "time:",toc() 571 toc = lambda: (datetime.datetime.now() - tic).total_seconds() 572 make(cylinder) 573 print "time:", toc() 574 574 575 575 def main(): … … 579 579 name = sys.argv[1] 580 580 import sasmodels.models 581 __import__('sasmodels.models.' +name)581 __import__('sasmodels.models.' + name) 582 582 model = getattr(sasmodels.models, name) 583 source, info = make(model); print source584 #print doc(model)583 source, _ = make(model) 584 print source 585 585 586 586 if __name__ == "__main__": -
sasmodels/models/parallelepiped.py
ra5d0d00 r33e91b1 9 9 ---------- 10 10 11 This model provides the form factor, *P(q)*, for a rectangular parallelepiped 12 (below) where the form factor is normalized by the volume of the 13 parallelepiped. If you need to apply polydispersity, see also the 11 This model provides the form factor, *P(q)*, for a rectangular parallelepiped 12 (below) where the form factor is normalized by the volume of the 13 parallelepiped. If you need to apply polydispersity, see also the 14 14 RectangularPrismModel_. 15 15 … … 20 20 P(Q) = {\text{scale} \over V} F^2(Q) + \text{background} 21 21 22 where the volume *V* = *A B C* and the averaging < > is applied over all 22 where the volume *V* = *A B C* and the averaging < > is applied over all 23 23 orientations for 1D. 24 24 … … 27 27 *Figure. Parallelepiped with the corresponding Definition of sides. 28 28 29 The edge of the solid must satisfy the condition that** *A* < *B* < *C*. 30 Then, assuming *a* = *A* / *B* < 1, *b* = *B* / *B* = 1, and 29 The edge of the solid must satisfy the condition that** *A* < *B* < *C*. 30 Then, assuming *a* = *A* / *B* < 1, *b* = *B* / *B* = 1, and 31 31 *c* = *C* / *B* > 1, the form factor is 32 32 33 33 .. math:: 34 34 35 P(q) = \frac{\textstyle{scale}}{V}\int_0^1 \phi(\mu \sqrt{1-\sigma^2},a) 35 P(q) = \frac{\textstyle{scale}}{V}\int_0^1 \phi(\mu \sqrt{1-\sigma^2},a) 36 36 [S(\mu c \sigma/2)]^2 d\sigma 37 37 … … 40 40 .. math:: 41 41 42 \phi(\mu,a) = \int_0^1 \{S[\frac{\mu}{2}\cos(\frac{\pi}{2}u)] 42 \phi(\mu,a) = \int_0^1 \{S[\frac{\mu}{2}\cos(\frac{\pi}{2}u)] 43 43 S[\frac{\mu a}{2}\sin(\frac{\pi}{2}u)]\}^2 du 44 44 45 45 S(x) = \frac{\sin x}{x} 46 46 47 47 \mu = qB 48 48 … … 53 53 \Delta\rho = \rho_{\textstyle p} - \rho_{\textstyle solvent} 54 54 55 The scattering intensity per unit volume is returned in units of |cm^-1|; 55 The scattering intensity per unit volume is returned in units of |cm^-1|; 56 56 ie, *I(q)* = |phi| *P(q)*\ . 57 57 58 NB: The 2nd virial coefficient of the parallelpiped is calculated based on 59 the averaged effective radius (= sqrt(*short_a* \* *short_b* / |pi|)) and 58 NB: The 2nd virial coefficient of the parallelpiped is calculated based on 59 the averaged effective radius (= sqrt(*short_a* \* *short_b* / |pi|)) and 60 60 length(= *long_c*) values, and used as the effective radius for 61 61 *S(Q)* when *P(Q)* \* *S(Q)* is applied. 62 62 63 To provide easy access to the orientation of the parallelepiped, we define 63 To provide easy access to the orientation of the parallelepiped, we define 64 64 three angles |theta|, |phi| and |bigpsi|. The definition of |theta| and |phi| 65 is the same as for the cylinder model (see also figures below). 66 The angle |bigpsi| is the rotational angle around the *long_c* axis against 67 the *q* plane. For example, |bigpsi| = 0 when the *short_b* axis is parallel 65 is the same as for the cylinder model (see also figures below). 66 The angle |bigpsi| is the rotational angle around the *long_c* axis against 67 the *q* plane. For example, |bigpsi| = 0 when the *short_b* axis is parallel 68 68 to the *x*-axis of the detector. 69 69 … … 83 83 ---------- 84 84 85 Validation of the code was done by comparing the output of the 1D calculation 86 to the angular average of the output of a 2D calculation over all possible 87 angles. The Figure below shows the comparison where the solid dot refers to 88 averaged 2D while the line represents the result of the 1D calculation (for 89 the averaging, 76, 180, 76 points are taken for the angles of |theta|, |phi|, 85 Validation of the code was done by comparing the output of the 1D calculation 86 to the angular average of the output of a 2D calculation over all possible 87 angles. The Figure below shows the comparison where the solid dot refers to 88 averaged 2D while the line represents the result of the 1D calculation (for 89 the averaging, 76, 180, 76 points are taken for the angles of |theta|, |phi|, 90 90 and |psi| respectively). 91 91 … … 96 96 *Figure. Comparison between 1D and averaged 2D.* 97 97 98 This model reimplements the form factor calculations implemented in a c-library 98 This model reimplements the form factor calculations implemented in a c-library 99 99 provided by the NIST Center for Neutron Research (Kline, 2006). 100 100 101 101 """ 102 102 103 from numpy import pi, inf 103 from numpy import pi, inf, sqrt 104 104 105 105 name = "parallelepiped" … … 108 108 P(q)= scale/V*integral from 0 to 1 of ... 109 109 phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma 110 110 111 111 phi(mu,a) = integral from 0 to 1 of .. 112 112 (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du 113 113 S(x) = sin(x)/x 114 114 mu = q*B 115 115 """ 116 116 category = "shape:parallelpiped" 117 117 118 118 parameters = [ 119 # [ "name", "units", default, [lower, upper], "type",120 # "description" ],121 [ "sld", "1e-6/Ang^2", 4, [-inf,inf], "",122 "Parallelepiped scattering length density"],123 [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "",124 "Solvent scattering length density"],125 [ "a_side", "Ang",35, [0, inf], "volume",126 "Shorter side of the parallelepiped"],127 [ "b_side", "Ang",75, [0, inf], "volume",128 "Second side of the parallelepiped"],129 [ "c_side", "Ang",400, [0, inf], "volume",130 "Larger side of the parallelepiped"],131 [ 132 "In plane angle"],133 [ 134 "Out of plane angle"],135 [ 136 "Rotation angle around its own c axis against q plane"],119 # [ "name", "units", default, [lower, upper], "type", 120 # "description" ], 121 ["sld", "1e-6/Ang^2", 4, [-inf, inf], "", 122 "Parallelepiped scattering length density"], 123 ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "", 124 "Solvent scattering length density"], 125 ["a_side", "Ang", 35, [0, inf], "volume", 126 "Shorter side of the parallelepiped"], 127 ["b_side", "Ang", 75, [0, inf], "volume", 128 "Second side of the parallelepiped"], 129 ["c_side", "Ang", 400, [0, inf], "volume", 130 "Larger side of the parallelepiped"], 131 ["theta", "degrees", 60, [-inf, inf], "orientation", 132 "In plane angle"], 133 ["phi", "degrees", 60, [-inf, inf], "orientation", 134 "Out of plane angle"], 135 ["psi", "degrees", 60, [-inf, inf], "orientation", 136 "Rotation angle around its own c axis against q plane"], 137 137 ] 138 138 139 source = [ "lib/J1.c", "lib/gauss76.c", "parallelepiped.c"]139 source = ["lib/J1.c", "lib/gauss76.c", "parallelepiped.c"] 140 140 141 141 def ER(a_side, b_side, c_side): … … 148 148 b = 0.5 * c_side 149 149 t1 = a * a * b 150 t2 = 1.0 + (b /a)*(1.0+a/b/2.0)*(1.0+pi*a/b/2.0)150 t2 = 1.0 + (b / a) * (1.0 + a / b / 2.0) * (1.0 + pi * a / b / 2.0) 151 151 ddd = 3.0 * t1 * t2 152 152 153 return 0.5 * (ddd) **(1./3.)153 return 0.5 * (ddd) ** (1. / 3.) 154 154 155 155 # parameters for demo … … 169 169 # For testing against the old sasview models, include the converted parameter 170 170 # names and the target sasview model name. 171 oldname ='ParallelepipedModel'172 oldpars =dict(theta='parallel_theta', phi='parallel_phi', psi='parallel_psi',173 a_side='short_a', b_side='short_b', c_side='long_c',174 sld='sldPipe', solvent_sld='sldSolv')171 oldname = 'ParallelepipedModel' 172 oldpars = dict(theta='parallel_theta', phi='parallel_phi', psi='parallel_psi', 173 a_side='short_a', b_side='short_b', c_side='long_c', 174 sld='sldPipe', solvent_sld='sldSolv') 175 175
Note: See TracChangeset
for help on using the changeset viewer.