text:
function fmri_design = spm_cfg_fmri_design
% SPM Configuration file for fMRI model specification (design only)
%_______________________________________________________________________
% Copyright (C) 2005-2016 Wellcome Trust Centre for Neuroimaging

% $Id: spm_cfg_fmri_design.m 6952 2016-11-25 16:03:13Z guillaume $

% ---------------------------------------------------------------------
% dir Directory
% ---------------------------------------------------------------------
dir         = cfg_files;
dir.tag     = 'dir';
dir.name    = 'Directory';
dir.help    = {'Select a directory where the SPM.mat file containing the specified design matrix will be written.'};
dir.filter = 'dir';
dir.ufilter = '.*';
dir.num     = [1 1];
% ---------------------------------------------------------------------
% units Units for design
% ---------------------------------------------------------------------
units         = cfg_menu;
units.tag     = 'units';
units.name    = 'Units for design';
units.help    = {'The onsets of events or blocks can be specified in either scans or seconds.'};
units.labels = {
                'Scans'
                'Seconds'
}';
units.values = {
                'scans'
                'secs'
}';
% ---------------------------------------------------------------------
% RT Interscan interval
% ---------------------------------------------------------------------
RT         = cfg_entry;
RT.tag     = 'RT';
RT.name    = 'Interscan interval';
RT.help    = {'Interscan interval, TR, (specified in seconds).  This is the time between acquiring a plane of one volume and the same plane in the next volume.  It is assumed to be constant throughout.'};
RT.strtype = 'r';
RT.num     = [1 1];
% ---------------------------------------------------------------------
% fmri_t Microtime resolution
% ---------------------------------------------------------------------
fmri_t         = cfg_entry;
fmri_t.tag     = 'fmri_t';
fmri_t.name    = 'Microtime resolution';
fmri_t.help    = {
                  'The microtime resolution, t, is the number of time-bins per scan used when building regressors. '
                  'If you have performed slice-timing correction, change this parameter to match the number of slices specified there; otherwise, you would typically not need to change this.'
                  ''
}';
fmri_t.strtype = 'n';
fmri_t.num     = [1 1];
fmri_t.def     = @(val)spm_get_defaults('stats.fmri.t', val{:});
% ---------------------------------------------------------------------
% fmri_t0 Microtime onset
% ---------------------------------------------------------------------
fmri_t0         = cfg_entry;
fmri_t0.tag     = 'fmri_t0';
fmri_t0.name    = 'Microtime onset';
fmri_t0.help    = {
                   'The microtime onset, t0, is the reference time-bin at which the regressors are resampled to coincide with data acquisition.'
                   'If you have performed slice-timing correction, you must change this parameter to match the reference slice specified there.'
                   'Otherwise, you might still want to change this if you have non-interleaved acquisition and you wish to sample the regressors so that they are appropriate for a slice in a particular part of the brain.'
                   'For example, if t0 = 1, then the regressors will be appropriate for the first slice; if t0=t, then the regressors will be appropriate for the last slice.'
                   'Setting t0 = t/2 is a good compromise if you are interested in slices at the beginning and end of the acquisition, or if you have interleaved data, or if you have 3D EPI data.'
                   ''
}';
fmri_t0.strtype = 'n';
fmri_t0.num     = [1 1];
fmri_t0.def     = @(val)spm_get_defaults('stats.fmri.t0', val{:});
% ---------------------------------------------------------------------
% timing Timing parameters
% ---------------------------------------------------------------------
timing         = cfg_branch;
timing.tag     = 'timing';
timing.name    = 'Timing parameters';
timing.val     = {units RT fmri_t fmri_t0 };
timing.help    = {
                  'Specify various timing parameters needed to construct the design matrix. This includes the units of the design specification and the interscan interval.'
                  ''
                  'Also, with longs TRs you may want to shift the regressors so that they are aligned to a particular slice.  This is effected by changing the microtime resolution and onset. '
}';
% ---------------------------------------------------------------------
% nscan Number of scans
% ---------------------------------------------------------------------
nscan         = cfg_entry;
nscan.tag     = 'nscan';
nscan.name    = 'Number of scans';
nscan.help    = {'Specify the number of scans for this session.The actual scans must be specified in a separate batch job ''fMRI data specification''.'};
nscan.strtype = 'n';
nscan.num     = [1 1];
% ---------------------------------------------------------------------
% name Name
% ---------------------------------------------------------------------
name         = cfg_entry;
name.tag     = 'name';
name.name    = 'Name';
name.help    = {'Condition Name'};
name.strtype = 's';
name.num     = [1 Inf];
% ---------------------------------------------------------------------
% onset Onsets
% ---------------------------------------------------------------------
onset         = cfg_entry;
onset.tag     = 'onset';
onset.name    = 'Onsets';
onset.help    = {'Specify a vector of onset times for this condition type. '};
onset.strtype = 'r';
onset.num     = [Inf 1];
% ---------------------------------------------------------------------
% duration Durations
% ---------------------------------------------------------------------
duration         = cfg_entry;
duration.tag     = 'duration';
duration.name    = 'Durations';
duration.help    = {'Specify the event durations. Epoch and event-related responses are modeled in exactly the same way but by specifying their different durations.  Events are specified with a duration of 0.  If you enter a single number for the durations it will be assumed that all trials conform to this duration. If you have multiple different durations, then the number must match the number of onset times.'};
duration.strtype = 'r';
duration.num     = [Inf 1];
% ---------------------------------------------------------------------
% tmod Time Modulation
% ---------------------------------------------------------------------
tmod         = cfg_menu;
tmod.tag     = 'tmod';
tmod.name    = 'Time Modulation';
tmod.help    = {
                'This option allows for the characterisation of linear or nonlinear time effects. For example, 1st order modulation would model the stick functions and a linear change of the stick function heights over time. Higher order modulation will introduce further columns that contain the stick functions scaled by time squared, time cubed etc.'
                ''
                'Interactions or response modulations can enter at two levels.  Firstly the stick function itself can be modulated by some parametric variate (this can be time or some trial-specific variate like reaction time) modeling the interaction between the trial and the variate or, secondly interactions among the trials themselves can be modeled using a Volterra series formulation that accommodates interactions over time (and therefore within and between trial types).'
}';
tmod.labels = {
               'No Time Modulation'
               '1st order Time Modulation'
               '2nd order Time Modulation'
               '3rd order Time Modulation'
               '4th order Time Modulation'
               '5th order Time Modulation'
               '6th order Time Modulation'
}';
tmod.values = {0 1 2 3 4 5 6};
tmod.val    = {0};
% ---------------------------------------------------------------------
% name Name
% ---------------------------------------------------------------------
name1         = cfg_entry;
name1.tag     = 'name';
name1.name    = 'Name';
name1.help    = {'Enter a name for this parameter.'};
name1.strtype = 's';
name1.num     = [1 Inf];
% ---------------------------------------------------------------------
% param Values
% ---------------------------------------------------------------------
param         = cfg_entry;
param.tag     = 'param';
param.name    = 'Values';
param.help    = {'Enter a vector of values, one for each occurence of the event.'};
param.strtype = 'r';
param.num     = [Inf 1];
% ---------------------------------------------------------------------
% poly Polynomial Expansion
% ---------------------------------------------------------------------
poly         = cfg_menu;
poly.tag     = 'poly';
poly.name    = 'Polynomial Expansion';
poly.help    = {'For example, 1st order modulation would model the stick functions and a linear change of the stick function heights over different values of the parameter. Higher order modulation will introduce further columns that contain the stick functions scaled by parameter squared, cubed etc.'};
poly.labels = {
               '1st order'
               '2nd order'
               '3rd order'
               '4th order'
               '5th order'
               '6th order'
}';
poly.values = {1 2 3 4 5 6};
% ---------------------------------------------------------------------
% pmod Parameter
% ---------------------------------------------------------------------
pmod         = cfg_branch;
pmod.tag     = 'pmod';
pmod.name    = 'Parameter';
pmod.val     = {name1 param poly };
pmod.help    = {
                'Model interactions with user specified parameters. This allows nonlinear effects relating to some other measure to be modelled in the design matrix.'
                ''
                'Interactions or response modulations can enter at two levels.  Firstly the stick function itself can be modulated by some parametric variate (this can be time or some trial-specific variate like reaction time) modeling the interaction between the trial and the variate or, secondly interactions among the trials themselves can be modeled using a Volterra series formulation that accommodates interactions over time (and therefore within and between trial types).'
}';
% ---------------------------------------------------------------------
% generic Parametric Modulations
% ---------------------------------------------------------------------
generic2         = cfg_repeat;
generic2.tag     = 'generic';
generic2.name    = 'Parametric Modulations';
generic2.help    = {'The stick function itself can be modulated by some parametric variate (this can be time or some trial-specific variate like reaction time) modeling the interaction between the trial and the variate. The events can be modulated by zero or more parameters.'};
generic2.values  = {pmod };
generic2.num     = [0 Inf];
% ---------------------------------------------------------------------
% porth Orthogonalise modulations
% ---------------------------------------------------------------------
porth         = cfg_menu;
porth.tag     = 'orth';
porth.name    = 'Orthogonalise modulations';
porth.help    = {'Orthogonalise regressors within trial types.'};
porth.labels  = {'Yes' 'No'};
porth.values  = {1 0};
porth.val     = {1};
% ---------------------------------------------------------------------
% cond Condition
% ---------------------------------------------------------------------
cond         = cfg_branch;
cond.tag     = 'cond';
cond.name    = 'Condition';
cond.val     = {name onset duration tmod generic2 porth};
cond.check   = @cond_check;
cond.help    = {'An array of input functions is contructed, specifying occurrence events or epochs (or both). These are convolved with a basis set at a later stage to give regressors that enter into the design matrix. Interactions of evoked responses with some parameter (time or a specified variate) enter at this stage as additional columns in the design matrix with each trial multiplied by the [expansion of the] trial-specific parameter. The 0th order expansion is simply the main effect in the first column.'};
% ---------------------------------------------------------------------
% generic Conditions
% ---------------------------------------------------------------------
generic1         = cfg_repeat;
generic1.tag     = 'generic';
generic1.name    = 'Conditions';
generic1.help    = {'You are allowed to combine both event- and epoch-related responses in the same model and/or regressor. Any number of condition (event or epoch) types can be specified.  Epoch and event-related responses are modeled in exactly the same way by specifying their onsets [in terms of onset times] and their durations.  Events are specified with a duration of 0.  If you enter a single number for the durations it will be assumed that all trials conform to this duration.For factorial designs, one can later associate these experimental conditions with the appropriate levels of experimental factors. '};
generic1.values  = {cond };
generic1.num     = [0 Inf];
% ---------------------------------------------------------------------
% multi Multiple conditions
% ---------------------------------------------------------------------
multi         = cfg_files;
multi.tag     = 'multi';
multi.name    = 'Multiple conditions';
multi.val     = {{''}};
multi.help    = {
                 'Select the *.mat file containing details of your multiple experimental conditions. '
                 ''
                 'If you have multiple conditions then entering the details a condition at a time is very inefficient. This option can be used to load all the required information in one go. You will first need to create a *.mat file containing the relevant information. '
                 ''
                 'This *.mat file must include the following cell arrays (each 1 x n): names, onsets and durations. eg. names=cell(1,5), onsets=cell(1,5), durations=cell(1,5), then names{2}=''SSent-DSpeak'', onsets{2}=[3 5 19 222], durations{2}=[0 0 0 0], contain the required details of the second condition. These cell arrays may be made available by your stimulus delivery program, eg. COGENT. The duration vectors can contain a single entry if the durations are identical for all events. Optionally, a (1 x n) cell array named orth can also be included, with a 1 or 0 for each condition to indicate whether parameteric modulators should be orthogonalised.'
                 ''
                 'Time and Parametric effects can also be included. For time modulation include a cell array (1 x n) called tmod. It should have a have a single number in each cell. Unused cells may contain either a 0 or be left empty. The number specifies the order of time modulation from 0 = No Time Modulation to 6 = 6th Order Time Modulation. eg. tmod{3} = 1, modulates the 3rd condition by a linear time effect.'
                 ''
                 'For parametric modulation include a structure array, which is up to 1 x n in size, called pmod. n must be less than or equal to the number of cells in the names/onsets/durations cell arrays. The structure array pmod must have the fields: name, param and poly.  Each of these fields is in turn a cell array to allow the inclusion of one or more parametric effects per column of the design. The field name must be a cell array containing strings. The field param is a cell array containing a vector of parameters. Remember each parameter must be the same length as its corresponding onsets vector. The field poly is a cell array (for consistency) with each cell containing a single number specifying the order of the polynomial expansion from 1 to 6.'
                 ''
                 'Note that each condition is assigned its corresponding entry in the structure array (condition 1 parametric modulators are in pmod(1), condition 2 parametric modulators are in pmod(2), etc. Within a condition multiple parametric modulators are accessed via each fields cell arrays. So for condition 1, parametric modulator 1 would be defined in  pmod(1).name{1}, pmod(1).param{1}, and pmod(1).poly{1}. A second parametric modulator for condition 1 would be defined as pmod(1).name{2}, pmod(1).param{2} and pmod(1).poly{2}. If there was also a parametric modulator for condition 2, then remember the first modulator for that condition is in cell array 1: pmod(2).name{1}, pmod(2).param{1}, and pmod(2).poly{1}. If some, but not all conditions are parametrically modulated, then the non-modulated indices in the pmod structure can be left blank. For example, if conditions 1 and 3 but not condition 2 are modulated, then specify pmod(1) and pmod(3). Similarly, if conditions 1 and 2 are modulated but there are 3 conditions overall, it is only necessary for pmod to be a 1 x 2 structure array.'
                 ''
                 'EXAMPLE:'
                 'Make an empty pmod structure: '
                 '  pmod = struct(''name'',{''''},''param'',{},''poly'',{});'
                 'Specify one parametric regressor for the first condition: '
                 '  pmod(1).name{1}  = ''regressor1'';'
                 '  pmod(1).param{1} = [1 2 4 5 6];'
                 '  pmod(1).poly{1}  = 1;'
                 'Specify 2 parametric regressors for the second condition: '
                 '  pmod(2).name{1}  = ''regressor2-1'';'
                 '  pmod(2).param{1} = [1 3 5 7]; '
                 '  pmod(2).poly{1}  = 1;'
                 '  pmod(2).name{2}  = ''regressor2-2'';'
                 '  pmod(2).param{2} = [2 4 6 8 10];'
                 '  pmod(2).poly{2}  = 1;'
                 ''
                 'The parametric modulator should be mean corrected if appropriate. Unused structure entries should have all fields left empty.'
}';
multi.filter = 'mat';
multi.ufilter = '.*';
multi.num     = [0 1];
% ---------------------------------------------------------------------
% name Name
% ---------------------------------------------------------------------
name         = cfg_entry;
name.tag     = 'name';
name.name    = 'Name';
name.help    = {'Enter name of regressor eg. First movement parameter'};
name.strtype = 's';
name.num     = [1 Inf];
% ---------------------------------------------------------------------
% val Value
% ---------------------------------------------------------------------
val         = cfg_entry;
val.tag     = 'val';
val.name    = 'Value';
val.help    = {'Enter the vector of regressor values'};
val.strtype = 'r';
val.num     = [Inf 1];
% ---------------------------------------------------------------------
% regress Regressor
% ---------------------------------------------------------------------
regress         = cfg_branch;
regress.tag     = 'regress';
regress.name    = 'Regressor';
regress.val     = {name val };
regress.help    = {'regressor'};
% ---------------------------------------------------------------------
% generic Regressors
% ---------------------------------------------------------------------
generic2         = cfg_repeat;
generic2.tag     = 'generic';
generic2.name    = 'Regressors';
generic2.help    = {'Regressors are additional columns included in the design matrix, which may model effects that would not be convolved with the haemodynamic response.  One such example would be the estimated movement parameters, which may confound the data.'};
generic2.values  = {regress };
generic2.num     = [0 Inf];
% ---------------------------------------------------------------------
% multi_reg Multiple regressors
% ---------------------------------------------------------------------
multi_reg         = cfg_files;
multi_reg.tag     = 'multi_reg';
multi_reg.name    = 'Multiple regressors';
multi_reg.val     = {{''}};
multi_reg.help    = {
                     'Select the *.mat/*.txt file(s) containing details of your multiple regressors. '
                     ''
                     'If you have multiple regressors eg. realignment parameters, then entering the details a regressor at a time is very inefficient. This option can be used to load all the required information in one go. '
                     ''
                     'You will first need to create a *.mat file containing a matrix R or a *.txt file containing the regressors. Each column of R will contain a different regressor. When SPM creates the design matrix the regressors will be named R1, R2, R3, ..etc.'
                     ''
                     'You can also select a PPI.mat file and SPM will automatically create regressors from fields PPI.ppi, PPI.Y and PPI.P.'
                     }';
multi_reg.filter = 'mat';
multi_reg.ufilter = '.*';
multi_reg.num     = [0 Inf];
% ---------------------------------------------------------------------
% hpf High-pass filter
% ---------------------------------------------------------------------
hpf         = cfg_entry;
hpf.tag     = 'hpf';
hpf.name    = 'High-pass filter';
hpf.help    = {'The default high-pass filter cutoff is 128 seconds.Slow signal drifts with a period longer than this will be removed. Use ''explore design'' to ensure this cut-off is not removing too much experimental variance. High-pass filtering is implemented using a residual forming matrix (i.e. it is not a convolution) and is simply to a way to remove confounds without estimating their parameters explicitly.  The constant term is also incorporated into this filter matrix.'};
hpf.strtype = 'r';
hpf.num     = [1 1];
hpf.def     = @(val)spm_get_defaults('stats.fmri.hpf', val{:});
% ---------------------------------------------------------------------
% sess Subject/Session
% ---------------------------------------------------------------------
sess         = cfg_branch;
sess.tag     = 'sess';
sess.name    = 'Subject/Session';
sess.val     = {nscan generic1 multi generic2 multi_reg hpf };
sess.check   = @sess_check;
sess.help    = {'The design matrix for fMRI data consists of one or more separable, session-specific partitions.  These partitions are usually either one per subject, or one per fMRI scanning session for that subject.'};
% ---------------------------------------------------------------------
% generic Data & Design
% ---------------------------------------------------------------------
generic         = cfg_repeat;
generic.tag     = 'generic';
generic.name    = 'Data & Design';
generic.help    = {
                   'The design matrix defines the experimental design and the nature of hypothesis testing to be implemented.  The design matrix has one row for each scan and one column for each effect or explanatory variable. (e.g. regressor or stimulus function).  '
                   ''
                   'This allows you to build design matrices with separable session-specific partitions.  Each partition may be the same (in which case it is only necessary to specify it once) or different.  Responses can be either event- or epoch related, where the latter model involves prolonged and possibly time-varying responses to state-related changes in experimental conditions.  Event-related response are modelled in terms of responses to instantaneous events.  Mathematically they are both modelled by convolving a series of delta (stick) or box-car functions, encoding the input or stimulus function. with a set of hemodynamic basis functions.'
}';
generic.values  = {sess };
generic.num     = [1 Inf];
% ---------------------------------------------------------------------
% name Name
% ---------------------------------------------------------------------
name         = cfg_entry;
name.tag     = 'name';
name.name    = 'Name';
name.help    = {'Name of factor, eg. ''Repetition'' '};
name.strtype = 's';
name.num     = [1 Inf];
% ---------------------------------------------------------------------
% levels Levels
% ---------------------------------------------------------------------
levels         = cfg_entry;
levels.tag     = 'levels';
levels.name    = 'Levels';
levels.help    = {'Enter number of levels for this factor, eg. 2'};
levels.strtype = 'n';
levels.num     = [Inf 1];
% ---------------------------------------------------------------------
% fact Factor
% ---------------------------------------------------------------------
fact         = cfg_branch;
fact.tag     = 'fact';
fact.name    = 'Factor';
fact.val     = {name levels };
fact.help    = {'Add a new factor to your experimental design'};
% ---------------------------------------------------------------------
% generic Factorial design
% ---------------------------------------------------------------------
generic1         = cfg_repeat;
generic1.tag     = 'generic';
generic1.name    = 'Factorial design';
generic1.help    = {
                    'If you have a factorial design then SPM can automatically generate the contrasts necessary to test for the main effects and interactions. '
                    ''
                    'This includes the F-contrasts necessary to test for these effects at the within-subject level (first level) and the simple contrasts necessary to generate the contrast images for a between-subject (second-level) analysis.'
                    ''
                    'To use this option, create as many factors as you need and provide a name and number of levels for each.  SPM assumes that the condition numbers of the first factor change slowest, the second factor next slowest etc. It is best to write down the contingency table for your design to ensure this condition is met. This table relates the levels of each factor to the conditions. '
                    ''
                    'For example, if you have 2-by-3 design  your contingency table has two rows and three columns where the the first factor spans the rows, and the second factor the columns. The numbers of the conditions are 1,2,3 for the first row and 4,5,6 for the second. '
}';
generic1.values  = {fact };
generic1.num     = [0 Inf];
% ---------------------------------------------------------------------
% derivs Model derivatives
% ---------------------------------------------------------------------
derivs         = cfg_menu;
derivs.tag     = 'derivs';
derivs.name    = 'Model derivatives';
derivs.help    = {'Model HRF Derivatives. The canonical HRF combined with time and dispersion derivatives comprise an ''informed'' basis set, as the shape of the canonical response conforms to the hemodynamic response that is commonly observed. The incorporation of the derivate terms allow for variations in subject-to-subject and voxel-to-voxel responses. The time derivative allows the peak response to vary by plus or minus a second and the dispersion derivative allows the width of the response to vary. The informed basis set requires an SPM{F} for inference. T-contrasts over just the canonical are perfectly valid but assume constant delay/dispersion. The informed basis set compares favourably with eg. FIR bases on many data sets. '};
derivs.labels = {
                 'No derivatives'
                 'Time derivatives'
                 'Time and Dispersion derivatives'
}';
derivs.values = {[0 0] [1 0] [1 1]};
derivs.val    = {[0 0]};
% ---------------------------------------------------------------------
% hrf Canonical HRF
% ---------------------------------------------------------------------
hrf         = cfg_branch;
hrf.tag     = 'hrf';
hrf.name    = 'Canonical HRF';
hrf.val     = {derivs };
hrf.help    = {'Canonical Hemodynamic Response Function. This is the default option. Contrasts of these effects have a physical interpretation and represent a parsimonious way of characterising event-related responses. This option is also useful if you wish to look separately at activations and deactivations (this is implemented using a t-contrast with a +1 or -1 entry over the canonical regressor). '};
% ---------------------------------------------------------------------
% length Window length
% ---------------------------------------------------------------------
length         = cfg_entry;
length.tag     = 'length';
length.name    = 'Window length';
length.help    = {'Post-stimulus window length (in seconds)'};
length.strtype = 'r';
length.num     = [1 1];
% ---------------------------------------------------------------------
% order Order
% ---------------------------------------------------------------------
order         = cfg_entry;
order.tag     = 'order';
order.name    = 'Order';
order.help    = {'Number of basis functions'};
order.strtype = 'n';
order.num     = [1 1];
% ---------------------------------------------------------------------
% fourier Fourier Set
% ---------------------------------------------------------------------
fourier         = cfg_branch;
fourier.tag     = 'fourier';
fourier.name    = 'Fourier Set';
fourier.val     = {length order };
fourier.help    = {'Fourier basis functions. This option requires an SPM{F} for inference.'};
% ---------------------------------------------------------------------
% length Window length
% ---------------------------------------------------------------------
length         = cfg_entry;
length.tag     = 'length';
length.name    = 'Window length';
length.help    = {'Post-stimulus window length (in seconds)'};
length.strtype = 'r';
length.num     = [1 1];
% ---------------------------------------------------------------------
% order Order
% ---------------------------------------------------------------------
order         = cfg_entry;
order.tag     = 'order';
order.name    = 'Order';
order.help    = {'Number of basis functions'};
order.strtype = 'n';
order.num     = [1 1];
% ---------------------------------------------------------------------
% fourier_han Fourier Set (Hanning)
% ---------------------------------------------------------------------
fourier_han         = cfg_branch;
fourier_han.tag     = 'fourier_han';
fourier_han.name    = 'Fourier Set (Hanning)';
fourier_han.val     = {length order };
fourier_han.help    = {'Fourier basis functions with Hanning Window - requires SPM{F} for inference.'};
% ---------------------------------------------------------------------
% length Window length
% ---------------------------------------------------------------------
length         = cfg_entry;
length.tag     = 'length';
length.name    = 'Window length';
length.help    = {'Post-stimulus window length (in seconds)'};
length.strtype = 'r';
length.num     = [1 1];
% ---------------------------------------------------------------------
% order Order
% ---------------------------------------------------------------------
order         = cfg_entry;
order.tag     = 'order';
order.name    = 'Order';
order.help    = {'Number of basis functions'};
order.strtype = 'n';
order.num     = [1 1];
% ---------------------------------------------------------------------
% gamma Gamma Functions
% ---------------------------------------------------------------------
gamma         = cfg_branch;
gamma.tag     = 'gamma';
gamma.name    = 'Gamma Functions';
gamma.val     = {length order };
gamma.help    = {'Gamma basis functions - requires SPM{F} for inference.'};
% ---------------------------------------------------------------------
% length Window length
% ---------------------------------------------------------------------
length         = cfg_entry;
length.tag     = 'length';
length.name    = 'Window length';
length.help    = {'Post-stimulus window length (in seconds)'};
length.strtype = 'r';
length.num     = [1 1];
% ---------------------------------------------------------------------
% order Order
% ---------------------------------------------------------------------
order         = cfg_entry;
order.tag     = 'order';
order.name    = 'Order';
order.help    = {'Number of basis functions'};
order.strtype = 'n';
order.num     = [1 1];
% ---------------------------------------------------------------------
% fir Finite Impulse Response
% ---------------------------------------------------------------------
fir         = cfg_branch;
fir.tag     = 'fir';
fir.name    = 'Finite Impulse Response';
fir.val     = {length order };
fir.help    = {'Finite impulse response - requires SPM{F} for inference.'};
% ---------------------------------------------------------------------
% bases Basis Functions
% ---------------------------------------------------------------------
bases         = cfg_choice;
bases.tag     = 'bases';
bases.name    = 'Basis Functions';
bases.val     = {hrf };
bases.help    = {'The most common choice of basis function is the Canonical HRF with or without time and dispersion derivatives. '};
bases.values  = {hrf fourier fourier_han gamma fir };
% ---------------------------------------------------------------------
% volt Model Interactions (Volterra)
% ---------------------------------------------------------------------
volt         = cfg_menu;
volt.tag     = 'volt';
volt.name    = 'Model Interactions (Volterra)';
volt.help    = {
                'Generalized convolution of inputs (U) with basis set (bf).'
                ''
                'For first order expansions the causes are simply convolved (e.g. stick functions) in U.u by the basis functions in bf to create a design matrix X.  For second order expansions new entries appear in ind, bf and name that correspond to the interaction among the orginal causes. The basis functions for these efects are two dimensional and are used to assemble the second order kernel. Second order effects are computed for only the first column of U.u.'
                'Interactions or response modulations can enter at two levels.  Firstly the stick function itself can be modulated by some parametric variate (this can be time or some trial-specific variate like reaction time) modeling the interaction between the trial and the variate or, secondly interactions among the trials themselves can be modeled using a Volterra series formulation that accommodates interactions over time (and therefore within and between trial types).'
}';
volt.labels = {
               'Do not model Interactions'
               'Model Interactions'
}';
volt.values = {1 2};
volt.val    = {1};
% ---------------------------------------------------------------------
% global Global normalisation
% ---------------------------------------------------------------------
xGlobal         = cfg_menu;
xGlobal.tag     = 'global';
xGlobal.name    = 'Global normalisation';
xGlobal.help    = {'Global intensity normalisation'};
xGlobal.labels = {
                  'Scaling'
                  'None'
}';
xGlobal.values = {
                  'Scaling'
                  'None'
}';
xGlobal.val     = {'None'};
%----------------------------------------------------------------------
% gMT Masking threshold
%----------------------------------------------------------------------
gMT         = cfg_entry;
gMT.tag     = 'mthresh';
gMT.name    = 'Masking threshold';
gMT.help    = {'Masking threshold, defined as proportion of globals.'};
gMT.strtype = 'r';
gMT.num     = [1 1];
gMT.def     = @(val)spm_get_defaults('mask.thresh', val{:});
% ---------------------------------------------------------------------
% cvi Serial correlations
% ---------------------------------------------------------------------
cvi         = cfg_menu;
cvi.tag     = 'cvi';
cvi.name    = 'Serial correlations';
cvi.help    = {
               'Serial correlations in fMRI time series due to aliased biorhythms and unmodelled neuronal activity can be accounted for using an autoregressive AR(1) model during Classical (ReML) parameter estimation.  '
               ''
               'This estimate assumes the same correlation structure for each voxel, within each session.  ReML estimates are then used to correct for non-sphericity during inference by adjusting the statistics and degrees of freedom appropriately.  The discrepancy between estimated and actual intrinsic (i.e. prior to filtering) correlations are greatest at low frequencies.  Therefore specification of the high-pass filter is particularly important. '
               ''
               'Serial correlation can be ignored if you choose the ''none'' option. Note that the above options only apply if you later specify that your model will be estimated using the Classical (ReML) approach. If you choose Bayesian estimation these options will be ignored. For Bayesian estimation, the choice of noisemodel (AR model order) is made under the estimation options. '
}';
cvi.labels  = {'none', 'AR(1)', 'FAST'};
cvi.values  = {'none', 'AR(1)', 'FAST'};
cvi.def     = @(val)spm_get_defaults('stats.fmri.cvi', val{:});
% ---------------------------------------------------------------------
% fmri_design fMRI model specification (design only)
% ---------------------------------------------------------------------
fmri_design         = cfg_exbranch;
fmri_design.tag     = 'fmri_design';
fmri_design.name    = 'fMRI model specification (design only)';
fmri_design.val     = {dir timing generic generic1 bases volt xGlobal gMT cvi };
fmri_design.help    = {
    'Specification of a General Linear Model for statistical analysis of fMRI data (design only).'
    ''
    'It comprises the following steps (1) specification of the GLM design matrix, fMRI data files and filtering (2) estimation of GLM paramaters using classical or Bayesian approaches and (3) interrogation of results using contrast vectors to produce Statistical Parametric Maps (SPMs) or Posterior Probability Maps (PPMs).'
    ''
    'The design matrix defines the experimental design and the nature of hypothesis testing to be implemented.  The design matrix has one row for each scan and one column for each effect or explanatory variable. (eg. regressor or stimulus function). You can build design matrices with separable session-specific partitions.  Each partition may be the same (in which case it is only necessary to specify it once) or different. '
    ''
    'Responses can be either event- or epoch related, the only distinction is the duration of the underlying input or stimulus function. Mathematically they are both modeled by convolving a series of delta (stick) or box functions (u), indicating the onset of an event or epoch with a set of basis functions.  These basis functions model the hemodynamic convolution, applied by the brain, to the inputs.  This convolution can be first-order or a generalized convolution modeled to second order (if you specify the Volterra option). The same inputs are used by the Hemodynamic model or Dynamic Causal Models which model the convolution explicitly in terms of hidden state variables. '
    ''
    'Basis functions can be used to plot estimated responses to single events once the parameters (i.e. basis function coefficients) have been estimated.  The importance of basis functions is that they provide a graceful transition between simple fixed response models (like the box-car) and finite impulse response (FIR) models, where there is one basis function for each scan following an event or epoch onset.  The nice thing about basis functions, compared to FIR models, is that data sampling and stimulus presentation does not have to be synchronized thereby allowing a uniform and unbiased sampling of peri-stimulus time.'
    ''
    'Event-related designs may be stochastic or deterministic.  Stochastic designs involve one of a number of trial-types occurring with a specified probability at successive intervals in time.  These probabilities can be fixed (stationary designs) or time-dependent (modulated or non-stationary designs).  The most efficient designs obtain when the probabilities of every trial type are equal. A critical issue in stochastic designs is whether to include null events If you wish to estimate the evoked response to a specific event type (as opposed to differential responses) then a null event must be included (even if it is not modeled explicitly).'
    ''
    'In SPM, analysis of data from multiple subjects typically proceeds in two stages using models at two ''levels''. The ''first level'' models are used to implement a within-subject analysis. Typically there will be as many first level models as there are subjects. Analysis proceeds as described using the ''Specify first level'' and ''Estimate'' options. The results of these analyses can then be presented as ''case studies''. More often, however, one wishes to make inferences about the population from which the subjects were drawn. This is an example of a ''Random-Effects (RFX) analysis'' (or, more properly, a mixed-effects analysis). In SPM, RFX analysis is implemented using the ''summary-statistic'' approach where contrast images from each subject are used as summary measures of subject responses. These are then entered as data into a ''second level'' model. '
    }';
fmri_design.prog = @spm_run_fmri_spec;
fmri_design.vout = @vout_stats;
fmri_design.modality = {'FMRI'};

%==========================================================================
function t = cond_check(job)
t   = {};
if (numel(job.onset) ~= numel(job.duration)) && (numel(job.duration)~=1)
    t = {sprintf('"%s": Number of event onsets (%d) does not match the number of durations (%d).',...
        job.name, numel(job.onset),numel(job.duration))};
end
for i=1:numel(job.pmod)
    if numel(job.onset) ~= numel(job.pmod(i).param)
        t = {t{:}, sprintf('"%s" & "%s":Number of event onsets (%d) does not equal the number of parameters (%d).',...
            job.name, job.pmod(i).name, numel(job.onset),numel(job.pmod(i).param))};
    end
end

%==========================================================================
function t = sess_check(sess)
t = {};
for i=1:numel(sess.regress)
    if sess.nscan ~= numel(sess.regress(i).val)
        t = {t{:}, sprintf('Num scans (%d) ~= Num regress[%d] (%d).',numel(sess.nscan),i,numel(sess.regress(i).val))};
    end
end

%==========================================================================
function dep = vout_stats(job)
dep(1)            = cfg_dep;
dep(1).sname      = 'SPM.mat File';
dep(1).src_output = substruct('.','spmmat');
dep(1).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});
%dep(2)            = cfg_dep;
%dep(2).sname      = 'SPM Variable';
%dep(2).src_output = substruct('.','spmvar');
%dep(2).tgt_spec   = cfg_findspec({{'strtype','e'}});

author: spm


text:
function M = Pythagor_tree(m,n,Colormap)
% function M = Pythagor_tree(m,n,Colormap)
% Compute Pythagoras_tree
% The Pythagoras Tree is a plane fractal constructed from squares.
% It is named after Pythagoras  because each triple of touching squares
% encloses a right triangle, in a configuration traditionally used to
% depict the Pythagorean theorem.
% http://en.wikipedia.org/wiki/Pythagoras_tree
%
% Input :
%       - m ( double m> 0) is the relative length of one of the side
%         right-angled triangle. The second side of the right-angle is
%         taken to be one.
%         To have a symmetric tree, m has to be 1.
%       - n ( integer ) is the level of recursion.
%         The number of elements of tree is equal to 2^(n+1)-1.
%         A reasonnable number for n is 10.
%       - Colormap: String used to generate color of the different levels
%         of the tree.
%       All these arguments are optional: the function can run with
%       argument.
% Output :
%       - Matrix M: Pyhagoras tree is stored in a matrix M.
%         This matrix has 5 columns.
%         Each row corresponds to the cordinate of each square of the tree
%         The two first columns give the bottom-left position of each
%         square. The third column corresponds to the orientation angle of
%         each square. The fourth column gives the size of each square. The
%         fifth column specifies the level of recursion of each square.
%         The first row corresponds to the root of the tree. It is always
%         the same
%         M(1,:) = [0 -1 0 1 1];
%         The leaf located at row i will give 2 leaves located at 2*i and
%         2*i+1.
%       - A svg file giving a vectorial display of the tree. The name of
%         file is generated from the parameter m,n,Colormap. The file is
%         stored in the current folder.
%
% 2010 02 29
% Guillaume Jacquenot
% guillaume dot jacquenot at gmail dot com

%% Check inputs
narg = nargin;
if narg <= 2
    % Colormap = 'jet';
    Colormap = 'summer';
    if narg <= 1
        n = 12; % Recursion level
        if nargin == 0
            m = 0.8;
        end
    end
end
if m <= 0
	error([mfilename ':e0'],'Length of m has to be greater than zero');
end
if rem(n,1)~=0
	error([mfilename ':e0'],'The number of level has to be integer');
end
if ~iscolormap(Colormap)
	error([mfilename ':e1'],'Input colormap is not valid');
end
%% Compute constants
d      = sqrt(1+m^2);                  %
c1     = 1/d;                          % Normalized length 1
c2     = m/d;                          % Normalized length 2
T      = [0 1/(1+m^2);1 1+m/(1+m^2)];  % Translation pattern
alpha1 = atan2(m,1);                   % Defines the first rotation angle
alpha2 = alpha1-pi/2;                  % Defines the second rotation angle
pi2    = 2*pi;                         % Defines pi2
nEle   = 2^(n+1)-1;                    % Number of elements (square)
M      = zeros(nEle,5);                % Matrice containing the tree
M(1,:) = [0 -1 0 1 1];                 % Initialization of the tree

%% Compute the level of each square contained in the resulting matrix
Offset = 0;
for i = 0:n
    tmp = 2^i;
    M(Offset+(1:tmp),5) = i;
    Offset = Offset + tmp;
end
%% Compute the position and size of each square wrt its parent
for i = 2:2:(nEle-1)
    j          = i/2;
    mT         = M(j,4) * mat_rot(M(j,3)) * T;
    Tx         = mT(1,:) + M(j,1);
    Ty         = mT(2,:) + M(j,2);
    theta1     = rem(M(j,3)+alpha1,pi2);
    theta2     = rem(M(j,3)+alpha2,pi2);
    M(i  ,1:4) = [Tx(1) Ty(1) theta1 M(j,4)*c1];
    M(i+1,1:4) = [Tx(2) Ty(2) theta2 M(j,4)*c2];
end
%% Display the tree
Pythagor_tree_plot(M,n);

%% Write results to an SVG file
Pythagor_tree_write2svg(m,n,Colormap,M);

function Pythagor_tree_write2svg(m,n,Colormap,M)
% Determine the bounding box of the tree with an offset
% Display_metadata = false;
Display_metadata = true;

nEle    = size(M,1);
r2      = sqrt(2);
LOffset = M(nEle,4) + 0.1;
min_x   = min(M(:,1)-r2*M(:,4)) - LOffset;
max_x   = max(M(:,1)+r2*M(:,4)) + LOffset;
min_y   = min(M(:,2)          ) - LOffset;  % -r2*M(:,4)
max_y   = max(M(:,2)+r2*M(:,4)) + LOffset;

% Compute the color of tree
ColorM = zeros(n+1,3);
eval(['ColorM = flipud(' Colormap '(n+1));']);
co   = 100;
Wfig = ceil(co*(max_x-min_x));
Hfig = ceil(co*(max_y-min_y));
filename = ['Pythagoras_tree_1_' strrep(num2str(m),'.','_') '_'...
             num2str(n) '_' Colormap '.svg'];
fid  = fopen(filename, 'wt');
fprintf(fid,'<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n');
if ~Display_metadata
    fprintf(fid,'<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"\n');
    fprintf(fid,'  "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">\n');
end
fprintf(fid,'<svg width="%d" height="%d" version="1.1"\n',Wfig,Hfig); %
% fprintf(fid,['<svg width="12cm" height="4cm" version="1.1"\n']); % Wfig,

% fprintf(fid,['<svg width="15cm" height="10cm" '...
%              'viewBox="0 0 %d %d" version="1.1"\n'],...
%              Wfig,Hfig);
if Display_metadata
    fprintf(fid,'\txmlns:dc="http://purl.org/dc/elements/1.1/"\n');
    fprintf(fid,'\txmlns:cc="http://creativecommons.org/ns#"\n');
    fprintf(fid,['\txmlns:rdf="http://www.w3.org/1999/02/22'...
                 '-rdf-syntax-ns#"\n']);
end
fprintf(fid,'\txmlns:svg="http://www.w3.org/2000/svg"\n');
fprintf(fid,'\txmlns="http://www.w3.org/2000/svg"\n');
fprintf(fid,'\txmlns:xlink="http://www.w3.org/1999/xlink">\n');

if Display_metadata
    fprintf(fid,'\t<title>Pythagoras tree</title>\n');
    fprintf(fid,'\t<metadata>\n');
    fprintf(fid,'\t\t<rdf:RDF>\n');
    fprintf(fid,'\t\t\t<cc:Work\n');
    fprintf(fid,'\t\t\t\trdf:about="">\n');
    fprintf(fid,'\t\t\t\t<dc:format>image/svg+xml</dc:format>\n');
    fprintf(fid,'\t\t\t\t<dc:type\n');
    fprintf(fid,'\t\t\t\t\trdf:resource="http://purl.org/dc/dcmitype/StillImage" />\n');
    fprintf(fid,'\t\t\t\t<dc:title>Pythagoras tree</dc:title>\n');
    fprintf(fid,'\t\t\t\t<dc:creator>\n');
    fprintf(fid,'\t\t\t\t\t<cc:Agent>\n');
    fprintf(fid,'\t\t\t\t\t\t<dc:title>Guillaume Jacquenot</dc:title>\n');
    fprintf(fid,'\t\t\t\t\t</cc:Agent>\n');
    fprintf(fid,'\t\t\t\t</dc:creator>\n');
    fprintf(fid,'\t\t\t\t<cc:license\n');
    fprintf(fid,'\t\t\t\t\t\trdf:resource="http://creativecommons.org/licenses/by-nc-sa/3.0/" />\n');
    fprintf(fid,'\t\t\t</cc:Work>\n');
    fprintf(fid,'\t\t\t<cc:License\n');
    fprintf(fid,'\t\t\t\trdf:about="http://creativecommons.org/licenses/by-nc-sa/3.0/">\n');
    fprintf(fid,'\t\t\t\t<cc:permits\n');
    fprintf(fid,'\t\t\t\t\trdf:resource="http://creativecommons.org/ns#Reproduction" />\n');
    fprintf(fid,'\t\t\t\t<cc:permits\n');
    fprintf(fid,'\t\t\t\t\trdf:resource="http://creativecommons.org/ns#Reproduction" />\n');
    fprintf(fid,'\t\t\t\t<cc:permits\n');
    fprintf(fid,'\t\t\t\t\trdf:resource="http://creativecommons.org/ns#Distribution" />\n');
    fprintf(fid,'\t\t\t\t<cc:requires\n');
    fprintf(fid,'\t\t\t\t\trdf:resource="http://creativecommons.org/ns#Notice" />\n');
    fprintf(fid,'\t\t\t\t<cc:requires\n');
    fprintf(fid,'\t\t\t\t\trdf:resource="http://creativecommons.org/ns#Attribution" />\n');
    fprintf(fid,'\t\t\t\t<cc:prohibits\n');
    fprintf(fid,'\t\t\t\t\trdf:resource="http://creativecommons.org/ns#CommercialUse" />\n');
    fprintf(fid,'\t\t\t\t<cc:permits\n');
    fprintf(fid,'\t\t\t\t\trdf:resource="http://creativecommons.org/ns#DerivativeWorks" />\n');
    fprintf(fid,'\t\t\t\t<cc:requires\n');
    fprintf(fid,'\t\t\t\t\trdf:resource="http://creativecommons.org/ns#ShareAlike" />\n');
    fprintf(fid,'\t\t\t</cc:License>\n');
    fprintf(fid,'\t\t</rdf:RDF>\n');
    fprintf(fid,'\t</metadata>\n');
end
fprintf(fid,'\t<defs>\n');
fprintf(fid,'\t\t<rect width="%d" height="%d" \n',co,co);
fprintf(fid,'\t\t\tx="0" y="0"\n');
fprintf(fid,'\t\t\tstyle="fill-opacity:1;stroke:#00d900;stroke-opacity:1"\n');
fprintf(fid,'\t\t\tid="squa"\n');
fprintf(fid,'\t\t/>	\n');
fprintf(fid,'\t</defs>\n');
fprintf(fid,'\t<g transform="translate(%d %d) rotate(180) " >\n',...
                round(co*max_x),round(co*max_y));
for i = 0:n
    fprintf(fid,'\t\t<g style="fill:#%s;" >\n',...
                generate_color_hexadecimal(ColorM(i+1,:)));
    Offset = 2^i-1;
    for j = 1:2^i
        k = j + Offset;
        fprintf(fid,['\t\t\t<use xlink:href="#squa" ',...
                     'transform="translate(%+010.5f %+010.5f)'...
                     ' rotate(%+010.5f) scale(%8.6f)" />\n'],...
                    co*M(k,1),co*M(k,2),M(k,3)*180/pi,M(k,4));
    end
    fprintf(fid,'\t\t</g>\n');
end
fprintf(fid,'\t</g>\n');
fprintf(fid,'</svg>\n');
fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function M = mat_rot(x)
c = cos(x);
s = sin(x);
M=[c -s; s c];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function H = Pythagor_tree_plot(D,ColorM)
if numel(ColorM) == 1
    ColorM = flipud(summer(ColorM+1));
end
H = figure('color','w');
hold on
axis equal
axis off
for i=1:size(D,1)
    cx    = D(i,1);
    cy    = D(i,2);
    theta = D(i,3);
    si    = D(i,4);
    M     = mat_rot(theta);
    x     = si*[0 1 1 0 0];
    y     = si*[0 0 1 1 0];
    pts   = M*[x;y];
    fill(cx+pts(1,:),cy+pts(2,:),ColorM(D(i,5)+1,:));
    % plot(cx+pts(1,1:2),cy+pts(2,1:2),'r');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Scolor = generate_color_hexadecimal(color)
Scolor = '000000';
for i=1:3
    c = dec2hex(round(255*color(i)));
    if numel(c)==1
        Scolor(2*(i-1)+1) = c;
    else
        Scolor(2*(i-1)+(1:2)) = c;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  res = iscolormap(cmap)
% This function returns true if 'cmap' is a valid colormap
LCmap = {...
    'autumn'
    'bone'
    'colorcube'
    'cool'
    'copper'
    'flag'
    'gray'
    'hot'
    'hsv'
    'jet'
    'lines'
    'pink'
    'prism'
    'spring'
    'summer'
    'white'
    'winter'
};

res = ~isempty(strmatch(cmap,LCmap,'exact'));

author: Sable

repo: mcbench-benchmarks


text:
function pcc=pro(pc);

test(1:100)=0;
l=round(100*pc);
test(1:l)=1;
n=round(rand*99)+1;
pcc=test(n);

author: HuangCongQing


text:
function b = dif2_rhs ( m, k )

%*****************************************************************************80
%
%% DIF2_RHS returns the DIF2 right hand side.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the row dimension.
%
%    Input, integer K, the column dimension ( should be 2).
%
%    Output, real B(M,K), the right hand side matrix.
%
  b = zeros ( m, k );

  b(1    ,1) = 1.0;
  b(2:m-1,1) = 0.0;
  b(  m,  1) = 1.0;

  b(1:m-1,2) = 0.0;
  b(  m,  2) = m + 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [err, hndl] = addsquare(Pll, Pur, C)
%
% [err, hndl] = addsquare(Pll, Pur, C)
%
%
%Purpose:
%  adds a rectangle on current plot
%
%
%
%Input Parameters:
%  Pll : [x y] of lower left point
%  Pur : [x y] of upper right point
%  C  : rgb triplet, between 0 and 1  rectangle color
%
%
%
%Output Parameters:
%   err : 0 No Problem
%       : 1 Mucho Problems
%   hndl : the handle to the rectangle object
%
%
%Key Terms:
%
%More Info :
%
%   add the square first, before making the plots,
%   otherwise, the patch will obscure the plots
%
%
%     Author : Ziad Saad
%     Date : Fri Oct 06 15:47:14 EDT 2000
%     LBC/NIMH/ National Institutes of Health, Bethesda Maryland

%Define the function name for easy referencing
FuncName = 'addsquare';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;

X = [Pll(1) Pur(1) Pur(1) Pll(1)]';
Y = [Pll(2) Pll(2) Pur(2) Pur(2)]';

hold on;
hndl = patch(X,Y, C);
set(hndl, 'EdgeColor', C);

err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
if (~exist('syms.m','file'))
    fprintf('The symbolic toolbox is not installed, skipping test\n');
    return;
end

%% Gradients
clear
clc

syms x1 x2
x = [0.1;0.1];

%OPTI Sym
fun = @(x) sin(x(1) + x(2)) + (x(1) - x(2))^2 - 1.5*x(1) + 2.5*x(2) + 1
[jac,jacp] = symJac(fun)
[hess,hessp] = symHess(fun)

%Symbolic
sfun = sin(x1+x2) + (x1-x2)^2 - 1.5*x1 + 2.5*x2 + 1
jac1 = jacobian(sfun)
hess1 = hessian(sfun)

%Eval
[jac(x);double(subs(jac1,{'x1','x2'},{x(1) x(2)}))]
full(jacp())
[hess(x);double(subs(hess1,{'x1','x2'},{x(1) x(2)}))]
full(hessp())

%OPTI Sym
fun = @(x) [x(1) + x(2) - 1;
            x(1)^2 + x(2)^2 - 1];
jac = symJac(fun)

%Symbolic
fun = [x1 + x2 - 1;
       x1^2 + x2^2 - 1];
jac1 = jacobian(fun)

%Eval
[jac(x);double(subs(jac1,{'x1','x2'},{x(1) x(2)}))]

%OPTI Sym
fun = @(x) 0.01*(x(1)-1)^2 + (x(2)-x(1)^2)^2;
[jac,jacp] = symJac(fun)
[hess,hessp] = symHess(fun)

%Symbolic
sfun = 0.01*(x1-1)^2 + (x2-x1^2)^2;
jac1 = jacobian(sfun)
hess1 = hessian(sfun)

%Eval
[jac(x);double(subs(jac1,{'x1','x2'},{x(1) x(2)}))]
full(jacp())
[hess(x);double(subs(hess1,{'x1','x2'},{x(1) x(2)}))]
full(hessp())

%% Symbolic Hessian of Lagrangian
clc
fun = @(x) x(1)*x(4)*(x(1) + x(2) + x(3)) + x(3);
% Nonlinear Constraints
 nlcon = @(x) [ x(1)*x(2)*x(3)*x(4);
                x(1)^2 + x(2)^2 + x(3)^2 + x(4)^2 ];

jac = symJac(fun)
hessObj = symHess(fun)

[hessLag,pattern] = symHessLag(fun,nlcon,4,true)

% Hessian
 H = @(x,sigma,lambda) sparse(sigma*[ 2*x(4)             0      0   0;
                                      x(4)               0      0   0;
                                      x(4)               0      0   0;
                                      2*x(1)+x(2)+x(3)  x(1)  x(1)  0 ] + ...
                      lambda(1)*[   0          0         0         0;
                                 x(3)*x(4)     0         0         0;
                                 x(2)*x(4) x(1)*x(4)     0         0;
                                 x(2)*x(3) x(1)*x(3) x(1)*x(2)     0  ] + ...
                      lambda(2)*diag([2 2 2 2]));

x = 0.1:0.1:0.4;
sigma = 0.3;
lambda = 1:2;

hessLag(x,sigma,lambda)
full(H(x,sigma,lambda))
full(pattern())

%% Application to an NLP, first no derivs
clc
% Objective
 fun = @(x) x(1)*x(4)*(x(1) + x(2) + x(3)) + x(3);

% Nonlinear Constraints
 nlcon = @(x) [ x(1)*x(2)*x(3)*x(4);
                x(1)^2 + x(2)^2 + x(3)^2 + x(4)^2 ];
 cl = [25;40];
 cu = [Inf;40];

% Bounds (lb <= x <= ub)
 lb = ones(4,1);
 ub = 5*ones(4,1);

% Initial Guess
 x0 = [1 5 5 1]';

opts = optiset('display','iter','solver','ipopt','derivCheck','on','solverOpts',ipoptset('derivative_test','first-order'));
Opt = opti('fun',fun,'nl',nlcon,cl,cu,'bounds',lb,ub,'x0',x0,'options',opts)

[x,f,e,i] = solve(Opt)

%% Application to an NLP, now with symJac, symHessLag
clc
% Objective
 fun = @(x) x(1)*x(4)*(x(1) + x(2) + x(3)) + x(3);

grad = symJac(fun);

% Nonlinear Constraints
 nlcon = @(x) [ x(1)*x(2)*x(3)*x(4);
                x(1)^2 + x(2)^2 + x(3)^2 + x(4)^2 ];

[nljac,nljacstr] = symJac(nlcon);

 cl = [25;40];
 cu = [Inf;40];

 %Hessian
 [hess,hstr] = symHessLag(fun,nlcon);

% Bounds (lb <= x <= ub)
 lb = ones(4,1);
 ub = 5*ones(4,1);

% Initial Guess
 x0 = [1 5 5 1]';

opts = optiset('display','iter','solver','ipopt','derivCheck','on','solverOpts',ipoptset('derivative_test','second-order'));
Opt = opti('fun',fun,'grad',grad,'nl',nlcon,cl,cu,'nljac',nljac,'nljacstr',nljacstr,'H',hess,'HStr',hstr,'bounds',lb,ub,'x0',x0,'options',opts)

[x,f,e,i] = solve(Opt)

%% Application to an NLP, analytical derivs supplied
clc
% Objective
 fun = @(x) x(1)*x(4)*(x(1) + x(2) + x(3)) + x(3);

% Gradient
 grad = @(x) [x(1)*x(4) + x(4)*sum(x(1:3)), x(1)*x(4),...
              x(1)*x(4) + 1,  x(1)*sum(x(1:3))];

% Nonlinear Constraints
 nlcon = @(x) [ x(1)*x(2)*x(3)*x(4);
                x(1)^2 + x(2)^2 + x(3)^2 + x(4)^2 ];

% Jacobian
 jac = @(x) sparse([prod(x')./x';
                    2.*x']);

% Jacobian Structure
 jacstr = @() sparse(ones(2,4));

 cl = [25;40];
 cu = [Inf;40];

 %Hessian
% Hessian
 H = @(x,sigma,lambda) sparse(sigma*[ 2*x(4)             0      0   0;
                                      x(4)               0      0   0;
                                      x(4)               0      0   0;
                                      2*x(1)+x(2)+x(3)  x(1)  x(1)  0 ] + ...
                      lambda(1)*[   0          0         0         0;
                                 x(3)*x(4)     0         0         0;
                                 x(2)*x(4) x(1)*x(4)     0         0;
                                 x(2)*x(3) x(1)*x(3) x(1)*x(2)     0  ] + ...
                      lambda(2)*diag([2 2 2 2]));

% Hessian Structure
 Hstr = @() sparse(tril(ones(4)));

% Bounds (lb <= x <= ub)
 lb = ones(4,1);
 ub = 5*ones(4,1);

% Initial Guess
 x0 = [1 5 5 1]';

opts = optiset('display','iter','solver','ipopt','derivCheck','on','solverOpts',ipoptset('derivative_test','second-order'));
Opt = opti('fun',fun,'grad',grad,'nl',nlcon,cl,cu,'nljac',nljac,'nljacstr',nljacstr,'H',H,'HStr',hstr,'bounds',lb,ub,'x0',x0,'options',opts)

[x,f,e,i] = solve(Opt)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function circle_segment_test05 ( )

%*****************************************************************************80
%
%% CIRCLE_SEGMENT_TEST05 tests the AREA and HEIGHT_FROM_AREA functions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_SEGMENT_TEST05\n' );
  fprintf ( 1, '  For circle segment with a given radius R,\n' );
  fprintf ( 1, '  CIRCLE_SEGMENT_AREA_FROM_HEIGHT computes the area A, given the height.\n' );
  fprintf ( 1, '  CIRCLE_SEGMENT_HEIGHT_FROM_AREA computes height H, given the area.\n' );
  fprintf ( 1, '  Check that these functions are inverses of each other\n' );
  fprintf ( 1, '  using random values of R, A, and H.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        R             H      =>     A    =>       H2\n' );
  fprintf ( 1, '\n' );

  seed = 123456789;

  for test = 1 : 5
    [ r, seed ] = r8_uniform_01 ( seed );
    r = 5.0 * r;
    [ h, seed ] = r8_uniform_01 ( seed );
    h = 2.0 * r * h;
    a = circle_segment_area_from_height ( r, h );
    h2 = circle_segment_height_from_area ( r, a );
    fprintf ( 1, '  %12f  %12f  %12f  %12f\n', r, h, a, h2 );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '        R             A      =>     H    =>       A2\n' );
  fprintf ( 1, '\n' );
  for test = 1 : 5
    [ r, seed ] = r8_uniform_01 ( seed );
    r = 5.0 * r;
    [ a, seed ] = r8_uniform_01 ( seed );
    a = pi * r * r * a;
    h = circle_segment_height_from_area ( r, a );
    a2 = circle_segment_area_from_height ( r, h );
    fprintf ( 1, '  %12f  %12f  %12f  %12f\n', r, a, h, a2 );
  end

  return
end


text:
function marker = extractScatterMarker(patch_data)

% EXTRACTS THE MARKER STYLE USED FOR MATLAB OBJECTS
% OF TYPE "PATCH". THESE OBJECTS ARE USED IN AREASERIES
% BARSERIES, CONTOURGROUP, SCATTERGROUP.

%-------------------------------------------------------------------------%

%-AXIS STRUCTURE-%
axis_data = get(ancestor(patch_data.Parent,'axes'));

%-FIGURE STRUCTURE-%
figure_data = get(ancestor(patch_data.Parent,'figure'));

%-INITIALIZE OUTPUT-%
marker = struct();

%-------------------------------------------------------------------------%

%-MARKER SIZEREF-%
marker.sizeref = 1;

%-------------------------------------------------------------------------%

%-MARKER SIZEMODE-%
marker.sizemode = 'area';

%-------------------------------------------------------------------------%

%-MARKER SIZE (STYLE)-%
marker.size = patch_data.SizeData;

%-------------------------------------------------------------------------%

%-MARKER SYMBOL (STYLE)-%
if ~strcmp(patch_data.Marker,'none')

    switch patch_data.Marker
        case '.'
            marksymbol = 'circle';
        case 'o'
            marksymbol = 'circle';
        case 'x'
            marksymbol = 'x-thin-open';
        case '+'
            marksymbol = 'cross-thin-open';
        case '*'
            marksymbol = 'asterisk-open';
        case {'s','square'}
            marksymbol = 'square';
        case {'d','diamond'}
            marksymbol = 'diamond';
        case 'v'
            marksymbol = 'triangle-down';
        case '^'
            marksymbol = 'triangle-up';
        case '<'
            marksymbol = 'triangle-left';
        case '>'
            marksymbol = 'triangle-right';
        case {'p','pentagram'}
            marksymbol = 'star';
        case {'h','hexagram'}
            marksymbol = 'hexagram';
    end

    marker.symbol = marksymbol;
end

%-------------------------------------------------------------------------%

%-MARKER LINE WIDTH (STYLE)-%
marker.line.width = patch_data.LineWidth;

%-------------------------------------------------------------------------%

%--MARKER FILL COLOR--%

%-figure colormap-%
colormap = figure_data.Colormap;

% marker face color
MarkerColor = patch_data.MarkerFaceColor;

filledMarkerSet = {'o','square','s','diamond','d',...
    'v','^', '<','>','hexagram','pentagram'};

filledMarker = ismember(patch_data.Marker,filledMarkerSet);

if filledMarker

    if isnumeric(MarkerColor)
        col = 255*MarkerColor;
        markercolor = ['rgb(' num2str(col(1)) ',' num2str(col(2)) ',' num2str(col(3)) ')'];
    else
        switch MarkerColor

            case 'none'

                markercolor = 'rgba(0,0,0,0)';

            case 'auto'

                if ~strcmp(axis_data.Color,'none')
                    col = 255*axis_data.Color;
                else
                    col = 255*figure_data.Color;
                end

                markercolor  = ['rgb(' num2str(col(1)) ',' num2str(col(2)) ',' num2str(col(3)) ')'];

            case 'flat'

                for n = 1:length(patch_data.CData)

                    capCD = max(min(patch_data.CData(1,n),axis_data.CLim(2)),axis_data.CLim(1));
                    scalefactor = (capCD - axis_data.CLim(1))/diff(axis_data.CLim);
                    col =  255*(colormap(1 + floor(scalefactor*(length(colormap)-1)),:));

                    markercolor{n} = ['rgb(' num2str(col(1)) ',' num2str(col(2)) ',' num2str(col(3)) ')'];

                end
        end
    end

    marker.color = markercolor;

end

%-------------------------------------------------------------------------%

%-MARKER LINE COLOR-%

% marker edge color
MarkerLineColor = patch_data.MarkerEdgeColor;

filledMarker = ismember(patch_data.Marker,filledMarkerSet);

if isnumeric(MarkerLineColor)
    col = 255*MarkerLineColor;
    markerlinecolor = ['rgb(' num2str(col(1)) ',' num2str(col(2)) ',' num2str(col(3)) ')'];
else
    switch MarkerLineColor

        case 'none'

            markerlinecolor = 'rgba(0,0,0,0)';

        case 'auto'

            EdgeColor = patch_data.EdgeColor;

            if isnumeric(EdgeColor)
                col = 255*EdgeColor;
                markerlinecolor = ['rgb(' num2str(col(1)) ',' num2str(col(2)) ',' num2str(col(3)) ')'];
            else

                switch EdgeColor

                    case 'none'

                        markerlinecolor = 'rgba(0,0,0,0)';

                    case {'flat', 'interp'}

                        for n = 1:length(patch_data.CData)

                            capCD = max(min(patch_data.CData(1,n),axis_data.CLim(2)),axis_data.CLim(1));
                            scalefactor = (capCD - axis_data.CLim(1))/diff(axis_data.CLim);
                            col =  255*(colormap(1 + floor(scalefactor*(length(colormap)-1)),:));

                            markerlinecolor{n} = ['rgb(' num2str(col(1)) ',' num2str(col(2)) ',' num2str(col(3)) ')'];

                        end

                end
            end

        case 'flat'

            for n = 1:length(patch_data.CData)

                capCD = max(min(patch_data.CData(1,n),axis_data.CLim(2)),axis_data.CLim(1));
                scalefactor = (capCD - axis_data.CLim(1))/diff(axis_data.CLim);
                col =  255*(colormap(1+floor(scalefactor*(length(colormap)-1)),:));

                markerlinecolor{n} = ['rgb(' num2str(col(1)) ',' num2str(col(2)) ',' num2str(col(3)) ')'];

            end
    end
end

if filledMarker
    marker.line.color = markerlinecolor;
else
    marker.color = markerlinecolor;
end

%-------------------------------------------------------------------------%

end

author: goGPS-Project


text:
function value = r8_degrees ( radians )

%*****************************************************************************80
%
%% R8_DEGREES converts an angle from radian to degree measure.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real RADIANS, the angle measurement in radians.
%
%    Output, real VALUE, the angle measurement in degrees.
%
  value = radians * 180.0 / pi;

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% hessian.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function h = mcs_hessian(i,k,x,x0,f,f0,g,G)
% computes the element G(i,k) of the Hessian of the local quadratic
% model
% Input:
% i, k		indices (k < i)
% x(1:n)	'neighbor' of x0 used for computing G(i,k); differs from
%		x0 only in the ith and kth component
% x0(1:n)	point around which the quadratic model is computed
% f		function value at x
% f0		function value at x0
% g(1:i)	components of the gradient of the local quadratic model %		that have already been computed
% G		components of the Hessian of the local quadratic model
%		that have already been computed
% Output:
% h = G(i,k)	newly computed nondiagonal element of the Hessian

function h = mcs_hessian(i,k,x,x0,f,f0,g,G)
h = f-f0-g(i)*(x(i)-x0(i))-g(k)*(x(k)-x0(k))-0.5*G(i,i)*(x(i)-x0(i))^2-0.5*G(k,k)*(x(k)-x0(k))^2;
h = h/(x(i)-x0(i))/(x(k)-x0(k));

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function p08_default ( dim_num )

%*****************************************************************************80
%
%% P08_DEFAULT sets default values for problem 08.
%
%  Discussion:
%
%    If a problem uses vector parameters, then the spatial dimension
%    DIM_NUM is needed as input, so that the vector parameters can be
%    properly dimensioned.
%
%    Every problem keeps a count of the number of function calls.  Calling
%    this routine causes that count to be reset to 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the problem.
%
  i4 = 0;
  p08_i4 ( 'D', '*', i4 );

  return
end


text:
%%
% Visual display of convolutive kernels.

addpath('../toolbox/');
addpath('../toolbox/img/');
addpath('../toolbox/export_fig/');

rep = 'results/rkhs/';
[~,~] = mkdir(rep);

normalize = @(x)x/sum(x(:));

if not(exist('kernel'))
kernel = 'gaussian-medium';
kernel = 'gaussian-large';
kernel = 'gaussian-small';
kernel = 'energy-dist';
end

% grid points for display
q = 512;
% with padding
q1 = q*8;

% input

rand('state', 1243);
N = 20;
mu = zeros(q,q);
I = randperm(q*q);
mu(I(1:N)) = 1/N;

nu = zeros(q,q);
I = randperm(q*q);
nu(I(1:N)) = 1/N;

mu = load_image('shape-1',q);
nu = load_image('shape-2',q);

mu = normalize(rescale(-mu));
nu = normalize(rescale(-nu));

xi = mu-nu;

% compute sqrt kernel
u = 2 * [0:q1/2, -q1/2+1:-1]' / q1;
[Y,X] = meshgrid(u,u);
D = sqrt(X.^2+Y.^2);

% soft max
s=1/5;
sm = @(t)s*log( 1+exp(t/s) );

switch kernel
    case 'energy-dist'
        r = .3;
        K = max(r-D,0);
        K = sm(r-D);
    case 'gaussian-small'
        sigma = .005;
        K = exp( -D.^2./(2*sigma^2) );
    case 'gaussian-medium'
        sigma = .02;
        K = exp( -D.^2./(2*sigma^2) );
    case 'gaussian-large'
        sigma = .05;
        K = exp( -D.^2./(2*sigma^2) );
end

Kh = fft2(K);
% Kh = max(real(Kh),0);

% for display
K1 = real( ifft2(sqrt(Kh)) );
K1 = fftshift(K1);
K1 = K1(end/2-q/2:end/2+q/2-1, end/2-q/2:end/2+q/2-1);

% zero padding
U = zeros(q1);
U(1:q,1:q) = xi;
xiC = real( ifft2( fft2(U) .* sqrt(Kh) ) );
xiC = xiC(1:q,1:q);

mu1 = mu/max(mu(:));
nu1 = nu/max(nu(:));
imagesc( cat(3, 2-nu1, 2-mu1-nu1, 2-mu1 )/2 );
axis image; axis off;
colormap jet(256);
% saveas(gcf, [rep 'input.png'], 'png');
export_fig([rep 'input.png'], '-m3');

disp_distrib(xiC);
axis image; axis off;
% saveas(gcf, [rep kernel '.png'], 'png');
export_fig([rep kernel '.png'], '-m3');

if strcmp(kernel, 'energy-dist')
    t = linspace(-1,1,q);
    [y,x] = meshgrid(t,t);
    r = .5*1e-3;
    A = -1 ./ sqrt(r + x.^2+y.^2); A = A/max(abs(A(:)));
    disp_distrib( A, 40 );
else
    disp_distrib(-K1, 40);
end
axis image; axis off;
% saveas(gcf, [rep kernel '-kernel.png'], 'png');
export_fig([rep kernel '-kernel.png'], '-m3');

author: optimaltransport


text:
function a = r8mat_test ( trans, lda, m, n )

%*****************************************************************************80
%
%  Purpose:
%
%    R8MAT_TEST sets up a test matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 February 2014
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, character TRANS, indicates whether matrix is to be
%    transposed.
%    'N', no transpose.
%    'T', transpose the matrix.
%
%    Input, integer LDA, the leading dimension of the matrix.
%
%    Input, integer M, N, the number of rows and columns of
%    the matrix.
%
%    Output, real A(LDA,*), the matrix.
%    if TRANS is 'N', then the matrix is stored in LDA*N entries,
%    as an M x N matrix;
%    if TRANS is 'T', then the matrix is stored in LDA*M entries,
%    as an N x M matrix.
%
  if ( trans == 'N' )

    a = zeros ( lda, n );

    for j = 1 : n
      for i = 1 : m
        a(i,j) = 10 * i + j;
      end
    end

  else

    a = zeros ( lda, m );

    for j = 1 : n
      for i = 1 : m
        a(j,i) = 10 * i + j;
      end
    end

  end

  return
end


text:
function value = r8_exponential_sample ( lambda )

%*****************************************************************************80
%
%% R8_EXPONENTIAL_SAMPLE samples the exponential PDF.
%
%  Discussion:
%
%    Note that the parameter LAMBDA is a multiplier.  In some formulations,
%    it is used as a divisor instead.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real LAMBDA, the parameter of the PDF.
%
%    Output, real VALUE, a sample of the PDF.
%
  r = r8_uniform_01_sample ( );

  value = - log ( r ) * lambda;

  return
end


text:
% Add baseline measure as Data Quality measure
% part of the ERPLAB Toolbox
% axs May 2019
%
% Format:
%  baseline_measure = dq_baseline(ERP);
%    or baseline_measure = dq_baseline(ERP,start_time, end_time, subtract_mean_flag)
%
% INPUT:      * - mandatory
%  * ERP  - an ERP structure
%
function [dq_struct, baseline_measure] = dq_baseline(ERP, start_ms, end_ms, subtract_mean_flag)

% Check input
try
    assert(isfield(ERP,'bindata'))
catch
    warning('Making Data Quality baseline measures requires an ERPSET')
    beep
    return
end

% Populate any empty args with defaults
if exist('start_ms','var') == 0 || isempty(start_ms)
    start_ms = ERP.times(1);
end
if exist('end_ms','var') == 0 || isempty(end_ms)
    end_ms = 0;
end
if exist('subtract_mean_flag','var') == 0 || isempty(subtract_mean_flag)
    subtract_mean_flag = 1;
end

% Check desired times are within availible times
good_start = start_ms >= ERP.times;
good_end = end_ms <= ERP.times;
try
    assert(any(good_start)==1)
    assert(any(good_end)==1)
catch
    beep
    warning('The times requested for Data Quality baseline are not available in this ERPset.')
    return
end

start_dp = find(abs(ERP.times-start_ms)==min(abs(ERP.times-start_ms))); % datapoint closest to desired ms
start_dp = start_dp(1); % if it's a tie, take the first element
end_dp = find(abs(ERP.times-end_ms)==min(abs(ERP.times-end_ms)));
end_dp = end_dp(1);

%

b_sd_1 =   std(ERP.bindata(1,start_dp:end_dp,1));

baseline_data = ERP.bindata(:,start_dp:end_dp,:);

if subtract_mean_flag == 1
    % Measure Standard Deviation, which removes mean
    baseline_measure = std(baseline_data,0,2);

elseif subtract_mean_flag == 0
    % use rms, not removing mean
    baseline_measure = rms(baseline_data,2);

elseif subtract_mean_flag == 2
    %Corrected SD
    nTimes = size(baseline_data,2);
    baseline_measure_biased = std(baseline_data,0,2);
    baseline_measure =  ((baseline_measure_biased * sqrt((nTimes-1)))/ sqrt((nTimes-(3/2)+(1/(8*(nTimes-1))))));

end

% quick plot of baseline sd
%imagesc(squeeze(baseline_measure))

% Prepare dq_struct
if subtract_mean_flag == 1
    dq_struct.type = 'Baseline Measure - SD';
elseif subtract_mean_flag == 0
    dq_struct.type = 'Baseline Measure - RMS';
elseif subtract_mean_flag == 2
    dq_struct.type = 'Baseline Measure - SD (Corrected)';
end

dq_struct.times = [start_ms end_ms];
dq_struct.data = baseline_measure;

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
load -ascii house.dat
houseD=house';
clear house

[N, m] = size(houseD)

class = 1

%rand('state',0); randn('state',0);
%houseD = houseD(:,randperm(m));

Napp = ceil(m*2/3);
Ntest = m-Napp

app  = houseD(:,1:Napp);size(app)
test = houseD(:,Napp+1:end);size(test)

unique(app(class,:))
unique(test(class,:))

ns = max(houseD')
clear houseD

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function pick = nms(boxes, overlap)
% top = nms(boxes, overlap)
% Non-maximum suppression. (FAST VERSION)
% Greedily select high-scoring detections and skip detections
% that are significantly covered by a previously selected
% detection.
%
% NOTE: This is adapted from Pedro Felzenszwalb's version (nms.m),
% but an inner loop has been eliminated to significantly speed it
% up in the case of a large number of boxes

% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

if isempty(boxes)
  pick = [];
  return;
end

x1 = boxes(:,1);
y1 = boxes(:,2);
x2 = boxes(:,3);
y2 = boxes(:,4);
if size(boxes,2)==4
  s = ones(1,size(boxes,1));
else
  s = boxes(:,end);
end

area = (x2-x1+1) .* (y2-y1+1);
[~, I] = sort(s);

pick = s*0;
counter = 1;
while ~isempty(I)
  last = length(I);
  i = I(last);
  pick(counter) = i;
  counter = counter + 1;

  xx1 = max(x1(i), x1(I(1:last-1)));
  yy1 = max(y1(i), y1(I(1:last-1)));
  xx2 = min(x2(i), x2(I(1:last-1)));
  yy2 = min(y2(i), y2(I(1:last-1)));

  w = max(0.0, xx2-xx1+1);
  h = max(0.0, yy2-yy1+1);

  inter = w.*h;
  o = inter ./ (area(i) + area(I(1:last-1)) - inter);

%   I = I(find(o<=overlap));
  I = I((o<=overlap));
end

pick = pick(1:(counter-1));

author: hbilen

repo: WSDDN

save_path: github-repos/MATLAB/hbilen-WSDDN


text:
function [ax, data] = scatterPlot(X, YLbls, markerSize, ax);

% SCATTERPLOT 2-D scatter plot of labelled points.
% FORMAT
% DESC plots a 2-D scatter plot of labelled points.
% ARG X : x points to plot.
% ARG lbls : the labels to plot.
% ARG markerSize : the size of the markers to use.
% ARG ax : the axes to plot on (default create new axes).
% RETURN ax : the axis handle of the plot.
% RETURN data : the handles of the data points.
%
% SEEALSO : plot, getSymbols
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2010

% NDLUTIL

  if nargin < 4
    ax = [];
  end
  if isempty(YLbls)
    symbol = [];
  else
    symbol = getSymbols(size(YLbls,2));
  end

  % Create the plot for the data
  clf
  if isempty(ax)
    ax = axes('position', [0.05 0.05 0.9 0.9]);
  else
    ax = axes(ax);
  end
  hold on

  data = twoDPlot(X, YLbls, symbol);
  for i = 1:2
    minX = min(X(:, i));
    maxX = max(X(:, i));
    xSpan = maxX - minX;
    xLim{i}(1) = minX - 0.1*xSpan;
    xLim{i}(2) = maxX + 0.1*xSpan;
  end
  if nargin>2
    for i = 1:length(data)
      set(data(i), 'markersize', markerSize);
    end
  end
  set(ax, 'xLim', xLim{1});
  set(ax, 'yLim', xLim{2});

  set(ax, 'fontname', 'arial');
  set(ax, 'fontsize', 20);
end

function returnVal = twoDPlot(X, label, symbol)

% GPLVMTWODPLOT Helper function for plotting the labels in 2-D.

  returnVal = [];

  if ~isempty(label)
    for i = 1:size(X, 1)
      labelNo = find(label(i, :));
      try
        returnVal = [returnVal; plot(X(i, 1), X(i, 2), symbol{labelNo})];
      catch
        if strcmp(lasterr, 'Index exceeds matrix dimensions.')
          error(['Only ' num2str(length(symbol)) ' labels supported (it''s easy to add more!)'])
        end
      end
    end
  else
    returnVal = plot(X(:, 1), X(:, 2), 'rx');
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function A = sdt_A(H, F)
% Calculate signal detection theory A statistic:
% a non-parametric estimate of sensitivity in ROC analysis
% per Zhang and Mueller (2005)
%
% :Usage:
% ::
%
%     A = sdt_A(H, F)
%
% :Inputs:
%
%   **H:**
%        hit rate
%
%   **F:**
%        false alarm rate
%
% :Outputs:
%
%   **A:**
%        A measure of sensitivity of a detector, essentially the area under
%        the ROC curve for either "normal" (H>=F) or "reverse skill" (F>H)
%        case (ranges from 0.5 to 1).
%
% Meaning of values where F > H ("reverse skill"):
%
% The meaningful measure of sensitivity in this case is A(F,H)
% as opposed to A(H,F).  Subject is assumed to still be sensitive,
% but to have inverted their interpretation (or "reverse skill").
%
% To distinguish from the normal case, such results are returned
% as (1 - sdt_A(F,H)), so the full range from 0 to 1 is used.
% For comparisons of sensitivity use abs(sdt_A(H,F) - 0.5).
%
% :Example:
%
% sdt_A(0.6,0.8) = 0.325. This indicates the same sensitivity as
% sdt_A(0.8,0.6) = 0.675, but with reversed skill.
%
% ..
%    Author: Joe Wielgosz 8/13/2009
% ..

if F <= H % Normal case

    if (F <= 0.5) && (0.5 <= H)
        A = 0.75 + (H-F)/4 - F*(1-H);
    elseif (F <= H) && (H <= 0.5)
        A = 0.75 + (H-F)/4 - F/(4*H);
    else % (0.5 < F) && (F <= H)
        A = 0.75 + (H-F)/4 - (1-H)/(4*(1-F));
    end

else % "Reverse skill" - more false positives than hits
    % Using Jason Buhle's approach
    A = 1 - sdt_A(F, H);
%    A = NaN;
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%  Figure 6.58      Feedback Control of Dynamic Systems, 5e
%                   Franklin, Powell, Emami
%

clear all
close all

num=9;
den=conv([1 0.5],[1 1]);
den=conv(den,[1 2]);
w=logspace(-1,1,500);
[mag,phas]=bode(num,den,w);
[OLGM,OLPM,OLWcg,OLWcp]=margin(mag,phas,w)

%Lead compensator, first guess
numl=[1 1];
denl=0.333*[1 3];
numc=conv(num,numl);
denc=conv(den,denl);
[magc,phasc]=bode(numc,denc,w);
[D1GM,D1PM,D1Wcg,D1Wcp]=margin(magc,phasc,w)
dencl=denc+[0 0 0 numc];
t=0:.1:20;
y=step(numc,dencl,t);
%Design Iteration, next guess
numl=[1 1.5];
denl=0.1*[1 15];
numcc=conv(num,numl);
dencc=conv(den,denl);
[magcc,phascc]=bode(numcc,dencc,w);
[D2GM,D2PM,D2Wcg,D2Wcp]=margin(magcc,phascc,w)
subplot(2,1,1)
loglog(w,mag,'-',w,magc,'--',w,magcc,'-.',w,ones(500,1),'-');
grid;
%xlabel('w (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.58 Bode Plot for lead-compensation design (a) magnitude');
subplot(2,1,2)
semilogx(w,phas,'-',w,phasc,'--',w,phascc,'-.',w,-180*ones(500,1));
grid;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
title('Fig. 6.58 (b) phase');


text:
function C = DiagCorr(a,b)
% e.g. dim of A: nFrames*nCells, like in Matlab 'corr' function
An=bsxfun(@minus,a,mean(a,1));
Bn=bsxfun(@minus,b,mean(b,1));
An=bsxfun(@times,An,1./sqrt(sum(An.^2,1)));
Bn=bsxfun(@times,Bn,1./sqrt(sum(Bn.^2,1)));
C=sum(An.*Bn,1);
end

author: xiuyechen

repo: FishExplorer


text:
%EX_TUMBLE Generate synthetic sensor data for a tumbling body
%
% Sets a number of workspace variables:
% t         time vector (1xN)
% w         angular velocity (3xN)
% truth     actual orientation as a unit quaternion (1xN)
% wm        simulated gyro output with bias (3xN)
% am        simulated accelerometer output with bias (3xN)
% mm        simulated magnetometer output with bias (3xN)
%
% Used on pages 81 and 89 of Robotics, Vision & Control, 2nd ed.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%% parameters (user adjustable)

dt = 0.05;          % sample interval

% make an asymmetric mass
J = diag([2 4 3]);
J(1,2) = -1;
J(2,1) = -1;
J(1,3) = -2;
J(3,1) = -2;
%eig(J)

% accelerometer
g0 = unit( [0, 0, 9.8]' );
gbias = 0.02*[2 -2 2]';  % bias 2% of norm

% magnetometer, use N E U data in nT  for Brisbane
m0 = unit( [28067.5, -5439.4, 44800.5]'*1e-9 );
mbias = 0.02*[-1 -1 2]';   % bias 2% of norm

% gyro
w0 = 0.2*[1 2 2]';
wbias = 0.1*[-1 2 -1]'; % bias 10% of max

%% simulation

% Solve Euler's rotational dynamic equation to get omega
[t,w] = ode45( @(t,w) -inv(J)*(cross(w, J*w)), [0:dt:20], w0);
w = w';  % one column per timestep

% Compute simulated sensor readings and true attitude
gm = zeros(3, numcols(w));
mm = zeros(3, numcols(w));

truth(1) = UnitQuaternion();

for k=1:numcols(w)-1
        iq = inv(truth(k));
        gm(:,k) = iq * g0 + gbias;  % sensor reading in body frame
        mm(:,k) = iq * m0 + mbias;  % sensor reading
        truth(k+1) = truth(k) .* UnitQuaternion.omega( w(:,k)*dt );
end

% add bias to measured angular velocity
wm = bsxfun(@plus, w, wbias);

author: petercorke

repo: robotics-toolbox-matlab


text:
function [minscale,crit,delta_list,lgd] = compute_minimum_scale(D, options)

% compute_minimum_scale - compute minimum deconvolution scale
%
%   D should be a (p,n) dictionary matrix. Each D(:,i) is an atom located
%   at some position i.
%
%   mscale(i) is a minimum scale as computed using a given criterion:
%       i=1 WERC criterion
%       i=2 ERC criterion
%       i=3 Fuchs criterion
%   one has mscale(i)>mscale(i+1) since the Fuch criterion is the finest
%   criterion (depends on sign).
%
%   options.mscale_type can be 'train' (spike train) or 'twodiracs'.
%   options.display=1 to display a graph of the criterions
%
%   Copyright (c) 2008 Gabriel Peyre

p = size(D,2); n = size(D,1);

options.null = 0;
mscale_type = getoptions(options, 'mscale_type', 'train');
delta_max = getoptions(options, 'delta_max', p/4);
disp = getoptions(options, 'display', 0);
delta_max = min(delta_max,p/2);
verb = getoptions(options, 'verb', 1);

subsampling = getoptions(options, 'subsampling', 1);

delta_list = 2:delta_max;
ntests = length(delta_list);

% record normalized gram
d = sqrt(sum(D.^2));
D1 = D ./ repmat(d, [size(D,1),1]);
G = abs(D1'*D1);

erc = [];
werc = [];
fuchs = [];
for i=1:ntests
    if verb
        progressbar(i,ntests);
    end
    delta = delta_list(i);
    % generate signal with spacing delta
    x = zeros(p,1);
    switch mscale_type
        case 'train'
            x(1:delta:p-delta+1) = 1;
        case 'twodiracs'
            x(round(end/2)) = 1;
            x(round(end/2)+delta) = 1;
        otherwise
            error('Unknown type');
    end
    % compute criterias
    werc(i) = compute_werc_criterion(D1,x,G);
    erc(i) = 1 - compute_erc_criterion(D1,x);
    fuchs(i) = compute_fuchs_criterion(D1,x);
end

crit = [werc; erc; fuchs];
lgd = {'werc', 'erc', 'fuchs'};
col = {'b', 'g', 'r'};
vmin = .7; vmax = 1.8;

if disp
    clf;
    hold on;
    aa = plot(delta_list*subsampling, crit); axis tight;
    set(aa, 'LineWidth', 2);
end
for i=1:3
    k = ntests;
    while crit(i,k)<1 && k>0
        k = k-1;
    end
    minscale(i) = delta_list(min(k+1,end))*subsampling;
    if disp
        aa = plot([1 1]*minscale(i), [vmin vmax], [col{i}(1) ':']);
        set(aa, 'LineWidth', 2);
    end
end
if disp
    axis([0 max(delta_list)*subsampling vmin vmax]);
    hold off; box on;
    legend(lgd);
end

author: gpeyre

repo: matlab-toolboxes


text:
function [out,pos] = max( in, plonk, dim )
%MAX for nested cells
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

% return scalar value
[out, ~] = flattenCellMatrix(in.data);
if(nargin == 1)
    out = max(cellfun(@(x) max(x(:)), out));
else
    [out, pos] = cellfun(@(x) max(x,[],dim), out, 'UniformOutput', false);
end

end

author: thomaskuestner


text:
function model = ml_traincov(varargin)
% Learn a linear predictive model using covariance-based classification.
% Model = ml_traincov(Trials, Targets, Options...)
%
% This method assumes that the given trials values represent covariance matrices and offers various
% methods, including information geometric approaches, to classify them. This implementation uses
% the covariance toolbox by Alexandre Barachant.
%
% In:
%   Trials       : training data, as in ml_train
%
%   Targets      : target variable, as in ml_train
%
% Out:
%   Model   : the predictive model; can be used with ml_predictcov
%
% Examples:
%   % learn a model using the defaults
%   model = ml_traincov(trials,targets);
%
% See also:
%   ml_predictcov
%
%                           Christian Kothe, Syntrogi
%                           2015-04-21
dp;

arg_define([0 3],varargin, ...
    arg_norep('trials'), ...
    arg_norep('targets'), ...
    arg({'classifier','ClassificationMethod'},'mdm',{'mdm','fgmdm','tslda'}, 'Classification method to use. The mdm method (minimum distance to mean) is based on minimum distance to class mean using a given distance metric and mean estimation metric. The fgmdm (filtered geodesic mdm) method additionally performs geodesic filtering. The tslda method (tangent-space linear discriminant analysis) is an LDA implementation that respects the Riemannian geometry of the data.'), ...
    arg({'mean_est','MeanEstimator'},'riemann',{'arithmetic','riemann','riemanndiag','riemanntrim','median','riemannmed','logeuclid','opttransp','ld','geodesic','harmonic','geometric'},'Method to average covariance matrices. Various methods are supported, including Riemannian mean/median/trimmed mean, log-euclidean mean, optimal transportation mean, log determinant mean, and geodesic iterative mean.'), ...
    arg({'distance_metric','DistanceMetric'},'riemann',{'euclid','riemann','kullback','logeuclid','opttransp','ld'},'Distance metric to use. This is only used for the mdm and fgmdm methods. Different distance metrics are supported, including the Euclidean metric, the Riemannian distance, Kullback-Leibler divergence, log-euclidean distance, optimal transportation distance, and log determinant distance'));

% find the class labels
classes = unique(targets);
if length(classes) == 1
    error('BCILAB:only_one_class','Your training data set has no trials for one of your classes; you need at least two classes to train a classifier.\n\nThe most likely reasons are that one of your target markers does not occur in the data, or that all your trials of a particular class are concentrated in a single short segment of your data (10 or 20 percent). The latter would be a problem with the experiment design.');
else
    % reformat feature shape to DxDxN
    if ndims(trials) == 2 %#ok<ISMAT,NODEF>
        [N,F] = size(trials);
        D = sqrt(F);
        if abs(D-round(D)) > 0
            error('The number of features in your trials must be a square of an integer.'); end
        trials = reshape(trials',[D,D,N]);
    else
        [U,V,N] = size(trials); %#ok<ASGLU>
        if U ~= V
            error('Your feature matrices are not square, i.e., cannot be covariance matrices.'); end
    end

    model.classes = classes;
    switch classifier
        case 'mdm'
            % minimum distance to mean: estimate class means
            for c=length(classes):-1:1
                model.C{c} = mean_covariances(trials(:,:,targets==classes(c)),mean_est); end
        case 'fgmdm'
            % filtered geodesic minimum distance to mean
            % geodesic filtering
            [model.W,model.Cg] = fgda(trials,targets,mean_est,{},'shcov',{});
            trials = geodesic_filter(trials,model.Cg,model.W(:,1:length(model.classes)-1));
            % estimate class means
            for c=length(classes):-1:1
                model.C{c} = mean_covariances(trials(:,:,targets==classes(c)),mean_est); end
        case 'tslda'
            error('This variant is not yet implemented!');
        otherwise
            error('Unsupported classification method: %s',hlp_tostring(classifier));
    end
end

model.classifier = classifier;
model.mean_est = mean_est;
model.distance_metric = distance_metric;

repo: REST