text:
function [Results, OPTIONS, HeadModel] = bst_wmne(HeadModel,OPTIONS)
% BST_WMNE: Compute the whitened and weighted minimum-norm operator (wMNE imaging kernel)
%
% USAGE:  [Results,OPTIONS] = bst_wmne(HeadModel, OPTIONS) : Compute mininum operator
%                   OPTIONS = bst_wmne()                   : Return default options
%
% DESCRIPTION:
%     This program computes the whitened and weighted minimum-norm operator,
%     (the wMNE imaging kernel), which is used to compute whitened
%     and weighted minimum-norm estimates (MNE).
%         (e.g., J=wMNEoperator*B; where B is the unwhitened data).
%     It can also compute the whitened and noise-normalized dynamic
%     statistical parametric mapping (dSPM) inverse operator, and/or the
%     whitened standardized low resolution brain electromagnetic tomography
%     (sLORETA) inverse operator, which are used to compute whitened source
%     activity dSPM and sLORETA maps.
%         (e.g., S_dSPM=dSPMoperator*B; where B is the unwhitened data).
%         (e.g., S_sLORETA=sLORETAoperator*B; where B is the unwhitened data).
%
%     The function was written with the goal of providing some of the same
%     functionality of the MNE software written by Matti Hamalainen, but no
%     guarantees are made that it performs all computations in the same
%     exact way. It also provides some functionalities not available in the
%     MNE software.
%
% INPUTS:
%    - HeadModel: Array of Brainstorm head model structures
%         |- Gain       : Forward field matrix for all the channels (unconstrained source orientations)
%         |- GridOrient : Dipole orientation matrix
%         |- area       : Vector with the areas (or possibly volumes) associated with the vertices of the source space.
%    - OPTIONS: structure
%         |- NoiseCov      : NoiseCov is the noise covariance matrix.
%         |- ChannelTypes  : Type of each channel (for each row of the Leadfield and the NoiseCov matrix)
%         |- InverseMethod : {'wmne', 'dspm', 'sloreta'}
%         |- SourceOrient  : String or a cell array of strings specifying the type of orientation constraints for each HeadModel (default: 'fixed')
%         |- SNR        : Signal-to noise ratio defined as in MNE (default: 3).
%         |- diagnoise  : Flag to discard off-diagonal elements of NoiseCov (assuming heteroscedastic uncorrelated noise) (default: 0)
%         |- loose      : Value that weights the source variances of the dipole components defining the tangent space of the cortical surfaces (default: []).
%         |- depth      : Flag to do depth weighting (default: 1).
%         |- weightexp  : Order of the depth weighting. {0=no, 1=full normalization, default=0.8}
%         |- weightlimit: Maximal amount depth weighting (default: 10).
%         |- magreg     : Amount of regularization of the magnetometer noise covariance matrix
%         |- gradreg    : Amount of regularization of the gradiometer noise covariance matrix.
%         |- eegreg     : Amount of regularization of the EEG noise covariance matrix.
%         |- ecogreg    : Amount of regularization of the ECOG noise covariance matrix.
%         |- seegreg    : Amount of regularization of the SEEG noise covariance matrix.
%         |- fMRI       : Vector of fMRI values are the source points.
%         |- fMRIthresh : fMRI threshold. The source variances of source points with OPTIONS.fMRI smaller
%         |               than fMRIthresh will be multiplied by OPTIONS.fMRIoff.
%         |- fMRIoff    : Weight assigned to non-active source points according to fMRI and fMRIthresh.
%
% OUTPUTS:
%    - Results : structure with the wMNE inverse operator and possibly the dSPM
%                and/or sLORETA inverse operators, and other information:

% NOTES:
%     - More leadfield matrices can be used: the solution will combine all
%       leadfield matrices appropriately.
%
%     - This leadfield structure allows to combine surface and volume
%       source spaces with and without dipole orientation constraints,and
%       with and without area or volumetric current density computations.
%       If using a single sphere headmodel, the silent radial
%       component could be eliminated using the SVD (e.g., use
%       bst_remove_silent.m).
%
%     - NoisCov: This should be computed from the pre-stimulus period for
%       averaged ERF data (e.g., using MNE), or from an empty room recording
%       for unaveraged spontaneous or induced data.
%
%     - Orientation constrains for dipoles (.SourceOrient field)
%         - "fixed"     : Dipoles constrained to point normal to cortical surfaces
%         - "free"      : No constraints, dipoles point in x,y,z directions
%         - "loose"     : Source variances of dipole components pointing tangentially to the cortical surfaces are multipled by OPTIONS.loose
%         - "truncated" : An SVD of the gain matrix for each source point is
%                         used to remove the dipole component with least variance, which for
%                         the Single Sphere Head Model, corresponds to the radialsilent component).
%       => For dealing with multiple source spaces with different types of orientation constraints use for example,
%          OPTION.SourceOrient{1}='fixed';
%          OPTION.SourceOrient{2}='free';
%          This has to correspond with the HeadModel Structure
%
%     - The output HeadModel param is used here in return to save LOTS of memory in the bst_wmne function,
%       event if it seems to be absolutely useless. Having a parameter in both input and output have the
%       effect in Matlab of passing them "by referece". So please, do NOT remove it from the function description
%
%     - sLORETA: Output values are multiplied by 1e12 for display in Brainstorm (time series and cortical maps).

% @=============================================================================
% 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.
% =============================================================================@
%
% Copyright (C) 2010 - Rey Rene Ramirez
%
% Authors:  Rey Rene Ramirez, Ph.D.   e-mail: rrramirez at mcw.edu
%           Francois Tadel, 2010-2013

%% ===== DEFINE DEFAULT OPTIONS =====
Def_OPTIONS.NoiseCov      = [];
Def_OPTIONS.InverseMethod = 'wmne';
Def_OPTIONS.ChannelTypes  = {};
Def_OPTIONS.SNR           = 3;
Def_OPTIONS.diagnoise     = 0;
%Def_OPTIONS.SourceOrient= {'free'};
Def_OPTIONS.SourceOrient= {'fixed'};
Def_OPTIONS.loose       = 0.2;
Def_OPTIONS.depth       = 1;
Def_OPTIONS.weightexp   = 0.5;
Def_OPTIONS.weightlimit = 10;
Def_OPTIONS.regnoise    = 1;
Def_OPTIONS.magreg      = .1;
Def_OPTIONS.gradreg     = .1;
Def_OPTIONS.eegreg      = .1;
Def_OPTIONS.ecogreg     = .1;
Def_OPTIONS.seegreg     = .1;
Def_OPTIONS.fMRI        = [];
Def_OPTIONS.fMRIthresh  = [];
Def_OPTIONS.fMRIoff     = 0.1;
Def_OPTIONS.pca         = 1;
% Return the default options
if (nargin == 0)
    Results = Def_OPTIONS;
    return
end
% Make the default for all the leadfields
numL = size(HeadModel,2);
Def_OPTIONS.SourceOrient = repmat(Def_OPTIONS.SourceOrient, [1 numL]);
% Copy default options to OPTIONS structure (do not replace defined values)
OPTIONS = struct_copy_fields(OPTIONS, Def_OPTIONS, 0);

%% ===== CHECK FOR INVALID VALUES =====
% Detect if the input noise covariance matrix is or should be diagonal
C_noise = OPTIONS.NoiseCov;
variances = diag(C_noise);
if isequal(C_noise, diag(variances))
    OPTIONS.diagnoise = 1;
    disp('wMNE> Detected diagonal noise covariance: setting diagnoise to 1');
end
% If OPTIONS.diagnoise is 1, then OPTIONS.pca=0
if OPTIONS.diagnoise
    OPTIONS.pca=0;
    disp('wMNE> If using diagonal noise covariance, PCA option should be off. Setting PCA option off.')
end
if isempty(OPTIONS.NoiseCov)
    error('You need to input the noise covariance in the NoiseCov field of OPTIONS.')
end
if (numL ~= length(OPTIONS.SourceOrient))
    error('The number of elements in the HeadModel structure should equal the length of the cell array OPTIONS.SourceOrient.')
end
if ~isempty(OPTIONS.loose) && (OPTIONS.loose>=1 || OPTIONS.loose<=0)
    error('loose value should be smaller than 1 and bigger than 0, or empty for no loose orientations.')
end
if OPTIONS.weightexp>1 || OPTIONS.weightexp<0
    error('weightexp should be a scalar between 0 and 1')
end
if OPTIONS.magreg>1 || OPTIONS.magreg<0
    error('magreg should be a scalar between 0 and 1')
end
if OPTIONS.eegreg>1 || OPTIONS.eegreg<0
    error('eegreg should be a scalar between 0 and 1')
end
if OPTIONS.ecogreg>1 || OPTIONS.ecogreg<0
    error('ecogreg should be a scalar between 0 and 1')
end
if OPTIONS.seegreg>1 || OPTIONS.seegreg<0
    error('seegreg should be a scalar between 0 and 1')
end

%% ===== NOISE COVARIANCE RANK =====
% Get indices of MEG and EEG channels
iMeg  = find(strncmpi(OPTIONS.ChannelTypes,'MEG',3));
iEeg  = find(strncmpi(OPTIONS.ChannelTypes,'EEG',3));
iEcog = find(strncmpi(OPTIONS.ChannelTypes,'ECOG',3));
iSeeg = find(strncmpi(OPTIONS.ChannelTypes,'SEEG',3));
% Diagonal noisecov
if OPTIONS.diagnoise
    C_noise = diag(variances);
    rnkC_noise_meg = length(iMeg);
    rnkC_noise_eeg = length(iEeg);
    rnkC_noise_ecog = length(iEcog);
    rnkC_noise_seeg = length(iSeeg);
    disp('wMNE> Setting off diagonal elements of the noise covariance to zero.');
    disp(['wMNE> Rank of noise covariance is ' num2str(size(C_noise,1))]);
% Full noisecov
else
    % Estimate noise covariance matrix rank separately for sensor types
    if ~isempty(iMeg)
        rnkC_noise_meg = rank(single(C_noise(iMeg,iMeg))); % Rey added this. Separate rank of MEG. 3/23/11
        disp(['wMNE> Rank of MEG part of noise covariance is ' num2str(rnkC_noise_meg)]);
    end
    if ~isempty(iEeg)
        rnkC_noise_eeg = rank(single(C_noise(iEeg,iEeg))); % Rey added this. Separate rank of EEG. 3/23/11
        disp(['wMNE> Rank of EEG part of noise covariance is ' num2str(rnkC_noise_eeg)]);
    end
    if ~isempty(iEcog)
        rnkC_noise_ecog = rank(single(C_noise(iEcog,iEcog)));  % FT added 21-Feb-13
        disp(['wMNE> Rank of ECOG part of noise covariance is ' num2str(rnkC_noise_ecog)]);
    end
    if ~isempty(iSeeg)
        rnkC_noise_seeg = rank(single(C_noise(iSeeg,iSeeg)));  % FT added 21-Feb-13
        disp(['wMNE> Rank of SEEG part of noise covariance is ' num2str(rnkC_noise_seeg)]);
    end
    % Sets off-diagonal terms to zero. Rey added this. 3/23/11
    C_noise_new = 0 * C_noise;
    if ~isempty(iMeg)
        C_noise_new(iMeg,iMeg) = C_noise(iMeg,iMeg);
    end
    if ~isempty(iEeg)
        C_noise_new(iEeg,iEeg) = C_noise(iEeg,iEeg);
    end
    if ~isempty(iEcog)
        C_noise_new(iEcog,iEcog) = C_noise(iEcog,iEcog);
    end
    if ~isempty(iSeeg)
        C_noise_new(iSeeg,iSeeg) = C_noise(iSeeg,iSeeg);
    end
    C_noise = C_noise_new;
end

%% ===== REGULARIZE NOISE COVARIANCE MATRIX =====
% Only if option is selected
if OPTIONS.regnoise
    listTypes = unique(OPTIONS.ChannelTypes);
    % Loop on all the required data types (MEG MAG, MEG GRAD, EEG)
    for iType = 1:length(listTypes)
        % Get channel indices
        iChan = find(strcmpi(OPTIONS.ChannelTypes, listTypes{iType}));
        % Regularize noise covariance matrix
        switch listTypes{iType}
            case 'MEG GRAD', reg = OPTIONS.gradreg;
            case 'MEG MAG',  reg = OPTIONS.magreg;
            case 'MEG',      reg = OPTIONS.gradreg;
            case 'EEG',      reg = OPTIONS.eegreg;
            case 'ECOG',     reg = OPTIONS.ecogreg;
            case 'SEEG',     reg = OPTIONS.seegreg;
        end
        C_noise(iChan,iChan) = C_noise(iChan,iChan) + (reg * mean(variances(iChan)) * eye(length(iChan)));
    end
end

% FT added 06-Nov-2014
% Make sure the noise covariance matrix is strictly symmetrical
% (the previous operations may cause rounding errors that make the matrix not exactly symmetrical)
if ~isequal(C_noise, C_noise')
    C_noise = tril(C_noise) + (tril(C_noise)' - diag(diag(C_noise)));
end

%% ===== WHITENING OPERATOR =====
% Rey added all of this, 3/23/11
% Modified FT 21-Feb-2013
% Whitening of each modality separately (MEG,EEG,ECOG,SEEG), which assumes
% zero covariance between them (i.e., a block diagonal noise covariance). This
% was recommended by Matti as EEG does not measure all the signals from the same
% environmental noise sources as MEG.
nChan = size(C_noise,2);
W = zeros(0, nChan);
if ~isempty(iMeg)
    W_meg = CalculateWhitener('MEG', C_noise, iMeg, rnkC_noise_meg, OPTIONS.pca);
    W_tmp = zeros(size(W_meg,1), nChan);
    W_tmp(:,iMeg) = W_meg;
    W = [W; W_tmp];
end
if ~isempty(iEeg)
    W_eeg = CalculateWhitener('EEG', C_noise, iEeg, rnkC_noise_eeg, OPTIONS.pca);
    W_tmp = zeros(size(W_eeg,1), nChan);
    W_tmp(:,iEeg) = W_eeg;
    W = [W; W_tmp];
end
if ~isempty(iEcog)
    W_ecog = CalculateWhitener('ECOG', C_noise, iEcog, rnkC_noise_ecog, OPTIONS.pca);
    W_tmp = zeros(size(W_ecog,1), nChan);
    W_tmp(:,iEcog) = W_ecog;
    W = [W; W_tmp];
end
if ~isempty(iSeeg)
    W_seeg = CalculateWhitener('SEEG', C_noise, iSeeg, rnkC_noise_seeg, OPTIONS.pca);
    W_tmp = zeros(size(W_seeg,1), nChan);
    W_tmp(:,iSeeg) = W_seeg;
    W = [W; W_tmp];
end
% Check for whitener integrity
if any(isnan(W(:))) || any(isinf(W(:)))
    error('Invalid noise covariance matrix.')
end
% Display rank of the whitener
rnkC_noise = size(W,1);
display(['wMNE> Total rank is ' num2str(rnkC_noise) '.'])

%% ===== PROCESSING LEAD FIELD MATRICES, WEIGHTS, AREAS & ORIENTATIONS =====
% Initializing.
spl = zeros(numL,1);
numdipcomp = spl;
for k = 1:numL
    sL = size(HeadModel(k).Gain, 2);
    switch OPTIONS.SourceOrient{k}
        case 'fixed',      numdipcomp(k)=1;
        case 'free',       numdipcomp(k)=3;
        case 'loose',      numdipcomp(k)=3;
        case 'truncated',  numdipcomp(k)=2;
    end
    spl(k) = (sL / 3) * numdipcomp(k); % This is a vector with the total number of dipole components per source space.
end
sspl = sum(spl); % This is the total number of dipole components across all source spaces.
L = zeros(rnkC_noise,sspl);
w = ones(sspl,1);
if isfield(HeadModel, 'area')
    areas = w;
else
    areas = [];
end
itangential = [];
start = 0;
Q_Cortex = cell(1,numL);
for k = 1:numL
    start = start + 1;
    endd = start + spl(k) - 1;
    Lk = HeadModel(k).Gain;
    HeadModel(k).Gain = [];
    %% ===== COMPUTE POWER =====
    szL = size(Lk);
    if OPTIONS.depth
        % display(['wMNE> Computing power of gain matrices at each source point for source space ' num2str(k) '.'])
        % Computing power
        wk = squeeze(sum(sum((reshape(Lk,[szL(1) 3 szL(2)/3])) .^2,1),2));
        wk = repmat(wk',[numdipcomp(k) 1]);
        wk = reshape(wk,[spl(k) 1]);
        w(start:endd) = wk;
        clear wk
    end
    switch OPTIONS.SourceOrient{k}
        case 'fixed'
            display('wMNE> Appying fixed dipole orientations.')
            Lk = bst_gain_orient(Lk,HeadModel(k).GridOrient);
        case 'free'
            display('wMNE> Using free dipole orientations. No constraints.')
        case 'loose'
            display('wMNE> Transforming lead field matrix to cortical coordinate system.')
            [Lk, Q_Cortex{k}] = bst_xyz2lf(Lk, HeadModel(k).GridOrient');
            % Getting indices for tangential dipoles.
            itangentialtmp = start:endd;
            itangentialtmp(1:3:end) = [];
            itangential = [itangential itangentialtmp];  %#ok<AGROW>
        case 'truncated'
            display('wMNE> Truncating the dipole component pointing in the direction with least variance (i.e., silent component for single sphere head model.')
            [Lk, Q_Cortex{k}] = bst_remove_silent(Lk);
    end
    %% ===== WHITEN LEAD FIELD MATRIX =====
    % Whiten lead field.
    % display(['wMNE> Whitening lead field matrix for source space ' num2str(k) '.'])
    Lk = W * Lk;
    if isfield(HeadModel(k),'area') && ~isempty(HeadModel(k).area)
        areav = HeadModel(k).area;
        areav = repmat(areav', [numdipcomp(k) 1]);
        areav = reshape(areav, [spl(k) 1]);
        areas(start:endd) = areav;
    end
    L(:,start:endd) = Lk;
    start = endd;
end
% Computing reciprocal of power.
w = 1 ./ w;
% Clear memory
clear Lk endd start itangentialtmp sL szL

%% ===== APPLY AREAS =====
if ~isempty(areas)
    display('wMNE> Applying areas to compute current source density.')
    areas = areas.^2;
    w = w .* areas;
end
clear areas

%% ===== APPLY DEPTH WEIGHTHING =====
if OPTIONS.depth
    % ===== APPLY WEIGHT LIMIT =====
    % Applying weight limit.
    % display('wMNE> Applying weight limit.')
    weightlimit2 = OPTIONS.weightlimit .^ 2;
    %limit=min(w(w>min(w)*weightlimit2));  % This is the Matti way.
    limit = min(w) * weightlimit2;  % This is the Rey way (robust to possible weight discontinuity).
    w(w>limit) = limit;

    % ===== APPLY WEIGHT EXPONENT =====
    % Applying weight exponent.
    % display('wMNE> Applying weight exponent.')
    w = w .^ OPTIONS.weightexp;
    clear limit weightlimit2
end

%% ===== APPLY LOOSE ORIENTATIONS =====
if ~isempty(itangential)
    display(['wMNE> Applying loose dipole orientations. Loose value of ' num2str(OPTIONS.loose) '.'])
    w(itangential) = w(itangential) * (OPTIONS.loose);
end

%% ===== APPLY fMRI PRIORS =====
% Apply fMRI Priors
if ~isempty(OPTIONS.fMRI)
    display('wMNE> Applying fMRI priors.')
    ifmri = (OPTIONS.fMRI < OPTIONS.fMRIthresh);
    w(ifmri) = w(ifmri) * OPTIONS.fMRIoff;
end

%% ===== ADJUSTING SOURCE COVARIANCE MATRIX =====
% Adjusting Source Covariance matrix to make trace of L*C_J*L' equal to number of sensors.
% display('wMNE> Adjusting source covariance matrix.')
C_J = speye(sspl, sspl);
C_J = spdiags(w, 0, C_J);
trclcl = trace(L * C_J * L');
C_J = C_J * (rnkC_noise / trclcl);
Rc = chol(C_J, 'lower');
LW = L * Rc;
clear C_J trclcl sspl itangential rnkC_noise

%% ===== SINGULAR VALUE DECOMPOSITION =====
% Set regularization parameter based on SNR
lambda2 = OPTIONS.SNR ^ (-2);
% Compute SVD.
% display('wMNE> Computing SVD of whitened and weighted lead field matrix.')
[U,S,V] = svd(LW,'econ');
s = diag(S);
ss = s ./ (s.^2 + lambda2);
display(sprintf('wMNE> Effective Max and Mean SNR: %.2f and %.2f', s(1)^2/lambda2, mean(s.^2)/lambda2));
clear LW lambda2 s

%% ===== WHITENED MNE IMAGING KERNEL =====
% Compute whitened MNE operator.
Kernel = Rc * V * diag(ss) * U';
clear Rc V ss U

%% ===== NORMALIZATIONS =====
start = 0;
for k = 1:numL
    start = start + 1;
    endd = start + spl(k) - 1;

    % ===== dSPM OPERATOR =====
    if strcmpi(OPTIONS.InverseMethod, 'dspm')
        dspmdiag = sum(Kernel(start:endd,:) .^2, 2);
        if (numdipcomp(k) == 1)
            dspmdiag = sqrt(dspmdiag);
        elseif (numdipcomp(k)==3 || numdipcomp(k)==2)
            dspmdiag = reshape(dspmdiag, [numdipcomp(k), spl(k)/numdipcomp(k)]);
            dspmdiag = sqrt(sum(dspmdiag)); % Taking trace and sqrt.
            dspmdiag = repmat(dspmdiag, [numdipcomp(k), 1]);
            dspmdiag = reshape(dspmdiag, [spl(k), 1]);
        end
        Kernel(start:endd,:) = bst_bsxfun(@rdivide, Kernel(start:endd,:), dspmdiag);

    % ===== sLORETA OPERATOR =====
    elseif strcmpi(OPTIONS.InverseMethod, 'sloreta')
        if (numdipcomp(k) == 1)
            sloretadiag = sqrt(sum(Kernel(start:endd,:) .* L(:,start:endd)', 2));
            Kernel(start:endd,:) = bst_bsxfun(@rdivide, Kernel(start:endd,:), sloretadiag);
        elseif (numdipcomp(k)==3 || numdipcomp(k)==2)
            for spoint = start:numdipcomp(k):endd
                R = Kernel(spoint:spoint+numdipcomp(k)-1,:) * L(:,spoint:spoint+numdipcomp(k)-1);
                SIR = sqrtm(pinv(R));
                Kernel(spoint:spoint+numdipcomp(k)-1,:) = SIR * Kernel(spoint:spoint+numdipcomp(k)-1,:);
            end
        end
    end

    % ===== LOOSE ORIENTATION: RE-ORIENT COMPONENTS =====
    % WARNING: Changing the orientations of the dipoles from MNE to Brainstorm
    % => Make "Unconstrained" and "loose"/"truncated" models directly comparable
    if ~isempty(Q_Cortex{k})
        % Creating a block diagonal matrix
        N = endd - start + 1;
        Nout = N / numdipcomp(k) * 3;
        iRow = reshape(repmat(reshape(1:Nout,3,[]), numdipcomp(k), 1), 1, []);
        iCol = reshape(repmat(1:N,3,1), 1, []);
        Q_Cortex{k} = sparse(iRow, iCol, Q_Cortex{k}(:));
        % Applying orientations
        Kernel(start:endd,:) = Q_Cortex{k} * Kernel(start:endd,:);
    end

    start=endd;
end
disp(' ');

%% ===== ASSIGN IMAGING KERNEL =====
% Multiply inverse operator and whitening matrix, so no need to whiten data.
Kernel = Kernel * W;
% Return results structure
Results.ImagingKernel = Kernel;
Results.ImageGridAmp  = [];
Results.Whitener      = W;
if (length(numdipcomp) > 1)
    Results.nComponents = 0;
else
    Results.nComponents = numdipcomp;
end

end

%% ==============================================================================
%  ===== HELPER FUNCTIONS =======================================================
%  ==============================================================================
function W = CalculateWhitener(Modality, C_noise, iChannel, rnkC_noise, isPca)
    [V,D] = eig(C_noise(iChannel,iChannel));
    % In some case, eig returns complex values
    if ~isreal(D)
        error('Complex values in the eigvalues...');
    end
    D = diag(D);
    [D,I] = sort(D,'descend');
    V = V(:,I);
    % No PCA case.
    if ~isPca
        display(['wMNE> Not doing PCA for ' Modality '.'])
        D = 1 ./ D;
        W = diag(sqrt(D)) * V';
    % Rey's approach. MNE has been changed to implement this.
    else
        display(['wMNE> Setting small ' Modality ' eigenvalues to zero.'])
        D = 1 ./ D;
        D(rnkC_noise+1:end) = 0;
        W = diag(sqrt(D)) * V';
        W = W(1:rnkC_noise,:); % This line will reduce the actual number of variables in data and leadfield to the true rank. This was not done in the original MNE C code.
    end
end


text:
function avw_img_write(avw, fileprefix, IMGorient, machine, verbose)

% avw_img_write - write Analyze image files (*.img)
%
% avw_img_write(avw,fileprefix,[IMGorient],[machine],[verbose])
%
% avw.img    - a 3D matrix of image data (double precision).
% avw.hdr    - a struct with image data parameters.  If
%              not empty, this function calls avw_hdr_write.
%
% fileprefix - a string, the filename without the .img
%              extension. If empty, may use avw.fileprefix
%
% IMGorient - optional int, force writing of specified
%             orientation, with values:
%
%   [],  if empty, will use avw.hdr.hist.orient field
%    0,  transverse/axial unflipped (default, radiological)
%    1,  coronal unflipped
%    2,  sagittal unflipped
%    3,  transverse/axial flipped, left to right
%    4,  coronal flipped, anterior to posterior
%    5,  sagittal flipped, superior to inferior
%
% This function will set avw.hdr.hist.orient and write the
% image data in a corresponding order.  This function is
% in alpha development, so it has not been exhaustively
% tested (07/2003). See avw_img_read for more information
% and documentation on the orientation option.
% Orientations 3-5 are NOT recommended!  They are part
% of the Analyze format, but only used in Analyze
% for faster raster graphics during movies.
%
% machine - a string, see machineformat in fread for details.
%           The default here is 'ieee-le'.
%
% verbose - the default is to output processing information to the command
%           window.  If verbose = 0, this will not happen.
%
% Tip: to change the data type, set avw.hdr.dime.datatype to:
%
%     1    Binary             (  1 bit  per voxel)
%     2    Unsigned character (  8 bits per voxel)
%     4    Signed short       ( 16 bits per voxel)
%     8    Signed integer     ( 32 bits per voxel)
%    16    Floating point     ( 32 bits per voxel)
%    32    Complex, 2 floats  ( 64 bits per voxel), not supported
%    64    Double precision   ( 64 bits per voxel)
%   128    Red-Green-Blue     (128 bits per voxel), not supported
%
% See also: avw_write, avw_hdr_write,
%           avw_read, avw_hdr_read, avw_img_read, avw_view
%

% $Revision$ $Date: 2009/01/14 09:24:45 $

% Licence:  GNU GPL, no express or implied warranties
% History:  05/2002, Darren.Weber@flinders.edu.au
%                    The Analyze format is copyright
%                    (c) Copyright, 1986-1995
%                    Biomedical Imaging Resource, Mayo Foundation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%------------------------------------------------------------------------
% Check inputs

if ~exist('avw','var')
    doc avw_img_write;
    ft_error('...no input avw.');
elseif isempty(avw)
    ft_error('...empty input avw.');
elseif ~isfield(avw,'img')
    ft_error('...empty input avw.img');
end

if ~exist('fileprefix','var')
    if isfield(avw,'fileprefix')
        if ~isempty(avw.fileprefix)
            fileprefix = avw.fileprefix;
        else
            fileprefix = [];
        end
    else
        fileprefix = [];
    end
end
if isempty(fileprefix)
    [fileprefix, pathname, filterindex] = uiputfile('*.hdr','Specify an output Analyze .hdr file');
    if pathname, cd(pathname); end
    if ~fileprefix
        doc avw_img_write;
        ft_error('no output .hdr file specified');
    end
end

if contains(fileprefix, '.hdr')
%    fprintf('AVW_IMG_WRITE: Removing .hdr extension from ''%s''\n',fileprefix);
    fileprefix = strrep(fileprefix,'.hdr','');
end
if contains(fileprefix, '.img')
%    fprintf('AVW_IMG_WRITE: Removing .img extension from ''%s''\n',fileprefix);
    fileprefix = strrep(fileprefix,'.img','');
end

if ~exist('IMGorient','var'), IMGorient = ''; end
if ~exist('machine','var'), machine = 'ieee-le'; end
if ~exist('verbose','var'), verbose = 1; end

if isempty(IMGorient), IMGorient = ''; end
if isempty(machine), machine = 'ieee-le'; end
if isempty(verbose), verbose = 1; end

%------------------------------------------------------------------------
% MAIN
if verbose
    version = '????????????????';
    fprintf('\nAVW_IMG_WRITE [v%s]\n',version(12:16));  tic;
end

fid = fopen(sprintf('%s.img',fileprefix),'w',machine);
if fid < 0
    ft_error('Cannot open file %s.img\n',fileprefix);
else
    avw = write_image(fid,avw,fileprefix,IMGorient,machine,verbose);
end

if verbose
    t=toc; fprintf('...done (%5.2f sec).\n\n',t);
end

% MUST write header after the image, to ensure any
% orientation changes during image write are saved
% in the header
avw_hdr_write(avw,fileprefix,machine,verbose);

return

function avw_hdr_write(avw, fileprefix, machine, verbose)

% AVW_HDR_WRITE - Write Analyze header file (*.hdr)
%
% avw_hdr_write(avw,[fileprefix],[machine],[verbose])
%
% eg, avw_hdr_write(avw,'test');
%
% avw        - a struct with .hdr field, which itself is a struct,
%              containing all fields of an Analyze header.
%              For details, see avw_hdr_read.m
%
% fileprefix - a string, the filename without the .hdr extension.
%              If empty, may use avw.fileprefix
%
% machine    - a string, see machineformat in fread for details.
%              The default here is 'ieee-le'.
%
% verbose - the default is to output processing information to the command
%           window.  If verbose = 0, this will not happen.
%
% See also, AVW_HDR_READ AVW_HDR_MAKE
%           AVW_IMG_READ AVW_IMG_WRITE
%

% $Revision$ $Date: 2009/01/14 09:24:45 $

% Licence:  GNU GPL, no express or implied warranties
% History:  05/2002, Darren.Weber@flinders.edu.au
%           02/2003, Bennett.Landman@ieee.org
%                    - more specific data history var sizes
%                    - 02/2003 confirmed, Darren
%
%                    The Analyze format and c code below is copyright
%                    (c) Copyright, 1986-1995
%                    Biomedical Imaging Resource, Mayo Foundation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('verbose','var'), verbose = 1; end
if ~exist('machine','var'), machine = 'ieee-le'; end

if verbose,
    version = '????????????????';
    fprintf('AVW_HDR_WRITE [v%s]\n',version(12:16));  tic;
end

%----------------------------------------------------------------------------
% Check inputs

if ~exist('avw','var')
    ft_warning('...no input avw - calling avw_hdr_make\n');
    avw = avw_hdr_make;
elseif isempty(avw)
    ft_warning('...empty input avw - calling avw_hdr_make\n');
    avw = avw_hdr_make;
elseif ~isfield(avw,'hdr')
    ft_warning('...empty input avw.hdr - calling avw_hdr_make\n');
    avw = avw_hdr_make;
end
if ~isequal(avw.hdr.hk.sizeof_hdr,348)
    ft_error('...avw.hdr.hk.sizeof_hdr must be 348!\n');
end

quit = 0;
if ~exist('fileprefix','var')
    if isfield(avw,'fileprefix')
        if ~isempty(avw.fileprefix)
            fileprefix = avw.fileprefix;
        else
            quit = 1;
        end
    else
        quit = 1;
    end
    if quit
        helpwin avw_hdr_write;
        ft_error('...no input fileprefix - see help avw_hdr_write\n\n');
        return;
    end
end

if contains(fileprefix, '.hdr')
%    fprintf('AVW_HDR_WRITE: Removing .hdr extension from ''%s''\n',fileprefix);
    fileprefix = strrep(fileprefix,'.hdr','');
end

%----------------------------------------------------------------------------
% MAIN

if verbose, tic; end

% % force volume to 4D if necessary; conforms to AVW standard.  i lifted this
% % code from mri_toolbox/avw_hdr_read.m and modified it just a little bit
% % (using minDim = 4; etc)
% minDim = 4;
% currDim = double(avw.hdr.dime.dim(1));
% if ( currDim < minDim )
%    % fprintf( 'Warning %s: Forcing %d dimensions in avw.hdr.dime.dim\n', ...
%    %    mfilename, minDim );
%    avw.hdr.dime.dim(1) = int16(minDim);
%    avw.hdr.dime.dim(currDim+2:minDim+1) = int16(1);
%    avw.hdr.dime.pixdim(1) = int16(minDim);
%    avw.hdr.dime.pixdim(currDim+2:minDim+1) = int16(1);
% end

fid = fopen(sprintf('%s.hdr',fileprefix),'w',machine);
if fid < 0
    ft_error('Cannot write to file %s.hdr\n',fileprefix);
else
    if verbose, fprintf('...writing %s Analyze header.\n',machine); end
    write_header(fid,avw,verbose);
end

if verbose, t=toc; fprintf('...done (%5.2f sec).\n\n',t); end

return

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

function avw = write_image(fid,avw,fileprefix,IMGorient,machine,verbose)

% short int bitpix;    /* Number of bits per pixel; 1, 8, 16, 32, or 64. */
% short int datatype   /* Datatype for this image set */
% /*Acceptable values for datatype are*/
% #define DT_NONE             0
% #define DT_UNKNOWN          0    /*Unknown data type*/
% #define DT_BINARY           1    /*Binary             ( 1 bit per voxel)*/
% #define DT_UNSIGNED_CHAR    2    /*Unsigned character ( 8 bits per voxel)*/
% #define DT_SIGNED_SHORT     4    /*Signed short       (16 bits per voxel)*/
% #define DT_SIGNED_INT       8    /*Signed integer     (32 bits per voxel)*/
% #define DT_FLOAT           16    /*Floating point     (32 bits per voxel)*/
% #define DT_COMPLEX         32    /*Complex,2 floats   (64 bits per voxel)/*
% #define DT_DOUBLE          64    /*Double precision   (64 bits per voxel)*/
% #define DT_RGB            128    /*A Red-Green-Blue datatype*/
% #define DT_ALL            255    /*Undocumented*/

switch double(avw.hdr.dime.datatype)
case   1
    avw.hdr.dime.bitpix = int16( 1); precision = 'bit1';
case   2
    avw.hdr.dime.bitpix = int16( 8); precision = 'uchar';
case   4
    avw.hdr.dime.bitpix = int16(16); precision = 'int16';
case   8
    avw.hdr.dime.bitpix = int16(32); precision = 'int32';
case  16
    avw.hdr.dime.bitpix = int16(32); precision = 'single';
case  32
    ft_error('...complex datatype not yet supported.\n');
case  64
    avw.hdr.dime.bitpix = int16(64); precision = 'double';
case 128
    ft_error('...RGB datatype not yet supported.\n');
otherwise
    ft_warning('...unknown datatype, using type 16 (32 bit floats).\n');
    avw.hdr.dime.datatype = int16(16);
    avw.hdr.dime.bitpix = int16(32); precision = 'single';
end

% write the .img file, depending on the .img orientation
if verbose,
    fprintf('...writing %s precision Analyze image (%s).\n',precision,machine);
end

fseek(fid,0,'bof');

% The standard image orientation is axial unflipped
if isempty(avw.hdr.hist.orient),
    msg = [ '...avw.hdr.hist.orient ~= 0.\n',...
            '   This function assumes the input avw.img is\n',...
            '   in axial unflipped orientation in memory.  This is\n',...
            '   created by the avw_img_read function, which converts\n',...
            '   any input file image to axial unflipped in memory.\n'];
    ft_warning(msg)
end

if isempty(IMGorient),
    if verbose,
        fprintf('...no IMGorient specified, using avw.hdr.hist.orient value.\n');
    end
    IMGorient = double(avw.hdr.hist.orient);
end

if ~isfinite(IMGorient),
    if verbose,
        fprintf('...IMGorient is not finite!\n');
    end
    IMGorient = 99;
end

switch IMGorient,

case 0, % transverse/axial unflipped

    % For the 'transverse unflipped' type, the voxels are stored with
    % Pixels in 'x' axis (varies fastest) - from patient right to left
    % Rows in   'y' axis                  - from patient posterior to anterior
    % Slices in 'z' axis                  - from patient inferior to superior

    if verbose, fprintf('...writing axial unflipped\n'); end

    avw.hdr.hist.orient = uint8(0);

    SliceDim = double(avw.hdr.dime.dim(4)); % z
    RowDim   = double(avw.hdr.dime.dim(3)); % y
    PixelDim = double(avw.hdr.dime.dim(2)); % x
    SliceSz  = double(avw.hdr.dime.pixdim(4));
    RowSz    = double(avw.hdr.dime.pixdim(3));
    PixelSz  = double(avw.hdr.dime.pixdim(2));

    x = 1:PixelDim;
    for z = 1:SliceDim,
        for y = 1:RowDim,
            fwrite(fid,avw.img(x,y,z),precision);
        end
    end

case 1, % coronal unflipped

    % For the 'coronal unflipped' type, the voxels are stored with
    % Pixels in 'x' axis (varies fastest) - from patient right to left
    % Rows in   'z' axis                  - from patient inferior to superior
    % Slices in 'y' axis                  - from patient posterior to anterior

    if verbose, fprintf('...writing coronal unflipped\n'); end

    avw.hdr.hist.orient = uint8(1);

    SliceDim = double(avw.hdr.dime.dim(3)); % y
    RowDim   = double(avw.hdr.dime.dim(4)); % z
    PixelDim = double(avw.hdr.dime.dim(2)); % x
    SliceSz  = double(avw.hdr.dime.pixdim(3));
    RowSz    = double(avw.hdr.dime.pixdim(4));
    PixelSz  = double(avw.hdr.dime.pixdim(2));

    x = 1:PixelDim;
    for y = 1:SliceDim,
        for z = 1:RowDim,
            fwrite(fid,avw.img(x,y,z),precision);
        end
    end

case 2, % sagittal unflipped

    % For the 'sagittal unflipped' type, the voxels are stored with
    % Pixels in 'y' axis (varies fastest) - from patient posterior to anterior
    % Rows in   'z' axis                  - from patient inferior to superior
    % Slices in 'x' axis                  - from patient right to left

    if verbose, fprintf('...writing sagittal unflipped\n'); end

    avw.hdr.hist.orient = uint8(2);

    SliceDim = double(avw.hdr.dime.dim(2)); % x
    RowDim   = double(avw.hdr.dime.dim(4)); % z
    PixelDim = double(avw.hdr.dime.dim(3)); % y
    SliceSz  = double(avw.hdr.dime.pixdim(2));
    RowSz    = double(avw.hdr.dime.pixdim(4));
    PixelSz  = double(avw.hdr.dime.pixdim(3));

    y = 1:PixelDim;
    for x = 1:SliceDim,
        for z = 1:RowDim,
            fwrite(fid,avw.img(x,y,z),precision);
        end
    end

case 3, % transverse/axial flipped

    % For the 'transverse flipped' type, the voxels are stored with
    % Pixels in 'x' axis (varies fastest) - from patient right to left
    % Rows in   'y' axis                  - from patient anterior to posterior*
    % Slices in 'z' axis                  - from patient inferior to superior

    if verbose,
        fprintf('...writing axial flipped (+Y from Anterior to Posterior)\n');
    end

    avw.hdr.hist.orient = uint8(3);

    SliceDim = double(avw.hdr.dime.dim(4)); % z
    RowDim   = double(avw.hdr.dime.dim(3)); % y
    PixelDim = double(avw.hdr.dime.dim(2)); % x
    SliceSz  = double(avw.hdr.dime.pixdim(4));
    RowSz    = double(avw.hdr.dime.pixdim(3));
    PixelSz  = double(avw.hdr.dime.pixdim(2));

    x = 1:PixelDim;
    for z = 1:SliceDim,
        for y = RowDim:-1:1, % flipped in Y
            fwrite(fid,avw.img(x,y,z),precision);
        end
    end

case 4, % coronal flipped

    % For the 'coronal flipped' type, the voxels are stored with
    % Pixels in 'x' axis (varies fastest) - from patient right to left
    % Rows in   'z' axis                  - from patient inferior to superior
    % Slices in 'y' axis                  - from patient anterior to posterior

    if verbose,
        fprintf('...writing coronal flipped (+Z from Superior to Inferior)\n');
    end

    avw.hdr.hist.orient = uint8(4);

    SliceDim = double(avw.hdr.dime.dim(3)); % y
    RowDim   = double(avw.hdr.dime.dim(4)); % z
    PixelDim = double(avw.hdr.dime.dim(2)); % x
    SliceSz  = double(avw.hdr.dime.pixdim(3));
    RowSz    = double(avw.hdr.dime.pixdim(4));
    PixelSz  = double(avw.hdr.dime.pixdim(2));

    x = 1:PixelDim;
    for y = 1:SliceDim,
        for z = RowDim:-1:1,
            fwrite(fid,avw.img(x,y,z),precision);
        end
    end

case 5, % sagittal flipped

    % For the 'sagittal flipped' type, the voxels are stored with
    % Pixels in 'y' axis (varies fastest) - from patient posterior to anterior
    % Rows in   'z' axis                  - from patient superior to inferior
    % Slices in 'x' axis                  - from patient right to left

    if verbose,
        fprintf('...writing sagittal flipped (+Z from Superior to Inferior)\n');
    end

    avw.hdr.hist.orient = uint8(5);

    SliceDim = double(avw.hdr.dime.dim(2)); % x
    RowDim   = double(avw.hdr.dime.dim(4)); % z
    PixelDim = double(avw.hdr.dime.dim(3)); % y
    SliceSz  = double(avw.hdr.dime.pixdim(2));
    RowSz    = double(avw.hdr.dime.pixdim(4));
    PixelSz  = double(avw.hdr.dime.pixdim(3));

    y = 1:PixelDim;
    for x = 1:SliceDim,
        for z = RowDim:-1:1, % superior to inferior
            fwrite(fid,avw.img(x,y,z),precision);
        end
    end

otherwise, % transverse/axial unflipped

    % For the 'transverse unflipped' type, the voxels are stored with
    % Pixels in 'x' axis (varies fastest) - from patient right to left
    % Rows in   'y' axis                  - from patient posterior to anterior
    % Slices in 'z' axis                  - from patient inferior to superior

    if verbose,
        fprintf('...unknown orientation specified, assuming default axial unflipped\n');
    end

    avw.hdr.hist.orient = uint8(0);

    SliceDim = double(avw.hdr.dime.dim(4)); % z
    RowDim   = double(avw.hdr.dime.dim(3)); % y
    PixelDim = double(avw.hdr.dime.dim(2)); % x
    SliceSz  = double(avw.hdr.dime.pixdim(4));
    RowSz    = double(avw.hdr.dime.pixdim(3));
    PixelSz  = double(avw.hdr.dime.pixdim(2));

    x = 1:PixelDim;
    for z = 1:SliceDim,
        for y = 1:RowDim,
            fwrite(fid,avw.img(x,y,z),precision);
        end
    end

end

fclose(fid);

% Update the header
avw.hdr.dime.dim(2:4) = int16([PixelDim,RowDim,SliceDim]);
avw.hdr.dime.pixdim(2:4) = single([PixelSz,RowSz,SliceSz]);

return

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

function write_header(fid,avw,verbose)

    header_key(fid,avw.hdr.hk);
    image_dimension(fid,avw.hdr.dime);
    data_history(fid,avw.hdr.hist);

    % check the file size is 348 bytes
    fbytes = ftell(fid);
    fclose(fid);
    if ~isequal(fbytes,348),
        msg = sprintf('...file size is not 348 bytes!\n');
        ft_warning(msg);
    end

return

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

function header_key(fid,hk)

	% Original header structures - ANALYZE 7.5
	% struct header_key                      /* header key      */
	%       {                                /* off + size      */
	%       int sizeof_hdr                   /*  0 +  4         */
	%       char data_type[10];              /*  4 + 10         */
	%       char db_name[18];                /* 14 + 18         */
	%       int extents;                     /* 32 +  4         */
	%       short int session_error;         /* 36 +  2         */
	%       char regular;                    /* 38 +  1         */
	%       char hkey_un0;                   /* 39 +  1         */
	%       };                               /* total=40 bytes  */

    fseek(fid,0,'bof');

    fwrite(fid, hk.sizeof_hdr(1),   'int32');  % must be 348!

    data_type = sprintf('%-10s',hk.data_type); % ensure it is 10 chars
    fwrite(fid, hk.data_type(1:10), 'uchar');

    db_name   = sprintf('%-18s',hk.db_name);   % ensure it is 18 chars
    fwrite(fid, db_name(1:18),      'uchar');

    fwrite(fid, hk.extents(1),      'int32');
    fwrite(fid, hk.session_error(1),'int16');

    regular   = sprintf('%1s',hk.regular);     % ensure it is 1 char
    fwrite(fid, regular(1),         'uchar');  % might be uint8

    %hkey_un0  = sprintf('%1s',hk.hkey_un0);    % ensure it is 1 char
    %fwrite(fid, hkey_un0(1),        'uchar');
    fwrite(fid, hk.hkey_un0(1),     'uint8');

    %    >Would you set hkey_un0 as char or uint8?
    %   Really doesn't make any difference.  As far as anyone here can remember,
    %   this was just to pad to an even byte boundary for that structure.  I guess
    %   I'd suggest setting it to a uint8 value of 0 (i.e, truly zero-valued) so
    %   that it doesn't look like anything important!
    %   Denny <hanson.dennis2@mayo.edu>

return

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

function image_dimension(fid,dime)

	%struct image_dimension
	%       {                                /* off + size      */
	%       short int dim[8];                /* 0 + 16          */
	%       char vox_units[4];               /* 16 + 4          */
	%       char cal_units[8];               /* 20 + 8          */
	%       short int unused1;               /* 28 + 2          */
	%       short int datatype;              /* 30 + 2          */
	%       short int bitpix;                /* 32 + 2          */
	%       short int dim_un0;               /* 34 + 2          */
	%       float pixdim[8];                 /* 36 + 32         */
	%			/*
	%				pixdim[] specifies the voxel dimensions:
	%				pixdim[1] - voxel width
	%				pixdim[2] - voxel height
	%				pixdim[3] - interslice distance
	%					..etc
	%			*/
	%       float vox_offset;                /* 68 + 4          */
	%       float roi_scale;                 /* 72 + 4          */
	%       float funused1;                  /* 76 + 4          */
	%       float funused2;                  /* 80 + 4          */
	%       float cal_max;                   /* 84 + 4          */
	%       float cal_min;                   /* 88 + 4          */
	%       int compressed;                  /* 92 + 4          */
	%       int verified;                    /* 96 + 4          */
	%       int glmax;                       /* 100 + 4         */
	%       int glmin;                       /* 104 + 4         */
	%       };                               /* total=108 bytes */

	fwrite(fid, dime.dim(1:8),      'int16');
	fwrite(fid, dime.vox_units(1:4),'uchar');
	fwrite(fid, dime.cal_units(1:8),'uchar');
	fwrite(fid, dime.unused1(1),    'int16');
	fwrite(fid, dime.datatype(1),   'int16');
	fwrite(fid, dime.bitpix(1),     'int16');
	fwrite(fid, dime.dim_un0(1),    'int16');
	fwrite(fid, dime.pixdim(1:8),   'float32');
	fwrite(fid, dime.vox_offset(1), 'float32');

    % Ensure compatibility with SPM (according to MRIcro)
    if dime.roi_scale == 0, dime.roi_scale = 0.00392157; end
	fwrite(fid, dime.roi_scale(1),  'float32');

	fwrite(fid, dime.funused1(1),   'float32');
	fwrite(fid, dime.funused2(1),   'float32');
	fwrite(fid, dime.cal_max(1),    'float32');
	fwrite(fid, dime.cal_min(1),    'float32');
	fwrite(fid, dime.compressed(1), 'int32');
	fwrite(fid, dime.verified(1),   'int32');
	fwrite(fid, dime.glmax(1),      'int32');
	fwrite(fid, dime.glmin(1),      'int32');

return

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

function data_history(fid,hist)

	% Original header structures - ANALYZE 7.5
	%struct data_history
	%       {                                /* off + size      */
	%       char descrip[80];                /* 0 + 80          */
	%       char aux_file[24];               /* 80 + 24         */
	%       char orient;                     /* 104 + 1         */
	%       char originator[10];             /* 105 + 10        */
	%       char generated[10];              /* 115 + 10        */
	%       char scannum[10];                /* 125 + 10        */
	%       char patient_id[10];             /* 135 + 10        */
	%       char exp_date[10];               /* 145 + 10        */
	%       char exp_time[10];               /* 155 + 10        */
	%       char hist_un0[3];                /* 165 + 3         */
	%       int views                        /* 168 + 4         */
	%       int vols_added;                  /* 172 + 4         */
	%       int start_field;                 /* 176 + 4         */
	%       int field_skip;                  /* 180 + 4         */
	%       int omax;                        /* 184 + 4         */
	%       int omin;                        /* 188 + 4         */
	%       int smax;                        /* 192 + 4         */
	%       int smin;                        /* 196 + 4         */
	%       };                               /* total=200 bytes */

    descrip     = sprintf('%-80s', hist.descrip);       % 80 chars
    aux_file    = sprintf('%-24s', hist.aux_file);      % 24 chars
    originator  = sprintf('%-10s', hist.originator);    % 10 chars
    generated   = sprintf('%-10s', hist.generated);     % 10 chars
    scannum     = sprintf('%-10s', hist.scannum);       % 10 chars
    patient_id  = sprintf('%-10s', hist.patient_id);    % 10 chars
    exp_date    = sprintf('%-10s', hist.exp_date);      % 10 chars
    exp_time    = sprintf('%-10s', hist.exp_time);      % 10 chars
    hist_un0    = sprintf( '%-3s', hist.hist_un0);      %  3 chars

    % ---
    % The following should not be necessary, but I actually
    % found one instance where it was, so this totally anal
    % retentive approach became necessary, despite the
    % apparently elegant solution above to ensuring that variables
    % are the right length.

    if length(descrip) < 80,
      paddingN = 80-length(descrip);
      padding = char(repmat(double(' '),1,paddingN));
      descrip = [descrip,padding];
    end
    if length(aux_file) < 24,
      paddingN = 24-length(aux_file);
      padding = char(repmat(double(' '),1,paddingN));
      aux_file = [aux_file,padding];
    end
    if length(originator) < 10,
      paddingN = 10-length(originator);
      padding = char(repmat(double(' '),1,paddingN));
      originator = [originator, padding];
    end
    if length(generated) < 10,
      paddingN = 10-length(generated);
      padding = char(repmat(double(' '),1,paddingN));
      generated = [generated, padding];
    end
    if length(scannum) < 10,
      paddingN = 10-length(scannum);
      padding = char(repmat(double(' '),1,paddingN));
      scannum = [scannum, padding];
    end
    if length(patient_id) < 10,
      paddingN = 10-length(patient_id);
      padding = char(repmat(double(' '),1,paddingN));
      patient_id = [patient_id, padding];
    end
    if length(exp_date) < 10,
      paddingN = 10-length(exp_date);
      padding = char(repmat(double(' '),1,paddingN));
      exp_date = [exp_date, padding];
    end
    if length(exp_time) < 10,
      paddingN = 10-length(exp_time);
      padding = char(repmat(double(' '),1,paddingN));
      exp_time = [exp_time, padding];
    end
    if length(hist_un0) < 10,
      paddingN = 10-length(hist_un0);
      padding = char(repmat(double(' '),1,paddingN));
      hist_un0 = [hist_un0, padding];
    end

    % -- if you thought that was anal, try this;
    % -- lets check for unusual ASCII char values!

    if find(double(descrip)>128),
      indexStrangeChar = find(double(descrip)>128);
      descrip(indexStrangeChar) = ' ';
    end
    if find(double(aux_file)>128),
      indexStrangeChar = find(double(aux_file)>128);
      aux_file(indexStrangeChar) = ' ';
    end
    if find(double(originator)>128),
      indexStrangeChar = find(double(originator)>128);
      originator(indexStrangeChar) = ' ';
    end
    if find(double(generated)>128),
      indexStrangeChar = find(double(generated)>128);
      generated(indexStrangeChar) = ' ';
    end
    if find(double(scannum)>128),
      indexStrangeChar = find(double(scannum)>128);
      scannum(indexStrangeChar) = ' ';
    end
    if find(double(patient_id)>128),
      indexStrangeChar = find(double(patient_id)>128);
      patient_id(indexStrangeChar) = ' ';
    end
    if find(double(exp_date)>128),
      indexStrangeChar = find(double(exp_date)>128);
      exp_date(indexStrangeChar) = ' ';
    end
    if find(double(exp_time)>128),
      indexStrangeChar = find(double(exp_time)>128);
      exp_time(indexStrangeChar) = ' ';
    end
    if find(double(hist_un0)>128),
      indexStrangeChar = find(double(hist_un0)>128);
      hist_un0(indexStrangeChar) = ' ';
    end

    % --- finally, we write the fields

    fwrite(fid, descrip(1:80),    'uchar');
    fwrite(fid, aux_file(1:24),   'uchar');

    %orient      = sprintf(  '%1s', hist.orient);        %  1 char
    %fwrite(fid, orient(1),        'uchar');
    fwrite(fid, hist.orient(1),   'uint8');     % see note below on char

    fwrite(fid, originator(1:10), 'uchar');
    fwrite(fid, generated(1:10),  'uchar');
    fwrite(fid, scannum(1:10),    'uchar');
    fwrite(fid, patient_id(1:10), 'uchar');
    fwrite(fid, exp_date(1:10),   'uchar');
    fwrite(fid, exp_time(1:10),   'uchar');
    fwrite(fid, hist_un0(1:3),    'uchar');

    fwrite(fid, hist.views(1),      'int32');
    fwrite(fid, hist.vols_added(1), 'int32');
    fwrite(fid, hist.start_field(1),'int32');
    fwrite(fid, hist.field_skip(1), 'int32');
    fwrite(fid, hist.omax(1),       'int32');
    fwrite(fid, hist.omin(1),       'int32');
    fwrite(fid, hist.smax(1),       'int32');
    fwrite(fid, hist.smin(1),       'int32');

return

% Note on using char:
% The 'char orient' field in the header is intended to
% hold simply an 8-bit unsigned integer value, not the ASCII representation
% of the character for that value.  A single 'char' byte is often used to
% represent an integer value in Analyze if the known value range doesn't
% go beyond 0-255 - saves a byte over a short int, which may not mean
% much in today's computing environments, but given that this format
% has been around since the early 1980's, saving bytes here and there on
% older systems was important!  In this case, 'char' simply provides the
% byte of storage - not an indicator of the format for what is stored in
% this byte.  Generally speaking, anytime a single 'char' is used, it is
% probably meant to hold an 8-bit integer value, whereas if this has
% been dimensioned as an array, then it is intended to hold an ASCII
% character string, even if that was only a single character.
% Denny  <hanson.dennis2@mayo.edu>

% See other notes in avw_hdr_read

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% polint1.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [g,G] = polint1(x,f)
% quadratic polynomial interpolation
% Input:
% x(1:3)  3 pairwise distinct support points
% f(1:3)  corresponding function values
% Output:
% g, G
% the interpolating polynomial is given by
% p(x) = f(1) + g(x - x(1)) + G/2(x - x(1))^2

function [g,G] = polint1(x,f)
f13 = (f(3)-f(1))/(x(3)-x(1));
f12 = (f(2)-f(1))/(x(2)-x(1));
f23 = (f(3)-f(2))/(x(3)-x(2));
g = f13 + f12 - f23;
G = 2*(f13-f12)/(x(3)-x(2));

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
% DeJong_f4.m
% De Jong's f4 function, ND, no noise
%
% described by Clerc in ...
% http://clerc.maurice.free.fr/pso/Semi-continuous_challenge/Semi-continuous_challenge.htm
%
% used to test optimization/global minimization problems
% in Clerc's "Semi-continuous challenge"
%
% f(x) = sum( [1:N].*(in.^4), 2)
%
% x = N element row vector containing [ x0, x1,..., xN ]
%   each row is processed independently,
%   you can feed in matrices of timeXN no prob
%
% example: cost = DeJong_f4([1,2;3,4;5,6])
% note minimum =0 @ x= all zeros

% Brian Birge
% Rev 1.0
% 9/12/04
function [out]=DeJong_f4(in)
 persistent D tlen d

% this speeds routine up a lot, if called from PSO these won't change from
% call to call (repmat is a cpu hog)
 Dx=length(in(1,:));
 tlenx=length(in(:,1));
 if isempty(D) | D~=Dx | tlen~=tlenx
   D=Dx; % dimension of prob
   tlen=tlenx; % how many separate states
   d=repmat([1:D],tlen,1); % needed to vectorize this
 end
 out = sum( d.*(in.^4), 2);

author: HuangCongQing


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
    = grw_run(fid, data, N, sigma, tc, opts)
% This file is the run core for the EG strategy.
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
%           = eg_run(fid, data, tc, opts)
%
% cum_ret: cumulative wealth achived at the end of a period.
% cumprod_ret: cumulative wealth achieved till the end each period.
% daily_ret: daily return achieved by a strategy.
% daily_portfolio: daily portfolio, achieved by the strategy
%
% data: market sequence vectors
% fid: handle for write log file
% gamma: swtiching parameter
% tc: transaction fee rate
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
%          = grw_run(fid, data, 10, 0.00005, 0, opts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n, m] = size(data);

% Variables for return, start with uniform weight
cum_ret = 1;
cumprod_ret = ones(n, 1);
daily_ret = ones(n, 1);
day_weight = ones(m, 1)/m;  %#ok<*NASGU>
day_weight_o = zeros(m, 1);
daily_portfolio = zeros(n, m);

% print file head
fprintf(fid, '-------------------------------------\n');
fprintf(fid, 'Parameters [N:%f, sigma:%f, tc:%f]\n', N, sigma, tc);
fprintf(fid, 'day\t Daily Return\t Total return\n');

fprintf(1, '-------------------------------------\n');
if(~opts.quiet_mode)
    fprintf(1, 'Parameters [N:%f, sigma:%f, tc:%f]\n', N, sigma, tc);
    fprintf(1, 'day\t Daily Return\t Total return\n');
end

% choose normal vectors
xi = normrnd(1/m, sigma^2, N, m);
gamma  = zeros(size(xi));

options = optimset('largescale','off', 'display', 'off');

for t = 1:1:n,
    % Calculate t's portfolio at the beginning of t-th trading day
    [day_weight] = grw_kernel(data(1, :), xi);

    % projection
    C = eye(m, m); d = day_weight;
    A = []; b=[];
    Aeq = ones(m, 1)'; beq = 1;
    lb = zeros(m, 1); ub = ones(m, 1);
    day_weight = lsqlin(C, d, A, b, Aeq, beq, lb, ub, day_weight, options);
    if or((day_weight < 0.0001), (day_weight'*ones(m, 1)>1.00001))
        fprintf(1, 't=%d, sum(day_weight)=%d, pause', t, day_weight'*ones(m, 1));
        pause;
    end

    % Normalize the constraint, always useless
    day_weight = day_weight./sum(day_weight);

    % Cal t's return and total return
    daily_ret(t, 1) = (data(t, :)*day_weight)*(1-tc/2*sum(abs(day_weight-day_weight_o)));
    cum_ret = cum_ret * daily_ret(t, 1);
    cumprod_ret(t, 1) = cum_ret;

    % Adjust weight(t, :) for the transaction cost issue
    day_weight_o = day_weight.*data(t, :)'/daily_ret(t, 1);

    % Monte Carlo update
    k = 1;
    while k <= N,   % Checks the size of our new sample
        u = rand(1); % Choose u from \mu[0, 1]
        j = randi([1, N]);  % Choose j randomly from {1, ... , N}

        if (u <= exp(xi(j, :))./sum(exp(xi(j, :)))*data(t, :)'/(max(data(t, :))) )  % accept-reject condition
            gamma(k, :) = xi(j, :);  k = k+1;   % Accepts xi if the above condition holds
        end
    end

    % Portfolio Update step
    h = normrnd(0, sigma^2, N, m);
    xi = gamma + h;

    % Debug information
    % Time consuming part, other way?
    fprintf(fid, '%d\t%f\t%f\n', t, daily_ret(t, 1), cumprod_ret(t, 1));
    if (~opts.quiet_mode)
        if (~mod(t, opts.display_interval)),
            fprintf(1, '%d\t%f\t%f\n', t, daily_ret(t, 1), cumprod_ret(t, 1));
        end
    end
end

% Debug Information
fprintf(fid, 'GRW(N:%d, sigma:%f, tc:%f), Final return: %f\n',...
    N, sigma, tc, cum_ret);
fprintf(fid, '-------------------------------------\n');
fprintf(1, 'GRW(N:%d, sigma:%f, tc:%f), Final return: %f\n',...
    N, sigma, tc, cum_ret);
fprintf(1, '-------------------------------------\n');

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function [pitch, roll, yaw] = q2att(qnb)
    q11 = qnb(1)*qnb(1); q12 = qnb(1)*qnb(2); q13 = qnb(1)*qnb(3); q14 = qnb(1)*qnb(4);
    q22 = qnb(2)*qnb(2); q23 = qnb(2)*qnb(3); q24 = qnb(2)*qnb(4);
    q33 = qnb(3)*qnb(3); q34 = qnb(3)*qnb(4);
    q44 = qnb(4)*qnb(4);
    C12=2*(q23-q14);
    C22=q11-q22+q33-q44;
    C31=2*(q24-q13); C32=2*(q34+q12); C33=q11-q22-q33+q44;

    pitch = asind(C32);
    roll = atan2d(-C31,C33);
    yaw = atan2d(C12,C22);
    yaw = yaw + (yaw<0)*360;
end

author: yandld


text:
function F = vertcat(F , G)
%VERTCAT   Vertical concatenation of SPHEREFUNV objects.
%   [F ; f] where F is a SPHEREFUNV with two components, and f is a SPHEREFUN or
%   scalar then returns a SPHEREFUNV with three components. The first and second
%   component remain unchanged and the third component is f.
%
%   [f ; F] where F is a SPHEREFUNV with two components, and f is a SPHEREFUN or
%   scalar then returns a SPHEREFUNV with three components. The first is f, and
%   the second and third are the first and second components of F.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

if ( isempty(F) || isempty(G) )
    F = spherefunv;
    return
end

if ( isa(G, 'double') )
    Fc = F.components;
    dom = Fc{1}.domain;
    G = spherefun(G, Fc{1}.domain);
elseif ( isa(F, 'double') )
    Gc = G.components;
    dom = Gc{1}.domain;
    F = spherefun(F, Gc{1}.domain);
elseif ( isa(G, 'spherefun') )
    if ( ~domainCheck(F.components{1}, G) )
        error('SPHEREFUN:SPHEREFUNV:vertcat:domain', 'Inconsistent domains.')
    end
    dom = G.domain;
else
    error('SPHEREFUN:SPHEREFUNV:vertcat:notSupported', ...
        'Vertical concatenation of these objects is not supported.')
end

if ( isa(F, 'spherefunv') )
    op = [ F.components, {G} ];
else
    op = [ {F}, G.components ];
end

F = spherefunv(op, dom);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%sim_lcmrawbasis.m
%Jamie Near, McGill University 2014.
%Modified (greatly simplified) for new FID-A spin system definitions, 2018.
%
% USAGE:
% [RF,out]=sim_lcmrawbasis(n,sw,Bfield,linewidth,metab,tau1,tau2,addref,makeraw,seq)
%
% DESCRIPTION:
% Generate an LCModel .RAW file to be used as an individual metabolite basis
% spectrum in an LCModel basis set.  The relevant characteristics of the
% acquisition can be specified (pulse sequence, number of points, spectral
% width, etc)
%
% INPUTS:
% n         = number of points in fid/spectrum
% sw        = desired spectral width in [Hz]
% Bfield    = main magnetic field strength in [T]
% linewidth = linewidth in [Hz]
% tau1      = first echo time in [ms] (if seq='st' or 'l', tau1 = TE)
% tau2      = second echo time in [ms].  (Used in Press, but not used in SE or LASER.
%             If seq='st', tau2=TM).
% addref    = add reference at 0ppm (for use in LCModel makebasis) ['y' or 'n']
% makeraw   = make output file for lcmodel ['y' or 'n']
% seq       = pulse sequence ['se' for Spin Echo, 'p' for Press, 'st' for Steam, or 'l' for LASER]
% metab     = one of the following choices
%   'H2O'    = Water
%   'Ala'    = Alanine
%   'Asp'    = Aspartate
%   'PCh'    = PhosphoCholine
%   'Cr'     = Creatine
%   'PCr'    = PhosphoCreatine
%   'GABA'   = Gamma-aminobutyric acid (kaiser)
%   'Gln'    = Glutamine
%   'Glu'    = Glutamate
%   'GSH'    = Glutathione
%   'Gly'    = Glycine
%   'Ins'    = Myo-inositol
%   'Lac'    = Lactate
%   'NAA'    = N-acetyl aspartate
%   'Scyllo' = Scyllo-inositol
%   'Tau'    = Taurine
%   'Asc'    = Ascorbate (Vitamin C)
%   'bHB'    = beta-Hydroxybutyrate
%   'bHG'    = beta-Hydroxyglutarate
%   'Glc'    = Glucose
%   'NAAG'   = N-acetyl aspartyl glutamate
%   'GPC'    = Glycero-phosphocholine
%   'PE'     = Phosphoryl ethanolamine
%   'Ser'    = Serine
%
% OUTPUTS:
% RF        = not used.
% out       = Simulated basis spectrum in FID-A structure format.

function [RF,out]=sim_lcmrawbasis(n,sw,Bfield,linewidth,metab,tau1,tau2,addref,makeraw,seq)

load('spinSystems.mat');

eval(['sys=sys' metab ';']);
spins=0;
for k=1:length(sys)
    spins=spins+length(sys(k).shifts);
end
disp(['simulating metabolite ' metab ' with ' num2str(spins) ' spins...  Please Wait...']);
switch seq
    case 'se'
        out = sim_spinecho(n,sw,Bfield,linewidth,sys,tau1);
    case 'p'
        out = sim_press(n,sw,Bfield,linewidth,sys,tau1,tau2);
    case 'st'
        out = sim_steam(n,sw,Bfield,linewidth,sys,tau1,tau2);
    case 'l'
        out = sim_laser(n,sw,Bfield,linewidth,sys,tau1);
    otherwise
        disp(['ERROR:  Sequence ' seq 'not recognized!!!']);
end

if addref=='y'||addref=='Y'
    addref=1;
elseif addref=='n'||addref=='N'
    addref=0;
end

if addref
    sysRef.J=0;
    sysRef.shifts=0;
    sysRef.name='Ref_0ppm';
    sysRef.scaleFactor=1;
    switch seq
        case 'se'
            ref = sim_spinecho(n,sw,Bfield,linewidth,sysRef,tau1);
        case 'p'
            ref = sim_press(n,sw,Bfield,linewidth,sysRef,tau1,tau2);
        case 'st'
            ref = sim_steam(n,sw,Bfield,linewidth,sysRef,tau1,tau2);
        case 'l'
            ref = sim_laser(n,sw,Bfield,linewidth,sysRef,tau1);
        otherwise
            disp(['ERROR:  Sequence ' seq 'not recognized!!!']);
    end
    out=op_addScans(out,ref);
end

if makeraw=='y'||makeraw=='Y'
    makeraw=1;
elseif makeraw=='n'||makeraw=='N'
    makeraw=0;
end

if makeraw
    RF=io_writelcmraw(out,[metab '.RAW'],metab);
else
    RF=[];
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function k = linardKernDiagCompute(kern, x)

% LINARDKERNDIAGCOMPUTE Compute diagonal of LINARD kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for the automatic relevance determination linear kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% computed at the given points.
%
% SEEALSO : linardKernParamInit, kernDiagCompute, kernCreate, linardKernCompute
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

scales = sparse(diag(sqrt(kern.inputScales)));
x = x*scales;

k = sum(x.*x, 2)*kern.variance;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% DISPLAY - display an EEG data class underlying structure
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Nov. 2008

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function display(obj);

    tmpMMO = memmapfile(obj.dataFile, 'writable', obj.writable, 'format', { 'single' obj.dimensions 'x' });
    if obj.transposed, disp('Warning: data does not display properly for memory mapped file which have been transposed'); end
    disp(tmpMMO.data.x);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [output,timing] = global_solve_upper(p,p_original,x,options,uppersolver,timing)

if ~p.solver.uppersolver.constraint.binary
    if ~isempty(p.binary_variables)
        local_gave_good = find(abs(x(p_original.binary_variables)-fix(x(p_original.binary_variables)))< options.bnb.inttol);
        p.lb(p.binary_variables(local_gave_good)) = fix(x(p.binary_variables(local_gave_good)));
        p.ub(p.binary_variables(local_gave_good)) = fix(x(p.binary_variables(local_gave_good)));
        for i = 1:3
            p = propagate_bounds_from_evaluations(p);
            p = propagate_bounds_from_equalities(p);
            p = update_monomial_bounds(p);
        end
    end

    if ~isempty(p_original.integer_variables)
        local_gave_good = find(abs(x(p.integer_variables)-fix(x(p.integer_variables)))< options.bnb.inttol);
        p.lb((p.integer_variables(local_gave_good))) = fix(x(p.integer_variables(local_gave_good)));
        p.ub((p.integer_variables(local_gave_good))) = fix(x(p.integer_variables(local_gave_good)));
        for i = 1:3
            p = propagate_bounds_from_evaluations(p);
            p = propagate_bounds_from_equalities(p);
            p = update_monomial_bounds(p);
        end
    end
end

% The bounds and relaxed solutions have been computed w.r.t to the relaxed
% bilinear model. We only need the original bounds and variables.
p.lb = p.lb(1:length(p_original.c));
p.ub = p.ub(1:length(p_original.c));

x = x(1:length(p_original.c));
%
p_upper = p_original;

p_upper = restoreQ(p_upper);

p_upper.lb = p.lb;
p_upper.ub = p.ub;

% KNITRO does not like fixed binary/integer variables, so we remove them
fixed_binary = find(p_upper.lb(p_upper.binary_variables) == p_upper.ub(p_upper.binary_variables));
fixed_integer = find(p_upper.lb(p_upper.integer_variables) == p_upper.ub(p_upper.integer_variables));
p_upper.binary_variables(fixed_binary) = [];if length(p_upper.binary_variables)==0;p_upper.binary_variables = [];end
p_upper.integer_variables(fixed_integer) = [];if length(p_upper.integer_variables)==0;p_upper.integer_variables = [];end

% Pick an initial point (this can be a bit tricky...)
lbinfbounds = find(isinf(p.lb));
ubinfbounds = find(isinf(p.ub));
switch p.options.bmibnb.localstart
    case 'relaxed'
        p_upper.x0 = x;
    case 'boxcenter'
        p_upper.x0 = x;
        p_upper.x0(lbinfbounds)=0;
        p_upper.x0(ubinfbounds)=0;
    case 'zero'
        p_upper.x0 = zeros(length(p.lb),1);
    otherwise
        error('Unknown local starting point option');
end

% Find x(a)*x(b) where x(a) or x(b) is fixed, and remove from Q
[a,b,kk] = find(triu(2*p_upper.Q));
for i = 1:length(a)
    if p.lb(a(i)) == p.ub(a(i))
        p_upper.c(b(i)) = p_upper.c(b(i)) + p.lb(a(i))*kk(i);
        p_upper.Q(a(i),b(i))=0;
        p_upper.Q(b(i),a(i))=0;
    elseif p.lb(b(i)) == p.ub(b(i))
        p_upper.c(a(i)) = p_upper.c(a(i)) + p.lb(b(i))*kk(i);
        p_upper.Q(a(i),b(i))=0;
        p_upper.Q(b(i),a(i))=0;
    end
end

% Save time
p_upper.options.saveduals = 0;

% Solve upper bounding problem
p_upper.options.warmstart = 1;
tstart = tic;
try
    output = feval(uppersolver,p_upper);
catch
    output.Primal = zeros(length(p_upper.lb),1);
    output.problem = -1;
end
if isempty(output.Primal)
     output.Primal = zeros(length(p_upper.lb),1);
end

timing.uppersolve = timing.uppersolve + toc(tstart);

% Project into the box (numerical issue)
output.Primal(output.Primal<p_upper.lb) = p_upper.lb(output.Primal<p_upper.lb);
output.Primal(output.Primal>p_upper.ub) = p_upper.ub(output.Primal>p_upper.ub);

function p = restoreQ(p)
% Sometimes the Q-term has been lifted and relaxed, but it might be the
% case that we are using a pure Q solver, and then we should really solve
% it as such
if nnz(p.Q)==0 && any(p.variabletype > 0)
    if p.solver.uppersolver.objective.quadratic.convex

        nonlinear = find(p.variabletype > 0 & p.variabletype <= 2);
        linear = find(p.variabletype == 0);
        % Find quadratic and linear terms
        used_in_c = find(p.c);
        quadraticterms = used_in_c(find(ismember(used_in_c,nonlinear)));
        Q = zeros(length(p.c),length(p.c));
        if ~isempty(quadraticterms)
            usedinquadratic = zeros(1,length(p.c));
            for i = 1:length(quadraticterms)
                Qij = p.c(quadraticterms(i));
                power_index = find(p.monomtable(quadraticterms(i),:));
                if length(power_index) == 1
                    Q(power_index,power_index) = Qij;
                else
                    Q(power_index(1),power_index(2)) = Qij/2;
                    Q(power_index(2),power_index(1)) = Qij/2;
                end
            end
        end
        p.Q = Q;
        p.c(nonlinear) = 0;
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function marginal = marginal_nodes(engine, b, nodes, t, add_ev, is_fam)
% function marginal = marginal_nodes(engine, b, nodes, t, add_ev, is_fam) (jtree_2TBN)

if nargin < 6, is_fam = 0; end
ss = engine.slice_size;

if ~is_fam & (t > 1) & all(nodes<=ss)
  nodes = nodes + ss;
end

if t==1
  c = clq_containing_nodes(engine.jtree_engine1, nodes, is_fam);
else
  c = clq_containing_nodes(engine.jtree_engine, nodes, is_fam);
end
if c == -1
  error(['no clique contains ' nodes])
end
bigpot = b.clpot{c};
pot = marginalize_pot(bigpot, nodes, engine.maximize);
marginal = pot_to_marginal(pot);

% we convert the domain to the unrolled numbering system
% so that add_ev_to_dmarginal (maybe called in update_ess) extracts the right evidence.
if t > 1
  marginal.domain = nodes+(t-2)*engine.slice_size;
end
assert(~add_ev);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [ panoLineImg ] = viewRegionHypothesis( img, AllCuboid, thres )
%VIEWREGIONHYPOTHESIS Summary of this function goes here
%   Detailed explanation goes here

AllLines = zeros(0,6);
for cid = 1:AllCuboid.count
    if AllCuboid.score(cid)<thres
        continue;
    end

    view = AllCuboid.views(cid,:);
    xyzBox = AllCuboid.xyzBox(cid,:);
    numRectangle = sum(view~=0);

    for rid = 1:numRectangle
        rect = reshape(xyzBox((rid*12-11):(rid*12)), 3, 4)';
        lines = lineFromTwoPoint(rect([1 2 3 4], :), rect([2 3 4 1], :));
        AllLines = [AllLines; lines];
    end
end

% img = imresize( im2double(img), [512 1024]);
% panoLineImg = paintParameterLine( AllLines, 1024, 512, img);
img = imresize( im2double(img), [1024 2048]);
panoLineImg = paintParameterLine( AllLines, 2048, 1024, img);

% imshow(panoLineImg);

end

author: zouchuhang

repo: LayoutNet


text:
% MP_SP_DRCHLT_L2_PROJ_UDOTN: assign the normal degrees of freedom trough an L2 projection for a multipatch geometry.
%  To be used with the 'RT' and 'NDL' spaces. The imposed condition reads   u \cdot n = h \cdot n
%
%   [vel, normal_dofs] = mp_sp_drchlt_l2_proj_udotn (space, msh, gnum, ornt, bnd_sides, bnd_func)
%
% INPUTS:
%
%    space:     space object (see sp_vector_div_transform)
%    msh:       mesh object (see msh_cartesian)
%    gnum:      global numbering of the degrees of freedom (see mp_interface_hdiv)
%    ornt:      global orientation of the degrees of freedom (see mp_interface_hdiv)
%  boundaries: array of structures containing the information for the boundaries (see mp_geo_load)
%    bnd_sides: boundary sides on which the Dirichlet condition is imposed
%    bnd_func:  the condition to be imposed (h in the equation)
%
% OUTPUT:
%
%     vel:         assigned value to the normal degrees of freedom
%     normal_dofs: global numbering of the normal basis functions
%
% Copyright (C) 2014, 2015 Rafael Vazquez
%
%    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 3 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, see <http://www.gnu.org/licenses/>.

function [u, normal_dofs] = mp_sp_drchlt_l2_proj_udotn (space, msh, gnum, ornt, boundaries, refs, bnd_func)

  if (isa (space, 'sp_multipatch'))
    warning ('For spaces of the class SP_MULTIPATCH, using the function SP_DRCHLT_PROJ_L2_UDOTN inside the class')
    [u, dofs] = sp_drchlt_proj_l2_udotn (sp, msh, refs, bnd_func);
    return
  end

  ndof = max ([gnum{:}]);
  M = spalloc (ndof, ndof, 3*ndof);
  rhs = zeros (ndof, 1);

  normal_dofs = [];
  for iref = refs
    for bnd_side = 1:boundaries(iref).nsides
      iptc = boundaries(iref).patches(bnd_side);
      iside = boundaries(iref).faces(bnd_side);

      if (strcmpi (space{iptc}.transform, 'div-preserving'))
        msh_side = msh_eval_boundary_side (msh{iptc}, iside);
        sp_side = sp_eval_boundary_side (space{iptc}, msh_side);
      else
        error ('The function only works with div-conforming spaces')
      end

      normal_dofs = union (normal_dofs, gnum{iptc}(sp_side.dofs));
      for idim = 1:msh{iptc}.rdim
        x{idim} = reshape (msh_side.geo_map(idim,:,:), msh_side.nqn, msh_side.nel);
      end
      g = bnd_func (x{:}, iref);

      M_loc = op_u_v (sp_side, sp_side, msh_side, ones(size(x{1})));
      rhs_loc = op_fdotn_v (sp_side, msh_side, g);

      global_dofs = gnum{iptc}(sp_side.dofs);
      M(global_dofs, global_dofs) = M(global_dofs, global_dofs) + M_loc;
      rhs(global_dofs) = rhs(global_dofs) + rhs_loc;
    end
  end

  u = M(normal_dofs, normal_dofs) \ rhs(normal_dofs);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function brs = ssr_brs_point_source(X,phi,xs,irs,conf)
%SSR_BRS_POINT_SOURCE binaural room scanning (BRS) set for use with the
%SoundScape Renderer
%
%   Usage: brs = ssr_brs_point_source(X,phi,xs,irs,conf)
%
%   Input parameters:
%       X       - listener position / m
%       phi     - azimuthal head orientation / rad
%       xs      - source position / m
%       irs     - impulse response data set for the second sources
%       conf    - configuration struct (see SFS_config)
%
%   Output parameters:
%       brs     - conf.N x 2*nangles matrix containing all impulse responses (2
%                 channels) for every angle of the BRS set
%
%   SSR_BRS_POINT_SOURCE(X,phi,xs,irs,conf) prepares a BRS set for a reference
%   source (single point source) for the given listener position. One way to use
%   this BRS set is using the SoundScapeRenderer (SSR), see
%   http://spatialaudio.net/ssr/
%
%   See also: ir_generic, ir_point_source, get_ir

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 5;
nargmax = 5;
narginchk(nargmin,nargmax);
isargposition(X);
isargxs(xs);
isargscalar(phi);
isargstruct(conf);

%% ===== Computation =====================================================
brs = ssr_brs(X,phi,[xs 0 -1 0 1],1,irs,conf);

author: sfstoolbox

repo: sfs-matlab


text:
function [vfN, vfR, vrN, vrR, vnetN, vnetR, uL, uC, u0L, u0C, lnxL, lnxC] = projectOntoSubspace(A, vf, vr, vnet, u, u0, lnx,printLevel,rowBool,colBool)
% Projects flux, net flux, potential and logarithmic concentration onto
% their respective subspaces of A using projection matrices generated either
% derived from SVD, or by using the Moore-Penrose pseudoinverse
%
% Optionally, a subset of the matrix A may be chosen by using A(rowBool,colBool)
% but then only the true rows of u, u0, lnx, and true columns of vf,vr,vnet are
% projected and the remaining rows and columns are not affected
%
% Let `M` denote the Moore-Penrose pseudoinverse of A and the subscripts are the following
% `_R` row space,
% `_N` nullspace,
% `_C` column space,
% `_L` left nullspace,
%
% Example for flux of net flux
%
% Let
%
% .. math::
%      vf   &= vf_R + vf_N \\
%      vf_R &= M A vf = PR vf \\
%      vf_N &= (I - M A) vf = PN vf
%
% Example for potential or logarithmic concentration
%
% Let
%
% .. math::
%      u   &= u_C + u_L \\
%      u_C &= A M u = PC u \\
%      u_L &= (I - A M) u = PL u
%
% USAGE:
%
%    [vfN, vfR, vrN, vrR, vnetN, vnetR, uL, uC, u0L, u0C, lnxL, lnxC] = projectOntoSubspace(modelT, vf, vr, vnet, u, u0, lnx)
%
% INPUTS:
%    A          `m x n` matrix
%    vf:        `n x 1` - forward flux
%    vr:        `n x 1` - reverse flux
%    vnet:      `n x 1` - net flux
%    u:         `m x 1` - chemical potential
%    u0:        `m x 1` - standard chemical potential
%    lnx:       `m x 1` - logarithmic concentration
% OPTIONAL INPUTS
%    rowBool    `m x 1` - boolean indicating the subset of rows of A
%    colBool    'n x 1' - boolean indicating the subset of cols of A
%
% OUTPUTS:
%    vfN:       forward flux - nullspace
%    vfR:       forward flux - row space
%    vrN:       reverse flux - nullspace
%    vrR:       reverse flux - row space
%    vnetN:     net flux - nullspace
%    vnetR:     net flux - row space
%    uL:        chemical potential - left nullspace
%    uC:        chemical potential - column space
%    lnxL:      logarithmic concentration - left nullspace
%    lnxC:      logarithmic concentration - column space

if ~isempty(vf)
    if any((vnet(colBool)- vf +vr)>1e-12) %sanity check
    error('Net flux does not equal the difference between forward and reverse flux')
    end
end

if ~isempty(lnx)
    if any((u - u0 - lnx)>1e-12)
        error('Chemical potential does not equal standard chemical potential plus logarithmic conc.')
    end
end

if ~exist('printLevel','var')
    printLevel=0;
end

[m,n]=size(A);

if ~exist('rowBool','var')
    rowBool=true(m,1);
end

if ~exist('colBool','var')
    colBool=true(n,1);
end

%generate fake outputs, or populate with unprojected vectors, part of which
%will be overwritten with the projected vectos further below.
if ~exist('u','var')
    uL=NaN*ones(m,1);
    uC=NaN*ones(m,1);
else
    uL=u;
    uC=u;
end
if ~exist('u0','var')
    u0L=NaN*ones(m,1);
    u0C=NaN*ones(m,1);
else
    u0L=u0;
    u0C=u0;
end
if ~exist('lnx','var')
    lnxL=NaN*ones(m,1);
    lnxC=NaN*ones(m,1);
else
    lnxL=lnx;
    lnxC=lnx;
end
if ~exist('vf','var')
    vfN=NaN*ones(n,1);
    vfR=NaN*ones(n,1);
else
    vfN=vf;
    vfR=vf;
end
if ~exist('vr','var')
    vrN=NaN*ones(n,1);
    vrR=NaN*ones(n,1);
else
    vrN=vr;
    vrR=vr;
end
if ~exist('vnet','var')
    vnetN=NaN*ones(n,1);
    vnetR=NaN*ones(n,1);
else
    vnetN=vnet;
    vnetR=vnet;
end

%generate projection matrices
sub_space='all';

[PR,PN,PC,PL]=subspaceProjector(A(rowBool,colBool),printLevel,sub_space);

%potential
if ~isempty(u)
    %concentration
    uC(rowBool)=PC*u(rowBool);
    uL(rowBool)=PL*u(rowBool);
end

if ~isempty(u0)
    %standard potential
    u0L(rowBool)=PL*u0(rowBool);
    u0C(rowBool)=PC*u0(rowBool);
end

if ~isempty(lnx)
    %concentration
    lnxC(rowBool)=PC*lnx(rowBool);
    lnxL(rowBool)=PL*lnx(rowBool);
end

%flux
if ~isempty(vnet)
    vnetN(colBool)=PN*vnet(colBool);
    vnetR(colBool)=PR*vnet(colBool);
end
if ~isempty(vf)
    vfN(colBool)=PN*vf(colBool);
    vrR(colBool)=PR*vr(colBool);
end
if ~isempty(vr)
    vrN(colBool)=PN*vr(colBool);
    vrR(colBool)=PR*vr(colBool);
end

author: opencobra

repo: cobratoolbox


text:
function [ fileList, fileFormat, fileFilter ] = java_getfile( DialogType, WindowTitle, DefaultDir, SelectionMode, FilesOrDir, Filters, defaultFilter)
% JAVA_GETFILE: Java-based file selection for opening and saving.
%
% USAGE: [fileList, fileFormat, fileFilter] = java_getfile(DialogType,WindowTitle,DefaultDir,SelectionMode, FilesOrDir,Filters,defaultFilter)
%
% INPUT :
%    - DialogType    : {'open', 'save'}
%    - WindowTitle   : String
%    - DefaultDir    : To ignore, set to []
%    - SelectionMode : {'single', 'multiple'}
%    - FilesOrDir    : {'files', 'dirs', 'files_and_dirs'}
%    - Filters       : {NbFilters x 2} cell array
%                      Filters(i,:) = {{'.ext1', '.ext2', '_tag1'...}, Description}
%    - defaultFilter : can be 1) the index of the default file filter
%                             2) the name of the filer to be used
% OUTPUT:
%    - fileList      : Cell-array of strings, full paths to the files that were selected
%    - fileFormat    : String that represents the format of the files that were selected
%    - fileFilter    : File filter that was selected when selecting the files

% @=============================================================================
% 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, 2008-2014

import org.brainstorm.file.*;

global GlobalData FileSelectorStatus;

%% ===== CONFIGURE DIALOG =====
% Initialize returned variables
fileList   = {};
fileFormat = '';
fileFilter = [];

% DialogType
if strcmpi(DialogType, 'save')
    DialogType = BstFileSelector.TYPE_SAVE;
else
    DialogType = BstFileSelector.TYPE_OPEN;
end
% SelectionMode
if strcmpi(SelectionMode, 'multiple')
    SelectionMode = BstFileSelector.SELECTION_MULTIPLE;
else
    SelectionMode = BstFileSelector.SELECTION_SINGLE;
end
% Files and/or directories
switch lower(FilesOrDir)
    case 'dirs'
        FilesOrDir = BstFileSelector.DIRECTORIES_ONLY;
    case 'files_and_dirs'
        FilesOrDir = BstFileSelector.FILES_AND_DIRECTORIES;
    otherwise
        FilesOrDir = BstFileSelector.FILES_ONLY;
end
% Default folder: user folder
if isempty(DefaultDir)
    DefaultDir = bst_get('UserDir');
end

% Filters
iDefaultFileFilter = 1;
for i=1:size(Filters, 1)
    % Filters cell array has the following format:
    %   - One row per filter,
    %   - A filter row can be {{'extensions_list'}, 'Description', 'FormatName'}
    %                      or {{'extensions_list'}, 'Description'}
    if (size(Filters, 2) == 3)
        fileFilters(i) = java_create('org.brainstorm.file.BstFileFilter', '[Ljava.lang.String;Ljava.lang.String;Ljava.lang.String;', Filters{i,1}, Filters{i,2}, Filters{i,3});
    else
        fileFilters(i) = java_create('org.brainstorm.file.BstFileFilter', '[Ljava.lang.String;Ljava.lang.String;', Filters{i,1}, Filters{i,2});
    end
    % If it is default file filter
    if ~isempty(defaultFilter) && (isnumeric(defaultFilter) && (i == defaultFilter)) || ...
            (ischar(defaultFilter) && (size(Filters, 2) == 3) && strcmpi(Filters{i,3}, defaultFilter))
        iDefaultFileFilter = i;
    end
end

% If a progress bar is displayed : hide it while displaying the file selector
pBar = GlobalData.Program.ProgressBar;
if ~isempty(pBar) && isfield(pBar, 'jWindow') && java_call(pBar.jWindow, 'isVisible')
    pBarHidden = 1;
    bst_progress('hide');
else
    pBarHidden = 0;
end

%% ===== HIDE MODAL WINDOWS =====
% Get brainstorm frame
jBstFrame = bst_get('BstFrame');
% If the frame is defined
if ~isempty(jBstFrame)
    jDialogModal = [];
    jDialogAlwaysOnTop = [];
    for i=1:length(GlobalData.Program.GUI.panels)
        panelContainer = get(GlobalData.Program.GUI.panels(i), 'container');
        panelContainer = panelContainer.handle{1};
        if isa(panelContainer, 'javax.swing.JDialog') && java_call(panelContainer, 'isModal')
            % A modal JDialog is found => Set it non non-modal
            jDialogModal = panelContainer;
            java_call(jDialogModal, 'setModal', 'Z', 0);
        end
        if (isa(panelContainer, 'javax.swing.JDialog') || isa(panelContainer, 'javax.swing.JFrame')) && java_call(panelContainer, 'isAlwaysOnTop')
            % An AlwaysOnTop frame is found => Remove always on top attribute
            jDialogAlwaysOnTop = panelContainer;
            java_call(jDialogAlwaysOnTop, 'setAlwaysOnTop', 'Z', 0);
        end
    end
end

%% ===== CREATE SELECTION DIALOG =====
% Initialize dialog status
FileSelectorStatus = 0;
% Create object
jBstSelector = java_create('org.brainstorm.file.BstFileSelector', 'ILjava.lang.String;Ljava.lang.String;II[Lorg.brainstorm.file.BstFileFilter;I', ...
                           DialogType, WindowTitle, DefaultDir, SelectionMode, FilesOrDir, fileFilters, iDefaultFileFilter - 1);
% Initialize a mutex (a figure that will be closed in the BstFileSelector close callback)
bst_mutex('create', 'FileSelector');
% Get JFileChooser dialog
jFileChooser = java_call(jBstSelector, 'getJFileChooser');
% Set dialog callback
java_setcb(jFileChooser, 'ActionPerformedCallback', @FileSelectorAction, ...
                         'PropertyChangeCallback',  @FileSelectorPropertyChanged);
drawnow;
% Display file selector
java_call(jBstSelector, 'showSameThread');
% MAGIC: Print something to the console output to get it to flush something, if not it crashes on MacOS 10.9.2
drawnow;
fprintf(1, ' \b');
% Wait for the file selector to be closed by the user
if ~isempty(bst_mutex('get', 'FileSelector'))
    bst_mutex('waitfor', 'FileSelector');
end

%% ===== PROCESS SELECTED FILES =====
% If user clicked OK after having selected a valid file
if FileSelectorStatus
    % Get file filter => file format
    fileFilter = java_call(jFileChooser, 'getFileFilter');
    fileFormat = char(java_call(fileFilter, 'getFormatName'));

    % If multiple selection
    if (SelectionMode == BstFileSelector.SELECTION_MULTIPLE)
        % Get selected files
        fs = java_call(jFileChooser, 'getSelectedFiles');
        % Convert them to a cell array of filenames
        fileList = cell(length(fs), 1);
        for i=1:length(fs)
            fileList{i} = char(java_call(fs(i), 'getAbsolutePath'));
        end
    % Else: single selection
    else
        % Get selected file
        fileList = char(java_call(jFileChooser, 'getSelectedFile'));

        % If SAVE dialog
        if (DialogType == BstFileSelector.TYPE_SAVE)
            % Get required extension
            suffix = fileFilter.getSuffixes();
            suffix = char(suffix(1));
            % Replace current extension with required extension (ONLY IF SUFFIX IS EXTENSION)
            if (suffix(1) == '.') && ~isequal(suffix, '.folder')
                [selPath, selBase, selExt] = bst_fileparts(fileList);
                fileList = bst_fullfile(selPath, [selBase, suffix]);
            end

            % If file already exist
            if file_exist(fileList) && ~isequal(suffix, '.folder')
                if ~java_dialog('confirm', sprintf('File already exist.\nDo you want to overwrite it?'), 'Save file')
                    fileList = [];
                    fileFormat = [];
                    fileFilter = [];
                end
            end
        end
    end
else
    fileList = [];
end

% Restore modal panels
if ~isempty(jBstFrame)
    if ~isempty(jDialogModal)
        java_call(jDialogModal, 'setModal', 'Z', 1);
    end
    if ~isempty(jDialogAlwaysOnTop)
        java_call(jDialogAlwaysOnTop, 'setAlwaysOnTop', 'Z', 1);
    end
end
% Restore progress bar
if pBarHidden
    bst_progress('show');
end

%% ===== CALLBACK FUNCTION =====
    function FileSelectorAction(h, ev)
        switch (char(java_call(ev, 'getActionCommand')))
            case 'ApproveSelection'
                FileSelectorStatus = 1;
            otherwise
                FileSelectorStatus = 0;
        end
        % Release mutex
        bst_mutex('release', 'FileSelector');
    end

    function FileSelectorPropertyChanged(h, ev)
        import org.brainstorm.file.*;
        % Only when saving
        if (DialogType == BstFileSelector.TYPE_SAVE)
            switch char(java_call(ev, 'getPropertyName'))
                case 'fileFilterChanged'
                    % Get new filter
                    newFilter = java_call(ev, 'getNewValue');
                    % New suffix
                    newSuffix = java_call(newFilter, 'getSuffixes');
                    newSuffix = char(newSuffix(1));
                    if isequal(newSuffix, '.folder')
                        newSuffix = '';
                    end
                    % Get old filename
                    selFile = java_call(jFileChooser, 'getSelectedFile');
                    if isempty(selFile)
                        oldFilename = DefaultDir;
                    else
                        jSelFile = java_call(jFileChooser, 'getSelectedFile');
                        oldFilename = char(java_call(jSelFile, 'getAbsolutePath'));
                    end

                    % Replace old extension by new one
                    [fPath, fBase, fExt] = bst_fileparts(oldFilename);
                    % Brainstorm or external file
                    if ~isempty(newSuffix) && (newSuffix(1) == '_')
                        fBase = strrep(fBase, ['_' newSuffix(2:end)], '');
                        fBase = strrep(fBase, [newSuffix(2:end) '_'], '');
                        fBase = [newSuffix(2:end) '_' fBase];
                        fExt  = '.mat';
                    else
                        fExt = newSuffix;
                    end
                    newFilename = bst_fullfile(fPath, [fBase, fExt]);
                    % Update default filename
                    java_call(jFileChooser, 'setSelectedFile', 'Ljava.io.File;', java_create('java.io.File', 'Ljava.lang.String;', newFilename));

                case 'directoryChanged'
                    DefaultDir = strrep(DefaultDir, char(java_call(ev, 'getOldValue')), char(java_call(ev, 'getNewValue')));
            end
        end
    end
end

author: brainstorm-tools


text:
function CaptureFigVid(ViewZ, FileName,OptionZ)
% CaptureFigVid(ViewZ, FileName,OptionZ)
% Captures a video of the 3D plot in the current axis as it rotates based
% on ViewZ and saves it as 'FileName.mpg'. Option can be specified.
%
% ViewZ:     N-rows with 2 columns, each row are the view angles in
%            degrees, First column is azimuth (pan), Second is elevation
%            (tilt) values outside of 0-360 wrap without error,
%            *If a duration is specified, angles are used as nodes and
%            views are equally spaced between them (other interpolation
%            could be implemented, if someone feels so ambitious).
%            *If only an initial and final view is given, and no duration,
%            then the default is 100 frames.
% FileName:  Name of the file of the produced animation. Because I wrote
%            the program, I get to pick my default of mpg-4, and the file
%            extension .mpg will be appended, even if the filename includes
%            another file extension. File is saved in the working
%            directory.
% (OptionZ): Optional input to specify parameters. The ones I use are given
%            below. Feel free to add your own. Any or all fields can be
%            used
% OptionZ.FrameRate: Specify the frame rate of the final video (e.g. 30;)
% OptionZ.Duration: Specify the length of video in seconds (overrides
%    spacing of view angles) (e.g. 3.5;)
% OptionZ.Periodic: Logical to indicate if the video should be periodic.
%    Using this removed the final view so that when the video repeats the
%    initial and final view are not the same. Avoids having to find the
%    interval between view angles. (e.g. true;)
%
% % % % Example (shown in published results, video attached) % % % %
% figure(171);clf;
% surf(peaks,'EdgeColor','none','FaceColor','interp','FaceLighting','phong')
% daspect([1,1,.3]);axis tight;
% OptionZ.FrameRate=15;OptionZ.Duration=5.5;OptionZ.Periodic=true;
% CaptureFigVid([-20,10;-110,10;-190,80;-290,10;-380,10],'WellMadeVid',OptionZ)
%
% Known issues: MPEG-4 video option only available on Windows machines. See
% fix where the VideoWriter is called.
%
% Getframe is used to capture image and current figure must be on monitor 1
% if multiple displays are used. Does not work if no display is used.
%
% Active windows that overlay the figure are captured in the movie.  Set up
% the current figure prior to calling the function. If you don't specify
% properties, such as tick marks and aspect ratios, they will likely change
% with the rotation for an undesirable effect.

% Cheers, Dr. Alan Jennings, Research assistant professor,
% Department of Aeronautics and Astronautics, Air Force Institute of Technology

%% preliminaries

% initialize optional argument
if nargin<3;     OptionZ=struct([]); end

% check orientation of ViewZ, should be two columns and >=2 rows
if size(ViewZ,2)>size(ViewZ,1); ViewZ=ViewZ.'; end
if size(ViewZ,2)>2
warning('AJennings:VidWrite',...
    'Views should have n rows and only 2 columns. Deleting extraneous input.');
ViewZ=ViewZ(:,1:2); %remove any extra columns
end

% Create video object
daObj=VideoWriter(FileName,'MPEG-4'); %my preferred format
% daObj=VideoWriter(FileName); %for default video format.
% MPEG-4 CANNOT BE USED ON UNIX MACHINES
% set values:
% Frame rate
if isfield(OptionZ,'FrameRate')
    daObj.FrameRate=OptionZ.FrameRate;
end
% Durration (if frame rate not set, based on default)
if isfield(OptionZ,'Duration') %space out view angles
    temp_n=round(OptionZ.Duration*daObj.FrameRate); % number frames
    temp_p=(temp_n-1)/(size(ViewZ,1)-1); % length of each interval
    ViewZ_new=zeros(temp_n,2);
    % space view angles, if needed
    for inis=1:(size(ViewZ,1)-1)
        ViewZ_new(round(temp_p*(inis-1)+1):round(temp_p*inis+1),:)=...
            [linspace(ViewZ(inis,1),ViewZ(inis+1,1),...
             round(temp_p*inis)-round(temp_p*(inis-1))+1).',...
             linspace(ViewZ(inis,2),ViewZ(inis+1,2),...
             round(temp_p*inis)-round(temp_p*(inis-1))+1).'];
    end
    ViewZ=ViewZ_new;
end
% space view angles, if needed
if length(ViewZ)==2 % only initial and final given
    ViewZ=[linspace(ViewZ(1,1),ViewZ(end,1)).',...
           linspace(ViewZ(1,2),ViewZ(end,2)).'];
end
% Periodicity
if isfield(OptionZ,'Periodic')&&OptionZ.Periodic==true
ViewZ=ViewZ(1:(end-1),:); %remove last sample
end
% open object, preparatory to making the video
open(daObj);

%% rotate the axis and capture the video
for kathy=1:size(ViewZ,1)
    view(ViewZ(kathy,:)); drawnow;
    writeVideo(daObj,getframe(gcf)); %use figure, since axis changes size based on view
end

%% clean up
close(daObj);

author: Sable

repo: mcbench-benchmarks