This example demonstates running a simple optimization with OPTESIM on a single machine. Below is the source file Demo1.m, followed by explanations:
0001 global global_opt_model global_spectra ... 0002 global_bvecs global_ese_routine 0003 0004 % Load Experimental Data 0005 filter_options.filter_range = 40:400; 0006 filter_options.data_range = 43:400; 0007 filter_options.var_range = 300:400; 0008 spectrum1 = filterRaw( 'WT304ns.mat', ... 0009 0.3090, filter_options ); 0010 0011 filter_options.filter_range = 20:400; 0012 filter_options.data_range = 13:400; 0013 spectrum2 = filterRaw( 'WT152ns.mat', ... 0014 0.3090, filter_options ); 0015 0016 global_spectra = {spectrum1, spectrum2}; 0017 0018 % Setup System Model 0019 optmodel = eseOptModel({'N14'},'1'); 0020 0021 optmodel = registerOptModelVariable(optmodel, ... 0022 'eeqQ', 1, 1, -10, 10, ... 0023 'eta', 1, 2, 0, 1); 0024 0025 optmodel = registerOptModelVariable(optmodel, ... 0026 'Aiso', 1, 3, -10, 10, ... 0027 'r', 1, 4, 0, 5); 0028 0029 optmodel = registerOptModelVariable(optmodel, ... 0030 'QEuler1', 1, 5, -360, 360, ... 0031 'QEuler2', 1, 6, -360, 360, ... 0032 'QEuler3', 1, 7, -360, 360); 0033 0034 global_opt_model = optmodel; 0035 0036 % Setup poweder average method 0037 % and field vectors for averaging 0038 global_ese_routine = @powderAverager; 0039 global_bvecs = gensphrser(4); 0040 0041 % Parameter Optimization 0042 x = [3.020, 0.598, 0.928, 3.415,20,20,20]; 0043 options = optimset('OutputFcn', @esePlot, 'Display','iter'); 0044 [ x_new, fval ] = fminsearch(@eseOptObjective, x, options); 0045 0046 % Estimate optimal parameter uncertainties 0047 bounds = eseOptCftInterval(x_new, fval, 0.99); 0048 0049 % Generate a report 0050 eseReport(x_new,fval,bounds);
Declare global variables. Those variables are required for OPTESIM
Load experimetal spectra. You may need to modify "filterRaw.m" to accommodate your own data structure. Please refer to the modification notes. Here we load two spectra of the same system under different experimental conditions. Line 0016 indicates that we are doing a global optimization of the two spectra simultaneously. filter_options provides parameters for processing the raw spectra, such as ranges for truncation, decay removal, and experimental noise calculation.
Setup an model for simulation. Here we define a single node model with one N14. The first parameter in eseOptModel.m is an cell array of nuclei, the second one is the combination rule of those nuclei. For example, for a system of one H2 and one N14 with a product rule, you can use:
optmodel = eseOptModel({'H2', 'N14'},'1*2');
Or, for a system with one H2 and another isotopic site of H2 or H1, you can use:
optmodel = eseOptModel({'H2', 'H2', 'H1'},'1*(2+3)');
Register varaibles for optimization. "registerOptModel"Variable function changes the "optmodel" by flagging varaibles to be subejcted to optimization. The flagging needs arguments in groups of five. The first one is the nuclear parameter. The second one tells which the nucleus this nuclear parameter belongs to. The third argument tells where to look for the value in the x vector we are providing in line 0042. You may use the same index for two or more different nuclear parameters, in that way, those nuclear parameters will be kept the same during optimization. The fourth and fiveth arguments are the lower and upper bound of this paramenter. In line 0034, we assign "optmodel" the global varaible "global_opt_model", so that it will be avaible to OPTESIM. Please note that you may change the proportions (weights) of isotopic nodes, too. In case you want to change a nuclear parameter, but do not wish to optimize on it, you may use the function "changeOptModelValue".
Indicate how we want to do the powder averaging. The options for "global_ese_routine" include "@powderAverager" for local computation, "@powderAveragerDistributed" for distributed computation, and "@powderAveragerDistributedNGP" for distributed computation with non-geometry preserving model. Line 0039 setups field orientations for powder averaging. The function "gensphrser" generates vectors distributed on a unit sphere. For crystal spectrum, you can specify only one orientation to global_bvecs.
Setup up initial value for optimization. Here we use MATLAB built-in fminsearch for the optimization. For practical use, you may want to repeat "fminsearch" serveral times to avoid local convergence.
Find out the uncertainties of the optimal parameters and generate a report.