text:
function [result,croppedImage] = isGoodPhoto(img)

minWidth = 480;
minHeight = 480;

%minWidth = 128;
%minHeight = 128;

croppedImage = [];

if isa(img,'char')
    img = imread(img);
end

if isa(img,'logical')
    result = 'binary';
    return;
end

if ndims(img)>3
    result = 'animation';
    return;
end

if size(img,1)<minHeight || size(img,2)<minWidth
    result = 'too small';
    return;
end

persistent emptyFlickrGIF;
persistent emptyFlickrPNG;
if isempty(emptyFlickrGIF)
    [emptyFlickrGIF,map]= imread('emptyFlickr.gif');
    if ~isempty( map )
        emptyFlickrGIF = ind2rgb( emptyFlickrGIF, map );
    end
end
if isempty(emptyFlickrPNG)
    emptyFlickrPNG = imread('emptyFlickr.png');
end

% check if it is flickr empty image
imgT = imresize(img,[size(emptyFlickrGIF,1), size(emptyFlickrGIF,2)]);
if size(imgT,3)==1
    diff= sum(reshape(abs(im2double(imgT)-rgb2gray(emptyFlickrGIF)),1,[]));
else
    diff= sum(reshape(abs(im2double(imgT)-emptyFlickrGIF),1,[]));
end
if diff<50
    result = 'flickr empty';
    return;
end

imgT = imresize(img,[size(emptyFlickrPNG,1), size(emptyFlickrPNG,2)]);
if size(imgT,3)==1
    diff= sum(reshape(abs(imgT-rgb2gray(emptyFlickrPNG)),1,[]));
else
    diff= sum(reshape(abs(imgT-emptyFlickrPNG),1,[]));
end
if diff<50
    result = 'flickr empty';
    return;
end

% check if it is line drawing
if size(img,3)==1
    h = imhist(img);
    if sum(h(240:end))/sum(h) > 0.5
        result = 'too white';
        return;
    elseif sum(h(1:10))/sum(h) > 0.5
        result = 'too black';
        return;
    end
else
    h = imhist(img(:,:,1)) + imhist(img(:,:,2)) + imhist(img(:,:,3));
    if sum(h(240:end))/sum(h) > 0.6
        result = 'too white';
        return;
    elseif sum(h(1:10))/sum(h) > 0.8
        result = 'too black';
        return;
    end
end

imgG = img;

G = fspecial('gaussian',[5 5],2);
%# Filter it
imgG = imfilter(imgG,G,'symmetric','same');

if size(imgG,3)>1
    imgG = rgb2gray(imgG);
end

tto = edge(imgG,'canny');

if sum(tto(:))/numel(tto) < 0.01
    result = 'too pure';
    return;
end

%{
% check if it is line drawing
imgG = img;
if size(imgG,3)>1
    imgG = rgb2gray(imgG);
end
h = imhist(imgG);
if sum(h(240:end))/sum(h) > 0.5
    result = 'too white';
    return;
elseif sum(h(1:10))/sum(h) > 0.5
    result = 'too black';
    return;
else
    [~,ind]=max(h);
    if sum(h(max(1,ind-7):min(256,ind+7)))/sum(h) > 0.7
        result = 'too pure';
        return;
    end
end
%}

% maybe useful http://www.mathworks.com/matlabcentral/fileexchange/25354-cropmat

if size(img,3)==1
    ttWhite = img>252;
else
    ttWhite = img(:,:,1)>252 & img(:,:,2)>252 & img(:,:,3)>252;
end

if size(img,3)==1
    ttBlack = img<5;
else
    ttBlack = img(:,:,1)<5 & img(:,:,2)<5 & img(:,:,3)<5;
end

wWhite=sum(ttWhite,1) > size(tto,1)*0.7;
hWhite=sum(ttWhite,2) > size(tto,2)*0.7;
wBlack=sum(ttBlack,1) > size(tto,1)*0.7;
hBlack=sum(ttBlack,2) > size(tto,2)*0.7;

anyEdgeW = sum(tto,1) > size(tto,1)*0.01;
anyEdgeH = sum(tto,2) > size(tto,2)*0.01;

minW=find(anyEdgeW, 1 )-1+2;   maxW=size(imgG,2)-find(anyEdgeW, 1, 'last' );
minH=find(anyEdgeH, 1 )-1+2;   maxH=size(imgG,1)-find(anyEdgeH, 1, 'last' );

try
    if minW<3 || maxW<3
        minW = 0;
        maxW = 0;
    else
        wMargin = wWhite | wBlack;
        if ~(any(wMargin(1:minW)) && any(wMargin(end-maxW:end)))
            minW = 0;
            maxW = 0;
        end
    end

    if minH<3 || maxH<3
        minH = 0;
        maxH = 0;
    else
        hMargin = hWhite | hBlack;
        if ~(any(hMargin(1:minH)) && any(hMargin(end-maxH:end)))
            minH = 0;
            maxH = 0;
        end
    end

    if minW>0 || maxW>0 || minH>0 || maxH>0
        croppedImage = img(minH+1:size(imgG,1)-maxH,minW+1:size(imgG,2)-maxW,:);
        result = 'crop';

        if size(croppedImage,1)<minHeight || size(croppedImage,2)<minWidth
            result = 'too small';
        end

        return;
    end
catch
end

%{

% automatic image cropping
% http://stackoverflow.com/questions/11121657/find-the-edges-of-image-and-crop-it-in-matlab
%# instead of "==" you can check for similarity within a tolerance
%tt=img(:,:,1)==img(:,:,2) & img(:,:,2) == img(:,:,3);

if size(img,3)==1
    tt = img>252;
else
    tt = img(:,:,1)>252 & img(:,:,2)>252 & img(:,:,3)>252;
end

%# invert tt so that it's 1 where there is signal
tt = ~tt;

%# clean up some of the smaller artifacts
tto = imopen(tt,strel('square',10));

%# get the areas and bounding box of the areas above threshold
%# as an additional criterion, you could also use excentricity

stats = regionprops(tto,'BoundingBox','Area');
if ~isempty(stats)
    area = cat(1,stats.Area);
    [~,maxAreaIdx] = max(area);
    bb = round(stats(maxAreaIdx).BoundingBox);

    %# note that regionprops switches x and y (it's a long story)
    croppedImage = img(bb(2):bb(2)+bb(4)-1,bb(1):bb(1)+bb(3)-1,:);

    if size(croppedImage,1)<size(img,1)-5 || size(croppedImage,2)<size(img,2)-5
        result = 'crop';
        return;
    end
end

%% black
if size(img,3)==1
    tt = img<5;
else
    tt = img(:,:,1)<5 & img(:,:,2)<5 & img(:,:,3)<5;
end

%# invert tt so that it's 1 where there is signal
tt = ~tt;

%# clean up some of the smaller artifacts
tto = imopen(tt,strel('square',10));

%# get the areas and bounding box of the areas above threshold
%# as an additional criterion, you could also use excentricity

stats = regionprops(tto,'BoundingBox','Area');
if ~isempty(stats)
    area = cat(1,stats.Area);
    [~,maxAreaIdx] = max(area);
    bb = round(stats(maxAreaIdx).BoundingBox);

    %# note that regionprops switches x and y (it's a long story)
    croppedImage = img(bb(2):bb(2)+bb(4)-1,bb(1):bb(1)+bb(3)-1,:);

    if size(croppedImage,1)<size(img,1)-5 || size(croppedImage,2)<size(img,2)-5
        result = 'crop';
        return;
    end
end

%}

result = 'good';
return;

author: jianxiongxiao

repo: ProfXkit


text:
function varargout = process_pac_dynamic_sur2( varargin )
% PROCESS_PAC_DYNAMIC: Compute the Time resolved Phase-Amplitude Coupling
%
% DOCUMENTATION

% @=============================================================================
% 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: Soheila Samiee, Francois Tadel 2013
%
% Updates:
%   - 1.0.4:  Soheila
%   - 2.0.0: Soheila, based on dpac v.1.1.3 and surrogate 3, June 2016
%   - 2.1:   SS, Based on dpac v2.1.0, Nov. 2016
%   - 3.1:   SS, block resampling (N = 2)
%   - 3.3:   SS, number of blocks for resampling could be more than 2
%            (default =5 ), May 2017
eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Surrogate tPAC - block resampling - version: 3.4';
    sProcess.FileTag     = '';
%     sProcess.Category    = 'Custom';
%     sProcess.SubGroup    = 'dPAC analysis';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Frequency','Time-resolved Phase-Amplitude Coupling'};
    sProcess.Index       = 660;
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw',      'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;

    % === TIME WINDOW
    sProcess.options.timewindow.Comment = 'Time:';
    sProcess.options.timewindow.Type    = 'timewindow';
    sProcess.options.timewindow.Value   = [];
    % === NESTING FREQ
    sProcess.options.nesting.Comment = 'Frequency for phase band (low):';
    sProcess.options.nesting.Type    = 'range';
    sProcess.options.nesting.Value   = {[8, 12], 'Hz', 2};
    % === NESTED FREQ
    sProcess.options.nested.Comment = 'Frequency for amplitude band (high):';
    sProcess.options.nested.Type    = 'range';
    sProcess.options.nested.Value   = {[40, 150], 'Hz', 2};
        % Band for fa
%     sProcess.options.label_fa.Comment = 'F_A frequency band:';
%     sProcess.options.label_fa.Type    = 'label';
    sProcess.options.fa_type.Comment = {'   Single band', ...
                                        '   More than one center frequencies (default: 20)' };
    sProcess.options.fa_type.Type    = 'radio';
    sProcess.options.fa_type.Value   = 1;

    % === WINDOW LENGTH
    sProcess.options.winLen.Comment = 'Length of sliding time window:';
    sProcess.options.winLen.Type    = 'value';
    sProcess.options.winLen.Value   = {1.10, 'S', 2};

    % === NUMBER OF SURROGATE DATASETS
    sProcess.options.Nsurrogate.Comment    = 'Number of surrogate datasets: ';
    sProcess.options.Nsurrogate.Type       = 'value';
    sProcess.options.Nsurrogate.Value      = {100, '', 0};

    % === SOURCES
    sProcess.options.label5.Comment = 'Sensors/sources to be investigated :';
    sProcess.options.label5.Type    = 'label';

    % === CLUSTERS
    sProcess.options.clusters.Comment = '';
    sProcess.options.clusters.Type    = 'scout_confirm';
    sProcess.options.clusters.Value   = {};
    sProcess.options.clusters.InputTypes = {'results'};

    % === SENSOR SELECTION
    sProcess.options.target_data.Comment    = 'Sensor types or names (empty=all): ';
    sProcess.options.target_data.Type       = 'text';
    sProcess.options.target_data.Value      = 'MEG, EEG';
    sProcess.options.target_data.InputTypes = {'data', 'raw'};
    % === SOURCE INDICES
    sProcess.options.target_res.Comment    = 'Source indices (empty=all): ';
    sProcess.options.target_res.Type       = 'text';
    sProcess.options.target_res.Value      = '';
    sProcess.options.target_res.InputTypes = {'results'};
    sProcess.options.label6.Comment = '(The indices will only be considered if the scouts are not selected)';
    sProcess.options.label6.Type    = 'label';

    % === ROW NAMES
    sProcess.options.target_tf.Comment    = 'Row names or indices (empty=all): ';
    sProcess.options.target_tf.Type       = 'text';
    sProcess.options.target_tf.Value      = '';
    sProcess.options.target_tf.InputTypes = {'timefreq', 'matrix'};

    % === LOOP METHOD
    sProcess.options.label1.Comment = 'Processing options [expert only]:';
    sProcess.options.label1.Type    = 'label';
    % === MAX_BLOCK_SIZE
    sProcess.options.max_block_size.Comment = 'Number of signals to process at once: ';
    sProcess.options.max_block_size.Type    = 'value';
    sProcess.options.max_block_size.Value   = {20, ' ', 0};

    % sProcess.options.filter_sensor.InputTypes = {'results'};
    % === AVERAGE OUTPUT FILES
    sProcess.options.label2.Comment = 'Output options:';
    sProcess.options.label2.Type    = 'label';
    sProcess.options.avgoutput.Comment = 'Save average PAC across trials';
    sProcess.options.avgoutput.Type    = 'checkbox';
    sProcess.options.avgoutput.Value   = 1;
%     % === SAVE PAC MAPS
%     sProcess.options.savefull.Comment = 'Save the full PAC maps';
%     sProcess.options.savefull.Type    = 'checkbox';
%     sProcess.options.savefull.Value   = 1;
end

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

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputsA) %#ok<DEFNU>
tic
    % Get options
    if isfield(sProcess.options, 'timewindow') && isfield(sProcess.options.timewindow, 'Value') && iscell(sProcess.options.timewindow.Value) && ~isempty(sProcess.options.timewindow.Value)
        OPTIONS.TimeWindow = sProcess.options.timewindow.Value{1};
    else
        OPTIONS.TimeWindow = [];
    end
    OPTIONS.FunctionVersion = sProcess.Comment;

    % Clusters
    if isfield(sProcess.options, 'clusters') && ~isempty(sProcess.options.clusters) && ~isempty(sProcess.options.clusters.Value)
        OPTIONS.Clusters = sProcess.options.clusters.Value;
    else
        OPTIONS.Clusters = [];
    end
    OPTIONS.BandNesting = sProcess.options.nesting.Value{1};
    OPTIONS.BandNested  = sProcess.options.nested.Value{1};
    OPTIONS.WinLen      = sProcess.options.winLen.Value{1};
    % Get target
    if  ~isempty(OPTIONS.Clusters)         % extract cluster
        OPTIONS.Target = OPTIONS.Clusters;
    elseif ismember(sInputsA(1).FileType, {'data','raw'}) && isfield(sProcess.options, 'target_data') && ~isempty(sProcess.options.target_data.Value)
        OPTIONS.Target = sProcess.options.target_data.Value;
    elseif strcmpi(sInputsA(1).FileType, 'results') && isfield(sProcess.options, 'target_res') && ~isempty(sProcess.options.target_res.Value)
        OPTIONS.Target = sProcess.options.target_res.Value;
    elseif ismember(sInputsA(1).FileType, {'timefreq', 'matrix'}) && isfield(sProcess.options, 'target_tf') && ~isempty(sProcess.options.target_tf.Value)
        OPTIONS.Target = sProcess.options.target_tf.Value;
    else
        OPTIONS.Target = [];
    end
    % All other options
    OPTIONS.MaxSignals   = sProcess.options.max_block_size.Value{1};
    if (strcmp(sInputsA(1).FileType,'data') && isempty(sProcess.options.target_data.Value)) || ...
            (strcmp(sInputsA(1).FileType,'results') && isempty(sProcess.options.target_res.Value))
        OPTIONS.isFullMaps   = 1;
    else
        OPTIONS.isFullMaps   = 0;
    end
    OPTIONS.isAvgOutput  = sProcess.options.avgoutput.Value;
    if (length(sInputsA) == 1)
        OPTIONS.isAvgOutput = 0;
    end
    OPTIONS.HighFreqs    = sProcess.options.fa_type.Value;
    OPTIONS.Nsur = sProcess.options.Nsurrogate.Value{1};

    % ===== INITIALIZE =====
    % Initialize output variables
    OutputFiles = {};
    sPAC_avg = [];
    nAvg = 0;
    % Initialize progress bar
    if bst_progress('isVisible')
        startValue = bst_progress('get');
    else
        startValue = 0;
    end
    % Options for LoadInputFile()
    if strcmpi(sInputsA(1).FileType, 'results')
        LoadOptions.LoadFull = 0;  % Load kernel-based results as kernel+data
    else
        LoadOptions.LoadFull = 1;  % Load the full file
    end
    LoadOptions.IgnoreBad   = 1;  % From raw files: ignore the bad segments
    LoadOptions.ProcessName = func2str(sProcess.Function);
    LoadOptions.TargetFunc = 'all';
    % Start the matlabpool for parallel processing in bst_pac

    % Loop over input files
    for iFile = 1:length(sInputsA)

        % ===== LOAD SIGNALS =====
        bst_progress('text', sprintf('PAC: Loading input file (%d/%d)...', iFile, length(sInputsA)));
        bst_progress('set', round(startValue + (iFile-1) / length(sInputsA) * 100));

        % Load input signals
        [sInput, nSignals, iRows] = bst_process('LoadInputFile', sInputsA(iFile).FileName, OPTIONS.Target, OPTIONS.TimeWindow, LoadOptions);
        if isempty(sInput) || isempty(sInput.Data)
            return;
        end

        % Get sampling frequency
        sRate = 1 / (sInput.Time(2) - sInput.Time(1));
        % Check the nested frequencies
        if (OPTIONS.BandNested(2) > sRate/3)
            % Warning
            strMsg = sprintf('Higher nesting frequency is too high (%d Hz) compared with sampling frequency (%d Hz): Limiting to %d Hz', round(OPTIONS.BandNested(2)), round(sRate), round(sRate/3));
            disp([10 'process_pac> ' strMsg]);
            bst_report('Warning', sProcess, [], strMsg);
            % Fix higher frequencyy
            OPTIONS.BandNested(2) = sRate/3;
        end
        % Check the extent of bandNested band
        if (OPTIONS.BandNested(2) <= OPTIONS.BandNested(1))
            bst_report('Error', sProcess, [], sprintf('Invalid frequency range: %d-%d Hz', round(OPTIONS.BandNested(1)), round(OPTIONS.BandNested(2))));
            continue;
        end

        % ===== COMPUTE PAC MEASURE =====
        % Number of blocks of signals
        MAX_BLOCK_SIZE = OPTIONS.MaxSignals;
        nBlocks = ceil(nSignals / MAX_BLOCK_SIZE);
        sPAC = [];
        % Display processing time
        disp(sprintf('Processing %d blocks of %d signals each.', nBlocks, MAX_BLOCK_SIZE));
        % Process each block of signals
        for iBlock = 1:nBlocks
%             tic
            bst_progress('text', sprintf('PAC: File %d/%d - Block %d/%d', iFile, length(sInputsA), iBlock, nBlocks));
            bst_progress('set', round(startValue + (iFile-1)/length(sInputsA)*100 + iBlock/nBlocks*100));

            % Indices of the signals
            iSignals = (iBlock-1)*MAX_BLOCK_SIZE+1 : min(iBlock*MAX_BLOCK_SIZE, nSignals);

            % Get signals to process
            if ~isempty(sInput.ImagingKernel)
                Fblock = sInput.ImagingKernel(iSignals,:) * sInput.Data;
            else
                Fblock = sInput.Data(iSignals,:);
            end

            % Defining the options
            PACoptions.doInterpolation = 1;%0               % Applying interpolation in frequency and time domain
            PACoptions.logCenters = 0; %1                   % Choose the center frequencies for f_A with log space in faBand
            PACoptions.overlap = 0.5;                       % Time window over lap (0<= value <1)
            if OPTIONS.HighFreqs ==1
                PACoptions.nHighFreqs = 1;                  % Number of high frequency centers
                PACoptions.doInterpolation = 0;
            else
                PACoptions.nHighFreqs = 20; %4               % Number of high frequency centers
            end
            PACoptions.nSur = OPTIONS.Nsur;
            OPTIONS.PACoptions = PACoptions;

            % Estimating DPAC
            sPACblock = Compute(Fblock, sRate,  OPTIONS.BandNested, OPTIONS.BandNesting, OPTIONS.WinLen, PACoptions);

            % Check for errors
            if isempty(sPACblock)
                return;
            end
            % Initialize output structure
            nTime = length(sPACblock.TimeOut);
            if isempty(sPAC)
                sPAC.ValPAC      = [];
                sPAC.NestingFreq = [];
                sPAC.NestedFreq  = [];
                sPAC.PhasePAC    = [];
                sPAC.DynamicPAC     = zeros(nSignals, nTime, length(sPACblock.HighFreqs), OPTIONS.Nsur);
                sPAC.DynamicNesting = zeros(nSignals, nTime, length(sPACblock.HighFreqs), OPTIONS.Nsur);
                sPAC.DynamicPhase   = zeros(nSignals, nTime, length(sPACblock.HighFreqs), OPTIONS.Nsur);
                sPAC.HighFreqs = sPACblock.HighFreqs;

                meanInputTime  = OPTIONS.TimeWindow(1)+(sPACblock.TimeOut(end)+OPTIONS.WinLen/2)/2;
                meanOutputTime = (sPACblock.TimeOut(1)+sPACblock.TimeOut(end))/2;
                sPAC.TimeOut   = sPACblock.TimeOut + (meanInputTime - meanOutputTime);
            end
            % Copy block results to output structure
            sPAC.DynamicPAC(iSignals,:,:,:)     = permute(sPACblock.DynamicPAC,[3,2,1,4]);
            sPAC.DynamicNesting(iSignals,:,:,:) = permute(sPACblock.DynamicNesting,[3,2,1,4]);
            sPAC.DynamicPhase(iSignals,:,:,:)   = permute(sPACblock.DynamicPhase,[3,2,1,4]);
        end

        % ===== APPLY SOURCE ORIENTATION =====
        if strcmpi(sInput.DataType, 'results') && (sInput.nComponents > 1)
            % Number of values per vertex
            switch (sInput.nComponents)
                case 2
                    sPAC.ValPAC         = (sPAC.ValPAC(1:2:end,:,:)           + sPAC.ValPAC(2:2:end,:,:))           / 2;
                    sPAC.NestingFreq    = (sPAC.NestingFreq(1:2:end,:,:)      + sPAC.NestingFreq(2:2:end,:,:))      / 2;
                    sPAC.NestedFreq     = (sPAC.NestedFreq(1:2:end,:,:)       + sPAC.NestedFreq(2:2:end,:,:))       / 2;
                    sPAC.PhasePAC       = (sPAC.PhasePAC(1:2:end,:,:)         + sPAC.PhasePAC(2:2:end,:,:))         / 2;
                    sPAC.DynamicPAC     = (sPAC.DynamicPAC(1:2:end,:,:,:)     + sPAC.DynamicPAC(2:2:end,:,:,:))     / 2;
                    sPAC.DynamicNesting = (sPAC.DynamicNesting(1:2:end,:,:,:) + sPAC.DynamicNesting(2:2:end,:,:,:)) / 2;
                    sInput.RowNames     = sInput.RowNames(1:2:end);
                    sPAC.DynamicPhase   = (sPAC.DynamicPhase(1:2:end,:,:,:)   + sPAC.DynamicPhase(2:2:end,:,:,:))   / 2;

                case 3
                    sPAC.ValPAC         = (sPAC.ValPAC(1:3:end,:,:)           + sPAC.ValPAC(2:3:end,:,:)           + sPAC.ValPAC(3:3:end,:,:))           / 3;
                    sPAC.NestingFreq    = (sPAC.NestingFreq(1:3:end,:,:)      + sPAC.NestingFreq(2:3:end,:,:)      + sPAC.NestingFreq(3:3:end,:,:))      / 3;
                    sPAC.NestedFreq     = (sPAC.NestedFreq(1:3:end,:,:)       + sPAC.NestedFreq(2:3:end,:,:)       + sPAC.NestedFreq(3:3:end,:,:))       / 3;
                    sPAC.PhasePAC       = (sPAC.PhasePAC(1:3:end,:,:)         + sPAC.PhasePAC(2:3:end,:,:)         + sPAC.PhasePAC(3:3:end,:,:))         / 3;
                    sPAC.DynamicPAC     = (sPAC.DynamicPAC(1:3:end,:,:,:)     + sPAC.DynamicPAC(2:3:end,:,:,:)     + sPAC.DynamicPAC(3:3:end,:,:,:))     / 3;
                    sPAC.DynamicNesting = (sPAC.DynamicNesting(1:3:end,:,:,:) + sPAC.DynamicNesting(2:3:end,:,:,:) + sPAC.DynamicNesting(3:3:end,:,:,:)) / 3;
                    sInput.RowNames     = sInput.RowNames(1:3:end);
                    sPAC.DynamicPhase   = (sPAC.DynamicPhase(1:3:end,:,:,:)   + sPAC.DynamicPhase(2:3:end,:,:,:)   + sPAC.DynamicPhase(3:3:end,:,:,:))   / 3;

            end
        end

        % ===== SAVE FILE =====
        % Detect incomplete lists of sources
        isIncompleteResult = strcmpi(sInput.DataType, 'results') && (length(sInput.RowNames) * sInput.nComponents < nSignals);
        % Comment
        Comment = 'Sur. DynamicPAC';
        if iscell(sInput.RowNames)
            % Find the scout name
            scoutName = sInput.RowNames{1};
            k = strfind(scoutName,'.');
            if isempty(k)
                Comment = [Comment, ': ' scoutName];
            else
                Comment = [Comment, ': ' scoutName(1:k-1)];
            end

        elseif (length(sInput.RowNames) == 1)
            Comment = [Comment, ': #', num2str(sInput.RowNames(1))];
        elseif isIncompleteResult
            Comment = [Comment, ': ', num2str(length(sInput.RowNames)), ' sources'];
        end

        if OPTIONS.isFullMaps
            Comment = [Comment, ' (Full)'];
        end
        % Output data type: if there are not all the sources, switch the datatype to "scout"
        if isIncompleteResult
            sInput.DataType = 'scout';
            % Convert source indices to strings
            if ~iscell(sInput.RowNames)
                sInput.RowNames = cellfun(@num2str, num2cell(sInput.RowNames), 'UniformOutput', 0);
            end
        end
        % Save each as an independent file
        if ~OPTIONS.isAvgOutput
            nAvg = 1;
            OutputFiles{end+1} = SaveFile(sPAC, sInput.iStudy, sInputsA(iFile).FileName, sInput, Comment, nAvg, OPTIONS);
        else
            % Compute online average of the connectivity matrices
            if isempty(sPAC_avg)
                sPAC_avg.ValPAC      = sPAC.ValPAC       ./ length(sInputsA);
                sPAC_avg.NestingFreq = sPAC.NestingFreq  ./ length(sInputsA);
                sPAC_avg.NestedFreq  = sPAC.NestedFreq   ./ length(sInputsA);
                sPAC_avg.PhasePAC(:,:,:,nAvg+1)       =  sPAC.PhasePAC;
                sPAC_avg.DynamicPAC(:,:,:,nAvg+1)     =  sPAC.DynamicPAC;
                sPAC_avg.DynamicNesting(:,:,:,nAvg+1) =  sPAC.DynamicNesting;
                sPAC_avg.DynamicPhase(:,:,:,nAvg+1)   =  sPAC.DynamicPhase;

                sPAC_avg.TimeOut     = sPAC.TimeOut;
                sPAC_avg.HighFreqs   = sPAC.HighFreqs;
            else
                sPAC_avg.ValPAC      = sPAC_avg.ValPAC      + sPAC.ValPAC      ./ length(sInputsA);
                sPAC_avg.NestingFreq = sPAC_avg.NestingFreq + sPAC.NestingFreq ./ length(sInputsA);
                sPAC_avg.NestedFreq  = sPAC_avg.NestedFreq  + sPAC.NestedFreq  ./ length(sInputsA);
                sPAC_avg.PhasePAC(:,:,:,nAvg+1)       =  sPAC.PhasePAC;
                sPAC_avg.DynamicPAC(:,:,:,nAvg+1)     =  sPAC.DynamicPAC;
                sPAC_avg.DynamicPhase(:,:,:,nAvg+1)   =  sPAC.DynamicPhase;
                sPAC_avg.DynamicNesting(:,:,:,nAvg+1) =  sPAC.DynamicNesting;
            end
            nAvg = nAvg + 1;
        end
    end

    % ===== SAVE AVERAGE =====
    if OPTIONS.isAvgOutput
        % Output study, in case of average
        [tmp, iOutputStudy] = bst_process('GetOutputStudy', sProcess, sInputsA);
        % Save file
        OutputFiles{1} = SaveFile(sPAC_avg, iOutputStudy, [], sInput, Comment, nAvg, OPTIONS);
    end
end

%% ========================================================================
%  ===== SUPPORT FUNCTIONS ================================================
%  ========================================================================

%% ===== SAVE FILE =====
function NewFile = SaveFile(sPAC, iOuptutStudy, DataFile, sInput, Comment, nAvg, OPTIONS)
    % ===== PREPARE OUTPUT STRUCTURE =====
    % Create file structure
    FileMat = db_template('timefreqmat');
    FileMat.TF        = sPAC.ValPAC;
    FileMat.Comment   = Comment;
    FileMat.Method    = 'dpac';
    FileMat.Measure   = 'maxpac';
    FileMat.DataFile  = file_win2unix(DataFile);
    FileMat.nAvg      = nAvg;
    FileMat.Freqs     = 0;
    % All the PAC fields
    FileMat.sPAC = rmfield(sPAC, 'ValPAC');
    % Time vector
    FileMat.Time = sPAC.TimeOut;

    % Output data type and Row names
    if isempty(OPTIONS.Target)
        FileMat.DataType = sInput.DataType;
        FileMat.RowNames = sInput.RowNames;
    elseif strcmpi(sInput.DataType, 'results') && ~isempty(OPTIONS.Target)
        FileMat.DataType = 'matrix';
        if isnumeric(sInput.RowNames)
        	FileMat.RowNames = cellfun(@num2str, num2cell(sInput.RowNames), 'UniformOutput', 0);
        else
            FileMat.RowNames = sInput.RowNames;
        end
    else
        FileMat.DataType = sInput.DataType;
        FileMat.RowNames = sInput.RowNames;
    end
    % Atlas
    if ~isempty(sInput.Atlas)
        FileMat.Atlas = sInput.Atlas;
    end
    if ~isempty(sInput.SurfaceFile)
        FileMat.SurfaceFile = sInput.SurfaceFile;
    end
    % History: Computation
    FileMat = bst_history('add', FileMat, 'compute', 'PAC measure (see the field "Options" for input parameters)');
    % Save options in the file
    FileMat.Options = OPTIONS;

    % ===== SAVE FILE =====
    % Get output study
    sOutputStudy = bst_get('Study', iOuptutStudy);
    % File tag
%     if OPTIONS.isFullMaps
        fileTag = 'timefreq_dpac_fullmaps';
%     else
%         fileTag = 'timefreq_dpac';
%     end
    % Output filename
    NewFile = bst_process('GetNewFilename', bst_fileparts(sOutputStudy.FileName), fileTag);
    % Save file
    bst_save(NewFile, FileMat, 'v6');
    % Add file to database structure
    db_add_data(iOuptutStudy, NewFile, FileMat);
    toc
end

% ===== COMPUTE PAC MEASURE =====
function sPAC = Compute(Xinput, sRate, faBand, fpBand, winLen, Options)
% USAGE:
%   sPAC = dPACestimate(Xinput, sRate, faBand, fpBand, winLen, PLvector)
%
% INPUTS:
%    - Xinput:        [nChannels,nTime] signal to process
%    - sRate:         Sampling frequency (Hz)
%    - faBand:        Nested Band: Minimum and maximum frequency for extraction of frequency for amplitude
%    - fpBand:        Nesting Band: Minimum and maximum frequency for extraction of frequency for phase (Hz)
%    - winLen:        Length of each time window for coupling estimation(S) (default: 1 Sec)
%    - isFullMaps:    If 1, save the full directPAC maps
%
% OUTPUTS:   sPAC structure [for each signal]
%    - TimeOut:        Output time vector (Sec)
%    - HighFreqs:       Frequency for amplitude vector
%    - ValPAC:         [nChannels, nTimeOut] Maximum PAC strength in each  time point
%    - NestedFreq:     [nChannels, nTimeOut] Fnested corresponding to maximum synchronization index in each time point
%    - NestingFreq:    [nChannels, nTimeOut] Fnesting corresponding to maximum synchronization index in each time point
%    - phasePAC:       [nChannels, nTimeOut] Phase corresponding to maximum synchronization index in each time point
%    - DynamicNesting: [nNestedCenters,nTimeOut,nChannels] Estimated nesting frequency for all times, channels and nested intervals
%    - DynamicPAC:     [nNestedCenters,nTimeOut,nChannels] full array of PAC
%
% DESCRIPTION:
%   Estimation of Phase Amplitude Coupling (PAC) with DPAC method.
%
% Author:  Soheila Samiee, 2013
%

if (nargin < 4) || isempty(fpBand)
    fpBand = [4, 8];
end
if (nargin < 5) || isempty(winLen)
    winLen = 1;
end
if ~isfield(Options, 'overlap')
    Options.overlap = 0.5;
end

if fpBand(2)>faBand(1)
    fpBand(2) = faBand(1)/2;
    error_msg = ['Maximum of Fp should be less than half of the minimum of Fa!' 10 10 ...
        'max{Fp} modified to ', num2str(fpBand(2))];
    bst_report('Error', sProcess, [], error_msg);
    disp(['Warning: ' error_msg]);
end

if winLen < 2*1/fpBand(1)
    error_msg = ['Window length is short for extracting this minimum fp!' 10 ...
        'Either increase window length to: ',num2str(2*1/fpBand(1)), ' or increase minimum fp to; ', num2str(2/winLen)];
    bst_report('Error', sProcess, [], error_msg);
    disp(['Warning: ' error_msg]);
    winLen = 2*1/fpBand(1);
end

% ===== SETTING THE PARAMETERS =====
tStep = winLen*(1-Options.overlap);     % Time step for sliding window on time (Sec) (Overlap: 50%)
margin = 2;%1                       % Margin (in time) for filtering (Sec) --- default: 2sec -> changed to 1 sec in May12,2016
hilMar = 1/5;                       % Percentage of margin for Hilber transform
nTS = size(Xinput,2);               % Number of temporal samples of data
bandNestingLen= max(3,1/(winLen+margin));% Length of band nesting -- considering the resolution in FFT domain with available window length
isMirror = 0;                       % Mirroring the data in filtering
isRelax  = 1;                       % Attenuation of the filter in the stopband (1 => 40 dB, 0 => 60 dB)
minExtracFreq = max(1/winLen, fpBand(1));  % minimum frequency that could be extracted as nestingFreq
sProcess = 'Process_pac_dynamic_sur2';

% Mode = Options.Mode;                % Mode of estimating the coupling (line or map)
doInterpolation = Options.doInterpolation;  % Applying interpolation in frequency and time domain
logCenters = Options.logCenters;    % Choose the center frequencies for f_A with log space in faBand
nHighFreqs = Options.nHighFreqs;    % Number of high frequency centers
% Mode = 'map';                       % Mode of estimating the coupling (line or map)
% doInterpolation = 1;                % Applying interpolation in frequency and time domain
% logCenters = 0; %0                  % Choose the center frequencies for f_A with log space in faBand
% nHighFreqs = 20;%4;                 % Number of high frequency centers
missedPcount = 0;                   % Number of intervals that do not have peak in their F_A envelope PSD
mirrorEffectSample = 40;            % Number of samples that can be affected due to mirroring effect

N = 5; % Number of blocks for shuffling

% ==== ADDING MARGING TO THE DATA => AVOID EDGE ARTIFACT (FILTERS AND HILBERT TRANSFORM) ====
nMargin = fix(margin*sRate);
nHilMar = fix(nMargin*hilMar);

if (nTS/sRate < 2*winLen)
    error_msg = 'Data length should be at least twice of window length';
    bst_report('Error', 'process_pac_dynamic', [], error_msg);
    disp(['Error: ' error_msg]);
    sPAC = [];
    return
% elseif (nTS < nMargin)
%     tmp = repmat([Xinput(:,end-1:-1:2),Xinput],1,ceil(nMargin/nTS));
%     Xinput = [tmp(:,end-nMargin+1-size(Xinput,2):end-size(Xinput,2)), Xinput, tmp(:,1:nMargin)];
%     clear tmp
% else
%     Xinput = [Xinput(:,nMargin+1:-1:2), Xinput, Xinput(:,end-1:-1:end-nMargin)];
end

% Zero-padding signal for the margin
Xinput = [zeros(size(Xinput,1),nMargin), Xinput, zeros(size(Xinput,1),nMargin)];

% ==== SETTING THE PARAMETERS OF THE FILTERS ====
if nHighFreqs > 1 %strcmp(Mode,'map')
%     Fconst = (faBand(end)-faBand(1))/nHighFreqs;                     % Frequency distant between successive fAs
    if logCenters
        nestedCenters = logspace(log10(faBand(1)),log10(faBand(end)),nHighFreqs);
    else
        nestedCenters = linspace(faBand(1),faBand(end),nHighFreqs);
    end
    Fstep = diff(nestedCenters)/2;  % the range of frequency around each nested center
    Fstep = [Fstep(1),Fstep,Fstep(end)];
    Fstep = max(Fstep, fpBand(2)/2);  % Minimum band width is defined to cover the whole interval between consecutive centre frequencies and at the same time consider all coupled frequencies to it in the range of interest.
else
    nestedCenters = mean(faBand);
    Fstep    = abs(faBand-nestedCenters);
end

fArolloff = [];
fProlloff = [];%1          % roll off frequency for filtering
sPAC.HighFreqs = nestedCenters;
nFa = length(nestedCenters);
nSources = size(Xinput,1);
isources = 1:nSources;
nTime = fix((nTS-fix(winLen*sRate))/fix(tStep*sRate))+1;
TimeOut = winLen/2 : tStep : winLen/2+(nTime-1)*tStep;        % Sec
nSur = Options.nSur;                % Number of surrogate itterations
PAC = zeros(nFa,nTime,nSources,nSur);               % PAC measure
nestingFreq = zeros(nFa,nTime,nSources,nSur);
DynamicPhase= zeros(nFa,nTime,nSources);               % PAC measure

% ===== MAIN LOOP ON FA ===== %
% Filtering in Fa band before cutting into smaller time windows
% (=> Higherfrequency resolution + faster process)

for ifreq=1:nFa
    % fA band
    bandNested = [nestedCenters(ifreq)-Fstep(ifreq),nestedCenters(ifreq)+Fstep(ifreq+1)];

    % Filtering in fA band
    Xnested = bst_bandpass_hfilter(Xinput, sRate,bandNested(1), bandNested(2), isMirror, isRelax, fArolloff);    % Filtering
    Xnested = Xnested(:,nMargin-nHilMar+1:end-nMargin+nHilMar);               % Removing part of the margin

    % Hilbert transform
    Z = hilbert(Xnested')';

    % Phase and envelope detection
    nestedEnv_total = abs(Z);             % Envelope of nested frequency rhythms
    nestedEnv_total = nestedEnv_total(:,nHilMar:end-nHilMar);              % Removing the margin

    % Loop on Time
    for iTime=1:nTime
        X = Xinput(:, (iTime-1)*fix(tStep*sRate)+[1:fix((2*margin+winLen)*sRate)]);
        nestedEnv = nestedEnv_total(:, (iTime-1)*fix(tStep*sRate)+[1:fix(winLen*sRate)]);

        % Time vector and number of samples
        nSample = size(nestedEnv,2);
        nFreq = 2^ceil(log2(nSample)+1);

        % Extraction of nesting frequency
        Ffft = abs(fft(nestedEnv-repmat(mean(nestedEnv,2),1,nSample),nFreq,2)).^2/nSample;
        freq = linspace(0,sRate,nFreq);
        x1 = X(:,nMargin+1:nMargin+fix(winLen*sRate));
        FfftSig = abs(fft(x1-repmat(mean(x1,2),1,nSample),nFreq,2)).^2/nSample;
        %%%

        % Finding the corresponding frequency component
        ind = bst_closest([minExtracFreq, fpBand(2)], freq);
        if freq(ind(1))<(minExtracFreq-diff(freq(1:2)))
            ind(1) = ind(1)+1;
        end
        if ind(2)>fpBand(2)
            ind(2) = ind(2)-1;
        end

                % Add previous and next point to the interval to give the algorithm
        % to find the local peaks even if they are in the first and last
        % point of interst in the spectrum
        if ind(1)>1
            ind(1) = ind(1)-1;
        end
        ind(2) = ind(2)+1;

        if freq(ind(1))<(minExtracFreq-diff(freq(1:2)))
            ind(1) = ind(1)+1;
        end

        indm = zeros(nSources,1);
        for iSource=1:nSources

            % Extracting the peak from envelope's PSD and then confirming
            % with a peak on the original signal
            [pks_env,locs_env] = findpeaks(Ffft(iSource,ind(1):ind(2)),'SORTSTR','descend');
            [pks_orig, locs_orig] = findpeaks(FfftSig(iSource,ind(1):ind(2)),'SORTSTR','descend');  % To check if a peak close to the coupled fp is available in the original signal
%             clear pks_env pks_orig

            % Ignore small peaks
            pks_orig = pks_orig/max(pks_orig);
            locs_orig = locs_orig(pks_orig>0.1);

            % Confirming the peak
            max_dist = max(1.5/winLen,1.5);     % maximum acceptable distance between peaks in evelope and the original signal's PSD
            count = 1;
            check_pks = 1;
            fp_loc = [];
            while check_pks && count<=length(locs_env)
                index = bst_closest(freq(locs_env(count)), freq(locs_orig));
                if abs(freq(locs_orig(index))-freq(locs_env(count)))<=max_dist
                    fp_loc = locs_env(count);
                    check_pks = 0;
                else
                    count = count+1;
                end
            end
            % If peak is not approved or no peak
            if isempty(fp_loc)
                fp_loc = 0;    % arbitrary value for fp  ==> will set the pac value to zero
                missedPcount = missedPcount +1;
            end

            indm(iSource) = fp_loc(1);
            clear pks_env locs_env
        end

        nestingFreq(ifreq,iTime,isources) = freq(ind(1)+indm-1);
        bandNesting = [max([squeeze(nestingFreq(ifreq,iTime,isources))-bandNestingLen/2,zeros(size(nestingFreq,3),1)],[],2),...
            squeeze(nestingFreq(ifreq,iTime,isources))+bandNestingLen/2];
        bandNesting(bandNesting<.15)=.15;

% Filtering in fP band
        if length(unique(bandNesting(:,1)))==1 && length(unique(bandNesting(:,2)))==1
            Xnesting = bst_bandpass_hfilter(X, sRate,bandNesting(1,1), bandNesting(1,2), isMirror, isRelax, fProlloff);    % Filtering
        else
            Xnesting = zeros(size(X));
            for i=1:length(isources)
                Xnesting(i,:) = bst_bandpass_hfilter(X(i,:), sRate, bandNesting(i,1), bandNesting(i,2),isMirror, isRelax, fProlloff);    % Filtering
            end
        end
        Xnesting = Xnesting(:,nMargin-nHilMar+1:fix((margin+winLen)*sRate)+nHilMar);              % Removing part of the margin
        % Hilbert transform
        Z = hilbert(Xnesting')';
        % Phase detection
        nestingPh = angle(Z-repmat(mean(Z,2),1,size(Z,2)));    % Phase of nesting frequency
        nestingPh = nestingPh(:,nHilMar:fix(winLen*sRate)+nHilMar-1);              % Removing the margin

        for ii=1:length(isources)
            iphase = find(diff(sign(nestingPh(ii,:) - nestingPh(ii,1)))==-2 | ...
                sign(nestingPh(ii,2:end)-nestingPh(ii,1))==0 | ...
                -(diff(sign(nestingPh(ii,:) - nestingPh(ii,1)))-1).*diff(nestingPh(ii,:)-nestingPh(ii,1)) >6 )-1;
            if isempty(iphase)
                iphase = length(nestingPh(ii,:));
            end
            amplitude = nestedEnv(ii,1:max(iphase));
            phase = nestingPh(ii,1:max(iphase));

            % Block resampling N=2
            numpoints=length(amplitude); %% number of sample points in raw signal

            if N==2
                minskip=numpoints/4; %% time lag must be at least this big
                maxskip=numpoints-numpoints/4; %% time lag must be smaller than this
                skip = ceil(minskip + (maxskip-minskip-1)*rand(2*nSur,1));
                skip(skip>maxskip)=[];
                skip(skip<minskip)=[];
                skip=skip(1:nSur,1);

                for iSur=1:nSur
                    surrogate_amplitude=[amplitude(skip(iSur):end) amplitude(1:skip(iSur)-1)];
                    PAC(ifreq,iTime,isources(ii),iSur) = sum(surrogate_amplitude.*exp(1i*phase),2)./max(iphase)./sqrt(mean(amplitude.^2,2));
                    DynamicPhase(ifreq,iTime,isources(ii),iSur) = angle(PAC(ifreq,iTime,isources(ii)));

                    if indm(ii)==0 % Fp not confirmed and arbitrary value for fp
                        PAC(ifreq,iTime,isources(ii),iSur) = 0;
                    end
                end

            elseif N>2 & N<numpoints/20
                % Block resampling N>2
%                 N = 10;
                numpoints=length(amplitude); %% number of sample points in raw signal
                block_len = fix(numpoints/N);
                remaining = numpoints - block_len*N;
                blocked_amp = reshape(amplitude(1:end-remaining),[block_len,N]);

                for iSur=1:nSur
                    order = randperm(N);
                    blocked_amp = blocked_amp(:,order);
                    surrogate_amplitude = [blocked_amp(:)', amplitude(end-remaining+1:end)];
                    PAC(ifreq,iTime,isources(ii),iSur) = sum(surrogate_amplitude.*exp(1i*phase),2)./max(iphase)./sqrt(mean(amplitude.^2,2));
                    DynamicPhase(ifreq,iTime,isources(ii),iSur) = angle(PAC(ifreq,iTime,isources(ii)));

                    if indm(ii)==0 % Fp not confirmed and arbitrary value for fp
                        PAC(ifreq,iTime,isources(ii),iSur) = 0;
                    end
                end

            else
                warning('Number of blocks should be at least 2, and less than 5% of data points')
            end
        end
    end
      disp(['iFreq: ', num2str(ifreq),' / ',num2str(nFa)]);
end

% ===== EXTRACTING THE PAC RELATED VALUES ===== %
[PACmax,maxInd] = max(abs(PAC),[],1);
% Fnested  = squeeze(nestedCenters(maxInd))';
% Sind     = repmat((1:nSources), nTime, 1);           % Source indices
% Tind     = repmat((1:nTime)', 1, nSources);              % Time indices
% linInd   = sub2ind(size(PAC),maxInd(:),Tind(:),Sind(:));
% Fnesting = reshape(nestingFreq(linInd),nTime,nSources)';
% phase    = reshape(angle(PAC(linInd)),nTime,nSources)'/pi*180;
% PACmax   = squeeze(PACmax)';

% ===== Interpolation in time domain for smoothing the results ==== %
if doInterpolation
    % Interpolation of PAC
    if nSources>1
        [X,Y,Z,W] = ndgrid(nestedCenters,TimeOut,[1:nSources], 1:nSur);
        ny = linspace(TimeOut(1), TimeOut(end), 2*length(TimeOut)-1);
        if logCenters
            nx = logspace(log10(nestedCenters(1)), log10(nestedCenters(end)), 2*nFa-1);
        else
            nx = linspace(nestedCenters(1), nestedCenters(end), 2*nFa-1);
        end
        [nX,nY,nZ,nW] = ndgrid(nx,ny,[1:nSources], 1:nSur);
        PAC = interpn(X,Y,Z,W,abs(PAC),nX,nY,nZ,nW,'linear',0);
        % make the notations similar
        tmp = nx; nx = ny;ny = tmp; clear tmp;

    else
        [X,Y,Z] = meshgrid(TimeOut,nestedCenters,[1:nSur]);
        nx = linspace(TimeOut(1), TimeOut(end), 2*nTime-1);
        if logCenters
            ny = logspace(log10(nestedCenters(1)), log10(nestedCenters(end)), 2*nFa-1);
        else
            ny = linspace(nestedCenters(1), nestedCenters(end), 2*nFa-1);
        end
        [nX,nY,nZ] = meshgrid(nx,ny,[1:nSur]);
        PAC = interp3(X,Y,Z,abs(squeeze(PAC)),nX,nY,nZ,'linear',0);
        PAC = permute(PAC,[1,2,4,3]);
    end
    TimeOut = nx;
    sPAC.HighFreqs = ny;
    clear nx nX nY nZ X Y Z

    % nestingFreq
    tmp = zeros(nFa*2-1, nTime, nSources, nSur);
    tmp(1:2:end,:,:,:) = nestingFreq;
    tmp(2:2:end,:,:,:) = nestingFreq(1:end-1,:,:,:);
    tmp2 = zeros(nFa*2-1, nTime*2-1, nSources, nSur);
    tmp2(:,1:2:end,:,:) = tmp;
    tmp2(:,2:2:end,:,:) = tmp(:,1:end-1,:,:);
    nestingFreq = tmp2;

    tmp = zeros(nFa*2-1, nTime, nSources, nSur);
    tmp(1:2:end,:,:,:) = DynamicPhase;
    tmp(2:2:end,:,:,:) = DynamicPhase(1:end-1,:,:,:);
    tmp2 = zeros(nFa*2-1, nTime*2-1, nSources, nSur);
    tmp2(:,1:2:end,:,:) = tmp;
    tmp2(:,2:2:end,:,:) = tmp(:,1:end-1,:,:);
    DynamicPhase = tmp2;
    clear tmp tmp2

end

Fnesting = [];
Fnested = [];
phase = [];

if missedPcount>0
disp(['Missed Peaks:',num2str(missedPcount),'/',num2str(nFa*nTime*nSources)])
end

% ===== OUTPUTS ===== %
if nTime >1
    sPAC.ValPAC = PACmax;
    sPAC.NestingFreq = Fnesting;
    sPAC.NestedFreq  = Fnested;
    sPAC.PhasePAC = phase;
    sPAC.TimeOut  = TimeOut;
    sPAC.DynamicPAC(:,:,1:nSources,:) = abs(PAC);
    sPAC.DynamicNesting(:,:,1:nSources,:)  = nestingFreq;
    sPAC.DynamicPhase(:,:,1:nSources,:)  = DynamicPhase;

        % == Generating two time points for Brainstorm structure ==
else
    sPAC.ValPAC = [PACmax(:), PACmax(:)];
    sPAC.NestingFreq = [Fnesting(:), Fnesting(:)];
    sPAC.NestedFreq  = [Fnested(:), Fnested(:)];
    sPAC.PhasePAC = [phase(:), phase(:)];
    sPAC.TimeOut  = [TimeOut, TimeOut+0.001];
    sPAC.DynamicPAC(:,1:2,1:nSources,:) = repmat(abs(PAC),[1,2,1]);
    sPAC.DynamicNesting(:,1:2,1:nSources,:)  = repmat(abs(nestingFreq),[1,2,1]);
    sPAC.DynamicPhase(:,1:2,1:nSources,:)  = repmat(abs(DynamicPhase),[1,2,1]);
end

end

author: brainstorm-tools


text:
function genPreLitTableLatex

p = path;
if exist('C:\Users\Luis\Dropbox\MATL', 'dir')
    baseFolder = 'C:\Users\Luis\Dropbox\MATL';
elseif exist('C:\Users\lmendo\Dropbox\MATL', 'dir')
    baseFolder = 'C:\Users\lmendo\Dropbox\MATL';
else
    error('Base folder not found')
end
addpath([baseFolder '\compiler'])

preLitTableFileName = [ baseFolder '\spec\preLitTable\preLitTable.tex' ];
preLitMatFileName = [ baseFolder '\compiler\preLit.mat' ];
preLitTxtFileName = [ baseFolder '\compiler\preLit.txt' ];

% Update preLit.mat, if needed; or load file
pLTxt = dir(preLitTxtFileName);
if ~isempty(pLTxt)
    pLTxt = datenum(pLTxt.date);
else
    pLTxt = -inf;
end
pLMat = dir(preLitMatFileName);
if ~isempty(pLMat)
    pLMat = datenum(pLMat.date);
else
    pLMat = -inf;
end
if pLMat < pLTxt % regeneration of `L.mat` required
    L = genPreLit(preLitTxtFileName, preLitMatFileName); % creates `preLit` struct array
    % definitions from funDefTxtFileName file, and saves it in funDefMatFileName file
else % load `preLit` variable from funDefMatFileName directly
    load(preLitMatFileName) % loads L variable
end

fn = fieldnames(L);
for n = 1:numel(fn)
    lines = {};
    %lines{end+1} = '\begin{tabular}{|l|l||l|l||l|l|}';
    %lines{end+1} = '\begin{tabular}{|p{.03\textwidth}|p{.3\textwidth}||p{.03\textwidth}|p{.3\textwidth}||p{.03\textwidth}|p{.3\textwidth}|}';
    %lines{end+1} = '\begin{longtable}{|l|l||l|l||l|l|}';
    lines{end+1} = '\begin{longtable}{|p{.03\textwidth}|p{.29\textwidth}||p{.03\textwidth}|p{.29\textwidth}||p{.03\textwidth}|p{.29\textwidth}|}';
    %lines{end+1} = '\esptab';
    lines{end+1} = '\hline';
    lines{end+1} = sprintf('\\multicolumn{6}{|c|}{Function \\matl{%s}} \\\\ \\hline\\hline', fn{n});
    %lines{end+1} = '\matl{99} & \matl{''CollapseDelimiters''} & \matl{99} & \matl{''CollapseDelimiters''} & \matl{99} & \matl{''CollapseDelimiters''} \kill';
    keys = L.(fn{n}).key;
    vals = L.(fn{n}).val;
    for m = 1:numel(keys)
        r = mod(m,3);
        if r
            lines{end+1} = sprintf('\\matl|%s| & \\matl|%s| &', num2str(keys(m)), vals{m}); % this assumes "|" is not be used in literals
        else
            lines{end+1} = sprintf('\\matl|%s| & \\matl|%s| \\\\ \\hline', num2str(keys(m)), vals{m});
        end
    end
    if r==1
        lines{end+1} = sprintf('& & & \\\\ \\hline');
    elseif r==2
        lines{end+1} = sprintf('& \\\\ \\hline');
    end
    lines{end+1} = '\end{longtable}';
    %lines{end+1} = '\ ';

    % Write file
    [p, f, e] = fileparts(preLitTableFileName);
    fname = [p filesep f fn{n} e];
    fid = fopen(fname, 'w');
    for k = 1:numel(lines)
        fwrite(fid, lines{k});
        fwrite(fid, [13 10]); % line break, Windows style
    end
    fclose(fid);

    % Convert "\matl" contents:
    convert_matl(fname, false);
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
%SF_DISC1 Linear discontinuous shape function (P1).
%
%   [ VBASE, NLDOF, XLDOF, SFUN ] = SF_DISC1( I_EVAL, N_SDIM, N_VERT, I_DOF, XI, AINVJAC, VBASE )
%   Evaluates discontinuous linear shape functions with with values
%   defined by a constant and derivatives in the cell interior.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       i_eval      scalar:  1             Evaluate function values
%                           >1             Evaluate values of derivatives
%       n_sdim      scalar: 1-3            Number of space dimensions
%       n_vert      scalar: 2-8            Number of vertices per cell
%       i_dof       scalar: 1-n_ldof       Local basis function to evaluate
%       xi          [n_sdim(+1)]           Local coordinates of evaluation point
%       aInvJac     [n,n_sdim(+1)*n_sdim]  Inverse of transformation Jacobian
%       vBase       [n]                    Preallocated output vector
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       vBase       [n]                    Evaluated function values
%       nLDof       [4]                    Number of local degrees of freedom on
%                                          vertices, edges, faces, and cell interiors
%       xLDof       [n_sdim,n_ldof]        Local coordinates of local dofs
%       sfun        string                 Function name of called shape function
%
%   See also SF_DISC0

% Copyright 2013-2022 Precise Simulation, Ltd.

% Check for evaluation on simplicies.

author: precise-simulation

repo: featool-multiphysics


text:
function [filteredPrice,filteredTime,actualTime] = realized_price_filter(price,time,timeType,samplingType,samplingInterval)
% Price filtering for computing realized variances and other
% high-frequency price based returns
%
% USAGE:
%   [FILTEREDPRICE,FILTEREDTIME,ACTUALTIME] =  realized_price_filter(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINTERVAL)
%
% INPUTS:
%   PRICE             - m by 1 column vector of prices
%   TIME              - m by 1 column vector of times corresponding to PRICE.
%                         Must be sorted and increasing.
%   TIMETYPE          - String describing the way times are measured
%                        'wall'    24-hour clock of the form HHMMSS, e.g. 101543 or 153217
%                        'seconds' Time measures in seconds past midnight.
%                          TIME must satisfy 0<=TIME<86400
%                        'unit'  Unit normalized date format, e.g. .1, .234, .9
%                          Unit normalized times are more general than the
%                          other types and can be applied to data from more
%                          than one calendar day
%   SAMPLINGTYPE      - String describing the type of sampling to use when
%                         filtering PRICE
%                         'CalendarTime' - Sample in calendar time using
%                           observations separated by SAMPLINGINTERVAL
%                           seconds. If TIMETYPE is 'unit',
%                           SAMPLINGINTERVAL must be between 0 and 1 and
%                           represents the fraction of the sample to skip
%                           when sampling.
%                         'CalendarUniform' - Sample in calendar time using
%                           SAMPLINGINTERVAL observations spread uniformly
%                           between TIME(1) and TIME(m)
%                         'BusinessTime' - Sample in business (tick) time
%                           using observation separated by SAMPLINGINTERVAL
%                           ticks
%                         'BusinessUniform' - Sample in business (tick)
%                           time using observations uniformly spaced in
%                           business time.
%                         'Fixed' - Sample at specific points in time. When
%                           using fixed, SAMPLINGINTERVAL must be a n by 1 vector
%                           of times with the same TIMETYPE as TIME (i.e.
%                           seconds if TIME is in seconds)
%   SAMPLINGINTERVAL   - Scalar integer or n by 1 vector whose meaning depends on the
%                          selected SAMPLINGTYPE
%
% OUTPUTS:
%   FILTEREDPRICE      - n by 1 vector of filtered prices.  n depends on
%                          SAMPLINGTYPE as well as the data
%   FILTEREDTIME       - n by 1 vector of sampling times corresponding to
%                          the chosen sampling scheme
%   ACTUALTIME         - n by 1 vector. ACTUALTIME(i) provides the actual
%                          observation time that was used when sampling at FILTEREDTIME(i)
%
%
% COMMENTS:
%   This is a helper function for REALIZED_KERNEL.  Last price
%   interpolation is always used.  If SAMPLINGTYPE is 'Fixed' then any
%   requested samples before the first observation will be back filled
%   with PRICE(1). An indicator for back filled values can be constructed
%   using ACTUALTIME>FILTEREDTIME
%
% EXAMPLES:
%   % 5-minute prices
%   FILTEREDPRICE =  realized_price_filter(PRICE,TIME,'wall','CalendarTime',300)
%
%   % 79 prices spread uniformly in calendar time
%   FILTEREDPRICE =  realized_price_filter(PRICE,TIME,'wall','CalendarUniform',79)
%
%   % 20-tick returns
%   FILTEREDPRICE = realized_price_filter(PRICE,TIME,'wall','CalendarUniform',79)
%
%   % 391 prices spread uniformly in business time
%   FILTEREDPRICE = realized_price_filter(PRICE,TIME,'wall','CalendarUniform',79)
%
%   % 5-minute prices coresponding to 9:30:00, 9:35:00, ...
%   sampleTimes =  seconds2wall(wall2seconds(93000):300:wall2seconds(160000))';
%   FILTEREDPRICE = realized_price_filter(PRICE,TIME,'wall','Fixed',sampleTimes)
%
%   % .05-calendar time sampling when using unit times
%   TIME = wall2unit(TIME,min(TIME),max(TIME));
%   FILTEREDPRICE = realized_price_filter(PRICE,TIME,'unit','CalendarTime',.05)
%
%
%  See also REALIZED_KERNEL, REALIZED_KERNEL_CORE, REALIZED_KERNEL_WEIGHTS
%  REALIZED_KERNEL_SELECT_LAG_LENGTH, REALIZED_VARIANCE

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 5/1/2008

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin~=5
    error('Five inputs required.')
end
if size(price,2)>size(price,1)
    price=price';
end
if size(price,2)>1
    error('PRICE must be a m by 1 vector.')
end
if size(time,2)>size(time,1)
    time=time';
end
if any(diff(time)<0)
    error('TIME must be sorted and increasing')
elseif any(diff(time)==0)
    warning('oxfordRealized:realizedPriceFilter','TIME contains multiple entries with the same value. This creates an ambiguity and FILTEREDPRICE will contain the last price if TIME does not only unique elements.')
end
if size(time,2)>1 || length(time)~=length(price)
    error('TIME must be a m by 1 vector.')
end
timeType=lower(timeType);
if ~ismember(timeType,{'wall','seconds','unit'})
    error('TIMETYPE must be one of ''wall'', ''seconds'' or ''unit''.');
end
% Inserted to protect against inputing integer times
time = double(time);

samplingType=lower(samplingType);
if ~ismember(samplingType,{'calendartime','calendaruniform','businesstime','businessuniform','fixed'})
    error('SAMPLINGTYPE must be one of ''CalendarTime'', ''CalendarUniform'', ''BusinessTime'', ''BusinessUniform'' or ''Fixed''.');
end

m=size(price,1);
t0Original=time(1);
tTOriginal=time(m);
if ismember(samplingType,{'calendartime','calendaruniform','businesstime','businessuniform'})
    % Must be a scalar integer, unless using unit times
    if (~isscalar(samplingInterval) || floor(samplingInterval)~=samplingInterval || samplingInterval<1) && ~strcmp(timeType,'unit')
        error('SAMPLINGINTERVAL must be a positive integer for the SAMPLINGTYPE selected.')
    end
else
    if size(samplingInterval,2)>size(samplingInterval,1)
        samplingInterval=samplingInterval';
    end
    if ~(any(samplingInterval>=t0Original) && any(samplingInterval<=tTOriginal))
        error('At least one sampling interval must be between min(TIME) and max(TIME) when using ''Fixed'' as SAMPLINGTYPE.')
    end
    if any(diff(samplingInterval)<=0)
        error('When using ''Fixed'' as SAMPLINGTYPE the vector of sampling times in SAMPLINGINTERVAL must be sorted and strictly increasing.')
    end
end

if strcmp(timeType,'unit') && strcmp(samplingType,'calendartime')
    % samplingInterval must be between 0 and 1
    if samplingInterval>1
        error('When TIMETYPE is ''unit'' and SAMPLINGTYPE is ''CalendarTime'', SAMPLINGINTERVAL must also be in ''unit'' terms, and so must be between 0 and 1')
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Get the size of the price
m=size(price,1);

% Convert times
if strcmp(samplingType,'fixed')
    time0 = min(samplingInterval);
    time1 = max(samplingInterval);
else
    time0 = min(time);
    time1 = max(time);
end
if strcmp(timeType,'wall')
    time = wall2unit(time,time0,time1);
elseif strcmp(timeType,'seconds')
    time = seconds2unit(time,time0,time1);
end
% If fixed convert these times as well
if strcmp(samplingType,'fixed')
    if strcmp(timeType,'wall')
        samplingInterval = wall2unit(samplingInterval,time0,time1);
    elseif strcmp(timeType,'seconds')
        samplingInterval = seconds2unit(samplingInterval,time0,time1);
    end
end

% Get the first and last time
t0=time(1);
tT=time(m);
% Switch based on the cases
switch samplingType
    case 'calendartime'
        % If calendar time
        % Convert sampling interval to be a fraction of a period
        % Unit strides
        if strcmp(timeType,'wall')
            timeBase = wall2seconds([time0 time1]);
            samplingInterval = samplingInterval/diff(timeBase);
        elseif strcmp(timeType,'seconds')
            samplingInterval = samplingInterval/(time1-time0);
        end
        filteredTime=(t0:samplingInterval:tT)';
        % Append the final observation
        if ~ismember(tT,filteredTime)
            filteredTime=[filteredTime; tT];
        end
        % Compute the filtered prices
        [filteredPrice,actualTime]=fasttimefilter(price,time,filteredTime);
    case 'calendaruniform'
        % Sampling interval contains the number of samples
        % Uniform spacing
        filteredTime=linspace(t0,tT,samplingInterval)';
        % Compute the filtered prices
        [filteredPrice,actualTime]=fasttimefilter(price,time,filteredTime);
    case 'businesstime'
        indices=(1:samplingInterval:m)';
        if ~ismember(m,indices)
            indices=[indices;m];
        end
        filteredPrice=price(indices);
        filteredTime=time(indices);
        actualTime = filteredTime;
    case 'businessuniform'
        % Sampling interval contains the number of samples
        indices=floor(linspace(1,m,samplingInterval));
        filteredPrice=price(indices);
        filteredTime=time(indices);
        actualTime = filteredTime;
    case 'fixed'
        filteredTime = samplingInterval;
        % Compute the filtered prices
        [filteredPrice,actualTime]=fasttimefilter(price,time,filteredTime);
    otherwise
        error('Unrecogniced SAMPLINGTYPE.')
end

% Finally convert filteredTime back to the original time format
if strcmp(timeType,'wall')
    filteredTime = unit2wall(filteredTime,time0,time1);
    actualTime = unit2wall(actualTime,time0,time1);
elseif strcmp(timeType,'seconds')
    filteredTime = unit2seconds(filteredTime,time0,time1);
    actualTime = unit2seconds(actualTime,time0,time1);
end

function [filteredPrice,actualTime]=fasttimefilter(price,time,filteredTime)

% Get the size of the inputs
m=size(price,1);
n=size(filteredTime,1);
% This uses a 2-index algorithm so that it is fast for any realistic size.
timeIndex=1;
filteredTimeIndex=1;
% Pl hold the index values.  Initializing to 1 makes the back fill easy
pl=ones(n,1);
while timeIndex<=m && filteredTimeIndex<=n
    if time(timeIndex)<=filteredTime(filteredTimeIndex)
        pl(filteredTimeIndex)=timeIndex;
        timeIndex=timeIndex+1;
    elseif filteredTimeIndex<n
        % Increment filteredTimeIndex until filteredTime>=time(timeIndex)
        while filteredTimeIndex<n && filteredTime(filteredTimeIndex)<time(timeIndex)
            filteredTimeIndex = filteredTimeIndex + 1;
            % Last price interploation for these since there is no new
            % information
            pl(filteredTimeIndex) = pl(filteredTimeIndex-1);
        end
        % Only assign if you aren't on the last one
        if filteredTimeIndex<n
            pl(filteredTimeIndex)=timeIndex;
        end
    elseif time(timeIndex)>filteredTime(n)
        % No more to assign, so break
        break
    end
end

% Clean up any trailing ones
starter = find(diff(pl)<0) + 1;
if ~isempty(starter)
    pl(starter:length(pl)) = pl(starter-1);
end
% Use pl the index the filteredPrice and filterTimeActual
filteredPrice = price(pl);
actualTime = time(pl);

author: bashtage

repo: mfe-toolbox


text:
function b = eq(f1,f2,varargin)
% ? sS1 == sS2
%
% Input
%  f1, f2 - @fibre
%
% Output
%  b - boolean
%
% Options
%  antipodal - include antipodal symmetry
%

if f1.CS ~= f2.CS || ...
    (length(f1)>1 && length(f2)>1 && any(size(f2)~=size(f1)))
  b = false;
  return
end

b = f1.h == f2.h & f1.r == f2.r;

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
classdef VideoMaker_Physical_Problem < VideoMaker

     properties

     end

      methods (Access = public)
        function Make_video_stress(obj,output_video_name,component)
            post = Postprocess_PhysicalProblem;
            field2print = post.stress_name;
            componentfield = [post.stress_component,component];
            obj.Make_video_standard_field(field2print,componentfield,output_video_name)
        end

        function Make_video_strain(obj,output_video_name,component)
            post = Postprocess_PhysicalProblem;
            field2print = post.strain_name;
            componentfield = [post.strain_component,component];
            obj.Make_video_standard_field(field2print,componentfield,output_video_name)
        end

        function Make_video_displacement(obj,output_video_name,component)
            post = Postprocess_PhysicalProblem;
            field2print = post.displ_name;
            componentfield = [post.displ_component,component];
            obj.Make_video_standard_field(field2print,componentfield,output_video_name)
         end

      end

      methods (Access = private)
        function Make_video_standard_field(obj,field2print,componentfield,output_video_name_in)
            file_tcl_name = 'tcl_gid.tcl';
            file_list = obj.create_file_list(obj.iterations_to_print,obj.file_name,obj.files_folder);
            [output_video_name] = obj.replace_special_character(output_video_name_in);
            [file_list] = obj.replace_special_character(file_list);

            file_tcl_name_with_path = fullfile(obj.files_folder,file_tcl_name);
            file_path_in = fullfile(pwd,'FEM','PostProcess','Make_Video_stress.tcl');
            filepath = obj.replace_special_character(file_path_in);

            fid = fopen(file_tcl_name_with_path,'w+');

            fprintf(fid,'GiD_Process PostProcess \n');
            fprintf(fid,'%s\n',['set arg1 "',file_list,'"']);
            fprintf(fid,'%s\n',['set arg2 "',output_video_name,'"']);
            fprintf(fid,'%s\n',['set arg3 "',field2print,'"']);
            fprintf(fid,'%s\n',['set arg4 "',componentfield,'"']);

            fprintf(fid,'%s\n',['source "',filepath,'"']);
            fprintf(fid,'%s\n',['Make_Video_stress $arg1 $arg2 $arg3 $arg4']);
            fprintf(fid,'%s\n',['GiD_Process Mescape Quit']);
            fclose(fid);
            obj.executeTclFiles(obj.gidPath,file_tcl_name_with_path)
        end

      end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [node,elem,indexMap] = rmisopoint(node,elem,gflag)
%% RMISOPOINT remove isolate points
%
% [node,elem] = rmisopoint(node,elem) rmove isolate points in the
% triangulation (node,elem). Isolate points are interiori vertices with
% valence 3 and boundary vertices with valence 2. Here the valence of a
% vertex is defined as the triangles containing that vertex.
%
% [node,elem,indexMap] = rmisopoint(node,elem) also return the index map
% between the input and output vertices. Since some vertices could be
% deleted in the procedure, the indexMap is important for the interpolation
% of functions defined on these vertices. See nodeinterpolate.
%
% Example
%   load airfoilperturb
%   [node,elem] = rmisopoint(node,elem);
%
% See also  bdsmoothing, edgeswap, nodeinterpolate
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

N = size(node,1); NT = size(elem,1);
%% Find boundary nodes
bdNode = findboundary(elem);
isBdNode = false(N,1);
isBdNode(bdNode) = true;
isIntNode = ~isBdNode;
intNode = find(~isBdNode);

%% Compute valence and find isolate points
valence = accumarray(elem(:),ones(3*NT,1),[N 1]);
isIsoNode = false(N,1);
isIsoNode(bdNode(valence(isBdNode)==2)) = true;
isIsoNode(intNode(valence(isIntNode)==3)) = true;
isoNode = find(isIsoNode);

%% Construct node star of isolated points
isIsoElem = isIsoNode(elem(:,1))|isIsoNode(elem(:,2))|isIsoNode(elem(:,3));
isoElem = find(isIsoElem);
tt = elem(isIsoElem,:);
q = simpqual(node,tt);
Ntt = size(tt,1);
tt2v = sparse([1:Ntt,1:Ntt,1:Ntt], tt(1:Ntt,:), 1, Ntt, N);

%% Remove isolate points
for i = 1:length(isoNode)
    pi = isoNode(i);
    ring = find(tt2v(:,pi));
    qi = min(q(ring));
    if qi>0.65  % don't remove this point
        isIsoNode(isoNode(i)) = false;
    else % remove this isolate point
        allpt = tt(ring,:);
        linkpt = setdiff(unique(allpt(:)),pi);
        if isIntNode(pi) % interior isolate points
            newt = fixorder(node,linkpt');
            tt(ring(1),:) = newt;
            tt(ring(2:end),:) = 0;
            % isoElem is used as an index mapping from tt to t
            elem(isoElem(ring(1)),:) = newt;
            elem(isoElem(ring(2:end)),1) = 0;
        elseif (isBdNode(pi) && isempty(find(linkpt == 0,1)))
        % Corner/feature points could be isolate nodes. Remove these
        % points will change the shape (and thus the area) of the mesh.
            newt = fixorder(node,linkpt');
            newarea = simplexvolume(node,newt);
            oldarea = sum(simplexvolume(node,tt(ring,:)));
            newq = simpqual(node,newt);
            if (abs(newarea-oldarea)/oldarea < 1e-3) && (newq>qi)
                tt(ring(1),:) = newt;
                q(ring(1)) = newq;
                tt(ring(2:end),:) = 0;
                elem(isoElem(ring(1)),:) = newt;
                elem(isoElem(ring(2:end)),1) = 0;
            else
                isIsoNode(isoNode(i)) = false;
            end
        else % linkpt contains 0, which means some triangles are deleted
            isIsoNode(isoNode(i)) = false;
        end
    end
end
%% Graph
if (nargin>2) && (gflag == 1)
    showmesh(node,elem);
    hold on; findnode(node,isIsoNode,'noindex','color','r')
end
%% Clean up
elem((elem(:,1) == 0),:) = [];
node(isIsoNode,:) = [];
indexMap = zeros(N,1);
indexMap(~isIsoNode)= 1:size(node,1);
elem = indexMap(elem);

repo: Aero-matlab


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% triple.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [xtrip,ftrip,g,G,x1,x2,nf] = triple(fcn,data,x,f,x1,x2,u,v,G)
% finding a local quadratic model by triple search
% Input:
% fcn = 'fun'	name of function fun(data,x), x an n-vector
% x(1:n)	start vector
% f		its function value
% x1(1:n), x2(1:n)  'neighbors' of x with x(i), x1(i), x2(i) pairwise
%		distinct for i = 1,...,n
% [u,v]		original box
% G		Hessian of the quadratic model around xtrip
% Output:
% xtrip(1:n)	best point found
% ftrip		its function value
% g, G		gradient and Hessian of the quadratic model around xtrip
% the quadratic model is given by
% q(x) = ftrip + g^T(x-xtrip) + 0.5(x-xtrip)G(x-xtrip)
% Uses the following m-files:
% hessian.m
% polint1.m

function [xtrip,ftrip,g,G,x1,x2,nf] = triple(fcn,data,x,f,x1,x2,u,v,hess,G)
nf = 0;
n = length(x);
g = zeros(n,1);
if nargin < 9
  hess = ones(n,n);
end
ind = find(u < x & x < v);
ind1 = find(x <= u | x >= v);
for j=1:length(ind1)
  g(ind1(j)) = 0;
  for k=1:n
    G(ind1(j),k) = 0;
    G(k,ind1(j)) = 0;
  end
end
if length(ind) <= 1
  xtrip = x;
  ftrip = f;
  if ~isempty(ind)
    g(ind) = 1;
    G(ind,ind) = 1;
  end
  return
end
if nargin < 10
  G = zeros(n,n);
end
xtrip = x;
ftrip = f;
xtripnew = x;
ftripnew = f;
for j=1:length(ind)
  i = ind(j);
  x = xtrip;
  f = ftrip;
  x(i) = x1(i);
  f1 = feval(fcn,data,x);
  x(i) = x2(i);
  f2 = feval(fcn,data,x);
  nf = nf + 2;
  [g(i),G(i,i)] = polint1([xtrip(i) x1(i) x2(i)],[f f1 f2]);
  if f1 <= f2
    if f1 < ftrip
      ftripnew = f1;
      xtripnew(i) = x1(i);
    end
  else
    if f2 < ftrip
      ftripnew = f2;
      xtripnew(i) = x2(i);
    end
  end
  if nargin < 10
    k1 = 0;
    if f1 <= f2
      x(i) = x1(i);
    else
      x(i) = x2(i);
    end
    for k=1:i-1
      if hess(i,k)
        if xtrip(k) > u(k) & xtrip(k) < v(k) & ~isempty(find(ind==k))
          q1 = ftrip + g(k)*(x1(k)-xtrip(k))+0.5*G(k,k)*(x1(k)-xtrip(k))^2;
          q2 = ftrip + g(k)*(x2(k)-xtrip(k))+0.5*G(k,k)*(x2(k)-xtrip(k))^2;
          if q1 <= q2
            x(k) = x1(k);
          else
            x(k) = x2(k);
          end
          f12 = feval(fcn,data,x);
          nf = nf + 1;
          G(i,k) = hessian(i,k,x,xtrip,f12,ftrip,g,G);
          G(k,i) = G(i,k);
          if f12 < ftripnew
            ftripnew = f12;
            xtripnew = x;
            k1 = k;
          end
          x(k) = xtrip(k);
        end
      else
        G(i,k) = 0;
        G(k,i) = 0;
      end
    end
  end
  if ftripnew < ftrip
    if x1(i) == xtripnew(i)
      x1(i) = xtrip(i);
    else
      x2(i) = xtrip(i);
    end
    if nargin < 10 & k1 > 0
      if xtripnew(k1) == x1(k1)
        x1(k1) = xtrip(k1);
      else
        x2(k1) = xtrip(k1);
      end
    end
    for k=1:i
      if ~isempty(find(ind==k))
        g(k) = g(k) + G(i,k)*(xtripnew(i) - xtrip(i));
        if nargin < 10 & k1 > 0
          g(k) = g(k) + G(k1,k)*(xtripnew(k1) - xtrip(k1));
        end
      end
    end
    xtrip = xtripnew;
    ftrip = ftripnew;
  end
end

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function rename_lowercase()
% renames all files in dir with lower case versions
% if all characters are capitalized.
%
% :Example:
% ::
%
%    d = dir('0*'); for i=1:length(d),cd(d(i).name),rename_lowercase, cd .., end
%
% ..
%    tor wager
% ..

d = dir; for i = 3:length(d),
    nm1 =(d(i).name);nm2=lower(nm1);
    str=['!mv ' nm1 ' ' nm2];

    if strcmp(nm1,'SPM.MAT'),
        str = ['!mv ' nm1 ' SPM.mat'];
        disp(str)
        eval(str)
    elseif strcmp(upper(nm1),nm1)
        disp(str)
        eval(str)
    end
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function saveTightFigure(h,outfilename)
% SAVETIGHTFIGURE(H,OUTFILENAME) Saves figure H in file OUTFILENAME without
%   the white space around it.
%
% by ``a grad student"
% http://tipstrickshowtos.blogspot.com/2010/08/how-to-get-rid-of-white-margin-in.html

% get the current axes
ax = get(h, 'CurrentAxes');

% make it tight
ti = get(ax,'TightInset');
set(ax,'Position',[ti(1) ti(2) 1-ti(3)-ti(1) 1-ti(4)-ti(2)]);

% adjust the papersize
set(ax,'units','centimeters');
pos = get(ax,'Position');
ti = get(ax,'TightInset');
set(h, 'PaperUnits','centimeters');
set(h, 'PaperSize', [pos(3)+ti(1)+ti(3) pos(4)+ti(2)+ti(4)]);
set(h, 'PaperPositionMode', 'manual');
set(h, 'PaperPosition',[0 0  pos(3)+ti(1)+ti(3) pos(4)+ti(2)+ti(4)]);

% save it
saveas(h,outfilename);

% Copyright (c) 2011, Emre Akbas
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * 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 OWNER 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.

author: alecjacobson

repo: gptoolbox


text:
function munu_l_i
%
%  Generates the data used in 'munu_l'. Data are stored in global
%  variables.
%
%  NOTE that 'munu_m_i' must be called before calling this routine.
%
%  Calling sequence:
%
%    munu_l_i
%
%  Remarks:
%
%    This routine has access to data, which is provided as global
%    data by the routine 'munu_m_i'.
%
%
%  LYAPACK 1.0 (Thilo Penzl, September 1999)

global LP_N LP_NL LP_NU

if ~length(LP_N)
  error('This routine needs global data which must be generated by calling ''munu_m_i'' first.');
end

[LP_NL,LP_NU] = lu(LP_N);

author: Sable

repo: mcbench-benchmarks


text:
function state = initstate(nvars,options,n,itr)
% Initialize swarm condition. Called by PSO.

% Initial particle velocities
state.Velocities = zeros(n,nvars) ;

% Initialize particle positions
state.Population = ...
    repmat(options.PopInitRange(1,:),n,1).*ones(n,nvars) + ...
    repmat((options.PopInitRange(2,:) - options.PopInitRange(1,:)),n,1) ...
    .*rand(n,nvars) ;

% Initialize the global and local fitness to the worst possible
state.fGlobalBest = ones(itr,1)*inf; % Global best fitness score
state.fLocalBests = ones(n,1)*inf ; % Individual best fitness score

% Initialize global and local best positions
state.xGlobalBest = ones(1,nvars)*inf ;
state.xLocalBests = ones(n,nvars)*inf ;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%%
clear;
close all;
%%
% load an image
x = uiuc_sample;
x = x(1:251, 1:256);

% compute filter bank
opt_filters = struct();
opt_filters.min_margin = [25, 800];
filters = morlet_filter_bank_2d(size(x), opt_filters);

%%
opt_wavelet = struct();
tic;
[x_phi, x_psi, meta_phi, meta_psi] = wavelet_2d(x, filters, opt_wavelet);
toc;
%%
% compute energy
energy_x = sum(x(:).^2);
energy_Wx = sum(x_phi(:).^2) + sum(cellfun(@(x)(sum(abs(x(:).^2))),x_psi));
assert(abs(energy_x-energy_Wx)/energy_x<1e-2);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function pt = bb_points(bb,numM,numN,margin)
% Generates numM x numN points on BBox.

bb(1:2) = bb(1:2)+margin;
bb(3:4) = bb(3:4)-margin;

if (numM == 1 && numN ==1)
    pt = bb_center(bb);
    return;
end

if (numM == 1 && numN > 1)
    c = bb_center(bb);
    stepW = (bb(3)-bb(1)) / (numN - 1);
    pt = ntuples(bb(1):stepW:bb(3),c(2));
    return;
end

if (numM > 1 && numN == 1)
    c = bb_center(bb);
    stepH = (bb(4)-bb(2)) / (numM - 1);
    pt = ntuples(c(1),(bb(2):stepH:bb(4)));
    return;
end

stepW = (bb(3)-bb(1)) / (numN - 1);
stepH = (bb(4)-bb(2)) / (numM - 1);

pt = ntuples(bb(1):stepW:bb(3),(bb(2):stepH:bb(4)));

if size(pt,2) < numM*numN
    count = numM * numN - size(pt, 2);
    app = repmat(pt(:,end), 1, count);
    pt = [pt, app];
end

author: yuxng


text:
function n = p37_n ( )

%*****************************************************************************80
%
%% P37_N returns the number of variables for problem 37.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = 2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [z varargout] = intersectm(varargin)

% INTERSECTM Find set intersection of several vectors
%
%   Z = INTERSECTM(A, B, C...) returns the values common to A, B, C... which can be numeric or character
%                              vectors or cell arrays of strings. Z is a vector sorted in ascending order.
%
%   Z = INTERSECTM(A, B, C...,'rows') when A, B, C... are matrices with the same number of columns returns
%                                     the rows common to the sets supplied.
%
%   [Z, IA, IB, IC...] = INTERSECTM(...) also returns column index vectors IA, IB, IC... such that Z = A(IA),
%                                        Z = B(IB), Z = C(IC)... (or Z = A(IA,:), Z = B(IB,:), Z = C(IC,:)...).
%
%
% EXAMPLES:
%   - Retrieve Z the set intersection of the elements in "a", "b" and "c" and the indexes IA, IB and IC
%
%       a =  1:5   ;
%       b = (3:9).';
%       c = -1:4   ;
%
%       [z, ia, ib, ic] = intersectm(a, b, c)
%       z =
%           3     4
%       ia =
%           3
%           4
%       ib =
%           1
%           2
%       ic =
%           5
%           6
%
% NOTE: IA, IB and IC have the same size as their corresponding inputs A, B and C while Z
%       has the size of A.
%
%   - Retrieve Z the set intersection of the rows "a", "b" and "c" and the indexes IA, IB and IC
%
%       a = [ 1:5 ;11:15];
%       b =  11:15;
%       c = [16:20;11:15];
%
%       [z, ia, ib, ic] = intersectm(a, b, c, 'rows')
%       z =
%           11    12    13    14    15
%       ia =
%           2
%       ib =
%           1
%       ic =
%           2
%
% Additional features:
% - FEX set functions page
%
% See also: INTERSECT, ISMEMBERM, SETDIFFM, SETXORM, UNIONM

% Author: Oleg Komarov (oleg.komarov@hotmail.it)
% Tested on R14SP3 (7.1) and on R2009b
% 23 jul 2010 - Created
% 30 jul 2010 - Added link to fex

% # inputs
error(nargchk(2,nargin,nargin))

% Check for 'rows' flag
IDXflag = strcmp('rows',varargin);
if any(IDXflag)
    % Extract flag
    flag     = varargin{ IDXflag};
    varargin = varargin(~IDXflag);
    % Check that all inputs have the same number of columns
    if logical(diff(cellfun('size',varargin,2)))
        error('intersectm:sameNumCols','A, B, C... must have the same number of columns.')
    end
else flag = [];
    % Check that all inputs are vectors
    if any(cellfun('prodofsize',varargin) ~= cellfun('length',varargin))
        error('intersectm:allVectors','A, B, C... must be all vectors or ''rows'' must be specified.')
    end
end

% Check same class
if ~cellfun('isclass',varargin(2:end), class(varargin{1}))
    error('intersectm:sameClass','A, B, C... should belong to the same class.')
end

% # of sets to intersect
numSets = numel(varargin);
if numSets == 1
    error('intersectm:noSet','Insufficient inputs: must supply at least two sets.')
end

% # outputs
error(nargoutchk(0,numSets+1,nargout))

% Find intersection set with ismemberm
if isempty(flag)
    z = varargin{1}(all(ismemberm(varargin{:}     ),2)  );
else
    z = varargin{1}(all(ismemberm(varargin{:},flag),2),:);
end

% Retrieve column index vectors
if nargout > 1
    varargout = cellfun(@(x) reshape(find(ismember(x,z,flag)),[],1),varargin(1:nargout-1),'un',0);
end
end


text:
function GSS = getWoLFScores(inputFile, kingdom)
% getWoLFScores
%   Call WoLF PSort to predict the sub-cellular localization of proteins.
%   The output can be used as input to predictLocalization. This function
%   is currently only available for Linux and requires Perl to be
%   installed. If one wants to use another predictor, see parseScores. The
%   function normalizes the scores so that the best score for each gene is
%   1.0.
%
%   Input:
%   inputFile	a FASTA file with protein sequences
%   kingdom     the kingdom of the organism, 'animal', 'fungi' or 'plant'
%
%   Output:
%   GSS         a gene scoring structure to be used in predictLocalization
%
%   Usage: GSS = getWoLFScores(inputFile, kingdom)

if ~isfile(inputFile)
    error('FASTA file %s cannot be found',string(inputFile));
end

kingdom=char(kingdom);
if ~any(strcmp(kingdom,{'animal','fungi','plant'}))
    EM='Allowed kingdoms are "animal", "fungi", and "plant"';
    dispEM(EM);
end

if ispc || ismac
    EM='This function currently runs only on Linux. Use parseScores if you want to use another predictor';
    dispEM(EM);
end

%Get the directory for RAVEN Toolbox
ravenPath=findRAVENroot();

%Temporary output name
outFile=tempname;
fid=fopen(outFile,'w');

%Do the prediction
[~, output]=unix(['perl "' ravenPath '/software/WoLFPSORT/bin/runWolfPsortSummary" ' kingdom ' < ' inputFile]);

%Save output and call the general parser
fprintf(fid,output);
fclose(fid);
GSS=parseScores(outFile,'wolf');

%Clean up
delete(outFile);
end

author: SysBioChalmers

repo: RAVEN


text:
function ef = cpf_vlim_event(cb_data, cx)
%CPF_VLIM_EVENT  Event function to detect bus voltage limit violations
%   EF = CPF_VLIM_EVENT(CB_DATA, CX)
%
%   CPF event function to detect bus voltage limits violations,
%   i.e. Vm <= Vmin or Vm >= Vmax.
%
%   Inputs:
%       CB_DATA : struct of data for callback functions
%       CX : struct containing info about current point (continuation soln)
%
%   Outputs:
%       EF : event function value

%   MATPOWER
%   Copyright (c) 2016-2017, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   and Ahmad Abubakar Sadiq, Federal University of Technology Minna, Nigeria
%   and Shrirang Abhyankar, Argonne National Laboratory
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% event function value is 2 nb x 1 vector equal to:
%%      [ Vmin - Vm ]
%%      [ Vm - Vmax ]

%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;

%% get updated MPC
d = cb_data;
mpc = cpf_current_mpc(d.mpc_base, d.mpc_target, ...
    d.Ybus, d.Yf, d.Yt, d.ref, d.pv, d.pq, cx.V, cx.lam, d.mpopt);

%% voltage magnitude violations
v_Vmin = mpc.bus(:, VMIN) - mpc.bus(:, VM);
v_Vmax = mpc.bus(:, VM) - mpc.bus(:, VMAX);

%% assemble event function value
ef = [v_Vmin;v_Vmax];

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [minFlux,maxFlux,Vmin,Vmax] = fluxVariability(model,optPercentage,osenseStr,rxnNameList,verbFlag, allowLoops)
%fluxVariability Performs flux variablity analysis
%
% [minFlux,maxFlux] = fluxVariability(model,optPercentage,osenseStr,rxnNameList,verbFlag, allowLoops)
%
%INPUT
% model             COBRA model structure
%
%OPTIONAL INPUTS
% optPercentage     Only consider solutions that give you at least a certain
%                   percentage of the optimal solution (Default = 100
%                   or optimal solutions only)
% osenseStr         Objective sense ('min' or 'max') (Default = 'max')
% rxnNameList       List of reactions for which FVA is performed
%                   (Default = all reactions in the model)
% verbFlag          Verbose output (opt, default false)
% allowLoops        Whether loops are allowed in solution. (Default = true)
%                   See optimizeCbModel for description
%
%OUTPUT
% minFlux           Minimum flux for each reaction
% maxFlux           Maximum flux for each reaction
%
%OPTIONAL OUTPUT
% Vmin          Matrix of column flux vectors, where each column is a
%               separate minimization.
% Vmax          Matrix of column flux vectors, where each column is a
%               separate maximization.
%

% Markus Herrgard  8/21/06 Original code.
% Ronan Fleming   01/20/10 Take the extremal flux from the flux vector,
%                          not from the objective since this is invariant
%                          to the value and sign of the coefficient
% Ronan Fleming   27/09/10 Vmin,Vmax

if (nargin < 2)
    optPercentage = 100;
end
if (nargin < 3)
    osenseStr = 'max';
end
if (nargin < 4)
    rxnNameList = model.rxns;
end
if (nargin < 5)
    verbFlag = false;
end
if (nargin < 6)
    allowLoops = true;
end
if (isempty(optPercentage))
    optPercentage = 100;
end
if (isempty(osenseStr))
    osenseStr = 'max';
end
if (isempty(rxnNameList))
    rxnNameList = model.rxns;
end

% LP solution tolerance
global CBT_LP_PARAMS
if (exist('CBT_LP_PARAMS', 'var'))
    if isfield(CBT_LP_PARAMS, 'objTol')
        tol = CBT_LP_PARAMS.objTol;
    else
        tol = 1e-6;
    end
    if isfield(CBT_LP_PARAMS, 'minNorm')
        minNorm = CBT_LP_PARAMS.minNorm;
    else
        minNorm = 0;
    end
else
    tol = 1e-6;
    minNorm = 0;
end

% Determine constraints for the correct space (0-100% of the full space)
if (sum(model.c ~= 0) > 0)
    hasObjective = true;
    optSol = optimizeCbModel(model,osenseStr, 0, allowLoops);
    if (optSol.stat > 0)
        objRxn = model.rxns(model.c~=0);
        if (strcmp(osenseStr,'max'))
            objValue = floor(optSol.f/tol)*tol*optPercentage/100;
        else
            objValue = ceil(optSol.f/tol)*tol*optPercentage/100;
        end
    else
        error('Infeasible problem - no optimal solution!');
    end
else
    hasObjective = false;
end

if (verbFlag == 1)
    h = waitbar(0,'Flux variability analysis in progress ...');
end
if (verbFlag > 1)
    fprintf('%4s\t%4s\t%10s\t%9s\t%9s\n','No','Perc','Name','Min','Max');
end

if (~isfield(model,'b'))
    model.b = zeros(size(model.S,1),1);
end
% Set up the general problem
[nMets,nRxns] = size(model.S);
rxnListFull = model.rxns;
LPproblem.c = model.c;
LPproblem.lb = model.lb;
LPproblem.ub = model.ub;
LPproblem.csense(1:nMets) = 'E';
LPproblem.csense = LPproblem.csense';
if hasObjective
    LPproblem.A = [model.S;columnVector(model.c)'];
    LPproblem.b = [model.b;objValue];
    if (strcmp(osenseStr,'max'))
        LPproblem.csense(end+1) = 'G';
    else
        LPproblem.csense(end+1) = 'L';
    end
else
    LPproblem.A = model.S;
    LPproblem.b = model.b;
end

% %solve to generate initial basis
LPproblem.osense = -1;
tempSolution = solveCobraLP(LPproblem);
LPproblem.basis = tempSolution.basis;

% Loop through reactions
maxFlux = zeros(length(rxnNameList), 1);
minFlux = zeros(length(rxnNameList), 1);

if length(minNorm)> 1 || minNorm > 0
    Vmin=zeros(nRxns,nRxns);
    Vmax=zeros(nRxns,nRxns);
    %minimizing the Euclidean norm gets rid of the loops, so there
    %is no need for a second slower MILP approach
    allowLoops=1;
else
    Vmin=[];
    Vmax=[];
end

solutionPool = zeros(length(model.lb), 0);
if ~exist('matlabpool') || (matlabpool('size') == 0) %aka nothing is active
    m = 0;
    for i = 1:length(rxnNameList)
        if mod(i,10) == 0, clear mex, end
        if (verbFlag == 1),fprintf('iteration %d.  skipped %d\n', i, round(m));end
        LPproblem.c = zeros(nRxns,1);
        rxnBool=strcmp(rxnListFull,rxnNameList{i});
        LPproblem.c(rxnBool) = 1; %no need to set this more than 1
        % do LP always
        LPproblem.osense = -1;
        LPsolution = solveCobraLP(LPproblem);
        %take the maximum flux from the flux vector, not from the obj -Ronan
        maxFlux(i) = LPsolution.full(LPproblem.c~=0);

        %minimise the Euclidean norm of the optimal flux vector to remove
        %loops -Ronan
        if length(minNorm)> 1 || minNorm > 0
            QPproblem=LPproblem;
            QPproblem.lb(LPproblem.c~=0)=maxFlux(i)-1e-12;
            QPproblem.ub(LPproblem.c~=0)=maxFlux(i)+1e12;
            QPproblem.c(:)=0;
            %Minimise Euclidean norm using quadratic programming
            if length(minNorm)==1
                minNorm=ones(nRxns,1)*minNorm;
            end
            QPproblem.F = spdiags(minNorm,0,nRxns,nRxns);
            %quadratic optimization
            solution = solveCobraQP(QPproblem);
            if isempty(solution.full)
                pause(eps)
            end
            Vmax(:,rxnBool)=solution.full(1:nRxns,1);
        end

        LPproblem.osense = 1;
        LPsolution = solveCobraLP(LPproblem);
        %take the maximum flux from the flux vector, not from the obj -Ronan
        minFlux(i) = LPsolution.full(LPproblem.c~=0);

        %minimise the Euclidean norm of the optimal flux vector to remove
        %loops
        %minimise the Euclidean norm of the optimal flux vector to remove
        %loops
        if length(minNorm)> 1 || minNorm > 0
            QPproblem=LPproblem;
            QPproblem.lb(LPproblem.c~=0)=maxFlux(i)-1e-12;
            QPproblem.ub(LPproblem.c~=0)=maxFlux(i)+1e12;
            QPproblem.c(:)=0;
            QPproblem.F = spdiags(minNorm,0,nRxns,nRxns);
            %Minimise Euclidean norm using quadratic programming
            if length(minNorm)==1
                minNorm=ones(nRxns,1)*minNorm;
            end
            QPproblem.F = spdiags(minNorm,0,nRxns,nRxns);
            %quadratic optimization
            solution = solveCobraQP(QPproblem);
            Vmin(:,rxnBool)=solution.full(1:nRxns,1);
        end

        if ~allowLoops
            if any( abs(LPproblem.c'*solutionPool - maxFlux(i)) < tol) % if any previous solutions are good enough.
                % no need to do anything.
                m = m+.5;
            else
                LPproblem.osense = -1;
                LPsolution = solveCobraMILP(addLoopLawConstraints(LPproblem, model));
                maxFlux(i) = LPsolution.obj/1000;
              end
            if any( abs(LPproblem.c'*solutionPool - minFlux(i)) < tol)
                m = m+.5;
                % no need to do anything.
            else
                LPproblem.osense = 1;
                LPsolution = solveCobraMILP(addLoopLawConstraints(LPproblem, model));
                minFlux(i) = LPsolution.obj/1000;
            end
        end
        if (verbFlag == 1)
            waitbar(i/length(rxnNameList),h);
        end
        if (verbFlag > 1)
            fprintf('%4d\t%4.0f\t%10s\t%9.3f\t%9.3f\n',i,100*i/length(rxnNameList),rxnNameList{i},minFlux(i),maxFlux(i));
        end
    end
else % parallel job.  pretty much does the same thing.
    parfor i = 1:length(rxnNameList)
        %if mod(i,10) == 0, clear mex, end
        %if (verbFlag == 1),fprintf('iteration %d.  skipped %d\n', i, round(m));end
        c = zeros(nRxns,1);
        c(strcmp(rxnListFull,rxnNameList{i})) = 1000;
        if allowLoops % do LP
            LPsolution = solveCobraLP(struct(...
                'A', LPproblem.A,...
                'b', LPproblem.b,...
                'lb', LPproblem.lb,...
                'ub', LPproblem.ub,...
                'csense', LPproblem.csense,...
                'c',c,...
                'osense',-1, ...
                'basis', LPproblem.basis ...
            ));
            %take the maximum flux from the flux vector, not from the obj -Ronan
            maxFlux(i) = LPsolution.full(c~=0);
            %LPproblemb.osense = 1;
            LPsolution = solveCobraLP(struct(...
                'A', LPproblem.A,...
                'b', LPproblem.b,...
                'lb', LPproblem.lb,...
                'ub', LPproblem.ub,...
                'csense', LPproblem.csense,...
                'c',c,...
                'osense',1, ... %only part that's different.
                'basis', LPproblem.basis ...
            ));
            minFlux(i) = LPsolution.full(c~=0);
        else
            LPsolution = solveCobraMILP(addLoopLawConstraints(struct(...
                'A', LPproblem.A,...
                'b', LPproblem.b,...
                'lb', LPproblem.lb,...
                'ub', LPproblem.ub,...
                'csense', LPproblem.csense,...
                'c',c,...
                'osense',-1 ...
            ), model));
            maxFlux(i) = LPsolution.obj/1000;

            LPsolution = solveCobraMILP(addLoopLawConstraints(struct(...
                'A', LPproblem.A,...
                'b', LPproblem.b,...
                'lb', LPproblem.lb,...
                'ub', LPproblem.ub,...
                'csense', LPproblem.csense,...
                'c',c,...
                'osense',1 ...
            ), model));%
            minFlux(i) = LPsolution.obj/1000;
        end
    end
end

if (verbFlag == 1)
	if ( regexp( version, 'R20') )
        	close(h);
	end
end

maxFlux = columnVector(maxFlux);
minFlux = columnVector(minFlux);

author: opencobra

repo: cobratoolbox