iFit: McStas/McXtrace simulation/optimization

McStas [www.mcstas.org]

  1. Building an instrument model
    1. Specifying which monitor to use
    2. Requirements
  2. Running a single simulation
  3. Getting instrument information and Displaying the instrument geometry
  4. Running series of simulations (scans)
  5. Optimizing instrument parameters
  6. 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

Building an instrument model

NOTE: the 'mcstas' routine was used up to iFit version 1.9.
It will be obsoleted in future versions, as the 'mccode' solution below is more robust.
You may still used it, and get help for it with: '>> help mcstas'.

The syntax for building a McCode (McStas or McXtrace) instrument model is:

>> model = mccode('instrument');	% uses the given instrument file (*.instr)
>> model = mccode('gui'); % lists all available instruments to select one
>> model = mccode(''); % pop-up a file selector to select an *.instr file
>> model = mccode('defaults'); % uses the templateDIFF instrument as example
which compile the instrument, and create an iFunc model.
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.

In addition, it is possible to select an instrument from the McStas installation with:

You may even assemble these models as iFunc arrays to iteratively perform series of simulations, or use arithmetic operations to combine models (see the iFunc documentation).

It is also possible to finely tune the instrument configuration using standard McStas/McXtrace options with the syntax:
>>  mccode('instrument', options)
and the following options:
options.dir directory where to store results, or set automatically (string) ;
the last simulation files are stored therein 'sim'.
options.ncount number of neutron events per iteration, e.g. 1e6 (double).
number of processors/cores to use with MPI on localhost or cluster (integer) ;
when MPI is available, and mpi options is not given, all cores are then used. Use options.mpi=1 for a serial calculation.
options.machines file name containing the list of machines/nodes to use (string), as often found in cluster job schedulers.
options.seed random number seed to use for each iteration (double)
options.gravitation 0 or 1 to set gravitation handling in neutron propagation (boolean)
options.monitor a single monitor name pattern to read, or left empty for the last (string). The whole monitors are available as iData objects in the model.UserData.monitors field. Wildcard expression can be used. See below for more advanced customisation.
options.mpirun set the executable path to 'mpirun', or automatically found when not set.
options.compile 0 or 1 to force re-compilation of the executable. Try that first if you can not create/use the object (an old executable may be used).
The options can be given as a structure or as a string, e.g. 'ncount=1e6; compile=1'.

NOTE: if the instrument can not be built, force the compilation with:
>>  mccode('instrument', 'compile=1')
as a previous executable may be in the way, and block the generation of the new model.

The scalar instrument parameters of type 'double' are used as model parameters.
Other parameters (e.g. of type string and int) are stored as a structure in:
The options ncount, seed, gravitation, monitor can be changed after the model creation in the UserData area of the object e.g.:

Specifying which monitor to use

The default behaviour is to leave options.monitor left unset at the model creation, so that all monitors are read when simulating instruments. The last monitor is then used as model value. This provides maximum information, but the processing overhead time is increased due to more data handling compared to a single/restricted monitor name selection.
In order to improve the computational efficiency, it is advised to specify which files have to be imported to define the model value by giving a file name or pattern (wildcard patterns can be used). This specification can be done either when building the model:
or afterwards:
which in both cases will select monitor files which match the given string as an exact file name or match the pattern when using a wildcard expression ( '*Theta*' ). In case the selection brings multiple files, the last one is used. When the monitor name selection does not match any file, or is given as empty, all monitors are imported, and the last one in the list is used.

All selected monitors are imported as iData object(s), in:
In practice, you will often either not specify the 'monitor' at the model creation to leave the system select the last one, or indicate a specific monitor name/pattern of your choice.

Customize the model value
In addition, you may customize how to use the monitor data either by adding expression after the initial filename/pattern, or by extending the model after its creation.

The monitor to use as optimization criteria (options.monitor) can be part of a more complex expression, which still must start with the file name/pattern to select the monitor(s) and then refer to itself with the 'signal' variable (which are iData object(s)), such as in the following example where we sum all selected monitors from the pattern (when more than one are found), and define the criteria as Amplitude/width˛ (for a single peak):
>> options.monitor='*Banana_Theta*; signal=plus(signal); signal=max(signal)/std(signal)^2;'
>> model=mccode('templateDIFF', options)
The options.monitor can be given in a compact way if the file pattern and the expression are enclosed in single or double quotes, such as in:
An other possibility (less recommended) is to extend the model after its creation:


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) or McXtrace. When used with the standalone version, no Matlab license is needed.

Running a single simulation

The model value is corresponds with the 'options.monitor' selected, i.e. the last one, or a specific name pattern. This selection can be changed any time without rebuilding the model (see above).

The instrument can be evaluated with given parameters and optionally axes:
>> results = feval(model);			% return the monitor, using current/default parameters
>> results = feval(model, parameters); % return the monitor as an array using guessed axes
>> results = feval(model, parameters, x,y,...); % return the monitor as an array, using specified axes
>> results = feval(model, parameters, nan); % return the raw monitor as an array
but it may be more convenient to manipulate the model value as iData objects:
>> result = iData(model, parameters);		% return monitor interpolated onto guessed axes
>> result = iData(model, parameters, nan); % return raw monitor as an iData object
>> result = iData(model, parameters, x,y,...); % return monitor interpolated onto axes x,y,...
>> plot(result);
>> plot(model.UserData.monitors)
where parameters is a structure, or string, or vector that holds parameter names and values, e.g.
or in a more compact way, giving parameters as a single string with members assigned with '=' and separated with the ';' character. An empty parameter value (e.g. []) will compute the model evaluation using current/default parameters.

In the following example, we plot the model value (raw monitor), and then all the available monitors in model.UserData.monitors
value =  iData (methods,doc,plot,plot(log),values,more...) 2D object:
    [Tag] [Dimension]                                     [Title] [Last command]          [Label/DisplayName]
  iD61442    [51 170]   'templateDIFF.instr McCode [... "signal"' iD61442=setalias(iD6... templateDI/templateDI...
model.UserData.monitors =  array [5  1] iData (methods,doc,plot,plot(log),values,more...) object:

Index     [Tag] [Dimension]                                     [Title] [Last command]          [Label/DisplayName]
    1   iD61449     [40 40] 'mccode.sim McCode sim file ... "Data ..."' iD61449=load(iData,'... Diff_Mono_XY
    2   iD61457      [40 1] 'mccode.sim McCode sim file ... "Inten..."' iD61457=load(iData,'... Diff_Mono_Lambda
    3   iD61459      [40 1] 'mccode.sim McCode sim file ... "Inten..."' iD61459=load(iData,'... Diff_Sample_Lambda
    4   iD61461     [340 1] 'mccode.sim McCode sim file ... "Inten..."' iD61461=load(iData,'... Diff_BananaTheta
    5   iD61463    [25 170] 'mccode.sim McCode sim file ... "Data ..."' iD61463=load(iData,'... Diff_BananaPSD
You can equally set the non numerical scalar parameters (e.g. strings and integers):
The unspecified instrument parameters are used with their default values (when defined).

The model value, and the raw monitors in model.UserData.monitors, in the form of iData objects have an additional Parameters field which holds the instrument parameters used for the simulation.
>> value.Parameters
>> get(model.UserData.monitors(1),'Parameters')
Usual data analysis (display, mathematics, fit, export) can be performed (refer to the iFit and iData documentation).

Getting instrument information and Displaying the instrument geometry

To display the instrument geometry, use the syntax:
>> mccode_display(model);
which open a geometrical representation of the instrument. The current model parameters are used. It can be given as well a set of parameter values:
>> mccode_display(model, parameters);
An empty parameter argument requests to use the current value. In addition, this routine returns the list of components, and can be given more options (3rd argument), as 'html','png', and 'x3d'.
>> mccode_display(model, [], 'html png x3d');
to generate hardcopy files and geometry rendering. The generated filenames are then displayed.

In order to get information about the instrument, especially its input parameters, you can access the model.UserData area:
>> model.UserData

ans =

instrument_source: [1x7011 char]
instrument_executable: [1299344x1 double]
instrument_exe: 'templateDIFF.out'
instrument_info: [1x960 char]
options: [1x1 struct]
Parameters: [1x1 struct]
Parameters_Constant: [1x1 struct]
monitors: [5x1 iData]
and especially the instrument_info and instrument_source, which you can edit with e.g:
or even in a separate window:
>> TextEdit(model.UserData.instrument_source)
>> TextEdit(model.UserData.instrument_info)

Running series of simulations (scans)

You may scan the model parameters when they are assigned as vector.
A more compact way can be used, using explicit scan values and axes:
The argument returned is an iData object array containing the model monitor as a function of the scanned parameters. You may compute the integral value with the 'sum' method:
>> integral = sum(scan, 0);
which 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 scans, use slice-o-matic for volume inspection
You may as well concatenate the scan steps into a higher dimensionality volume using e.g. the iData/cat method:
>> catenated = cat(scan, 2);
iFit: MsCtas          scan: ploting the integral monitor value
Multidimensional scans are also possible. In the following example, we launch a quick 2D scan with parameters specified as a structure:
>> model = mccode('templateDIFF','monitor=Diff_BananaTheta');
>> scan = iData(model, struct('RV',[0.7 0.9 1.2],'L2',[1 1.3 1.5]), nan);
Then we may plot, as a 2D surface, the integral value for the model value monitor:
>> plot(iData(sum(scan,0)));
iFit: McStas:
          evolution of a monitorYou may look at the instrument parameters used for the scans, e.g.
>> get(scan,'Parameters.RV')
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.

You may also scan the non numerical parameters with e.g. for the 'powder' string parameter:

Optimizing instrument parameters

The optimization of instrument parameters is performed as simply as a single simulation, using the fmax method to maximize the monitor value (you may similarly use fmin to minimize the model value).
>> best_parameters = fmax(model);		% maximize with the current/default parameters
>> best_parameters = fmax(model, parameters); % maximize with given initial parameters
The parameters, as a vector, string or named structure, are initiated to their starting value, from which the optimization proceeds. An empty value requests to use the current/default values. 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.

You can tune the optimizer configuration, as explained in the iFunc documentation, with syntax:
>> best_parameters = fmax(model, parameters, options);
as well as any other optimizer configuration fields such as
The options can be given as a structure or assembled into a string. The default optimizer configuration is used when not specified (see the Optimizers page for more details). The optimizer routine is automatically guessed, but can be forced (see example below):
Eventhough shown as 'fminXXX' methods, the iFunc/fmax method in use here will make sure to maximize the monitor value by using the specified optimiser.

It is highly recommended to request the opimisation using the raw monitors. As indicated above, the axes are specified as additional arguments after the options:
>> best_parameters = fmax(model, parameters, options, x, y, ...);
and in particular you should probably use:
>> best_parameters = fmax(model, parameters, options, nan);	% use raw monitors
The optimization returns the best parameter guess, as well as the monitor value for this solution.

In the following example, a McStas model is assembled, specifying which monitor to use as model value. Then all parameters are fixed except RV and L2 which are bound (see section below). Last, the optimizer is explicitly indicated (else an automatic guess is made) as well as a request for a plot during the optimisation. In addition, similarly with other iFit optimization methods, the optimizer final status and optional returned information can be obtained with e.g. 
>> model = mccode('templateDIFF','monitor=*Theta*');	% monitor contains 'Theta'
>> fix(model, 'all'); % fix all but RV L2
>> model.RV='free'; model.RV=[0.5 0.8 1.5]; % set RV bounds and value in a single assignment
>> model.L2='free'; model.L2=[1 3 4];
>> [parameters, fval, status, output]=fmax(model,[], 'optimizer=fminimfil; OutputFcn=fminplot', nan);
where output stores most of the optimization process data, including uncertainty on the best parameter values, output.parsHistoryUncertainty and output.parsHessianUncertainty.

To abort the optimization, just close the optimization plot window (or click the STOP button). The actual end occurs after a few more iterations.

You can then plot the monitor(s) for the optimised parameters with:  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='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 parameter search to specified ranges, and fix some of the parameters. As the model is an iFunc object, you may fix some parameter with any of:
or alternatively clear (free) them with:
You may as well bound a parameter with:iFit: McStas wrapper:
        dynamic optimization evolution or even bound a parameter and set its value in a single line:
In the following example we optimise the monochromator radius of curvature in the templateDIF diffractometer, in order to get the maximum flux in the diffractogram:
>> model = mccode('templateDIFF','monitor=*Theta*');	% monitor contains 'Theta'
>> fix(model, 'all'); % fix all
>> model.RV='free'; model.RV=[0.5 1 1.5]; % set RV free, bounds and value in a single assignment
>> [parameters, fval, status, output]=fmax(model,[], ...
'optimizer=fminpso; OutputFcn=fminplot;TolFun =5%;TolX=5%;ncount=1e5;MaxFunEvals=100', nan);
You can then plot the monitor(s) for the optimised parameters with: 
We strongly recommend to limit the parameter search within boundaries during optimizations.

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