Merge pull request fitbenchmarking#1405 from fitbenchmarking/add_spin…

…w_1d_slices

Added support for SpinW 1D cuts of 2D powder data

docs/source/users/BenchmarkProblems.rst

@@ -372,3 +372,15 @@ conditions for the parameters.
 The data for testing the `Gaussian` fitting can be found within
 the `gaussian` subfolder. It contains a dataset with 16 different starting 
 conditions for the parameters. 
+
+SpinW 2D Powder Data
+====================
+
+**Download** :download:`.zip <https://fitbenchmarking.github.io/assets/datasets/spinw_powder_data.zip>`
+or :download:`.tar.gz <https://fitbenchmarking.github.io/assets/datasets/spinw_powder_data.tar.gz>`
+
+This problem (also found in the folder `examples/benchmark_problems/SpinW_powder_data`)
+contains 1D cuts of 2D powder data simulated using SpinW, using the approach outlined
+in `this tutorial <https://spinw.org/tutorials/39tutorial>`_ .
+
+This problem has 8 unknown parameters and 186 data points.

docs/source/users/install_instructions/externals.rst

@@ -53,6 +53,13 @@ website <https://pace-neutrons.github.io/Horace/3.6.0/Download_and_setup.html>`_
 In addition, MATLAB and the MATLAB engine must be installed following the
 :ref:`instructions given below<matlab-install>`.
 
+SpinW
+-----
+
+SpinW can be installed by following the instructions `on the SpinW website
+<https://spinw.org/IntroToSpinW/#/install1>`__. In addition, MATLAB and the MATLAB
+engine must be installed following the :ref:`instructions given below<matlab-install>`.
+
 .. _levmar-install:
 
 Levmar

docs/source/users/problem_definition_files/horace.rst

@@ -15,6 +15,8 @@ Examples of horace problems are:
 
 .. literalinclude:: ../../../../examples/benchmark_problems/Horace/PCSMO_at_001_data.txt
 
+.. literalinclude:: ..../../../examples/benchmark_problems/SpinW_powder_data/tri_AFM_powder.txt
+
 .. note::
   The Horace file format requires you to have run the benchmark problem in Horace 
   using :code:`fit()` and :code:`simulate()` successfully. Relevant links on 
@@ -166,6 +168,14 @@ Examples of the simulate_function:
 
 .. literalinclude:: ../../../../examples/benchmark_problems/Horace/m_scripts/simulate_functions/fb_simulate_pcsmo_test.m
 
+plot_type
+  For SpinW 2D powder fitting problems, the plot type must be specified. Currently, support is only available for plotting
+  1D cuts and so plot_type should be set to '1D_cuts'.
+
+q_cens
+  For SpinW 2D powder fitting problems, the values of Q at which the 1D cuts have been taken.
+  These values should be provided as a comma-separated string.
+
 .. note::
    All the functions needed in the fitting must be in the subdirectory of the benchmark problem.
 

examples/benchmark_problems/SpinW_powder_data/m_scripts/functions/tri_AFM_powder.m

@@ -0,0 +1,3 @@
+function y = tri_AFM_powder(x, p)
+
+    y = @(x, p) matparser(x, 'param', p, 'mat', {'J_1'}, 'init', true);

examples/benchmark_problems/SpinW_powder_data/m_scripts/simulate_functions/fb_simulate_tri_AFM.m

@@ -0,0 +1,43 @@
+function y = fb_simulate_IX_1D_test1(w,fitpars,msk)
+% simulate loop to solve for the parameters 
+
+fit_func = @tri_AFM_powder;
+model_param = reshape(fitpars(1), 1, []);
+
+Ei = 20; % meV
+% Resolution for MARI with Gd chopper running at 200Hz and Ei=20meV (from PyChop)
+eres = @(en) 2.1750e-04*sqrt((Ei-en).^3 .* ( (32.27217*(0.168+0.400*(Ei./(Ei-en)).^1.5)).^2 + (14.53577*(1.168+0.400*(Ei./(Ei-en)).^1.5)).^2) );
+% Q-resolution (parameters for MARI)
+e2k = @(en) sqrt( en .* (2*1.67492728e-27*1.60217653e-22) )./1.05457168e-34./1e10;
+L1 = 11.8;  % Moderator to Sample distance in m
+L2 = 4.0;   % Sample to detector distance in m
+ws = 0.05;  % Width of sample in m
+wm = 0.12;  % Width of moderator in m
+wd = 0.025; % Width of detector in m
+ki = e2k(Ei);
+a1 = ws/L1; % Angular width of sample seen from moderator
+a2 = wm/L1; % Angular width of moderator seen from sample
+a3 = wd/L2; % Angular width of detector seen from sample
+a4 = ws/L2; % Angular width of sample seen from detector
+dQ = 2.35 * sqrt( (ki*a1)^2/12 + (ki*a2)^2/12 + (ki*a3)^2/12 + (ki*a4)^2/12 );
+
+tri = sw_model('triAF',1);
+fitpow = sw_fitpowder(tri, w, fit_func, model_param, 'independent');
+
+% set parameters passed to powspec
+fitpow.powspec_args.dE = eres; % energy resolution
+fitpow.powspec_args.fastmode = true;
+fitpow.powspec_args.neutron_output = true;
+fitpow.powspec_args.nRand = 2e2; % low for speed (typically want > 1e3)
+fitpow.powspec_args.hermit = true;
+% set parameters passsed to sw_instrument
+fitpow.sw_instrument_args = struct('dQ', dQ, 'ThetaMin', 3.5, 'Ei', Ei);
+
+fitpow.powspec_args.dE = eres;
+fitpow.powspec_args.fastmode = true;
+fitpow.powspec_args.neutron_output = true;
+
+[y, bg] = fitpow.calc_spinwave_spec(fitpars)
+y = y'
+
+end

examples/benchmark_problems/SpinW_powder_data/m_scripts/wye_functions/fb_wye_tri_AFM.m

@@ -0,0 +1,9 @@
+function [w, y, e, msk] = fb_wye_tri_AFM(datafile, path)
+% Gets w , x, y ,e and msk from the object
+
+w = load(datafile).w;
+y = [w(1).y; w(2).y; w(3).y];
+e = [w(1).e; w(2).e; w(3).e];
+msk = []
+
+end

examples/benchmark_problems/SpinW_powder_data/tri_AFM_powder.txt

@@ -0,0 +1,10 @@
+# FitBenchmark Problem
+software = 'Horace'
+name = 'Triangular AFM Powder Spectrum'
+description = 'Simulated powder spectrum of a triangular AFM with constant background and normally distributed noise, fit to 1D cuts rather than 2D data'
+input_file = 'tri_AFM_powder.mat'
+function = 'foreground=m_scripts/functions/tri_AFM_powder.m,J1=0.85,b1=0,b2=1.4,b3=0,b4=1.4,b5=0,b6=1.4,b7=29.2'
+wye_function = 'matlab_script=m_scripts/wye_functions/fb_wye_tri_AFM.m'
+simulate_function = 'matlab_script=m_scripts/simulate_functions/fb_simulate_tri_AFM.m'
+plot_type = '1D_cuts'
+q_cens = '0.8,1.2,1.6'

fitbenchmarking/parsing/horace_parser.py

@@ -23,9 +23,9 @@ def horace_on():
         horace_location = os.environ["HORACE_LOCATION"]
         spinw_location = os.environ["SPINW_LOCATION"]
         eng.evalc("restoredefaultpath")
-        eng.evalc(f"addpath('{horace_location}')")
+        eng.evalc(f"addpath(genpath('{horace_location}'))")
         eng.evalc("horace_on")
-        eng.evalc(f"addpath('{spinw_location}')")
+        eng.evalc(f"addpath(genpath('{spinw_location}'))")
     elif (
         "HORACE_LOCATION" not in os.environ and "SPINW_LOCATION" in os.environ
     ):
@@ -181,7 +181,6 @@ def _create_function(self) -> typing.Callable:
 
         def fit_function(x, *p):
             # Assume, for efficiency, matching shape => matching values
-            # print(*p)
             if x.shape != self._horace_x.shape:
                 return np.ones(x.shape)
             eng.evalc(f"global {self._horace_msk}")
@@ -232,4 +231,41 @@ def set_persistent_vars(path):
             print(eng.evalc(f"load('{path}')"))
 
         self.fitting_problem.set_persistent_vars = set_persistent_vars
+
+        plot_type_options = ["1d_cuts"]
+
+        if "plot_type" in self._entries:
+            if self._entries["plot_type"].lower() in plot_type_options:
+                self.fitting_problem.additional_info["plot_type"] = (
+                    self._entries["plot_type"].lower()
+                )
+            else:
+                raise ParsingError(
+                    "The plot type should be one of these "
+                    f"options {plot_type_options}"
+                )
+
+        if self.fitting_problem.additional_info["plot_type"] == "1d_cuts":
+            eng.evalc(f"ebin_cens = {self._horace_w}(1).x")
+            self.fitting_problem.additional_info["ebin_cens"] = np.array(
+                eng.workspace["ebin_cens"], dtype=np.float64
+            )[0]
+            if "q_cens" in self._entries:
+                self.fitting_problem.additional_info["q_cens"] = self._entries[
+                    "q_cens"
+                ].split(",")
+                self.fitting_problem.additional_info["n_cuts"] = len(
+                    self.fitting_problem.additional_info["q_cens"]
+                )
+            else:
+                raise ParsingError("q_cens are required for plotting 1D cuts")
+            if not float(
+                len(self.fitting_problem.data_y)
+                / self.fitting_problem.additional_info["n_cuts"]
+            ).is_integer():
+                raise ParsingError(
+                    "Number of data points must be divisible "
+                    "by number of q_cens"
+                )
+
         return super()._set_additional_info()