text:
% SP_BSPLINE_FLUID: Construct different pair of B-Splines spaces on the physical domain for fluid problems.
%
%   [spv, spp] = sp_bspline_fluid (elem_name, knots, nsub, ...
%                                               degreep, regularity, msh)
%
% INPUTS:
%
%   elem_name:  the name of the element. Right now 'TH' (Taylor-Hood),
%                'NDL' (Nedelec, 2nd family), 'RT' (Raviart-Thomas) and
%                'SG' (SubGrid) are supported.
%   knots:      knot vector of the coarse geometry.
%   nsub:       number of subdivisions of each interval.
%   degreep:    degree of the pressure space along each parametric direction
%   regularity: continuity of the pressure space along each parametric direction
%   msh:        msh object containing (in the field msh.qn) the points
%                along each parametric direction in the parametric
%                domain at which to evaluate, i.e. quadrature points
%                or points for visualization (see msh_cartesian).
%
% OUTPUT:
%
%   spv: object representing the discrete velocity function space (see sp_vector)
%   spp: object representing the discrete pressure function space (see sp_scalar)
%
%   For more details, see:
%      A. Buffa, C. de Falco, G. Sangalli,
%      IsoGeometric Analysis: Stable elements for the 2D Stokes equation
%      Internat. J. Numer. Methods Fluids, 2011
%
%      A. Bressan, G. Sangalli,
%      Isogeometric discretizations of the Stokes problem: stability
%       analysis by the macroelement technique
%      IMA J. Numer. Anal., 2013.
%
% Copyright (C) 2009, 2010, 2011 Carlo de Falco
% Copyright (C) 2011, 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 [spv, spp] = sp_bspline_fluid (element_name, ...
                   knots, nsub_p, degree_p, regularity_p, msh)

% Construction of the knot vectors and the discrete space for the velocity
switch (lower (element_name))
  case {'th'}
    knotsp = kntrefine (knots, nsub_p-1, degree_p, regularity_p);
    spp = sp_bspline (knotsp, degree_p, msh);

    degree_v = degree_p + 1;
    regularity_v = regularity_p;
    nsub_v = nsub_p;
    knots_v = kntrefine (knots, nsub_v-1, degree_v, regularity_v);
    scalar_space = sp_bspline (knots_v, degree_v, msh);
    for idim = 1:msh.ndim
      scalar_spaces{idim} = scalar_space;
    end
    spv = sp_vector (scalar_spaces, msh);
    clear scalar_spaces scalar_space

  case {'sg'}
    knotsp = kntrefine (knots, nsub_p-1, degree_p, regularity_p);
    spp = sp_bspline (knotsp, degree_p, msh);

    degree_v = degree_p + 1;
    regularity_v = regularity_p+1;
    nsub_v = 2*nsub_p;
    knots_v = kntrefine (knots, nsub_v-1, degree_v, regularity_v);
    scalar_space = sp_bspline (knots_v, degree_v, msh);
    for idim = 1:msh.ndim
      scalar_spaces{idim} = scalar_space;
    end
    spv = sp_vector (scalar_spaces, msh);
    clear scalar_spaces scalar_space

  case {'ndl'}
% In this case the regularity is assigned first in the velocity space
    degree_h1 = degree_p + 1;
    regularity_h1 = regularity_p + 1;
    knots_h1 = kntrefine (knots, nsub_p-1, degree_h1, regularity_h1);
    [knots_hdiv, degree_hdiv] = knt_derham (knots_h1, degree_h1, 'Hdiv');

    degree_v = degree_h1;
    for idim = 1:msh.ndim
      knots_v{idim} = knots_hdiv{idim}{idim};
      for jdim = setdiff (1:msh.ndim, idim)
        knots_v{jdim} = sort ([knots_hdiv{idim}{jdim}, unique(knots_hdiv{idim}{jdim})]);
      end
      scalar_spaces{idim} = sp_bspline (knots_v, degree_v, msh);
    end
    spv = sp_vector (scalar_spaces, msh, 'div-preserving');
    clear scalar_spaces

    [knotsp, degp] = knt_derham (knots_h1, degree_h1, 'L2');
    spp = sp_bspline (knotsp, degp, msh, 'integral-preserving');

  case {'rt'}
% In this case the regularity is assigned first in the velocity space
    degree_h1 = degree_p + 1;
    regularity_h1 = regularity_p + 1;
    knots_h1 = kntrefine (knots, nsub_p-1, degree_h1, regularity_h1);

    [knots_v, degree_v] = knt_derham (knots_h1, degree_h1, 'Hdiv');
    for idim = 1:msh.ndim
      scalar_spaces{idim} = sp_bspline (knots_v{idim}, degree_v{idim}, msh);
    end
    spv = sp_vector (scalar_spaces, msh, 'div-preserving');
    clear scalar_spaces

    [knotsp, degp] = knt_derham (knots_h1, degree_h1, 'L2');
    spp = sp_bspline (knotsp, degp, msh, 'integral-preserving');

%     if (nargout == 3)
%       PI = b2nst__ (spp, knotsp, degree_p, msh);
%     end
  otherwise
    error ('sp_bspline_fluid: unknown element type')
end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [out_val] = per_graph(in_vec)
% per_graph.m - draws PER graph
global numPkts PER_snr PER_exp PER_bitPayload;

% PER graph settings (actual values are user-selected (GUI))
%%%%% numPkts = 200;
%%%%% PER_snr = 29:33;

% Input parameters
Clk        = in_vec(1);
bitsPerBlk = in_vec(2:9);
bmode      = in_vec(10);
dbits      = in_vec(11:end);
txbits     = dbits(1:end/2);
rxbits     = dbits(end/2+1:end);

% Viterbi trace-back depth (hard-coded for now)
link_delay = 128; %%34;

% Reset PER graph data
if (Clk==0)  PER_exp = zeros(size(PER_snr));  end

% Count number of error packets
if (floor(Clk/numPkts)<length(PER_snr))
    idx = floor(Clk/numPkts)+1;
    % Compare tx, rx bits
    txb = txbits( 1 : PER_bitPayload-link_delay );
    rxb = rxbits( link_delay+1 : PER_bitPayload );

    if (1) % for PER
        if ~all(txb==rxb)  PER_exp(idx) = PER_exp(idx) + 1;  end
    else   % for BER
        PER_exp(idx) = PER_exp(idx) + sum(xor(txb,rxb))/(PER_bitPayload-link_delay);
    end

    % Display PER information
    disp( ['Pkt:' num2str(mod(Clk,numPkts)+1) ', Idx=' num2str(idx) ...
           ', PER_list: ' num2str(PER_exp)] );

    % save to file (every tenth packet)
    if (mod(Clk,10)==9)
        fid = fopen('..\per_results.txt','w');
        fprintf(fid, 'Pkt: %d, Idx=%d,', mod(Clk,numPkts)+1, idx);
        fprintf(fid, 'SNR_list: [%s]\n', num2str(PER_snr));
        fprintf(fid, 'Pkt: %d, Idx=%d,', mod(Clk,numPkts)+1, idx);
        fprintf(fid, 'PER_list: [%s]\n', num2str(PER_exp));
        fclose(fid);
    end
end

if (Clk==length(PER_snr)*numPkts)
    PER_fig = figure;
    % Draw PER Curves
    PER_exp = PER_exp/numPkts;
    semilogy(PER_snr, PER_exp, 'g');
    %axis([0 25 1/1000 1]);
    axis([0 50 1/1000 1]);

    % Label graph...
    title('PER Curve');
    xlabel('SNR (dB)');
    ylabel('PER');
    legend('Simulation Result');
    grid on;

    % Display end message
    disp('PER test complete.');
end

% Output percentage done
out_val = min([100 100*Clk/(length(PER_snr)*numPkts)]);

author: Sable

repo: mcbench-benchmarks


text:
function dist_signed = quad_point_dist_signed_2d ( q, p )

%*****************************************************************************80
%
%% QUAD_POINT_DIST_SIGNED_2D: signed distanct ( quadrilateral, point ) in 2D.
%
%  Discussion:
%
%    The quadrilateral must be convex.  DIST_SIGNED is actually the maximum
%    of the signed distances from the point to each of the four lines that
%    make up the quadrilateral.
%
%    Essentially, if the point is outside the convex quadrilateral,
%    only one of the signed distances can be positive, or two can
%    be positive and equal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real Q(2,4), the vertices of the quadrilateral.
%
%    Input, real P(2,1), the point which is to be checked.
%
%    Output, real DIST_SIGNED, the signed distance from the
%    point to the convex quadrilateral.  If DIST_SIGNED is
%    0.0, the point is on the boundary;
%    negative, the point is in the interior;
%    positive, the point is in the exterior.
%

%
%  Compare the signed distance from each line segment to the point,
%  with the signed distance to the midpoint of the opposite line.
%
%  The signed distances should all be negative if the point is inside.
%
%  Side 12
%
  dis12 = line_exp_point_dist_signed_2d ( q(1:2,1), q(1:2,2), p );

  pm(1:2,1) = 0.5 * ( q(1:2,3) + q(1:2,4) );

  dis = line_exp_point_dist_signed_2d ( q(1:2,1), q(1:2,2), pm );

  if ( 0.0 < dis )
    dis = -dis;
    dis12 = -dis12;
  end
%
%  Side 23
%
  dis23 = line_exp_point_dist_signed_2d ( q(1:2,2), q(1:2,3), p );

  pm(1:2,1) = 0.5 * ( q(1:2,4) + q(1:2,1) );

  dis = line_exp_point_dist_signed_2d ( q(1:2,2), q(1:2,3), pm );

  if ( 0.0 < dis )
    dis = -dis;
    dis23 = -dis23;
  end
%
%  Side 34
%
  dis34 = line_exp_point_dist_signed_2d ( q(1:2,3), q(1:2,4), p );

  pm(1:2,1) = 0.5 * ( q(1:2,1) + q(1:2,2) );

  dis = line_exp_point_dist_signed_2d ( q(1:2,3), q(1:2,4), pm );

  if ( 0.0 < dis )
    dis = -dis;
    dis34 = -dis34;
  end
%
%  Side 41
%
  dis41 = line_exp_point_dist_signed_2d ( q(1:2,4), q(1:2,1), p );

  pm(1:2,1) = 0.5 * ( q(1:2,2) + q(1:2,3) );

  dis = line_exp_point_dist_signed_2d ( q(1:2,4), q(1:2,1), pm );

  if ( 0.0 < dis )
    dis = -dis;
    dis41 = -dis41;
  end

  dist_signed = max ( dis12, max ( dis23, max ( dis34, dis41 ) ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = scaleNoiseOut(noise, mu, varSigma)

% SCALENOISEOUT A simple noise model that scales and centres the data.

% NOISE

y = zeros(size(mu));
for i = 1:size(mu, 2)
  y(:, i) = noise.bias(i) + mu(:, i)*noise.scale(i);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function varargout = process_inverse_2016( varargin )
% PROCESS_INVERSE_2016: Compute an inverse model.

% @=============================================================================
% 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, 2012-2019

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % ===== PROCESS =====
    % Description the process
    sProcess.Comment     = 'Compute sources [2016]';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Sources';
    sProcess.Index       = 326;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/SourceEstimation';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'raw'};
    sProcess.OutputTypes = {'results', 'results'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Option: Output
    sProcess.options.output.Comment = {'Kernel only: shared', 'Kernel only: one per file', 'Full results: one per file'};
    sProcess.options.output.Type    = 'radio';
    sProcess.options.output.Value   = 1;
    % Options: MNE options
    sProcess.options.inverse.Comment = {'panel_inverse_2016', 'Source estimation options: '};
    sProcess.options.inverse.Type    = 'editpref';
    sProcess.options.inverse.Value   = bst_inverse_linear_2016();
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    OutputFiles = {};

    % ===== GET OPTIONS =====
    % Default inverse options
    OPTIONS = Compute();
    % Get options edited by the user
    OPTIONS = struct_copy_fields(OPTIONS, sProcess.options.inverse.Value, 1);
    % Output
    switch (sProcess.options.output.Value)
        % Kernel only: shared
        case 1
            [sChannels, iStudies] = bst_get('ChannelForStudy', unique([sInputs.iStudy]));
            iStudies = unique(iStudies);
            iDatas   = [];
            OPTIONS.ComputeKernel = 1;
        % Kernel only: one per file
        case 2
            iStudies = [sInputs.iStudy];
            iDatas   = [sInputs.iItem];
            OPTIONS.ComputeKernel = 1;
        % Full results: one per file
        case 3
            iStudies = [sInputs.iStudy];
            iDatas   = [sInputs.iItem];
            OPTIONS.ComputeKernel = 0;
    end
    % No messages
    OPTIONS.DisplayMessages = 0;

    % ===== START COMPUTATION =====
    % Call head modeler
    [AllFiles, errMessage] = Compute(iStudies, iDatas, OPTIONS);
    % Report errors
    if isempty(AllFiles) && ~isempty(errMessage)
        bst_report('Error', sProcess, sInputs, errMessage);
        return;
    elseif ~isempty(errMessage)
        bst_report('Warning', sProcess, sInputs, errMessage);
    end
    % For shared kernels: Return only the source files corresponding to the recordings that were in input
    if isempty(iDatas)
        % Loop on the output files (all links): find the ones that match the input
        for iFile = 1:length(AllFiles)
            % Resolve link: get data file
            [ResFile, DataFile] = file_resolve_link(AllFiles{iFile});
            % Find one that matches the inputs
            iInput = find(file_compare({sInputs.FileName}, file_short(DataFile)));
            % If founf: add to the output files
            if ~isempty(iInput)
                OutputFiles{end+1} = AllFiles{iFile};
            end
        end
    else
        OutputFiles = AllFiles;
    end
end

%% ===== COMPUTE INVERSE SOLUTION =====
% USAGE:      OPTIONS = Compute()
%         OutputFiles = Compute(iStudies, iDatas, OPTIONS)
%
% Authors: Sylvain Baillet, October 2002
%          Esen Kucukaltun-Yildirim, 2004
%          Syed Ashrafulla, John Mosher, Rey Ramirez, 2009-2012
%          Francois Tadel, John Mosher, 2009-2016
%
function [OutputFiles, errMessage] = Compute(iStudies, iDatas, OPTIONS)
    % Initialize returned variables
    OutputFiles = {};
    errMessage = [];
    % Default options settings
    Def_OPTIONS = struct(...
        'InverseMethod',       'minnorm', ... % A string that specifies the imaging method
        'InverseMeasure',      'dspm', ...
        'SourceOrient',        'fixed', ...
        'DataTypes',           [], ...     % Cell array of strings: list of modality to use for the reconstruction (MEG, MEG GRAD, MEG MAG, EEG)
        'Comment',             '', ...     % Inverse solution description (optional)
        'DisplayMessages',     1, ...
        'ComputeKernel',       1);         % If 1, compute MN kernel to be applied subsequently to data instead of full ImageGridAmp array
    % Return default options
    if (nargin < 2)
        OutputFiles = Def_OPTIONS;
        return;
    end
    % Use default options
    if (nargin < 3) || isempty(OPTIONS)
        OPTIONS = Def_OPTIONS;
    else
        % Check field names of passed OPTIONS and fill missing ones with default values
        OPTIONS = struct_copy_fields(OPTIONS, Def_OPTIONS, 0);
    end

    %% ===== GET INPUT INFORMATION =====
    isShared = isempty(iDatas);
    % Get all the study structures
    sStudies = bst_get('Study', unique(iStudies));
    % Get channel studies
    if isShared
        sChanStudies = sStudies;
    else
        [tmp, iChanStudies] = bst_get('ChannelForStudy', unique(iStudies));
        sChanStudies = bst_get('Study', iChanStudies);
    end
    % Check that there are channel files available
    if any(cellfun(@isempty, {sChanStudies.Channel}))
        errMessage = 'No channel file available.';
        return;
    end
    % Check head model
    if any(cellfun(@isempty, {sChanStudies.HeadModel}))
        errMessage = 'No head model available.';
        return;
    end
    % Check noise covariance
    if any(cellfun(@isempty, {sChanStudies.NoiseCov}))
        errMessage = 'No noise covariance matrix available.';
        return;
    end
    % Loop through all the channel files to find the available modalities and head model types
    AllMod = {};
    HeadModelType = 'surface';
    MEGMethod = [];
    nSamplesNoise = [];
    nSamplesData  = [];
    for i = 1:length(sChanStudies)
        AllMod = union(AllMod, sChanStudies(i).Channel.DisplayableSensorTypes);
        if isempty(sChanStudies(i).HeadModel(sChanStudies(i).iHeadModel).MEGMethod)
            AllMod = setdiff(AllMod, {'MEG GRAD','MEG MAG','MEG'});
        end
        if isempty(sChanStudies(i).HeadModel(sChanStudies(i).iHeadModel).EEGMethod)
            AllMod = setdiff(AllMod, {'EEG'});
        end
        if isempty(sChanStudies(i).HeadModel(sChanStudies(i).iHeadModel).ECOGMethod)
            AllMod = setdiff(AllMod, {'ECOG'});
        end
        if isempty(sChanStudies(i).HeadModel(sChanStudies(i).iHeadModel).SEEGMethod)
            AllMod = setdiff(AllMod, {'SEEG'});
        end
        if ~strcmpi(sChanStudies(i).HeadModel(sChanStudies(i).iHeadModel).HeadModelType, 'surface')
            HeadModelType = sChanStudies(i).HeadModel(sChanStudies(i).iHeadModel).HeadModelType;
        end
        if ~isempty(sChanStudies(i).HeadModel(sChanStudies(i).iHeadModel).MEGMethod) && isempty(MEGMethod)
            MEGMethod = sChanStudies(i).HeadModel(sChanStudies(i).iHeadModel).MEGMethod;
        end
        % First file only: Load the number of samples from the covariance files
        if (i == 1)
            % Noise covariance
            if (length(sChanStudies(i).NoiseCov) >= 1) && ~isempty(sChanStudies(i).NoiseCov(1).FileName)
                covMat = load(file_fullpath(sChanStudies(i).NoiseCov(1).FileName), 'nSamples');
                if isfield(covMat, 'nSamples') && ~isempty(covMat.nSamples)
                    nSamplesNoise = covMat.nSamples;
                end
            end
            % Data covariance
            if (length(sChanStudies(i).NoiseCov) >= 2) && ~isempty(sChanStudies(i).NoiseCov(2).FileName)
                covMat = load(file_fullpath(sChanStudies(i).NoiseCov(2).FileName), 'nSamples');
                if isfield(covMat, 'nSamples') && ~isempty(covMat.nSamples)
                    nSamplesData = covMat.nSamples;
                end
            end
        end
    end
    % Keep only MEG and EEG
    if any(ismember(AllMod, {'MEG GRAD','MEG MAG'}))
        AllMod = intersect(AllMod, {'MEG GRAD', 'MEG MAG', 'EEG', 'ECOG', 'SEEG'});
    else
        AllMod = intersect(AllMod, {'MEG', 'EEG', 'ECOG', 'SEEG'});
    end
    % Check that at least one modality is available
    if isempty(AllMod)
        errMessage = 'No valid sensor types to estimate sources: please calculate an appropriate headmodel.';
        return;
    end

    %% ===== SELECT INVERSE METHOD =====
    % Select method
    if OPTIONS.DisplayMessages
        % Options dialog window
        sMethod = gui_show_dialog('Compute sources', @panel_inverse_2016, 1, [], AllMod, isShared, HeadModelType, nSamplesNoise, nSamplesData);
        if isempty(sMethod)
            return;
        end
        % Override default options
        OPTIONS = struct_copy_fields(OPTIONS, sMethod, 1);

        % === BRAINENTROPY MEM ===
        % Display additional option windows
        if strcmpi(OPTIONS.InverseMethod, 'mem')
            % No data files found
            if isShared
                errMessage = 'Cannot compute shared kernels with this method.';
                return
            end
            % Install/load brainentropy plugin
            [isInstalled, errMessage] = bst_plugin('Install', 'brainentropy', 1);
            if ~isInstalled
                return;
            end
            % Default options
            MethodOptions = be_main();
            % Interface to edit options
            MethodOptions = gui_show_dialog('MEM options', @panel_brainentropy, [], [], MethodOptions);
            % Add fields that are not defined by the options of the MEM interface
            if ~isempty(MethodOptions)
                switch (HeadModelType)
                    case {'surface', 'ImageGrid'}
                        MethodOptions.SourceOrient{1} = 'fixed';
                    case 'volume'
                        MethodOptions.SourceOrient{1} = 'free';
                        MethodOptions.flagSourceOrient = [0 0 2 0];
                end
            end
            % Canceled by user
            if isempty(MethodOptions)
                return
            end
            % Add options to list
            OPTIONS = struct_copy_fields(OPTIONS, MethodOptions, 1);
        end
    end
    % If no MEG and no EEG selected
    if isempty(OPTIONS.DataTypes)
        errMessage = 'Please select at least one modality.';
        return;
    end
    % Tags corresponding to the different methods
    methodTag = panel_inverse_2016('GetMethodComment', OPTIONS.InverseMethod, OPTIONS.InverseMeasure);

    %% ===== COMMENT =====
    % Base comment: "METHOD: MODALITIES"
    if isempty(OPTIONS.Comment)
        OPTIONS.Comment = [methodTag, ': ' GetModalityComment(OPTIONS.DataTypes)];
    end
    % Add source orientation option string
    strOptions = '';
    if isempty(OPTIONS.SourceOrient)
        strOptions = 'Mixed';
    elseif ~strcmpi(OPTIONS.InverseMethod, 'mem')
        switch (OPTIONS.SourceOrient{1})
            case 'fixed',      strOptions = 'Constr';
            case 'loose',      strOptions = 'Loose';
            case 'free',       strOptions = 'Unconstr';
        end
    end
    % Add Kernel/Full option string
    if ~OPTIONS.ComputeKernel
        if ~isempty(strOptions)
            strOptions = [',' strOptions];
        end
        strOptions = ['Full', strOptions];
    end
    % Final comment
    if ~isempty(strOptions)
        OPTIONS.Comment = [OPTIONS.Comment, '(', strOptions, ')'];
    end

    %% ===== LOOP ON INPUT FILES =====
    % Initializations
    initOPTIONS = OPTIONS;
    % Display progress bar
    bst_progress('start', 'Compute sources', 'Initialize...', 0, 3*length(iStudies) + 1);
    % Process each input
    for iEntry = 1:length(iStudies)
        OPTIONS = initOPTIONS;

        % ===== LOAD CHANNEL FILE =====
        bst_progress('text', 'Reading channel information...');
        % Get study structure
        iStudy = iStudies(iEntry);
        sStudy = bst_get('Study', iStudy);
        % Check if default study
        isDefaultStudy = strcmpi(sStudy.Name, bst_get('DirDefaultStudy'));
        % Get channel file for study
        [sChannel, iStudyChannel] = bst_get('ChannelForStudy', iStudy);
        ChannelFile = sChannel.FileName;
        % Load channel file
        ChannelMat = in_bst_channel(ChannelFile, 'Channel', 'Projector');

        % ===== LOAD DATA FILES =====
        bst_progress('text', 'Getting bad channels...');
        % Single inverse file
        if ~isShared
            % Get only one file
            DataFile = sStudy.Data(iDatas(iEntry)).FileName;
            % Load data file info (only 'mem' requires the recordings to be loaded here)
            if strcmpi(OPTIONS.InverseMethod, 'mem')
                DataMat = in_bst_data(DataFile, 'ChannelFlag', 'Time', 'nAvg', 'Leff', 'F');
            else
                DataMat = in_bst_data(DataFile, 'ChannelFlag', 'Time', 'nAvg', 'Leff');
            end
            ChannelFlag = DataMat.ChannelFlag;
            nAvg        = DataMat.nAvg;
            Leff        = DataMat.Leff;
            Time        = DataMat.Time;
            % Is it a Raw file?
            isRaw = strcmpi(sStudy.Data(iDatas(iEntry)).DataType, 'raw');
        % Shared inverse kernel
        else
            % Get all the dependent data files
            [iRelatedStudies, iRelatedData] = bst_get('DataForStudy', iStudy);
            % List all the data files
            nAvgAll     = zeros(1,length(iRelatedStudies));
            LeffAll     = zeros(1,length(iRelatedStudies));
            BadChannels = [];
            nChannels   = [];
            for i = 1:length(iRelatedStudies)
                % Get data file
                sStudyRel = bst_get('Study', iRelatedStudies(i));
                DataFull = file_fullpath(sStudyRel.Data(iRelatedData(i)).FileName);
                % Read bad channels and nAvg
                DataMat = load(DataFull, 'ChannelFlag', 'nAvg', 'Leff');
                if isfield(DataMat, 'nAvg') && ~isempty(DataMat.nAvg)
                    nAvgAll(i) = DataMat.nAvg;
                else
                    nAvgAll(i) = 1;
                end
                if isfield(DataMat, 'Leff') && ~isempty(DataMat.Leff)
                    LeffAll(i) = DataMat.Leff;
                else
                    LeffAll(i) = 1;
                end
                % Count number of times the channe is bad
                if isempty(BadChannels)
                    BadChannels = double(DataMat.ChannelFlag < 0);
                else
                    BadChannels = BadChannels + (DataMat.ChannelFlag < 0);
                end
                % Channel number
                if isempty(nChannels)
                    nChannels = length(DataMat.ChannelFlag);
                elseif (nChannels ~= length(DataMat.ChannelFlag))
                    errMessage = 'All data files must have the same number of channels.';
                    continue;
                end
            end
            % Get list of sensors selected for inversion
            iChanInv = good_channel(ChannelMat.Channel, [], OPTIONS.DataTypes);
            % Mark all the channels that are not selected here as good
            BadChannels(setdiff(1:length(BadChannels), iChanInv)) = 0;

            % === CHECK nAVG ===
            % if ~isempty(iRelatedStudies) && any(nAvgAll ~= nAvgAll(1)) && isFirstWarnAvg
            %     % Display a warning in a dialog window
            %     if OPTIONS.DisplayMessages
            %         isConfirm = java_dialog('confirm', ...
            %             ['Warning: You should estimate separataley the sources for the averages and the single trials.', 10 ...
            %             'The level of noise in the files might be different, this may cause inaccurate results.' 10 10 ...
            %             ' - For several averages: compute sources separately for each file.' 10 ...
            %             ' - For single trials: compute a shared solution, and move the averages in another condition.' 10 10 ...
            %             'Ignore this warning and compute sources ?'], 'Compute sources');
            %         if ~isConfirm
            %             return
            %         end
            %     % Return a warning message
            %     else
            %         errMessage = [errMessage 'Mixing averages and single trials. Result might be inaccurate.' 10];
            %     end
            %     isFirstWarnAvg = 0;
            % end
            nAvg = min([nAvgAll 1]);
            Leff = min([LeffAll 1]);

            % === BAD CHANNELS ===
            if any(BadChannels)
                % Display a warning in a dialog window
                if OPTIONS.DisplayMessages
                    % Build list of bad channels
                    strBad = '';
                    iBad = find(BadChannels);
                    for i = 1:length(iBad)
                        strBad = [strBad sprintf('%d(%d)   ', iBad(i), BadChannels(iBad(i)))];
                        if (mod(i,6) == 0)
                            strBad = [strBad 10];
                        end
                    end
                    % Ask user confirmation
                    [res, isCancel] = java_dialog('input', ...
                        ['Some channels are bad in at least one file (total ' num2str(length(iRelatedStudies)) ' files): ' 10 ...
                         '(in parentheses, the number of files for which the channel is bad)' 10 10 ...
                         strBad 10 10 ...
                         'The following channels will be considered as BAD for' 10 ...
                         'all the recordings and excluded from the source estimation:' 10 10], ...
                        'Exclude bad channels', [], sprintf('%d ', iBad));
                    if isCancel
                        continue;
                    elseif (~isempty(res) && isempty(str2num(res)))
                        errMessage = 'Invalid bad channel list of indices.';
                        continue;
                    end
                    % Get bad channels
                    BadChannels = 0 * BadChannels;
                    if ~isempty(res)
                        iBad = str2num(res);
                        BadChannels(iBad) = 1;
                    end
                % Return a warning message
                else
                    errMessage = [errMessage 'Bad channels for all the trials: ' sprintf('%d ', find(BadChannels)) 10];
                end
            end
            % Build a resulting ChannelFlag
            ChannelFlag = ones(length(ChannelMat.Channel), 1);
            ChannelFlag(BadChannels > 0) = -1;
            % No data loaded
            DataFile = [];
            Time = [];
        end

        % ===== CHANNEL FLAG =====
        % Get the list of good channels
        GoodChannel = good_channel(ChannelMat.Channel, ChannelFlag, OPTIONS.DataTypes);
        if isempty(GoodChannel)
            errMessage = [errMessage 'No good channels available.' 10];
            break;
        end

        % ===== LOAD NOISE COVARIANCE =====
        % Get channel study
        sStudyChannel = bst_get('Study', iStudyChannel);
        % Load NoiseCov file
        NoiseCovMat = load(file_fullpath(sStudyChannel.NoiseCov(1).FileName));
        % Check for NaN values in the noise covariance
        if ~isempty(NoiseCovMat.NoiseCov) && (nnz(isnan(NoiseCovMat.NoiseCov(GoodChannel, GoodChannel))) > 0)
            errMessage = [errMessage 'The noise covariance contains NaN values. Please re-calculate it after tagging correctly the bad channels in the recordings.' 10];
            break;
        end
        % Divide noise covariance by number of trials
        if ~isempty(nAvg) && (nAvg > 1)
            NoiseCovMat.NoiseCov = NoiseCovMat.NoiseCov ./ nAvg;
        end

        % ===== LOAD DATA COVARIANCE =====
        % Load DataCov file
        if (length(sStudyChannel.NoiseCov) >= 2) && ~isempty(sStudyChannel.NoiseCov(2).FileName)
            DataCovMat = load(file_fullpath(sStudyChannel.NoiseCov(2).FileName));
            % Check for NaN values in the noise covariance
            if ~isempty(DataCovMat.NoiseCov) && (nnz(isnan(DataCovMat.NoiseCov(GoodChannel, GoodChannel))) > 0)
                errMessage = [errMessage 'The data covariance contains NaN values. Please re-calculate it after tagging correctly the bad channels in the recordings.' 10];
                break;
            end
            % Divide data covariance by number of trials
            if isempty(nAvg) && (nAvg > 1)
                DataCovMat.NoiseCov = DataCovMat.NoiseCov ./ nAvg;
            end
        else
            DataCovMat = [];
        end
        % Beamformers: Require a data covariance matrix
        if strcmpi(OPTIONS.InverseMethod, 'lcmv') && isempty(DataCovMat)
            errMessage = [errMessage 'You need to calculate a data covariance before using the "beamformer" option.' 10];
            break;
        end
        % Shrinkage: Require the FourthMoment matrix
        if strcmpi(OPTIONS.NoiseMethod, 'shrink') && ...
                ((~isempty(DataCovMat)  && (~isfield(DataCovMat, 'FourthMoment')  || isempty(DataCovMat.FourthMoment))) || ...
                 (~isempty(NoiseCovMat) && (~isfield(NoiseCovMat, 'FourthMoment') || isempty(NoiseCovMat.FourthMoment))))
            errMessage = [errMessage 'Please recalculate the noise and data covariance matrices for using the "automatic shrinkage" option.' 10];
            break;
        end

        % ===== LOAD HEAD MODEL =====
        bst_progress('text', 'Loading head model...');
        bst_progress('inc', 1);
        % Get headmodel file
        HeadModelFile = sStudyChannel.HeadModel(sStudyChannel.iHeadModel).FileName;
        % Load head model
        HeadModel = in_bst_headmodel(HeadModelFile, 0, 'Gain', 'GridLoc', 'GridOrient', 'GridAtlas', 'SurfaceFile', 'MEGMethod', 'EEGMethod', 'ECOGMethod', 'SEEGMethod', 'HeadModelType');
        % Apply current SSP projectors
        if ~isempty(ChannelMat.Projector)
            % Rebuild projector in the expanded form (I-UUt)
            Proj = process_ssp2('BuildProjector', ChannelMat.Projector, [1 2]);
            % Apply projectors
            if ~isempty(Proj)
                % Get all sensors for which the gain matrix was successfully computed
                iGainSensors = find(sum(isnan(HeadModel.Gain), 2) == 0);
                % Apply projectors to gain matrix
                HeadModel.Gain(iGainSensors,:) = Proj(iGainSensors,iGainSensors) * HeadModel.Gain(iGainSensors,:);
                % Apply SSPs on both sides of the noise covariance matrix
                NoiseCovMat.NoiseCov = Proj * NoiseCovMat.NoiseCov * Proj';
                if ~isempty(DataCovMat)
                    DataCovMat.NoiseCov = Proj * DataCovMat.NoiseCov * Proj';
                end
            end
        end
        % Select only good channels
        HeadModel.Gain = HeadModel.Gain(GoodChannel, :);
        % Apply average reference: separately SEEG, ECOG, EEG
        if any(ismember(unique({ChannelMat.Channel.Type}), {'EEG','ECOG','SEEG'}))
            % Create average reference montage
            sMontage = panel_montage('GetMontageAvgRef', [], ChannelMat.Channel(GoodChannel), ChannelFlag(GoodChannel), 0);
            HeadModel.Gain = sMontage.Matrix * HeadModel.Gain;
            % Apply average reference operator on both sides of the noise covariance matrix
            NoiseCovMat.NoiseCov(GoodChannel, GoodChannel) = sMontage.Matrix * NoiseCovMat.NoiseCov(GoodChannel, GoodChannel) * sMontage.Matrix';
            if ~isempty(DataCovMat)
                DataCovMat.NoiseCov(GoodChannel, GoodChannel) = sMontage.Matrix * DataCovMat.NoiseCov(GoodChannel, GoodChannel) * sMontage.Matrix';
            end
        end
        % Copy initial head model
        HeadModelInit = HeadModel;
        % Get number of sources
        nSources =  size(HeadModelInit.Gain,2) / 3;

        % ===== MIXED HEADMODEL =====
        if strcmpi(HeadModelInit.HeadModelType, 'mixed') && ~isempty(HeadModel.GridAtlas) && ~isempty(HeadModel.GridAtlas(1).Scouts)
            % Only supported for wMNE
            if ~ismember(OPTIONS.InverseMethod, {'minnorm', 'gls', 'lcmv'})
                errMessage = [errMessage 'The mixed headmodel is currently only supported for the following inverse solutions: Minimum norm, dipole fitting, beamformer.' 10];
                break;
            end
            % Initialize variable
            HeadModel.Gain      = [];
            HeadModel.GridAtlas = [];
            HeadModel           = repmat(HeadModel, 1, length(HeadModelInit.GridAtlas(1).Scouts));
            iVert2Grid = [];
            iAllGrid   = [];
            iAllSource = [];
            iOffset    = 0;
            % Split the head model in multiple entries
            for iScout = 1:length(HeadModelInit.GridAtlas(1).Scouts)
                % Get indices in the Gain matrix
                sScout = HeadModelInit.GridAtlas(1).Scouts(iScout);
                iGainRows = sort([3*sScout.GridRows-2, 3*sScout.GridRows-1, 3*sScout.GridRows]);
                % Create the headmodel structure for the current region
                HeadModel(iScout).Gain       = HeadModelInit.Gain(:, iGainRows);
                HeadModel(iScout).GridLoc    = HeadModelInit.GridLoc(sScout.GridRows, :);
                HeadModel(iScout).GridOrient = HeadModelInit.GridOrient(sScout.GridRows, :);
                switch (sScout.Region(3))
                    case 'C',  OPTIONS.SourceOrient{iScout} = 'fixed';  nComp = 1;
                    case 'U',  OPTIONS.SourceOrient{iScout} = 'free';   nComp = 3;
                    case 'L',  OPTIONS.SourceOrient{iScout} = 'loose';  nComp = 3;
                end
                % In the case of a surface region, add the match of the vertices in the cortex surface and the GridLoc matrix
                if strcmpi(sScout.Region(2), 'S')
                    iVert2Grid = [iVert2Grid; sScout.Vertices(:), sScout.GridRows(:)];
                end
                % Add to the scout definition the indices in the ImageGrid
                iAllGrid   = [iAllGrid,   reshape(repmat(sScout.GridRows,nComp,1), 1, [])];
                iAllSource = [iAllSource, iOffset + (1:nComp*length(sScout.GridRows))];
                iOffset = iOffset + nComp*length(sScout.GridRows);
            end
            % Create sparse conversion matrices between indices
            if ~isempty(iVert2Grid)
                HeadModelInit.GridAtlas(1).Vert2Grid = logical(sparse(iVert2Grid(:,2), iVert2Grid(:,1), ones(size(iVert2Grid,1),1)));
            else
                HeadModelInit.GridAtlas(1).Vert2Grid = [];
            end
            if ~isempty(iAllSource)
                HeadModelInit.GridAtlas(1).Grid2Source = logical(sparse(iAllSource, iAllGrid, ones(size(iAllSource))));
            else
                HeadModelInit.GridAtlas(1).Grid2Source = [];
            end
        end

        %% ===== COMPUTE INVERSE SOLUTION =====
        bst_progress('text', 'Estimating sources...');
        bst_progress('inc', 1);
        % NoiseCov: keep only the good channels
        OPTIONS.NoiseCovMat = NoiseCovMat;
        OPTIONS.NoiseCovMat.NoiseCov = OPTIONS.NoiseCovMat.NoiseCov(GoodChannel, GoodChannel);
        if isfield(OPTIONS.NoiseCovMat, 'FourthMoment') && ~isempty(OPTIONS.NoiseCovMat.FourthMoment)
            OPTIONS.NoiseCovMat.FourthMoment = OPTIONS.NoiseCovMat.FourthMoment(GoodChannel, GoodChannel);
        end
        if isfield(OPTIONS.NoiseCovMat, 'nSamples') && ~isempty(OPTIONS.NoiseCovMat.nSamples)
            OPTIONS.NoiseCovMat.nSamples = OPTIONS.NoiseCovMat.nSamples(GoodChannel, GoodChannel);
        end
        % DataCov: keep only the good channels
        if ~isempty(DataCovMat)
            OPTIONS.DataCovMat = DataCovMat;
            OPTIONS.DataCovMat.NoiseCov     = OPTIONS.DataCovMat.NoiseCov(GoodChannel, GoodChannel);
            OPTIONS.DataCovMat.FourthMoment = OPTIONS.DataCovMat.FourthMoment(GoodChannel, GoodChannel);
            OPTIONS.DataCovMat.nSamples     = OPTIONS.DataCovMat.nSamples(GoodChannel, GoodChannel);
        end
        % Get channels types
        OPTIONS.ChannelTypes = {ChannelMat.Channel(GoodChannel).Type};
        % Switch depending on the selected inverse method
        switch( OPTIONS.InverseMethod )
            case {'minnorm', 'gls', 'lcmv'}
                % Call John's wmne function
                % NOTE: The output HeadModel param is used here in return to save LOTS of memory in the bst_inverse_linear_2016 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 reference".
                try
                    [Results, OPTIONS] = bst_inverse_linear_2016(HeadModel, OPTIONS);
                catch e
                    if bst_get('MatlabVersion') == 904
                        errMsg = ['Note: Matlab 2018a changed the behavior of the SVD() function. ' ...
                            10 'If issues arise, we recommend using another version.'];
                        e = MException(e.identifier, [e.message 10 10 errMsg]);
                        throw(e);
                    else
                        rethrow(e);
                    end
                end
            case 'mem'
                % Add options needed by the MEM functions
                OPTIONS.DataFile      = DataFile;
                OPTIONS.DataTime      = Time;
                OPTIONS.Channel       = ChannelMat.Channel(GoodChannel);
                OPTIONS.Data          = DataMat.F(GoodChannel,:);
                OPTIONS.ChannelFlag   = ChannelFlag(GoodChannel);
                OPTIONS.ResultFile    = [];
                OPTIONS.HeadModelFile = HeadModelFile;
                OPTIONS.GoodChannel   = GoodChannel;
                OPTIONS.FunctionName  = 'mem';
                % Call the mem solver
                [Results, OPTIONS] = be_main(HeadModel, OPTIONS);
                Results.nComponents = round(max(size(Results.ImageGridAmp,1),size(Results.ImagingKernel,1)) / nSources);
                % Get outputs
                DataFile = OPTIONS.DataFile;
                Time     = OPTIONS.DataTime;
            otherwise
                error('Unknown method');
        end
        % Error handling
        if isempty(Results)
            errMessage = [errMessage 'The inverse function returned an empty structure.' 10];
            break;
        end
        % Copy outputs to a standard results structure
        ResultsMat = db_template('resultsmat');
        ResultsMat = struct_copy_fields(ResultsMat, Results, 1);

        % ===== COMPUTE FULL RESULTS =====
        % Full results
        if (OPTIONS.ComputeKernel == 0) && ~isempty(ResultsMat.ImagingKernel) && ~isempty(DataFile)
            % Load data
            DataMat = in_bst_data(DataFile, 'F');
            % Incompatible options: Full results + raw files (impossible to view after)
            if isstruct(DataMat.F)
                errMessage = [errMessage 'Cannot compute full results for raw files: import the files first or compute an inversion kernel only.' 10];
                break;
            end
            % Multiply inversion kernel with the recordings
            ResultsMat.ImageGridAmp = ResultsMat.ImagingKernel * DataMat.F(GoodChannel, :);
            ResultsMat.ImagingKernel = [];
        end
        % If kernel only on raw data: do not save the Time vector
        if ~isempty(ResultsMat.ImagingKernel) && ~isShared && isRaw
            Time = [];
        % Check if no time dimension
        elseif ~isempty(ResultsMat.ImageGridAmp) && (size(ResultsMat.ImageGridAmp, 2) == 1)
            % If the output is not a timecourse: replicate the results to get two time instants
            ResultsMat.ImageGridAmp = repmat(ResultsMat.ImageGridAmp, [1,2]);
            % Keep only the first and last time instants
            Time = [Time(1), Time(end)];
        end

        %% ===== SAVE RESULTS FILE =====
        bst_progress('text', 'Saving results...');
        bst_progress('inc', 1);
        % ===== OUTPUT FILENAME =====
        % Add method name
        strMethod = methodTag;
        % Add modality
        for i = 1:length(OPTIONS.DataTypes)
            strMethod = [strMethod, '_', file_standardize(OPTIONS.DataTypes{i})];
        end
        % Add kernel tag
        if OPTIONS.ComputeKernel
            strMethod = [strMethod, '_KERNEL'];
        end
        % Output folder
        if isempty(DataFile)
            OutputDir = bst_fileparts(file_fullpath(ChannelFile));
        else
            OutputDir = bst_fileparts(file_fullpath(DataFile));
        end
        % Output filename
        ResultFile = bst_process('GetNewFilename', OutputDir, ['results_', strMethod]);

        % ===== CREATE FILE STRUCTURE =====
        ResultsMat.Comment       = [OPTIONS.Comment ' 2016'];
        ResultsMat.Function      = OPTIONS.FunctionName;
        ResultsMat.Time          = Time;
        ResultsMat.DataFile      = DataFile;
        ResultsMat.HeadModelFile = HeadModelFile;
        ResultsMat.HeadModelType = HeadModelInit.HeadModelType;
        ResultsMat.ChannelFlag   = ChannelFlag;
        ResultsMat.GoodChannel   = GoodChannel;
        ResultsMat.SurfaceFile   = file_short(HeadModelInit.SurfaceFile);
        switch lower(ResultsMat.HeadModelType)
            case 'volume'
                ResultsMat.GridLoc    = HeadModelInit.GridLoc;
                % ResultsMat.GridOrient = [];
            case 'surface'
                ResultsMat.GridLoc    = [];
                % ResultsMat.GridOrient = [];    % THE ORIENTATION CAN BE RETURNED BY THE INVERSE METHOD ('optim')
            case 'mixed'
                ResultsMat.GridLoc    = HeadModelInit.GridLoc;
                ResultsMat.GridOrient = HeadModelInit.GridOrient;
        end
        ResultsMat.GridAtlas = HeadModelInit.GridAtlas;
        ResultsMat.nAvg      = nAvg;
        ResultsMat.Leff      = Leff;
        ResultsMat.Options   = OPTIONS;
        % History
        ResultsMat = bst_history('add', ResultsMat, 'compute', ['Source estimation: ' OPTIONS.InverseMethod]);
        % Make file comment unique
        if ~isempty(sStudy.Result)
            ResultsMat.Comment = file_unique(ResultsMat.Comment, {sStudy.Result.Comment});
        end
        % Save new file structure
        bst_save(ResultFile, ResultsMat, 'v6');

        % ===== REGISTER NEW FILE =====
        % Create new results structure
        newResult = db_template('results');
        newResult.Comment       = ResultsMat.Comment;
        newResult.FileName      = file_short(ResultFile);
        newResult.DataFile      = DataFile;
        newResult.isLink        = 0;
        newResult.HeadModelType = ResultsMat.HeadModelType;
        % Add new entry to the database
        iResult = length(sStudy.Result) + 1;
        sStudy.Result(iResult) = newResult;
        % Update Brainstorm database
        bst_set('Study', iStudy, sStudy);

        % ===== UPDATE DISPLAY =====
        % Update tree
        panel_protocols('UpdateNode', 'Study', iStudy);
        % Update links
        if isShared
            if isDefaultStudy
                % If added to a 'default_study' node: need to update results links
                OutputLinks = db_links('Subject', sStudy.BrainStormSubject);
                % Update whole tree display
                panel_protocols('UpdateTree');
            else
                % Update links to the new results file
                OutputLinks = db_links('Study', iStudy);
                % Update display of the study node
                panel_protocols('UpdateNode', 'Study', iStudy);
            end
            % Find in the links the ones that are based on the node that was just calculated
            isNewLink = ~cellfun(@(c)isempty(strfind(c, newResult.FileName)), OutputLinks);
            OutputFiles = cat(2, OutputFiles, OutputLinks(isNewLink));
        else
            % Store output filename
            OutputFiles{end+1} = newResult.FileName;
        end
        % Expand data node
        panel_protocols('SelectNode', [], newResult.FileName);
    end
    % Save database
    db_save();
    % Hide progress bar
    bst_progress('stop');
end

%% ===== GET MODALITY COMMENT =====
function Comment = GetModalityComment(Modalities)
    % Replace "MEG GRAD+MEG MAG" with "MEG ALL"
    if all(ismember({'MEG GRAD', 'MEG MAG'}, Modalities))
        Modalities = setdiff(Modalities, {'MEG GRAD', 'MEG MAG'});
        Modalities{end+1} = 'MEG ALL';
    end
    % Loop to build comment
    Comment = '';
    for im = 1:length(Modalities)
        if (im >= 2)
            Comment = [Comment, '+'];
        end
        Comment = [Comment, Modalities{im}];
    end
end

author: brainstorm-tools


text:
function [F1,Se,PPV,Nb] = run_sqi(refqrs,testqrs,thres,margin,windowlen,fs)
% sqi = run_sqi(refqrs,testqrs,thres,margin,windowlen,fs)
% compare two sets of annotation with one as the reference (refqrs) and one
% as the test (testqrs)
%
% inputs
%     refqrs:       reference qrs annotation (in sec)
%     testqrs:      test qrs annotations (in sec)
%     thres:        threshold (in sec,default 0.05s)
%     margin:       margin time not include in comparison (in sec,default 2s)
%     windowlen:    length of the comparison window (in sec,default 60s)
%     fs:           sampling frequency
%
% output
%     sqi: match proportion according to some criteria you can change
%     depending on what you are looking for (can be Se, PPV or F1 measure).
%     See at the end of the function.
%
% When using this work, then please cite [1] and [2]:
%     [1] Behar Joachim, Oster Julien, Qiao Li, Clifford Gari D. Signal Quality
%     During Arrhythmia and its Application to False Alarm Reduction.
%     IEEE Transactions on Biomedical Engineering. 60(6). 1660-6. 2013.
%
%     [2] Li, Qiao, Roger G. Mark, and Gari D. Clifford. "Robust heart rate estimation
%     from multiple asynchronous noisy sources using signal quality indices and
%     a Kalman filter." Physiological measurement 29.1 (2008): 15.
%
% PCinCC2014, version 1.0, June 2014
% Released under the GNU General Public License
%
% Copyright (C) 2014  Joachim Behar
% Oxford university, Intelligent Patient Monitoring Group - Oxford 2014
% joachim.behar@gmail.com
%
% Updates:
% 07-02-2014
% JB- testes with Octave -> running OK
%
% 02-10-2014
% Bug fix: Dealing with annotations close to the border of the search
% window. Lines 70 - 100
% David Springer
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.

% == managing inputs
if nargin<2; error('bsqi: wrong number of input arguments \n'); end;
if nargin<3; thres=0.05; end;
if nargin<4; margin=2; end;
if nargin<5; windowlen=60; end;
if nargin<6; fs=1000; end;

if size(refqrs,1)>size(refqrs,2); refqrs=refqrs';end
if size(testqrs,1)>size(testqrs,2); testqrs=testqrs';end

start = margin*fs;
stop = (windowlen-margin)*fs;
refqrs = refqrs.*fs; % convert into samples from time
testqrs = testqrs.*fs; % convert into samples from time

try
    refqrs  = refqrs(refqrs>start & refqrs<stop)'; % reference annotations
    testqrs = testqrs(testqrs>start & testqrs<stop)'; % test annotations

    if ~isempty(refqrs)

        NB_REF  = length(refqrs);
        NB_TEST = length(testqrs);

        % == removing borders problems
        indbord = find(refqrs<thres*fs | refqrs>(windowlen-thres)*fs); % reference QRS at the border
        if ~isempty(indbord)
            [IndMatchBord,DistQRSbord] = dsearchn(testqrs,refqrs(indbord));
            Indeces_below_threshold = DistQRSbord<thres*fs; %Added line to find indices of annotation of interest (refqrs) below threshold (02-10-14)
            IndMatchBord = IndMatchBord(Indeces_below_threshold);
            NB_QRS_BORD = length(indbord); % QRS at the border of the window (0,1,2)
            NB_MATCHING = length(IndMatchBord); % nb of corresponding mathing QRS (0,1,2)
            if isempty(IndMatchBord);
                refqrs(indbord) = [];
            elseif NB_MATCHING<NB_QRS_BORD
                refqrs(indbord(~Indeces_below_threshold)) = []; %Removing other indices not below threshold (02-10-14)
            end
        end

        indbord = find(testqrs<thres*fs | testqrs>(windowlen-thres)*fs); % reference QRS at the border
        if ~isempty(indbord)
            [IndMatchBord,DistQRSbord] = dsearchn(refqrs,testqrs(indbord));
            Indeces_below_threshold = DistQRSbord<thres*fs; %Added line to find indices of annotation of interest (testqrs) below threshold (02-10-14)
            IndMatchBord = IndMatchBord(Indeces_below_threshold);
            NB_QRS_BORD = length(indbord); % QRS at the border of the window (0,1,2)
            NB_MATCHING = length(IndMatchBord); % nb of corresponding mathing QRS (0,1,2)
            if isempty(IndMatchBord);
                testqrs(indbord) = [];
            elseif NB_MATCHING<NB_QRS_BORD
                testqrs(indbord(~Indeces_below_threshold)) = []; %Removing other indices not below threshold (02-10-14)
            end
        end

        % == core function
        [IndMatch,Dist] = dsearchn(refqrs,testqrs);         % closest ref for each point in test qrs
        IndMatchInWindow = IndMatch(Dist<thres*fs);         % keep only the ones within a certain window
        NB_MATCH_UNIQUE = length(unique(IndMatchInWindow)); % how many unique matching
        TP = NB_MATCH_UNIQUE;                               % number of identified ref QRS
        FN = NB_REF-TP;                                     % number of missed ref QRS
        FP = NB_TEST-TP;                                    % how many extra detection?
        Se  = TP/(TP+FN);
        PPV = TP/(FP+TP);
        F1 = 2*Se*PPV/(Se+PPV);                             % accuracy measure

        Nb.TP = TP;
        Nb.FN = FN;
        Nb.FP = FP;
    else
        F1=[];Se=[];PPV=[];Nb=[];
    end
catch
    F1=[];Se=[];PPV=[];Nb=[];
end

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function cadj=ref_spreadadj(coef);
%REF_SPREADADJ  Symbol of adjoint preading function.
%   Usage: cadj=ref_spreadadj(c);
%
%   cadj=SPREADADJ(c) will compute the symbol cadj of the spreading
%   operator that is the adjoint of the spreading operator with symbol c.
%
%   The algorithm converts the symbol to the matrix representation,
%   adjoints its, and finds its spreading function.

L=size(coef,1);

T=tfmat('spread',coef);
cadj=spreadfun(T');

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
addpath('mex');

% we provide two sequences "car" and "table"
% example = 'table';
example = 'car';

% load the two frames
im1 = im2double(imread([example '1.jpg']));
im2 = im2double(imread([example '2.jpg']));

% im1 = imresize(im1,0.5,'bicubic');
% im2 = imresize(im2,0.5,'bicubic');

% set optical flow parameters (see Coarse2FineTwoFrames.m for the definition of the parameters)
alpha = 0.012;
ratio = 0.75;
minWidth = 20;
nOuterFPIterations = 7;
nInnerFPIterations = 1;
nSORIterations = 30;

para = [alpha,ratio,minWidth,nOuterFPIterations,nInnerFPIterations,nSORIterations];

% this is the core part of calling the mexed dll file for computing optical flow
% it also returns the time that is needed for two-frame estimation
tic;
[vx,vy,warpI2] = Coarse2FineTwoFrames(im1,im2,para);
toc

figure;imshow(im1);figure;imshow(warpI2);

% output gif
clear volume;
volume(:,:,:,1) = im1;
volume(:,:,:,2) = im2;
if exist('output','dir')~=7
    mkdir('output');
end
frame2gif(volume,fullfile('output',[example '_input.gif']));
volume(:,:,:,2) = warpI2;
frame2gif(volume,fullfile('output',[example '_warp.gif']));

% visualize flow field
clear flow;
flow(:,:,1) = vx;
flow(:,:,2) = vy;
imflow = flowToColor(flow);

figure;imshow(imflow);
imwrite(imflow,fullfile('output',[example '_flow.jpg']),'quality',100);

author: qingsenyangit

repo: AHDRNet


text:
function fixPSlinestyle(varargin)

%FIXPSLINESTYLE Fix line styles in exported post script files
%
% FIXPSLINESTYLE(FILENAME) fixes the line styles in the postscript file
%   FILENAME. The file will be over-written. This takes a .PS or .EPS file
%   and fixes the dotted and dashed line styles to be a little bit more
%   esthetically pleasing. It fixes the four default line styles (line,
%   dotted, dashed, dashdot).
%
% FIXPSLINESTYLE(FILENAME, NEWFILENAME) creates a new file NEWFILENAME.
%
%   This is meant to be used with postscript files created by MATLAB
%   (print, export).
%
% Example:
%   x  = 1:.1:10;
%   y1 = sin(x);
%   y2 = cos(x);
%   h  = plot(x, y1, '--', x, y2, '-.');
%   set(h, 'LineWidth', 2);
%   grid on;
%   legend('line 1', 'line2');
%
%   print -depsc test.eps
%   fixPSlinestyle('test.eps', 'fixed_test.eps');
%
% See also PRINT.

% Copyright 2005-2010 The MathWorks, Inc.

% Error checking
error(nargchk(1, 2, nargin));
if ~ischar(varargin{1}) || (nargin == 2 && ~ischar(varargin{2}))
  error('Input arguments must be file names (char).');
end

% Make sure the files specified are postscript files
[p1, n1, e1] = fileparts(varargin{1});
if isempty(e1) || ~ismember(lower(e1), {'.ps', '.eps'})
  error('The extension has to be .ps or .eps');
end

% Open file and read it in
fid = fopen(varargin{1}, 'r');
str = fread(fid);
str = char(str');
fclose(fid);

% Find where the line types are defined
id   = findstr(str, '% line types:');
str1 = str(1:id-1);
[line1     , remline ] = strtok(str(id:end), '/');
[replacestr, remline2] = strtok(remline    , '%');

% Define the new line styles
solidLine   = sprintf('/SO { [] 0 setdash } bdef\n');
dotLine     = sprintf('/DO { [3 dpi2point mul 3 dpi2point mul] 0 setdash } bdef\n');
dashedLine  = sprintf('/DA { [6 dpi2point mul] 0 setdash } bdef\n');
dashdotLine = sprintf('/DD { [2 dpi2point mul 2 dpi2point mul 6 dpi2point mul 2 dpi2point mul] 0 setdash } bdef\n');

% Construct the new file with the new line style definitions
newText     = [str1, line1, solidLine, dotLine, dashedLine, dashdotLine, remline2];

% Check for output file name
if nargin == 2
  [p2, n2, e2] = fileparts(varargin{2});
  if isempty(e2)
    fname = fullfile(p2, [n2, e1]);
  else
    if strcmpi(e1, e2)
      fname = varargin{2};
    else
      error('Output file must have same file extension.');
    end
  end
else % if not defined, over-write
  fname = varargin{1};
end

% Write out to file
fid = fopen(fname, 'w');
fprintf(fid, '%s', newText);
fclose(fid);

author: masumhabib

repo: PlotPub

save_path: github-repos/MATLAB/masumhabib-PlotPub


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function decentr_state = readDecentrStateFromOptG2oFiles(...
    g2o_dir, decentr_state, suffix)
% Inspired by Luca Carlone's
% https://bitbucket.org/lucacarlone/pgo3d-duality-opencode/src/ebb6e1b8cebaad7f2aaf581b1d0c0bad737faebb/lib/readG2oDataset3D.m?at=master&fileviewer=file-view-default

nr_robots = numel(decentr_state);

for robot_i = 1:nr_robots
    file_id = fopen(...
        [g2o_dir '/' num2str(robot_i - 1) suffix '.g2o'], 'r');

    while true
        line = fgets(file_id);
        data = textscan(line, '%s %d64 %f %f %f %f %f %f %f');
        if ~strcmp(data{1}, 'VERTEX_SE3:QUAT')
            break
        end
        [robot_i_val, frame_i] = gtsamFrameIdToIndices(data{2});
        assert(robot_i_val == robot_i);

        x = data{3}; y = data{4}; z = data{5};
        qx = data{6}; qy = data{7}; qz = data{8}; qw = data{9};

        Sim_W_C = eye(4);
        Sim_W_C(1:3, 4) = [x y z]';
        q = [qw, qx qy, qz]';
        if(abs(norm(q)-1) > 1e-3)
            norm(q)
            error('Quaternion has not unit norm');
        else
            q = q/norm(q); % we normalize anyway
        end

        Sim_W_C(1:3, 1:3) = fixR(quat2rot(q));

        decentr_state{robot_i}.Sim_O_C{frame_i} = ...
            Sim_W_C;
    end
end

end

author: uzh-rpg


text:
function [fh yMin, yMax] = tapas_uniqc_plot_abs_image(Y,iDynSli,fh, yMin, yMax, ...
    colorMap, colorBar)
%simple plotting routine for one dynamic within a 3D-dataset
%
%   [fh yMin, yMax] = tapas_uniqc_plot_image_diagnostics(Y,iDynSli,fh, yMin, yMax)
%
% IN
%   Y           [nVoxelX nVoxelY, nSlices, nVolumes] real valued data matrix
%               OR
%               [nSamples, nCoils, nDynSlis] real-valued data matrix
%   iDynSli     index of which dynamic/slice shall be plotted
%   fh          figure handle to be plotted into; if empty or missing, new figure is
%               created
%   yMin        [nPlots,1] min value in ylim, one for each set of coils, if
%               empty or missing, newly created from data limits
%   yMax        [nPlots,1] max value in ylim, one for each set of coils, if
%               empty or missing, newly created from data limits
%
% OUT
%   fh          figure handle plotted into
%   yMin        [nPlots,2] min value in ylim, one for each set of coils,
%               abs and phase
%   yMax        [nPlots,2] max value in ylim, one for each set of coils
%               abs and phase
% EXAMPLE
%   [fh yMin, yMax] = tapas_uniqc_plot_image_diagnostics(Y,iDynSli,fh, yMin, yMax, coilPlots)
%
%   See also plotTrajDiagnostics guiTrajDiagnostics

% Author: Lars Kasper
% Created: 2013-01-06
% Copyright 2013 Institute for Biomedical Engineering, ETH/Uni Zurich.

if nargin < 3 || isempty(fh)
    fh = figure('Name','Video of ImageDiagnostics', 'WindowStyle', 'normal');
else
    figure(fh);
end

% determine plot limits
if nargin < 5 || isempty(yMin)
    validIndices = find(~isinf(Y) & ~isnan(Y));
    yMin = min(Y(validIndices));
    yMax = max(Y(validIndices));
end

if nargin < 6
    colorMap = 'gray';
end

if nargin < 7
    colorBar = 'off';
end

doPlotColorBar = strcmpi(colorBar, 'on');

% plot abs data always
imagesc(Y(:,:,iDynSli));
%colormap gray; axis image;
colormap(colorMap); axis image;

caxis([yMin, yMax]);

if doPlotColorBar
    colorbar;
end

stringTitle = sprintf('abs, iDynSli = %d', iDynSli);
if exist('suptitle', 'builtin')
    suptitle(stringTitle);
else
    title(stringTitle);
end

end

author: translationalneuromodeling

repo: tapas


text:
function ov = regionOverlap(region1, region2, area)
% compute intersection over union of two regions
% ov = regionOverlap(region1, region2, area)

area = area(:);
r1 = false(numel(area), 1);
r2 = false(numel(area), 1);
r1(region1) = true;
r2(region2) = true;
ov = sum(area(r1 & r2)) / sum(area(r1 | r2));

author: Cloud-CV

repo: object-proposals


text:
function [ result ] = get_subwindow( input_image, index, resx, resy, spacingx, spacingy )
% index is 0-N
%GET_SCRIPT Summary of this function goes here
%   Detailed explanation goes here

	maxPerRow = floor((size(input_image,2)-resx+spacingx) / spacingx);
    maxRows = floor((size(input_image,1) - resy+spacingy )/ spacingy);
    currentRow = floor(index / (maxPerRow));

    if (currentRow < maxRows)

        cutX = (mod(index,maxPerRow)*spacingx);
        cutY = (currentRow * spacingy);

        result = input_image(cutY+1:(cutY+resy),cutX+1:(cutX+resx),:);

    else
        result = 'error';
        %error('index too high, image overflow');
    end

end


text:
%% |plotboxpos.m|: Plot box (visible) position
% Author: Kelly Kearney
%
% This repository includes the code for the |plotboxpos.m| Matlab function.
%
% This function returns the position of the plotted region of an axis,
% which may differ from the actual axis position, depending on the axis
% limits, data aspect ratio, and plot box aspect ratio.  The position is
% returned in the same units as the those used to define the axis itself.
% This function can only be used for a 2D plot.
%
%% Getting started
%
% *Prerequisites*
%
% This function requires Matlab R14 or later.
%
% *Downloading and installation*
%
% This code can be downloaded from <https://github.com/kakearney/plotboxpos-pkg/ Github>
% or the
% <http://www.mathworks.com/matlabcentral/fileexchange/9615-plotboxpos
% MatlabCentral File Exchange>.  The File Exchange entry is updated daily
% from the GitHub repository.
%
% *Matlab Search Path*
%
% The following folders need to be added to your Matlab Search path (via
% |addpath|, |pathtool|, etc.):
%
%   plotboxpos-pkg/plotboxpos

%% Syntax
%
% |pos = plotboxpos(h)| returns the 1 x 4 position vector |pos| for
% the axis with handle |h|.  The units of |pos| will match those of the
% axis' parent object (typically the figure).

%% Examples
%
% We start by plotting a circle, changing the axis aspect ratio to be
% 1:1 and the axis limits to be tight to the data.

h.ax(1) = axes;
h.ex = rectangle('position', [0 0 1 1], 'curvature', [1 1]);
axis tight equal;
box on;

%%
% The axis |'Position'| property still corresponds to the full potential axis
% space, discounting the axis modifications we just made.

pos1 = get(h.ax(1), 'position');
annotation('rectangle', pos1, 'edgecolor', 'b', 'linestyle', '--');

%%
% The |plotboxpos| function returns the position the axis is actually using
% with its current axis ratio settings;

pos2 = plotboxpos(h.ax(1));
annotation('rectangle', pos2, 'edgecolor', 'r', 'linestyle', '-.');

%% Contributions
%
% Community contributions to this package are welcome!
%
% To report bugs, please submit
% <https://github.com/kakearney/plotboxpos-pkg/issues an issue> on GitHub and
% include:
%
% * your operating system
% * your version of Matlab and all relevant toolboxes (type |ver| at the Matlab command line to get this info)
% * code/data to reproduce the error or buggy behavior, and the full text of any error messages received
%
% Please also feel free to submit enhancement requests, or to send pull
% requests (via GitHub) for bug fixes or new features.
%
% I do monitor the MatlabCentral FileExchange entry for any issues raised
% in the comments, but would prefer to track issues on GitHub.
%

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function varargout = ellipsoidMesh(elli, varargin)
%ELLIPSOIDMESH Convert a 3D ellipsoid to face-vertex mesh representation
%
%   [V, F] = ellipsoidMesh(ELLI)
%   ELLI is given by:
%   [XC YC ZC  A B C  PHI THETA PSI],
%   where (XC, YC, ZC) is the ellipsoid center, A, B and C are the half
%   lengths of the ellipsoid main axes, and PHI THETA PSI are Euler angles
%   representing ellipsoid orientation, in degrees.
%
%
%   See also
%   meshes3d, drawEllipsoid, sphereMesh, inertiaEllipsoid
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-03-12,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

%% Default values

% number of meridians
nPhi    = 32;

% number of parallels
nTheta  = 16;

%% Extract input arguments

% Parse the input (try to extract center coordinates and radius)
if nargin == 0
    % no input: assumes ellipsoid with default shape
    elli = [0 0 0 5 4 3 0 0 0];
end

% default set of options for drawing meshes
options = {'FaceColor', 'g', 'linestyle', 'none'};

while length(varargin) > 1
    switch lower(varargin{1})
        case 'nphi'
            nPhi = varargin{2};

        case 'ntheta'
            nTheta = varargin{2};

        otherwise
            % assumes this is drawing option
            options = [options varargin(1:2)]; %#ok<AGROW>
    end

    varargin(1:2) = [];
end

%% Parse numerical inputs

% Extract ellipsoid parameters
xc  = elli(:,1);
yc  = elli(:,2);
zc  = elli(:,3);
a   = elli(:,4);
b   = elli(:,5);
c   = elli(:,6);
k   = pi / 180;
ellPhi   = elli(:,7) * k;
ellTheta = elli(:,8) * k;
ellPsi   = elli(:,9) * k;

%% Coordinates computation

% convert unit basis to ellipsoid basis
sca     = createScaling3d(a, b, c);
rotZ    = createRotationOz(ellPhi);
rotY    = createRotationOy(ellTheta);
rotX    = createRotationOx(ellPsi);
tra     = createTranslation3d([xc yc zc]);

% concatenate transforms
trans   = tra * rotZ * rotY * rotX * sca;

%% parametrisation of ellipsoid

% spherical coordinates
theta   = linspace(0, pi, nTheta+1);
phi     = linspace(0, 2*pi, nPhi+1);

% convert to cartesian coordinates
sintheta = sin(theta);
x = cos(phi') * sintheta;
y = sin(phi') * sintheta;
z = ones(length(phi),1) * cos(theta);

% transform mesh vertices
[x, y, z] = transformPoint3d(x, y, z, trans);

% convert to FV mesh
[vertices, faces] = surfToMesh(x, y, z, 'xPeriodic', false, 'yPeriodic', true);

% format output
varargout = formatMeshOutput(nargout, vertices, faces);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [ y, j, f, ierror ] = yjf_check_common ( y, j, f )

%*****************************************************************************80
%
%% YJF_CHECK_COMMON normalizes a Common YJF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, J, real F, the YJF date.
%
%    Output, integer IERROR, nonzero if there was an error.
%
  [ y, j, ierror ] = yj_check_common ( y, j );

  if ( ierror ~= 0 )
    return
  end
%
%  Force the fraction to lie between 0 and 1.
%
  while ( f < 0.0 )
    f = f + 1.0;
    j = j - 1;
  end

  while ( 1.0 <= f )
    f = f - 1.0;
    j = j + 1;
  end

  return
end


text:
function varargout = lvd_EditTankToTankConnsGUI(varargin)
% LVD_EDITTANKTOTANKCONNSGUI MATLAB code for lvd_EditTankToTankConnsGUI.fig
%      LVD_EDITTANKTOTANKCONNSGUI, by itself, creates a new LVD_EDITTANKTOTANKCONNSGUI or raises the existing
%      singleton*.
%
%      H = LVD_EDITTANKTOTANKCONNSGUI returns the handle to a new LVD_EDITTANKTOTANKCONNSGUI or the handle to
%      the existing singleton*.
%
%      LVD_EDITTANKTOTANKCONNSGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LVD_EDITTANKTOTANKCONNSGUI.M with the given input arguments.
%
%      LVD_EDITTANKTOTANKCONNSGUI('Property','Value',...) creates a new LVD_EDITTANKTOTANKCONNSGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before lvd_EditTankToTankConnsGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to lvd_EditTankToTankConnsGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help lvd_EditTankToTankConnsGUI

% Last Modified by GUIDE v2.5 21-Jan-2019 11:37:06

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_EditTankToTankConnsGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_EditTankToTankConnsGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before lvd_EditTankToTankConnsGUI is made visible.
function lvd_EditTankToTankConnsGUI_OpeningFcn(hObject, eventdata, handles, varargin)
    % This function has no output args, see OutputFcn.
    % hObject    handle to figure
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    % varargin   command line arguments to lvd_EditTankToTankConnsGUI (see VARARGIN)

    % Choose default command line output for lvd_EditTankToTankConnsGUI
    handles.output = hObject;

    lvdData = varargin{1};
    setappdata(hObject,'lvdData',lvdData);

    populateGUI(handles, lvdData);

    % Update handles structure
    guidata(hObject, handles);

    % UIWAIT makes lvd_EditTankToTankConnsGUI wait for user response (see UIRESUME)
    uiwait(handles.lvd_EditTankToTanksConnsGUI);

function populateGUI(handles, lvdData)
    set(handles.connsListBox,'String',lvdData.launchVehicle.getTankToTankConnectionsListBoxStr());

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditTankToTankConnsGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = [];
    else
        close(handles.lvd_EditTankToTanksConnsGUI);
    end

% --- Executes on button press in closeButton.
function closeButton_Callback(hObject, eventdata, handles)
% hObject    handle to closeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    uiresume(handles.lvd_EditTankToTanksConnsGUI);

% --- Executes on selection change in connsListBox.
function connsListBox_Callback(hObject, eventdata, handles)
% hObject    handle to connsListBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns connsListBox contents as cell array
%        contents{get(hObject,'Value')} returns selected item from connsListBox
    if(strcmpi(get(handles.lvd_EditTankToTanksConnsGUI,'SelectionType'),'open'))
        lvdData = getappdata(handles.lvd_EditTankToTanksConnsGUI,'lvdData');
        lv = lvdData.launchVehicle;

        [~,t2TConns] = lvdData.launchVehicle.getTankToTankConnectionsListBoxStr();

        if(~isempty(t2TConns))
            selT2TConn = get(handles.connsListBox,'Value');
            t2TConn = lv.getTankToTankForInd(selT2TConn);

            lvd_EditTankToTankConnGUI(t2TConn);

            lvState = lvdData.initStateModel.lvState;
            connState = lvState.getTank2TankConnStateForConn(t2TConn);
            connState.flowRate = t2TConn.initFlowRate;

            set(handles.connsListBox,'String',lvdData.launchVehicle.getTankToTankConnectionsListBoxStr());
        end
    end

% --- Executes during object creation, after setting all properties.
function connsListBox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to connsListBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in addConnButton.
function addConnButton_Callback(hObject, eventdata, handles)
% hObject    handle to addConnButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditTankToTanksConnsGUI,'lvdData');

    tempSrcTank = lvdData.launchVehicle.getTankForInd(1);
    tempTgtTank = LaunchVehicleTank.empty(1,0);

    t2TConn = TankToTankConnection(tempSrcTank, tempTgtTank);
    useT2TConn = lvd_EditTankToTankConnGUI(t2TConn);

    if(useT2TConn)
        lvdData.launchVehicle.addTankToTankConnection(t2TConn);

        lvState = lvdData.initStateModel.lvState;
        connState = TankToTankConnState(t2TConn);
        connState.flowRate = t2TConn.initFlowRate;
        lvState.addT2TConnState(connState);

        set(handles.connsListBox,'String',lvdData.launchVehicle.getTankToTankConnectionsListBoxStr());
    end

% --- Executes on button press in removeConnButton.
function removeConnButton_Callback(hObject, eventdata, handles)
% hObject    handle to removeConnButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditTankToTanksConnsGUI,'lvdData');
    lv = lvdData.launchVehicle;

    selConn = get(handles.connsListBox,'Value');
    conn = lv.getTankToTankForInd(selConn);

    tf = conn.isInUse();

    if(tf == false)
        lvState = lvdData.initStateModel.lvState;
        lvState.removeT2TConnStateForConn(conn);

        lv.removeTankToTankConnection(conn);

        listBoxStr = lvdData.launchVehicle.getTankToTankConnectionsListBoxStr();
        set(handles.connsListBox,'String',listBoxStr);

        numConns = length(listBoxStr);
        if(numConns == 0)
            handles.connsListBox.String = ' ';
            set(handles.connsListBox,'Value',1);
        elseif(selConn > numConns)
            set(handles.connsListBox,'Value',numConns);
        end
    else
        warndlg(sprintf('Could not delete the "%s" connection because it is in use as part of an event termination condition, event action, objective function, or constraint.  Remove the dependencies before attempting to delete the connection.', conn.getName()),'Cannot Delete Connection','modal');
    end

% --- Executes on key press with focus on lvd_EditTankToTankConnsGUI or any of its controls.
function lvd_EditTankToTanksConnsGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditTankToTankConnsGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            uiresume(handles.lvd_EditTankToTanksConnsGUI);
        case 'enter'
            uiresume(handles.lvd_EditTankToTanksConnsGUI);
        case 'escape'
            uiresume(handles.lvd_EditTankToTanksConnsGUI);
    end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%--- help for sym/collect ---
%
% COLLECT Collect coefficients.
%    COLLECT(S,v) regards each element of the symbolic matrix S as a
%    polynomial in v and rewrites S in terms of the powers of v.
%    COLLECT(S) uses the default variable determined by SYMVAR.
%    COLLECT(S, 'f') uses all symbolic calls to the function f as variables.
%    COLLECT(S, {'f1', ..., 'fk'}) uses all symbolic calls to any of the
%    functions f1, ..., fk as variables.
%
%    Examples:
%       syms x y
%
%       collect(x^2*y + y*x - x^2 - 2*x)  returns (y - 1)*x^2 + (y - 2)*x
%
%       f = -1/4*x*exp(-2*x)+3/16*exp(-2*x)
%       collect(f,exp(-2*x))  returns -(x/4 - 3/16)/exp(2*x)
%
%       f = x*sin(x) + sin(x) + x*sin(2*x) - sin(2*x)
%       collect(f, 'sin')  returns (x + 1)*sin(x) + (x - 1)*sin(2*x)
%
%    See also SYM/SIMPLIFY, SYM/FACTOR, SYM/EXPAND, SYM/SYMVAR.
%
%    Other functions named collect
%
%       ts/collect
%

author: jmaih


text:
function [xu,yu,xl,yl,xc,yc]=GenerateNACASeries4Airfoil(afid,dotcount,xmod)

%--------------------------------------------------------------------------
%GenerateNACASeries4Airfoil
%Version 1.20
%Created by Stepen (zerocross_raptor@yahoo.com)
%Created 20 November 2010
%Last modified 30 November 2011
%--------------------------------------------------------------------------
%GenerateNACASeries4Airfoil generates the airfoil vertexes' coordinate of
%the given NACA Series 4 Airfoil. The math equation used to generates the
%airfoil coordinates is based on Theory of Wing Section Chapter 6 by Abbott
%and Doenhoff.
%--------------------------------------------------------------------------
%Syntax:
%[xu,yu,xl,yl,xc,yc]=GenerateNACASeries4Airfoil(afid,dotcount,xmod)
%Input argument:
%- afid (1 x 4 str) specifies NACA Series 4 Airfoil identifier.
%- dotcount (1 x 1 int) specifies the number of vertexes to be generated on
%  the airfoil's camber/mean line.
%- xmod (str) specifies the mode of airfoil vertex distribution. Enter
%  'Uniform' to create uniform vertex distribution or 'Cosine' to create
%  vertex distribution based on cosine function (More vertex on the leading
%  edge region).
%Output argument:
%- xu (i x 1 num) specifies the x axis location of airfoil's upper surface
%  vertexes in fraction of chord. The airfoil's upper surface vertex are
%  arranged from leading edge (the first element of xu) to the trailing
%  edge (the last element of xu).
%- yu (i x 1 num) specifies the y axis location of airfoil's upper surface
%  vertexes in fraction of chord. The airfoil's upper surface vertex are
%  arranged from leading edge (the first element of yu) to the trailing
%  edge (the last element of yu).
%- xl (i x 1 num) specifies the x axis location of airfoil's lower surface
%  vertexes in fraction of chord. The airfoil's lower surface vertex are
%  arranged from leading edge (the first element of xl) to the trailing
%  edge (the last element of xl).
%- yl (i x 1 num) specifies the y axis location of airfoil's lower surface
%  vertexes in fraction of chord. The airfoil's lower surface vertex are
%  arranged from leading edge (the first element of yl) to the trailing
%  edge (the last element of yl).
%- xc (i x 1 num) specifies the x axis location of airfoil's camber line
%  vertexes in fraction of chord. The airfoil's camber line vertex are
%  arranged from leading edge (the first element of xc) to the trailing
%  edge (the last element of xc).
%- yc (i x 1 num) specifies the y axis location of airfoil's camber line
%  vertexes in fraction of chord. The airfoil's camber line vertex are
%  arranged from leading edge (the first element of yc) to the trailing
%  edge (the last element of yc).
%--------------------------------------------------------------------------

%CodeStart-----------------------------------------------------------------
%Checking input afid
    if ~ischar(afid)
        error('Airfoil identifier must be a string!')
    end
    if numel(afid)~=4
        error('Airfoil identifier must be a 4 digit number!')
    end
    if isempty(str2double(afid))
        error('Airfoil identifier must be a 4 digit number!')
    end
%Checking input dotcount
    if numel(dotcount)~=1
        error('Number of vertex must be scalar!')
    end
    if (mod(dotcount,1~=0))||(dotcount<0)
        error('Number of vertex must be positive integer!')
    end
%Checking input xmod
    if nargin<3
        xmod='Cosine';
    end
    if (~strcmpi(xmod,'Uniform'))&&(~strcmpi(xmod,'Cosine'))
        error('Vertex distribution input must be Uniform or Cosine!')
    end
%Assigning identifier to equation coefficient
    id1=str2double(afid(1));
    id2=str2double(afid(2));
    id3=str2double(afid([3,4]));
    m=id1*(1/100);            %Maximum camber
    p=id2*(10/100);           %Maximum camber location
    t=id3*(1/100);            %Maximum thickness
    if t==0
        warning(['Zero thickness airfoil!',...
                 ' Airfoil will be just a camber line!'])
    end
%Calculating x-axis location of camber line vertexes
    panelcount=dotcount-1;
    if strcmpi(xmod,'Uniform')
        panellength=1/panelcount;
        xc=(0:panellength:1)';
    elseif strcmpi(xmod,'Cosine')
        deltadeg=90/panelcount;
        xc=1-cosd(0:deltadeg:90)';
    end
%Preallocating array for speed
    yc=zeros(dotcount,1);
    gc=zeros(dotcount,1);
    yt=zeros(dotcount,1);
    xu=zeros(dotcount,1);
    xl=zeros(dotcount,1);
    yu=zeros(dotcount,1);
    yl=zeros(dotcount,1);
%Calculating y-axis location of camber line vertexes and camber gradient
    if m~=0
        for i=1:1:dotcount
            if xc(i)<=p
                yc(i)=(m/(p^2))*((2*p*xc(i))-(xc(i)^2));
            elseif xc(i)>p
                yc(i)=(m/((1-p)^2))*(1-(2*p)+(2*p*xc(i))-(xc(i)^2));
            end
            gc(i)=(m/p^2)*(2*p+2*xc(i));
        end
    end
%Converting camber gradient to camber slope and normal
    sc=atand(gc);
%Calculating thickness distribution
    for i=1:1:dotcount
        yt(i)=5*t*((0.29690*(xc(i)^0.5))-...
                   (0.12600*xc(i))-...
                   (0.35160*(xc(i)^2))+...
                   (0.28430*(xc(i)^3))-...
                   (0.10150*(xc(i)^4)));
    end
%Generating airfoil vertexes
    for i=1:1:dotcount
        xu(i)=xc(i)-yt(i)*sind(sc(i));
        yu(i)=yc(i)+yt(i)*cosd(sc(i));
        xl(i)=xc(i)+yt(i)*sind(sc(i));
        yl(i)=yc(i)-yt(i)*cosd(sc(i));
    end
%CodeEnd-------------------------------------------------------------------

end


text:
% The COBRAToolbox: testGenerateRules.m
%
% Purpose:
%     - testGenerateRules tests generateRules
%
% Authors:
%     - Initial Version: Uri David Akavia - August 2017

global CBTDIR

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testGenerateRules'));
cd(fileDir);

if 1
    modelsToTry = {'Acidaminococcus_intestini_RyC_MR95.mat', 'Acidaminococcus_sp_D21.mat', 'Recon1.0model.mat', 'ecoli_core_model.mat', 'modelReg.mat'};
else
    modelsToTry = {'Recon2.v05.mat'};% TODO fix why this is not working
end
for i=1:length(modelsToTry)
    model = getDistributedModel(modelsToTry{i});
    fprintf('Beginning model %s\n', modelsToTry{i});

    model2 = generateRules(model);
    model.rules = strrep(model.rules, '  ', ' ');
    model3 = model;
    model3.genes = []; % Empty genes field
    model3 = generateRules(model3, 0);
    model4 = model;
    model4 = rmfield(model4, 'genes');
    model4 = generateRules(model4, 0);

    fp = FormulaParser();
    % fix for Recon2
    if strcmp(modelsToTry{i}, 'Recon2.v05.mat')
        model.rules(2240) = {'(x(2)) | (x(4)) | (x(3))'}; % '(26.1) or (314.1) or (314.2)'
        model.rules(2543) = {'(x(2)) | (x(4)) | (x(1)) | (x(3))'}; % '(26.1) or (314.1) or (8639.1) or (314.2)'
        model.rules(2750) = {'(x(2)) | (x(4)) | (x(1)) | (x(3))'}; % '(26.1) or (314.1) or (8639.1) or (314.2)'
        model.rules(2940) = {'(x(4)) | (x(60)) | (x(2)) | (x(61)) | (x(3))'}; % '(314.1) or (4128.1) or (26.1) or (4129.1) or (314.2)'
        model.rules(3133) = {'(x(2)) | (x(4)) | (x(1)) | (x(3))'}; % '(26.1) or (314.1) or (314.2)'
    end
    %also, introduce spaces around the individual genes:
    model.rules = regexprep(model.rules,'([^\s])(x\([0-9]+\))','$1 $2');
    model.rules = regexprep(model.rules,'(x\([0-9]+\))([^\s])','$1 $2');

    % Create a model with sorted genes
    sortedModel = model;
    sortedRules = sortedModel.rules;
    [sortedModel.genes, newInd] = sort(model.genes);
    [~, newInd] = sort(newInd);
    for j=1:length(newInd)
        sortedRules = strrep(sortedRules, ['x(', num2str(j), ')'], ['x(New', num2str(newInd(j)), ')']);
    end
    sortedRules = strrep(sortedRules, 'New', '');
    sortedModel.rules = sortedRules;

    for rule = 1:numel(model.rules)
        if isempty(model.rules{rule})
            %assert that both formulas are empty (and thus equal).
            assert(isequal(model.rules{rule},model2.rules{rule}))
            assert(isequal(sortedModel.rules{rule},model3.rules{rule}))
            assert(isequal(sortedModel.rules{rule},model4.rules{rule}))
            continue;
        elseif strcmp(model.rules{rule},model2.rules{rule})
            %If the strings are identical than we don't need to check
            %equivalence.
            assert(isequal(model.rules{rule},model2.rules{rule}))
            assert(isequal(sortedModel.rules{rule},model3.rules{rule}))
            assert(isequal(sortedModel.rules{rule},model4.rules{rule}))
            continue
        end

        head1 = fp.parseFormula(model.rules{rule});
        head2 = fp.parseFormula(model2.rules{rule});
        headSort = fp.parseFormula(sortedModel.rules{rule});
        head3 = fp.parseFormula(model3.rules{rule});
        head4 = fp.parseFormula(model4.rules{rule});

        %Assert that the formulas are equal.
        assert(head1.isequal(head2));
        assert(headSort.isequal(head3));
        assert(headSort.isequal(head4));
    end
    fprintf('Succesfully completed model %s\n', modelsToTry{i});
end

author: opencobra

repo: cobratoolbox


text:
function err = mne_write_w_file1(filename, w)
% mne_write_w_file1(filename, w)
%
%  writes a binary 'w' file
%
%  filename - name of file to write to
%  w        - a structure with the following fields:
%
% vertices - vector of vertex indices (1-based)
% data     - vector of data values
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
% Revision 1.5  2006/04/23 15:29:41  msh
% Added MGH to the copyright
%
% Revision 1.4  2006/04/10 23:26:54  msh
% Added fiff reading routines
%
% Revision 1.3  2005/12/05 20:23:21  msh
% Added fiff_save_evoked. Improved error handling.
%
% Revision 1.2  2005/11/21 03:19:12  msh
% Improved error handling
%
% Revision 1.1  2005/11/21 02:15:51  msh
% Added more routines
%

me='MNE:mne_write_w_file1';
if(nargin ~= 2)
    error(me,'usage: mne_write_w_file1(filename, w)');
    return;
end

try
    w.vertices = w.vertices - 1;
    mne_write_w_file(filename,w);
catch
    w.vertices = w.vertices - 1;
    error(me,'%s',mne_omit_first_line(lasterr));
end

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip