Skip to content

Remove langevin from C library #624

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
krzywon merged 2 commits into from
Feb 11, 2025
Merged

Remove langevin from C library #624

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3b9b17e
Remove langevin from C library. Fixes #612
Jan 17, 2025
2711eba
Merge branch 'master' into 612_remove_langevin
pkienzle Jan 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 20 additions & 3 deletions explore/precision.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -413,6 +413,7 @@ def add_function(name, mp_function, np_function, ocl_function,
ocl_function=make_ocl("return fmod(q, 2*M_PI);", "sas_fmod"),
)

# TODO: move to sas_special
def sas_langevin(x):
scalar = np.isscalar(x)
if scalar:
Expand Down Expand Up @@ -451,14 +452,30 @@ def sas_langevin_x(x):
mp_function=lambda x: (1/mp.tanh(x) - 1/x),
np_function=sas_langevin,
#ocl_function=make_ocl("return q < 0.7 ? q*(1./3. + q*q*(-1./45. + q*q*(2./945. + q*q*(-1./4725.) + q*q*(2./93555.)))) : 1/tanh(q) - 1/q;", "sas_langevin"),
ocl_function=make_ocl("return q < 1e-5 ? q/3. : 1/tanh(q) - 1/q;", "sas_langevin"),
#ocl_function=make_ocl("return q < 1e-5 ? q/3. : 1/tanh(q) - 1/q;", "sas_langevin"),
ocl_function=make_ocl("""
#if FLOAT_SIZE>4 // DOUBLE_PRECISION
# define LANGEVIN_CUTOFF 0.1
#else
# define LANGEVIN_CUTOFF 1.0
#endif
const double qsq = q*q;
return (q < LANGEVIN_CUTOFF) ? q / (3. + qsq / (5. + qsq/(7. + qsq/(9.)))) : 1/tanh(q) - 1/q;
""", "sas_langevin"),
)
add_function(
name="langevin(x)/x",
mp_function=lambda x: (1/mp.tanh(x) - 1/x)/x,
#np_function=lambda x: sas_langevin(x)/x, # Note: need to test for x=0
np_function=sas_langevin_x,
ocl_function=make_ocl("return q < 1e-5 ? 1./3. : (1/tanh(q) - 1/q)/q;", "sas_langevin_x"),
ocl_function=make_ocl("""
#if FLOAT_SIZE>4 // DOUBLE_PRECISION
# define LANGEVIN_CUTOFF 0.1
#else
# define LANGEVIN_CUTOFF 1.0
#endif
const double qsq = q*q;
return (q < LANGEVIN_CUTOFF) ? 1. / (3. + qsq / (5. + qsq/(7. + qsq/(9.)))) : (1/tanh(q) - 1/q)/q;
""", "sas_langevin_x"),
)
add_function(
name="gauss_coil",
Expand Down
28 changes: 0 additions & 28 deletions sasmodels/models/lib/magnetic_functions.c
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,34 +1,6 @@
//These functions are required for magnetic analysis models. They are copies
//from sasmodels/kernel_iq.c, to enables magnetic parameters for 1D and 2D models.


static double langevin(
double x) {
// cotanh(x)-1\x

if (x < 0.00001) {
// avoid dividing by zero
return 1.0/3.0*x-1.0/45.0 * pow(x, 3) + 2.0/945.0 * pow(x, 5) - 1.0/4725.0 * pow(x, 7);
} else {
return 1.0/tanh(x)-1/x;
}
}

static double langevinoverx(
double x)
{
// cotanh(x)-1\x

if (x < 0.00001) {
// avoid dividing by zero
return 1.0/3.0-1.0/45.0 * pow(x, 2) + 2.0/945.0 * pow(x, 4) - 1.0/4725.0 * pow(x, 6);
} else {
return langevin(x)/x;
}
}



//weighting of spin resolved cross sections to reconstruct partially polarised beam with imperfect optics using up_i/up_f.
static void set_weights(double in_spin, double out_spin, double weight[8]) //from kernel_iq.c
{
Expand Down
Loading