3 dpolarisation

[dehoni]

This branch contains the changes for a generalized 3D description of magnetic SANS (ticket #432). The necessary adaption in SASview are found under SasView/sasview#1714. For consistency, the orientation/directions of magnetization and polarization should be similar to the particle orientation.

[pkienzle]

I tried comparing the results against the existing SasView and found that the differ when up_phi=0. Is this intended? I may have messed up on the tests.

Looking at the code in kernel_iq.c, you dividing q by |q|^2 rather than |q|, which is suspicious. I think this only affects the imaginary portion of du/ud since the other calculations are using ORTH_VEC which does it's own norm.

I have some code efficiency comments once you are satisfied that the code is working as intended.

The code in explore/realspace.py does not run. The purpose of this code is to cross check the analytic solution against Monte Carlo sampling from a shape.

You should be able to run it using, for example:

python explore/realspace.py -q0.1 -d2 sphere up_phi=30 up_angle=45 rho_m=3 theta_m=-10 phi_m=-20

I can get it to run with the following change to magnetic_sld:

    Mperp = Mvector-norm*np.dot(Mvector.T, qvector)[None, :]*qvector[:, None]
    MperpP = Mperp-(np.dot(Mperp.T, Pvector)[None, :]*Pvector[:, None])
    MperpPperpQ = MperpP- norm*np.dot(MperpP.T, qvector)[None, :]*qvector[:, None]
    return [
        rho - np.dot(Pvector,Mperp),   # dd => sld - D M_perpx
        np.sqrt(np.sum(MperpPperpQ**2, axis=0)) - 1j*np.dot(MperpP.T,qvector), # du => -D (M_perpy + j M_perpz)
        np.sqrt(np.sum(MperpPperpQ**2, axis=0)) + 1j*np.dot(MperpP.T,qvector), # ud => -D (M_perpy - j M_perpz)
        rho + np.dot(Pvector,Mperp),   # uu => sld + D M_perpx
    ]

but the result is nowhere close to matching the new magnetic sphere model. You appear to be applying the norm twice in this code compared to once in kernel_iq.c, so not surprising.

Following the code in kernel_iq.c, which defines ortho(a,b) = a - (a.b)/(b.b) b you can rewrite realspace.py using:

def orth(A, b_hat): # A = 3 x n, and b_hat unit vector
    return A - np.sum(A*b_hat[:, None], axis=0)[None, :] * b_hat[:, None]

def magnetic_sld(qx, qy, up_angle, up_phi, rho, rho_m):
    """
    Compute the complex sld for the magnetic spin states.
    Returns effective rho for spin states [dd, du, ud, uu].
    """
    q_norm = 1/sqrt(qx**2 + qy**2) if qx != 0. or qy != 0. else 0.
    cos_spin, sin_spin = cos(-radians(up_angle)), sin(-radians(up_angle))
    cos_phi, sin_phi = cos(radians(up_phi)), sin(radians(up_phi))
    M = rho_m
    p_hat = np.array([sin_spin*cos_phi, sin_spin*sin_phi, cos_spin])
    q_hat = np.array([qx, qy, 0.])*q_norm
    M_perp = orth(M, q_hat)
    M_perpP = orth(M_perp, p_hat)
    M_perpP_perpQ = orth(M_perpP, q_hat)
    perpx = np.dot(p_hat, M_perp)
    perpy = np.sqrt(np.sum(M_perpP_perpQ**2, axis=0))
    perpz = np.dot(q_hat, M_perpP)
    return [
        rho - perpx,       # dd
        perpy - 1j*perpz,  # du
        perpy + 1j*perpz,  # ud
        rho + perpx,       # uu
    ]

This is not the same as what you have, specifically regarding the q norm.

[wpotrzebowski]

sasmodels ready for testing on OSX

[dehoni]

Thanks for the review. The sugested changes in realspace.py are done and the odd pieces resolved. It seems to reproduce nicely the "theory" for a number of known but also odd cases. I think that is the definition of magic: Me caveman has seen a burning stick.

[butlerpd]

LOL ... lots of magic in sasmodels 😄 very cool that it works!

[butlerpd]

As agreed at the fortnightly meeting of March 16, this PR is ready to merge pending verification of the build "failures." Assuming they are the usual docs building issues this should be merged

[butlerpd]

The failures seem to stem from the model unit tests failing. In particular cylinders and ellipsoids etc. This suggest that the changes to kernel and/or model info have broken at least the nuclear scattering in some way? These need to be fixed before this can be merged.

[pkienzle]

The failing tests are all of the form P@S. That's because product.py needs the following:

diff --git a/sasmodels/product.py b/sasmodels/product.py
index 704c148..2f5b34c 100644
--- a/sasmodels/product.py
+++ b/sasmodels/product.py
@@ -418,7 +418,7 @@ class ProductKernel(Kernel):
         # they will only be for form factor P, not structure factor S.
         first_mag = last_s + have_beta_mode + have_er_mode
         mag_pars = 3*p_info.parameters.nmagnetic
-        last_mag = first_mag + (mag_pars + 3 if mag_pars else 0)
+        last_mag = first_mag + (mag_pars + 4 if mag_pars else 0)
         self._magentic_slice = slice(first_mag, last_mag)
 
     def Iq(self, call_details, values, cutoff, magnetic):

[pkienzle]

Some trivial code style comments.

So other than the fix to product.py it is okay to merge.

I'm very surprised that moving things into vector functions doesn't slow down the code, at least with the compilers that I'm using. Maybe it is worth testing in a couple more operating environments (windows, nvidia) to be sure.

[dehoni]

Thanks for the suggestions. The code has been changed to reflect the coding style. I see the first successful tests coming in.

[pkienzle]

BTW, lots of lines in the doc strings are still longer than 80 characters in magnetism.rst. Since you are in there anyway you may as well shorten them.

[pkienzle]

The latest spurious doc build warning is:

2021-03-18T07:17:15.8980166Z /home/runner/work/sasmodels/sasmodels/sasmodels/modelinfo.py:docstring of sasmodels.modelinfo.ModelInfo.tests:: WARNING: py:class reference target not found: TestCondition

I don't understand why this happens. Worse, I can't reproduce on my local machine using (cd doc && make html); maybe a sphinx version issue?

Anyway, I suppress it by adding a line to doc/conf.py:

nitpick_ignore = [
    ...
    ('py:class', 'SourceModule'),
    ('py:class', 'TestCondition'),
]

[pkienzle]

Note: as painful as these spurious warnings may be, I appreciate the fact that doc build failures are no longer silently ignored. Thanks, @llimeht!

[pkienzle]

After the last doc nits are picked can this be merged?

Or are there changes needed to sasview first?

[butlerpd]

Yes. According to the discussion on this PR at the last fortnightly meeting, this can be merged as soon as it is ready. The claim is that it is independent of the sasview PR in the sense that the GUI won't have access necessarily to everything but it should not break things and will be available via scripting. The next step would be to finalize the GUI PR. Of course if that turns out to be wrong we will have to revisit but this should not affect the release of 5.0.4 which is being built off of an earlier tagged version of sasmodels

[dehoni]

Sasview does not need changes to pick up the new magnetism definition in sasmodels on the fitpage. The magnetic parameters are added with sasmodels/modelinfo.py.

[pkienzle]

I meant the '...' to represent what was already there. So:

nitpick_ignore = [
    ('py:class', 'argparse.Namespace'),
    ('py:class', 'bumps.parameter.Parameter'),
    ('py:class', 'collections.OrderedDict'),
    ('py:class', 'cuda.Context'),
    ('py:class', 'cuda.Function'),
    ('py:class', 'np.dtype'),
    ('py:class', 'numpy.dtype'),
    ('py:class', 'np.ndarray'),
    ('py:class', 'numpy.ndarray'),
    ('py:class', 'pyopencl.Program'),
    ('py:class', 'pyopencl._cl.Context'),
    ('py:class', 'pyopencl._cl.CommandQueue'),
    ('py:class', 'pyopencl._cl.Device'),
    ('py:class', 'pyopencl._cl.Kernel'),
    ('py:class', 'QWebView'),
    ('py:class', 'unittest.suite.TestSuite'),
    ('py:class', 'wx.Frame'),
    # autodoc and namedtuple is completely broken
    ('py:class', 'integer -- return number of occurrences of value'),
    ('py:class', 'integer -- return first index of value.'),
    # autodoc doesn't handle these type definitions
    ('py:class', 'Data'),
    ('py:class', 'Data1D'),
    ('py:class', 'Data2D'),
    ('py:class', 'Kernel'),
    ('py:class', 'ModelInfo'),
    ('py:class', 'module'),
    ('py:class', 'SesansData'),
    ('py:class', 'SourceModule'),
    ('py:class', 'TestCondition'),
]

Sorry for the confusion.

[butlerpd]

finally clean! I note that the only travis error now is the python 2.7 build. Is this because of 4.x? Should we drop the travis check of 2.7 now?

[pkienzle]

It seems we can no longer use anaconda to set up a python environment for 2.7.

I suggest we stop using travisCI now, but that will be a separate PR (ticket #446).

[pkienzle]

BTW, 4.x runs on python 3.