text:
function [predictedY,predictedYD,YD,stats] = predictPhenotype_CVHMM (Yin,Fin,T,...
    options_HMM,options_prediction,varargin)
%
% Kernel ridge regression or nearest-neighbour estimation using
% a distance matrix using (stratified) LOO.
% The difference with predictPhenotype.m is that the HMM is run also within
% the cross-validation scheme (run on training data, reapplied on test), so it's slower.
%
% Since the HMM is purely unsupervised (does not use Yin), it is
% safe to run it just once, and have only the prediction within the
% cross-validation loop (whether this is actually OK depends on the specific application).
% For this, use computeDistMatrix() and predictPhenotype()
%
% INPUT
% Yin       (no. subjects by no. phenotypes) matrix of phenotypic values to predict,
%           where each element can be continuous or binary. If a multiclass variable
%           is to be predicted, then Yin should be encoded by a
%           (no. subjects by no. classes) matrix, with zeros or ones
%           indicator entries.
% Fin       (no. subjects by 1) cell with (i) files of subject data, or (ii) just the data
% options   Struct with the prediction options, with fields:
%   + alpha - for method='KRR', a vector of weights on the L2 penalty on the regression
%           By default: [0.0001 0.001 0.01 0.1 0.4 0.7 1.0 10 100]
%   + sigmafact - for method='KRR', a vector of parameters for the kernel; in particular,
%           this is the factor by which we will multiply a data-driven estimation
%           of the kernel parameter. By default: [1/5 1/3 1/2 1 2 3 5];
%   + CVscheme - vector of two elements: first is number of folds for model evaluation;
%             second is number of folds for the model selection phase (0 in both for LOO)
%   + CVfolds - prespecified CV folds for the outer loop
%   + biascorrect - whether we correct for bias in the estimation
%                   (Smith et al. 2019, NeuroImage)
%   + verbose -  display progress?
% cs        optional (no. subjects X no. subjects) dependency structure matrix with
%           specifying possible relations between subjects (e.g., family
%           structure), or a (no. subjects X 1) vector defining some
%           grouping, with (1...no.groups) or 0 for no group
% confounds     (no. subjects by  no. of confounds) matrix of features that
%               potentially influence the phenotypes and we wish to control for
%               (optional)
%
% OUTPUT
% predictedY    predicted response,in the original (non-decounfounded) space
% predictedYD    predicted response,in the decounfounded space
% YD    response,in the decounfounded space
% stats         structure, with fields
%   + pval - permutation-based p-value, if permutation is run;
%            otherwise, correlation-based p-value
%   + cod - coeficient of determination
%   + corr - correlation between predicted and observed Y
%   + baseline_corr - baseline correlation between predicted and observed Y for null model
%   + sse - sum of squared errors
%   + baseline_sse - baseline sum of squared errors
%   PLUS: All of the above +'_deconf' in the deconfounded space, if counfounds were specified
%
% Author: Diego Vidaurre, OHBA, University of Oxford

options = options_prediction;
[N,q] = size(Yin);

if ~iscell(Fin), error('Argument Fin must be a cell'); end

which_nan = false(N,1);
if q == 1
    which_nan = isnan(Yin);
    if any(which_nan)
        Yin = Yin(~which_nan);
        Fin = Fin(~which_nan,:);
        warning('NaN found on Yin, will remove...')
    end
    N = size(Yin,1);
end

if isempty(options), options = struct(); end
if nargin<3, options = struct(); end

if ~isfield(options,'alpha')
    alpha = [0.0001 0.001 0.01 0.1 0.4 0.7 1.0 10 100];
else
    alpha = options.alpha;
end
if ~isfield(options,'sigmafact')
    sigmafact = [1/5 1/3 1/2 1 2 3 5];
else
    sigmafact = options.sigmafact;
end
if ~isfield(options,'K')
    K = 1:min(50,round(0.5*N));
else
    K = options.K;
end

if ~isfield(options,'CVscheme'), CVscheme = [10 10];
else, CVscheme = options.CVscheme; end
if ~isfield(options,'CVfolds'), CVfolds = [];
else, CVfolds = options.CVfolds; end
% if ~isfield(options,'biascorrect'), biascorrect = 0;
% else, biascorrect = options.biascorrect; end
if ~isfield(options,'verbose'), verbose = 1;
else, verbose = options.verbose; end

% check correlation structure
allcs = [];
if (nargin>=6) && ~isempty(varargin{1})
    cs = varargin{1};
    if ~isempty(cs)
        is_cs_matrix = (size(cs,2) == size(cs,1));
        if size(cs,2)>1 % matrix format
            if any(which_nan)
                cs = cs(~which_nan,~which_nan);
            end
            [allcs(:,2),allcs(:,1)]=ind2sub([length(cs) length(cs)],find(cs>0));
        else
            if any(which_nan)
                cs = cs(~which_nan);
            end
            allcs = find(cs > 0);
        end
    end
else, cs = [];
end

% get confounds
if (nargin>=7) && ~isempty(varargin{2})
    confounds = varargin{2};
    confounds = confounds - repmat(mean(confounds),N,1);
    deconfounding = 1;
    if any(which_nan)
        confounds = confounds(~which_nan,:);
    end
else
    confounds = []; deconfounding = 0;
end

Ymean = zeros(N,q);
YD = zeros(N,q); % deconfounded signal
YmeanD = zeros(N,q); % mean in deconfounded space
predictedY = zeros(N,q);
if deconfounding, predictedYD = zeros(N,q); end

% create the inner CV structure - we can't to stratified because it's the
% same fold structure for all variables
if isempty(CVfolds)
    if CVscheme(1)==1
        folds = {1:N};
    elseif q == 1
        Yin_copy = Yin; Yin_copy(isnan(Yin)) = realmax;
        folds = cvfolds(Yin_copy,CVscheme(1),allcs);
    else  % no stratification
        folds = cvfolds(randn(size(Yin,1),1),CVscheme(1),allcs);
    end
else
    folds = CVfolds;
end

for ifold = 1:length(folds)

    if verbose, fprintf('CV iteration %d \n',ifold); end

    J = folds{ifold}; % test
    if isempty(J), continue; end
    if length(folds)==1
        ji = J;
    else
        ji = setdiff(1:N,J); % train
    end

    hmm = hmmmar(Fin(ji),T(ji),options_HMM);
    HMMs_dualregr = cell(N,1);

    % dual-estimation
    for j = 1:N
        if ischar(Fin{j})
            fsub = Fin{j};
            loadfile_sub;
        else
            X = Fin{j};
        end
        HMMs_dualregr{j} = hmmdual(X,T{j},hmm);
        HMMs_dualregr{j}.state = rmfield(HMMs_dualregr{j}.state,'prior');
    end

    disp('HMM done')

    Din = zeros(N,N);
    parfor n1 = 1:N-1
        din = zeros(1,N);
        for n2 = n1+1:N
            % FO is contained in TPC; TPC is contained in HMM
             din(n2) = (hmm_kl(HMMs_dualregr{n1},HMMs_dualregr{n2}) ...
                + hmm_kl(HMMs_dualregr{n2},HMMs_dualregr{n1}))/2;
        end
        Din(n1,:) = din;
    end; Din = Din' + Din;

    D = Din(ji,ji);
    Y = Yin(ji,:);
    D2 = Din(J,ji);

    % family structure for this fold
    Qallcs=[];
    if (~isempty(cs))
        if is_cs_matrix
            [Qallcs(:,2),Qallcs(:,1)] = ...
                ind2sub([length(cs(ji,ji)) length(cs(ji,ji))],find(cs(ji,ji)>0));
        else
            Qallcs = find(cs(ji) > 0);
        end
    end

    parfor ii = 1:q

        Dii = D; Yii = Y;

        ind = find(~isnan(Y(:,ii)));
        Yii = Yii(ind,ii);
        QDin = Dii(ind,ind);
        QN = length(ind);

        Qfolds = cvfolds(Yii,CVscheme(2),Qallcs); % we stratify

        % deconfounding business
        if deconfounding
            Cii = confounds(ji,:); Cii = Cii(ind,:);
            [betaY,interceptY,Yii] = deconfoundPhen(Yii,Cii);
        end

        % centering response
        my = mean(Yii);
        Yii = Yii - repmat(my,size(Yii,1),1);
        Ymean(J,ii) = my;
        QYin = Yii;

        Dev = Inf(length(alpha),length(sigmafact));

        for isigm = 1:length(sigmafact)

            sigmf = sigmafact(isigm);

            QpredictedY = Inf(QN,length(alpha));
            QYinCOMPARE = QYin;

            % Inner CV loop
            for Qifold = 1:length(Qfolds)

                QJ = Qfolds{Qifold}; Qji = setdiff(1:QN,QJ);
                QD = QDin(Qji,Qji);
                QY = QYin(Qji); Qmy = mean(QY); QY = QY-Qmy;
                Nji = length(Qji);
                QD2 = QDin(QJ,Qji);

                sigmabase = auto_sigma(QD);
                sigma = sigmf * sigmabase;

                K = gauss_kernel(QD,sigma);
                K2 = gauss_kernel(QD2,sigma);
                I = eye(Nji);
                ridg_pen_scale = mean(diag(K));

                for ialph = 1:length(alpha)
                    alph = alpha(ialph);
                    beta = (K + ridg_pen_scale * alph * I) \ QY;
                    QpredictedY(QJ,ialph) = K2 * beta + repmat(Qmy,length(QJ),1);
                end
            end

            Dev(:,isigm) = (sum(( QpredictedY - ...
                repmat(QYinCOMPARE,1,length(alpha))).^2) / QN)';

        end

        [~,m] = min(Dev(:)); % Pick the one with the lowest deviance
        [ialph,isigm] = ind2sub(size(Dev),m);
        sigmf = sigmafact(isigm);
        sigmabase = auto_sigma(D);
        sigma = sigmf * sigmabase;
        alph = alpha(ialph);

        Dii = D(ind,ind); D2ii = D2(:,ind);

        K = gauss_kernel(Dii,sigma);
        K2 = gauss_kernel(D2ii,sigma);
        Nji = length(ind);
        I = eye(Nji);

        ridg_pen_scale = mean(diag(K));
        beta = (K + ridg_pen_scale * alph * I) \ Yii;

        % predict the test fold
        predictedY(J,ii) = K2 * beta + my; % some may be NaN actually

        % predictedYD and YD in deconfounded space; Yin and predictedY are confounded
        predictedYD(J,ii) = predictedY(J,ii);
        YD(J,ii) = Yin(J,ii);
        YmeanD(J,ii) = Ymean(J,ii);
        if deconfounding % in order to later estimate prediction accuracy in deconfounded space
            [~,~,YD(J,ii)] = deconfoundPhen(YD(J,ii),confounds(J,:),betaY,interceptY);
            % original space
            predictedY(J,ii) = confoundPhen(predictedY(J,ii),confounds(J,:),betaY,interceptY);
            Ymean(J,ii) = confoundPhen(YmeanD(J,ii),confounds(J,:),betaY,interceptY);
        end

%         if biascorrect % we do this in the original space
%             Yhattrain = K * beta + my;
%             if deconfounding
%                 Yhattrain = confoundPhen(Yhattrain,confounds(ji,:),betaY,interceptY);
%                 Ytrain = [confoundPhen(QYin,confounds(ji,:),betaY,interceptY) ...
%                     ones(size(QYin,1),1)];
%             else
%                 Ytrain = [QYin ones(size(QYin,1),1)];
%             end
%             b = pinv(Ytrain) * Yhattrain;
%             predictedY(J,ii) = (predictedY(J,ii) - b(2)) / b(1);
%         end

        if rem(ii,100)==0, disp(['Variable ' num2str(ii) ]); end

    end

    disp(['Fold ' num2str(ifold) ])

end

stats = struct();
stats.sse = zeros(q,1);
stats.cod = zeros(q,1);
stats.corr = zeros(q,1);
stats.baseline_corr = zeros(q,1);
stats.pval = zeros(q,1);
if deconfounding
    stats.sse_deconf = zeros(q,1);
    stats.cod_deconf = zeros(q,1);
    stats.corr_deconf = zeros(q,1);
    stats.baseline_corr_deconf = zeros(q,1);
    stats.pval_deconf = zeros(q,1);
end

for ii = 1:q
    ind = find(~isnan(Yin(:,ii)));
    stats.sse(ii) = sum((Yin(ind,ii)-predictedY(ind,ii)).^2);
    nullsse = sum((Yin(ind,ii)-Ymean(ind,ii)).^2);
    stats.cod(ii) = 1 - stats.sse(ii) / nullsse;
    stats.corr(ii) = corr(Yin(ind,ii),predictedY(ind,ii));
    stats.baseline_corr(ii) = corr(Yin(ind,ii),Ymean(ind,ii));
    [~,pv] = corrcoef(Yin(ind,ii),predictedY(ind,ii)); % original space
    if corr(Yin(ind,ii),predictedY(ind,ii))<0, stats.pval(ii) = 1;
    else, stats.pval(ii) = pv(1,2);
    end
    if deconfounding
        stats.sse_deconf(ii) = sum((YD(ind,ii)-predictedYD(ind,ii)).^2);
        nullsse_deconf = sum((YD(ind,ii)-YmeanD(ind,ii)).^2);
        stats.cod_deconf(ii) = 1 - stats.sse_deconf(ii) / nullsse_deconf;
        stats.corr_deconf(ii) = corr(YD(ind,ii),predictedYD(ind,ii));
        stats.baseline_corr_deconf(ii) = corr(YD(ind,ii),YmeanD(ind,ii));
        [~,pv] = corrcoef(YD(ind,ii),predictedYD(ind,ii)); % original space
        if corr(YD(ind,ii),predictedYD(ind,ii))<0, stats.pval_deconf(ii) = 1;
        else, stats.pval_deconf(ii) = pv(1,2);
        end
    end
end

end

function K = gauss_kernel(D,sigma)
% Gaussian kernel
D = D.^2; % because distance is sqrt-ed
K = exp(-D/(2*sigma^2));
end

function sigma = auto_sigma (D)
% gets a data-driven estimation of the kernel parameter
D = D(triu(true(size(D,1)),1));
sigma = median(D);
end

function [betaY,my,Y] = deconfoundPhen(Y,confX,betaY,my)
if nargin<3, betaY = []; end
if isempty(betaY)
    my = mean(Y);
    Y = Y - repmat(my,size(Y,1),1);
    betaY = (confX' * confX + 0.00001 * eye(size(confX,2))) \ confX' * Y;
end
Y = Y - confX*betaY;
end

function Y = confoundPhen(Y,conf,betaY,my)
Y = Y + conf*betaY + repmat(my,size(Y,1),1);
end

author: OHBA-analysis

repo: HMM-MAR


text:
function g = dftcorr(f, w)
%DFTCORR 2-D correlation in the frequency domain.
%   G = DFTCORR(F, W) performs the correlation of a mask, W, with
%   image F. The output, G, is the correlation image, of class
%   double. The output is of the same size as F. When, as is
%   generally true in practice, the mask image is much smaller than
%   G, wraparound error is negligible if W is padded to size(F).

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.4 $  $Date: 2003/10/26 23:20:44 $

[M, N] = size(f);
f = fft2(f);
w = conj(fft2(w, M, N));
g = real(ifft2(w.*f));

repo: weiyanmin


text:
function title = p01_title ( )

%*****************************************************************************80
%
%% P01_TITLE returns a title for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The Fletcher-Powell helical valley function.';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%kxvortho 'Visualize and navigate by orthogonal planes in a volume'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros xvortho.pane file
%
% Parameters:
% InputFile: i ' Input Image ', optional: 'Input image filename'
% InputFile: cmap 'Colormap Image', optional: 'input alternate colormap file'
% Integer: x 'X Placement', default: -1: 'x location for GUI autoplacement'
% Integer: y 'Y Placement', default: -1: 'y location for GUI autoplacement'
% Double: update 'Update Time', default: 0: 'initial update time'
%
% Example: kxvortho({i, cmap}, {'i','';'cmap','';'x',-1;'y',-1;'update',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% xvortho - Interactively navigate by orthogonal planes in a volume
%
%  DESCRIPTION
% "xvortho" is an interactive tool which allows navigation by any input
% object.
% The features provided by "xvortho" include:
% - An input file history list
% - Normalization of colormaps
% - Control over which columns of the map data are displayed as red, green, and blue
% - Application of functions to the map data values to produce new colormaps
% - Pre-defined colormap installation
%
% The master form of "xvortho" displays the input object in a workspace.
% The object is displayed as three orthogonal planes, named planes X, Y
% and Z, as the input object was a 3D volume.
% Plane X corresponds to a constant dimention w and it is a function of d and h.
% Plane Y corresponds to a constant dimention h and it is a function of w and d.
% Plane Z corresponds to a constant dimention d and it is a function of w and h.
%
% The initial configuration shows plane X with w=0, plane Y with h=0 and
% plane Z with d=0. Initial time and element dimentions are also zero. Clicking
% in one of the three displayed planes is equivalent to changing the choosen
% positions w, h or d, depending on which plane the cursor was pointing on
% during the clicking. For instance, if we are in the position w=10, h=20 and
% d=30, and the cursor is in the plane Z, on the
% position (50,100) when the button is pressed, then the new positions in w and
% h will be 50 and 100 respectively, which will cause the planes X and Y to be
% updated. Changing of time and element dimentions are also possible using
% scroll bars on the bottom of the workspace.
%
% Above the workspace are "License", "Help", and "Quit" buttons, from which one
% may access the license agreement, display online help pages, and quit the
% program, respectively.  To the left of the image workspace are two subform
% buttons labeled "Files" and "Options".
%
% Use the "Files" subform button to display the "Input and Output Files" subform.
% From this subform, you may input a new image. You may also output any of the
% three displayed
% images, either as they currently appears, or as they really are in the input
% image.
%
% A number of features are available from the "Options" subform, which
% contains three panes: "Image List", "Display" and "Colormaps". Those were
% adapted from editimage.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% editimage
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993, 1994, 1995 UNICAMP, R A Lotufo,  All rights reserved.
%

function varargout = kxvortho(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kxvortho(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'cmap', '__input';'x', -1;'y', -1;'update', 0};
maxval={1,1,1000,1000,1};
minval={1,1,-1,-1,1};
istoggle=[1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','Integer','Integer','Double'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'xvortho"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function setup_ilsvrcanimalpart()
%% download images
if(~exist('../data/detanimalpart','dir'))
    system('rm -r ../data/detanimalpart');
    system('git clone https://github.com/zqs1022/detanimalpart.git');
    system('mv detanimalpart ../data/');
end

%% settings
conf.data.imgdir='../data/detanimalpart/';
conf.data.readCode='./data_input/data_input_DET-2015-subset-ILSVRC2013train/';
conf.data.minArea=50^2;

conf.output.dir='./mat/';

Name_batch={'n01443537','n01503061','n01639765','n01662784','n01674464','n01882714','n01982650','n02084071','n02118333','n02121808','n02129165','n02129604','n02131653','n02324045','n02342885','n02355227','n02374451','n02391049','n02395003','n02398521','n02402425','n02411705','n02419796','n02437136','n02444819','n02454379','n02484322','n02503517','n02509815','n02510455'};
for i=1:numel(Name_batch)
    mkdir([conf.output.dir,Name_batch{i}]);
    save([conf.output.dir,Name_batch{i},'/conf.mat'],'conf');
end

addpath(genpath('./tool'));
addpath(conf.data.readCode);

repo: interpretableCNN


text:
function [samples, energies, diagn] = metrop2(f, x, opt, gradf, varargin)
%METROP2	Markov Chain Monte Carlo sampling with Metropolis algorithm.
%
%	Description
%	SAMPLES = METROP(F, X, OPT) uses the Metropolis algorithm to
%	sample from the distribution P ~ EXP(-F), where F is the first
%	argument to METROP.   The Markov chain starts at the point X and each
%	candidate state is picked from a Gaussian proposal distribution and
%	accepted or rejected according to the Metropolis criterion.
%
%	SAMPLES = METROP(F, X, OPT, [], P1, P2, ...) allows additional
%	arguments to be passed to F().  The fourth argument is ignored, but
%	is included for compatibility with HMC and the optimizers.
%
%	[SAMPLES, ENERGIES, DIAGN] = METROP(F, X, OPT) also returns a log
%	of the energy values (i.e. negative log probabilities) for the
%	samples in ENERGIES and DIAGN, a structure containing diagnostic
%	information (position and acceptance threshold) for each step of the
%	chain in DIAGN.POS and DIAGN.ACC respectively.  All candidate states
%	(including rejected ones) are stored in DIAGN.POS.
%
%	S = METROP('STATE') returns a state structure that contains the state
%	of the two random number generators RAND and RANDN. These are
%	contained in fields randstate,  randnstate.
%
%	METROP('STATE', S) resets the state to S.  If S is an integer, then
%	it is passed to RAND and RANDN. If S is a structure returned by
%	METROP('STATE') then it resets the generator to exactly the same
%	state.
%
%       See METROP2_OPT for the optional parameters in the OPTIONS
%       structure.
%
%	See also
%	HMC, METROP2_OPT
%

%	Copyright (c) Christopher M Bishop, Ian T Nabney (1996, 1997)
%	Copyright (c) 1998-2000 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% Global variable to store state of momentum variables: set by set_state
% Used to initialize variable if set
global HMC_MOM
if nargin <= 2
  if ~strcmp(f, 'state')
    error('Unknown argument to metrop2');
  end
  switch nargin
   case 1
    samples = get_state(f);
    return;
   case 2
    set_state(f, x);
    return;
  end
end

% Set empty omptions to default values
opt=metrop2_opt(opt);
% Refrence to structures is much slower, so...
opt_nsamples=opt.nsamples;
opt_display =opt.display;
stddev = opt.stddev;

% Set up string for evaluating potential function.
%f = fcnchk(f, length(varargin));

nparams = length(x);
samples = zeros(opt_nsamples, nparams);		% Matrix of returned samples.
if nargout >= 2
  en_save = 1;
  energies = zeros(opt_nsamples, 1);
else
  en_save = 0;
end
if nargout >= 3
  diagnostics = 1;
  diagn_pos = zeros(opt_nsamples, nparams);
  diagn_acc = zeros(opt_nsamples, 1);
else
  diagnostics = 0;
end

% Main loop.
k = - opt.nomit + 1;
Eold = f(x, varargin{:});	% Evaluate starting energy.
nreject = 0;				% Initialise count of rejected states.
while k <= opt_nsamples

  xold = x;
  % Sample a new point from the proposal distribution
  x = xold + randn(1, nparams)*stddev;

  % Now apply Metropolis algorithm.
  Enew = f(x, varargin{:});	% Evaluate new energy.
  a = exp(Eold - Enew);			% Acceptance threshold.
  if (diagnostics & k > 0)
    diagn_pos(k,:) = x;
    diagn_acc(k,:) = a;
  end
  if (opt_display > 1)
    fprintf(1, 'New position is\n');
    disp(x);
  end

  if a > rand(1)	% Accept the new state.
    Eold = Enew;
    if (opt_display > 0)
      fprintf(1, 'Finished step %4d  Threshold: %g\n', k, a);
    end
  else			% Reject the new state
    if k > 0
      nreject = nreject + 1;
    end
    x = xold;	% Reset position
    if (opt_display > 0)
      fprintf(1, '  Sample rejected %4d.  Threshold: %g\n', k, a);
    end
  end
  if k > 0
    samples(k,:) = x;			% Store sample.
    if en_save
      energies(k) = Eold;		% Store energy.
    end
  end
  k = k + 1;
end

if (opt_display > 0)
  fprintf(1, '\nFraction of samples rejected:  %g\n', ...
          nreject/(opt_nsamples));
end

if diagnostics
  diagn.pos = diagn_pos;
  diagn.acc = diagn_acc;
end

% Return complete state of the sampler.
function state = get_state(f)

state.randstate = rand('state');
state.randnstate = randn('state');
return

% Set state of sampler, either from full state, or with an integer
function set_state(f, x)

if isnumeric(x)
  rand('state', x);
  randn('state', x);
else
  if ~isstruct(x)
    error('Second argument to metrop must be number or state structure');
  end
  if (~isfield(x, 'randstate') | ~isfield(x, 'randnstate'))
    error('Second argument to metrop must contain correct fields')
  end
  rand('state', x.randstate);
  randn('state', x.randnstate);
end
return

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [cfg] = ft_singleplotER(cfg, varargin)

% FT_SINGLEPLOTER plots the event-related fields or potentials of a single
% channel or the average over multiple channels. Multiple datasets can be
% overlayed.
%
% Use as
%   ft_singleplotER(cfg, data)
% or
%   ft_singleplotER(cfg, data1, data2, ..., datan)
%
% The data can be an erp/erf produced by FT_TIMELOCKANALYSIS, a power
% spectrum produced by FT_FREQANALYSIS or connectivity spectrum produced by
% FT_CONNECTIVITYANALYSIS.
%
% The configuration can have the following parameters:
%   cfg.parameter     = field to be plotted on y-axis, for example 'avg', 'powspctrm' or 'cohspctrm' (default is automatic)
%   cfg.maskparameter = field in the first dataset to be used for masking of data; this is not supported when
%                       computing the mean over multiple channels, or when giving multiple input datasets (default = [])
%   cfg.maskstyle     = style used for masking of data, 'box', 'thickness' or 'saturation' (default = 'box')
%   cfg.maskfacealpha = mask transparency value between 0 and 1
%   cfg.xlim          = 'maxmin' or [xmin xmax] (default = 'maxmin')
%   cfg.ylim          = 'maxmin', 'maxabs', 'zeromax', 'minzero', or [ymin ymax] (default = 'maxmin')
%   cfg.channel       = Nx1 cell-array with selection of channels (default = 'all'), see FT_CHANNELSELECTION for details
%   cfg.title         = string, title of plot
%   cfg.showlegend    = 'yes' or 'no', show the legend with the colors (default = 'no')
%   cfg.refchannel    = name of reference channel for visualising connectivity, can be 'gui'
%   cfg.baseline      = 'yes', 'no' or [time1 time2] (default = 'no'), see ft_timelockbaseline
%   cfg.baselinetype  = 'absolute' or 'relative' (default = 'absolute')
%   cfg.trials        = 'all' or a selection given as a 1xn vector (default = 'all')
%   cfg.fontsize      = font size of title (default = 8)
%   cfg.hotkeys       = enables hotkeys (leftarrow/rightarrow/uparrow/downarrow/m) for dynamic zoom and translation (ctrl+) of the axes
%   cfg.interactive   = interactive plot 'yes' or 'no' (default = 'yes')
%                       in a interactive plot you can select areas and produce a new
%                       interactive plot when a selected area is clicked. multiple areas
%                       can be selected by holding down the shift key.
%   cfg.figure        = 'yes' or 'no', whether to open a new figure. You can also specify a figure handle from FIGURE, GCF or SUBPLOT. (default = 'yes')
%   cfg.position      = location and size of the figure, specified as [left bottom width height] (default is automatic)
%   cfg.renderer      = string, 'opengl', 'zbuffer', 'painters', see RENDERERINFO (default is automatic, try 'painters' when it crashes)
%   cfg.linestyle     = linestyle/marker type, see options of the PLOT function (default = '-')
%                       can be a single style for all datasets, or a cell-array containing one style for each dataset
%   cfg.linewidth     = linewidth in points (default = 0.5)
%   cfg.linecolor     = color(s) used for plotting the dataset(s). The default is defined in LINEATTRIBUTES_COMMON, see
%                       the help of this function for more information
%   cfg.directionality = '', 'inflow' or 'outflow' specifies for
%                       connectivity measures whether the inflow into a
%                       node, or the outflow from a node is plotted. The
%                       (default) behavior of this option depends on the dimor
%                       of the input data (see below).
%   cfg.select        = 'intersect' or 'union' (default = 'intersect')
%                       with multiple input arguments determines the
%                       pre-selection of the data that is considered for
%                       plotting.
%   cfg.showlocations = 'no' (default), or 'yes'. plot a small spatial layout of all sensors, highlighting the specified subset
%   cfg.layouttopo    = filename, or struct (see FT_PREPARE_LAYOUT) used for showing the locations with cfg.showlocations = 'yes'
%
% The following options for the scaling of the EEG, EOG, ECG, EMG, MEG and NIRS channels
% is optional and can be used to bring the absolute numbers of the different
% channel types in the same range (e.g. fT and uV). The channel types are determined
% from the input data using FT_CHANNELSELECTION.
%   cfg.eegscale       = number, scaling to apply to the EEG channels prior to display
%   cfg.eogscale       = number, scaling to apply to the EOG channels prior to display
%   cfg.ecgscale       = number, scaling to apply to the ECG channels prior to display
%   cfg.emgscale       = number, scaling to apply to the EMG channels prior to display
%   cfg.megscale       = number, scaling to apply to the MEG channels prior to display
%   cfg.gradscale      = number, scaling to apply to the MEG gradiometer channels prior to display (in addition to the cfg.megscale factor)
%   cfg.magscale       = number, scaling to apply to the MEG magnetometer channels prior to display (in addition to the cfg.megscale factor)
%   cfg.nirsscale      = number, scaling to apply to the NIRS channels prior to display
%   cfg.mychanscale    = number, scaling to apply to the channels specified in cfg.mychan
%   cfg.mychan         = Nx1 cell-array with selection of channels
%   cfg.chanscale      = Nx1 vector with scaling factors, one per channel specified in cfg.channel
%
% For the plotting of directional connectivity data the cfg.directionality
% option determines what is plotted. The default value and the supported
% functionality depend on the dimord of the input data. If the input data
% is of dimord 'chan_chan_XXX', the value of directionality determines
% whether, given the reference channel(s), the columns (inflow), or rows
% (outflow) are selected for plotting. In this situation the default is
% 'inflow'. Note that for undirected measures, inflow and outflow should
% give the same output. If the input data is of dimord 'chancmb_XXX', the
% value of directionality determines whether the rows in data.labelcmb are
% selected. With 'inflow' the rows are selected if the refchannel(s) occur in
% the right column, with 'outflow' the rows are selected if the
% refchannel(s) occur in the left column of the labelcmb-field. Default in
% this case is '', which means that all rows are selected in which the
% refchannel(s) occur. This is to robustly support linearly indexed
% undirected connectivity metrics. In the situation where undirected
% connectivity measures are linearly indexed, specifying 'inflow' or
% 'outflow' can result in unexpected behavior.
%
% to facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% if you specify this option the input data will be read from a *.mat
% file on disk. this mat files should contain only a single variable named 'data',
% corresponding to the input structure.
%
% See also FT_SINGLEPLOTTFR, FT_MULTIPLOTER, FT_MULTIPLOTTFR, FT_TOPOPLOTER, FT_TOPOPLOTTFR

% Copyright (C) 2003-2006, Ole Jensen
% Copyright (C) 2006-2022, Donders Centre for Cognitive Neuroimaging
%
% 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$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEVELOPERS NOTE: This code is organized in a similar fashion for multiplot/singleplot/topoplot
% and for ER/TFR and should remain consistent over those 6 functions.
% Section 1: general cfg handling that is independent from the data
% Section 2: data handling, this also includes converting bivariate (chan_chan and chancmb) into univariate data
% Section 3: select the data to be plotted and determine min/max range
% Section 4: do the actual plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Section 1: general cfg handling that is independent from the data

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar varargin
ft_preamble provenance varargin

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

Ndata = numel(varargin);
for i=1:Ndata
  % check if the input data is valid for this function
  varargin{i} = ft_checkdata(varargin{i}, 'datatype', {'timelock', 'freq'});
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels', 'trial'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'unused',     {'cohtargetchannel'});
cfg = ft_checkconfig(cfg, 'renamedval', {'zlim', 'absmax', 'maxabs'});
cfg = ft_checkconfig(cfg, 'renamedval', {'directionality', 'feedforward', 'outflow'});
cfg = ft_checkconfig(cfg, 'renamedval', {'directionality', 'feedback',    'inflow'});
cfg = ft_checkconfig(cfg, 'renamed',    {'matrixside',     'directionality'});
cfg = ft_checkconfig(cfg, 'renamed',    {'channelindex',   'channel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'channelname',    'channel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'cohrefchannel',  'refchannel'});
cfg = ft_checkconfig(cfg, 'renamed',	  {'zparam',         'parameter'});
cfg = ft_checkconfig(cfg, 'renamed',    {'graphcolor',     'linecolor'});
cfg = ft_checkconfig(cfg, 'deprecated', {'xparam'});
cfg = ft_checkconfig(cfg, 'renamed',    {'newfigure',      'figure'});

% set the defaults
cfg.baseline        = ft_getopt(cfg, 'baseline',      'no');
cfg.trials          = ft_getopt(cfg, 'trials',        'all', 1);
cfg.xlim            = ft_getopt(cfg, 'xlim',          'maxmin');
cfg.ylim            = ft_getopt(cfg, 'ylim',          'maxmin');
cfg.zlim            = ft_getopt(cfg, 'zlim',          'maxmin');
cfg.comment         = ft_getopt(cfg, 'comment',        strcat([date '\n']));
cfg.axes            = ft_getopt(cfg, ' axes',         'yes');
cfg.fontsize        = ft_getopt(cfg, 'fontsize',       8);
cfg.interpreter     = ft_getopt(cfg, 'interpreter',   'none');  % none, tex or latex
cfg.hotkeys         = ft_getopt(cfg, 'hotkeys',       'yes');
cfg.interactive     = ft_getopt(cfg, 'interactive',   'yes');
cfg.maskparameter   = ft_getopt(cfg, 'maskparameter',  []);
cfg.colorgroups     = ft_getopt(cfg, 'colorgroups',   'condition'); % this is the only supported option
cfg.linecolor       = ft_getopt(cfg, 'linecolor',     []);
cfg.linestyle       = ft_getopt(cfg, 'linestyle',     '-');
cfg.linewidth       = ft_getopt(cfg, 'linewidth',      0.5);
cfg.maskstyle       = ft_getopt(cfg, 'maskstyle',     'box');
cfg.maskfacealpha   = ft_getopt(cfg, 'maskfacealpha', 1);
cfg.channel         = ft_getopt(cfg, 'channel',       'all');
cfg.title           = ft_getopt(cfg, 'title',          []);
cfg.directionality  = ft_getopt(cfg, 'directionality', []);
cfg.figurename      = ft_getopt(cfg, 'figurename',     []);
cfg.preproc         = ft_getopt(cfg, 'preproc',        []);
cfg.frequency       = ft_getopt(cfg, 'frequency',     'all'); % needed for frequency selection with TFR data
cfg.latency         = ft_getopt(cfg, 'latency',       'all'); % needed for latency selection with TFR data, FIXME, probably not used
cfg.showlegend      = ft_getopt(cfg, 'showlegend',    'no');
cfg.renderer        = ft_getopt(cfg, 'renderer',       []); % let MATLAB decide on the default
cfg.select          = ft_getopt(cfg, 'select',         'intersect'); % for ft_selectdata
cfg.showlocations   = ft_getopt(cfg, 'showlocations',  'no');
cfg.layouttopo      = ft_getopt(cfg, 'layouttopo');

% this is needed for the figure title and correct labeling of linecolor later on
if isfield(cfg, 'dataname') && ~isempty(cfg.dataname)
  dataname = cfg.dataname;
elseif isfield(cfg, 'inputfile') && ~isempty(cfg.inputfile)
  dataname = cfg.inputfile;
elseif nargin>1
  dataname = arrayfun(@inputname, 2:nargin, 'UniformOutput', false);
else
  dataname = {};
end

%% Section 2: data handling, this also includes converting bivariate (chan_chan and chancmb) into univariate data

for i=1:Ndata
  dtype{i}    = ft_datatype(varargin{i});
  hastime(i)  = isfield(varargin{i}, 'time');
  hasfreq(i)  = isfield(varargin{i}, 'freq');
end

% check if the input has consistent datatypes
if ~all(strcmp(dtype, dtype{1})) || ~all(hastime==hastime(1)) || ~all(hasfreq==hasfreq(1))
  ft_error('different datatypes are not allowed as input');
else
  dtype   = dtype{1};
  hastime = hastime(1);
  hasfreq = hasfreq(1);
end

% Set x/y/parameter according to datatype and dimord
switch dtype
  case 'timelock'
    xparam = 'time';
    if isfield(varargin{1}, 'trial')
      cfg.parameter = ft_getopt(cfg, 'parameter', 'trial');
    elseif isfield(varargin{1}, 'individual')
      cfg.parameter = ft_getopt(cfg, 'parameter', 'individual');
    elseif isfield(varargin{1}, 'avg')
      cfg.parameter = ft_getopt(cfg, 'parameter', 'avg');
    end
  case 'freq'
    if hastime && hasfreq
      xparam = 'time';   % average over selected frequencies
    else
      xparam = 'freq';
    end
    cfg.parameter = ft_getopt(cfg, 'parameter', 'powspctrm');
  case 'comp'
    % not supported
  otherwise
    % not supported
end

% check whether rpt/subj is present and remove if necessary
dimord = getdimord(varargin{1}, cfg.parameter);
dimtok = tokenize(dimord, '_');
hasrpt = any(ismember(dimtok, {'rpt' 'subj'}));

if ~hasrpt
  assert(isequal(cfg.trials, 'all') || isequal(cfg.trials, 1), 'incorrect specification of cfg.trials for data without repetitions');
else
  assert(~isempty(cfg.trials), 'empty specification of cfg.trials for data with repetitions');
end

% parse cfg.channel
if isfield(cfg, 'channel') && isfield(varargin{1}, 'label')
  cfg.channel = ft_channelselection(cfg.channel, varargin{1}.label);
elseif isfield(cfg, 'channel') && isfield(varargin{1}, 'labelcmb')
  cfg.channel = ft_channelselection(cfg.channel, unique(varargin{1}.labelcmb(:)));
end

% apply baseline correction
if ~strcmp(cfg.baseline, 'no')
  tmpcfg = keepfields(cfg, {'baseline', 'baselinetype', 'baselinewindow', 'demean', 'parameter', 'channel'});
  for i=1:Ndata
    % keep mask-parameter if it is set
    if ~isempty(cfg.maskparameter)
      tempmask = varargin{i}.(cfg.maskparameter);
    end
    if strcmp(dtype, 'timelock') && strcmp(xparam, 'time')
      varargin{i} = ft_timelockbaseline(tmpcfg, varargin{i});
    elseif strcmp(dtype, 'freq') && strcmp(xparam, 'time')
      varargin{i} = ft_freqbaseline(tmpcfg, varargin{i});
    elseif strcmp(dtype, 'freq') && strcmp(xparam, 'freq')
      ft_error('baseline correction is not supported for spectra without a time dimension');
    else
      ft_warning('baseline correction not applied, please set xparam');
    end
    % put mask-parameter back if it is set
    if ~isempty(cfg.maskparameter)
      varargin{i}.(cfg.maskparameter) = tempmask;
    end
  end
end

% channels should NOT be selected and averaged here, since a topoplot might follow in interactive mode
tmpcfg = keepfields(cfg, {'trials', 'select', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
if hasrpt
  tmpcfg.avgoverrpt = 'yes';
else
  tmpcfg.avgoverrpt = 'no';
end
if hastime && hasfreq
  tmpcfg.avgoverfreq = 'yes';       % average over selected frequencies
  tmpcfg.frequency = cfg.frequency; % not to be confused with cfg.xlim or cfg.ylim
  tmpcfg.keepfreqdim = 'no';
else
  tmpcfg.avgoverfreq = 'no';
end
tmpvar = varargin{1};
[varargin{:}] = ft_selectdata(tmpcfg, varargin{:});
% restore the provenance information and put back cfg.channel
tmpchannel  = cfg.channel;
[cfg, varargin{:}] = rollback_provenance(cfg, varargin{:});
cfg.channel = tmpchannel;

if isfield(tmpvar, cfg.maskparameter) && ~isfield(varargin{1}, cfg.maskparameter)
  % the mask parameter is not present after ft_selectdata, because it is
  % not included in all input arguments. Make the same selection and copy
  % it over
  tmpvar = ft_selectdata(tmpcfg, tmpvar);
  varargin{1}.(cfg.maskparameter) = tmpvar.(cfg.maskparameter);
end

clear tmpvar tmpcfg dimord dimtok hastime hasfreq hasrpt

% ensure that the preproc specific options are located in the cfg.preproc
% substructure, but also ensure that the field 'refchannel' remains at the
% highest level in the structure. This is a little hack by JM because the field
% refchannel can relate to connectivity or to an EEg reference.

if isfield(cfg, 'refchannel'), refchannelincfg = cfg.refchannel; cfg = rmfield(cfg, 'refchannel'); end
cfg = ft_checkconfig(cfg, 'createsubcfg',  {'preproc'});
if exist('refchannelincfg', 'var'), cfg.refchannel  = refchannelincfg; end

if ~isempty(cfg.preproc)
  % preprocess the data, i.e. apply filtering, baselinecorrection, etc.
  fprintf('applying preprocessing options\n');
  if ~isfield(cfg.preproc, 'feedback')
    cfg.preproc.feedback = cfg.interactive;
  end
  for i=1:Ndata
    varargin{i} = ft_preprocessing(cfg.preproc, varargin{i});
  end
end

% Handle the bivariate case
dimord = getdimord(varargin{1}, cfg.parameter);
if startsWith(dimord, 'chan_chan_') || startsWith(dimord, 'chancmb_')
  % convert the bivariate data to univariate and call this plotting function again
  cfg.originalfunction = 'ft_singleplotER';
  cfg.trials = 'all'; % trial selection has been taken care off
  bivariate_common(cfg, varargin{:});
  return
end

% Apply channel-type specific scaling
fn = fieldnames(cfg);
fn = setdiff(fn, {'skipscale', 'showscale', 'gridscale'}); % these are for the layout and plotting, not for CHANSCALE_COMMON
fn = fn(endsWith(fn, 'scale') | startsWith(fn, 'mychan') | strcmp(fn, 'channel') | strcmp(fn, 'parameter'));
tmpcfg = keepfields(cfg, fn);
if ~isempty(tmpcfg)
  for i=1:Ndata
    varargin{i} = chanscale_common(tmpcfg, varargin{i});
  end
  % remove the scaling fields from the configuration, to prevent them from being called again in interactive mode
  % but keep the parameter and channel field
  cfg = removefields(cfg, setdiff(fn, {'parameter', 'channel'}));
else
  % do nothing
end

%% Section 3: select the data to be plotted and determine min/max range

if istrue(cfg.showlocations)
  % Read or create the layout that will be used for plotting, if specified
  tmpcfg = keepfields(cfg, {'rows', 'columns', 'commentpos', 'scalepos', 'projection', 'viewpoint', 'rotate', 'width', 'height', 'elec', 'grad', 'opto', 'layouttopo', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  tmpcfg.skipcomnt = 'yes';
  tmpcfg.skipscale = 'yes';
  tmpcfg.pointcolor = cfg.linecolor; % switch of name for ft_prepare_layout
  if ~isempty(cfg.layouttopo)
    tmpcfg.layout = cfg.layouttopo;
  elseif isfield(cfg, 'layout') && ~isempty(cfg.layout)
    tmpcfg.layout = cfg.layout;
  else
    ft_warning('no explicit layout specified, attempting to create one from the data');
  end
  cfg.layouttopo = ft_prepare_layout(tmpcfg, varargin{1});
end

% Take the desided subselection of channels, this is the same in all datasets
[selchan] = match_str(varargin{1}.label, cfg.channel);

% Get physical min/max range of x, i.e. time or frequency
if strcmp(cfg.xlim, 'maxmin')
  % Find maxmin throughout all varargins:
  xmin = [];
  xmax = [];
  for i=1:Ndata
    xmin = min([xmin varargin{i}.(xparam)]);
    xmax = max([xmax varargin{i}.(xparam)]);
  end
else
  xmin = cfg.xlim(1);
  xmax = cfg.xlim(2);
end

% Get the index of the nearest bin, this is the same in all datasets
xminindx = nearest(varargin{1}.(xparam), xmin);
xmaxindx = nearest(varargin{1}.(xparam), xmax);
xmin = varargin{1}.(xparam)(xminindx);
xmax = varargin{1}.(xparam)(xmaxindx);
selx = xminindx:xmaxindx;
xval = varargin{1}.(xparam)(selx);

% get physical y-axis range, i.e. parameter to be plotted
if ~isnumeric(cfg.ylim)
  % find maxmin throughout all varargins
  ymin = +inf;
  ymax = -inf;
  for i=1:Ndata
    % select the channels in the data that match with the layout and that are selected for plotting
    dat = nanmean(varargin{i}.(cfg.parameter)(selchan, selx), 1); % mean over channels

    ymin = min(ymin, min(dat(:)));
    ymax = max(ymax, max(dat(:)));
  end
  switch cfg.ylim
    case 'maxmin'
      % keep them as they are
    case 'maxabs'
      ymax = max(abs(ymax), abs(ymin));
      ymin = -ymax;
    case 'zeromax'
      ymin = 0;
    case 'minzero'
      ymax = 0;
    otherwise
      ft_error('invalid specification of cfg.ylim');
  end
else
  ymin = cfg.ylim(1);
  ymax = cfg.ylim(2);
end

% gather the data from all input data structures
datamatrix = zeros(Ndata, length(selx));
for i=1:Ndata
  datamatrix(i,:) = mean(varargin{i}.(cfg.parameter)(selchan, selx), 1); % mean over channels
end

% gather the mask from the first input data structure
if ~isempty(cfg.maskparameter)
  maskmatrix = mean(varargin{1}.(cfg.maskparameter)(selchan, selx), 1); % mean over channels
else
  % create an Nx0 matrix
  maskmatrix = zeros(length(selchan), 0);
end

%% Section 4: do the actual plotting

% determine the coloring of channels/conditions
[linecolor, linestyle, linewidth] = lineattributes_common(cfg, varargin{:});
linecolor = mean(linecolor(selchan, :, :), 1);
linestyle = linestyle(selchan(1), :);
linewidth = mean(linewidth(selchan, :), 1);

% open a new figure, or add it to the existing one
open_figure(keepfields(cfg, {'figure', 'position', 'visible', 'renderer', 'figurename', 'title'}));

yval = datamatrix;
mask = maskmatrix;

if strcmp(cfg.maskstyle, 'difference')
  % combine the conditions in a single plot, highlight the difference
  ft_plot_vector(xval, yval, 'color', permute(linecolor, [3 2 1]), 'style', linestyle(1,:), 'linewidth', linewidth(1,:), 'highlight', mask, 'highlightstyle', cfg.maskstyle, 'hlim', [xmin xmax], 'vlim', [ymin ymax], 'facealpha', cfg.maskfacealpha);
else
  % loop over the conditions, plot them on top of each other
  for i=1:Ndata
    ft_plot_vector(xval, yval(i,:), 'color', linecolor(1,:,i), 'style', linestyle{1,i}, 'linewidth', linewidth(1,i), 'highlight', mask, 'highlightstyle', cfg.maskstyle, 'hlim', [xmin xmax], 'vlim', [ymin ymax], 'facealpha', cfg.maskfacealpha);
  end
end

if ischar(linecolor)
  set(gca, 'ColorOrder', char2rgb(linecolor))
elseif isnumeric(linecolor)
  set(gca, 'ColorOrder', shiftdim(linecolor(1,:,:),1)');
end

% show the legend with the colors of the conditions
if istrue(cfg.showlegend) && Ndata>1
  if strcmp(cfg.maskstyle, 'difference')
    colorLabels = {'difference'};
  else
    colorLabels = {};
  end
  for i=1:Ndata
    if ischar(linecolor)
      colorLabels{end+1} = [dataname{i} '='         linecolor(i)     ];
    elseif isnumeric(linecolor)
      colorLabels{end+1} = [dataname{i} '=' num2str(linecolor(1, :, i)) ];
    end
  end
  legend(colorLabels)
end

% set xlim and ylim
if xmin~=xmax
  xlim([xmin xmax]);
end
if ymin~=ymax
  ylim([ymin ymax]);
end

% adjust mask box extents to ymin/ymax
if ~isempty(cfg.maskparameter)
  ptchs = findobj(gcf, 'type', 'patch');
  for i = 1:length(ptchs)
    YData = get(ptchs(i), 'YData');
    YData(YData == min(YData)) = ymin;
    YData(YData == max(YData)) = ymax;
    set(ptchs(i), 'YData',YData);
  end
end

% Set callback to adjust axes
if strcmp('yes', cfg.hotkeys)
  %  attach data and cfg to figure and attach a key listener to the figure
  set(gcf, 'KeyPressFcn', {@key_sub, xmin, xmax, ymin, ymax})
end

% Create axis title containing channel name(s) and channel number(s):
if ~isempty(cfg.title)
  t = cfg.title;
else
  if length(cfg.channel) == 1
    t = [char(cfg.channel) ' / ' num2str(selchan) ];
  else
    t = sprintf('mean(%0s)', join_str(', ', cfg.channel));
  end
end
title(t, 'fontsize', cfg.fontsize, 'interpreter', cfg.interpreter);

% set the figure window title, add channel labels if number is small
if isempty(get(gcf, 'Name'))
  if length(selchan) < 5
    chans = join_str(', ', cfg.channel);
  else
    chans = '<multiple channels>';
  end
  if ~isempty(cfg.figurename)
    set(gcf, 'name', cfg.figurename);
    set(gcf, 'NumberTitle', 'off');
  elseif ~isempty(dataname)
    set(gcf, 'Name', sprintf('%d: %s: %s (%s)', double(gcf), mfilename, join_str(', ', dataname), chans));
    set(gcf, 'NumberTitle', 'off');
  else
    set(gcf, 'Name', sprintf('%d: %s (%s)', double(gcf), mfilename, chans));
    set(gcf, 'NumberTitle', 'off');
  end
end

if istrue(cfg.showlocations)
  hpos = xmin+(xmax-xmin)*0.1;
  vpos = ymin+(ymax-ymin)*0.9;
  h    = 0.2*(ymax-ymin);
  w    = 0.2*(xmax-xmin);
  pointcolor  = zeros(numel(cfg.layouttopo.label),3);
  pointsize   = ones(numel(cfg.layouttopo.label),2);
  pointsize(selchan) = 4;
  pointsymbol = 'o';
  ft_plot_layout(cfg.layouttopo, 'box', 'no', 'label', 'off', 'pointsize', pointsize, 'pointcolor', pointcolor, 'pointsymbol', pointsymbol, 'hpos', hpos, 'vpos', vpos, 'width', w, 'height', h);
end

% make the figure interactive
if strcmp(cfg.interactive, 'yes')
  % add the cfg/data/channel information to the figure under identifier linked to this axis
  ident                  = ['axh' num2str(round(sum(clock.*1e6)))]; % unique identifier for this axis
  set(gca, 'tag', ident);
  info                   = guidata(gcf);
  info.(ident).cfg       = cfg;
  info.(ident).varargin  = varargin;
  info.(ident).dataname  = dataname;
  if exist('linecolor', 'var')
    info.(ident).linecolor   = linecolor;
  end
  guidata(gcf, info);
  set(gcf, 'windowbuttonupfcn',     {@ft_select_range, 'multiple', false, 'yrange', false, 'callback', {@select_topoplotER}, 'event', 'windowbuttonupfcn'});
  set(gcf, 'windowbuttondownfcn',   {@ft_select_range, 'multiple', false, 'yrange', false, 'callback', {@select_topoplotER}, 'event', 'windowbuttondownfcn'});
  set(gcf, 'windowbuttonmotionfcn', {@ft_select_range, 'multiple', false, 'yrange', false, 'callback', {@select_topoplotER}, 'event', 'windowbuttonmotionfcn'});
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous varargin
ft_postamble provenance
ft_postamble savefig

% add a menu to the figure, but only if the current figure does not have subplots
menu_fieldtrip(gcf, cfg, false);

if ~ft_nargout
  % don't return anything
  clear cfg
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION which is called after selecting a time range
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function select_topoplotER(range, varargin)

% fetch cfg/data based on axis indentifier given as tag
ident    = get(gca, 'tag');
info     = guidata(gcf);
cfg      = info.(ident).cfg;
varargin = info.(ident).varargin;
if ~isempty(range)
  cfg = removefields(cfg, 'inputfile');   % the reading has already been done and varargin contains the data
  cfg = removefields(cfg, 'showlabels');  % this is not allowed in topoplotER
  cfg.baseline = 'no';                    % make sure the next function does not apply a baseline correction again
  cfg.dataname = info.(ident).dataname;   % put data name in here, this cannot be resolved by other means
  cfg.channel = 'all';                    % make sure the topo displays all channels, not just the ones in this singleplot
  cfg.comment = 'auto';
  cfg.trials = 'all';                     % trial selection has already been taken care of
  cfg.xlim = range(1:2);
  % if user specified a ylim, copy it over to the zlim of topoplot
  if isfield(cfg, 'ylim')
    cfg.zlim = cfg.ylim;
    cfg = rmfield(cfg, 'ylim');
  end
  fprintf('selected cfg.xlim = [%f %f]\n', cfg.xlim(1), cfg.xlim(2));
  % ensure that the new figure appears at the same position
  cfg.figure = 'yes';
  cfg.position = get(gcf, 'Position');
  ft_topoplotER(cfg, varargin{:});
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION which handles hot keys in the current plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function key_sub(handle, eventdata, varargin)
xlimits = xlim;
ylimits = ylim;
incr_x = abs(xlimits(2) - xlimits(1)) /10;
incr_y = abs(ylimits(2) - ylimits(1)) /10;

if length(eventdata.Modifier) == 1 && strcmp(eventdata.Modifier{:}, 'control')
  % TRANSLATE by 10%
  switch eventdata.Key
    case 'leftarrow'
      xlim([xlimits(1)+incr_x xlimits(2)+incr_x])
    case 'rightarrow'
      xlim([xlimits(1)-incr_x xlimits(2)-incr_x])
    case 'uparrow'
      ylim([ylimits(1)-incr_y ylimits(2)-incr_y])
    case 'downarrow'
      ylim([ylimits(1)+incr_y ylimits(2)+incr_y])
  end % switch
else
  % ZOOM by 10%
  switch eventdata.Key
    case 'leftarrow'
      xlim([xlimits(1)-incr_x xlimits(2)+incr_x])
    case 'rightarrow'
      xlim([xlimits(1)+incr_x xlimits(2)-incr_x])
    case 'uparrow'
      ylim([ylimits(1)-incr_y ylimits(2)+incr_y])
    case 'downarrow'
      ylim([ylimits(1)+incr_y ylimits(2)-incr_y])
    case 'm'
      xlim([varargin{1} varargin{2}])
      ylim([varargin{3} varargin{4}])
  end % switch
end % if


text:
classdef AffineGridGenerator < dagnn.Layer
%DAGNN.AFFINEGRIDGENERATIOR  Generate an affine grid for bilinear resampling
%   This layer maps 1 x 1 x 6 x N affine transforms to 2 x Ho x Wo x N
%   sampling grids compatible with dagnn.BlilinearSampler.

% (c) 2016 Ankush Gupta

 properties
     Ho = 0;
     Wo = 0;
 end

  properties (Transient)
    % the grid (normalized \in [-1,1]) --> this is cached
    % has the size: [HoWo x 2]
    xxyy ;
  end

  methods
    function outputs = forward(obj, inputs, ~)
      % input is a 1x1x6xN TENSOR corresponding to:
      % [ c1 c2 c5 ]
      % [ c3 c4 c6 ]
      % [  0  0  1 ]
      % i.e., [d1] = [c1 c2]  * [d1] + [c5]
      %       [d2]   [c3 c4]    [d2]   [c6]
      % where, di is the i-th dimension.
      %
      % OUTPUT is a 2xHoxWoxN grid which corresponds to applying
      % the above affine transform to the [-1,1] normalized x,y
      % coordinates.

      %fprintf('affineGridGenerator forward\n');
      useGPU = isa(inputs{1}, 'gpuArray');

      % reshape the tfm params into matrices:
      A = inputs{1};
      nbatch = size(A,4);
      A = reshape(A, 2,3,nbatch);
      L = A(:,1:2,:);
      L = reshape(L,2,2*nbatch); % linear part

      % generate the grid coordinates:
      if isempty(obj.xxyy)
        obj.initGrid(useGPU);
      end

      % transform the grid:
      t = A(:,3,:); % translation
      t = reshape(t,1,2*nbatch);
      g = bsxfun(@plus, obj.xxyy * L, t); % apply the transform
      g = reshape(g, obj.Wo,obj.Ho,2,nbatch);

      % cudnn compatibility:
      g = permute(g, [3,2,1,4]);

      outputs = {g};
    end

    function [derInputs, derParams] = backward(obj, inputs, ~, derOutputs)

      useGPU = isa(derOutputs{1}, 'gpuArray');
      dY = derOutputs{1};
      nbatch = size(dY,4);

      % cudnn compatibility:
      dY = permute(dY, [3,2,1,4]);

      % create the gradient buffer:
      dA = zeros([2,3,nbatch], 'single');
      if useGPU, dA = gpuArray(dA); end

      dY = reshape(dY, obj.Ho*obj.Wo, 2*nbatch);
      % gradient wrt the linear part:
      dL = obj.xxyy' * dY;
      dL = reshape(dL,2,2,nbatch);
      dA(:,1:2,:) = dL;

      % gradient wrt translation (or bias):
      dt = reshape(sum(dY,1),2,1,nbatch);
      dA(:,3,:) = dt;

      dA = reshape(dA, size(inputs{1}));
      derInputs = {dA};
      derParams = {};
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      nBatch = inputSizes{1}(4);
      outputSizes = {[2, obj.Ho, obj.Wo, nBatch]};
    end

    function obj = AffineGridGenerator(varargin)
      obj.load(varargin) ;
      % get the output sizes:
      obj.Ho = obj.Ho ;
      obj.Wo = obj.Wo ;
      obj.xxyy = [] ;
    end

    function obj = reset(obj)
      reset@dagnn.Layer(obj) ;
      obj.xxyy = [] ;
    end

    function initGrid(obj, useGPU)
      % initialize the grid:
      % this is a constant
      xi = linspace(-1, 1, obj.Ho);
      yi = linspace(-1, 1, obj.Wo);

      [yy,xx] = meshgrid(xi,yi);
      xxyy = [yy(:), xx(:)] ; % Mx2
      if useGPU
        xxyy = gpuArray(xxyy);
      end
      obj.xxyy = xxyy ;
    end

  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function L = cmpndLikelihood(noise, mu, varsigma, y)

% CMPNDLIKELIHOOD Likelihood of data under compound noise model.

% NOISE

% NOISE

L = zeros(size(mu));
for i = 1:length(noise.comp)
  L(:, i) = noiseLikelihood(noise.comp{i},...
                mu(:, i), ...
                varsigma(:, i), ...
                y(:, i));
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% Build and initialize the computational graph for extracting log Mel
% filterbank features
%
function [layer, para] = BuildFbankExtractionNet()
para.output = 'tmp';
para.IO.nStream = 1;
para.NET.sequential = 1;
para.cost_func.layer_idx = [];

para = ConfigBasicSTFT(para);
layer = genNetworkFbankExtraction(para.topology);     % generate the network graph

% generating the scaling factor for the input, as we will need to use a
% small constant in the logarithm. We need to make sure that the power of
% speech are larger than this constant most of the time.
scale = 1e4;        % we hard code the scale to be a constant so that all network will use the same number
scale = scale/2^16; % note that we are using int16 to store waveform samples, so need to scale down
layer = InitWavScaleLayer(layer, scale);

% set Mel filterbank linear transform
layer = InitMelLayer(layer, para);

para.out_layer_idx = length(layer);

para = ParseOptions2(para);

end

author: singaxiong

repo: SignalGraph


text:
function title = p19_title ( title )

%*****************************************************************************80
%
%% P19_TITLE returns the title of problem 19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'The jumping cosine';

  return
end


text:
function varargout=voicebox(varargin)
%
% For calling details please see v_voicebox.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_voicebox.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: voicebox.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_voicebox(varargin{:});
else
    v_voicebox(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [t,z] = plotDNLS(prob,opts,xb,confStats,tspan)
%plotDNLS Plot Parameter Estimation Problem

%   Copyright (C) 2013 Jonathan Currie (IPL)

if(nargin < 5), tspan = []; end
if(nargin < 4), confStats = []; end

%Measurement Plot Color
measC = [0.4 0.4 0.4];
%Measurement Plot Style
measS = 'o';
%Initial Condition plot style
icS = 'sq';

%Plot confidence lines if present
if(~isempty(confStats))
    plot(NaN,NaN); %hack to prevent patch removing outer borders
    plotConfReg(confStats);
    hold on;
end

% Insert initial conditions if also solved
ind = isnan(prob.odez0);
if(any(ind))
    len = sum(ind);
    estZ0 = xb(end-len+1:end);
    prob.odez0(ind) = estZ0;
    xb = xb(1:end-len);
end

% Generate smooth plot
dopts = optidynset(opts.dynamicOpts);
if(isempty(tspan))
    tspan = [prob.xdata(1) prob.xdata(end)];
end
ode = @(t,z) prob.ode(t,z,xb);
switch(dopts.integrator)
    case 'ode45'
        [t,z] = ode45(ode,tspan,prob.odez0,dopts.odeOpts);
    case 'ode15s'
        [t,z] = ode15s(ode,tspan,prob.odez0,dopts.odeOpts);
    case 'ode23'
        [t,z] = ode23(ode,tspan,prob.odez0,dopts.odeOpts);
    case 'ode113'
        [t,z] = ode113(ode,tspan,prob.odez0,dopts.odeOpts);
    case 'ode23t'
        [t,z] = ode23t(ode,tspan,prob.odez0,dopts.odeOpts);
    case 'ode23tb'
        [t,z] = ode23tb(ode,tspan,prob.odez0,dopts.odeOpts);
    case 'ode23s'
        [t,z] = ode23s(ode,tspan,prob.odez0,dopts.odeOpts);
    case 'ode15i'
        [t,z] = ode15i(ode,tspan,prob.odez0,dopts.odeOpts);
    otherwise
        optiwarn('OPTI:Plot','Unknown integrator chosen! Using ode45.');
        [t,z] = ode45(ode,tspan,prob.odez0,dopts.odeOpts);
end
if(nargout > 1), return; end %just for evaluating above

%Plot smooth sim
plot(t,z);
title('Dynamic Parameter Estimation Solution');
xlabel('time'); ylabel('z');
hold on;

%Plot Solved Initial Conditions
if(any(ind))
    for i = 1:sum(ind)
        plot(tspan(1),estZ0(i),icS,'Color',measC);
    end
end

%If we have xdata_old, then should be in cell format, plot each cell
if(isfield(prob.misc,'xdata_orig') && ~isempty(prob.misc.xdata_orig))
    if(iscell(prob.misc.xdata_orig) && iscell(prob.misc.ydata_orig))
        for i = 1:length(prob.misc.xdata_orig)
            plot(prob.misc.xdata_orig{i},prob.misc.ydata_orig{i},measS,'Color',measC);
        end
    elseif(~iscell(prob.misc.xdata_orig) && ~iscell(prob.misc.ydata_orig))
        for i = 1:size(prob.misc.ydata_orig,2)
            plot(prob.misc.xdata_orig,prob.misc.ydata_orig(:,i),measS,'Color',measC);
        end
    else
        error('Expected misc.xdata_orig and misc.ydata_orig to be cell arrays');
    end
else
    %Otherwise reshape ydata (always a column in OPTI) and plot
    if(~isempty(dopts.stateIndex))
        nstates = length(prob.odez0(dopts.stateIndex));
    else
        nstates = length(prob.odez0);
    end
    ydata = reshape(prob.ydata,length(prob.xdata),nstates);
    %Plot
    plot(prob.xdata,ydata,measS,'Color',measC);
end

hold off;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef MIndenter < handle
    % Performs code indenting. Should not be used directly but only by
    % MBeautify.

    properties (Constant)
        Delimiters = {' ', '\f', '\n', '\r', '\t', '\v', ...
            ','};
        KeywordsIncrease = {'function', 'classdef', 'properties', ...
            'methods', 'if', 'for', 'parfor', 'switch', 'try', 'while', ...
            'arguments', 'enumeration'};
        KeywordsSandwich = {'else', 'elseif', 'case', 'otherwise', ...
            'catch'};
        KeywordsDecrease = {'end', 'end;'};
    end

    properties (Access = private)
        Configuration;
    end

    methods
        function obj = MIndenter(configuration)
            % Creates a new formatter using the passed configuration.
            obj.Configuration = configuration;
        end

        function indentedSource = performIndenting(obj, source)
            % indentation strategy:
            % allfunctions (default) = AllFunctionIndent
            % nestedfunctions = MixedFunctionIndent (MATLAB default)
            % noindent = ClassicFunctionIndent
            strategy = obj.Configuration.specialRule('Indentation_Strategy').Value;

            % determine indent string
            indentationCharacter = obj.Configuration.specialRule('IndentationCharacter').Value;
            indentationCount = obj.Configuration.specialRule('IndentationCount').ValueAsDouble;
            if strcmpi(indentationCharacter, 'white-space')
                indent = ' ';
                for i = 2:indentationCount
                    indent = [' ', indent];
                end
            elseif strcmpi(indentationCharacter, 'tab')
                indent = '\t';
            else
                warning('MBeautifier:IllegalSetting:IndentationCharacter', 'MBeautifier: The indentation character must be set to "white-space" or "tab". MBeautifier using MATLAB defaults.');
                indent = '    ';
            end

            % TODO
            makeBlankLinesEmpty = obj.Configuration.specialRule('Indentation_TrimBlankLines').ValueAsDouble;

            % currently in continuation mode (line before ended with ...)?
            continuationMode = 0;
            % layer of indentatino (next line)
            layerNext = 0;
            % this stack keeps track of the keywords
            stack = {};

            % start indenting
            newLine = MBeautifier.Constants.NewLine;
            lines = regexp(source, newLine, 'split');
            for linect = 1:numel(lines)
                % layer of indentation (current line)
                layer = layerNext;

                % remove existing indentation and whitespace
                lines{linect} = strtrim(lines{linect});

                % split line in words
                pattern = ['[', obj.joinString(obj.Delimiters, '|'), ']'];
                words = regexp(lines{linect}, pattern, 'split');
                % ignore empty lines and comments
                if (~isempty(lines{linect}) && (lines{linect}(1) ~= '%'))
                    % find keywords and adjust indent
                    for wordct = 1:numel(words)
                        % detect end of line comments
                        if (strcmp(words{wordct}, '%'))
                            break;
                        end

                        % look for keywords that increase indent
                        if (sum(strcmp(words{wordct}, obj.KeywordsIncrease)))
                            layerNext = layerNext + 1;
                            % push keyword onto stack
                            stack = [stack, words{wordct}];

                            % correction for function keywords according to
                            % configuration
                            if (strcmp(stack{end}, 'function'))
                                switch lower(strategy)
                                    case 'nestedfunctions'
                                        if (numel(stack) == 1)
                                            % top level function
                                            layerNext = layerNext - 1;
                                        else
                                            if (strcmp(stack{end-1}, 'function'))
                                                % nested function
                                                layer = layer + 1;
                                                layerNext = layerNext + 1;
                                            else
                                                % class method
                                                % do nothing
                                            end
                                        end
                                    case 'noindent'
                                        layerNext = layerNext - 1;
                                    otherwise
                                end
                            end

                            % correction for switch
                            if (strcmp(stack{end}, 'switch'))
                                layerNext = layerNext + 1;
                            end
                        end

                        % look for sandwich keywords
                        if (sum(strcmp(words{wordct}, obj.KeywordsSandwich)))
                            if (wordct == 1)
                                % at the beginning, decrease only current indent
                                layer = layer - 1;
                            end
                        end

                        % look for end that decreases the indent
                        if (sum(strcmp(words{wordct}, obj.KeywordsDecrease)))
                            if (wordct == 1)
                                % end at the beginning decreases indent of this line
                                layer = layer - 1;
                            end
                            % inline end may alter the indent of the next line
                            layerNext = layerNext - 1;

                            if isempty(stack)
                               continue
                            end

                            % correction for function keywords according to
                            % configuration
                            if (strcmp(stack{end}, 'function'))
                                switch lower(strategy)
                                    case 'nestedfunctions'
                                        if (numel(stack) == 1)
                                            % top level function
                                            layerNext = layerNext + 1;
                                        else
                                            if (strcmp(stack{end-1}, 'function'))
                                                % nested function
                                                layer = layer - 1;
                                            else
                                                % class method
                                                % do nothing
                                            end
                                        end
                                    case 'noindent'
                                        if (wordct == 1)
                                            % end at the beginning decreases indent of this line
                                            layer = layer + 1;
                                        end
                                        layerNext = layerNext + 1;
                                    otherwise
                                        % do nothing
                                end
                            end

                            % correction for switch
                            if (strcmp(stack{end}, 'switch'))
                                if (wordct == 1)
                                    % end at the beginning decreases indent of this line
                                    layer = layer - 1;
                                end
                                layerNext = layerNext - 1;
                            end

                            % pop keyword
                            stack(end) = [];
                        end
                    end

                    % look for continuation lines
                    if (strcmp(words{end}, '...'))
                        if (~continuationMode)
                            continuationMode = 1;
                            layerNext = layerNext + 1;
                        end
                    else
                        if (continuationMode)
                            continuationMode = 0;
                            layerNext = layerNext - 1;
                        end
                    end
                end

                % add correct indentation
                for ict = 1:layer
                    if ~makeBlankLinesEmpty || ~isempty(lines{linect})
                        lines{linect} = [indent, lines{linect}];
                    end
                end
            end

            indentedSource = obj.joinString(lines, MBeautifier.Constants.NewLine);
        end
    end

    methods (Access = private, Static)
        % TODO: Create a public utility function
       function outStr = joinString(cellStr, delim)
            outStr = '';
            for i = 1:numel(cellStr)
                outStr = [outStr, cellStr{i}, delim];
            end

            outStr(end-numel(delim)+1:end) = '';
        end
    end
end

author: davidvarga

repo: MBeautifier


text:
% [edges] = csaAssign(n,graph)
%
% Compute min-cost assignment with non-negative integral edge weights
% using Andrew Goldberg's CSA package (precise costs version).
%
% INPUT
%	n	Number of nodes in the bipartite graph (must be even).
%	graph	3xm matrix describing graph.
%
% OUTPUT
%	edges	3xn matrix of edges in assignment.
%
% You must ensure that an assignment involving all nodes exists, else
% the code may hang.  This is a feature of the CSA package.  If your
% problem does not necessarily provide such an assignment, then you
% should overlay a high-cost perfect match as a safety net.
%
% Both graph and edges matrices have the same structure.  Each column
% gives a graph edge e.  The two nodes are given by e(1) and e(2):
%
%	e(1) < e(2)
%	1 <= e(1) <= n/2
%	n/2 < e(2) <= n
%
% The edge weight is given by e(3).
%
% Since the output edge matrix should contain one reference to each
% node, sum(sum(edges(1:2,:))) == n*(n+1)/2.
%
% David Martin <dmartin@eecs.berkeley.edu>
% January, 2003

author: Cloud-CV

repo: object-proposals


text:
function [cc] = yd32cc(yd3)
% Convert volume from cubic yards to cubic centimeters*.
% Chad Greene 2012
cc = yd3*764554.85798;

% *Not to be confused with ancient Egyptian cubic cubits.

author: Sable

repo: mcbench-benchmarks


text:
function [xm xcov] = UT(Xi, W)
%
%
[n, kmax] = size(Xi);

xm = 0;
for k=1:kmax
  xm = xm + W(k)*Xi(:, k);
end

xcov = zeros(n, n);
for k=1:kmax
  xcov = xcov + W(k)*(Xi(:, k) - xm)*(Xi(:, k) - xm)';
end
xcov = xcov;

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
function hostname = hlp_hostname()
% Find the hostname in a platform-independent fashion.

% Copyright (C) Christian Kothe, SCCN, 2010, christian@sccn.ucsd.edu
%
% 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

import java.net.InetAddress;
hostname = char(InetAddress.getLocalHost());
hostname = hostname(1:find(hostname=='/',1)-1);

repo: REST


text:
function minDist = distancePointPolygon(point, poly)
%DISTANCEPOINTPOLYGON Shortest distance between a point and a polygon
%
%   DIST = distancePointPolygon(POINT, POLYGON)
%   Computes the shortest distance between the point POINT and the polygon
%   given by POLYGON. POINT is a 1-by-2 row vector, and POLYGON is a N-by-2
%   array containing vertex coordinates.
%   The distance is computed as the minimal distance to the boundary edges.
%
%   Example
%     % Computes the distance between a point and a square
%     square = [0 0; 10 0;10 10;0 10];
%     p0 = [16 3];
%     distancePointPolygon(p0, square)
%     ans =
%          6
%
%   See also
%   polygons2d, points2d, distancePointPolyline, distancePointEdge,
%   projPointOnPolyline
%

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2009-04-30,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

% eventually copy first point at the end to ensure closed polygon
if sum(poly(end, :) == poly(1,:)) ~= 2
    poly = [poly; poly(1,:)];
end

% call to distancePointPolyline
minDist = distancePointPolyline(point, poly);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% eeg_interp_scalp3d_test - Script to test 3D interpolation of scalp potentials

clear

p = eeg_toolbox_defaults;
p = eeg_open(p);
p = elec_open(p);
p = mesh_open(p);

p = mesh_plot(p);

close all

scalpvert = p.mesh.data.vertices{3};

Mean = mean(scalpvert);
STD  = std(scalpvert);

Min = Mean - 3*STD;
Max = Mean + 3*STD;

delta = 0.05;
Xd = Min(1):delta:Max(1);
Yd = Min(2):delta:Max(2);
Zd = Min(3):delta:Max(3);

[x0,y0,z0] = meshgrid(Xd,Yd,Zd);

X0 = [x0(:), y0(:), z0(:)];

X = p.mesh.data.vertices{4};
v = p.volt.data(1,:)';

v0 = griddatan(X,v,X0);
v0 = reshape(v0, size(x0));

V = v0(:);
Vfinite = find(isfinite(V));
X0finite = X0(Vfinite,:);

trisurf(convhulln(X0finite),X0finite(:,1),X0finite(:,2),X0finite(:,3),Vfinite,'EdgeColor','none','FaceColor','interp');

fprintf('\ndone\n');
return

%X0scalpindices = dsearchn(X0,scalpvert);
X0scalpindices = dsearchn(X0finite,scalpvert);

newscalpvert(:,1) = x0(X0scalpindices);
newscalpvert(:,2) = y0(X0scalpindices);
newscalpvert(:,3) = z0(X0scalpindices);
newvoltage = v0(X0scalpindices);

trisurf(convhulln(newscalpvert),newscalpvert(:,1),newscalpvert(:,2),newscalpvert(:,3),newvoltage,'EdgeColor','none','FaceColor','interp');

return

Vplot = mean(v);

Hp = patch(isosurface(x0,y0,z0,v0,Vplot));
isonormals(x0,y0,z0,v0,Hp);
set(Hp,'FaceColor','red','EdgeColor','none');
view(3);
camlight;
lighting phong
%axis equal
title('Interpolated isosurface from scattered data')

author: PatternRecognition

repo: OpenBMI


text:
function r = bisect_characteristic ( x0, theta, characteristic )

%*****************************************************************************80
%
%% BISECT_CHARACTERISTIC: characteristic function transition surface.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X0(M,1), the coordinates of the base point inside the surface.
%
%    Input, real THETA, the coordinates of a direction from the base point.
%
%    Input, real y = CHARACTERISTIC ( m, n, x ), the handle for a function which
%    evaluates the characteristic for the object at N M-dimensional points X,
%    returning a 1 for points inside the object, and 0 for points outside.
%
%    Output, real R, the distance from X0 to the surface in direction THETA.
%

  x0 = x0(:);

  m = length ( x0 );
  n = 1;
%
%  Initially, X1 = X0, and so Y1 should be 1.
%
  x1(1:m,1) = x0(1:m,1);
  y1 = characteristic ( m, n, x1 );
%
%  Seek an exterior point.
%
  x2 = exterior_point_characteristic ( m, x0, theta, characteristic );
%
%  Carry out the bisection search in the interval [X1,X2].
%
  it_num = 0;

  while ( 1.0E-3 < sqrt ( norm ( x2 - x1 ) ) )

    x3 = ( x1 + x2 ) / 2.0;
    y3 = characteristic ( m, n, x3 );

    if ( y3 == 0 )
      x2 = x3;
    else
      x1 = x3;
    end

    it_num = it_num + 1;

    if ( 1000 < it_num )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'BISECT_CHARACTERISTIC - Fatal error!\n' );
      fprintf ( 1, '  Too many iterations.\n' );
      error ( 'BISECT_CHARACTERISTIC - Fatal error!' )
    end

  end
%
%  Measure the distance from X0 to the transition.
%  We estimate the transition to occur at the average of X1 and X2.
%
  x3 = ( x1 + x2 ) / 2.0;
  r = norm ( x3 - x0 );

  return
end


text:
%KNNR Trainable Nearest Neighbor Regression
%
%    Y = KNNR(X,K)
%    Y = X*KNNR([],K)
%    Y = X*KNNR(K)
%
% INPUT
%   X    Regression dataset, used for training
%   K    number of neighbors (default K=3)
%
% OUTPUT
%   Y    k-nearest neighbor regression
%
% DESCRIPTION
% Define a k-Nearest neighbor regression on dataset X.
%
% SEE ALSO (PRTools Guide)
% LINEARR, TESTR, PLOTR

% Copyright: D.M.J. Tax, D.M.J.Tax@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function y = knnr(varargin)

  mapname = 'KNN regression';
  argin = shiftargin(varargin,'scalar');
  argin = setdefaults(argin,[],3);

  if mapping_task(argin,'definition')

    y = define_mapping(argin,'untrained',mapname);

	elseif mapping_task(argin,'training')			% Train a mapping.

		[x,k] = deal(argin{:});
    [n,d] = size(x);
    W.x = +x;
    W.y = gettargets(x);
    W.k = k;
    y = prmapping(mfilename,'trained',W,1,d,1);
    y = setname(y,'k-nearest neighbor regression');

  else                                      % Evaluation

    [x,v] = deal(argin{1:2});
    w = getdata(v);
    [n,d] = size(x);
    D = distm(+x,w.x);
    [sD,I] = sort(D,2);
    if n==1
      out = mean(w.y(I(:,1:w.k)));
    else
      out = mean(w.y(I(:,1:w.k)),2);
    end
    y = setdat(x,out);

  end

author: marianux

repo: ecg-kit

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


text:
function varargout = eig_yalmip_internal(varargin)

switch class(varargin{1})

    case 'double'
        X = varargin{1}(:);
        n = sqrt(length(X));
        X = reshape(X,n,n);
        varargout{1} = eig((X+X')/2);

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.

        x = varargin{1};
        y = yalmip('definemulti',mfilename,x,[size(x,1) 1]);
        varargout{1} = y;

    case 'char' % YALMIP send 'graph' when it wants the epigraph or hypograph
        %   disp('The EIG operator is not supported in optimization problems.')
        %   error

        operator = CreateBasicOperator('callback');

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

        % Inofficial way to model several nonlinear variables in
        % one call
        % varargout{2}.models = getvariables(t):getvariables(t)+length(X)-1;

    otherwise
        error([upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [J, grad] = costFunction(theta, X, y)
%COSTFUNCTION Compute cost and gradient for logistic regression
%   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
%   parameter for logistic regression and the gradient of the cost
%   w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Note: grad should have the same dimensions as theta
%

J = (-1/m)*((y'*log(sigmoid(X*theta))) + ((1-y)'*log(1 - sigmoid(X*theta))));

grad = (1/m)*(X' * (sigmoid(X*theta) - y));

% =============================================================

end

author: UtkarshPathrabe