iFit: McStas simulation/optimization

McStas [www.mcstas.org]

  1. Running a single simulation
  2. Getting instrument information and Displaying the instrument geometry
  3. Running series of simulations (scans)
  4. Optimizing instrument parameters
  5. Optimizing instrument parameters with constraints

Commands we use in this page: mcstas (mcxtrace)

In this document, we demonstrate how iFit/Optimizers and iData can all be used transparently to simulate and optimize McStas instrument models.

To import a McStas simulation results, just use
>> results = iData('single_detector_file')
>> results = iData('directory') % import everything in the directory, including sim files (files may imported more than once)
>> results = iData('scan_directory/mcstas.dat') % for scans
>> results = iData('scan_directory/mcstas.sim') % for all files from the directory, uniquely.
which returns a single or an array of iData object(s), that you can plot with e.g. subplot.
The resulting McStas data has an additional Parameters field which holds the instrument parameters used for the simulation.
>> results(1).Parameters

Running a single simulation

As iFit contains powerful import routines, that support McStas data files, it is straight forward to launch a McStas instrument simulation from Matlab and get back the simulation results. To do so, you will obviously need to have iFit installed (refer to the Install page), as well as McStas (refer to the McStas web site). An easy solution is to install our ILL/CS Live DVD ready-to-run system with McStas pre-installed, and then put on top Matlab (requires a valid license).

The syntax for running a McStas instrument model is:
>> results = mcstas('instrument', parameters, options);
where parameters is a structure that holds parameter names and values, e.g.

parameters.RV=1
parameters.lambda=2.36
parameters.powder='Na2Ca3Al2F14.laz'

and both numerical values and strings are supported. Numerical values should be single scalars, whereas strings are used as-is, that is they may assign DEFINITION instrument parameters such as file names, arrays and other types that are sent to the commands line. The unspecified instrument parameters are used with their default values (when defined).

The parameters can also be given as a single string with members separated with the ';' character ans assigned with '=':
>> results = mcstas('instrument', 'RV=1; lambda=2.36; powder=Na2Ca3Al2F14.laz');
When the instrument is not found in the current directory, it is searched in the McStas library (may be specified with the MCSTAS and MCXTRACE environment variables), and eventually copied locally, which triggers its compilation.

If McStas/McXtrace executables (e.g. 'mcrun') are not found, you may have to extend the PATH with e.g.:
The last argument options controls simulation configuration.
In the following examples, we define required structures on the fly, run the templateDIFF instrument (the file must be accessible from Matlab current path) and get back monitors:
>> monitors=mcstas('templateDIFF', struct('RV',1), struct('ncount',1e6))
>> monitors=mcstas('templateDIFF', 'RV=1', 'ncount=1e6'); % same as above
monitors = array [1 5] iData object:

Index [Tag] [Dimension] [Title] [Last command] [Label]
1 id445199 [40 40] 'File Diff_Mono_XY_1299773321...' id445199=load(iData,... Diff_Mono_XY
2 id445212 [40 1] 'File Diff_Mono_Lambda_129977...' id445212=load(iData,... Diff_Mono_...
3 id445224 [40 1] 'File Diff_Sample_Lambda_1299...' id445224=load(iData,... Diff_Sampl...
4 id445236 [340 1] 'File Diff_BananaTheta_129977...' id445236=load(iData,... Diff_Banan...
5 id445248 [170 25] 'File Diff_BananaPSD_12997733...' id445248=load(iData,... Diff_Banan...
>> subplot(monitors)
The last line displays the simulation results. Usual data analysis (display, mathematics, fit, export) can be performed (refer to the iFit and iData documentation).
The resulting McStas data has an additional Parameters field which holds the instrument parameters used for the simulation.
>> monitors(1).Parameters
The 'monitors' option can be specified to limit the number of monitor files to load after each simulation. Specifying a monitor file name list highly speeds-up the importation of files after each iteration.
It should contain one token string or a cellstr {...} of tokens to search for. These tokens are compared with the monitor file names, and only matches are imported. If no match can be identified, all monitors are read, the monitor list can be reduced by searching tokens in the data sets contents.

Getting instrument information and Displaying the instrument geometry

To display the instrument geometry, use the syntax:
>> results = mcstas('templateDIFF', 'RV=1','mode=display');
which will launch 'mcdisplay' in the background, and open the resulting figure. No neutron trajectory is shown though, only the instrument geometry.

In order to just get information about the instrument, especially its input parameters, you can run
>> results = mcstas('instrument', 'parameters=values...','mode=info');
>> result = mcstas('instrument', '--info'); % same as above

Running series of simulations (scans)

It is very simple to perform parameter scans, even in multi-dimensional spaces. For this, the parameter values just need to be vectors of numeric or vector cells, such as
Scans with string parameters are also possible, such as with
and then executes for instance:
>> [integral,monitors]=mcstas('templateDIFF', parameters, struct('ncount',1e4));
A live display of the on-going scan is shown when specifying OutputFcn='fminplot' in the last options argument (struct).

The integral argument returned is an iData object containing the integral of monitors as a function of the scanned parameters. It can be directly plotted with e.g.
>> plot(integral);     % simple plot, or median surface for multidimensional scans
>> scatter3(integral); % coloured points (possibly in 2D,3D)
>> plot3(integral); % coloured lines in 2D, or volume for multidimensional scans
>> slice(integral); % only for 3D, use slice-o-matic for volume inspection
iFit: MsCtas
          scan: ploting the integral monitor valueThe monitors argument returned is an array containing individual scan monitors. You may extract and catenate a selection (slices) or it by e.g. the iData/cat method.

In the following example, we launch a quick 2D scan:
>> [p,m]=mcstas('templateDIFF', struct('RV',[0.7 0.9 1.2],'L2',[1 1.3 1.5]),struct('ncount',1e4));
Then we may plot, as a 2D surface, the integral value for the last (5th) monitor:
>> plot(p(:,:,5));
iFit: McStas:
          evolution of a monitorYou may look at the instrument parameters used for the scans, e.g.
>> get(m,'RV')
which are also available in the Parameters field of the objects.
Let's finally look at the evolution of the 4-th monitor, which is a 1D diffractogram, along an 'RV' parameter scan:
>> a = squeeze(m(2,:,:));           % extract an 'RV' line, with L2=1.3
>> b = a(:,4); % get 4th monitor on this line
>> setaxis(b,2,'RV'); % define a new axis along the line RV
>> c = cat(2, b); % assemble into a single iData object
>> plot(c,'tight interp'); % plot
This methodology also works with 2D monitors assembled along, e.g. a scan line, which creates a 3D volume to explore with slice, or plot  or plot3.

Warning: Beware of the number of iterations required to scan large dimensionality spaces. In practice you should avoid spaces above 2 or 3 scanned parameters. Exceeding this will undoubtedly require very long computation times, and large memory storage.

Optimizing instrument parameters

The optimization of instrument parameters is performed as simply as a single simulation, just by e.g. setting options.mode to 'optimize'.
>> [solution, monitors] = mcstas('instrument', parameters, options);
The parameters, as a named structure, are initiated to their starting value, from which the optimization proceeds. Only numerical parameters can be optimized. Also, it is highly encouraged to use bounded parameters (that is define constraints, see below), else the optimizer can get lost, produce un-meaningful results, or return NaN parameter values.

More specifically, the options structure may contain, in addition to the optimization mode, the following members:
as well as any other optimizer configuration fields such as
The default optimizer configuration is used when not specified (see the Optimizers page for more details), as well as the interactive plotting and verbose information during the optimization.

The optimization returns the best parameter guess, as well as the monitor value for this solution.
In addition, similarly with other iFit optimization methods, the optimizer final status and optional returned information can be obtained with e.g.
>> [parameters, monitors, status, output]=mcstas('templateDIFF',struct('RV',1), ...
struct('monitors','Banana','mode','optimize','optimizer','fminimfil'));
where output stores most of the optimization process data, including uncertainty on the best parameter values, output.parsHistoryUncertainty and output.parsHessianUncertainty.

iFit: McStas wrapper:
        dynamic optimization evolutionIn the following example, we search for an optimal monochromator curvature (parameter RV) on the templateDIFF instrument model, in order to maximize the Banana monitors (as 'monitors' option):
>> [parameters, monitors]=mcstas('templateDIFF',struct('RV',1), ...
struct('monitors','Banana','optimizer','fminimfil'));
parameters =
RV: 0.75
The monitors to use as optimization criteria (options.monitors) can be part of a more complex expression, which still must start with the file name of the monitor, and then refer to itself with the 'this' variable (which is an iData object), such as in the following example where we define the criteria as Amplitude/width˛ (for a single peak)
options.monitors='Banana; this=max(this)/std(this)^2'
options.monitors='Banana/std(this)^4'
When the monitor filter selects more than one data set (e.g. more than one file match the monitor selection name), 'this' is an iData array. Specifying a monitor list highly speeds-up the importation of files after each iteration.
For simpler definitions, the options can be given as a string, similarly to the starting parameters:
mcstas('templateDIFF','RV=[1 3]', 'monitors=Banana; mode=optimize;')
During the optimization, a dynamic plot of the criteria evolution (the sum of the Banana detector counts) and the parameters evolution is shown. The number of optimized parameters can be large, and only the first 3 ones will be shown on the dynamic plot.
The default stopping condition is met when the integral monitor change is lower than 1%.

To abort the optimization, just close the optimization plot window. The actual end occurs after a few more iterations.

Warning: as McStas uses a Monte-Carlo technique (which means 'random'), the criteria is a noisy function. The stop conditions must then usually be quite relaxed for the optimization to end, or the number of ncounts must be large enough to minimize the statistical uncertainty. TolFun='0.1%' and TolX='0.1%' are fair choices. Also, we recommend to constrain parameter ranges (see below). Also, increasing the 'ncount' (which is 1e5 by default) may help the optimizer to search in a noisy criteria landscape.

Optimizing instrument parameters with constraints

It is quite usual the restrict the parameters search to specified ranges, and to fix some of the parameters. We have seen that an optimization can be initiated by defining the starting parameter value with e.g. parameters.RV=1, together with e.g. options.mode='optimize'.

To fix a parameter value during the optimization, use a string parameter value, such as :
To restrict a parameter search to a given range, set its minimum and maximum values :
If you wish to restrict tne search range as well as define the starting parameter value, use a 3 element vector [min start max ]:
We strongly recommend to limit the parameter search within boundaries during optimizations.
>> [parameters, monitors]=mcstas('templateDIFF',struct('RV',[0.7 0.8 1.2]), ...
struct('monitors','Banana','mode','optimize','mpi',4,'compile',1));


E. Farhi - iFit/McStas - Mar. 22, 2017 1.9 - back to Main iFit Page ILL,
            Grenoble, France <www.ill.eu>