diff --git a/sasmodels/compare.py b/sasmodels/compare.py index 081e2589..57370c75 100755 --- a/sasmodels/compare.py +++ b/sasmodels/compare.py @@ -80,7 +80,7 @@ def __call__(self, **par: float) -> np.ndarray: ... === model parameters === -preset*/-random[=seed] preset or random parameters - -sets=n generates n random datasets with the seed given by -random=seed + -sets=n generates n random datasets using -seed=seed -pars/-nopars* prints the parameter set or not -sphere[=150] set up spherical integration over theta/phi using n points -mono*/-poly suppress or allow polydispersity on generated parameters @@ -1023,7 +1023,7 @@ def plot_models(opts, result, limits=None, setnum=0): '2d', '1d', 'sesans', # Parameter set - 'preset', 'random', 'random=', 'sets=', + 'preset', 'random', 'random=', 'sets=', 'seed=', 'nopars', 'pars', 'sphere', 'sphere=', # integrate over a sphere in 2d with n points 'poly', 'mono', @@ -1227,6 +1227,7 @@ def parse_opts(argv): elif arg.startswith('-res='): opts['res'] = arg[5:] elif arg.startswith('-noise='): opts['noise'] = float(arg[7:]) elif arg.startswith('-sets='): opts['sets'] = int(arg[6:]) + elif arg.startswith('-seed='): opts['seed'] = int(arg[6:]) elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:] elif arg.startswith('-cutoff='): opts['cutoff'] = arg[8:] elif arg.startswith('-title='): opts['title'] = arg[7:] diff --git a/sasmodels/direct_model.py b/sasmodels/direct_model.py index 7662e17a..2fad3406 100644 --- a/sasmodels/direct_model.py +++ b/sasmodels/direct_model.py @@ -453,14 +453,18 @@ def test_reparameterize(): except Exception: pass -def _direct_calculate(model, data, pars): +def _direct_calculate(model, data, pars, ngauss=0): from .core import build_model, load_model_info + from .generate import set_integration_size + model_info = load_model_info(model) + if ngauss != 0: + set_integration_size(model_info, ngauss) kernel = build_model(model_info) calculator = DirectModel(data, kernel) return calculator(**pars) -def Iq(model, q, dq=None, ql=None, qw=None, **pars): +def Iq(model, q, dq=None, ql=None, qw=None, ngauss=0, **pars): """ Compute I(q) for *model*. Resolution is *dq* for pinhole or *ql* and *qw* for slit geometry. Use 0 or None for infinite slits. @@ -498,16 +502,16 @@ def broadcast(v): else np.full(len(q), v) if np.isscalar(v) else _as_numpy(v)) data.dxl, data.dxw = broadcast(ql), broadcast(qw) - return _direct_calculate(model, data, pars) + return _direct_calculate(model, data, pars, ngauss=ngauss) -def Iqxy(model, qx, qy, dqx=None, dqy=None, **pars): +def Iqxy(model, qx, qy, dqx=None, dqy=None, ngauss=0, **pars): """ Compute I(qx, qy) for *model*. Resolution is *dqx* and *dqy*. See :func:`Iq` for details on model and parameters. """ from .data import Data2D data = Data2D(x=qx, y=qy, dx=dqx, dy=dqy) - return _direct_calculate(model, data, pars) + return _direct_calculate(model, data, pars, ngauss=ngauss) def Gxi(model, xi, **pars): """ @@ -528,6 +532,8 @@ def main(): if len(sys.argv) < 3: print("usage: python -m sasmodels.direct_model modelname (q|qx,qy) par=val ...") sys.exit(1) + + ngauss = 0 model = sys.argv[1] call = sys.argv[2].upper() pars = dict((k, (float(v) if not k.endswith("_pd_type") else v)) @@ -542,13 +548,13 @@ def main(): dq = dqw = dql = None #dq = [q*0.05] # 5% pinhole resolution #dqw, dql = [q*0.05], [1.0] # 5% horizontal slit resolution - print(Iq(model, [q], dq=dq, qw=dqw, ql=dql, **pars)[0]) + print(Iq(model, [q], dq=dq, qw=dqw, ql=dql, ngauss=ngauss, **pars)[0]) #print(Gxi(model, [q], **pars)[0]) elif len(values) == 2: qx, qy = values dq = None #dq = [0.005] # 5% pinhole resolution at q = 0.1 - print(Iqxy(model, [qx], [qy], dqx=dq, dqy=dq, **pars)[0]) + print(Iqxy(model, [qx], [qy], dqx=dq, dqy=dq, ngauss=ngauss, **pars)[0]) else: print("use q or qx,qy") sys.exit(1) diff --git a/sasmodels/generate.py b/sasmodels/generate.py index 7db1c3ec..55f395c6 100644 --- a/sasmodels/generate.py +++ b/sasmodels/generate.py @@ -291,13 +291,27 @@ def set_integration_size(info, n): Note: this really ought to be a method in modelinfo, but that leads to import loops. """ - if info.source and any(lib.startswith('lib/gauss') for lib in info.source): - from .gengauss import gengauss - path = joinpath(MODEL_PATH, "lib", "gauss%d.c"%n) - if not exists(path): - gengauss(n, path) - info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss') - else lib for lib in info.source] + from .gengauss import gengauss + + if not info.source: + return + + # Generate the integration points + path = joinpath(MODEL_PATH, "lib", f"gauss{n}.c") + if not exists(path): + # print(f"building Gaussian integration points of size {n} in {str(path)}") + gengauss(n, path) + + # Replace adaptive.c or lib/gauss.c + try: + index = info.source.index("lib/adaptive.c") + info.source[index:index+1] = [f"lib/gauss{n}.c", "lib/nonadaptive.c"] + except ValueError: + for index in range(len(info.source)-1, -1, -1): + if info.source[index].startswith("lib/gauss"): + info.source[index] = f"lib/gauss{n}.c" + break + # print("info.source is now", info.source) def format_units(units): # type: (str) -> str diff --git a/sasmodels/gengauss.py b/sasmodels/gengauss.py index a9ce027f..2ef06313 100755 --- a/sasmodels/gengauss.py +++ b/sasmodels/gengauss.py @@ -21,19 +21,18 @@ def gengauss(n, path): array_size = n with open(path, "w") as fid: - fid.write("""\ -// Generated by sasmodels.gengauss.gengauss(%d) + fid.write(f"""\ +// Generated by sasmodels.gengauss.gengauss({n}) #ifdef GAUSS_N # undef GAUSS_N # undef GAUSS_Z # undef GAUSS_W #endif -#define GAUSS_N %d -#define GAUSS_Z Gauss%dZ -#define GAUSS_W Gauss%dWt - -"""%(n, n, n, n)) +#define GAUSS_N {n} +#define GAUSS_Z Gauss{n}Z +#define GAUSS_W Gauss{n}Wt +""") if array_size != n: fid.write("// Note: using array size %d so that it is a multiple of 4\n\n"%array_size) diff --git a/sasmodels/model_test.py b/sasmodels/model_test.py index ef9b079b..4d0ca340 100755 --- a/sasmodels/model_test.py +++ b/sasmodels/model_test.py @@ -579,6 +579,31 @@ def _build_test(test): for test in tests: yield _build_test(test) +def _generate_target_values(modelname, ngauss=0): + from .generate import set_integration_size + + model_info = load_model_info(modelname) + if ngauss != 0: + set_integration_size(model_info, ngauss) + model = build_model(model_info, platform="dll", dtype="d") + + for pars, q, Iq in model_info.tests: + qin = q + if isinstance(Iq, float): + q, Iq = [q], [Iq] + if isinstance(q[0], tuple): + qx, qy = zip(*q) + q_vectors = [np.array(qx), np.array(qy)] + else: + q_vectors = [np.array(q)] + kernel = model.make_kernel(q_vectors) + target = np.array(Iq) + actual = call_kernel(kernel, pars) + if True or (actual != target).any(): + print("Test:", modelname, pars) + print(" q = ", qin) + print(f" current => [{', '.join(f'{v:.15g}' for v in target)}]") + print(f" ngauss={ngauss} => [{', '.join(f'{v:.15g}' for v in actual)}]") def main(): # type: () -> int @@ -601,6 +626,11 @@ def main(): help="Engines on which to run the test. " "Valid values are opencl, cuda, dll, and all. " "Defaults to all if no value is given") + parser.add_argument("-t", "--targets", action="store_true", + help="Generate target values for test.") + parser.add_argument("--ngauss", type=int, default=10000, + help="Number of gauss points to use in integration for " + "target values. Warning: this is very slow the first time.") parser.add_argument("models", nargs="*", help="The names of the models to be tested. " "If the first model is 'all', then all but the listed " @@ -630,9 +660,13 @@ def main(): print("unknown engine " + opts.engine) return 1 - runner = TestRunner(verbosity=opts.verbose, **test_args) - result = runner.run(make_suite(loaders, opts.models)) - return 1 if result.failures or result.errors else 0 + if opts.targets: + for model in opts.models: + _generate_target_values(model, ngauss=opts.ngauss) + else: + runner = TestRunner(verbosity=opts.verbose, **test_args) + result = runner.run(make_suite(loaders, opts.models)) + return 1 if result.failures or result.errors else 0 if __name__ == "__main__": diff --git a/sasmodels/models/barbell.c b/sasmodels/models/barbell.c index 87f1553e..c1344d0f 100644 --- a/sasmodels/models/barbell.c +++ b/sasmodels/models/barbell.c @@ -19,13 +19,18 @@ _bell_kernel(double qab, double qc, double h, double radius_bell, const double m = radius_bell*qc; // cos argument slope const double b = (half_length+h)*qc; // cos argument intercept const double qab_r = radius_bell*qab; // Q*R*sin(theta) + + const double qr_max = fmax(qab_r, m); + constant double *w, *z; + int n = gauss_weights(qr_max, &w, &z); + double total = 0.0; - for (int i = 0; i < GAUSS_N; i++){ - const double t = GAUSS_Z[i]*zm + zb; + for (int i = 0; i < n; i++){ + const double t = z[i]*zm + zb; const double radical = 1.0 - t*t; const double bj = sas_2J1x_x(qab_r*sqrt(radical)); const double Fq = cos(m*t + b) * radical * bj; - total += GAUSS_W[i] * Fq; + total += w[i] * Fq; } // translate dx in [-1,1] to dx in [lower,upper] const double integral = total*zm; @@ -110,21 +115,25 @@ Fq(double q,double *F1, double *F2, double sld, double solvent_sld, const double h = sqrt(square(radius_bell) - square(radius)); const double half_length = 0.5*length; + const double qr_max = q*fmax(radius, half_length); + constant double *w, *z; + int n = gauss_weights(qr_max, &w, &z); + // translate a point in [-1,1] to a point in [0, pi/2] const double zm = M_PI_4; const double zb = M_PI_4; double total_F1 = 0.0; double total_F2 = 0.0; - for (int i = 0; i < GAUSS_N; i++){ - const double theta = GAUSS_Z[i]*zm + zb; + for (int i = 0; i < n; i++){ + const double theta = z[i]*zm + zb; double sin_theta, cos_theta; // slots to hold sincos function output SINCOS(theta, sin_theta, cos_theta); const double qab = q*sin_theta; const double qc = q*cos_theta; const double Aq = _fq(qab, qc, h, radius_bell, radius, half_length); // scale by sin_theta for spherical coord integration - total_F1 += GAUSS_W[i] * Aq * sin_theta; - total_F2 += GAUSS_W[i] * Aq * Aq * sin_theta; + total_F1 += w[i] * Aq * sin_theta; + total_F2 += w[i] * Aq * Aq * sin_theta; } // translate dx in [-1,1] to dx in [lower,upper] const double form_avg = total_F1 * zm; diff --git a/sasmodels/models/barbell.py b/sasmodels/models/barbell.py index d7f3d307..4f3fbcb7 100644 --- a/sasmodels/models/barbell.py +++ b/sasmodels/models/barbell.py @@ -117,7 +117,7 @@ ] # pylint: enable=bad-whitespace, line-too-long -source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"] +source = ["lib/polevl.c", "lib/sas_J1.c", "lib/adaptive.c", "barbell.c"] valid = "radius_bell >= radius" have_Fq = True radius_effective_modes = [ diff --git a/sasmodels/models/capped_cylinder.c b/sasmodels/models/capped_cylinder.c index c792574c..80c23a18 100644 --- a/sasmodels/models/capped_cylinder.c +++ b/sasmodels/models/capped_cylinder.c @@ -26,13 +26,18 @@ _cap_kernel(double qab, double qc, double h, double radius_cap, const double m = radius_cap*qc; // cos argument slope const double b = (half_length+h)*qc; // cos argument intercept const double qab_r = radius_cap*qab; // Q*R*sin(theta) + + const double qr_max = fmax(qab_r, m); + constant double *w, *z; + int n = gauss_weights(qr_max, &w, &z); + double total = 0.0; - for (int i=0; i [0, 1] const double m = 0.5; const double b = 0.5; double total_F1 = 0.0; //initialize intergral double total_F2 = 0.0; //initialize intergral - for(int i=0;i