text:
function varargout = process_cohere1n_2021( varargin )
% PROCESS_COHERE1N_2021: Compute the coherence between all the pairs of signals, in one file.
%
% USAGE:   OutputFiles = process_cohere1n_2021('Run', sProcess, sInputA)
%                        process_cohere1n_2021('Test')

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Coherence NxN [2021]';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Connectivity';
    sProcess.Index       = 656;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;

    % === CONNECT INPUT
    sProcess = process_corr1n('DefineConnectOptions', sProcess, 1);
    % === REMOVE EVOKED REPONSE
    sProcess.options.removeevoked.Comment = 'Remove evoked response from each trial';
    sProcess.options.removeevoked.Type    = 'checkbox';
    sProcess.options.removeevoked.Value   = 0;
    sProcess.options.removeevoked.Group   = 'input';
    % === COHERENCE METHOD
    sProcess.options.cohmeasure.Comment = {...
        ['Magnitude-squared Coherence
' ...
        '|C|^2 = |Cxy|^2/(Cxx*Cyy)'], ...
        ['Imaginary Coherence (2019)
' ...
        'IC    = |imag(C)|'], ...
        ['Lagged Coherence (2019)
' ...
        'LC    = |imag(C)|/sqrt(1-real(C)^2)'], ...
        [' Imaginary Coherence (before 2019)
' ...
        ' IC    = imag(C)^2 / (1-real(C)^2) ']; ...
        'mscohere', 'icohere2019', 'lcohere2019', 'icohere'};
    sProcess.options.cohmeasure.Type    = 'radio_label';
    sProcess.options.cohmeasure.Value   = 'mscohere';
    % === WINDOW LENGTH
    sProcess.options.win_length.Comment = 'Window length for PSD estimation:';
    sProcess.options.win_length.Type    = 'value';
    sProcess.options.win_length.Value   = {1, 's', []};
    % === OVERLAP
    sProcess.options.overlap.Comment = 'Overlap for PSD estimation:' ;
    sProcess.options.overlap.Type    = 'value';
    sProcess.options.overlap.Value   = {50, '%', []};
    % === HIGHEST FREQUENCY OF INTEREST
    sProcess.options.maxfreq.Comment = 'Highest frequency of interest:';
    sProcess.options.maxfreq.Type    = 'value';
    sProcess.options.maxfreq.Value   = {60,'Hz',2};
    % === OUTPUT MODE 2021
    sProcess.options.outputmode.Comment = {'Save individual results (one output file per input file)', 'Average cross-spectra of input files (one output file)'; ...
                                           'input', 'avgcoh'};
    sProcess.options.outputmode.Type    = 'radio_label';
    sProcess.options.outputmode.Value   = 'input';
    sProcess.options.outputmode.Group   = 'output';
end

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

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputA) %#ok<DEFNU>
    % Input options
    OPTIONS = process_corr1n('GetConnectOptions', sProcess, sInputA);
    if isempty(OPTIONS)
        OutputFiles = {};
        return
    end

    % Metric options
    OPTIONS.Method = 'cohere';
    OPTIONS.RemoveEvoked  = sProcess.options.removeevoked.Value;
    OPTIONS.WinLen        = sProcess.options.win_length.Value{1};
    OPTIONS.MaxFreq       = sProcess.options.maxfreq.Value{1};
    OPTIONS.CohOverlap    = 0.50;  % First pre-define the overlap
    OPTIONS.CohMeasure    = sProcess.options.cohmeasure.Value;

    % Change the overlap if it is specified
    if isfield(sProcess.options, 'overlap') && isfield(sProcess.options.overlap, 'Value') && ...
       iscell(sProcess.options.overlap.Value) && ~isempty(sProcess.options.overlap.Value) && ~isempty(sProcess.options.overlap.Value{1})
        OPTIONS.CohOverlap = sProcess.options.overlap.Value{1}/100 ;
    end

    % Compute metric
    OutputFiles = bst_connectivity({sInputA.FileName}, [], OPTIONS);
end

%% ===== TEST FUNCTION =====
function Test() %#ok<DEFNU>
    % Start a new report
    bst_report('Start');
    % Get test datasets
    sFile = process_simulate_ar('Test'); % Fs = 1200 Hz
    % NOTES:
    % bst_cohn.m (2019) uses 2^nextpow2(round(Fs / MaxFreqRes)) samples
    % of DATA for the FFT, there is no zero padding
    %
    % bst_cohn_2021.m uses nWinLen = round(WinLen * Fs) samples of DATA,
    % that are zero padded to 2^nextpow2(nWinLen * 2) for the FFT
    %
    % To get the similar results with bst_cohn and bst_cohn_2021:
    % 1. Select a MaxFreqRes that leads to a power-of-2 number of samples,
    % thus, we can be sure that no extra data is used in the FFT.
    % 2. Use the duration associated to that MaxFreqRes as WinLen parameter
    % for bst_cohn_2021, with this we will compute coherence on the same data

    Fs = 1200;      % Default Fs for process_simulate_ar('Test')
    nSamples = 512; % Desired number of samples
    MaxFreqRes = Fs/nSamples; % Hz ==> round(Fs / MaxFreqRes) = 512 samples
    WinLen = 1 / MaxFreqRes;  % s  ==> 512 samples zero-padded to 1024

    % Coherence process with bst_cohn.m (2019)
    tic;
    sTmp = bst_process('CallProcess', 'process_cohere1n', sFile, [], ...
        'timewindow',   [], ...          % All the time in input
        'cohmeasure',   'mscohere', ...  % 1=Magnitude-squared, 2=Imaginary
        'overlap',      50, ...          % 50%
        'maxfreqres',   MaxFreqRes, ...  % VARIES
        'maxfreq',      [], ...          % No maximum frequency
        'outputmode',   1);              % Save individual results (one file per input file)
    t = toc;
    % Execution time
    bst_report('Info', 'process_cohere1n', sFile, sprintf('Execution time: %1.6f seconds', t));
    % Add tag
    bst_process('CallProcess', 'process_add_tag', sTmp.FileName, [], 'tag', '(2019)' );
    % Snapshot: spectrum
    bst_process('CallProcess', 'process_snapshot', sTmp, [], ...
        'target',       11, ...  % Connectivity matrix (image)
        'modality',     1, 'orient', 1, 'time', 0, 'contact_time', [-40, 110], 'contact_nimage', 16, ...
        'Comment',      [sTmp.Comment, ': (2019)']);

    % Coherence process with bst_cohn_2021.m
    tic;
    sTmp = bst_process('CallProcess', 'process_cohere1n_2021', sFile, [], ...
        'timewindow',   [], ...          % All the time in input
        'includebad',   1, ...
        'removeevoked', 0, ...
        'cohmeasure',   'mscohere', ...  % Magnitude-squared Coherence|C|^2 = |Cxy|^2/(Cxx*Cyy)
        'win_length',   WinLen, ...
        'overlap',      50, ...          % 50%
        'maxfreq',      [], ...
        'outputmode',   'input');        % Save individual results (one output file per input file)
    t = toc;
    % Execution time
    bst_report('Info', 'process_cohere1n_2021', sFile, sprintf('Execution time: %1.6f seconds', t));
    % Snapshot: spectrum
    bst_process('CallProcess', 'process_snapshot', sTmp, [], ...
        'target',         11, ...  % Connectivity matrix (image)
        'modality',       1, ...
        'orient',         1, ...
        'contact_nimage', 16, ...
        'Comment',      [sTmp.Comment, ': (2021)']);

    % Save and display report
    ReportFile = bst_report('Save', sTmp);
    bst_report('Open', ReportFile);
end

author: brainstorm-tools


text:
classdef HomogVoigtPlaneStressHomogenizer < SequentialLaminateHomogenizer

    properties (Access = private)
    end

    methods (Access = public)

        function obj = HomogVoigtPlaneStressHomogenizer(c0,c1,dir,mi,theta)
            obj.init(c0,c1,dir,mi,theta);
            obj.computeAnisotropicTensors();
            obj.homogenize()
            obj.transformToVoigt()
            obj.makeHomogenizedTensorPlaneStress();
        end
    end

    methods (Access = protected)

        function transformToVoigt(obj)
            t = obj.homogTensor;
            t = Tensor2VoigtConverter.convert(t);
            obj.homogTensor = t;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ lines ] = findlines( imge,map,theta,rho,minLen,maxLineNum,maxPeakNum,resol,maxGap )
%FINDLINE Find line features in the image using Hough transform
%   LINES = FINDLINES(IMGE,MAP,THETA,RHO,MINLEN,MAXLINENUM,MAXPEAKNUM,RESOL,MAXGAP)
%	IMGE is the binary edge image;
%	MAP,THETA,RHO is the output of HOUGH;
%	MINLEN is the minimum length of lines that will be found;
%	No more MAXLINENUM lines will be returned;
%	MAXPEAKNUM is the maximum peaks	that will be detected on the map;
%	MAXGAP is the maximum gap length that will be merged;
%	RESOL is the resolution of the detected lines, the smaller, the lines
%	must be further away from each other. A suggest value is 100.
%	LINES is a struct array with members: theta,rho,point1([row,col]),point2.
%
%	Yan Ke @ THUEE, 20110123, xjed09@gmail.com

% imge = edge(img);
% [map theta rho] = hough1(imge,[]);
[pri pti] = findpeaks(map,minLen,maxPeakNum,resol);
peakNum = length(pri);
prho = rho(pri); % peak points
pthe = theta(pti)*pi/180;
rhoThres = abs(rho(2)-rho(1))/2;

fitPts = cell(1,peakNum);
fitPtsNum = zeros(1,peakNum);
for p = 1:peakNum
	 % store points which are on the lines
	fitPts{p} = zeros(2,ceil(map(pri(p),pti(p))*1.2));
end
[Y X] = find(imge);
X = X-1; % The origin is the bottom-left corner.
M = size(imge,1);
Y = M-Y;

% find the points that's on the selected lines
for p = 1:length(X)
	rho1 = X(p)*sin(pthe)+Y(p)*cos(pthe);
	isOnPeak = find(abs(rho1-prho)<=rhoThres);
	fitPtsNum(isOnPeak) = fitPtsNum(isOnPeak)+1;
	for q = isOnPeak
		fitPts{q}(:,fitPtsNum(q)) = [X(p);Y(p)];
	end
end

% track the points on each line, find the start points and the end points,
lines = [];
tempLine = struct;
for p = 1:peakNum
	t1 = [-cos(pthe(p)),sin(pthe(p))]; % the unit directional vector of the line
	dist1 = t1*fitPts{p}(:,1:fitPtsNum(p)); % dist along the line
	[ordDist ordIdx] = sort(dist1,'ascend');
	ordPts = fitPts{p}(:,ordIdx);

	% handle with the gaps
	dist2 = diff(ordDist); % dist with each other
	gapPos = [0 find(dist2>maxGap) length(ordDist)];
	lineLen = diff(gapPos);
	for q = find(lineLen>=minLen)
		tempLine.theta = pthe(p);
		tempLine.rho = prho(p);
		tempLine.point1 = [M-ordPts(2,gapPos(q)+1),ordPts(1,gapPos(q)+1)+1];
		tempLine.point2 = [M-ordPts(2,gapPos(q+1)),ordPts(1,gapPos(q+1))+1];
		lines = [lines,tempLine];
		if length(lines) == maxLineNum, break; end
	end
	if length(lines) == maxLineNum, break; end
end

end

function [rows cols] = findpeaks(map,minLen,maxPeakNum,resol)
%	Find the at most maxLinNum points that's no less than minLen in map, meanwhile not
%	8-adjacent to each other, in MAP. Rows and cols are returned.

if max(max(map)) < minLen, return; end
rows = zeros(1,maxPeakNum);
cols = rows;
sz = size(map);
sup = ceil(sz/resol);
for p = 1:maxPeakNum
	[V,I] = max(map(:));
	if V<minLen
		rows = rows(1:p-1);
		cols = cols(1:p-1);
		break;
	end
	[rows(p),cols(p)] = ind2sub(sz,I);

	% non-maximal suppression
	left = max(1,cols(p)-sup(2));
	lr = min(sup(2)*2+1,sz(2)-left);
	up = max(1,rows(p)-sup(1));
	ud = min(sup(1)*2+1,sz(1)-up);
	map(up:up+ud,left:left+lr) = 0;
end

% if max(max(map)) < minLen, return; end
% map(map<minLen) = 0;
% mapMax = ordfilt2(map,25,ones(5));
% mapPeak = (mapMax == map) & map; % non-maximal suppression
% pt = find(mapPeak);
% pv = map(pt);
% [pvo,idx] = sort(pv,'descend');
% pt = pt(idx);
% if length(pvo) > maxPeakNum, pt = pt(1:maxPeakNum); end
% [rows cols] = ind2sub(size(map),pt);

end

author: Sable

repo: mcbench-benchmarks


text:
function A = metric_03 ( p )

%*****************************************************************************80
%
%% METRIC_03 evaluates metric #3 at any point.
%
%  Discussion:
%
%    This routine evaluates the matrix that determines the metric
%    at a point.
%
%    This particular matrix exaggerates distances in the Y direction.
%
%    It is diagonal, and it is spatially constant.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 May 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P(2), the point at which the metric matrix is to
%    be evaluated.
%
%    Output, real A[2,2], the metric matrix.
%
  A = [ 1.0, 0.0; 0.0, 100.0 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function maximise( varargin )

%MAXIMISE Specifiies a concave (or affine) objective to be maximized.

global cvx___
prob = evalin( 'caller', 'cvx_problem', '[]' );
if ~isa( prob, 'cvxprob' ),
    error( 'No CVX model exists in this scope.' );
elseif isempty( cvx___.problems ) || cvx___.problems( end ).self ~= prob,
    error( 'Internal CVX data corruption. Please CLEAR ALL and rebuild your model.' );
elseif nargin < 1,
    error( 'Objective expression missing.' );
elseif iscellstr( varargin ),
    x = evalin( 'caller', sprintf( '%s ', varargin{:} ) );
elseif nargin > 1,
    error( 'Too many input arguments.' );
else
    x = varargin{1};
end
try
    newobj( prob, 'maximize', x );
catch exc
    rethrow( exc )
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
%  OPTIMZE_LOGISTIC -- finds the optimal parameters of a logistic function
%
%  ::
%
%
%    ab=optimize_logistic(x1,p1,x2,p2)
%
%    ab=optimize_logistic(x1,p1,x2,p2,op)
%
%  Args:
%
%     - **x1,x2** [scalars] : points at which to evaluate the logistic
%
%     - **p1,p2** [scalars] : probability values for x1 and x2
%
%     - **op** [{'+'}|'-'] : decides whether b is to be multiplied by 1 or by
%       -1 in the evaluation the logistic function.
%
%  Returns:
%     :
%
%     - **ab** [vector] : optimized parameters of the logistic function
%
%  Note:
%
%       The function has the form f(x,a,b)=a/(a+exp(b*x)) so that f(x1,a,b)=p1
%       and f(x2,a,b)=p2
%
%  Example:
%
%     See also:
%

author: jmaih


text:
function [summary] = tapas_sem_posterior_summary(posterior, time)
%% Generate a summary of the posterior estimates.
%
% Input
%       posterior       -- Estimate from the model.
% Output
%       summary         -- Summary structue computed from the posterior.
%

% aponteeduardo@gmail.com
% copyright (C) 2019
%

n = 1;

n = n + 1;
if nargin < n
    time = linspace(0, 10.0, 1200);
end

% Plot the data.
data = posterior.data;
try
    model = posterior.model.graph{1}.htheta.model;
catch err
    try
        model = posterior.ptheta;
    catch
        error('tapas:sem:display_posterior', 'Improper input')
    end
end

ns = numel(data);

summary = struct('summaries', [], 'fits', []);

summaries = cell(ns, 1);
fits = cell(ns, 1);

for i = 1:ns
    samples = posterior.ps_theta(i, :);
    fits{i} = tapas_sem_generate_fits(data(i), samples, model, time);
    summaries{i} = tapas_sem_generate_summary(data(i), samples, model, i);
end

summary.fits = fits;
summary.summaries = vertcat(summaries{:});

end

author: translationalneuromodeling

repo: tapas


text:
function bessel_k0_values_test ( )

%*****************************************************************************80
%
%% BESSEL_K0_VALUES_TEST demonstrates the use of BESSEL_K0_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BESSEL_K0_VALUES_TEST:\n' );
  fprintf ( 1, '  BESSEL_K0_VALUES stores values of \n' );
  fprintf ( 1, '  the Bessel K0 function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X            K0(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = bessel_k0_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end


text:
function UpdatedDA=UpdateDA(CA,DA,Q,W)
% Update DA
% CA is the Archive that has  been updated.
% DA is the Archive that has not been updated
% Q  is the set of offspring
% W  is the set of weight vectors

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    S=[];                                                       % S is the set used for output
    Hd=[DA,Q];
    N=size(W,1);
    [~,Region_Hd] = max(1-pdist2(real(Hd.objs),W,'cosine'),[],2);
    [~,Region_CA] = max(1-pdist2(CA.objs,W,'cosine'),[],2);     % associat the individuals in Hd and CA with their corresponding subregions
    itr=1;
    while length(S)<N
        for i=1:N                                               % here i denotes the order of the subregion
            current_c=find(Region_CA==i);                       % current_c denotes that the current_c_th individual is/are already in the ith region
            if length(current_c)<itr
                for j=1:itr-length(current_c)                   % j denotes the number of solutions from Hd that need to join into the region(i)
                    current_d=find(Region_Hd==i);
                    if isempty(current_d)~=1
                        [FrontNO,~]=NDSort(Hd(current_d).objs,inf);
                        O=Hd(current_d(FrontNO==1));            % O is the set of nondominated solutions from region(i) in Hd
                        [~,Region_O] = max(1-pdist2(real(O.objs),W,'cosine'),[],2);
                        Z = min(O.objs,[],1);
                        g_tch=max(abs(O.objs-repmat(Z,length(O),1))./W(Region_O,:),[],2);
                        [~,order]=min(g_tch);
                        x_best=O(order);

                        Hd(current_d(Hd(current_d)==x_best))=[];% update Region_Hd
                        if isempty(Hd)==1
                            Region_Hd=[];
                        else
                            [~,Region_Hd] = max(1-pdist2(real(Hd.objs),W,'cosine'),[],2);
                        end
                        if length(S)<N                          % add the best individual into S
                            S=cat(2,S,x_best);
                        end
                    else
                        break;
                    end
               end
            end
            if length(S)==N
                break;
            end
        end
        itr=itr+1;
    end
    UpdatedDA=S;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [t,x_new,f_new,g_new,funEvals,H] = ArmijoBacktrack(...
    x,t,d,f,fr,g,gtd,c1,LS_interp,LS_multi,progTol,debug,doPlot,saveHessianComp,funObj,options,varargin)
% [t,x_new,f_new,g_new,funEvals,H] = ArmijoBacktrack(...
%    x,t,d,f,fr,g,gtd,c1,LS_interp,LS_multi,progTol,debug,doPlot,saveHessianComp,funObj,varargin)
%
% Backtracking linesearch to satisfy Armijo condition
%
% Inputs:
%   x: starting location
%   t: initial step size
%   d: descent direction
%   f: function value at starting location
%   fr: reference function value (usually funObj(x))
%   gtd: directional derivative at starting location
%   c1: sufficient decrease parameter
%   debug: display debugging information
%   LS_interp: type of interpolation
%   progTol: minimum allowable step length
%   doPlot: do a graphical display of interpolation
%   funObj: objective function
%   varargin: parameters of objective function
%
% Outputs:
%   t: step length
%   f_new: function value at x+t*d
%   g_new: gradient value at x+t*d
%   funEvals: number function evaluations performed by line search
%   H: Hessian at initial guess (only computed if requested)
%
% recet change: LS changed to LS_interp and LS_multi

% Evaluate the Objective and Gradient at the Initial Step
if nargout == 6
    [f_new,g_new,H] = funObj(x + t*d,varargin{:});
else

  x_struct = vector2struct(x+t*d, options.TEMPLATE);
  [f_new,g_new_struct] = feval(funObj, x_struct,varargin{:});
  g_new = struct2vector(g_new_struct);

end
funEvals = 1;

while f_new > fr + c1*t*gtd || ~isLegal(f_new)
    temp = t;

    if LS_interp == 0 || ~isLegal(f_new)
        % Ignore value of new point
        if debug
            fprintf('Fixed BT\n');
        end
        t = 0.5*t;
    elseif LS_interp == 1 || ~isLegal(g_new)
        % Use function value at new point, but not its derivative
        if funEvals < 2 || LS_multi == 0 || ~isLegal(f_prev)
            % Backtracking w/ quadratic interpolation based on two points
            if debug
                fprintf('Quad BT\n');
            end
            t = polyinterp([0 f gtd; t f_new sqrt(-1)],doPlot,0,t);
        else
            % Backtracking w/ cubic interpolation based on three points
            if debug
                fprintf('Cubic BT\n');
            end
            t = polyinterp([0 f gtd; t f_new sqrt(-1); t_prev f_prev sqrt(-1)],doPlot,0,t);
        end
    else
        % Use function value and derivative at new point

        if funEvals < 2 || LS_multi == 0 || ~isLegal(f_prev)
            % Backtracking w/ cubic interpolation w/ derivative
            if debug
                fprintf('Grad-Cubic BT\n');
            end
            t = polyinterp([0 f gtd; t f_new g_new'*d],doPlot,0,t);
        elseif ~isLegal(g_prev)
            % Backtracking w/ quartic interpolation 3 points and derivative
            % of two
            if debug
                fprintf('Grad-Quartic BT\n');
            end
            t = polyinterp([0 f gtd; t f_new g_new'*d; t_prev f_prev sqrt(-1)],doPlot,0,t);
        else
            % Backtracking w/ quintic interpolation of 3 points and derivative
            % of two
            if debug
                fprintf('Grad-Quintic BT\n');
            end
            t = polyinterp([0 f gtd; t f_new g_new'*d; t_prev f_prev g_prev'*d],doPlot,0,t);
         end
    end

    % Adjust if change in t is too small/large
    if t < temp*1e-3
        if debug
            fprintf('Interpolated Value Too Small, Adjusting\n');
        end
        t = temp*1e-3;
    elseif t > temp*0.6
        if debug
            fprintf('Interpolated Value Too Large, Adjusting\n');
        end
        t = temp*0.6;
    end

    % Store old point if doing three-point interpolation
    if LS_multi
        f_prev = f_new;
        t_prev = temp;
        if LS_interp == 2
            g_prev = g_new;
        end
    end

    if ~saveHessianComp && nargout == 6
        [f_new,g_new,H] = funObj(x + t*d,varargin{:});
    else

      x_struct = vector2struct(x+t*d, options.TEMPLATE);
      [f_new,g_new_struct] = feval(funObj, x_struct,varargin{:});
      g_new = struct2vector(g_new_struct);

%         [f_new,g_new] = funObj(x + t*d,varargin{:});
    end
    funEvals = funEvals+1;

    % Check whether step size has become too small
    if max(abs(t*d)) <= progTol
        if debug
            fprintf('Backtracking Line Search Failed\n');
        end
        t = 0;
        f_new = f;
        g_new = g;
        break;
    end
end

% Evaluate Hessian at new point
if nargout == 6 && funEvals > 1 && saveHessianComp

  [f_new,g_new,H] = funObj(x + t*d,varargin{:});
  funEvals = funEvals+1;

end

x_new = x + t*d;

end

repo: CategoryShapes


text:
function [overallFitness,mostFit,f] = getOverallFitness(cbalColinPowerWeights,fitnessMatrix)
% function [overallFitness,mostFit,f] = getOverallFitness(cbalColinPowerWeights,fitnessMatrix)
%
%
% Tor Wager, 11/17/01

   	% Determine overall fitness
   	for i = 1:size(cbalColinPowerWeights,2)				% convert to z scores
      if not(sum(fitnessMatrix(i,:)) == 0)              % only for rows with > 0 weight
         if std(fitnessMatrix(i,:)) == 0
            % disp(['Warning: All fitness scores for row ' num2str(i) ' are the same. Sum ' num2str(sum(fitnessMatrix(i,:)))])
            zfitnessMatrix(i,:) = (fitnessMatrix(i,:) - mean(fitnessMatrix(i,:))) / 1;
         else
      		zfitnessMatrix(i,:) = (fitnessMatrix(i,:) - mean(fitnessMatrix(i,:))) / std(fitnessMatrix(i,:));
         end
      end
   	end
   zfitnessMatrix((cbalColinPowerWeights == 0),:) = 0;
   overallFitness = (cbalColinPowerWeights * zfitnessMatrix) / sum(cbalColinPowerWeights);  % weighted avg sum z-scores

   % determine the most fit organism overall
   mostFit = find(overallFitness == max(overallFitness));

   % avoid duplicates
   mostFit = mostFit(1);

   % return fitness scores for most fit
   f = fitnessMatrix(:,mostFit);

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function ziggurat_test ( )

%*****************************************************************************80
%
%% ZIGGURAT_TEST tests the ZIGGURAT library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2008
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ZIGGURAT_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the ZIGGURAT library.\n' );
%
%  Make sure that SEED controls the sequence and can restart it.
%
  ziggurat_test001 ( );
  ziggurat_test02 ( );
  ziggurat_test03 ( );
  ziggurat_test04 ( );
%
%  Measure the time it takes to generate a lot of variables.
%
  sample_num = 100000;

  ziggurat_test05 ( sample_num );
  ziggurat_test06 ( sample_num );
  ziggurat_test07 ( sample_num );
  ziggurat_test08 ( sample_num );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ZIGGURAT_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function y = all(x)

% ALL True if all elements of a vector are nonzero.
%
% For vectors, ALL(V) returns logical 1 (TRUE) if none of the elements
% of the vector are zero.  Otherwise it returns logical 0 (FALSE).  For
% matrices, ALL(X) operates on the columns of X, returning a row vector
% of logical 1's and 0's. For N-D arrays, ALL(X) operates on the first
% non-singleton dimension.

% Copyright (C) 2012, Donders Centre for Cognitive Neuroimaging, Nijmegen, NL
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if nargin>1
  error('this implementation is only supported with one input argument');
end

siz = size(x);
if numel(siz)>2
  error('this implementation is only supported with vector or matrix input');
end

if siz(1)==1 || siz(2)==1
  y = true;
  for i=1:prod(siz)
    if ~x(i)
      y = false;
      break
    end
  end

else
  y = true(1,siz(2));
  for j=1:siz(2)
    for i=1:siz(1)
      if ~x(i,j)
        y(1,j) = false;
        break
      end
    end
  end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%demoL2iPotts_Deconv
% Reconstruction of a blurred jump-sparse signal from incomplete measurements under
% Gaussian noise using the inverse L2-Potts functional

% load signal
groundTruth = loadPcwConst('sampleDec');
n = numel(groundTruth);

% create Gaussian kernel
K = convkernel('gaussian', 51, 6);
Afull = spconvmatrix(K, numel(groundTruth));

% select random measurements
idx = sort(randidx(n, 0.5)) ;
A = Afull(idx, :);

% create blurred and noisy signal (Gaussian noise)
fBlurry = A * groundTruth(:);
fNoisy = fBlurry + 0.05 * randn(size(fBlurry));

% reconstruction using the inverse L2-Potts problem
gamma = 0.03;
[u, dataError, nJumps, energy] = minL2iPotts(fNoisy, gamma, A);

% show result
showPotts(fNoisy, u, groundTruth, 'L^2-iPotts')

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
function [flowx,flowy] = specific_wind(x,y,nel)
%vertical_wind   Reference problem 3.1 convective wind
%   [flowx,flowy] = specific_wind(x,y,nel);
%   input
%          x          x coordinate vector
%          y          y coordinate vector
%          nel        number of elements
%   specifies constant vertical wind.
%   IFISS function: DJS; 5 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      flowx =  zeros(nel,1);  flowy =  ones(nel,1);
      return


text:
function imageOrigin = getImageOrigin(MRSIStruct, dimension)
    switch(lower(dimension))
        case 'x'
            imageOrigin = MRSIStruct.imageOrigin(1);
        case 'y'
            imageOrigin = MRSIStruct.imageOrigin(2);
        case 'z'
            imageOrigin = MRSIStruct.imageOrigin(3);
        otherwise
            error('voxSize not found for dimension %s', dimension);
    end
end

author: CIC-methods

repo: FID-A

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


text:
function NDFUN
%NDFUN Matrix operations on N-D matrices
%   NDFUN treats an N-D matrix of double precision values as a set of pages
%   of 2D matrices, and performs various matrix operations on those pages.
%   The BLAS and LAPACK routines compiled into MATLAB are used for all these
%   operations, so results will be very close, and are usually identical, to
%   "native MATLAB results".  Available commands are:
%
%   C = NDFUN('mult', A, B)
%   C = NDFUN('backslash', A, B)
%   C = NDFUN('inv', A)
%   C = NDFUN('eig', A)
%   [C, D] = NDFUN('eig', A)
%   C = NDFUN('version')
%
%   The two-argument commands perform operations equivalent to:
%       for i=1:N
%           C(:,:,i) = A(:,:,i) * B(:,:,i);
%       end
%   The one-argument command
%       for i=1:N
%           C(:,:,i) = inv(A(:,:,i));
%       end
%
%   Any number of dimensions is supported, but dimensions > 2 must match:
%       C = ndfun('mult', rand(4,3,2,2,2), rand(3,1,2,2,2))
%   C will have size = [4 1 2 2 2]
%
%   NDFUN will reuse 2D arguments when needed, much like scalar
%   operations.  A single 2D matrix can be multiplied (or solved with)
%   each 2D page of the other argument.  For instance:
%       A = rand(4,3);  B = rand(3,10,100);
%       C = ndfun('mult', A, B);
%   is equivalent to:
%       for i=1:100
%           C(:,:,i) = A * B(:,:,i);
%       end
%   The reverse also works.  These types of operations are especially
%   efficient for the backslash operator.
%
%
%   Author: Peter Boettcher <boettcher@ll.mit.edu>
%   Source: www.mit.edu/~pwb/matlab/ndfun/
%
% Last modified: <Mon Nov 25 14:31:39 2002 by pwb>
%
% NOTES:
%   1) The option 'backslash' has given a segmentation fault.
%   2) Avoid NaNs. NaNs are not handled well in general.
%   3) Unpredictable behavior: the options 'inv' and 'eig' will
%       sometimes not give an answer. Chances are better if output
%       variables are one capital eltter and do not exist (?!)
%

  error('MEX file not found.  Try ''mex ndfun.c'' to compile.');

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function izakt
%IZAKT	Unit test for the function izak.

%	O. Lemoine - February 1996.

N=256;

% Perfect reconstruction
sig=noisecg(N);
DZT=zak(sig);
sigr=izak(DZT);
errors=find(any(abs(sig-sigr)>sqrt(eps)));
if length(errors)~=0,
  error('izak test 1 failed');
end

sig=noisecg(N);
DZT=zak(sig,8,32);
sigr=izak(DZT);
errors=find(any(abs(sig-sigr)>sqrt(eps)));
if length(errors)~=0,
  error('izak test 2 failed');
end

sig=noisecg(N);
DZT=zak(sig,128,2);
sigr=izak(DZT);
errors=find(any(abs(sig-sigr)>sqrt(eps)));
if length(errors)~=0,
  error('izak test 3 failed');
end

N=315;

% Perfect reconstruction
sig=noisecg(N);
DZT=zak(sig);
sigr=izak(DZT);
errors=find(any(abs(sig-sigr)>sqrt(eps)));
if length(errors)~=0,
  error('izak test 4 failed');
end

sig=noisecg(N);
DZT=zak(sig,9,35);
sigr=izak(DZT);
errors=find(any(abs(sig-sigr)>sqrt(eps)));
if length(errors)~=0,
  error('izak test 5 failed');
end

sig=noisecg(N);
DZT=zak(sig,3,105);
sigr=izak(DZT);
errors=find(any(abs(sig-sigr)>sqrt(eps)));
if length(errors)~=0,
  error('izak test 6 failed');
end

author: HeLiangHIT


text:
%--------------------------------------------------------------------------------------------------------
% The system is created based on the principles described in the following paper
% Jimmy SJ. Ren, Li Xu, Qiong Yan, Wenxiu Sun, "Shepard Convolutional Neural Networks",
% Advances in Neural Information Processing Systems (NIPS 2015)
% email: jimmy.sj.ren@gmail.com
%--------------------------------------------------------------------------------------------------------
addpath applications/Shepard_CNN/Shepard_super_res/
addpath utils/
addpath cuda/
addpath mem/
addpath layers/
addpath layers_adapters/
addpath optimization/
addpath pipeline/
addpath data/

clearvars -global config;
clearvars -global mem;
clear gen_mask_patch_cat_idx_for_super_res;
clear;
global config mem;
shepard_sr_x4_configure();
init(0);

load('data/Shepard_CNN/Shepard_super_res/x4/val_1ch/val_1');

images_t = reshape(images_t, size(images_t,1), size(images_t,2), 1, size(images_t,3));
labels_t = reshape(labels_t, size(labels_t,1), size(labels_t,2), 1, size(labels_t,3));

perm = randperm(size(images_t, 4));
images_t = images_t(:,:,:,perm);
labels_t = labels_t(:,:,:,perm);
test_samples = config.NEW_MEM(images_t(:,:,:,1:1000));
test_labels = config.NEW_MEM(labels_t(:,:,:,1:1000));

mask = config.NEW_MEM([1 0 0 0;0 0 0 0;0 0 0 0;0 0 0 0]);
mask = repmat(mask, config.input_size(1)/4, config.input_size(2)/4, config.chs);
%mask = config.NEW_MEM([1 0; 0 0]);
%mask = repmat(mask, config.input_size(1)/2, config.input_size(2)/2, config.chs);
mask = repmat(mask, 1,1,1,config.batch_size);

count = 0;
cost_avg = 0;
epoc = 0;
points_seen = 0;
display_points = 5000;
save_points = 50000;
max_grad = 1;
fprintf('%s\n', datestr(now, 'dd-mm-yyyy HH:MM:SS FFF'));
for pass = 1:10
    for p = 1:50
        load(strcat('data/Shepard_CNN/Shepard_super_res/x4/train_1ch/patches_', num2str(p), '.mat'));

        images = reshape(images, size(images,1), size(images,2), 1, size(images,3));
        labels = reshape(labels, size(labels,1), size(labels,2), 1, size(labels,3));

        perm = randperm(20000);
        images = images(:,:,:,perm);
        labels = labels(:,:,:,perm);
        train_imgs = config.NEW_MEM(images);
        train_labels = config.NEW_MEM(labels);

        for i = 1:size(train_labels, 4) / config.batch_size
            points_seen = points_seen + config.batch_size;
            in = train_imgs(:,:,:,(i-1)*config.batch_size+1:i*config.batch_size);
            %in = train_labels(:,:,:,(i-1)*config.batch_size+1:i*config.batch_size);
            out = train_labels(:,:,:,(i-1)*config.batch_size+1:i*config.batch_size);
            out = out((size(in, 1) - config.output_size(1)) / 2 + 1:(size(in, 1) - config.output_size(1)) / 2 + config.output_size(1), ...
                      (size(in, 2) - config.output_size(2)) / 2 + 1:(size(in, 2) - config.output_size(2)) / 2 + config.output_size(2), :, :);

            % make the mask list
            mask_li = {};
            mask_li{1} = mask;
            % operate the training pipeline
            op_train_pipe_with_mask(in.*mask, mask_li, out);
            % update the weights
            config.UPDATE_WEIGHTS();

            if(cost_avg == 0)
                cost_avg = config.cost;
            else
                cost_avg = (cost_avg + config.cost) / 2;
            end

            % display point
            if(mod(points_seen, display_points) == 0)
                count = count + 1;
                fprintf('%d ', count);
            end
            % save point
            if(mod(points_seen, save_points) == 0)
                fprintf('\n%s', datestr(now, 'dd-mm-yyyy HH:MM:SS FFF'));
                epoc = epoc + 1;
                test_cost = 0;
                for t = 1:size(test_samples, 4) / config.batch_size
                    t_label = test_labels(:,:,:,(t-1)*config.batch_size+1:t*config.batch_size);
                    t_label = config.NEW_MEM(t_label((size(in, 1) - config.output_size(1)) / 2 + 1:(size(in, 1) - config.output_size(1)) / 2 + config.output_size(1), ...
                                            (size(in, 2) - config.output_size(2)) / 2 + 1:(size(in, 2) - config.output_size(2)) / 2 + config.output_size(2), :));

                    op_test_pipe_with_mask(test_samples(:,:,:,(t-1)*config.batch_size+1:t*config.batch_size).*mask, mask_li, t_label);

                    test_out = gather(mem.output);
                    test_cost = test_cost + config.cost;
                end
                test_cost = test_cost / (size(test_samples, 4) / config.batch_size);
                fprintf('\nepoc %d, training avg cost: %f, test avg cost: %f\n', epoc, cost_avg, test_cost);

                save_weights(strcat('applications/Shepard_CNN/Shepard_super_res/results/shepard_layer_x4/epoc', num2str(epoc), '.mat'));

                cost_avg = 0;
            end
        end
    end
end

author: jimmy-ren


text:
function [varargout]=padLinDim(varargin)

% function [XP,indOriginal]=padLinDim(X,numPad,padDim)
% -----------------------------------------------------------------------
% Pad arrays (including vectors) allong dimension padDim with numPad
% entries which are linearly extrapolated allong that direction.
%
% -----------------------------------------------------------------------

%% Parse input

switch nargin
    case 1
        X=varargin{1};
        numPad=1;
        padDim=1;
        dirOpt='both';
    case 2
        X=varargin{1};
        numPad=varargin{2};
        padDim=1;
        dirOpt='both';
    case 3
        X=varargin{1};
        numPad=varargin{2};
        padDim=varargin{3};
        dirOpt='both';
    case 4
        X=varargin{1};
        numPad=varargin{2};
        padDim=varargin{3};
        dirOpt=varargin{4};
end

switch dirOpt
    case 'both'
        [XP,indOriginal]=padLinDimStart(X,numPad,padDim);
        XP=flip(XP,padDim);
        XP=padLinDimStart(XP,numPad,padDim);
        XP=flip(XP,padDim);
    case 'start'
        [XP,indOriginal]=padLinDimStart(X,numPad,padDim);
    case 'end'
        X=flip(X,padDim);
        [XP,indOriginal]=padLinDimStart(X,numPad,padDim);
        XP=flip(XP,padDim);
end

varargout{1}=XP;
varargout{2}=indOriginal;

end

%%

function [varargout]=padLinDimStart(X,numPad,padDim)
numDims=ndims(X);
siz=size(X);

Q=1:1:numDims;
Q=Q(Q~=padDim);
n=prod(siz(Q));
A=zeros(n,numDims);
A(:,padDim)=1;
for q=Q(:)'
    s=siz(q);
    A(:,q)=repmat(1:1:s,1,n/s);
end
B=A;
B(:,padDim)=2;

indDiff_1=sub2indn(siz,A);
indDiff_2=sub2indn(siz,B);

sizD=siz;
sizD(padDim)=1;

X1=reshape(X(indDiff_1),sizD);
X2=reshape(X(indDiff_2),sizD);
dX=X2-X1;

sizn=siz;
sizn(padDim)=sizn(padDim)+numPad;

IND=1:1:prod(sizn);
A=ind2subn(sizn,IND);
logicOld=(A(:,padDim)-numPad)>0;
indOriginal=IND(logicOld);

XP=zeros(sizn);
XP(indOriginal)=X;

A=zeros(n,numDims);
for qq=Q(:)'
    s=sizn(qq);
    A(:,qq)=repmat(1:1:s,1,n/s);
end

for q=1:1:numPad
    A(:,padDim)=q;
    indSet=sub2indn(sizn,A);
    XP(indSet)=((numPad-q)+1);
    XP(indSet)=X1-(((numPad-q)+1).*dX);
end

varargout{1}=XP;
varargout{2}=indOriginal;
end
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% Compute Hoyer sparsity of x

function spx = sp(x,w)

r = length(x);

spx = 0;
for i = 1 : r
    if x{i} == 0
        spx = 1;
    else
        ni = length(x{i});
        if nargin <= 1
            spx = spx + (sqrt(ni)-norm(x{i},1)/norm(x{i},2))/(sqrt(ni)-1);
        else
            nw = norm(w{i},2);
            spx = spx + (nw-w{i}'*abs(x{i})/norm(x{i},2))/(nw-min(w{i}));
        end
    end
end
spx = spx/r;

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [bool, storageType, remotePath, fileName] = isLocal(variable)
%"isLocal"
%   Returns false if passed variable has been stored outside of the planC.
%
% JRA 10/14/04
%
%Usage:
%   [bool, storageType, remotePath] = isLocal(variable);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

bool = 1;
if isstruct(variable) && isfield(variable, 'remotePath') & isfield(variable, 'info') && strcmpi(variable.info, 'Remotely stored data.  Use getRemoteVariable to access.')
    bool = 0;
    storageType = variable.storageType;
    remotePath  = variable.remotePath;
    fileName    = variable.filename;
else
    storageType = [];
    remotePath  = [];
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
neuron.A(:, IDs_last) = neuron_bk.A(:, IDs_last);
neuron.C(IDs_last,:) = neuron_bk.C(IDs_last,:);
neuron.C_raw(IDs_last,:) = neuron_bk.C_raw(IDs_last,:);
ind_del(IDs_last) = false;

axes(ax_merged); cla;
axes(ax_merged_trace); cla;

IDs = IDs_last;

author: zhoupc


text:
function [chansel, trlsel, cfg] = rejectvisual_summary(cfg, data)

% SUBFUNCTION for ft_rejectvisual

% Copyright (C) 2005-2006, Markus Bauer, Robert Oostenveld
% Copyright (C) 2006-2022, Robert Oostenveld, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if ft_platform_supports('matlabversion', '2014b', inf)
  % make the figure large enough to hold stuff and place it centered at the top of the
  % screen where it would appear by default
  pos = get(groot, 'DefaultFigurePosition'); % groot was introduced in R2014b
  w = 800;
  h = 600;
  x = pos(1) + pos(3)/2 - w/2;
  y = pos(2) + pos(4)/2 - h/2;
  cfg.position = ft_getopt(cfg, 'position', [x y w h]);
end

% determine the initial selection of trials
ntrl = length(data.trial);
if isequal(cfg.trials, 'all') % support specification like 'all'
  cfg.trials = 1:ntrl;
elseif isempty(cfg.trials)
  cfg.trials = 1:ntrl;
elseif isnumeric(cfg.trials)
  % use the selection as it is
end
trlsel = false(1, ntrl);
trlsel(cfg.trials) = true;

% determine the initial selection of channels
nchan = length(data.label);
cfg.channel = ft_channelselection(cfg.channel, data.label); % support specification like 'all'
chansel = false(1, nchan);
chansel(match_str(data.label, cfg.channel)) = true;

% compute the sampling frequency from the first two timepoints
fsample = 1/mean(diff(data.time{1}));

% % select the specified latency window from the data
% % here it is done BEFORE filtering and metric computation
% for i=1:ntrl
%   begsample = nearest(data.time{i}, cfg.latency(1));
%   endsample = nearest(data.time{i}, cfg.latency(2));
%   data.time{i} = data.time{i}(begsample:endsample);
%   data.trial{i} = data.trial{i}(:, begsample:endsample);
% end
if ischar(cfg.latency)
  cfg.latency(1) = min(cellfun(@min, data.time));
  cfg.latency(2) = max(cellfun(@max, data.time));
end

% compute the offset from the time axes
offset = zeros(ntrl, 1);
for i=1:ntrl
  offset(i) = time2offset(data.time{i}, fsample);
end

% set up guidata info
info                = [];
info.data           = data;
info.cfg            = cfg;
info.metric         = cfg.metric;
info.previousmetric = 'none';
info.level          = nan(nchan, ntrl);
info.ntrl           = ntrl;
info.nchan          = nchan;
info.trlsel         = trlsel;
info.chansel        = chansel;
info.fsample        = fsample;
info.offset         = offset;
info.quit           = 0;

if isfield(cfg, 'neighbours')
  % prepare the neighbours, load from disk, make channel selection, etc.
  tmpcfg = keepfields(cfg, {'neighbours'});
  info.cfg.neighbours = ft_prepare_neighbours(tmpcfg, data);
  % creates a NxN Boolean matrix that describes whether channels are connected as neighbours
  info.cfg.connectivity = channelconnectivity(info.cfg, data);
else
  info.cfg.connectivity = [];
end

h = create_figure(info);

% Compute initial metric...
compute_metric(h);

while ishandle(h)
  redraw(h);
  info = guidata(h);
  if info.quit == 0
    uiwait(h);
  else
    chansel = info.chansel;
    trlsel  = info.trlsel;
    cfg     = info.cfg;
    delete(h);
    break;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = create_figure(info)
% Creates the GUI and plots the data for the first time.
% All additional changes are done by manipulating the properties of the
% specific graphical objects
h = open_figure(keepfields(info.cfg, {'figure', 'position', 'visible', 'renderer', 'figurename', 'title'}));

% define three axes
info.axes(1) = axes('position', [0.100 0.650 0.375 0.300]);  % summary
info.axes(2) = axes('position', [0.575 0.650 0.375 0.300]);  % channels
info.axes(3) = axes('position', [0.100 0.250 0.375 0.300]);  % trials
% callback function (for toggling trials/channels) is set later, so that
% the user cannot try to toggle trials while nothing is present on the
% plots

% set up radio buttons for choosing metric
bgcolor = get(h, 'color');

metriclist = {'var' 'std' 'min' 'max' 'maxabs' 'range' 'kurtosis' '1/var' 'zvalue' 'maxzvalue' 'neighbexpvar' 'neighbcorr' 'neighbstdratio' 'db'};

uicontrol(h, 'Units', 'normalized', 'position', [0.575 0.50 0.24 0.04], 'Style', 'text', 'HorizontalAlignment', 'left', 'backgroundcolor', get(h, 'color'), 'string', 'Metric used:'); % text string for metric
uicontrol(h, 'Units', 'normalized', 'position', [0.725 0.51 0.20 0.04], 'Style', 'popupmenu', 'backgroundcolor', bgcolor, 'string', metriclist,       'HandleVisibility', 'off', 'callback', @change_metric); % popup menu for metric
uicontrol(h, 'Units', 'normalized', 'position', [0.575 0.46 0.24 0.04], 'Style', 'text', 'HorizontalAlignment', 'left', 'backgroundcolor', get(h, 'color'), 'string', 'Plot trial:');
info.plottrltxt = uicontrol(h, 'Units', 'normalized', 'position', [0.725 0.47 0.20 0.04], 'Style', 'edit', 'HorizontalAlignment', 'left', 'backgroundcolor', [1 1 1], 'callback', @plot_trials); % editbox for trial plotting

% editboxes for manually specifying which channels/trials to toggle on or off
uicontrol(h, 'Units', 'normalized', 'position', [0.575 0.42 0.24 0.04], 'Style', 'text', 'HorizontalAlignment', 'left', 'backgroundcolor', get(h, 'color'), 'string', 'Toggle trial:');
uicontrol(h, 'Units', 'normalized', 'position', [0.725 0.43 0.20 0.04], 'Style', 'edit', 'HorizontalAlignment', 'left', 'backgroundcolor', [1 1 1], 'callback', @toggle_trials, 'tag', 'edit_toggle_trials');
uicontrol(h, 'Units', 'normalized', 'position', [0.575 0.38 0.24 0.04], 'Style', 'text', 'HorizontalAlignment', 'left', 'backgroundcolor', get(h, 'color'), 'string', 'Toggle channel:');
uicontrol(h, 'Units', 'normalized', 'position', [0.725 0.39 0.20 0.04], 'Style', 'edit', 'HorizontalAlignment', 'left', 'backgroundcolor', [1 1 1], 'callback', @toggle_channels, 'tag', 'edit_toggle_channels');

info.excludetrllbl  = uicontrol(h, 'Units', 'normalized', 'position', [0.575 0.34 0.24 0.04], 'Style', 'text', 'HorizontalAlignment', 'left', 'backgroundcolor', get(h, 'color'), 'string', sprintf('Rejected trials: %i/%i', sum(info.trlsel==0), info.ntrl));
info.excludetrltxt  = uicontrol(h, 'Units', 'normalized', 'position', [0.605 0.30 0.45 0.04], 'Style', 'text', 'HorizontalAlignment', 'left', 'backgroundcolor', get(h, 'color'), 'Enable', 'Inactive');
info.excludechanlbl = uicontrol(h, 'Units', 'normalized', 'position', [0.575 0.26 0.28 0.04], 'Style', 'text', 'HorizontalAlignment', 'left', 'backgroundcolor', get(h, 'color'), 'string', sprintf('Rejected channels: %i/%i', sum(info.chansel==0), info.nchan));
info.excludechantxt = uicontrol(h, 'Units', 'normalized', 'position', [0.605 0.22 0.45 0.04], 'Style', 'text', 'HorizontalAlignment', 'left', 'backgroundcolor', get(h, 'color'), 'Enable', 'Inactive');

% instructions
instructions = sprintf('Drag the mouse over the channels or trials you wish to exclude');
uicontrol(h, 'Units', 'normalized', 'position', [0.575 0.125 0.35 0.08], 'Style', 'text', 'HorizontalAlignment', 'left', 'backgroundcolor', get(h, 'color'), 'string', instructions, 'FontWeight', 'bold', 'ForegroundColor', 'r');

% logbox
info.output_box = uicontrol(h, 'Units', 'normalized', 'position', [0.00 0.00 1.00 0.1], 'Style', 'edit', 'HorizontalAlignment', 'left', 'Max', 3, 'Min', 1, 'Enable', 'inactive', 'FontName', get(0, 'FixedWidthFontName'), 'FontSize', 9, 'ForegroundColor', [0 0 0], 'BackgroundColor', [1 1 1]);

% quit button
uicontrol(h, 'Units', 'normalized', 'position', [0.825 0.11 0.10 0.04], 'string', 'quit', 'callback', @quit);

% store the updated info in the figure
guidata(h, info);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function compute_metric(h)
info = guidata(h);
tmp = info.level(info.chansel, info.trlsel);
if isequal(info.metric, info.previousmetric) && all(~isnan(tmp(:)))
  % there is no reason to recompute the metric
  return
end

if contains(info.metric, 'zvalue')
  % cellmean and cellstd (see FT_DENOISE_PCA) would work instead of for-loops, but they are too memory-intensive
  runsum = zeros(info.nchan, 1);
  runss  = zeros(info.nchan, 1);
  runnum = 0;
  for i=1:info.ntrl
    % not entirely sure whether info.data.time{i} is correct, so I am constructing it here on the fly
    dat = preproc(info.data.trial{i}, info.data.label, offset2time(info.offset(i), info.fsample, size(info.data.trial{i}, 2)), info.cfg.preproc);
    runsum = runsum + nansum(dat, 2);
    runss  = runss  + nansum(dat.^2, 2);
    runnum = runnum + sum(isfinite(dat), 2);
  end
  mval = runsum./runnum;
  sd   = sqrt(runss./runnum - (runsum./runnum).^2);
else
  mval = [];
  sd   = [];
end

update_log(info.output_box, sprintf('Computing %s...', info.metric));
ft_progress('init', info.cfg.feedback, sprintf('computing %s', info.metric));
level = nan(info.nchan, info.ntrl);
for i=1:info.ntrl
  ft_progress(i/info.ntrl, 'computing %s %d of %d\n', i, info.metric, info.ntrl);
  % not entirely sure whether info.data.time{i} is correct, so I am constructing it here on the fly
  dat = preproc(info.data.trial{i}, info.data.label, offset2time(info.offset(i), info.fsample, size(info.data.trial{i}, 2)), info.cfg.preproc);
  val = artifact_level(dat, info.metric, mval, sd, info.cfg.connectivity);
  if isvector(val)
    level(:,i) = val;
  else
    level(:,i) = nanmedian(val);
  end
end

ft_progress('close');
update_log(info.output_box, 'Done.');

info.level = level;
info.previousmetric = info.metric;
guidata(h, info);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%s%%%%
function redraw(h)
info  = guidata(h);
% work with a copy of the data
level = info.level;

[maxperchan, maxpertrl, maxperchan_all, maxpertrl_all] = set_maxper(level, info.chansel, info.trlsel, info.metric);

% make the three figures
if gcf~=h, figure(h); end
% datacursormode on;

set(h, 'CurrentAxes', info.axes(1))
cla(info.axes(1));
switch info.cfg.viewmode
  case {'remove', 'toggle'}
    tmp = level;
    tmp(~info.chansel, :) = nan;
    tmp(:, ~info.trlsel)  = nan;
    imagesc(tmp, 'AlphaData', ~isnan(tmp));
    caxis([min(tmp(:)) max(tmp(:))]);
  case 'hide'
    imagesc(level(info.chansel==1, info.trlsel==1));
    if ~all(info.trlsel)
      set(info.axes(1), 'Xtick', []);
    end
    if ~all(info.chansel)
      set(info.axes(1), 'Ytick', []);
    end
end % switch
axis ij;
% colorbar;
title(info.cfg.method);
ylabel('channel number');
xlabel('trial number');

set(h, 'CurrentAxes', info.axes(2))
cla(info.axes(2));
switch info.cfg.viewmode
  case 'remove'
    plot(maxperchan(info.chansel==1),     find(info.chansel==1), '.');
    xmin = min(maxperchan);
    xmax = max(maxperchan);
    ymin = 1;
    ymax = info.nchan;
  case 'toggle'
    plot(maxperchan_all(info.chansel==1), find(info.chansel==1), '.');
    hold on;
    plot(maxperchan_all(info.chansel==0), find(info.chansel==0), 'o');
    hold off;
    xmin = min(maxperchan_all);
    xmax = max(maxperchan_all);
    ymin = 1;
    ymax = info.nchan;
  case 'hide'
    xmin = min(maxperchan);
    xmax = max(maxperchan);
    ymin = 1;
    ymax = sum(info.chansel==1);
    plot(maxperchan(info.chansel==1), 1:ymax, '.');
    if ~all(info.chansel)
      set(info.axes(2), 'Ytick', []);
    end
end % switch
% don't try to rescale the axes if they are empty
if any(info.chansel) && any(info.trlsel)
  % ensure that the horizontal and vertical range increase, also when negative
  % see https://github.com/fieldtrip/fieldtrip/issues/1150
  axis(fixrange([0.8*xmin 1.2*xmax ymin-0.5 ymax+0.5]));
end
axis ij;
set(info.axes(2), 'ButtonDownFcn', @toggle_visual);  % needs to be here; call to axis resets this property
ylabel('channel number');

set(h, 'CurrentAxes', info.axes(3))
cla(info.axes(3));
switch info.cfg.viewmode
  case 'remove'
    plot(find(info.trlsel==1), maxpertrl(info.trlsel==1), '.');
    xmin = 1;
    xmax = info.ntrl;
    ymin = min(maxpertrl);
    ymax = max(maxpertrl);
  case 'toggle'
    plot(find(info.trlsel==1), maxpertrl_all(info.trlsel==1), '.');
    hold on;
    plot(find(info.trlsel==0), maxpertrl_all(info.trlsel==0), 'o');
    hold off;
    xmin = 1;
    xmax = info.ntrl;
    ymin = min(maxpertrl_all);
    ymax = max(maxpertrl_all);
  case 'hide'
    xmin = 1;
    xmax = sum(info.trlsel==1);
    ymin = min(maxpertrl);
    ymax = max(maxpertrl);
    plot(1:xmax, maxpertrl(info.trlsel==1), '.');
    if ~all(info.trlsel)
      set(info.axes(3), 'Xtick', []);
    end
end % switch
% don't try to rescale the axes if they are empty
if any(info.chansel) && any(info.trlsel)
  % ensure that the horizontal and vertical range increase, also when negative
  % see https://github.com/fieldtrip/fieldtrip/issues/1150
  axis(fixrange([xmin-0.5 xmax+0.5 0.8*ymin 1.2*ymax]));
end
set(info.axes(3), 'ButtonDownFcn', @toggle_visual);  % needs to be here; call to axis resets this property
xlabel('trial number');

% put excluded trials/channels in their respective edit boxes
set(info.excludechanlbl, 'string', sprintf('Channels to exclude: %i/%i', sum(info.chansel==0), info.nchan));
set(info.excludetrllbl, 'string', sprintf('Trials to exclude: %i/%i', sum(info.trlsel==0), info.ntrl));

if ~all(info.trlsel)
  excludetrltxt = sprintf('%d, ', find(~info.trlsel));
  excludetrltxt = excludetrltxt(1:end-2);
  set(info.excludetrltxt, 'String', excludetrltxt, 'FontAngle', 'normal', 'ButtonDownFcn', @toggle_trials_fill);
else
  set(info.excludetrltxt, 'String', 'No trials to exclude', 'FontAngle', 'italic');
end

if ~all(info.chansel)
  if false % isfield(info.data, 'label')
    excludechantxt = sprintf('%s, ', info.data.label{~info.chansel});
  else
    excludechantxt = sprintf('%d, ', find(~info.chansel));
  end
  excludechantxt = excludechantxt(1:end-2);
  set(info.excludechantxt, 'String', excludechantxt, 'FontAngle', 'normal', 'ButtonDownFcn', @toggle_channels_fill);
else
  set(info.excludechantxt, 'String', 'No channels to exclude', 'FontAngle', 'italic');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function toggle_trials(h, eventdata)
info = guidata(h);

% process input from the "toggle trials" textbox
text = get(h, 'string');
set(h, 'string', '');
trls = handle_edit_input(text);

if isempty(trls)
  update_log(info.output_box, sprintf('Please enter one or more trials'));
  uiresume
  return
elseif any(trls<1)
  update_log(info.output_box, sprintf('ERROR: Trial value too small!'));
  uiresume
  return
elseif any(trls>length(info.trlsel))
  update_log(info.output_box, sprintf('ERROR: Trial value too large!'));
  uiresume
  return
else
  info.trlsel(trls) = ~info.trlsel(trls);
  % recalculate the metric
  compute_metric(h)
  guidata(h, info);
end
uiresume

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function toggle_trials_fill(h, eventdata)
info = guidata(h);
excludetrltxt = sprintf('%d, ', find(~info.trlsel));
excludetrltxt = excludetrltxt(1:end-2);
set(findobj('tag', 'edit_toggle_trials'), 'string', excludetrltxt)
set(info.plottrltxt, 'string', excludetrltxt)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function toggle_channels(h, eventdata)
info = guidata(h);

% process input from the "toggle channels" textbox
text = get(h, 'string');
set(h, 'string', '');
chans = handle_edit_input(text);                      % this can be 'all' or 'megref'
chans = ft_channelselection(chans, info.data.label);  % this is a cell-array of strings
chans = match_str(info.data.label, chans);            % this is a list of numbers

if isempty(chans)
  update_log(info.output_box, sprintf('Please enter one or more channels'));
  uiresume
  return
elseif any(chans<1)
  update_log(info.output_box, sprintf('ERROR: Channel value too small!'));
  uiresume
  return
elseif any(chans>length(info.data.label))
  update_log(info.output_box, sprintf('ERROR: Channel value too large!'));
  uiresume
  return
else
  info.chansel(chans) = ~info.chansel(chans);
  % recalculate the metric
  compute_metric(h)
  guidata(h, info);
end
uiresume

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function toggle_channels_fill(h, eventdata)
info = guidata(h);
excludechantxt = sprintf('%d, ', find(~info.chansel));
excludechantxt = excludechantxt(1:end-2);
set(findobj('tag', 'edit_toggle_channels'), 'string', excludechantxt)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function toggle_visual(h, eventdata)
% copied from FT_SELECT_BOX, but without the waitforbuttonpress command since here it is triggered by the ButtonDown event
point1 = get(gca, 'CurrentPoint');    % button down detected
rbbox;                                % this draws the rubber-band box
point2 = get(gca, 'CurrentPoint');    % button up detected
point1 = point1(1, 1:2);              % extract x and y
point2 = point2(1, 1:2);
x = sort([point1(1) point2(1)]);
y = sort([point1(2) point2(2)]);

g     = get(gca, 'Parent');
info  = guidata(g);

[maxperchan, maxpertrl, maxperchan_all, maxpertrl_all] = set_maxper(info.level, info.chansel, info.trlsel, info.metric);

switch gca
  case info.axes(1)
    % visual selection in the summary plot is not supported

  case info.axes(2)
    % the visual selection was made in the channels plot
    switch info.cfg.viewmode
      case 'toggle'
        chanlabels = 1:info.nchan;
        toggle = ...
          chanlabels >= y(1) & ...
          chanlabels <= y(2) & ...
          maxperchan_all(chanlabels)' >= x(1) & ...
          maxperchan_all(chanlabels)' <= x(2);
        info.chansel(toggle) = ~info.chansel(toggle);

      case 'remove'
        chanlabels     = 1:info.nchan;
        toggle = ...
          chanlabels >= y(1) & ...
          chanlabels <= y(2) & ...
          maxperchan(chanlabels)' >= x(1) & ...
          maxperchan(chanlabels)' <= x(2);
        info.chansel(toggle) = false;

      case 'hide'
        chanlabels = 1:sum(info.chansel==1);
        [junk, origchanlabels] = find(info.chansel==1);
        toggle = ...
          chanlabels >= y(1) & ...
          chanlabels <= y(2) & ...
          maxperchan(origchanlabels)' >= x(1) & ...
          maxperchan(origchanlabels)' <= x(2);
        info.chansel(origchanlabels(toggle)) = false;
    end

  case info.axes(3)
    % the visual selection was made in the trials plot
    switch info.cfg.viewmode
      case 'toggle'
        trllabels = 1:info.ntrl;
        toggle = ...
          trllabels >= x(1) & ...
          trllabels <= x(2) & ...
          maxpertrl_all(trllabels) >= y(1) & ...
          maxpertrl_all(trllabels) <= y(2);
        info.trlsel(toggle) = ~info.trlsel(toggle);

      case 'remove'
        trllabels = 1:info.ntrl;
        toggle = ...
          trllabels >= x(1) & ...
          trllabels <= x(2) & ...
          maxpertrl(trllabels) >= y(1) & ...
          maxpertrl(trllabels) <= y(2);
        info.trlsel(toggle) = false;

      case 'hide'
        trllabels = 1:sum(info.trlsel==1);
        [junk, origtrllabels] = find(info.trlsel==1);
        toggle = ...
          trllabels >= x(1) & ...
          trllabels <= x(2) & ...
          maxpertrl(origtrllabels) >= y(1) & ...
          maxpertrl(origtrllabels) <= y(2);
        info.trlsel(origtrllabels(toggle)) = false;
    end

end % switch gca

% recalculate the metric
compute_metric(h);
guidata(h, info);
uiresume

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function quit(h, eventdata)
info = guidata(h);
info.quit = 1;
guidata(h, info);
uiresume

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function change_metric(h, eventdata)
info = guidata(h);
%info.metric = get(eventdata.NewValue, 'string');
info.metric = h.String{h.Value};
guidata(h, info);
compute_metric(h);
uiresume

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function update_log(h, new_text)
new_text        = [datestr(now, 13) '# ' new_text];
curr_text       = get(h, 'string');
size_curr_text  = size(curr_text, 2);
size_new_text   = size(new_text, 2);
if size_curr_text > size_new_text
  new_text        = [new_text blanks(size_curr_text-size_new_text)];
else
  curr_text   = [curr_text repmat(blanks(size_new_text-size_curr_text), size(curr_text, 1), 1)];
end
set(h, 'String', [new_text; curr_text]);
drawnow

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [maxperchan, maxpertrl, maxperchan_all, maxpertrl_all] = set_maxper(level, chansel, trlsel, metric)
if ismember(metric, {'var', 'std', 'max', 'maxabs', 'range', 'kurtosis', '1/var', 'zvalue', 'maxzvalue', 'neighbstdratio'})
  % take the maximum
  extreme = @max;
elseif ismember(metric, {'min', 'neighbexpvar', 'neighbcorr'})
  % take the minimum
  extreme = @min;
end
% determine the extreme value
maxperchan_all = extreme(level, [], 2);
maxpertrl_all  = extreme(level, [], 1);
% determine the extreme value over the remaining selection
level(~chansel, :) = nan;
level(:, ~trlsel)  = nan;
maxperchan = extreme(level, [], 2);
maxpertrl  = extreme(level, [], 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_trials(h, eventdata)
info = guidata(h);

% process input from the "plot trials" textbox
text = get(h, 'string');
set(h, 'string', '');
trls = handle_edit_input(text);

if isempty(trls)
  update_log(info.output_box, sprintf('Please enter one or more trials'));
  uiresume
  return
elseif any(trls<1)
  update_log(info.output_box, sprintf('ERROR: Trial value too small!'));
  uiresume
  return
elseif any(trls>length(info.trlsel))
  update_log(info.output_box, sprintf('ERROR: Trial value too large!'));
  uiresume
  return
else
  update_log(info.output_box, 'Making multiplot of individual trials ...');
end

cfg_mp = [];
% disable hashing of input data (speeds up things)
cfg_mp.trackcallinfo = 'no';
cfg_mp.layout   = info.cfg.layout;
cfg_mp.channel  = info.data.label(info.chansel);
cfg_mp.dataname = info.cfg.dataname;
cfg_mp.ylim     = info.cfg.ylim;
currfig = gcf;
for n = 1:length(trls)
  % ft_multiplotER should be able to make the selection, but fails due to http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2978
  % that bug is hard to fix, hence it is solved here with a work-around
  cfg_sd = [];
  cfg_sd.trials = trls(n);
  tmpdata = ft_selectdata(cfg_sd, info.data);

  figure()
  cfg_mp.interactive = 'yes';
  ft_multiplotER(cfg_mp, tmpdata);
  title(sprintf('Trial %i', trls(n)));
end
figure(currfig);
update_log(info.output_box, 'Done.');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function range = fixrange(range)
% ensure that the horizontal and vertical range always increase
% also when both are negative, or both are the same, or both are or zero
% see https://github.com/fieldtrip/fieldtrip/issues/1150
if range(1)>range(2)
  % swap them
  range([1 2]) = range([2 1]);
elseif range(1)==range(2)
  if range(1)==0
    % move them a little bit apart from each other
    range([1 2]) = range([1 2]) + [-eps +eps];
  else
    % move them a little bit apart from each other, scaled with the value of interest
    range([1 2]) = range([1 2]) + [-eps +eps]*range(1);
  end
end
if range(3)>range(4)
  % swap them
  range([3 4]) = range([4 3]);
elseif range(3)==range(4)
  if range(3)==0
    % move them a little bit apart from each other
    range([3 4]) = range([3 4]) + [-eps +eps];
  else
    % move them a little bit apart from each other, scaled with the value of interest
    range([3 4]) = range([3 4]) + [-eps +eps]*range(3);
  end
end


text:
% efficiencyOFexpDesigns2 : OVERLAPPING EVENTS
% ER-fMRI data analysis
% script:
% throughout we assume TR=2.

loadDataYes=1; % =1 -> loads existing data file

pluginCorrYes=0;
daleYes=1; % 1-usual Dale efficiency. 0-Fisher

% event matrices are defined based on vectors with the average
% energy removed

cutEndPad=0; %cutting reduces efficiency!
nReps=10000; %n-2;
n2=24; % assumed # TRs to cover HDR fn
nVals=1;
pwRange=6:10; %7

if ~loadDataYes,

rndEffAll=[];
mEffAll=[];

tic
for pwr=pwRange,
% creating event vectors
ms=m2bin([mseq(2,pwr,0,1)])'; % fixing to the length of 2^n
n=length(ms)  % scan duration in TRs
if pluginCorrYes
  b=[0.406;0.8825]; % parameters from SPG fMRI data
  fittedACorr=autocorrFnct(b,1:n+(n2-1)*rem(cutEndPad+1,2));
  Cninv=inv(toeplitz(fittedACorr));
else
  Cninv=eye(n+(n2-1)*rem(cutEndPad+1,2));
end;

mEff=zeros(1,n-1);

shift=2^(pwr-1);
%shift1=2^(pwr-2)*3;
%XXXXXXXXXXXXXXXXXXXXX
for cycle=1:(2^pwRange(1))-2,
  mEvent=ms; %row 1
  mEvent=[mEvent; [ms(shift+1:end) ms(1:shift)]];%row 3
  mEvent=[mEvent; [ms(cycle+1:end) ms(1:cycle)]];%row 2

  % defining event matrix as convolution matrix
  %eventMatrix=makeEventMtrx(mEvent-(ones(length(mEvent),1)*...
  %sum(mEvent'))'/length(mEvent),n2);

  eventMatrix=makeEventMtrx(mEvent,n2);

  if cutEndPad, eventMatrix=eventMatrix(1:n,:); end;
  eventMatrix=eventMatrix-ones(size(eventMatrix,1),1)* ...
    sum(eventMatrix)/size(eventMatrix,1);

  if daleYes,
    designEff=1/trace(inv(eventMatrix'*Cninv*eventMatrix));
  else
    designEff=trace(eventMatrix'*Cninv*eventMatrix);
  end;

  mEff(cycle)=designEff;
end;

%mProbab=sum(mEvent'); mProbab=mean(mProbab)/n;
plot(mEff,'k','LineWidth',2);
set(gca,'FontSize',12);
xlabel('Cyclical shift of event vector #3','FontSize',16)
ylabel('Efficiency','FontSize',16)
%axis([0 n 0 1])
set(gca,'Position',[.2,.15,.7,.7]);
set(gcf,'PaperPosition',[1,1,5,4]);
%print -dpsc2 cyclingOverlapping3ev
drawnow;
pause

foo=max(mEff);
mEffAll=[mEffAll, foo(1)];

% for simply randomized designs:
nVals=size(mEvent,1);
nOnes=0.5;
overlapYes=1;

rndEff=[];
for k=1:nReps;
   %nOnes=k/(n-1);

	rndEvent=balancedRnd(n,nVals,nOnes,overlapYes);
	%if k==1, sum(rndEvent'), end;
        eventMatrix=makeEventMtrx(rndEvent-(ones(length(rndEvent),1)*sum(rndEvent'))'/length(rndEvent),n2);
   if cutEndPad, eventMatrix=eventMatrix(1:n,:); end;
        if daleYes,
	  designEff=1/trace(inv(eventMatrix'*Cninv*eventMatrix));
	else
	  designEff=trace(eventMatrix'*Cninv*eventMatrix);
	end;

	rndEff=[rndEff, designEff];
end;
rndEffAll=[rndEffAll;rndEff];
end;
toc;

else
  eval(['load dataEff',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nh']);
end %fi loadDataYes

nn=2.^[pwRange]-1;

figure(1); clf;
%[b,a]=hist(rndEff,20);
%bar(a,b/nReps); hold on;
%plot(max(mEff),0,'k*');

loglog(nn,mEffAll,'r*');hold on;

% theoretical max:
%theoMax=(nn+1)/n2/4;
%loglog(nn,theoMax,'g+');hold on;
p99=max(rndEffAll'); %prctile(rndEffAll',99.9);
p00=min(rndEffAll'); %prctile(rndEffAll',0.1);
med=median(rndEffAll');
loglog(nn,med,'k.-','LineWidth',2);
loglog([nn;nn],[p00;p99],'k')
loglog([nn*.93;nn*1.1],[p00;p00],'k')
loglog([nn*.93;nn*1.1],[p99;p99],'k')

axis([50 10^3*1.2 10^-.5 10^1*3]);
%set(gca,'YTick',[0.003 .01 .1 .3 1 3 10 30]);
set(gca,'XTick',[2^6 2^7 2^8 2^9 2^10]);
set(gca,'LineWidth',2,'FontSize',12);

xlabel('Sequence length','FontSize',12);
ylabel('Efficiency','FontSize',12);
title(['With Overlap: n_e=',num2str(nVals),' p=',num2str(nOnes),' n_h=',num2str(n2)]);

set(gca,'Position',[.2,.15,.7,.7]);
set(gcf,'PaperPosition',[1,1,4,3]);
eval(['print -dpsc2 eff',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nh_over']);
disp(['print -dpsc2 eff',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nh_over']);

figure(2); clf;
semilogx(nn,mEffAll./max(rndEffAll'),'-ok','LineWidth',2); hold on;
semilogx(nn,mEffAll./med,':+k','LineWidth',2)
set(gca,'Position',[.2,.15,.7,.7]);
xlabel('Sequence length','FontSize',12);
ylabel('m-seq/random efficiency ratio','FontSize',12);
title(['With overlap: n_e=',num2str(nVals),' p=',num2str(nOnes),' n_h=',num2str(n2)]);
set(gca,'Position',[.2,.15,.7,.7]);
set(gcf,'PaperPosition',[1,1,4,3]);
set(gca,'XTick',[2^6 2^7 2^8 2^9 2^10]);
set(gca,'LineWidth',2,'FontSize',12);
axis([50 10^3*1.2 .9 2.5]);
eval(['print -dpsc2 eff',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nhRATIO']);
disp(['print -dpsc2 eff',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nhRATIO']);

eval(['save dataEff',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nh rndEffAll p99 p00 med mEffAll']);
disp(['save dataEff',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nh rndEffAll p99 p00 med mEffAll']);


text:
% Function that optimizes Pose given an initial solution of the rotation only
%
% IMPORTANT: THIS FUNCTION ONLY DEALS WITH THE ORTHOGRAPHIC CASE RIGHT NOW
%
% This function optimizes the pose estimation of one camera with respect to
% another given only the 2 corresponding views.  The data is supposed to be
% centered and only the rotation has to be estimated. I it is here
% performed using simple gradient descent.
%
% The following code explains how the code to compute the
% error/gradient/hessian was created using matlab symbolic toolbox.
%
% Vincent's Structure From Motion Toolbox      Version 3.0
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

%%%%%  Infimum in TSFM   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms a b c cc d dd real

R = [ a^2+b^2-c^2-d^2 2*b*c-2*a*d 2*a*c+2*b*d; ...
 2*a*d+2*b*c a^2-b^2+c^2-d^2 2*c*d-2*a*b; ...
 2*b*d-2*a*c 2*a*b+2*c*d a^2-b^2-c^2+d^2 ]/(a^2+b^2+c^2+d^2);
R = R(1:2,:);
dRa = simple(diff(R,a));
dRb = simple(diff(R,b));
dRc = simple(diff(R,c));
dRd = simple(diff(R,d));
RtR = R'*R;
RtdRa = R'*simple(diff(R,a));
RtdRb = R'*simple(diff(R,b));
RtdRc = R'*simple(diff(R,c));
RtdRd = R'*simple(diff(R,d));

%%%%%  Rest of the code    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l = [ R(:); RtR(:) ]
l = [ R(:); dRa(:); dRb(:); dRc(:); dRd(:); RtR(:); RtdRa(:); RtdRb(:); RtdRc(:); RtdRd(:) ];
l = simple( l );

maple restart;
com = 'res := [ ';
k = 1;
for j = 1 : 6
  com = [ com 'r[' num2str(j) ']=' char(l(k)) ', '];
  k = k + 1;
end
kk=1
for i = 1 : 4
  for j = 1 : 6
	com = [ com 'dR[' num2str(kk) ']=' char(l(k)) ', '];
	k = k + 1;
	kk = kk + 1;
  end
end
for j = 1 : 9
  com = [ com 'rtR[' num2str(j) ']=' char(l(k)) ', '];
  k = k + 1;
end
kk = 1;
for i = 1 : 4
  for j = 1 : 9
    com = [ com 'rtDR[' num2str(kk) ']=' char(l(k)) ', '];
	k = k + 1;
	kk = kk + 1;
	end
end
com = [ com(1:end-2) '];' ];

%%%%%  Below is the code to generate C-code    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Mupad try, but it sucks ...
%  reset(symengine);
%  out = evalin(symengine, [ 'opts := generate::optimize([R=expr(' char(Rp') ') ]);' ] );
%
%  out = evalin(symengine, [ 'opts := generate::optimize([R=expr(' char(Rp') '), dR=expr(' char(l(:,3:10)) '), ddR=expr(' char(l(:,11:end)) ')]);' ] );
%  evalin(symengine, [ ':=rhs(opts[-1]);' ] );
%
%
%
%  out = [ char(evalin(symengine, [ 'generate::C(opts[1..-2]);' ] )) ...
%    char( evalin(symengine, [ 'generate::C(rhs(opts[-1]))' ] )) ];
%  out = strrep( out, '"', '' );
%  out = strrep( out, [ var '[0][' ], [ var '[' ] );

%%% now, save your crazy array into the res variable
maple( com );

%  %%% cost of the assigments with no optimization
%  maple('codegen[cost](res)')
%  % maple( 'opt1 := [codegen[optimize](l)]' );
%  %%% cost of the assigment with the normal optimization
%  in = [ 'l:=' char( l ) ]; maple( in ); maple('codegen[cost](codegen[optimize](l))')

%%% code for better optimization and corresponding cost
maple( 'opt2 := [codegen[optimize](res,tryhard)]' );
maple('codegen[cost](opt2)')

%%% Convert to C code
out = maple( 'codegen[C](opt2)' );

%%% clean the C code and convert it to Matlab code (replace [] by ()
out = strrep( out, ';   ', ';\n' );
out = strrep( out, '~', '' );

tout = ''
existing = cell(1,0);
tmp = regexp( out,'(t\d*)', 'tokens' )
for i = 1 : length(tmp)
  ii = tmp(i);
  ii = ii{1}{1};
  doExist = false;
  for j = 1 : size(existing,2)
    if strcmp(existing{j},ii)
	  doExist = true;
	  break;
	end
  end
  if ~doExist
    existing{1,end+1} = ii;
    tout = [ tout, ', ', ii ];
  end
end
out = [ 'double ', tout(2:end), ';\n', out ];
fprintf(out)

while 1
  tok = regexp( out,'\[(\d*)\]\[(\d*)\]', 'tokens' );
  tokStart = regexp( out,'\[(\d*)\]\[(\d*)\]', 'start' );
  tokEnd = regexp( out,'\[(\d*)\]\[(\d*)\]', 'end' );

  if isempty(tok); break; end

  out = [ out( 1:tokStart(1)-1 ) '(' num2str(str2double(tok{1}{1})+1) ...
      ',' num2str(str2double(tok{1}{2})+1) ')' out( tokEnd(1)+1 : end ) ];
end

fprintf( out );

author: vrabaud


text:
% Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
% Copyright (C) 2003 Doug Stewart <dastew@sympatico.ca>
% Copyright (C) 2011 Alexander Klein <alexander.klein@math.uni-giessen.de>
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, see <http://www.gnu.org/licenses/>.

% Generate a butterworth filter.
% Default is a discrete space (Z) filter.
%
% [b,a] = butter(n, Wc)
%    low pass filter with cutoff pi*Wc radians
%
% [b,a] = butter(n, Wc, 'high')
%    high pass filter with cutoff pi*Wc radians
%
% [b,a] = butter(n, [Wl, Wh])
%    band pass filter with edges pi*Wl and pi*Wh radians
%
% [b,a] = butter(n, [Wl, Wh], 'stop')
%    band reject filter with edges pi*Wl and pi*Wh radians
%
% [z,p,g] = butter(...)
%    return filter as zero-pole-gain rather than coefficients of the
%    numerator and denominator polynomials.
%
% [...] = butter(...,'s')
%     return a Laplace space filter, W can be larger than 1.
%
% [a,b,c,d] = butter(...)
%  return  state-space matrices
%
% References:
%
% Proakis & Manolakis (1992). Digital Signal Processing. New York:
% Macmillan Publishing Company.

function [a, b, c, d] = oc_butter (n, W, varargin)

if (nargin>4 || nargin<2) || (nargout>4 || nargout<2)
    print_usage;
end

% interpret the input parameters
if (~(length(n)==1 && n == round(n) && n > 0))
    error ('butter: filter order n must be a positive integer');
end

stop = 0;
digital = 1;
for i=1:length(varargin)
    switch varargin{i}
        case 's', digital = 0;
        case 'z', digital = 1;
        case { 'high', 'stop' }, stop = 1;
        case { 'low',  'pass' }, stop = 0;
        otherwise,  error ('butter: expected [high|stop] or [s|z]');
    end
end

[r, c]=size(W);
if (~(length(W)<=2 && (r==1 || c==1)))
    error ('butter: frequency must be given as w0 or [w0, w1]');
elseif (~(length(W)==1 || length(W) == 2))
    error ('butter: only one filter band allowed');
elseif (length(W)==2 && ~(W(1) < W(2)))
    error ('butter: first band edge must be smaller than second');
end

if ( digital && ~all(W >= 0 & W <= 1))
    error ('butter: critical frequencies must be in (0 1)');
elseif ( ~digital && ~all(W >= 0 ))
    error ('butter: critical frequencies must be in (0 inf)');
end

% Prewarp to the band edges to s plane
if digital
    T = 2;       % sampling frequency of 2 Hz
    W = 2/T*tan(pi*W/T);
end

% Generate splane poles for the prototype butterworth filter
% source: Kuc
C = 1; % default cutoff frequency
pole = C*exp(1i*pi*(2*(1:n) + n - 1)/(2*n));
if mod(n,2) == 1, pole((n+1)/2) = -1; end  % pure real value at exp(i*pi)
zero = [];
gain = C^n;

% splane frequency transform
[zero, pole, gain] = oc_sftrans(zero, pole, gain, W, stop);

% Use bilinear transform to convert poles to the z plane
if digital
    [zero, pole, gain] = oc_bilinear(zero, pole, gain, T);
end

% convert to the correct output form
if nargout==2,
    a = real(gain*poly(zero));
    b = real(poly(pole));
elseif nargout==3,
    a = zero;
    b = pole;
    c = gain;
else
    % output ss results
    % [a, b, c, d] = oc_zp2ss(zero, pole, gain);
    error('Not supported');
end

end

author: brainstorm-tools


text:
function yi = finterp1(x,y,xi,outOfRangeVal)
%"finterp1"
%   Fast 1-D linear interpolation of regularly spaced vector, to regularly
%   spaced vector/point.
%
% This function does the following
% transformation:	t = (xi-x(1))/(x(end)-x(1))
%                   yi = y(1) + t ( y(end) - y(1) )
%
%   x & y must be a monotonically increasing vector/array;
%
% Written DK 09-27-06
%
% Usage:
%       yi = finterp1(x,y,xi,outOfRangeVal);
%
% See also FINTERP2, FINTERP3, FINTERP3NOMESH ,INTERP1, INTERP2, INTERP3
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('outOfRangeVal')
    outOfRangeVal = NaN;
end

n = length(x);
% check if all the parameters passed are of the same class
superiorfloat(x,y,xi);

try
    h = diff(x);
    [ignore,k] = histc(xi,x);
    k(xi<x(1) | ~isfinite(xi)) = 1;
    k(xi>=x(n)) = n-1;
    t = (xi - x(k))./h(k);
catch % if the vector is monotonically decreasing
    t = (xi-x(1))/(x(end)-x(1));
    k = y(1) + t*(y(end) - y(1));
    k = round(k);
    k(xi<x(1) | ~isfinite(xi)) = 1;
    k(xi>=x(n)) = n-1;
end

yi = y(k)+t.*(y(k+1) - y(k));

% APA Check
% x = x(:);
% xi = xi(:);
% y = y(:);
% dxV = [diff(x); 1]; %DUMMY 1
% if all(dxV(1:end-1)<0)
%     x = flipud(x);
%     y = flipud(y);
%     dxV = [diff(x); 1];
% end
% dyV = [diff(y); 1]; %DUMMY 1
% indNaN = xi < min(x) | xi> max(x);
% xi(indNaN) = x(1); %assign DUMMY value
% [jnk,binIndex] = histc(xi,x);
% slopeV = dyV./dxV;
% slopeV(binIndex==length(x)) = 0;
% yi = y(binIndex) + slopeV(binIndex).*(xi-x(binIndex));
% yi(indNaN) = outOfRangeVal;
%
% return;