Merge pull request #592 from SasView/insert-marketplacemodel-bulk-fer…

…romagnets

Insert marketplacemodel bulk ferromagnets

.github/workflows/test.yml

@@ -29,11 +29,7 @@ jobs:
     - name: setup apt dependencies for Linux
       if: ${{ matrix.os == 'ubuntu-latest' }}
       run: |
-        wget https://repo.radeon.com/amdgpu-install/6.1.1/ubuntu/jammy/amdgpu-install_6.1.60101-1_all.deb
-        sudo apt install ./amdgpu-install_6.1.60101-1_all.deb
         sudo apt update
-        sudo apt install amdgpu-dkms
-        sudo apt install rocm
 
     - name: Install Python dependencies
       run: |
@@ -56,9 +52,11 @@ jobs:
         choco install opencl-intel-cpu-runtime
         python -m pip install --only-binary=pyopencl --find-links http://www.silx.org/pub/wheelhouse/ --trusted-host www.silx.org pyopencl
 
-    - name: Test with pytest
+    - name: Test with pytest (only on Windows for now since PoCL is failing on Ubuntu)
+
       env:
         PYOPENCL_COMPILER_OUTPUT: 1
+        SAS_OPENCL: none
       run: |
         # other CI uses the following, but `setup.py test` is a deprecated way
         # of running tests
@@ -68,6 +66,8 @@ jobs:
 
     - name: check that the docs build (linux only)
       if: ${{ matrix.os == 'ubuntu-latest' }}
+      env:
+        SAS_OPENCL: none
       run: |
         make -j 4 -C doc SPHINXOPTS="-W --keep-going -n" html
 

CHANGES.rst

@@ -1,10 +1,12 @@
 Release notes
 =============
 
-v1.0.8 2024-06-??
+v1.0.8 2024-09-26
 -----------------
+* New model: Bulk ferromagnets model from marketplace
 * Doc update: Archive built docs on Github
 * Doc update: Display math correctly
+* Fix error in FCC paracrystaline models
 * Fix parameter name checking in kernel call
 
 v1.0.7 2023-03-23

doc/guide/plugin.rst

@@ -851,6 +851,8 @@ Some non-standard constants and functions are also provided:
         $x^2$
     cube(x):
         $x^3$
+    clip(a, a_min, a_max):
+        $\min(\max(a, a_\text{min}), a_\text{max})$, or NaN if $a$ is NaN.
     sas_sinx_x(x):
         $\sin(x)/x$, with limit $\sin(0)/0 = 1$.
     powr(x, y):

explore/precision.py

@@ -190,13 +190,13 @@ def plotdiff(x, target, actual, label, diff):
     if diff == "relative":
         err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd')
         #err = np.clip(err, 0, 1)
-        pylab.loglog(x, err, '-', label=label)
+        pylab.loglog(x, err, '-', label=label, alpha=0.7)
     elif diff == "absolute":
         err = np.array([abs((t-a)) for t, a in zip(target, actual)], 'd')
-        pylab.loglog(x, err, '-', label=label)
+        pylab.loglog(x, err, '-', label=label, alpha=0.7)
     else:
         limits = np.min(target), np.max(target)
-        pylab.semilogx(x, np.clip(actual, *limits), '-', label=label)
+        pylab.semilogx(x, np.clip(actual, *limits), '-', label=label, alpha=0.7)
 
 def make_ocl(function, name, source=[]):
     class Kernel(object):
@@ -412,6 +412,54 @@ def add_function(name, mp_function, np_function, ocl_function,
     np_function=lambda x: np.fmod(x, 2*np.pi),
     ocl_function=make_ocl("return fmod(q, 2*M_PI);", "sas_fmod"),
 )
+
+def sas_langevin(x):
+    scalar = np.isscalar(x)
+    if scalar:
+        x = np.array([x]) # should inherit dtype for single if given single
+    f = np.empty_like(x)
+    cutoff = 0.1 if f.dtype == np.float64 else 1.0
+    #cutoff *= 10
+    index = x < cutoff
+    xp = x[index]
+    xpsq = xp*xp
+    f[index] = xp / (3. + xpsq / (5. + xpsq/(7. + xpsq/(9.))))
+    # 4 terms gets to 1e-7 single, 1e-14 double. Can get to 1e-15 double by adding
+    # another 4 terms and setting cutoff at 1.0. Not worthwhile. Instead we would
+    # need an expansion about x somewhere between 1 and 10 for the interval [0.1, 100.]
+    #f[index] = xp / (3. + xpsq / (5. + xpsq/(7. + xpsq/(9. + xpsq/(11.0 + xpsq/(13. + xpsq/(15. + xpsq/17.)))))))
+    xp = x[~index]
+    f[~index] = 1/np.tanh(xp) - 1/xp
+    return f[0] if scalar else f
+
+def sas_langevin_x(x):
+    scalar = np.isscalar(x)
+    if scalar:
+        x = np.array([x]) # should inherit dtype for single if given single
+    f = np.empty_like(x)
+    cutoff = 0.1 if f.dtype == np.float64 else 1.0
+    index = x < cutoff
+    xp = x[index]
+    xpsq = xp*xp
+    f[index] = 1. / (3. + xpsq / (5. + xpsq/(7. + xpsq/(9.))))
+    xp = x[~index]
+    f[~index] = (1/np.tanh(xp) - 1/xp)/xp
+    return f[0] if scalar else f
+
+add_function(
+    name="langevin(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"),
+)
+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"),
+)
 add_function(
     name="gauss_coil",
     mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4,

sasmodels/kernel_header.c

@@ -187,8 +187,54 @@
 #endif
 inline double square(double x) { return x*x; }
 inline double cube(double x) { return x*x*x; }
+// clip() follows numpy.clip() semantics, returning (x < low ? low : x > high ? high : x)
+// OpenCL/CUDA clamp() returns fmin(fmax(x, low), high)
+// C++(17) clamp() matches numpy.clip()
+// If x is NaN numpy.clip() returns NaN but OpenCL clamp() returns low.
+inline double clip(double x, double low, double high) { return x < low ? low : x > high ? high : x; }
 inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; }
 
+// vector algebra
+static void SET_VEC(double *vector, double v0, double v1, double v2)
+{
+    vector[0] = v0;
+    vector[1] = v1;
+    vector[2] = v2;
+}
+
+static void SCALE_VEC(double *vector, double a)
+{
+    vector[0] = a*vector[0];
+    vector[1] = a*vector[1];
+    vector[2] = a*vector[2];
+}
+
+static void ADD_VEC(double *result_vec, double *vec1, double *vec2)
+{
+    result_vec[0] = vec1[0] + vec2[0];
+    result_vec[1] = vec1[1] + vec2[1];
+    result_vec[2] = vec1[2] + vec2[2];
+}
+
+static double SCALAR_VEC( double *vec1, double *vec2)
+{
+    return vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2];
+}
+
+static double MAG_VEC( double *vec)
+{
+    return sqrt(SCALAR_VEC(vec,vec));
+}
+
+
+static void ORTH_VEC(double *result_vec, double *vec1, double *vec2)
+{
+    double scale =  SCALAR_VEC(vec1,vec2) / SCALAR_VEC(vec2,vec2);
+    result_vec[0] = vec1[0] - scale * vec2[0];
+    result_vec[1] = vec1[1] - scale * vec2[1];
+    result_vec[2] = vec1[2] - scale * vec2[2];
+}
+
 // CRUFT: support old style models with orientation received qx, qy and angles
 
 // To rotate from the canonical position to theta, phi, psi, first rotate by

sasmodels/kernel_iq.c

@@ -70,52 +70,8 @@ typedef union {
 #if defined(MAGNETIC) && NUM_MAGNETIC > 0
 // ===== Helper functions for magnetism =====
 
-// vector algebra
-void SET_VEC(double *vector, double v0, double v1, double v2)
-{
-    vector[0] = v0;
-    vector[1] = v1;
-    vector[2] = v2;
-}
-
-void SCALE_VEC(double *vector, double a)
-{
-    vector[0] = a*vector[0];
-    vector[1] = a*vector[1];
-    vector[2] = a*vector[2];
-}
 
-void ADD_VEC(double *result_vec, double *vec1, double *vec2)
-{
-    result_vec[0] = vec1[0] + vec2[0];
-    result_vec[1] = vec1[1] + vec2[1];
-    result_vec[2] = vec1[2] + vec2[2];
-}
-
-static double SCALAR_VEC( double *vec1, double *vec2)
-{
-    return vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2];
-}
 
-static double MAG_VEC( double *vec)
-{
-    return sqrt(SCALAR_VEC(vec,vec));
-}
-
-void ORTH_VEC(double *result_vec, double *vec1, double *vec2)
-{
-    double scale =  SCALAR_VEC(vec1,vec2) / SCALAR_VEC(vec2,vec2);
-    result_vec[0] = vec1[0] - scale * vec2[0];
-    result_vec[1] = vec1[1] - scale * vec2[1];
-    result_vec[2] = vec1[2] - scale * vec2[2];
-}
-
-
-// Return value restricted between low and high
-static double clip(double value, double low, double high)
-{
-  return (value < low ? low : (value > high ? high : value));
-}
 
 // Compute spin cross sections given in_spin and out_spin
 // To convert spin cross sections to sld b:
@@ -178,11 +134,10 @@ static double mag_sld(
 )
 {
   double Mvector[3];
-  double Pvector[3]; 
+  double Pvector[3];
   double qvector[3];
-  double rhom[3]; 
   double perpy[3];
-  double perpz[3];  
+  double perpz[3];
   double Mperp[3];
 
   const double qsq = sqrt(qx*qx + qy*qy); 

sasmodels/models/lib/magnetic_functions.c

@@ -0,0 +1,130 @@
+//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
+{
+ double norm=out_spin;
+
+
+   in_spin = clip(sqrt(square(in_spin)), 0.0, 1.0);//opencl has ambiguities for abs()
+   out_spin = clip(sqrt(square(out_spin)), 0.0, 1.0);
+
+   if (out_spin < 0.5){norm=1-out_spin;}
+
+// The norm is needed to make sure that the scattering cross sections are
+//correctly weighted, such that the sum of spin-resolved measurements adds up to
+// the unpolarised or half-polarised scattering cross section. No intensity weighting
+// needed on the incoming polariser side assuming that a user has normalised
+// to the incoming flux with polariser in for SANSPOl and unpolarised beam, respectively.
+
+
+   weight[0] = (1.0-in_spin) * (1.0-out_spin) / norm; // dd.real
+   weight[1] = weight[0]; // dd.imag
+   weight[2] = in_spin * out_spin / norm;             // uu.real
+   weight[3] = weight[2];             // uu.imag
+   weight[4] = (1.0-in_spin) * out_spin / norm;       // du.real
+   weight[5] = weight[4];       // du.imag 
+   weight[6] = in_spin * (1.0-out_spin) / norm;       // ud.real
+   weight[7] = weight[6];       // ud.imag
+ }
+
+
+
+//transforms scattering vector q in polarisation/magnetisation coordinate system
+ static void set_scatvec(double *qrot, double q,
+  double cos_theta, double sin_theta,
+  double alpha, double beta)
+{
+  double cos_alpha, sin_alpha;
+  double cos_beta, sin_beta;
+
+  SINCOS(alpha*M_PI/180, sin_alpha, cos_alpha);
+  SINCOS(beta*M_PI/180, sin_beta, cos_beta);
+  //field is defined along (0,0,1), orientation of detector
+  //is precessing in a cone around B with an inclination of theta
+
+  qrot[0] = q*(cos_alpha * cos_theta);
+  qrot[1] = q*(cos_theta * sin_alpha*sin_beta + 
+  cos_beta * sin_theta);
+  qrot[2] = q*(-cos_beta * cos_theta* sin_alpha + 
+  sin_beta * sin_theta);
+}
+
+
+
+//Evaluating the magnetic scattering vector (Halpern Johnson vector) for general orientation of q and collecting terms for the spin-resolved (POLARIS) cross sections. Mz is along the applied magnetic field direction, which is also the polarisation direction.
+static void mag_sld(
+   // 0=dd.real, 1=dd.imag, 2=uu.real, 3=uu.imag,  4=du.real, 5=du.imag,  6=ud.real, 7=ud.imag
+ double qx, double qy, double qz,
+  double mxreal, double mximag, double myreal,  double myimag, double mzreal,double mzimag, double nuc, double sld[8])
+{
+  double vector[3];
+  //The (transversal) magnetisation and hence the magnetic scattering sector is here a complex quantity. The spin-flip (magnetic) scattering amplitude is given with
+  //  MperpPperpQ \pm i MperpP  (Moon-Riste-Koehler Phys Rev 181, 920, 1969) with Mperp and MperpPperpQ the
+  //magnetisation scattering vector components perpendicular to the polarisation/field direction. Collecting terms in SF that are real (MperpPperpQreal + SCALAR_VEC(MperpPimag,qvector) )  and imaginary (MperpPperpQimag \pm SCALAR_VEC(MperpPreal,qvector) )
+  double Mvectorreal[3];
+  double Mvectorimag[3];
+  double Pvector[3];
+  double perpy[3];
+  double perpx[3]; 
+
+  double Mperpreal[3];
+  double Mperpimag[3];
+
+  const double q = sqrt(qx*qx + qy*qy + qz*qz);
+  SET_VEC(vector, qx/q, qy/q, qz/q); 
+
+  //Moon-Riste-Koehler notation choose z as pointing along field/polarisation axis
+  //totally different to what is used in SASview (historical reasons)
+  SET_VEC(Mvectorreal, mxreal, myreal, mzreal);
+  SET_VEC(Mvectorimag, mximag, myimag, mzimag);
+
+  SET_VEC(Pvector, 0, 0, 1);
+  SET_VEC(perpy, 0, 1, 0);  
+  SET_VEC(perpx, 1, 0, 0);   
+   //Magnetic scattering vector Mperp could be simplified like in Moon-Riste-Koehler
+   //leave the generic computation just to check
+  ORTH_VEC(Mperpreal, Mvectorreal, vector);
+  ORTH_VEC(Mperpimag, Mvectorimag, vector);
+
+
+   sld[0] = nuc - SCALAR_VEC(Pvector,Mperpreal); // dd.real => sld - D Pvector \cdot Mperp 
+   sld[1] = +SCALAR_VEC(Pvector,Mperpimag); //dd.imag  = nuc_img - SCALAR_VEC(Pvector,Mperpimg); nuc_img only exist for noncentrosymmetric nuclear structures; 
+   sld[2] = nuc + SCALAR_VEC(Pvector,Mperpreal);              // uu => sld + D Pvector \cdot Mperp
+   sld[3] = -SCALAR_VEC(Pvector,Mperpimag); //uu.imag
+
+   sld[4] = SCALAR_VEC(perpy,Mperpreal)+SCALAR_VEC(perpx,Mperpimag);       // du.real => real part along y + imaginary part along x
+   sld[5] = SCALAR_VEC(perpy,Mperpimag)-SCALAR_VEC(perpx,Mperpreal);       // du.imag => imaginary component along y - i *real part along x
+   sld[6] = SCALAR_VEC(perpy,Mperpreal)-SCALAR_VEC(perpx,Mperpimag);      // ud.real =>  real part along y - imaginary part along x
+   sld[7] = SCALAR_VEC(perpy,Mperpimag)+SCALAR_VEC(perpx,Mperpreal);         // ud.imag => imaginary component along y + i * real part along x
+
+ }