text:
function events = in_events_ricoh(sFile, EventFile)
% IN_EVENTS_RICOH: Read the events descriptions for a RICOH MEG file.
%
% USAGE:  events = in_events_ricoh(sFile, EventFile)
%
% This function is based on the Ricoh MEG reader toolbox version 1.0.
% For copyright and license information and software documentation,
% please refer to the contents of the folder brainstorm3/external/ricoh

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2018

% Read file info using Ricoh functions
header.events     = getRHdrEvent(EventFile);       % Get information about trigger events.
header.annotation = getRHdrAnnotation(EventFile);  % Get information about annotations.

% Initialize returned structure
events = repmat(db_template('event'), 0);
% Triggers
if ~isempty(header.events)
    % All all the events types
    uniqueNames = unique({header.events.name});
    % Create events structures: one per category of event
    for i = 1:length(uniqueNames)
        % Add a new event category
        iEvt = length(events) + 1;
        % Find all the occurrences of event #iEvt
        iMrk = find(strcmpi({header.events.name}, uniqueNames{i}));
        % Get all samples
        allSamples = [header.events(iMrk).sample_no];
        % Get the epoch numbers (only ones for continuous and averaged files)
        if (sFile.header.acq.acq_type == 1) || (sFile.header.acq.acq_type == 2)
            iEpochs = ones(1, length(iMrk));
        else
            iEpochs = floor(allSamples / sFile.header.acq.frame_length) + 1;
            allSamples = allSamples - (iEpochs-1) .* sFile.header.acq.frame_length;
        end
        % Add event structure
        events(iEvt).label      = uniqueNames{i};
        events(iEvt).epochs     = iEpochs;
        events(iEvt).times      = allSamples ./ sFile.prop.sfreq;
        events(iEvt).reactTimes = [];
        events(iEvt).select     = 1;
        events(iEvt).channels   = [];
        events(iEvt).notes      = [];
    end
end
% Annotations
if ~isempty(header.annotation)
    % Get annotation label
    annotLabel = cell(1,length(header.annotation));
    for i = 1:length(header.annotation)
        if ~isempty(strtrim(header.annotation(i).comment))
            annotLabel{i} = strtrim(header.annotation(i).comment);
            annotLabel{i}(annotLabel{i} == 0) = [];
        elseif ~isempty(header.annotation(i).label) && isnumeric(header.annotation(i).label)
            annotLabel{i} = ['annot_', num2str(header.annotation(i).label)];
        else
            annotLabel{i} = 'Unknown';
        end
    end

    % All all the events types
    uniqueNames = unique(annotLabel);
    % Create events structures: one per category of event
    for i = 1:length(uniqueNames)
        % Add a new event category
        iEvt = length(events) + 1;
        % Find all the occurrences of event #iEvt
        iMrk = find(strcmpi(annotLabel, uniqueNames{i}));
        % Get all samples
        allSamples = [header.annotation(iMrk).sample_no];
        % Get the epoch numbers (only ones for continuous and averaged files)
        if (sFile.header.acq.acq_type == 1) || (sFile.header.acq.acq_type == 2)
            iEpochs = ones(1, length(iMrk));
        else
            iEpochs = floor(allSamples / sFile.header.acq.frame_length) + 1;
            allSamples = allSamples - (iEpochs-1) .* sFile.header.acq.frame_length;
        end
        % Add event structure
        events(iEvt).label      = uniqueNames{i};
        events(iEvt).epochs     = iEpochs;
        events(iEvt).times      = allSamples ./ sFile.prop.sfreq;
        events(iEvt).reactTimes = [];
        events(iEvt).select     = 1;
        events(iEvt).channels   = [];
        events(iEvt).notes      = [];
    end
end

author: brainstorm-tools


text:
function best_weights = GetBestMixingWeight(stand_dev_of_variables)
% Best mixing of random variables with common mean
% and built-in constraint
%   best_weights = GetBestMixingWeight(stand_dev_of_variables)
%
% Input argument
%   stand_dev_of_variables: standard deviation of each variable
%
% Output argument
%   best_weights: best set of weights for mixing random variables

% Copyright 2016 Google Inc. All Rights Reserved.
% Author: hidekik@google.com (Hideki Kawahara)

narginchk(1, 1);
n_variables = length(stand_dev_of_variables);
H = ones(n_variables - 1, n_variables - 1) * ...
  stand_dev_of_variables(n_variables) ^ 2;
H = H + diag(stand_dev_of_variables(1 : n_variables - 1) .^ 2);
v = ones(n_variables - 1, 1) * stand_dev_of_variables(n_variables) ^ 2;
a = H \ v;
best_weights = [a; 1 - sum(a)];
end

author: google


text:
function a = lehmer ( n )

%*****************************************************************************80
%
%% LEHMER returns the Lehmer matrix.
%
%  Discussion:
%
%    The Lehmer matrix is a symmetric positive definite N by N
%    matrix with
%
%      A(i,j) = min ( i, j ) / max ( i, j )
%
%    A is totally nonnegative.  The inverse of A is tridiagonal,
%    and explicit formulas are known for its entries.
%
%    The condition number of A satisfies the bounds:
%
%      N <= condition ( A ) <= 4*N*N.
%
%  Reference:
%
%    M. Newman and J. Todd,
%    The evaluation of matrix inversion programs,
%    Journal of the Society for Industrial and Applied Mathematics,
%    Volume 6, 1958, pages 466-476.
%
%    Solutions to problem E710 (proposed by D H Lehmer):
%    The inverse of a matrix,
%    American Mathematical Monthly,
%    Volume 53, 1946, pages 534-535.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(N,N), the matrix.
%
  a = ones ( n, 1 ) * ( 1:n );
  a = a ./ a';
  a = tril ( a ) + tril ( a, -1 )';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%test_PID
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR 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 General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

figure

colormap jet

s = 12;

PIDM = getPID(s, 1, planC, 'transverse', 'max');

subplot(1,3,1), imagesc(PIDM), axis image, title('transverse')

PIDM = getPID(s, 1, planC, 'saggital', 'max');

subplot(1,3,2), imagesc(PIDM), axis image, title('saggital')

PIDM = getPID(s, 1, planC, 'coronal', 'max');

subplot(1,3,3), imagesc(PIDM), axis image, colorbar, title('coronal')

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function aps_powerlaw(start_step,end_step,save_path)
% aps_powerlaw(start_step,end_step)
% Program to compute the powerlaw tropsoheric delays.
%   STEP 1 = Estimate powerlaw delay and maximum height from sounding data
%            (when this is available), or use defaults or user set values
%   STEP 2 = Rotate dataset to minimize interpolation effects, Powerlaw scaling and interpolate to regular grid
%   STEP 3 = 1D/2D bandfiltering in the frequency domain and interpolate
%            back to local grid
%   STEP 4 = Spatial local estimation of the powerlaw scaling coefficient,
%            estimate final values from reliable frequency band, extrapolate
%            to all point locations and compute powerlaw tropospheric delay.
%
%   To change the processing parameters use setparm_aps and getparm_aps
%
% **** when using this estimation method please cite:
%       D.P.S. Bekaert, A.J. Hooper and T.J. Wright, A spatially-variable power-law tropospheric correction
% 	  	 technique for InSAR data, JGR, doi:10.1029/2014JB011558 *****
%
%
%     Copyright (C) 2015  Bekaert David - University of Leeds
%     Email: eedpsb@leeds.ac.uk or davidbekaert.com
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program 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 General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% By David Bekaert - University of Leeds
%
% modifications:
% 04/2013   DB:     Incorportate a parm_aps list with all the processing
%                   options and variables.
% 05/2013   DB:     Include powerlaw coefficient sensitivity analysis
% 05/2013   DB:     Allow for an interferogram based estimation
% 12/2013   DB:     Put the flag of ifg correction on prior to powerlaw computation
% 02/2014   DB:     Integrate the plotting with the running of the sounding data
% 02/2014   DB:     Include a mountain ridge code, user interaction needed
% 07/2014   DB:     Have the option to plot the ridges
% 07/2014   DB:     Include a cropping option
% 11/2014   DB:     Fix bug and include paper reference
% 03/2015   DB:     Ask for the reference technique in step 5
% 01/2016   DB:     Fix ll.lonlat, ph_uw, hgt for non-stamps data, replace
%                   save calls with aps_save, in case no -append flag is included
% 02/2016   DB:     Expand option 5 to plot with the unwrapped phase too.

if nargin<3
   save_path = [pwd filesep];
end

currdir = pwd;		% current working directory
plot_flag = 0;
hydro =1;
wet=1;

fprintf(['\n********************************************************************************************************** \n','D.P.S. Bekaert, A.J. Hooper and T.J. Wright,\n', 'A spatially-variable power-law tropospheric correction technique for InSAR data, JGR,  doi:10.1029/2014JB011558\n', '**********************************************************************************************************\n\n'])

%% PART 1: Define the powerlaw coefficents, from the user or by using sounding data.
sounding_data = getparm_aps('sounding_data',1);
stamps_processed = getparm_aps('stamps_processed',1);
powerlaw_ridge_constraint = getparm_aps('powerlaw_ridge_constraint',1);

if strcmp(stamps_processed,'y')
    load psver
else
    psver = 2;
end

% file names of the output data
apsname = [save_path filesep 'tca' num2str(psver) '.mat'];
apssbname = [save_path filesep 'tca_sb' num2str(psver) '.mat'];
apsbandsname = [save_path filesep 'tca_bands' num2str(psver) '.mat'];
apsbandssbname = [save_path filesep 'tca_bands_sb' num2str(psver) '.mat'];

% saving part of the data in a subfolder
save_path = [save_path filesep 'aps_p'];
if exist(save_path,'dir')~=7
    mkdir(save_path);
end

if start_step==0 && strcmp(powerlaw_ridge_constraint,'y')

    mountain_ridge=[];
    if exist('tca_support.mat')==2
        load('tca_support.mat','powerlaw_ridges')
        if exist('powerlaw_ridges','var')==1
            mountain_ridge = powerlaw_ridges.mountain_ridge;
        end
    end
    if isempty(mountain_ridge)==1
        str_repr='y';
        str_vis = 'n';
    end

    % data has been computed before - reporcess? - visualize?
    if length(mountain_ridge)>0
        str_repr = '';
        while strcmpi(str_repr,'y')~=1 && strcmpi(str_repr,'n')~=1
            str_repr = input(['Ridges have been defined before, reprocess the data? [y/n] \n'],'s');
        end

        % no reprocessing, do you want to visualize it?
        if strcmpi(str_repr,'n')
            str_vis = '';
            while strcmpi(str_vis,'y')~=1 && strcmpi(str_vis,'n')~=1
                str_vis = input(['Do you want to visualize the data? [y/n] \n'],'s');
            end
        else
            str_vis = 'n';
        end

    end

    if strcmpi(str_repr,'y')
        [mountain_ridge] = aps_powerlaw_watershed;
    end
    if strcmpi(str_vis,'y')
        aps_support_plot(1);
    end
end

if start_step==1 && strcmp(sounding_data,'y')
    fprintf('\n\nStep 1: Powerlaw coefficients from the sounding data\n')

    % getting the variables from the parm_aps file
    look_angle = getparm_aps('look_angle',1);
    lambda = getparm_aps('lambda',1);
    sounding_h0= getparm_aps('sounding_h0',1);
    sounding_h_alpha_thres= getparm_aps('sounding_h_alpha_thres',1);
    sounding_start_date= getparm_aps('sounding_start_date',1);
    sounding_end_date= getparm_aps('sounding_end_date',1);
    sounding_time_stamp= getparm_aps('sounding_time_stamp',1);
    sounding_dir= getparm_aps('sounding_dir',1);
    sounding_error_promp= getparm_aps('sounding_error_promp',1);
    sounding_sensitivity= getparm_aps('sounding_sensitivity',1);
    n_months = getparm_aps('sounding_months',1);
    time_stamp = getparm_aps('sounding_time_stamp',1);
    sounding_ifg_dates = getparm_aps('sounding_ifg_dates',1);

    %%% check if the file is already exisiting - ask if reprocessing is
    %%% needed, if not load exisitng file and plot results
    time_stamp_str = [];
    for k=1:size(time_stamp,1)
        if k>1
            time_stamp_str = [time_stamp_str '_' time_stamp(k,:)];
        else
            time_stamp_str = [time_stamp(k,:)];
        end
    end

    if strcmp(sounding_sensitivity,'y')
        fprintf('Sensitivity analyses of sounding data \n\n')

        % check first if the output file is already exisiting otherwize
        % ask if it needs to be reprocessed or visualised.

        if strcmp(sounding_ifg_dates,'y')
            fprintf('Coefficient estimation for each interferogram \n\n')
            if hydro==1 && wet==0
                save_name = [sounding_dir filesep 'Powerlaw' filesep 'Powerlaw_sensitivity_hydro_SAR_dates_1month_' time_stamp_str 'Hr.mat'  ];
            elseif hydro==0 && wet==1
                save_name = [sounding_dir filesep 'Powerlaw' filesep 'Powerlaw_sensitivity_wet_SAR_dates_1month_' time_stamp_str 'Hr.mat'  ];
            else
                save_name = [sounding_dir filesep 'Powerlaw' filesep 'Powerlaw_sensitivity_SAR_dates_1month_' time_stamp_str 'Hr.mat'  ];
            end
            n_months =1;
        else
            fprintf('Coefficient estimation based on average \n\n')
            % putting the variables in the right set-up
            start_year = str2num(sounding_start_date(1:4));
            end_year = str2num(sounding_end_date(1:4));
            start_str = sounding_start_date(5:6);
            end_str = sounding_end_date(5:6);

            % the file name to be loaded
            if hydro==1 && wet==0
                save_name = [sounding_dir filesep 'Powerlaw' filesep 'Powerlaw_sensitivity_hydro_' num2str(n_months) 'month_' time_stamp_str 'Hr_' num2str(start_year) start_str '_' num2str(end_year) end_str '.mat'  ];
            elseif hydro==0 && wet==1
                save_name = [sounding_dir filesep 'Powerlaw' filesep 'Powerlaw_sensitivity_wet_' num2str(n_months) 'month_' time_stamp_str 'Hr_' num2str(start_year) start_str '_' num2str(end_year) end_str '.mat'  ];
            else
                save_name = [sounding_dir filesep 'Powerlaw' filesep 'Powerlaw_sensitivity_' num2str(n_months) 'month_' time_stamp_str 'Hr_' num2str(start_year) start_str '_' num2str(end_year) end_str '.mat'  ];
            end
        end

        if exist(save_name,'file')==2
            str = '';
            while strcmpi(str,'y')~=1 && strcmpi(str,'n')~=1
                str = input(['This file already exist, do you want to re-process the data? [y/n] \n'],'s');
            end

            if strcmpi(str,'y')
                sounding_powerlaw_sens(hydro,wet)
            else
                sounding_powerlaw_sens_display
            end
        else
            % estimate the coefficients for individual interferograms
            sounding_powerlaw_sens(hydro,wet)
        end

    else
        fprintf('Based on specified sounding period \n\n')
        % When having sounding data estimate the powerlaw and b coefficient
        [alpha_log_all,alpha,h0,n_soundings] = sounding;

        % setting the estimated parameters in the parameter file
        fprintf('Updating the parm_aps list \n')
        setparm_aps('powerlaw_h0',h0)
        setparm_aps('powerlaw_alpha',alpha)

        % saving the data
        save([save_path filesep 'aps_p_step1.mat'],'alpha','h0','look_angle','lambda','sounding_h0','sounding_h_alpha_thres','sounding_start_date','sounding_end_date','sounding_time_stamp')

    end
end
if start_step==1 && strcmp(sounding_data,'n')
     fprintf('\n\n Step 1: powerlaw coefficents set by user \n\n')

     % as set by the user
     h0 = getparm_aps('powerlaw_h0');
     alpha = getparm_aps('powerlaw_alpha');
end

%% PART 2: Rotate the dataset to get mimumal white space around it
% Scaling and interpolation of the data to a regular grid
if start_step<=2 && end_step >=2
    fprintf('\n\nStep 2: Rotate the dataset \n')

    % getting the variables from the parm_aps file
    heading = getparm_aps('heading',1);
    crop_flag =  getparm_aps('crop_flag',1);

    % loading the lonlat information
    ll_matfile = getparm_aps('ll_matfile',1);
    ll = load(ll_matfile);
    ll = ll.lonlat;

    % retrieving the coefficients
    h0 = getparm_aps('powerlaw_h0',1);
    alpha = getparm_aps('powerlaw_alpha',1);
    getparm_aps('powerlaw_xy_res',1);
    phuw_matfile = getparm_aps('phuw_matfile',1);
    hgt_matfile = getparm_aps('hgt_matfile',1);
    DEM_corr = getparm_aps('powerlaw_DEM_corr',1);
    bperp_matfile = getparm_aps('bperp_matfile',1);

    % rotating the data
    [xy_rotatelocal_full,xy_local_full,local_origin_full] = ll2rotatelocal(ll,heading);

    % cropping if needed
    crop = [];
    if strcmpi(crop_flag,'y')
        if exist('area_ex.mat','file')==2
            crop = load('area_ex.mat');
        elseif exist('../area_ex.mat','file')==2
            crop = load('../area_ex.mat');
        else
            crop = [];
            fprintf(['There is no area_ex.mat file found, no data is cropped... \n'])
        end
    end
    if ~isempty(crop)
        ix_remove =  inpolygon(ll(:,1),ll(:,2),crop.lonlat(:,1),crop.lonlat(:,2));
    else
        ix_remove = [];
    end
    ll(ix_remove,:) =[];

    % get the InSAR convexhull -- not needed at this stage!
    ix_convexhull = convhull(ll(:,1),ll(:,2));
    InSAR_convexhull = ll(ix_convexhull,:);
    if exist('tca_support.mat','file')==2
      save('tca_support.mat','-append','InSAR_convexhull')
    else
      save('tca_support.mat','InSAR_convexhull')
    end
    clear InSAR_convexhull

    % rotating the data using all data including the region cropped.
    % This will be used in the last stage of the power-law code.
    [xy_rotatelocal,xy_local,local_origin] = ll2rotatelocal(ll,heading);

    fprintf('        Scaling and interpolation to a regular grid \n\n')
    % loading the unwrapped interferograms
    ph = load(phuw_matfile);
    ph = ph.ph_uw;

    % keep track of those pixels which are NaN
    ix_phase_nan = isnan(ph);

    % remove crop if needed
    ph(ix_remove,:)=[];

    % loading the topography information
    hgt = load(hgt_matfile);
    hgt = hgt.hgt;
    % remove crop if needed
    hgt(ix_remove,:)=[];

    % geting the number of interferograms and points from the phase matrix
    n_interferograms= size(ph,2);
    n_points = size(ph,1);

    % Scaling of the topography according to the powerlaw
    % ifg based powerlaw correction?
    if length(alpha)==n_interferograms
        hgt_max = max(hgt);
        hgt_scaled=repmat(hgt,1,size(ph,2));
        ifg_based_correction='y';

        for k=1:length(alpha)
            if h0(k)*1000-hgt_max<=0
                error('myApp:argChk', ['Power law is not valid above h0 heights,... \nAbort,... \n'])
            end
            hgt_scaled(:,k)=(h0(k)*1000-hgt_scaled(:,k)).^(alpha(k));
        end
    else
        % Scaling of the topography according to the powerlaw
        ifg_based_correction='n';
        ix = find(h0*1000-hgt < 0);
        if isempty(ix)~=1
            error('myApp:argChk', ['Power law is not valid above h0 heights,... \nAbort,... \n'])
        end
        hgt_scaled = (h0*1000-hgt).^(alpha);
    end

    % define a new dataset were the height is actually the first column and
    % the interferograms are in the other columns
    % optional remove an estimate for the DEM errors
    if strcmp(DEM_corr,'y') && n_interferograms>5;
       % these are erros that scale with perpendicular baseline
       % loading the perpendicualr baseline information
       bperp = load(bperp_matfile);
        if strcmp(stamps_processed,'y')
          bperp = bperp.bperp;
       end
       % checking the size of bperp
       if size(bperp,2) > 1
          bperp = bperp';
          if size(bperp,2) >1
              error('myApp:argChk', ['bperp is not a vector,... \nAbort,... \n'])
          end
       end
       % estimating the correlated errors
       DEM_corr_e = lscov(bperp,ph')';

       % removing DEM correlated errors
       ph_temp = ph-repmat(bperp',n_points,1).*repmat(DEM_corr_e,1,n_interferograms);
       dataset = [hgt_scaled ph_temp];
       clear ph_temp A
    else
        if strcmp(DEM_corr,'y') && n_interferograms<=5
            fprintf('Not enough interferograms to make a reliable estimate for the DEM error \n')
            DEM_corr = 'n';
        end
        dataset = [hgt_scaled ph];
        DEM_corr_e = zeros([n_points 1]);
    end

    %% do an initial reference correction for each interferogram.
    % This is to scope for regions that are small and which can have bandfiltering
    % arctifacts introduced. What happens is that each interferogram is
    % shifted such that the phase is zero at h0. Note that this is a long
    % wavelength shift that will be filtered out latter on. But in case of
    % irregularities it will reduce introduced arctifacts. To cope with
    % spatial varying signals this offset estimation is done over multiple
    % windows, but the constant offset is the same for all windows.
    % get the position of some subwindows

    [iterate] = window_generation(xy_rotatelocal,1);

    counter=0;
    for k=1:length(iterate.window_ix)
        counter = counter + size(iterate.window_ix{k},1);
    end

    for k=1:n_interferograms
        % separate between ifg based correction or using a single alpha and h0
        if length(alpha)==n_interferograms
            height_data = dataset(:,k);
            scaling = 1./nanmean(abs(dataset(:,k)));
            phase_data = dataset(:,n_interferograms+k);
        else
            scaling = 1./nanmean(abs(dataset(:,1)));
            height_data = dataset(:,1);
            phase_data = dataset(:,1+k);
        end

        % estimating the reference, forced to be the same for all windows
        A_temp = zeros([counter length(iterate.window_ix)]);
        data_temp =  zeros([counter size(phase_data,2)]);
        counter =1;
        for kk=1:length(iterate.window_ix)
            A_temp(counter:counter+size(iterate.window_ix{kk},1)-1,kk) = [height_data(iterate.window_ix{kk},1)]*scaling;
            data_temp(counter:counter+size(iterate.window_ix{kk},1)-1,:) = [phase_data(iterate.window_ix{kk},1)];
            counter =  counter + size(iterate.window_ix{kk},1);

        end
        A_temp = [A_temp ones([size(data_temp,1) 1])];

        % allowing for NaN in interferogram {DB}
%        ix = find(isnan(data_temp(:,1))~=1);
        ix = ~((isnan(data_temp(:,1)) + sum(isnan(A_temp),2))>=1);
        coeff = lscov(A_temp(ix,:),data_temp(ix,1));

        % plot the interferogram before and after the reference shift.
        if plot_flag ==1
           figure('name',['Correction of the reference for ifg ', num2str(k)])
           subplot(2,1,1)
           plot(height_data,phase_data,'k.')
           if length(alpha)==n_interferograms
                temp =getparm_aps('powerlaw_alpha');
                xlabel(['(h_0-h)^{' num2str(temp(k)) '}'])
           else
                xlabel(['(h_0-h)^{' num2str(getparm_aps('powerlaw_alpha')) '}'])
           end
           ylabel('Phase')
        end

       if plot_flag ==1
           for kk=1:length(iterate.window_ix)
                y_temp = [[height_data(iterate.window_ix{kk},1)]*scaling ones(size([height_data(iterate.window_ix{kk},1)]*scaling))]*[coeff(kk) coeff(end)]';
                hold on
                plot([height_data(iterate.window_ix{kk},1)],y_temp,'r-')
           end
           title('power-law plot with ifg reference to be the mean')
       end

       % correct the data for the reference
       data_temp = data_temp - coeff(end);
       % separate between ifg based correction or using a single alpha and h0
       if length(alpha)==n_interferograms
           dataset(:,n_interferograms+k) = dataset(:,n_interferograms+k) - coeff(end);
       else
           dataset(:,k+1) = dataset(:,k+1) - coeff(end);
       end

       % plot the corrected result
       if plot_flag ==1
           subplot(2,1,2)
           if length(alpha)==n_interferograms
               plot(dataset(:,k),dataset(:,n_interferograms+k),'k.')
               temp =getparm_aps('powerlaw_alpha');
               xlabel(['(h_0-h)^{' num2str(temp(k)) '}'])
           else
              plot(dataset(:,1),dataset(:,k+1),'k.')
              xlabel(['(h_0-h)^{' num2str(getparm_aps('powerlaw_alpha')) '}'])
           end
           ylabel('Phase shifted to h0 reference')

           % plot line on top
           coeff = lscov(A_temp,data_temp(:,1));
           for kk=1:length(iterate.window_ix)
                if length(alpha)==n_interferograms
                    y_temp = [[dataset(iterate.window_ix{kk},k)]*scaling ones(size([dataset(iterate.window_ix{kk},k)]*scaling))]*[coeff(kk) coeff(end)]';
                else
                    y_temp = [[dataset(iterate.window_ix{kk},1)]*scaling ones(size([dataset(iterate.window_ix{kk},1)]*scaling))]*[coeff(kk) coeff(end)]';
                end
                hold on
                plot([dataset(iterate.window_ix{kk},k)],y_temp,'r-')

           end
           title('power law plot with h0 as reference')
       end
    end

    %% Extrapolate such the convexhull contains the full grid
    xy_res = getparm_aps('powerlaw_xy_res');
    [xy_local_box,dataset_box] = extrapolate_local_new(xy_rotatelocal,dataset, xy_res(1), xy_res(2));

    % Interpolate to a regular grid
    [X_regular,Y_regular,Z_regular] = interpolate_regular(xy_local_box,dataset_box, xy_res(1), xy_res(2));

    % Saving the results of this step
    aps_save([save_path filesep 'aps_p_step2.mat'],xy_rotatelocal,xy_rotatelocal_full,ifg_based_correction,xy_local,xy_local_full,heading,ll,local_origin,local_origin_full,X_regular,Y_regular,Z_regular,xy_rotatelocal,xy_res,DEM_corr,DEM_corr_e,ix_remove,ix_phase_nan)

end

%% PART 3: Bandfiltering and interpolation back to a local grid
if start_step<=3 && end_step >=3
    fprintf('Step 3: Bandfiltering of the regular grid and converting back to a local grid \n\n')
    % Loading the data from the previous step
    if start_step==3
        load([save_path filesep 'aps_p_step2.mat'])
    end

     % loading the lonlat information
    ll_matfile = getparm_aps('ll_matfile');
    temp = load(ll_matfile);
    if strcmp(stamps_processed,'y')
        temp = temp.lonlat;
    end

    n_points_original = size(temp,1);
    clear temp;

    % getting the extend of the spatial bandfilters
    spatial_bands = getparm_aps('powerlaw_spatial_bands');

    % Bandfiltering
    bandfiltering(Z_regular,xy_res(1),xy_res(2),spatial_bands,save_path,ifg_based_correction)

    % Interpolate back to local grid
    n_datasets = size(Z_regular,3);
    % calling the function for each bandfiltered dataset
    fprintf(['*Interpolate the band filtered data to the local grid: \n'])
    if strcmp(ifg_based_correction,'y')
        h_ifg_number = n_datasets/2;
    else
        h_ifg_number = 1;
    end
    for k=1:n_datasets

        % the file names, vary depending if the correction is varying for each interferogram
        if k<=h_ifg_number && h_ifg_number~=1
            % thse are interferograms
            load_name = ['bandfilter_regular_hgt_ifg_' num2str(k) '.mat'];
            save_name = ['bandfilter_local_hgt_ifg' num2str(k) '.mat'];
        elseif k<=h_ifg_number && h_ifg_number==1
            % these are the heights
            load_name = 'bandfilter_regular_hgt.mat';
            save_name = 'bandfilter_local_hgt.mat';
        else
            % thse are interferograms
             load_name = ['bandfilter_regular_ifg_' num2str(k-h_ifg_number) '.mat'];
             save_name = ['bandfilter_local_ifg_' num2str(k-h_ifg_number) '.mat'];
        end

        dataset_band = load([save_path filesep load_name]);

        % Storing the information what type of filtering was done for which band
        dimension_filter = dataset_band.dimension_filter;      % 1 for 1D, 2 for 2D
        % loading the data for each dataset
        dataset_band = dataset_band.data_band_out;

        % computing the bandfiltered data at the local grid
        % note that this is xy_rotatelocal and does not include the cropped
        % region!
        [z_local_out_temp] = interpolate2local(X_regular,Y_regular,dataset_band,xy_rotatelocal);
        % in case data has been cropped then path the cropped region with NaNs
        z_local_out = NaN([n_points_original size(z_local_out_temp,2)]);

        if isempty(ix_remove)
            z_local_out = z_local_out_temp;
        else
            z_local_out(~ix_remove,:) = z_local_out_temp;
        end

        % Check in case some of the orginal phase values were nan.
        % if so set them to nan here too.
        if k<=h_ifg_number && h_ifg_number~=1
            % these are interferograms heights not the phase. No need to
            % modify this here.
        elseif k<=h_ifg_number && h_ifg_number==1
            % these are the heights, nothing needs to be done as the
            % heights remain the same for each interferogram.
        else
            % these are interferograms
            fprintf('Allowing for NaN''s in the interferograms \n')
            z_local_out(ix_phase_nan(:,k-h_ifg_number),:)=NaN;
        end

        % saving the data for each dataset
        aps_save([save_path filesep save_name],z_local_out,spatial_bands,dimension_filter);
        % outputting the progress to the user
        fprintf(['Progress: ' num2str(k) '/' num2str(n_datasets) ' done \n'])
    end

    % Saving the results of this step
    aps_save([save_path filesep 'aps_p_step3.mat'],ifg_based_correction,n_datasets,local_origin,local_origin_full,xy_rotatelocal,xy_rotatelocal_full,DEM_corr,DEM_corr_e,dimension_filter,ix_remove,ix_phase_nan)
end

%% PART 4: Estimating the scaling coefficient of the powerlaw for the local dataset
% and computing the corresponding troposheric delay
if start_step<=4 && end_step >=4

    fprintf('\n\nStep 4: Estimating scaling coefficient of the powerlaw locally \n')
    fprintf('and estimating the tropospheric interferometric phase delay \n\n')

    % Loading the data from the previous step
    if start_step==4
        load([save_path filesep 'aps_p_step3.mat'])
    end

    % loading information from the parm file
    hgt_matfile = getparm_aps('hgt_matfile');
    powerlaw_all_bands = getparm_aps('powerlaw_all_bands');
    hgt = load(hgt_matfile);
    hgt = hgt.hgt;

    % checking if its an interferogram based correction or not:
    if strcmp(ifg_based_correction,'y')
        n_interferograms = n_datasets/2;
    else
        n_interferograms = n_datasets-1;
    end

    % initialization
    n_points = size(hgt,1);
    ph_tropo_powerlaw = NaN([n_points n_interferograms]);
    K_tropo_powerlaw = NaN([n_points n_interferograms]);
    perc_cons_K_powerlaw= NaN([ n_interferograms 1]);
    if strcmp(powerlaw_all_bands,'y')
        ph_tropo_powerlaw_bands = NaN([n_points n_interferograms size(getparm_aps('powerlaw_spatial_bands'),1)]);
        K_tropo_powerlaw_bands = NaN([n_points n_interferograms size(getparm_aps('powerlaw_spatial_bands'),1)]);
    end

    % Linear estimation for multiple windows based on bandfitlered data
    % final slope value is estimated for reliable bands only
    iterate = [];
    % getting the patch information from the parm_aps file

    % processing only those images that are not dropped in case of stamps
    ix_ifgs = 1:n_interferograms;
    if strcmp(stamps_processed,'y')
       ix_drop_ifgs = getparm('drop_ifg');
    else
        ix_drop_ifgs = [];
    end
    ix_ifgs(ix_drop_ifgs)=[];

    for k=1:length(ix_ifgs)

        ix_interferogram = ix_ifgs(k);

        % the file names, vary depending if the correction is varying for each interferogram
        if strcmp(ifg_based_correction,'y') && n_interferograms~=1
            % the first dataset was set the be the topography
            topo_data = load([save_path filesep 'bandfilter_local_hgt_ifg' num2str(ix_interferogram) '.mat']);
            topo_data = topo_data.z_local_out;
            % the other dataset were the phase
            phase_data = load([save_path filesep 'bandfilter_local_ifg_' num2str(ix_interferogram) '.mat']);
            phase_data = phase_data.z_local_out;
            ifg_number=k;
        else
            % the first dataset was set the be the topography
            topo_data = load([save_path filesep 'bandfilter_local_hgt.mat']);
            topo_data = topo_data.z_local_out;
            % the other dataset were the phase
            phase_data = load([save_path filesep 'bandfilter_local_ifg_' num2str(ix_interferogram) '.mat']);
            phase_data = phase_data.z_local_out;
            ifg_number=1;
        end

        % keep the figures of the outlier rejection as validation
        if k==0
            save_path_outlier_windows = [save_path filesep 'ifg_' num2str(ix_interferogram) '_outlier_window'];
        else
            save_path_outlier_windows=[];
        end

        if k==10000
            debug_fig_linear=1;
        else
            debug_fig_linear=0;
        end

        [ph_tropo_powerlaw_temp,ph_tropo_powerlaw_band_temp,iterate,K_temp,K_band,perc_cons_K] = aps_powerlaw_linear_local(xy_rotatelocal_full,local_origin_full,topo_data,phase_data,iterate,dimension_filter,save_path_outlier_windows,ifg_number,ix_phase_nan(:,ix_interferogram));

        % storing of the data
        % Computation of the tropospheric delay for all ifgs
        ph_tropo_powerlaw(:,ix_interferogram) = ph_tropo_powerlaw_temp;
        K_tropo_powerlaw(:,ix_interferogram) = K_temp;
        perc_cons_K_powerlaw(ix_interferogram,1) = perc_cons_K;
        if strcmp(powerlaw_all_bands,'y')
            % keeping all bands with in the columns the ifgs and the third
            % dimention the bands
            ph_tropo_powerlaw_bands(:,ix_interferogram,:) = ph_tropo_powerlaw_band_temp;
            K_tropo_powerlaw_bands(:,ix_interferogram,:)=K_band;

            % keeping the interferogram and the respective bands.
            eval(['ph_tropo_powerlaw_ifg_' num2str(ix_interferogram) '= ph_tropo_powerlaw_band_temp;']);
            eval(['K_tropo_powerlaw_ifg_' num2str(ix_interferogram) '= K_band;']);

            clear ph_tropo_powerlaw_band_temp
            % saving the data from this step
            if strcmp(stamps_processed,'y')
                % This is StaMPS
                if strcmp(getparm('small_baseline_flag'),'y')
                    if exist(apsbandssbname,'file')==2
                        save(apsbandssbname,'-append',['ph_tropo_powerlaw_ifg_' num2str(ix_interferogram)],['K_tropo_powerlaw_ifg_' num2str(ix_interferogram)])
                    else
                        save(apsbandssbname,['ph_tropo_powerlaw_ifg_' num2str(ix_interferogram)],['K_tropo_powerlaw_ifg_' num2str(ix_interferogram)])
                    end
                else
                    if exist(apsbandsname,'file')==2
                        save(apsbandsname,'-append',['ph_tropo_powerlaw_ifg_' num2str(ix_interferogram)],['K_tropo_powerlaw_ifg_' num2str(ix_interferogram)])
                    else
                        save(apsbandsname,['ph_tropo_powerlaw_ifg_' num2str(ix_interferogram)],['K_tropo_powerlaw_ifg_' num2str(ix_interferogram)])
                    end
                end
            else
                % This is not StaMPS
                if exist(apsbandsname,'file')==2
                    save(apsbandsname,'-append',['ph_tropo_powerlaw_ifg_' num2str(ix_interferogram)],['K_tropo_powerlaw_ifg_' num2str(ix_interferogram)])
                else
                    save(apsbandsname,['ph_tropo_powerlaw_ifg_' num2str(ix_interferogram)],['K_tropo_powerlaw_ifg_' num2str(ix_interferogram)])
                end
            end
            eval(['clear ph_tropo_powerlaw_ifg_' num2str(ix_interferogram) ' K_tropo_powerlaw_ifg_' num2str(ix_interferogram)  ';']);
        end

        % give output to the screen
        fprintf(['Progress ifgs: ' num2str(k) '/' num2str(length(ix_ifgs)) ' done \n'])
    end
    % saving the data from this step
    aps_save([save_path filesep 'aps_p_step4.mat'],ph_tropo_powerlaw,DEM_corr,DEM_corr_e,iterate,K_tropo_powerlaw)

    % saving the data into the final variable togehter were other
    % correction values are saved.
    ph_tropo_powerlaw_or= ph_tropo_powerlaw;
    K_tropo_powerlaw_or=K_tropo_powerlaw;
    if strcmp(stamps_processed,'y')
        % This is StaMPS
        if strcmp(getparm('small_baseline_flag'),'y')
            if exist(apssbname,'file')==2
                save(apssbname,'-append','ph_tropo_powerlaw','K_tropo_powerlaw','perc_cons_K_powerlaw')
            else
                save(apssbname,'ph_tropo_powerlaw','K_tropo_powerlaw','perc_cons_K_powerlaw')
            end

            % saving the band information
            if strcmpi(powerlaw_all_bands,'y')
                if exist(apsbandssbname,'file')==2
                    save(apsbandssbname,'-append','ph_tropo_powerlaw_or','K_tropo_powerlaw_or','ph_tropo_powerlaw_bands','K_tropo_powerlaw_bands','perc_cons_K_powerlaw')
                else
                    save(apsbandssbname,'ph_tropo_powerlaw_or','K_tropo_powerlaw_or','ph_tropo_powerlaw_bands','K_tropo_powerlaw_bands','perc_cons_K_powerlaw')
                end
            end
        else
            % this is single master
            if exist(apsname,'file')==2
                save(apsname,'-append','ph_tropo_powerlaw','K_tropo_powerlaw','perc_cons_K_powerlaw')
            else
                save(apsname,'ph_tropo_powerlaw','K_tropo_powerlaw','perc_cons_K_powerlaw')
            end

            % saving the band information
            if strcmpi(powerlaw_all_bands,'y')
                if exist(apsbandsname,'file')==2
                    save(apsbandsname,'-append','ph_tropo_powerlaw_or','K_tropo_powerlaw_or','ph_tropo_powerlaw_bands','K_tropo_powerlaw_bands','perc_cons_K_powerlaw')
                else
                    save(apsbandsname,'ph_tropo_powerlaw_or','K_tropo_powerlaw_or','ph_tropo_powerlaw_bands','K_tropo_powerlaw_bands','perc_cons_K_powerlaw')
                end
            end
        end
    else
        % This is not StaMPS
        if exist(apsname,'file')==2
            save(apsname,'-append','ph_tropo_powerlaw','K_tropo_powerlaw','perc_cons_K_powerlaw')
        else
            save(apsname,'ph_tropo_powerlaw','K_tropo_powerlaw','perc_cons_K_powerlaw')
        end

         % saving the band information
        if strcmpi(powerlaw_all_bands,'y')
            if exist(apsbandsname,'file')==2
                save(apsbandsname,'-append','ph_tropo_powerlaw_or','K_tropo_powerlaw_or','ph_tropo_powerlaw_bands','K_tropo_powerlaw_bands','perc_cons_K_powerlaw')
            else
                save(apsbandsname,'ph_tropo_powerlaw_or','K_tropo_powerlaw_or','ph_tropo_powerlaw_bands','K_tropo_powerlaw_bands','perc_cons_K_powerlaw')
            end
        end
    end
    setparm_aps('powerlaw_kept',[0]);

end

if start_step<=5 && end_step >=5 &&  exist('aps_RMSE_comparison')==2
    fprintf('\n\nStep 5: Powerlaw band comparison with reference\n')

    if strcmpi(getparm_aps('stamps_processed'),'y')
        % check which technique it needs to be compared with
        technique_number = -10;
        while isnumeric(technique_number)~=1 | (technique_number<1 | technique_number>7)
            technique_number = input(['*Select the technique number you want to compare with: \n---- MERIS (1), MODIS (2), MODIS Recal (3), ERA-I (4), WRF (5), u(sb) (6), u(sb)-d (7) : \n'],'s');
            technique_number = floor(str2num(technique_number));
        end
        % the technique component as a string
        if technique_number==1
            comp_str = 'a_mi';
        elseif technique_number==2
            comp_str = 'a_MI';
        elseif technique_number==3
            comp_str = 'a_RMI';
        elseif technique_number==4
            comp_str = 'a_e';
        elseif technique_number==5
            comp_str = 'a_w';
        elseif technique_number==6
            if strcmpi(getparm('small_baseline_flag'),'y')
                comp_str = 'usb';
            else
                comp_str = 'u';
            end
        elseif technique_number==7
            if strcmpi(getparm('small_baseline_flag'),'y')
                comp_str = 'usb-d';
            else
                comp_str = 'u-d';
            end
        end

        % check if we should add the hydrostatic delays too in case this is
        % possible
        if technique_number==1 | technique_number==2 | technique_number==3
            hydro_comp = inf;
            while isnumeric(hydro_comp)~=1 | (hydro_comp<0 | hydro_comp>2)
                hydro_comp = input(['*Do you want to add a hydrostatic component:\n---- None (0), ERA-I hydr (1), WRF hydro (2): \n'],'s');
                hydro_comp = floor(str2num(hydro_comp));
            end
        else
            hydro_comp = 0;
        end
        % the hydrostatic component as a string
        hydro_str = '';
        if hydro_comp==1
            hydro_str = '+a_eh';
        elseif hydro_comp==2
            hydro_str = '+a_wh';
        end

        % The reference technique is combination of both
        tech_ref = [comp_str hydro_str];
        clear hydro_str comp_str technique_number hydro_comp

        % running the RMSE script
        aps_RMSE_comparison(tech_ref,'a_pbands')
        str = '';
        while strcmpi(str,'y')~=1 && strcmpi(str,'n')~=1
            str = input(['Do you want to update ''powerlaw_kept'' to a different band? [y/n] \n'],'s');
        end
        if strcmpi(str,'y')
            bandnumber = -10;
            while (bandnumber<0 | bandnumber>size(getparm_aps('powerlaw_spatial_bands'),1))
                bandnumber = input(['Which the band? (row number of the band to keep, or 0 for all bands combined) \n'],'s');
                bandnumber = str2num(bandnumber);
                if isempty(bandnumber)
                    bandnumber = -10;
                end
            end
            setparm_aps('powerlaw_kept',[bandnumber]);
        end
    else
       fprintf('This option uses the StaMPS ps_plot function, and therefore only works if your data is stamps processed \n')
    end

end

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
echo on
% This script shows how to use the ga using a float representation.
% You should see the demos for
% more information as well. gademo1, gademo2, gademo3
global bounds

% Setting the seed back to the beginning for comparison sake
rand('seed',0)

% Crossover Operators
xFns = 'simpleXover';
xOpts = [.4];

% Mutation Operators
mFns = 'binaryMutation';

mOpts = [0.005];

% Termination Operators
termFns = 'maxGenTerm';
termOps = [200]; % 200 Generations

% Selection Function
selectFn = 'roulette'
selectOps = [];

% Evaluation Function
evalFn = 'gaMichEval';
evalOps = [];

type gaMichEval

% Bounds on the variables
bounds = [-3 12.1; 4.1 5.8];

% GA Options [epsilon float/binar display]
gaOpts=[1e-6 0 1];

% Generate an intialize population of size 20
startPop = initializega(20,bounds,'gaMichEval',[],[1e-6 0]);

% Lets run the GA
% Hit a return to continue
pause

[x endPop bestPop trace]=ga(bounds,evalFn,evalOps,startPop,gaOpts,...
    termFns,termOps,selectFn,selectOps,xFns,xOpts,mFns,mOpts);

% x is the best solution found
x
% Hit a return to continue
pause

% endPop is the ending population
endPop
% Hit a return to continue
pause

% trace is a trace of the best value and average value of generations
trace
% Hit a return to continue
pause

% Plot the best over time
clf
plot(trace(:,1),trace(:,2));
% Hit a return to continue
pause

% Add the average to the graph
hold on
plot(trace(:,1),trace(:,3));
% Hit a return to continue
pause

% Lets increase the population size by running the defaults
%
rand('seed',0)
termOps=[100];
[x endPop bestPop trace]=ga(bounds,evalFn,evalOps,[],gaOpts,termFns,termOps,...
    selectFn,selectOps);

% x is the best solution found
x
% Hit a return to continue
pause

% endPop is the ending population
endPop
% Hit a return to continue
pause

% trace is a trace of the best value and average value of generations
trace
% Hit a return to continue
pause

% Plot the best over time
clf
plot(trace(:,1),trace(:,2));
% Hit a return to continue
pause

% Add the average to the graph
hold on
plot(trace(:,1),trace(:,3));

echo off

author: HuangCongQing


text:
function r = subsref(a,s)
%SUBSREF      Implements subscripted references for Hessians
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/03/12     S.M. Rump  INTLAB_HESSIAN_DERIV_ERROR removed
% modified 10/06/12     S.M. Rump  internal use
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  N = getappdata(0,'INTLAB_HESSIAN_NUMVAR');

  while 1
    if ~isa(a,'hessian')                 % index reference a.x(i) etc.
      r = subsref(a,s(1));
    elseif strcmp(s(1).type,'()')        % index reference a(i)
      r.x = a.x(s(1).subs{:});
      index = reshape(1:prod(size(a.x)),size(a.x));
      index = index(s(1).subs{:});
      r.dx = a.dx(:,index(:));
      r.hx = a.hx(:,index(:));
      % avoid Matlab 6.5f bug:
      % a = sparse([],[],[],1,1); a = reshape(a,1,1); abs(a)
      % produces  9.6721e-317  or similar number in underflow range
      if prod(size(r.x))==1
        r.x = full(r.x);
      end
      if prod(size(r.dx))==1
        r.dx = full(r.dx);
      end
      if prod(size(r.hx))==1
        r.hx = full(r.hx);
      end
      r = class(r,'hessian');
    elseif strcmp(s(1).type,'.')         % index reference a.x, a.dx or a.hx
      if strcmp(s(1).subs,'x')
        r = a.x;
      elseif strcmp(s(1).subs,'dx')
        sizeax = size(a.x);
        if ( length(sizeax)==2 ) & ( sizeax(2)==1 )
          sizeax = sizeax(1);            % row gradient for column vector a.x
        end
        if issparse(a.dx) & ( length(sizeax)>1 )
          error('access of .dx sparse hessian with more than one column, see hessianinit')
        end
        r = reshape(transpose(a.dx),[sizeax N]);
      elseif strcmp(s(1).subs,'ddx')
        r = a.dx;
      elseif strcmp(s(1).subs,'hx')
        if issparse(a.hx) & ( prod(size(a.x))>1 )
          error('access of .hx of non-scalar sparse hessian, see hessianinit')
        end
        index = reshape(1:N*N,N,N)';
        r = a.hx + a.hx(index(:),:);     % Hessian is .hx + transpose(.hx)
        sizeax = size(a.x);
        if prod(sizeax)==1
          r = reshape(r,N,N);
        else
          r = reshape(r,[sizeax N N]);
        end
      elseif strcmp(s(1).subs,'hhx')
        index = reshape(1:N*N,N,N)';
        r = a.hx + a.hx(index(:),:);     % Hessian is .hx + transpose(.hx)
      elseif strcmp(s(1).subs,'mid')
        r = mid(a);
      else
        error('invalid subscript reference for hessian')
      end
    else
      error('invalid index reference for hessian')
    end
    if length(s)==1
      if rndold
        setround(rndold)
      end
      return
    end
    s = s(2:end);
    a = r;
  end

repo: OpenMAS


text:
function [p,S] = polyfit(x,y,n)
    %POLYFIT Polynomial curve fitting.
    %	POLYFIT(x,y,n) finds the coefficients of a polynomial p(x) of
    %	degree n that fits the data, p(x(i)) ~= y(i), in a least-squares sense.
    %
    %	[p,S] = POLYFIT(x,y,n) returns the polynomial coefficients p and a
    %	matrix S for use with POLYVAL to produce error estimates on predictions.
    %	If the errors in the data, y, are independent normal with constant
    %	variance, POLYVAL will produce error bounds which contain at least 50%
    %	of the predictions.
    %
    %	See also POLY, POLYVAL, ROOTS.

    %	J.N. Little 4-21-85, 8-23-86; CBM, 12-27-91 BAJ, 5-7-93.
    %	Copyright (c) 1984-94 by The MathWorks, Inc.

    % The regression problem is formulated in matrix format as:
    %
    %    y = V*p    or
    %
    %          3  2
    %    y = [x  x  x  1] [p3
    %                      p2
    %                      p1
    %                      p0]
    %
    % where the vector p contains the coefficients to be found.  For a
    % 7th order polynomial, matrix V would be:
    %
    % V = [x.^7 x.^6 x.^5 x.^4 x.^3 x.^2 x ones(size(x))];

    %report_this_filefun(mfilename('fullpath'));

    if any(size(x) ~= size(y))
        error('X and Y vectors must be the same size.')
    end

    x = x(:);
    y = y(:);

    % Construct Vandermonde matrix.
    V(:,n+1) = ones(length(x),1);
    for j = n:-1:1
        V(:,j) = x.*V(:,j+1);
    end

    % Solve least squares problem.
    [Q,R] = eval('qr(V,0)','qr(V)');

    % The current PC version does not have the two-argument form of qr
    [rows, cols] = size(R);
    if rows ~= cols
        R = R(1:cols,:);
        Q = Q(:,1:cols);
    end

    p = R\(Q'*y);    % Same as p = V\y;
    r = y - V*p;
    p = p';          % Polynomial coefficients are row vectors by convention.

    % S is a structure containing three elements: the Cholesky factor of the
    % Vandermonde matrix, the degrees of freedom and the norm of the residuals.

    df = length(y) - (n+1);
    S = [R; [df zeros(1,n)]; [norm(r) zeros(1,n)]];

author: CelsoReyes


text:
function [ t_start, y_start ] = p29_start ( neqn )

%*****************************************************************************80
%
%% P29_START returns the starting point for problem p29.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_START, Y_START(NEQN,1), the initial data.
%
  y_start = zeros ( neqn, 1 );

  t_start = 0.0;
  y_start = 1.0;

  return
end


text:
function value = r4_cos ( x )

%*****************************************************************************80
%
%% R4_COS evaluates the cosine of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the cosine of X.
%

  persistent ntsn
  persistent pi2rec
  persistent pihi
  persistent pilo
  persistent pirec
  persistent sincs
  persistent xmax
  persistent xsml
  persistent xwarn
%
%  pihi + pilo = pi.  pihi is exactly representable on all machines
%  with at least 8 bits of precision.  whether it is exactly
%  represented depends on the compiler.  this routine is more
%  accurate if it is exactly represented.
%
  pi2 = 1.57079632679489661923;
  pi2rec = 0.636619772367581343;
  pihi = 3.140625;
  pilo = 9.6765358979323846E-04;
  pirec = 0.31830988618379067;

  if ( isempty ( ntsn ) )

    sincs = [ ...
      -0.374991154955873175840, ...
      -0.181603155237250201864, ...
      +0.005804709274598633559, ...
      -0.000086954311779340757, ...
      +0.000000754370148088851, ...
      -0.000000004267129665056, ...
      +0.000000000016980422945, ...
      -0.000000000000050120579, ...
      +0.000000000000000114101, ...
      -0.000000000000000000206 ]';

    ntsn = r4_inits ( sincs, 10, 0.1 * r4_mach ( 3 ) );
    xsml = sqrt ( 2.0 * r4_mach ( 3 ) );
    xmax = 1.0 / r4_mach ( 4 );
    xwarn = sqrt ( xmax );

  end

  absx = abs ( x );
  y = absx + pi2;

  if ( xmax < y )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_COS - Warning!\n' );
    fprintf ( 1, '  No precision because |X| is big.\n' );
    value = 0.0;
    return
  end

  if ( xwarn < y )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_COS - Warning!\n' );
    v '  Answer < half precision because |X| is big.\n' );
  end

  value = 1.0;

  if ( absx < xsml )
    return
  end

  xn = r4_aint ( y * pirec + 0.5 );
  n2 = r4_aint ( mod ( xn, 2.0 ) + 0.5 );
  xn = xn - 0.5;
  f = ( absx - xn * pihi ) - xn * pilo;

  xn = 2.0 * ( f * pi2rec ) * ( f * pi2rec ) - 1.0;
  value = f + f * r4_csevl ( xn, sincs, ntsn );

  if ( n2 ~= 0 )
    value = - value;
  end

  if ( value < - 1.0 )
    value = - 1.0;
  elseif ( 1.0 < value )
    value = + 1.0;
  end

  return
end


text:
function tutorial_fem(tutorial_dir)
% TUTORIAL_FEM: Script that reproduces the online tutorial "FEM tutorial: MEG/EEG Median nerve stimulation"
%
% REFERENCE:
%     https://neuroimage.usc.edu/brainstorm/Tutorials/FemMedianNerve
%
% INPUTS:
%     tutorial_dir: Directory where the sample_fem.zip file has been unzipped

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Francois Tadel, 2021

% ===== FILES TO IMPORT =====
% You have to specify the folder in which the tutorial dataset is unzipped
if (nargin == 0) || isempty(tutorial_dir) || ~file_exist(tutorial_dir)
    error('The first argument must be the full path to the tutorial dataset folder.');
end
% Build the path of the files to import
T1Nii  = fullfile(tutorial_dir, 'sample_fem', 'sub-fem01', 'ses-mri', 'anat', 'sub-fem01_ses-mri_T1w.nii.gz');
T2Nii  = fullfile(tutorial_dir, 'sample_fem', 'sub-fem01', 'ses-mri', 'anat', 'sub-fem01_ses-mri_T2w.nii.gz');
DwiNii = fullfile(tutorial_dir, 'sample_fem', 'sub-fem01', 'ses-mri', 'dwi', 'sub-fem01_ses-mri_dwi.nii.gz');
DwiBval = fullfile(tutorial_dir, 'sample_fem', 'sub-fem01', 'ses-mri', 'dwi', 'sub-fem01_ses-mri_dwi.bval');
DwiBvec = fullfile(tutorial_dir, 'sample_fem', 'sub-fem01', 'ses-mri', 'dwi', 'sub-fem01_ses-mri_dwi.bvec');
FifFile = fullfile(tutorial_dir, 'sample_fem', 'sub-fem01', 'ses-meg', 'meg', 'sub-fem01_ses-meg_task-mediannerve_run-01_proc-tsss_meg.fif');
% Check if the folder contains the required files
if ~file_exist(T1Nii) || ~file_exist(FifFile)
    error(['The folder ' tutorial_dir ' does not contain the folder from the file sample_fem.zip.']);
end
% Subject name
SubjectName = 'Subject01';
% Latency of interest
Latency = 0.022;

% ===== CHECK SOFTWARE DEPENDENCIES =====
% Start brainstorm without the GUI
if ~brainstorm('status')
    brainstorm nogui
end
% SimNIBS
status = system('headreco --version');
if (status ~= 0)
    error('SimNIBS is not installed or not added to the system path: the command "headreco" could not be found.');
end
% BrainSuite
if ~file_exist(bst_fullfile(bst_get('BrainSuiteDir'), 'bin'))
    error('BrainSuite is not configured in the Brainstorm preferences.');
end
% SPM12
[isInstalled, errMsg] = bst_plugin('Install', 'spm12', 0);
if ~isInstalled
    error(['Could not install plugin: spm12' 10 errMsg]);
end
% CAT12
[isInstalled, errMsg] = bst_plugin('Install', 'cat12', 0);
if ~isInstalled
    error(['Could not install plugin: cat12' 10 errMsg]);
end
% Iso2mesh
[isInstalled, errMsg] = bst_plugin('Install', 'iso2mesh', 0);
if ~isInstalled
    error(['Could not install plugin: iso2mesh' 10 errMsg]);
end

% ===== CREATE PROTOCOL =====
% The protocol name has to be a valid folder name (no spaces, no weird characters...)
ProtocolName = 'TutorialFem';
% Delete existing protocol
gui_brainstorm('DeleteProtocol', ProtocolName);
% Create new protocol
gui_brainstorm('CreateProtocol', ProtocolName, 0, 1);
% Start a new report
bst_report('Start');

%% ===== IMPORT ANATOMY =====
% ===== IMPORT MRI VOLUMES =====
% Create subject
[sSubject, iSubject] = db_add_subject(SubjectName, [], 0, 0);
% Import T1 MRI
T1File = import_mri(iSubject, T1Nii, 'ALL', 0, 0);
% Compute the MNI normalization
bst_normalize_mni(T1File);
% Import T2 MRI
T2File = import_mri(iSubject, T2Nii, 'ALL', 0, 0);
% Volumes are not registered: Register with SPM
mri_coregister(T2File, T1File, 'spm', 1);
% Delete the non-registered T2
file_delete(file_fullpath(T2File), 1);
db_reload_subjects(iSubject);

% ===== IMPORT DTI =====
% Process: Convert DWI to DTI (BrainSuite)
bst_process('CallProcess', 'process_dwi2dti', [], [], ...
    'subjectname', SubjectName, ...
    'dwifile',     {DwiNii, 'DWI-NII'}, ...
    'bvalfile',    {DwiBval, 'DWI-BVAL'}, ...
    'bvecfile',    {DwiBvec, 'DWI-BVEC'});

% ===== RUN SIMNIBS =====
% Process: Generate FEM mesh
bst_process('CallProcess', 'process_fem_mesh', [], [], ...
    'subjectname',   SubjectName, ...
    'method',        'simnibs', ...  % SimNIBS:Call SimNIBS to segment and mesh the T1 (and T2) MRI.
    'vertexdensity', 0.5, ...
    'nvertices',     15000, ...
    'zneck',         0);
% Select default cortex: Low-resolution
[sSubject, iSubject, iCortex] = bst_get('SurfaceFile', bst_fullfile(SubjectName, 'tess_cortex_pial_low.mat'));
db_surface_default(iSubject, 'Cortex', iCortex);

% ===== REMESH WITH ISO2MESH =====
% Part skipped because:
% 1) Not needed: DUNEuro works with cavities in this case
% 2) No process call for: Extract surfaces (import_femlayers)
% 3) Remeshing with iso2mesh doesn't work...

% ===== COMPUTE FEM TENSORS =====
% Process: Compute FEM tensors
bst_process('CallProcess', 'process_fem_tensors', [], [], ...
    'subjectname', SubjectName, ...
    'femcond',     struct(...
         'FemCond',         [0.14, 0.33, 1.79, 0.008, 0.43], ...
         'isIsotropic',     [0, 1, 1, 1, 1], ...
         'AnisoMethod',     'ema+vc', ...
         'SimRatio',        10, ...
         'SimConstrMethod', 'wolters'));

% ===== COMPUTE BEM SURFACES =====
% Process: Generate BEM surfaces
bst_process('CallProcess', 'process_generate_bem', [], [], ...
    'subjectname', SubjectName, ...
    'nscalp',      1922, ...
    'nouter',      1922, ...
    'ninner',      1922, ...
    'thickness',   4);

%% ===== ACCESS THE RECORDINGS =====
% Process: Create link to raw file
sFilesRaw = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
    'subjectname',    SubjectName, ...
    'datafile',       {FifFile, 'FIF'}, ...
    'channelreplace', 0, ...
    'channelalign',   1, ...    % Automatic registration with digitized head points
    'evtmode',        'value');

% Process: Events: Read from channel
sFilesRaw = bst_process('CallProcess', 'process_evt_read', sFilesRaw, [], ...
    'stimchan',  'STI101', ...
    'trackmode', 1, ...  % Value: detect the changes of channel value
    'zero',      0);

% Process: Project electrodes on scalp
bst_process('CallProcess', 'process_channel_project', sFilesRaw, []);

% Process: Snapshot: Sensors/MRI registration
bst_process('CallProcess', 'process_snapshot', sFilesRaw, [], ...
    'target',   1, ...  % Sensors/MRI registration
    'modality', 1, ...  % MEG (All)
    'orient',   1, ...  % left
    'Comment',  'MEG/MRI Registration');
% Process: Snapshot: Sensors/MRI registration
bst_process('CallProcess', 'process_snapshot', sFilesRaw, [], ...
    'target',   1, ...  % Sensors/MRI registration
    'modality', 4, ...  % EEG
    'orient',   1, ...  % left
    'Comment',  'EEG/MRI Registration');

%% ===== PRE-PROCESSING =====
% Process: Power spectrum density (Welch)
sFilesPsd = bst_process('CallProcess', 'process_psd', sFilesRaw, [], ...
    'timewindow',  [18, 148.999], ...
    'win_length',  5, ...
    'win_overlap', 50, ...
    'units',       'physical', ...  % Physical: U2/Hz
    'sensortypes', 'MEG, EEG', ...
    'win_std',     0, ...
    'edit',        struct(...
         'Comment',         'Power', ...
         'TimeBands',       [], ...
         'Freqs',           [], ...
         'ClusterFuncTime', 'none', ...
         'Measure',         'power', ...
         'Output',          'all', ...
         'SaveKernel',      0));

% Process: Snapshot: Frequency spectrum
bst_process('CallProcess', 'process_snapshot', sFilesPsd, [], ...
    'target',   10, ...  % Frequency spectrum
    'Comment',  'Power spectrum density');

% Process: Band-pass:20Hz-250Hz
sFilesBand = bst_process('CallProcess', 'process_bandpass', sFilesRaw, [], ...
    'sensortypes', 'MEG, EEG', ...
    'highpass',    20, ...
    'lowpass',     250, ...
    'tranband',    0, ...
    'attenuation', 'strict', ...  % 60dB
    'ver',         '2019', ...  % 2019
    'mirror',      0, ...
    'read_all',    0);

% Process: Notch filter: 60Hz 120Hz 180Hz
sFilesClean = bst_process('CallProcess', 'process_notch', sFilesBand, [], ...
    'sensortypes', 'MEG, EEG', ...
    'freqlist',    [60, 120, 180], ...
    'cutoffW',     2, ...
    'useold',      0, ...
    'read_all',    0);

% Process: Delete selected files
bst_process('CallProcess', 'process_delete', [sFilesRaw, sFilesBand], [], ...
    'target', 2);  % Delete selected folders

% Process: Re-reference EEG
bst_process('CallProcess', 'process_eegref', sFilesClean, [], ...
    'eegref',      'AVERAGE', ...
    'sensortypes', 'EEG');

%% ===== IMPORT RECORDINGS =====
% Process: Import MEG/EEG: Events
sFilesEpochs = bst_process('CallProcess', 'process_import_data_event', sFilesClean, [], ...
    'subjectname', SubjectName, ...
    'condition',   '', ...
    'eventname',   '2', ...
    'timewindow',  [18, 148.999], ...
    'epochtime',   [-0.1, 0.2], ...
    'createcond',  0, ...
    'ignoreshort', 1, ...
    'usectfcomp',  1, ...
    'usessp',      1, ...
    'freq',        [], ...
    'baseline',    []);

% Process: Average: By trial group (folder average)
sFilesAvg = bst_process('CallProcess', 'process_average', sFilesEpochs, [], ...
    'avgtype',       5, ...  % By trial group (folder average)
    'avg_func',      1, ...  % Arithmetic average:  mean(x)
    'weighted',      0, ...
    'keepevents',    0);

% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
    'target',   5, ...  % Recordings time series
    'modality', 4, ...  % EEG
    'Comment',  'EEG ERP');
% Process: Snapshot: Recordings topography (contact sheet)
bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
    'target',   6, ...  % Recordings topography (ont time)
    'modality', 4, ...  % EEG
    'time',     Latency, ...
    'Comment',  'EEG ERP');
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
    'target',   5, ...  % Recordings time series
    'modality', 1, ...  % MEG (All)
    'Comment',  'MEG ERF');
% Process: Snapshot: Recordings topography (contact sheet)
bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
    'target',   6, ...  % Recordings topography (ont time)
    'modality', 1, ...  % MEG (All)
    'time',     Latency, ...
    'Comment',  'MEG ERF');

%% ===== NOISE COVARIANCE =====
% Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ...
    'baseline',       [-0.1, -0.01], ...
    'sensortypes',    'MEG, EEG', ...
    'target',         1, ...  % Noise covariance     (covariance over baseline time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'replacefile',    1);  % Replace

%% ===== FORWARD: FEM EEG DTI =====
% Process: Compute head model
bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ...
    'Comment',     'DUNEuro FEM EEG DTI', ...
    'sourcespace', 1, ...  % Cortex surface
    'meg',         1, ...  %
    'eeg',         4, ...  % DUNEuro FEM
    'ecog',        1, ...  %
    'seeg',        1, ...  %
    'duneuro',     struct(...
         'FemCond',             [], ...
         'FemSelect',           [1, 1, 1, 1, 1], ...
         'UseTensor',           1, ...
         'Isotropic',           1, ...
         'SrcShrink',           0, ...
         'SrcForceInGM',        0, ...
         'FemType',             'fitted', ...
         'SolverType',          'cg', ...
         'GeometryAdapted',     0, ...
         'Tolerance',           1e-08, ...
         'ElecType',            'normal', ...
         'MegIntorderadd',      0, ...
         'MegType',             'physical', ...
         'SolvSolverType',      'cg', ...
         'SolvPrecond',         'amg', ...
         'SolvSmootherType',    'ssor', ...
         'SolvIntorderadd',     0, ...
         'DgSmootherType',      'ssor', ...
         'DgScheme',            'sipg', ...
         'DgPenalty',           20, ...
         'DgEdgeNormType',      'houston', ...
         'DgWeights',           1, ...
         'DgReduction',         1, ...
         'SolPostProcess',      1, ...
         'SolSubstractMean',    0, ...
         'SolSolverReduction',  1e-10, ...
         'SrcModel',            'venant', ...
         'SrcIntorderadd',      0, ...
         'SrcIntorderadd_lb',   2, ...
         'SrcNbMoments',        3, ...
         'SrcRefLen',           20, ...
         'SrcWeightExp',        1, ...
         'SrcRelaxFactor',      6, ...
         'SrcMixedMoments',     1, ...
         'SrcRestrict',         1, ...
         'SrcInit',             'closest_vertex', ...
         'BstSaveTransfer',     0, ...
         'BstEegTransferFile',  'eeg_transfer.dat', ...
         'BstMegTransferFile',  'meg_transfer.dat', ...
         'BstEegLfFile',        'eeg_lf.dat', ...
         'BstMegLfFile',        'meg_lf.dat', ...
         'UseIntegrationPoint', 1, ...
         'EnableCacheMemory',   0, ...
         'MegPerBlockOfSensor', 0), ...
    'channelfile', '');
% Process: Compute sources [2018]
sSrcEegFemDti = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
    'output',  1, ...  % Kernel only: shared
    'inverse', struct(...
         'Comment',        'Dipoles: EEG FEM DTI', ...
         'InverseMethod',  'gls', ...
         'InverseMeasure', 'performance', ...
         'SourceOrient',   {{'free'}}, ...
         'Loose',          0.2, ...
         'UseDepth',       1, ...
         'WeightExp',      0.5, ...
         'WeightLimit',    10, ...
         'NoiseMethod',    'median', ...
         'NoiseReg',       0.1, ...
         'SnrMethod',      'rms', ...
         'SnrRms',         1e-06, ...
         'SnrFixed',       3, ...
         'ComputeKernel',  1, ...
         'DataTypes',      {{'EEG'}}));
% Process: Dipole scanning
sDipEegFemDti = bst_process('CallProcess', 'process_dipole_scanning', sSrcEegFemDti, [], ...
    'timewindow', [Latency, Latency], ...
    'scouts',     {});

%% ===== FORWARD: FEM MEG DTI =====
% Process: Compute head model
bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ...
    'Comment',     'DUNEuro FEM MEG DTI', ...
    'sourcespace', 1, ...  % Cortex surface
    'meg',         5, ...  % DUNEuro FEM
    'eeg',         1, ...  %
    'ecog',        1, ...  %
    'seeg',        1, ...  %
    'duneuro',     struct(...
         'FemCond',             [], ...
         'FemSelect',           [1, 1, 1, 0, 0], ...
         'UseTensor',           1, ...
         'Isotropic',           1, ...
         'SrcShrink',           0, ...
         'SrcForceInGM',        0, ...
         'FemType',             'fitted', ...
         'SolverType',          'cg', ...
         'GeometryAdapted',     0, ...
         'Tolerance',           1e-08, ...
         'ElecType',            'normal', ...
         'MegIntorderadd',      0, ...
         'MegType',             'physical', ...
         'SolvSolverType',      'cg', ...
         'SolvPrecond',         'amg', ...
         'SolvSmootherType',    'ssor', ...
         'SolvIntorderadd',     0, ...
         'DgSmootherType',      'ssor', ...
         'DgScheme',            'sipg', ...
         'DgPenalty',           20, ...
         'DgEdgeNormType',      'houston', ...
         'DgWeights',           1, ...
         'DgReduction',         1, ...
         'SolPostProcess',      1, ...
         'SolSubstractMean',    0, ...
         'SolSolverReduction',  1e-10, ...
         'SrcModel',            'venant', ...
         'SrcIntorderadd',      0, ...
         'SrcIntorderadd_lb',   2, ...
         'SrcNbMoments',        3, ...
         'SrcRefLen',           20, ...
         'SrcWeightExp',        1, ...
         'SrcRelaxFactor',      6, ...
         'SrcMixedMoments',     1, ...
         'SrcRestrict',         1, ...
         'SrcInit',             'closest_vertex', ...
         'BstSaveTransfer',     0, ...
         'BstEegTransferFile',  'eeg_transfer.dat', ...
         'BstMegTransferFile',  'meg_transfer.dat', ...
         'BstEegLfFile',        'eeg_lf.dat', ...
         'BstMegLfFile',        'meg_lf.dat', ...
         'UseIntegrationPoint', 1, ...
         'EnableCacheMemory',   0, ...
         'MegPerBlockOfSensor', 0), ...
    'channelfile', '');
% Process: Compute sources [2018]
sSrcMegFemDti = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
    'output',  1, ...  % Kernel only: shared
    'inverse', struct(...
         'Comment',        'Dipoles: MEG FEM DTI', ...
         'InverseMethod',  'gls', ...
         'InverseMeasure', 'performance', ...
         'SourceOrient',   {{'free'}}, ...
         'Loose',          0.2, ...
         'UseDepth',       1, ...
         'WeightExp',      0.5, ...
         'WeightLimit',    10, ...
         'NoiseMethod',    'median', ...
         'NoiseReg',       0.1, ...
         'SnrMethod',      'rms', ...
         'SnrRms',         1e-06, ...
         'SnrFixed',       3, ...
         'ComputeKernel',  1, ...
         'DataTypes',      {{'MEG GRAD', 'MEG MAG'}}));
% Process: Dipole scanning
sDipMegFemDti = bst_process('CallProcess', 'process_dipole_scanning', sSrcMegFemDti, [], ...
    'timewindow', [Latency, Latency], ...
    'scouts',     {});

%% ===== FORWARD: FEM EEG ISO =====
% Get updated subject structure
sSubject = bst_get('Subject', iSubject);
% Get FEM file
FemFile = file_fullpath(sSubject.Surface(sSubject.iFEM).FileName);
% Remove tensors
process_fem_tensors('ClearTensors', FemFile);
% Process: Compute head model
bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ...
    'Comment',     'DUNEuro FEM EEG ISO', ...
    'sourcespace', 1, ...  % Cortex surface
    'meg',         1, ...  %
    'eeg',         4, ...  % DUNEuro FEM
    'ecog',        1, ...  %
    'seeg',        1, ...  %
    'duneuro',     struct(...
         'FemCond',             [0.14, 0.33, 1.79, 0.008, 0.43], ...
         'FemSelect',           [1, 1, 1, 1, 1], ...
         'UseTensor',           0, ...
         'Isotropic',           1, ...
         'SrcShrink',           0, ...
         'SrcForceInGM',        0, ...
         'FemType',             'fitted', ...
         'SolverType',          'cg', ...
         'GeometryAdapted',     0, ...
         'Tolerance',           1e-08, ...
         'ElecType',            'normal', ...
         'MegIntorderadd',      0, ...
         'MegType',             'physical', ...
         'SolvSolverType',      'cg', ...
         'SolvPrecond',         'amg', ...
         'SolvSmootherType',    'ssor', ...
         'SolvIntorderadd',     0, ...
         'DgSmootherType',      'ssor', ...
         'DgScheme',            'sipg', ...
         'DgPenalty',           20, ...
         'DgEdgeNormType',      'houston', ...
         'DgWeights',           1, ...
         'DgReduction',         1, ...
         'SolPostProcess',      1, ...
         'SolSubstractMean',    0, ...
         'SolSolverReduction',  1e-10, ...
         'SrcModel',            'venant', ...
         'SrcIntorderadd',      0, ...
         'SrcIntorderadd_lb',   2, ...
         'SrcNbMoments',        3, ...
         'SrcRefLen',           20, ...
         'SrcWeightExp',        1, ...
         'SrcRelaxFactor',      6, ...
         'SrcMixedMoments',     1, ...
         'SrcRestrict',         1, ...
         'SrcInit',             'closest_vertex', ...
         'BstSaveTransfer',     0, ...
         'BstEegTransferFile',  'eeg_transfer.dat', ...
         'BstMegTransferFile',  'meg_transfer.dat', ...
         'BstEegLfFile',        'eeg_lf.dat', ...
         'BstMegLfFile',        'meg_lf.dat', ...
         'UseIntegrationPoint', 1, ...
         'EnableCacheMemory',   0, ...
         'MegPerBlockOfSensor', 0), ...
    'channelfile', '');
% Process: Compute sources [2018]
sSrcEegFemIso = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
    'output',  1, ...  % Kernel only: shared
    'inverse', struct(...
         'Comment',        'Dipoles: EEG FEM ISO', ...
         'InverseMethod',  'gls', ...
         'InverseMeasure', 'performance', ...
         'SourceOrient',   {{'free'}}, ...
         'Loose',          0.2, ...
         'UseDepth',       1, ...
         'WeightExp',      0.5, ...
         'WeightLimit',    10, ...
         'NoiseMethod',    'median', ...
         'NoiseReg',       0.1, ...
         'SnrMethod',      'rms', ...
         'SnrRms',         1e-06, ...
         'SnrFixed',       3, ...
         'ComputeKernel',  1, ...
         'DataTypes',      {{'EEG'}}));
% Process: Dipole scanning
sDipEegFemIso = bst_process('CallProcess', 'process_dipole_scanning', sSrcEegFemIso, [], ...
    'timewindow', [Latency, Latency], ...
    'scouts',     {});

%% ===== FORWARD: FEM MEG ISO =====
% Process: Compute head model
bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ...
    'Comment',     'DUNEuro FEM MEG ISO', ...
    'sourcespace', 1, ...  % Cortex surface
    'meg',         5, ...  % DUNEuro FEM
    'eeg',         1, ...  %
    'ecog',        1, ...  %
    'seeg',        1, ...  %
    'duneuro',     struct(...
         'FemCond',             [0.14, 0.33, 1.79, 0.008, 0.43], ...
         'FemSelect',           [1, 1, 1, 0, 0], ...
         'UseTensor',           0, ...
         'Isotropic',           1, ...
         'SrcShrink',           0, ...
         'SrcForceInGM',        0, ...
         'FemType',             'fitted', ...
         'SolverType',          'cg', ...
         'GeometryAdapted',     0, ...
         'Tolerance',           1e-08, ...
         'ElecType',            'normal', ...
         'MegIntorderadd',      0, ...
         'MegType',             'physical', ...
         'SolvSolverType',      'cg', ...
         'SolvPrecond',         'amg', ...
         'SolvSmootherType',    'ssor', ...
         'SolvIntorderadd',     0, ...
         'DgSmootherType',      'ssor', ...
         'DgScheme',            'sipg', ...
         'DgPenalty',           20, ...
         'DgEdgeNormType',      'houston', ...
         'DgWeights',           1, ...
         'DgReduction',         1, ...
         'SolPostProcess',      1, ...
         'SolSubstractMean',    0, ...
         'SolSolverReduction',  1e-10, ...
         'SrcModel',            'venant', ...
         'SrcIntorderadd',      0, ...
         'SrcIntorderadd_lb',   2, ...
         'SrcNbMoments',        3, ...
         'SrcRefLen',           20, ...
         'SrcWeightExp',        1, ...
         'SrcRelaxFactor',      6, ...
         'SrcMixedMoments',     1, ...
         'SrcRestrict',         1, ...
         'SrcInit',             'closest_vertex', ...
         'BstSaveTransfer',     0, ...
         'BstEegTransferFile',  'eeg_transfer.dat', ...
         'BstMegTransferFile',  'meg_transfer.dat', ...
         'BstEegLfFile',        'eeg_lf.dat', ...
         'BstMegLfFile',        'meg_lf.dat', ...
         'UseIntegrationPoint', 1, ...
         'EnableCacheMemory',   0, ...
         'MegPerBlockOfSensor', 0), ...
    'channelfile', '');
% Process: Compute sources [2018]
sSrcMegFemIso = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
    'output',  1, ...  % Kernel only: shared
    'inverse', struct(...
         'Comment',        'Dipoles: MEG FEM ISO', ...
         'InverseMethod',  'gls', ...
         'InverseMeasure', 'performance', ...
         'SourceOrient',   {{'free'}}, ...
         'Loose',          0.2, ...
         'UseDepth',       1, ...
         'WeightExp',      0.5, ...
         'WeightLimit',    10, ...
         'NoiseMethod',    'median', ...
         'NoiseReg',       0.1, ...
         'SnrMethod',      'rms', ...
         'SnrRms',         1e-06, ...
         'SnrFixed',       3, ...
         'ComputeKernel',  1, ...
         'DataTypes',      {{'MEG GRAD', 'MEG MAG'}}));
% Process: Dipole scanning
sDipMegFemIso = bst_process('CallProcess', 'process_dipole_scanning', sSrcMegFemIso, [], ...
    'timewindow', [Latency, Latency], ...
    'scouts',     {});

%% ===== FORWARD: BEM EEG =====
% Process: Compute head model
bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ...
    'Comment',     'OpenMEEG BEM EEG', ...
    'sourcespace', 1, ...  % Cortex surface
    'meg',         1, ...  %
    'eeg',         3, ...  % OpenMEEG BEM
    'ecog',        1, ...  %
    'seeg',        1, ...  %
    'openmeeg',    struct(...
         'BemSelect',    [1, 1, 1], ...
         'BemCond',      [1, 0.0125, 1], ...
         'BemNames',     {{'Scalp', 'Skull', 'Brain'}}, ...
         'BemFiles',     {{}}, ...
         'isAdjoint',    0, ...
         'isAdaptative', 1, ...
         'isSplit',      0, ...
         'SplitLength',  4000), ...
    'channelfile', '');
% Process: Compute sources [2018]
sSrcEegBem = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
    'output',  1, ...  % Kernel only: shared
    'inverse', struct(...
         'Comment',        'Dipoles: EEG BEM', ...
         'InverseMethod',  'gls', ...
         'InverseMeasure', 'performance', ...
         'SourceOrient',   {{'free'}}, ...
         'Loose',          0.2, ...
         'UseDepth',       1, ...
         'WeightExp',      0.5, ...
         'WeightLimit',    10, ...
         'NoiseMethod',    'median', ...
         'NoiseReg',       0.1, ...
         'SnrMethod',      'rms', ...
         'SnrRms',         1e-06, ...
         'SnrFixed',       3, ...
         'ComputeKernel',  1, ...
         'DataTypes',      {{'EEG'}}));
% Process: Dipole scanning
sDipEegBem = bst_process('CallProcess', 'process_dipole_scanning', sSrcEegBem, [], ...
    'timewindow', [Latency, Latency], ...
    'scouts',     {});

%% ===== FORWARD: OS MEG =====
% Process: Compute head model
bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ...
    'Comment',     'OS MEG', ...
    'sourcespace', 1, ...  % Cortex surface
    'meg',         3, ...  % Overlapping spheres
    'eeg',         1, ...  %
    'ecog',        1, ...  %
    'seeg',        1, ...  %
    'channelfile', '');
% Process: Compute sources [2018]
sSrcMegOs = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
    'output',  1, ...  % Kernel only: shared
    'inverse', struct(...
         'Comment',        'Dipoles: MEG OS', ...
         'InverseMethod',  'gls', ...
         'InverseMeasure', 'performance', ...
         'SourceOrient',   {{'free'}}, ...
         'Loose',          0.2, ...
         'UseDepth',       1, ...
         'WeightExp',      0.5, ...
         'WeightLimit',    10, ...
         'NoiseMethod',    'median', ...
         'NoiseReg',       0.1, ...
         'SnrMethod',      'rms', ...
         'SnrRms',         1e-06, ...
         'SnrFixed',       3, ...
         'ComputeKernel',  1, ...
         'DataTypes',      {{'MEG GRAD', 'MEG MAG'}}));
% Process: Dipole scanning
sDipMegOs = bst_process('CallProcess', 'process_dipole_scanning', sSrcMegOs, [], ...
    'timewindow', [Latency, Latency], ...
    'scouts',     {});

%% ===== MERGE DIPOLES =====
% sDipEegFemDti, sDipMegFemDti, sDipEegFemIso, sDipMegFemIso, sDipEegBem, sDipMegOs
% FEM DTI: EEG vs MEG
% sDipFemDti = dipoles_merge({sDipEegFemDti.FileName, sDipMegFemDti.FileName});

% Save and display report
ReportFile = bst_report('Save', []);
bst_report('Open', ReportFile);


text:
function x = p20_start ( option, nvar )

%*****************************************************************************80
%
%% P20_START returns a starting point for problem 20.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Output, real X(NVAR), the starting point.
%
  hold_index = p20_i4_get ( 'hold_index' );
  hold_value = p20_r8_get ( 'hold_value' );

  if ( hold_index == 1 )
    l = hold_value;
    theta = pi / 8.0;
  elseif ( hold_index == 2 )
    l = 0.25;
    theta = hold_value;
  else
    l = 0.25;
    theta = pi / 8.0;
  end

  [ lambda, mu ] = p20_setup ( l, theta );

  x = [ l, theta, lambda, mu ]';

  return
end


text:
function D = calc_dist_matrix(M, samples)

if nargin==1 || isempty(samples)
    samples = 1:M.n;
end

march = fastmarchmex('init', int32(M.TRIV-1), double(M.VERT(:,1)), double(M.VERT(:,2)), double(M.VERT(:,3)));

D = zeros(length(samples));

for i=1:length(samples)
%     fprintf('(%d/%d)\n', i, length(samples));
    source = inf(M.n,1);
    source(samples(i)) = 0;
    d = fastmarchmex('march', march, double(source));
    D(:,i) = d(samples);
end

fastmarchmex('deinit', march);

D = 0.5*(D+D');

end

author: OshriHalimi


text:
function varargout = drawVector(pos, vect, varargin)
%DRAWVECTOR Draw vector at a given position
%
%   drawVector(POS, VECT)
%   POS should be a N-by-2 or N-by-3 array containing position of vector
%   origins, and VECT should be a N-by-2 or N-by-3 array containing the
%   direction of the vectors.
%
%   Example
%     figure; hold on;
%     drawVector([1 2], [3 2]);
%     drawVector([1 2], [-2 3]);
%     axis equal;
%
%   See also
%     quiver, drawVector3d
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2013-03-18,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2013 INRA - Cepia Software Platform.

% check input dimension
nd = size(pos, 2);
if size(vect, 2) ~= nd
    error('input vector and position must have same dimension');
end

if nd == 2
    % Display 2D vectors
    h = quiver(pos(:, 1), pos(:, 2), vect(:, 1), vect(:, 2), 0, varargin{:});

elseif nd == 3
    % Display 3D vectors
    h = quiver3(pos(:, 1), pos(:, 2), pos(:, 3), ...
        vect(:, 1), vect(:, 2), vect(:, 3), 0, varargin{:});

else
    error('Can not display vectors of dimension > 3');
end

% format output
if nargout > 0
    varargout{1} = h;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [Rtrunc,q,u,v]=nmt_svdtrunc(R,signalspace)
% [Rtrunc,q,u,v] = nmt_svdtrunc(R,signalspace)
%
% Allows user to reject undesired SVD components of an arbitrary matrix R.
% It can be used to, e.g., reduce the matrix's rank.
%
% signalspace should be vector of which components to include, e.g. [1 2]
% for the equivalent of cfg.reducerank=2

[u,q,v]=svd(R,'econ');

notsignalspace = setdiff(1:length(q),signalspace);
q(notsignalspace,notsignalspace)=0;
Rtrunc = u*q*v';

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef DataReader < handle
methods
    function data = next(this)
        data = this.nextImpl;
    end
    function data = previous(this)
        data = this.previousImpl;
    end
    function data = current(this)
        data = this.currentImpl;
    end
end

methods(Abstract, Access=protected)
    data = nextImpl(this);
    data = previousImpl(this);
    data = currentImpl(this);
end
end

repo: ekfmonoslam


text:
function [ cfg ] = ft_trialfun_balert( cfg )
%ft_definetrialBAlert: Extract trials B-Alert data from CSV file.
%   FieldTrip cannot yet interpret the event markers from B-Alert data.
%   Therefore, it is necessary to have B-Alert LAB. This is (paid) software
%   from Advanced Brain Monitoring, in which you extract the evetmakers
%   using the function: readevents(*.Events.edf,*.Signals.Raw.edf) . This
%   function returns a *.csv file in which, for every epoch (second)
    fid = fopen(cfg.dataset,'rt');
    frmt = '%d%d%s%d%d%s%d%d%d';
    events = textscan(fid,frmt,'Delimiter',',','CollectOutput',true, 'HeaderLines',1);
    len = size(events{1,1},1);
    eventMarkers = str2num(cell2mat(events{1,4}(2:len-1,:))); % event marker data
    epochTime = events{1,1}(2:len-1,:); %epoched time
    % Resample to 1024 Hz and save as sample numbers
    cfg.event = struct(...
        'type',repmat('STATUS',size(epochTime,1)),...
        'sample',epochTime(:,1)*1024 + round(epochTime(:,2)/256*1024),...
        'value',eventMarkers);

    % Find stimulus
    stimIndex = zeros(size(eventMarkers,1),1);
    for i = 1:length(cfg.eventvalue)
        stimIndex = stimIndex + (eventMarkers==cfg.eventvalue(i));
    end
    stimIndex = logical(stimIndex);
    cfg.trl(:,1) = cfg.event.sample(stimIndex)-cfg.prestim*cfg.newfs;
    cfg.trl(:,2) = cfg.event.sample(stimIndex)+cfg.poststim*cfg.newfs;
    cfg.trl(:,3) = eventMarkers(stimIndex);
end

author: spm


text:
%CONCATM Fixed-cell mapping concatenating cell array of datasets or mappings
%
%   W = V*CONCATM
%   B = A*CONCATM
%
% INPUT
%   V  Cell array of mappings
%   A  Cell array of datasets
%
% OUTPUT
%   W  Concatenation of mappings: W = [V{:}]
%   A  Concatenation of datasets: B = [A{:}]
%
% DESCRIPTION
% This routine is an exception of the general treatment of cell arrays and
% mappings in PRTools. As a rule the elements of a cell array are processed
% one by one by a mapping. This routine takes the entire cell array and
% concatenates its elements.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com

function out = concatm(argin)

if nargin == 0
  out = prmapping(mfilename,'fixed_cell');
else
  if ~iscell(argin)
    error('Cell array expected')
  else
    out = [];
    for j=1:size(argin,1)
      out = [out;[argin{j,:}]];
    end
  end
end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function pt = steinerPoint(varargin)
%STEINERPOINT Compute steiner point (weighted centroid) of a polygon
%
%   PT = steinerPoint(POINTS);
%   PT = steinerPoint(PTX, PTY);
%   Computes steiner point of a polygon defined by POINTS. POINTS is a
%   [N*2] array of double.
%
%   The steiner point is computed the same way as the polygon centroid,
%   except that a weight depending on the angle is given to each vertex.
%
%   See also:
%   polygons2d, polygonArea, polygonCentroid, drawPolygon
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 11/11/2004.
%

if nargin==1
    var = varargin{1};
    px = var(:,1);
    py = var(:,2);
elseif nargin==2
    px = varargin{1};
    py = varargin{2};
end

% Algorithme P. Bourke
sx = 0;
sy = 0;
N = length(px);
for i=1:N-1
    sx = sx + (px(i)+px(i+1))*(px(i)*py(i+1) - px(i+1)*py(i));
    sy = sy + (py(i)+py(i+1))*(px(i)*py(i+1) - px(i+1)*py(i));
end
sx = sx + (px(N)+px(1))*(px(N)*py(1) - px(1)*py(N));
sy = sy + (py(N)+py(1))*(px(N)*py(1) - px(1)*py(N));

pt = [sx sy]/6/polygonArea(px, py);


text:
function test_cases = testSubfunctions
%testSubfunctions Contains two passing subfunction tests

%   Steven L. Eddins
%   Copyright 2008 The MathWorks, Inc.

findSubfunctionTests;

function testSub1

function testSub2

author: ricedsp