Changeset d2bb604 in sasmodels
- Timestamp:
- Apr 5, 2016 12:34:30 PM (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:
- 21b116f
- Parents:
- 1e2a1ba
- Location:
- sasmodels
- Files:
-
- 1 deleted
- 21 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/compare.py
rea05c87 rd2bb604 38 38 from . import core 39 39 from . import kerneldll 40 from . import product41 40 from .data import plot_theory, empty_data1D, empty_data2D 42 41 from .direct_model import DirectModel -
sasmodels/convert.json
rec45c4f rd2bb604 623 623 "SphereSLDModel", 624 624 { 625 "n": "n_shells", 625 626 "radius_core": "rad_core", 626 627 "sld_solvent": "sld_solv" -
sasmodels/core.py
r1e2a1ba rd2bb604 178 178 from the *pars* dictionary for parameter value and parameter dispersion. 179 179 """ 180 value = values.get(parameter.name, parameter.default)180 value = float(values.get(parameter.name, parameter.default)) 181 181 relative = parameter.relative_pd 182 182 limits = parameter.limits … … 200 200 """ 201 201 value, weight = zip(*pars) 202 weight = [w if w else [1.] for w in weight] 202 203 value = [v.flatten() for v in meshgrid(*value)] 203 204 weight = np.vstack([v.flatten() for v in meshgrid(*weight)]) … … 279 280 else: 280 281 vol_pars = [get_weights(parameter, values) 281 for parameter in model_info['parameters'] 282 for parameter in model_info['parameters'].call_parameters 282 283 if parameter.type == 'volume'] 283 284 value, weight = dispersion_mesh(vol_pars) 284 285 individual_radii = ER(*value) 285 #print(values[0].shape, weights.shape, fv.shape)286 286 return np.sum(weight*individual_radii) / np.sum(weight) 287 287 … … 297 297 else: 298 298 vol_pars = [get_weights(parameter, values) 299 for parameter in model_info['parameters'] 299 for parameter in model_info['parameters'].call_parameters 300 300 if parameter.type == 'volume'] 301 301 value, weight = dispersion_mesh(vol_pars) -
sasmodels/generate.py
r1e2a1ba rd2bb604 486 486 # Define the function calls 487 487 if partable.form_volume_parameters: 488 refs = _call_pars(" v.", partable.form_volume_parameters)489 call_volume = "#define CALL_VOLUME( v) form_volume(%s)" % (",".join(refs))488 refs = _call_pars("_v.", partable.form_volume_parameters) 489 call_volume = "#define CALL_VOLUME(_v) form_volume(%s)" % (",".join(refs)) 490 490 else: 491 491 # Model doesn't have volume. We could make the kernel run a little … … 495 495 source.append(call_volume) 496 496 497 refs = [" q[i]"] + _call_pars("v.", partable.iq_parameters)498 call_iq = "#define CALL_IQ( q,i,v) Iq(%s)" % (",".join(refs))497 refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters) 498 call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs)) 499 499 if _have_Iqxy(user_code): 500 500 # Call 2D model 501 refs = ["q[2*i]", "q[2*i+1]"] + _call_pars(" v.", partable.iqxy_parameters)502 call_iqxy = "#define CALL_IQ( q,i,v) Iqxy(%s)" % (",".join(refs))501 refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("_v.", partable.iqxy_parameters) 502 call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs)) 503 503 else: 504 504 # Call 1D model with sqrt(qx^2 + qy^2) 505 505 warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 506 506 # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 507 pars_sqrt = ["sqrt( q[2*i]*q[2*i]+q[2*i+1]*q[2*i+1])"] + refs[1:]508 call_iqxy = "#define CALL_IQ( q,i,v) Iq(%s)" % (",".join(pars_sqrt))507 pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:] 508 call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 509 509 510 510 # Fill in definitions for numbers of parameters -
sasmodels/kernelcl.py
rba32cdd rd2bb604 55 55 56 56 try: 57 raise NotImplementedError("OpenCL not yet implemented for new kernel template") 57 58 import pyopencl as cl 58 59 # Ask OpenCL for the default context so that we know that one exists -
sasmodels/model_test.py
r4d76711 rd2bb604 89 89 90 90 if is_py: # kernel implemented in python 91 continue # TODO: re-enable python tests 91 92 test_name = "Model: %s, Kernel: python"%model_name 92 93 test_method_name = "test_%s_python" % model_name -
sasmodels/modelinfo.py
rfb5914f rd2bb604 144 144 # support for the form 145 145 # demo(thickness=0, thickness2=50) 146 for k in pars[name].length:146 for k in range(1, pars[name].length+1): 147 147 key = name+str(k) 148 148 if key not in scalars: -
sasmodels/models/flexible_cylinder.c
r43b7eea rd2bb604 1 double form_volume(double length, double kuhn_length, double radius); 2 double Iq(double q, double length, double kuhn_length, double radius, 3 double sld, double solvent_sld); 4 double Iqxy(double qx, double qy, double length, double kuhn_length, 5 double radius, double sld, double solvent_sld); 6 double flexible_cylinder_kernel(double q, double length, double kuhn_length, 7 double radius, double sld, double solvent_sld); 8 9 10 double form_volume(double length, double kuhn_length, double radius) 1 static double 2 form_volume(length, kuhn_length, radius) 11 3 { 12 13 return 0.0; 4 return 1.; 14 5 } 15 6 16 double flexible_cylinder_kernel(double q, 17 double length, 18 double kuhn_length, 19 double radius, 20 double sld, 21 double solvent_sld) 7 static double 8 Iq(double q, 9 double length, 10 double kuhn_length, 11 double radius, 12 double sld, 13 double solvent_sld) 22 14 { 23 24 const double cont = sld-solvent_sld; 25 const double qr = q*radius; 26 //const double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr); 27 const double crossSect = sas_J1c(qr); 28 double flex = Sk_WR(q,length,kuhn_length); 29 flex *= crossSect*crossSect; 30 flex *= M_PI*radius*radius*length; 31 flex *= cont*cont; 32 flex *= 1.0e-4; 33 return flex; 15 const double contrast = sld - solvent_sld; 16 const double cross_section = sas_J1c(q*radius); 17 const double volume = M_PI*radius*radius*length; 18 const double flex = Sk_WR(q, length, kuhn_length); 19 return 1.0e-4 * volume * square(contrast*cross_section) * flex; 34 20 } 35 36 double Iq(double q,37 double length,38 double kuhn_length,39 double radius,40 double sld,41 double solvent_sld)42 {43 44 double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);45 return result;46 }47 48 double Iqxy(double qx, double qy,49 double length,50 double kuhn_length,51 double radius,52 double sld,53 double solvent_sld)54 {55 double q;56 q = sqrt(qx*qx+qy*qy);57 double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);58 59 return result;60 } -
sasmodels/models/hardsphere.py
rec45c4f rd2bb604 149 149 """ 150 150 151 Iqxy = """152 // never called since no orientation or magnetic parameters.153 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);154 """155 156 151 # ER defaults to 0.0 157 152 # VR defaults to 1.0 -
sasmodels/models/hayter_msa.py
rec45c4f rd2bb604 87 87 return 1.0; 88 88 """ 89 Iqxy = """90 // never called since no orientation or magnetic parameters.91 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);92 """93 89 # ER defaults to 0.0 94 90 # VR defaults to 1.0 -
sasmodels/models/lamellar.py
rec45c4f rd2bb604 82 82 """ 83 83 84 Iqxy = """85 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);86 """87 88 84 # ER defaults to 0.0 89 85 # VR defaults to 1.0 -
sasmodels/models/lamellar_hg.py
rec45c4f rd2bb604 101 101 """ 102 102 103 Iqxy = """104 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);105 """106 107 103 # ER defaults to 0.0 108 104 # VR defaults to 1.0 -
sasmodels/models/lamellar_hg_stack_caille.py
rec45c4f rd2bb604 120 120 """ 121 121 122 Iqxy = """123 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);124 """125 126 122 # ER defaults to 0.0 127 123 # VR defaults to 1.0 -
sasmodels/models/lamellar_stack_caille.py
rec45c4f rd2bb604 104 104 """ 105 105 106 Iqxy = """107 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);108 """109 110 106 # ER defaults to 0.0 111 107 # VR defaults to 1.0 -
sasmodels/models/lamellar_stack_paracrystal.py
rec45c4f rd2bb604 132 132 """ 133 133 134 Iqxy = """135 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);136 """137 138 134 # ER defaults to 0.0 139 135 # VR defaults to 1.0 -
sasmodels/models/rpa.c
r13ed84c rd2bb604 1 1 double Iq(double q, double case_num, 2 double Na, double Phia, double va, double a_sld, double ba, 3 double Nb, double Phib, double vb, double b_sld, double bb, 4 double Nc, double Phic, double vc, double c_sld, double bc, 5 double Nd, double Phid, double vd, double d_sld, double bd, 2 double N[], double Phi[], double v[], double L[], double b[], 6 3 double Kab, double Kac, double Kad, 7 4 double Kbc, double Kbd, double Kcd 8 5 ); 9 6 10 double Iqxy(double qx, double qy, double case_num,11 double Na, double Phia, double va, double a_sld, double ba,12 double Nb, double Phib, double vb, double b_sld, double bb,13 double Nc, double Phic, double vc, double c_sld, double bc,14 double Nd, double Phid, double vd, double d_sld, double bd,15 double Kab, double Kac, double Kad,16 double Kbc, double Kbd, double Kcd17 );18 19 double form_volume(void);20 21 double form_volume(void)22 {23 return 1.0;24 }25 26 7 double Iq(double q, double case_num, 27 double Na, double Phia, double va, double La, double ba, 28 double Nb, double Phib, double vb, double Lb, double bb, 29 double Nc, double Phic, double vc, double Lc, double bc, 30 double Nd, double Phid, double vd, double Ld, double bd, 8 double N[], double Phi[], double v[], double L[], double b[], 31 9 double Kab, double Kac, double Kad, 32 10 double Kbc, double Kbd, double Kcd … … 36 14 #if 0 // Sasview defaults 37 15 if (icase <= 1) { 38 N a=Nb=1000.0;39 Phi a=Phib=0.0000001;16 N[0]=N[1]=1000.0; 17 Phi[0]=Phi[1]=0.0000001; 40 18 Kab=Kac=Kad=Kbc=Kbd=-0.0004; 41 L a=Lb=1e-12;42 v a=vb=100.0;43 b a=bb=5.0;19 L[0]=L[1]=1e-12; 20 v[0]=v[1]=100.0; 21 b[0]=b[1]=5.0; 44 22 } else if (icase <= 4) { 45 Phi a=0.0000001;23 Phi[0]=0.0000001; 46 24 Kab=Kac=Kad=-0.0004; 47 L a=1e-12;48 v a=100.0;49 b a=5.0;25 L[0]=1e-12; 26 v[0]=100.0; 27 b[0]=5.0; 50 28 } 51 29 #else 52 30 if (icase <= 1) { 53 N a=Nb=0.0;54 Phi a=Phib=0.0;31 N[0]=N[1]=0.0; 32 Phi[0]=Phi[1]=0.0; 55 33 Kab=Kac=Kad=Kbc=Kbd=0.0; 56 L a=Lb=Ld;57 v a=vb=vd;58 b a=bb=0.0;34 L[0]=L[1]=L[3]; 35 v[0]=v[1]=v[3]; 36 b[0]=b[1]=0.0; 59 37 } else if (icase <= 4) { 60 N a= 0.0;61 Phi a=0.0;38 N[0] = 0.0; 39 Phi[0]=0.0; 62 40 Kab=Kac=Kad=0.0; 63 L a=Ld;64 v a=vd;65 b a=0.0;41 L[0]=L[3]; 42 v[0]=v[3]; 43 b[0]=0.0; 66 44 } 67 45 #endif 68 46 69 const double Xa = q*q*b a*ba*Na/6.0;70 const double Xb = q*q*b b*bb*Nb/6.0;71 const double Xc = q*q*b c*bc*Nc/6.0;72 const double Xd = q*q*b d*bd*Nd/6.0;47 const double Xa = q*q*b[0]*b[0]*N[0]/6.0; 48 const double Xb = q*q*b[1]*b[1]*N[1]/6.0; 49 const double Xc = q*q*b[2]*b[2]*N[2]/6.0; 50 const double Xd = q*q*b[3]*b[3]*N[3]/6.0; 73 51 74 52 // limit as Xa goes to 0 is 1 … … 98 76 #if 0 99 77 const double S0aa = icase<5 100 ? 1.0 : N a*Phia*va*Paa;78 ? 1.0 : N[0]*Phi[0]*v[0]*Paa; 101 79 const double S0bb = icase<2 102 ? 1.0 : N b*Phib*vb*Pbb;103 const double S0cc = N c*Phic*vc*Pcc;104 const double S0dd = N d*Phid*vd*Pdd;80 ? 1.0 : N[1]*Phi[1]*v[1]*Pbb; 81 const double S0cc = N[2]*Phi[2]*v[2]*Pcc; 82 const double S0dd = N[3]*Phi[3]*v[3]*Pdd; 105 83 const double S0ab = icase<8 106 ? 0.0 : sqrt(N a*va*Phia*Nb*vb*Phib)*Pa*Pb;84 ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 107 85 const double S0ac = icase<9 108 ? 0.0 : sqrt(N a*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);86 ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 109 87 const double S0ad = icase<9 110 ? 0.0 : sqrt(N a*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);88 ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 111 89 const double S0bc = (icase!=4 && icase!=7 && icase!= 9) 112 ? 0.0 : sqrt(N b*vb*Phib*Nc*vc*Phic)*Pb*Pc;90 ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 113 91 const double S0bd = (icase!=4 && icase!=7 && icase!= 9) 114 ? 0.0 : sqrt(N b*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);92 ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 115 93 const double S0cd = (icase==0 || icase==2 || icase==5) 116 ? 0.0 : sqrt(N c*vc*Phic*Nd*vd*Phid)*Pc*Pd;94 ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 117 95 #else // sasview equivalent 118 //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N c,Phic,vc,Pcc);119 double S0aa = N a*Phia*va*Paa;120 double S0bb = N b*Phib*vb*Pbb;121 double S0cc = N c*Phic*vc*Pcc;122 double S0dd = N d*Phid*vd*Pdd;123 double S0ab = sqrt(N a*va*Phia*Nb*vb*Phib)*Pa*Pb;124 double S0ac = sqrt(N a*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);125 double S0ad = sqrt(N a*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);126 double S0bc = sqrt(N b*vb*Phib*Nc*vc*Phic)*Pb*Pc;127 double S0bd = sqrt(N b*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);128 double S0cd = sqrt(N c*vc*Phic*Nd*vd*Phid)*Pc*Pd;96 //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc); 97 double S0aa = N[0]*Phi[0]*v[0]*Paa; 98 double S0bb = N[1]*Phi[1]*v[1]*Pbb; 99 double S0cc = N[2]*Phi[2]*v[2]*Pcc; 100 double S0dd = N[3]*Phi[3]*v[3]*Pdd; 101 double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 102 double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 103 double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 104 double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 105 double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 106 double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 129 107 switch(icase){ 130 108 case 0: … … 311 289 // Note: 1e-13 to convert from fm to cm for scattering length 312 290 const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13; 313 const double Lad = icase<5 ? 0.0 : (L a/va - Ld/vd)*sqrt_Nav;314 const double Lbd = icase<2 ? 0.0 : (L b/vb - Ld/vd)*sqrt_Nav;315 const double Lcd = (L c/vc - Ld/vd)*sqrt_Nav;291 const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav; 292 const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav; 293 const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav; 316 294 317 295 const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33 … … 321 299 322 300 } 323 324 double Iqxy(double qx, double qy,325 double case_num,326 double Na, double Phia, double va, double a_sld, double ba,327 double Nb, double Phib, double vb, double b_sld, double bb,328 double Nc, double Phic, double vc, double c_sld, double bc,329 double Nd, double Phid, double vd, double d_sld, double bd,330 double Kab, double Kac, double Kad,331 double Kbc, double Kbd, double Kcd332 )333 {334 double q = sqrt(qx*qx + qy*qy);335 return Iq(q,336 case_num,337 Na, Phia, va, a_sld, ba,338 Nb, Phib, vb, b_sld, bb,339 Nc, Phic, vc, c_sld, bc,340 Nd, Phid, vd, d_sld, bd,341 Kab, Kac, Kad,342 Kbc, Kbd, Kcd);343 } -
sasmodels/models/spherical_sld.py
rec45c4f rd2bb604 170 170 # pylint: disable=bad-whitespace, line-too-long 171 171 # ["name", "units", default, [lower, upper], "type", "description"], 172 parameters = [["n _shells","", 1, [0, 9], "", "number of shells"],172 parameters = [["n", "", 1, [0, 9], "", "number of shells"], 173 173 ["radius_core", "Ang", 50.0, [0, inf], "", "intern layer thickness"], 174 174 ["sld_core", "1e-6/Ang^2", 2.07, [-inf, inf], "", "sld function flat"], … … 192 192 193 193 demo = dict( 194 n _shells=4,194 n=4, 195 195 scale=1.0, 196 196 solvent_sld=1.0, -
sasmodels/models/squarewell.py
rec45c4f rd2bb604 123 123 """ 124 124 125 Iqxy = """126 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);127 """128 129 125 # ER defaults to 0.0 130 126 # VR defaults to 1.0 -
sasmodels/models/stickyhardsphere.py
rec45c4f rd2bb604 171 171 """ 172 172 173 Iqxy = """174 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);175 """176 177 173 # ER defaults to 0.0 178 174 # VR defaults to 1.0 -
sasmodels/resolution.py
raaf75b6 rd2bb604 502 502 from scipy.integrate import romberg 503 503 504 par_set = set([p.name for p in form.info['parameters'] ])504 par_set = set([p.name for p in form.info['parameters'].call_parameters]) 505 505 if any(k not in par_set for k in pars.keys()): 506 506 extra = set(pars.keys()) - par_set 507 raise ValueError("bad parameters: [%s] not in [%s]"% 508 (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 507 raise ValueError("bad parameters: [%s] not in [%s]" 508 % (", ".join(sorted(extra)), 509 ", ".join(sorted(pars.keys())))) 509 510 510 511 if np.isscalar(width): … … 556 557 from scipy.integrate import romberg 557 558 558 par_set = set([p.name for p in form.info['parameters'] ])559 par_set = set([p.name for p in form.info['parameters'].call_parameters]) 559 560 if any(k not in par_set for k in pars.keys()): 560 561 extra = set(pars.keys()) - par_set 561 raise ValueError("bad parameters: [%s] not in [%s]" %562 (", ".join(sorted(extra)),563 ", ".join(sorted(pars.keys()))))562 raise ValueError("bad parameters: [%s] not in [%s]" 563 % (", ".join(sorted(extra)), 564 ", ".join(sorted(pars.keys())))) 564 565 565 566 _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) -
sasmodels/sasview_model.py
r1e2a1ba rd2bb604 346 346 """ 347 347 if self._model is None: 348 self._model = core.build_model(self._model_info , platform='dll')348 self._model = core.build_model(self._model_info) 349 349 q_vectors = [np.asarray(q) for q in args] 350 350 kernel = self._model.make_kernel(q_vectors) … … 352 352 for p in self._model_info['parameters'].call_parameters] 353 353 details, weights, values = core.build_details(kernel, pairs) 354 re turnkernel(details, weights, values, cutoff=self.cutoff)354 result = kernel(details, weights, values, cutoff=self.cutoff) 355 355 kernel.q_input.release() 356 356 kernel.release()
Note: See TracChangeset
for help on using the changeset viewer.