text:
classdef MOEADDE < ALGORITHM
% <multi/many> <real/integer>
% MOEA/D based on differential evolution
% delta --- 0.9 --- The probability of choosing parents locally
% nr    ---   2 --- Maximum number of solutions replaced by each offspring

%------------------------------- Reference --------------------------------
% H. Li and Q. Zhang, Multiobjective optimization problems with complicated
% Pareto sets, MOEA/D and NSGA-II, IEEE Transactions on Evolutionary
% Computation, 2009, 13(2): 284-302.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [delta,nr] = Algorithm.ParameterSet(0.9,2);

            %% Generate the weight vectors
            [W,Problem.N] = UniformPoint(Problem.N,Problem.M);
            T = ceil(Problem.N/10);

            %% Detect the neighbours of each solution
            B = pdist2(W,W);
            [~,B] = sort(B,2);
            B = B(:,1:T);

            %% Generate random population
            Population = Problem.Initialization();
            Z = min(Population.objs,[],1);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                % For each solution
                for i = 1 : Problem.N
                    % Choose the parents
                    if rand < delta
                        P = B(i,randperm(end));
                    else
                        P = randperm(Problem.N);
                    end

                    % Generate an offspring
                    Offspring = OperatorDE(Problem,Population(i),Population(P(1)),Population(P(2)));

                    % Update the ideal point
                    Z = min(Z,Offspring.obj);

                    % Update the solutions in P by Tchebycheff approach
                    g_old = max(abs(Population(P).objs-repmat(Z,length(P),1)).*W(P,:),[],2);
                    g_new = max(repmat(abs(Offspring.obj-Z),length(P),1).*W(P,:),[],2);
                    Population(P(find(g_old>=g_new,nr))) = Offspring;
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [img] = hyperConvert3d(img, h, w, numBands)
% HYPERCONVERT2D Converts an 2D matrix to a 3D data cube
% Converts a 2D matrix (p x N) to a 3D data cube (m x n x p)
% where N = m * n
%
% Usage
%   [M] = hyperConvert3d(M)
% Inputs
%   M - 2D data matrix (p x N)
% Outputs
%   M - 3D data cube (m x n x p)

if (ndims(img) ~= 2)
    error('Input image must be p x N.');
end

[numBands, N] = size(img);

if (1 == N)
    img = reshape(img, h, w);
else
    img = reshape(img.', h, w, numBands);
end

return;

author: BehnoodRasti


text:
function [conf, tx] = useBoundaryClassifier(bndinfo, X, classifier, featureset)
% [conf, tx] = useBoundaryClassifier(bndinfo, X, classifier, featureset)

tx = double(getBoundaryClassifierFeatures(bndinfo, X, (1:bndinfo.ne)));

if exist('featureset', 'var') && ~isempty(featureset)
    tx = getPartialFeatures(tx, [], featureset);
end

conf = test_boosted_dt_mc(classifier, tx);
conf = 1 ./ (1+exp(-conf));
conf = conf ./ repmat(sum(conf, 2), [1 size(conf, 2)]);

author: Cloud-CV

repo: object-proposals


text:
function updatefly(h,x,y,theta,maj,min)

% draw an isosceles triangle with center (x,y)
% with rotation theta
% with height maj*4
% with base min*4

if isstruct(x),
  trx = x;
  t = y;
  x = trx.x(t);
  y = trx.y(t);
  theta = trx.theta(t);
  maj = trx.a(t);
  min = trx.b(t);
end

if 0,

  ellipseupdate(h,maj*2,min*2,x,y,theta);

else

% isosceles triangle not yet rotated or centered
pts = [-maj*2,-min*2,
       -maj*2,min*2,
       maj*2,0];

% rotate
costheta = cos(theta);
sintheta = sin(theta);
R = [costheta,sintheta;-sintheta,costheta];
pts = pts*R;

% translate
pts(:,1) = pts(:,1) + x;
pts(:,2) = pts(:,2) + y;

% plot
set(h,'xdata',pts([1:3,1],1),'ydata',pts([1:3,1],2));

end

author: kristinbranson

repo: JAABA


text:
function op = smooth_logLLogistic(y)
% SMOOTH_LOGLLOGISTIC Log-likelihood function of a logistic: sum_i( y_i mu_i - log( 1+exp(mu_i) ) )
%   OP = SMOOTH_LOGLLOGISTIC( Y )
%   returns a function that computes the log-likelihood function
%   in a standard logistic regression model with independent entries. There
%   are two classes y_i = 0 and y_i = 1 with
%
%   prob(y_i = 1) = exp(mu_i)/(1 + exp(mu_i)
%
%   so that the log-likelihood is given by
%
%       log-likelihood(mu) = sum_i ( y_i mu_i - log(1+ exp(mu_i)) )
%
%   where mu is the parameter of the distribution (this is unknown,
%   so it is the variable), and Y is a vector of observations.

error(nargchk(1,1,nargin));
op = tfocs_smooth( @smooth_logLlogistic_impl);

function [ v, g ] = smooth_logLlogistic_impl( mu )

  if length(mu) == 1,
      mu = mu * ones(size(y));
  elseif size(mu) ~= size(y),
      error('Parameters and data must be of the same size'),
  end

  aux = 1 + exp(-abs(mu));
  v = tfocs_dot(y-1,mu.*(mu > 0)) ...
         + tfocs_dot(y,mu.*(mu < 0)) ...
              - tfocs_dot(ones(size(y)), log(aux));
  if nargout > 1,
      g = y - ((mu > 0) + (mu <= 0).*exp(mu))./aux;
  end
end

end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [u] = spm_erp_u(t,P,M)
% returns the [scalar] input for EEG models (Gaussian function)
% FORMAT [u] = spm_erp_u(t,P,M)
% t      - PST (seconds)
% P      - parameter structure
%   P.R  - scaling of [Gaussian] parameters
%
% u   - stimulus-related (subcortical) input
%
% See spm_fx_erp.m and spm_erp_priors.m
%__________________________________________________________________________
%
% David O, Friston KJ (2003) A neural mass model for MEG/EEG: coupling and
% neuronal dynamics. NeuroImage 20: 1743-1755
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_erp_u.m 7679 2019-10-24 15:54:07Z spm $

% preliminaries - check durations (ms)
%--------------------------------------------------------------------------
try
    if length(M.dur) ~= length(M.ons)
        M.dur = M.dur(1) + M.ons - M.ons;
    end
catch
    M.dur = 32 + M.ons - M.ons;
end

% check sustained input (0,1)
%--------------------------------------------------------------------------
try
    if length(M.sus) ~= length(M.ons)
        M.sus = M.sus(1) + M.ons - M.ons;
    end
catch
    M.sus = 0 + M.ons - M.ons;
end

% stimulus - Gaussian (subcortical) impulse
%--------------------------------------------------------------------------
nu    = length(M.ons);
u     = sparse(length(t),nu);
t     = t*1000;
for i = 1:nu

    % Gaussian bump function
    %----------------------------------------------------------------------
    delay  = M.ons(i) + 128*P.R(i,1);
    scale  = M.dur(i) * exp(P.R(i,2));
    U      = exp(-(t - delay).^2/(2*scale^2));

    % sustained inputs
    %----------------------------------------------------------------------
    try
        prop = M.sus(i)*exp(P.R(i,3));
    catch
        prop = M.sus(i);
    end
    U      = prop*cumsum(U)/sum(U) + U*(1 - prop);
    u(:,i) = 32*U;
end

author: spm


text:
function [valMIC] = applyMIC(pathMINE,variable,matData,jobName)
% -------------------------------------------------------------------------
% function [valMIC] = applyMIC(pathMINE,variable,matData,jobName)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes the Maximal Information Coefficient (MIC) between
% the feature contained in the vector 'variable' and all the features
% contained in the matrix 'matData'. MIC is computed using the executable
% MINE.jar, which can be downloaded at: <http://www.exploredata.net/>.
% -------------------------------------------------------------------------
% INPUTS:
% - pathMINE: Full path to the MINE.jar executable.
% - variable: Column cector of size [nInst X 1], where 'nInst' refers to the
%             number of instances for the given feature tested.
% - matData: Matrix of size [nInst X nFeat], where 'nFeat' is the number of
%            features to be tested against 'variable'.
% - jobName: String specifying the name of the job to be sent to MINE.jar.
% -------------------------------------------------------------------------
% OUTPUTS:
% -  valMIC: Column vector of size [nFeat X 1], specifying the MIC between
%            each feature in 'matData' and the 'variable' features. Entry
%            number in 'valMIC' corresponds to the column number in
%            'matData'.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: May 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
%
%    _______________________________________________________________
%
% MINE version 1.0.1d
% Copyright 2011 by David Reshef and Yakir Reshef.
%
% This application is licensed under a Creative Commons
% Attribution-NonCommercial-NoDerivs 3.0 Unported License.
%
% See http://creativecommons.org/licenses/by-nc-nd/3.0/ for more information.
% -------------------------------------------------------------------------

startpath=pwd;
cd(pathMINE)

% Writing the MIC input file
if isempty(strfind(jobName,'.csv'))
    jobName = [jobName,'.csv'];
end
id = (0:size(matData,2))';
data = [variable,matData]';
csvwrite(jobName,[id,data],0,0)

% Executing MIC
commandSys = ['java -jar MINE.jar ',jobName,' -masterVariable 0'];
[~,~] = system(commandSys);

% Reading the output from MIC
commandRead = [jobName,',mv=0,cv=0.0,B=n^0.6,Results.csv'];
micOutput = csvread(commandRead,1,1,[1 1 size(matData,2) 2]);

% Sorting the results
[val,ind] = sort(micOutput(:,1));
valMIC = val;
valMIC(1:end) = micOutput(ind(1:end),2);

% Cleaning up
delete(jobName)
delete(commandRead)
commandLast=[jobName,',mv=0,cv=0.0,B=n^0.6,Status.txt'];
delete(commandLast)

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [ b, m1, m2, r11, r12, r21, r22, seed ] = p36_p_val ( ...
  action, b, m1, m2, r11, r12, r21, r22, seed )

%*****************************************************************************80
%
%% P36_P_VAL sets or gets parameters for problem 36.
%
%  Discussion:
%
%    If ACTION is 'SET', the parameters are input values, and set
%    the internal values.
%
%    If ACTION is 'GET', the parametrs are output values, copied from
%    the internal values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION, specifies the action desired.
%    'SET' sets the parameters;
%    'GET' gets the current values.
%
%    Input/output, real B, M1, M2, R11, R12, R21, R22,
%    problem parameters.
%
%    Input/output, integer SEED, a seed for the random number generator.
%
  persistent b_save
  persistent m1_save
  persistent m2_save
  persistent r11_save
  persistent r12_save
  persistent r21_save
  persistent r22_save
  persistent seed_save

  if ( action(1) == 'G' || action(1) == 'g' )

    b = b_save;
    m1 = m1_save;
    m2 = m2_save;
    r11 = r11_save;
    r12 = r12_save;
    r21 = r21_save;
    r22 = r22_save;
    seed = seed_save;

  elseif ( action(1) == 'S' || action(1) == 's' )

    b_save = b;
    m1_save = m1;
    m2_save = m2;
    r11_save = r11;
    r12_save = r12;
    r21_save = r21;
    r22_save = r22;
    seed_save = seed;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function cnn_imagenet_minimal()
% CNN_IMAGENET_MINIMAL   Minimalistic demonstration of how to run an ImageNet CNN model

% Setup MatConvNet.
run(fullfile(fileparts(mfilename('fullpath')), ...
  '..', '..', 'matlab', 'vl_setupnn.m')) ;

% Download a pre-trained CNN from the web.
if ~exist('imagenet-vgg-f.mat', 'file')
  fprintf('Downloading the VGG-F model ... this may take a while\n') ;
  urlwrite('http://www.vlfeat.org/matconvnet/models/imagenet-vgg-f.mat', ...
    'imagenet-vgg-f.mat') ;
end

% Load the model and upgrade it to MatConvNet current version.
net = load('imagenet-vgg-f.mat') ;
net = vl_simplenn_tidy(net) ;

% Obtain and preprocess an image.
im = imread('peppers.png') ;
im_ = single(im) ; % note: 255 range
im_ = imresize(im_, net.meta.normalization.imageSize(1:2)) ;
im_ = im_ - net.meta.normalization.averageImage ;

% Run the CNN.
res = vl_simplenn(net, im_) ;

% Show the classification result.
scores = squeeze(gather(res(end).x)) ;
[bestScore, best] = max(scores) ;
figure(1) ; clf ; imagesc(im) ;
title(sprintf('%s (%d), score %.3f',...
   net.meta.classes.description{best}, best, bestScore)) ;

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function value = year_length_coptic ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_COPTIC returns the number of days in a Coptic year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of
%    days in the year.
%
  if ( year_is_leap_coptic ( y ) )
    value = 366;
  else
    value = 365;
  end

  return
end


text:
%   This directory contains the test suite for the mUnit test framework.  Before
%   running the test suite, do the following:
%
%   1. Make sure the mUnit test framework directory is on your path
%   2. Make sure the helper_classes subdirectory of the test directory is on
%      your path.
%   3. Make the test directory your current directory.
%
%   To run the test suite:
%
%       run(TestSuite())

%   Steven L. Eddins
%   Copyright 2008 The MathWorks

help Readme

author: ricedsp


text:
function error_frobenius = inverse_error ( n, a, b )

%*****************************************************************************80
%
%% INVERSE_ERROR determines the error in an inverse matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(N,N), the matrix.
%
%    Input, real B(N,N), the inverse.
%
%    Output, real ERROR_FROBENIUS, the Frobenius norm
%    of (A*B-I) + (B*A-I).
%
  error_frobenius = norm ( a * b - eye ( n ), 'fro' ) ...
                  + norm ( b * a - eye ( n ), 'fro' );

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   demo script to convert a closed surface to a binary image
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% preparation
% user must add the path of iso2mesh to matlab path list
% addpath('../');

%% load the sample data
load rat_head.mat

% first, generate a surface from the original image
% similar to demo_shortcuts_ex1.m

[node,face,regions,holes]=v2s(volimage,0.5,3);

node=sms(node,face(:,1:3),3,0.5); % apply 3 mesh smoothing

mdim=ceil(max(node)+1);
dstep=0.25;
zslice=15;
xrange=0:dstep:mdim(1);
yrange=0:dstep:mdim(2);
zrange=0:dstep:mdim(3);
img=surf2vol(node,face(:,1:3),xrange,yrange,zrange);

imagesc(squeeze(img(:,:,zslice))); % z=10

hold on

z0=zslice*dstep;
plane=[min(node(:,1)) min(node(:,2)) z0
       min(node(:,1)) max(node(:,2)) z0
       max(node(:,1)) min(node(:,2)) z0];

% run qmeshcut to get the cross-section information at z=mean(node(:,1))
% use the x-coordinates as the nodal values

[bcutpos,bcutvalue,bcutedges]=qmeshcut(face(:,1:3),node,node(:,1),plane);
[bcutpos,bcutedges]=removedupnodes(bcutpos,bcutedges);
bcutloop=extractloops(bcutedges);
bcutloop(isnan(bcutloop))=[]; % there can be multiple loops, remove the separators
plot(bcutpos(bcutloop,2)*(1/dstep),bcutpos(bcutloop,1)*(1/dstep),'w');

if(isoctavemesh)
  if(~exist('bwfill'))
    error('you need to install octave-image toolbox first');
  end
  img2=zeros(size(img),'uint8');
  for i=1:size(img,3)
    img2(:,:,i)=bwfill(img(:,:,i),'holes');
  end
  img2=img2+img;
else
  img2=imfill(img,'holes')+img;
end
figure;
imagesc(squeeze(img2(:,:,zslice))); % z=10
hold on;
plot(bcutpos(bcutloop,2)*(1/dstep),bcutpos(bcutloop,1)*(1/dstep),'y--');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [x, infos] = acc_online_mu_nmf(V, rank, in_options)
% Accelerated online non-negative matrix factorization (ONMF) algorithm.
%
% Inputs:
%       matrix      V
%       rank        rank
%       options     options
% Output:
%       w           solution of w
%       infos       information
%
% References:
%       S. S. Bucak, B. Gunsel,
%       "Incremental Subspace Learning via Non-negative Matrix Factorization,"
%       Pattern Recognition, 2009.
%
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai and H.Sakai on Feb. 12, 2017
%
% Change log:
%
%       Feb. 12, 2017 (Hiroyuki Kasai): Fixed algorithm.
%
%       May. 20, 2019 (Hiroyuki Kasai): Added initialization module.
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.rep_mode = 'fix';
    local_options.w_repeat = 1;
    local_options.h_repeat = 1;
    local_options.alpha    = 2;
    local_options.delta    = 0.1;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    Wt = init_factors.W;
    H = init_factors.H;

    % initialize
    method_name = 'ACC-Online-NMF';
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    %At = zeros(m, rank);
    %Bt = zeros(rank, rank);

    % initialize for this algorithm
    if strcmp(options.rep_mode, 'adaptive')
        K = m*n;
        rhoh = 1+(K+m*rank)/(n*(rank+1));
    end

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, Wt, H, [], options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('%s: Epoch = 0000, cost = %.16e, optgap = %.4e\n', method_name, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main outer loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % Reset sufficient statistic
        At = zeros(m, rank);
        Bt = zeros(rank, rank);

        % main inner loop
        for t = 1 : options.batch_size : n - 1

            % retrieve vt and ht
            vt = V(:, t:t+options.batch_size-1);
            ht = H(:, t:t+options.batch_size-1);

%             % uddate ht
%             Wtv = Wt.' * vt;
%             WtW = Wt.' * Wt;
%             for iii=1:h_repeat
%                 ht = ht .* (Wtv) ./ (WtW * ht);
%                 ht = ht + (ht<eps) .* eps;
%             end

            Wtv = Wt.' * vt;
            WtW = Wt.' * Wt;
            if strcmp(options.rep_mode, 'adaptive')
                gamma = 1;
                eps0 = 1;
                j = 1;
                rhoh_alpha = rhoh * options.alpha;

                %while j <= floor(1+rhoh*alpha) &&  gamma >= delta*eps0
                while j <= rhoh_alpha && gamma >= options.delta * eps0
                    ht0 = ht;
                    ht = ht .* (Wtv) ./ (WtW * ht);
                    ht = ht + (ht<eps) .* eps;
                    if j == 1
                        eps0 = norm(ht0 - ht);
                    end
                    gamma = norm(ht0 - ht);
                    j = j+1;
                end
            else
                for iii = 1 : options.h_repeat
                    ht = ht .* (Wtv) ./ (WtW * ht);
                    ht = ht + (ht<eps) .* eps;
                end
            end

            % update sufficient statistics
            At = At + vt *  ht';
            Bt = Bt + ht *  ht';

            % update W
            for iii = 1 : options.w_repeat
                Wt = Wt .* At ./ (Wt * Bt);
                Wt = Wt + (Wt<eps) .* eps;
            end

            % store new h
            H(:, t:t+options.batch_size-1) = ht;

            grad_calc_count = grad_calc_count + m * options.batch_size;
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, Wt, H, [], options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = Wt;
    x.H = H;
    x.R = [];

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [roi, roiNot] = dtiRoiClip(roi, rlClip, apClip, siClip)
%
% [roi, roiNot] = dtiRoiClip(roi, rlClip, apClip, siClip)
%
% Clips a dti ROI.
%
% HISTORY:
% 2005.01.05 RFD wrote it.

if(nargin<4) siClip = []; end
if(nargin<3) apClip = []; end
if(nargin==1)
    newName = [roi.name '_clip'];
    prompt = {'Left (-80) Right (+80) clip (blank for none):',...
          'Posterior (-120) Anterior (+80) clip (blank for none):',...
          'Inferior (-50) Superior (+90) clip (blank for none):',...
          'New ROI name:'};
    defAns = {'','','',newName};
    resp = inputdlg(prompt,'Clip Current ROI',1,defAns);
    if(isempty(resp))
        disp('User cancelled clip.');
        return;
    end
    rlClip = str2num(resp{1});
    apClip = str2num(resp{2});
    siClip = str2num(resp{3});
    roi.name = resp{4};
end

keep = ones(size(roi.coords,1),1);
if(~isempty(rlClip))
    keep = keep & (roi.coords(:,1)<rlClip(1) | roi.coords(:,1)>rlClip(2));
end
if(~isempty(apClip))
    keep = keep & (roi.coords(:,2)<apClip(1) | roi.coords(:,2)>apClip(2));
end
if(~isempty(siClip))
    keep = keep & (roi.coords(:,3)<siClip(1) | roi.coords(:,3)>siClip(2));
end
coords = roi.coords;
roi.coords = coords(keep,:);
if(nargout>1)
    roiNot = roi;
    roiNot.coords = coords(~keep,:);
    roiNot.name = [roi.name '_NOT'];
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function linplus_test48 ( )

%*****************************************************************************80
%
%% TEST48 tests R8PBU_ML.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  mu = 3;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST48\n' );
  fprintf ( 1, '  R8PBU_ML computes A*x\n' );
  fprintf ( 1, '    where A has been factored by R8PBU_FA.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
  fprintf ( 1, '  Upper bandwidth MU = %d\n', mu );
%
%  Set the matrix.
%
  [ a, seed ] = r8pbu_random ( n, mu, seed );
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
  b = r8pbu_mxv ( n, mu, a, x );
%
%  Factor the matrix.
%
  [ a_lu, info ] = r8pbu_fa ( n, mu, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Fatal error!\n' );
    fprintf ( 1, '  R8PBU_FA declares the matrix is singular!\n' );
    fprintf ( 1, '  The value of INFO is %d\n', info );
    return
  end
%
%  Now multiply factored matrix times solution to get right hand side again.
%
  b2 = r8pbu_ml ( n, mu, a_lu, x );

  r8vec2_print_some ( n, b, b2, 10, '  A*x and PLU*x' );

  return
end


text:
 function [alpha, beta, ok] = nufft_alpha(N, J, K, alpha, beta)
%	Determine alpha and beta, as associated with the scaling factors
%	for min-max interpolation, from arguments.
%	in
%		alpha	[L,1]	Fourier series coefficients of scaling factors
%				or a string (see below)
%		beta		scale gamma=2pi/K by this in Fourier series
%				typically is K/N (me) or 0.5 (Liu)
%	out
%		alpha,beta
%
%	Copyright 2001-3-30	Jeff Fessler	The University of Michigan

%	if no arguments
if nargin < 3
	help(mfilename)
	error(mfilename)
return
end

if ~isvar('alpha') | isempty(alpha)
	alpha = [1];	% default Fourier series coefficients of scaling factors
end
if ~isvar('beta') | isempty(beta)
	beta = 0.5;	% default is Liu version for now
end

%
%	see if 'best' alpha is desired
%
if ischar(alpha)
	if streq(alpha, 'uniform')
		alpha = [1];
		beta = 0.5;
	else
		if streq(alpha, 'best')
			L = 0;
		elseif streq(alpha, 'best,L=1')
			L = 1;
		elseif streq(alpha, 'best,L=2')
			L = 2;
		else
			error 'unknown alpha argument'
		end
		[alpha, beta, ok] = nufft_best_alpha(J, L, K/N);
		if ~ok
			tmp = 'optimal alpha unknown for J=%d, K/N=%g, L=%d';
			warning(sprintf(tmp, J, K/N, L))
		end
	end
end

author: thomaskuestner


text:
function test_suite = test_crossvalidate
% tests for test_crossvalidate
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_crossvalidate_basics
    classifier=@cosmo_classify_nn;
    randint=@()ceil(rand()*5+5);

    ds=cosmo_synthetic_dataset('ntargets',randint(),...
                               'nchunks',randint(),...
                               'nreps',randint(),...
                               'seed',0);    % random data
    nsamples=size(ds.samples,1);
    nfolds=randint();

    partitions=struct();
    partitions.train_indices=cell(nfolds,1);
    partitions.test_indices=cell(nfolds,1);

    train_size=ceil(nsamples*(rand()*.5+.25));

    pred=NaN(nsamples,nfolds);

    for fold=1:nfolds
        all_idx=randperm(nsamples);
        train_idx=all_idx(1:train_size);
        test_idx=all_idx((train_size+1):end);

        partitions.train_indices{fold}=train_idx;
        partitions.test_indices{fold}=test_idx;

        pred(test_idx,fold)=classifier(ds.samples(train_idx,:),...
                                        ds.sa.targets(train_idx),...
                                        ds.samples(test_idx,:));
    end

    pred_msk=~isnan(pred);
    is_correct=bsxfun(@eq,ds.sa.targets,pred) & pred_msk;
    acc=sum(is_correct(:))/sum(pred_msk(:));

    opt=struct();
    opt.check_partitions=false;

    [res_pred,res_acc]=cosmo_crossvalidate(ds,classifier,partitions,opt);
    assertEqual(res_pred,pred);
    assertElementsAlmostEqual(res_acc,acc);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function h = drawLine3d(lin, varargin)
%DRAWLINE3D Draw a 3D line clipped by the current axes.
%
%   drawLine3d(LINE) draws the line LINE on the current axis, by clipping
%   with the current axis.
%
%   drawLine3d(LINE, PARAM, VALUE) accepts parameter/value pairs, like
%   for plot function. Color of the line can also be given as a single
%   parameter.
%
%   drawLine3d(AX,...) plots into AX instead of GCA.
%
%   H = drawLine3d(...) returns a handle to the created line object.
%   If the line is not clipped by the axis, function returns -1.
%
%   Example
%     % draw a sphere together with the three main axes
%     figure; hold on;
%     drawSphere([40 30 20 30]);
%     view(3); axis equal; axis([-20 80 -20 80 -20 80])
%     drawLine3d([0 0 0   1 0 0], 'k');
%     drawLine3d([0 0 0   0 1 0], 'k');
%     drawLine3d([0 0 0   0 0 1], 'k');
%     light;
%
%
%   See also:
%     lines3d, createLine3d, clipLine3d
%

% ---------
% author : David Legland
% INRA - TPV URPOI - BIA IMASTE
% created the 17/02/2005.

% Parse and check inputs
isLine3d = @(x) validateattributes(x,{'numeric'},...
    {'nonempty','nonnan','real','finite','size',[nan,6]});
defOpts.Color = 'b';
[hAx, lin, varargin] = ...
    parseDrawInput(lin, isLine3d, 'line', defOpts, varargin{:});

% extract limits of the bounding box
box = [get(gca, 'xlim') get(gca, 'ylim') get(gca, 'zlim')];

% clip the line with the limits of the current axis
edge = clipLine3d(lin, box);

% draw the clipped line
if sum(isnan(edge)) == 0
    hh = drawEdge3d(hAx, edge);
    if ~isempty(varargin)
        set(hh, varargin{:});
    end
else
    hh = [];
end

% process output
if nargout > 0
    h = hh;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function net= addLayers(net, opts, dbTrain)

    methodOpts= strsplit(opts.method, '_');
    [~, sz]= relja_netOutputDim(net);
    D= sz(3);

    if ismember('preL2', methodOpts)
        % normalize feature-wise
        net.layers{end+1}= struct('type', 'normalize', 'name', 'preL2', ...
            'param', [2*D, 1e-12, 1, 0.5], 'precious', 0);
        methodOpts= removeOpt(methodOpts, 'preL2');
        doPreL2= true;
    else
        doPreL2= false;
    end

    if ismember('max', methodOpts)
        methodOpts= removeOpt(methodOpts, 'max');
        net.layers{end+1}= layerTotalMaxPool('max:core');

    elseif ismember('avg', methodOpts)
        methodOpts= removeOpt(methodOpts, 'avg');
        net.layers{end+1}= layerTotalAvgPool('avg:core');

    elseif any( ismember( {'vlad', 'vladv2'}, methodOpts) )

        if doPreL2
            L2str= '_preL2';
        else
            L2str= '';
        end

        whichDesc= sprintf('%s_%s%s', opts.netID, opts.layerName, L2str);

        k= 64;
        paths= localPaths();
        trainDescFn= sprintf('%s%s_%s_traindescs.mat', paths.initData, dbTrain.name, whichDesc);
        clstFn= sprintf('%s%s_%s_k%03d_clst.mat', paths.initData, dbTrain.name, whichDesc, k);

        clsts= getClusters(net, opts, clstFn, k, dbTrain, trainDescFn);

        load( trainDescFn, 'trainDescs');
        load( clstFn, 'clsts');
        net.meta.sessionID= sprintf('%s_%s', net.meta.sessionID, dbTrain.name);

        % --- VLAD layer

        if ismember('vladv2', methodOpts)
            methodOpts= removeOpt(methodOpts, 'vladv2');

            % set alpha for sparsity
            [~, dsSq]= yael_nn(clsts, trainDescs, 2); clear trainDescs;
            alpha= -log(0.01)/mean( dsSq(2,:)-dsSq(1,:) ); clear dsSq;

            net.layers{end+1}= layerVLADv2('vlad:core');
            net.layers{end}= net.layers{end}.constructor({alpha*2*clsts, -alpha*sum(clsts.^2,1), -clsts});

        elseif ismember('vlad', methodOpts)
            % see comments on vladv2 vs vlad in the README_more.md

            methodOpts= removeOpt(methodOpts, 'vlad');

            % set alpha for sparsity
            clstsAssign= relja_l2normalize_col(clsts);
            dots= sort(clstsAssign'*trainDescs, 1, 'descend'); clear trainDescs;
            alpha= -log(0.01)/mean( dots(1,:) - dots(2,:) ); clear dots;

            net.layers{end+1}= layerVLAD('vlad:core');
            net.layers{end}= net.layers{end}.constructor({alpha*clstsAssign, clsts});

        else
            error('Unsupported method "%s"', opts.method);
        end

        if ismember('intra', methodOpts)
            % --- intra-normalization
            net.layers{end+1}= struct('type', 'normalize', 'name', 'vlad:intranorm', ...
                'param', [2*D, 1e-12, 1, 0.5], 'precious', 0);
            methodOpts= removeOpt(methodOpts, 'intra');
        end

    else
        error('Unsupported method "%s"', opts.method);
    end

    % --- final normalization
    net.layers{end+1}= layerWholeL2Normalize('postL2');

    % --- check if all options are used
    if ~isempty(methodOpts)
        error('Unsupported options (method=%s): %s', opts.method, strjoin(methodOpts, ', '));
    end

    net.meta.sessionID= sprintf('%s_%s', net.meta.sessionID, opts.method);
    net.meta.epoch= 0;

end

function opts= removeOpt(opts, optName)
    opts(ismember(opts, optName))= [];
end

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
% Author: Ricardo Baptista and Matthias Poloczek
% Date:   June 2018
%
% See LICENSE.md for copyright information
%

function KL = KL_divergence_ising(Theta_P, moments, x)
% KL_divergence_ising: Function evaluates the KL divergence objective
% for Ising Models

n_vars = size(Theta_P,1);

% Generate all binary vectors
bin_vals = dec2bin(0:2^n_vars-1)-'0';
bin_vals(bin_vals == 0) = -1;
n_vectors = size(bin_vals,1);

% Compute normalizing constant for P
P_vals = zeros(n_vectors,1);
for i=1:n_vectors
	P_vals(i) = exp(bin_vals(i,:)*Theta_P*bin_vals(i,:)');
end
Zp = sum(P_vals);

% Run computation for each x
n_xvals = size(x,1);
KL = zeros(n_xvals,1);

for j=1:n_xvals

    % Apply elementwise masking to Theta
    Theta_Q = tril(Theta_P,-1);
    nnz_Q   = find(Theta_Q);
    Theta_Q(nnz_Q) = Theta_Q(nnz_Q).*(x(j,:))';
    Theta_Q = Theta_Q + Theta_Q';

    % Compute normalizing constant for Q
    Q_vals = zeros(n_vectors,1);
    for i=1:n_vectors
        Q_vals(i) = exp(bin_vals(i,:)*Theta_Q*bin_vals(i,:)');
    end
    Zq = sum(Q_vals);

    % compute KL
    KL(j) = sum(sum((Theta_P - Theta_Q).*moments)) + log(Zq) - log(Zp);

end

end

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
function plotData(X, y)
%PLOTDATA Plots the data points X and y into a new figure
%   PLOTDATA(x,y) plots the data points with + for the positive examples
%   and o for the negative examples. X is assumed to be a Mx2 matrix.

% Create New Figure
figure; hold on;

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the positive and negative examples on a
%               2D plot, using the option 'k+' for the positive
%               examples and 'ko' for the negative examples.
%

% Find Indices of Positive and Negative Examples
pos = find(y==1); neg = find(y == 0);

% Plot Examples
plot(X(pos, 1), X(pos, 2), 'k+', 'LineWidth', 2, 'MarkerSize', 7);
plot(X(neg, 1), X(neg, 2), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7);

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

hold off;

end

author: benoitvallon

repo: coursera-machine-learning


text:
function [ok,msg] = checkOptiSol(optObj,tol)
%CHECKOPTISOL  Check the solution to an optimization for errors
%
%   Called by OPTI / checkSol

%   Copyright (C) 2011 Jonathan Currie (I2C2)

ok = 1;
msg = sprintf('Solver Status:');

if(nargin < 2), tol = 1e-6; end

lintol = tol;
quadtol = tol;
nllintol = tol;
sdtol = tol;
inttol = 1e-5;

%Check Exit Flag
switch(optObj.ef)
    case 1
        msg = sprintf('%s Solved\n',msg);
    case 0
        msg = sprintf('%s Iterations / Time / Nodes Exceeded\n',msg);
        ok = 0;
    case -1
        msg = sprintf('%s Infeasible\n',msg);
        ok = 0;
    otherwise
        msg = sprintf('%s Solver Error\n',msg);
        ok = 0;
end

%Get Problem Properties
prob = optObj.prob;
x = optObj.sol; [r,c] = size(x);
if(c > r), x = x'; end

%Convert SeDuMi problem to OPTI if specified
if(~isempty(prob.sdcone) && isstruct(prob.sdcone))
    prob = sedumi2opti(prob);
end

%Check Linear Constraints
if(~isempty(prob.rl))
    %Check Row Bounds
    if(~isempty(prob.A))
        cval = prob.A*x;
        [msg,ok] = checkRowCon(cval,prob.rl,prob.ru,'Linear',lintol,msg,ok);
    end
else
    %Check Inequality
    if(~isempty(prob.A))
        err = prob.A*x - prob.b;
        chk = err > lintol;
        if(any(chk))
            ok = 0;
            msg = sprintf('%s\nLinear <= Inequality Constraint(s) Broken: [tol %g]\n',msg,lintol);
            msg = solErrMsg(msg,true(size(chk)),chk,err);
        end
    end
    %Check Equality
    if(~isempty(prob.Aeq))
        err = abs(prob.Aeq*x - prob.beq);
        chk = err > lintol;
        if(any(chk))
            ok = 0;
            msg = sprintf('%s\nLinear == Equality Constraint(s) Broken: [tol %g]\n',msg,lintol);
            msg = solErrMsg(msg,true(size(chk)),chk,err);
        end
    end
end
%Check Bounds
if(~isempty(prob.lb))
    err = abs(x - prob.lb);
    chk = x < (prob.lb-lintol);
    if(any(chk))
        ok = 0;
        msg = sprintf('%s\nDecision Variable Lower Bound(s) Broken: [tol %g]\n',msg,lintol);
        msg = solErrMsg(msg,true(size(chk)),chk,err);
    end
end
if(~isempty(prob.ub))
    err = abs(x - prob.ub);
    chk = x > (prob.ub+lintol);
    if(any(chk))
        ok = 0;
        msg = sprintf('%s\nDecision Variable Upper Bound(s) Broken: [tol %g]\n',msg,lintol);
        msg = solErrMsg(msg,true(size(chk)),chk,err);
    end
end

%Check Quadratic Constraints
if(~isempty(prob.Q))
    if(iscell(prob.Q)) %multiple
        %Evaluate each constraint
        cval = zeros(size(prob.qrl));
        for i = 1:length(prob.qrl)
           cval(i) = x'*prob.Q{i}*x + prob.l(:,i)'*x;
        end
        %Check all constraints
        [msg,ok] = checkRowCon(cval,prob.qrl,prob.qru,'Quadratic',quadtol,msg,ok);
    else %single
        cval = x'*prob.Q*x + prob.l'*x;
        [msg,ok] = checkRowCon(cval,prob.qrl,prob.qru,'Quadratic',quadtol,msg,ok);
    end
end

%Check Integer Constraints
if(any(prob.int.str == 'I'))
    idx = prob.int.str == 'I';
    err = abs(x - round(x));
    chk = err > inttol;
    if(any(chk(idx)))
        ok = 0;
        msg = sprintf('%s\nDecision Variable Integer Constraint(s) Broken: [tol %g]\n',msg,inttol);
        msg = solErrMsg(msg,idx,chk,err);
    end
end
%Check Binary Constraints
if(any(prob.int.str == 'B'))
    idx = prob.int.str == 'B';
    err = abs(max([x - round(x), x - 1, 0 - x],[],2));
    chk = err > inttol;
    if(any(chk(idx)))
        ok = 0;
        msg = sprintf('%s\nDecision Variable Binary Constraint(s) Broken: [tol %g]\n',msg,inttol);
        msg = solErrMsg(msg,idx,chk,err);
    end
end

%Check Semidefinite Constraints
if(~isempty(prob.sdcone))
    if(iscell(prob.sdcone))
        err = zeros(length(prob.sdcone),1);
        chk = zeros(length(prob.sdcone),1);
        for i = 1:length(prob.sdcone)
            cval = evalSDCone(prob.sdcone{i},x);
            err(i) = -min(eig(cval));
            chk(i) = err(i) > sdtol;
        end
    else
        cval = evalSDCone(prob.sdcone,x);
        err = -min(eig(cval));
        chk = err > sdtol;
    end
    if(any(chk))
        ok = 0;
        msg = sprintf('%s\nSemidefinite Constraint(s) Broken: [tol %g]\n',msg,sdtol);
        msg = solErrMsg(msg,true(size(chk)),chk,err);
    end
end

%Check Nonlinear Constraints
if(~isempty(prob.nlcon))
    %Convert general constraints to row (if not already)
    if(~isempty(prob.nlrhs)), prob = nmix2row(prob); end
    %Get cval
    cval = prob.nlcon(x);
    %Run check
    [msg,ok] = checkRowCon(cval,prob.cl,prob.cu,'Nonlinear',nllintol,msg,ok);
end

function [msg,ok] = checkRowCon(cval,rl,ru,str,tol,msg,ok)
%Indices
eq = rl == ru; neq = ~eq;
ile = isfinite(ru) & neq;
ige = isfinite(rl) & neq;
if(any(ile))
    err = (cval - ru);
    chk = err > tol;
    if(any(chk(ile)))
        ok = 0;
        msg = sprintf('%s\n%s <= Inequality Constraint(s) Broken: [tol %g]\n',msg,str,tol);
        msg = solErrMsg(msg,ile,chk,err);
    end
end
if(any(ige))
    err = (rl - cval);
    chk = err > tol;
    if(any(chk(ige)))
        ok = 0;
        msg = sprintf('%s\n%s >= Inequality Constraint(s) Broken: [tol %g]\n',msg,str,tol);
        msg = solErrMsg(msg,ige,chk,err);
    end
end
if(any(eq))
    err = abs(cval - rl);
    chk = err > tol;
    if(any(chk(eq)))
        ok = 0;
        msg = sprintf('%s\n%s == Equality Constraint(s) Broken: [tol %g]\n',msg,str,tol);
        msg = solErrMsg(msg,eq,chk,err);
    end
end

function msg = solErrMsg(msg,indx,chk,err)
for i = 1:length(chk)
    if(indx(i) == true && chk(i) == true)
        msg = sprintf('%s #%-5d Error: %-10g\n',msg,i,err(i));
    end
end

function cval = evalSDCone(cone,x)
dim = sqrt(size(cone,1));
cval = -cone(:,1);
for i = 1:length(x)
    cval = cval + x(i).*cone(:,i+1);
end
cval = reshape(cval,dim,dim);


text:
function NAMES = fieldnamesr(S,varargin)
%FIELDNAMESR Get structure field names in recursive manner.
%
%   NAMES = FIELDNAMESR(S) returns a cell array of strings containing the
%   structure field names associated with s, the structure field names
%   of any structures which are fields of s, any structures which are
%   fields of fields of s, and so on.
%
%   NAMES = FIELDNAMESR(S,DEPTH) is an optional field which allows the depth
%   of the search to be defined. Default is -1, which does not limit the
%   search by depth. A depth of 1 will return only the field names of s
%   (behaving like FIELDNAMES). A depth of 2 will return those field
%   names, as well as the field names of any structure which is a field of
%   s.
%
%   NAMES = FIELDNAMESR(...,'full') returns the names of fields which are
%   structures, as well as the contents of those structures, as separate
%   cells. This is unlike the default where a structure-tree is returned.
%
%   NAMES = FIELDNAMESR(...,'prefix') returns the names of the fields
%   prefixed by the name of the input structure.
%
%   NAMES = FIELDNAMESR(...,'struct') returns only the names of fields
%   which are structures.
%
%   See also FIELDNAMES, ISFIELD, GETFIELD, SETFIELD, ORDERFIELDS, RMFIELD.

%   Developed in MATLAB 7.13.0.594 (2011b).
%   By AATJ (adam.tudorjones@pharm.ox.ac.uk). 2011-10-11. Released under
%   the BSD license.

%Set optional arguments to default settings, if they are not set by the
%caller.
if size(varargin,1) == 0
    optargs = {-1} ;
    optargs(1:length(varargin)) = varargin ;
    depth = optargs{:} ;

    full = [0] ; [prefix] = [0] ; [struct] = [0] ;
else
    if any(strcmp(varargin,'full') == 1)
        full = 1 ;
    else
        full = 0 ;
    end

    if any(strcmp(varargin,'prefix') == 1)
        prefix = 1 ;
    else
        prefix = 0 ;
    end

    if any(strcmp(varargin,'struct') == 1)
        struct = 1 ;
    else
        struct = 0 ;
    end

    if any(cellfun(@isnumeric,varargin) == 1)
        depth = varargin{find(cellfun(@isnumeric,varargin))} ;
    else
        depth = -1 ;
    end
end

%Return fieldnames of input structure, prefix these with "S.".
NAMES = cellfun(@(x) strcat('S.',x),fieldnames(S),'UniformOutput',false) ;

%Set state counters to initial values (k terminates recursive loop, g makes
%recursive loop behave in a different way.
k = 1 ; g = 0 ;

fndstruct = {} ;

%k is ~0 while all fields have not been searched to exhaustion or to
%specified depth.
while k ~= 0
    fndtemp = {} ;

    k = length(NAMES) ;

    %g set to 1 prevents fieldnames from being added to output NAMES.
    %Useful when determining whether fields at the lowest specified depth
    %are structures, without adding their child fieldnames (i.e. at
    %specified depth + 1) to NAMES output.
    if depth == 1
        g = 1 ;
    end

    for i = 1:length(NAMES)
        %If the current fieldname is a structure, find its child
        %fieldnames, add to NAMES if not at specified depth (g = 0). Add to
        %fndstruct (list of structures).
        if isstruct(eval(NAMES{i})) == 1
            if g ~= 1
                fndtemp2 = fieldnames(eval(NAMES{i})) ;
                fndtemp2 = cellfun(@(x) strcat(sprintf('%s.'            ...
                    ,NAMES{i}),x),fndtemp2,'UniformOutput',false) ;
                fndtemp = cat(1,fndtemp,fndtemp2) ;
            elseif g == 1
                fndtemp = cat(1,fndtemp,NAMES{i}) ;
                k = k - 1 ;
            end
            fndstruct = cat(1,fndstruct,NAMES{i}) ;
        else
            fndtemp = cat(1,fndtemp,NAMES{i}) ;
            k = k - 1 ;
        end
    end

    NAMES = fndtemp ;

    %If we have reached depth, stop recording children fieldnames to NAMES
    %output, but determine whether fields at final depth are structures or
    %not (g). After this, terminate loop by setting k to 0.
    if depth ~= -1                                                      ...
            && any(cellfun(@(x) size(find(x == '.'),2),NAMES) == depth)
        g = 1 ;
    elseif depth ~= -1                                                  ...
            && any(cellfun(@(x) size(find(x == '.'),2),NAMES) > depth)
        k = 0 ;
    end
end

%Return names of fields which are structures, as well as fields which are
%not structures, if 'full' optional argument is set.
if full == 1
    for i = 1:length(fndstruct)

        %If depth is specified, add structure names to appropriate depth to
        %output NAMES.
        if depth == -1 || size(find(fndstruct{i} == '.'),2) <= depth-1
            structi = find(cellfun(@(x) isempty(x) == 0                 ...
                ,strfind(NAMES,fndstruct{i})),1) ;

            %If the current structure name is related to the first field of
            %NAMES, add structure name to NAMES in a particular way, else
            %add it in another way.
            if structi > 1
                NAMES = cat(1,NAMES(1:(structi-1)),fndstruct{i}         ...
                    ,NAMES(structi:end)) ;
            elseif isempty(structi)
                error ('AATJ:fieldnamesr,NotFindStructure'              ...
                    ,'Could not find structure name in NAMES') ;
            else
                NAMES = cat(1,fndstruct{i},NAMES) ;
            end
        end
    end
end

%Return only fields which are structures, if optional argument is defined.
%If 'full' is not set, do not include both parent and child structures.
if struct == 1
    if full == 0
        fndstruct2 = {} ;
        for i = 1:length(fndstruct)
            if isempty(cell2mat(strfind(fndstruct,strcat(fndstruct{i},'.')))) == 1
                fndstruct2(end+1,1) = fndstruct(i) ;
            end
        end
        fndstruct = fndstruct2 ;
    end
    NAMES = fndstruct ;
end

%Prefix input structure name on all fields of output cell array NAME if
%user set this option.
if prefix == 1
    structname = inputname(1) ;
    NAMES = cellfun(@(x) strcat(sprintf('%s.',structname),x(3:end))     ...
        ,NAMES,'UniformOutput',false) ;
else
    NAMES = cellfun(@(x) sprintf('%s',x(3:end)),NAMES,'UniformOutput'   ...
        ,false) ;
end

author: translationalneuromodeling

repo: tapas


text:
function alpha = p19_alpha ( )

%*****************************************************************************80
%
%% P19_ALPHA returns ALPHA for problem 19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real ALPHA, the value of ALPHA.
%
  alpha = 0.0;

  return
end


text:
function [cfg, artifact] = ft_artifact_tms(cfg, data)

% FT_ARTIFACT_TMS reads the data segments of interest from file and identifies artefacts in
% EEG recordings that were done during TMS stimulation.
%
% Use as
%   [cfg, artifact] = ft_artifact_tms(cfg)
% with the configuration options
%   cfg.dataset     = string with the filename
% or
%   cfg.headerfile  = string with the filename
%   cfg.datafile    = string with the filename
% and optionally
%   cfg.headerformat
%   cfg.dataformat
%
% Alternatively you can use it as
%   [cfg, artifact] = ft_artifact_tms(cfg, data)
% where the input data is a structure as obtained from FT_PREPROCESSING.
%
% In both cases the configuration should also contain
%   cfg.trl         = structure that defines the data segments of interest, see FT_DEFINETRIAL
%   cfg.continuous  = 'yes' or 'no' whether the file contains continuous data (default = 'yes')
% and
%   cfg.method      = 'detect' or 'marker', see below.
%   cfg.prestim     = scalar, time in seconds prior to onset of detected event to mark as artifactual (default = 0.005 seconds)
%   cfg.poststim    = scalar, time in seconds post onset of detected even to mark as artifactual (default = 0.010 seconds)
%
% The different methods are described in detail below.
%
% With cfg.method='detect', TMS-artifact are detected on basis of transient
% high-amplidude gradients that are typical for TMS-pulses. The data is preprocessed
% (again) with the following settings, which are optimal for identifying TMS-pulses.
% Artifacts are identified by means of thresholding the z-transformed value of the
% preprocessed data. This method acts as a wrapper around FT_ARTIFACT_ZVALUE.
%   cfg.artfctdef.tms.derivative  = 'yes'
%   cfg.artfctdef.tms.channel     = Nx1 cell-array with selection of channels, see FT_CHANNELSELECTION for details
%   cfg.artfctdef.tms.cutoff      = z-value at which to threshold (default = 4)
%   cfg.artfctdef.tms.trlpadding  = 0.1
%   cfg.artfctdef.tms.fltpadding  = 0.1
%   cfg.artfctdef.tms.artpadding  = 0.01
% Be aware that if one artifact falls within this specified range of another artifact, both
% artifact will be counted as one. Depending on cfg.prestim and cfg.poststim you may not mark
% enough data as artifactual.
%
% With cfg.method='marker', TMS-artifact onsets and offsets are based on markers/triggers that
% are written into the EEG dataset. This method acts as a wrapper around FT_DEFINETRIAL to
% determine on- and offsets of TMS pulses by reading markers in the EEG.
%   cfg.trialfun            = function name, see below (default = 'ft_trialfun_general')
%   cfg.trialdef.eventtype  = 'string'
%   cfg.trialdef.eventvalue = number, string or list with numbers or strings
% The cfg.trialfun option is a string containing the name of a function that you wrote
% yourself and that FT_ARTIFACT_TMS will call. The function should take the cfg-structure as
% input and should give a NxM matrix with M equal to or larger than 3) in the same format as
% "trl" as the output. You can add extra custom fields to the configuration structure to
% pass as arguments to your own trialfun.  Furthermore, inside the trialfun you can use the
% FT_READ_EVENT function to get the event information from your data file.
%
% The output argument "artifact" is a Nx2 matrix comparable to the
% "trl" matrix of FT_DEFINETRIAL. The first column of which specifying the
% beginsamples of an artifact period, the second column contains the
% endsamples of the artifactperiods.
%
% To facilitate data-handling and distributed computing with the peer-to-peer
% module, this function has the following option:
%   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_REJECTARTIFACT, FT_ARTIFACT_CLIP, FT_ARTIFACT_ECG, FT_ARTIFACT_EOG,
% FT_ARTIFACT_JUMP, FT_ARTIFACT_MUSCLE, FT_ARTIFACT_THRESHOLD, FT_ARTIFACT_ZVALUE

% Copyright (C) 2003-2011, Jan-Mathijs Schoffelen & Robert Oostenveld
%
% 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$

% 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 provenance is not needed because just a call to ft_artifact_zvalue
% ft_preamble loadvar data is not needed because ft_artifact_zvalue will do this

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed',    {'datatype', 'continuous'});
cfg = ft_checkconfig(cfg, 'renamedval', {'continuous', 'continuous', 'yes'});
cfg = ft_checkconfig(cfg, 'required', 'method');
cfg = ft_checkconfig(cfg, 'allowedval', {'method', 'detect', 'marker'});

% set default rejection parameters
if ~isfield(cfg, 'artfctdef'),                       cfg.artfctdef                   = [];        end
if ~isfield(cfg, 'method'),                          cfg.method                      = 'detect';  end
if ~isfield(cfg.artfctdef, 'tms'),                   cfg.artfctdef.tms               = [];        end
if ~isfield(cfg, 'prestim'),                         cfg.prestim                     = 0.005;     end
if ~isfield(cfg, 'poststim'),                        cfg.poststim                    = 0.010;     end

if isfield(cfg.artfctdef.tms, 'artifact')
  fprintf('tms artifact detection has already been done, retaining artifacts\n');
  artifact = cfg.artfctdef.tms.artifact;
  return
end

switch cfg.method
  case 'detect'
    % settings for preprocessing
    if ~isfield(cfg.artfctdef.tms, 'derivative'),   cfg.artfctdef.tms.derivative   = 'yes';    end
    % settings for the zvalue subfunction
    if ~isfield(cfg.artfctdef.tms, 'method'),     cfg.artfctdef.tms.method      = 'zvalue';    end
    if ~isfield(cfg.artfctdef.tms, 'channel'),    cfg.artfctdef.tms.channel     = 'all';       end
    if ~isfield(cfg.artfctdef.tms, 'trlpadding'), cfg.artfctdef.tms.trlpadding  = 0.1;         end
    if ~isfield(cfg.artfctdef.tms, 'fltpadding'), cfg.artfctdef.tms.fltpadding  = 0.1;         end
    if ~isfield(cfg.artfctdef.tms, 'artpadding'), cfg.artfctdef.tms.artpadding  = 0.01;        end
    if ~isfield(cfg.artfctdef.tms, 'cutoff'),     cfg.artfctdef.tms.cutoff      = 4;           end
    % construct a temporary configuration that can be passed onto artifact_zvalue
    tmpcfg                  = [];
    tmpcfg.trl              = cfg.trl;
    tmpcfg.artfctdef.zvalue = cfg.artfctdef.tms;
    if isfield(cfg, 'continuous'),   tmpcfg.continuous       = cfg.continuous;    end
    if isfield(cfg, 'dataformat'),   tmpcfg.dataformat       = cfg.dataformat;    end
    if isfield(cfg, 'headerformat'), tmpcfg.headerformat     = cfg.headerformat;  end
    % call the zvalue artifact detection function

    % the data is either passed into the function by the user or read from file with cfg.inputfile
    hasdata = exist('data', 'var');

    if hasdata
      % read the header
      cfg = ft_checkconfig(cfg, 'forbidden', {'dataset', 'headerfile', 'datafile'});
      fsample = data.fsample;
      [tmpcfg, artifact] = ft_artifact_zvalue(tmpcfg, data);
    else
      cfg = ft_checkconfig(cfg, 'dataset2files', 'yes');
      cfg = ft_checkconfig(cfg, 'required', {'headerfile', 'datafile'});
      hdr = ft_read_header(cfg.headerfile);
      fsample = hdr.Fs;
      tmpcfg.datafile    = cfg.datafile;
      tmpcfg.headerfile  = cfg.headerfile;
      [tmpcfg, artifact] = ft_artifact_zvalue(tmpcfg);
    end
    cfg.artfctdef.tms = tmpcfg.artfctdef.zvalue;

    % adjust artifact definition so that Nx2 matrix contains detected TMS
    % events with user-specified pre- and post stimulus period included.
    % The reason for this is that ft_artifact_zvalue centers the period
    % marked as artifactual around the detected event. In the case of a TMS
    % pulse the window you would like to mark as artifactual is not
    % symmetrical around the onset of the pulse.

    % get values and express in samples
    prestim = round(cfg.prestim * fsample);
    poststim = round(cfg.poststim * fsample);

    % adjust Nx2 artifact matrix to be centered non-symmetrically around
    % detected TMS-pulse
    artifact(:,1) = (artifact(:,1)+artifact(:,2))./2 - prestim;
    artifact(:,2) = artifact(:,1) + poststim;
    cfg.artfctdef.tms.artifact = artifact;

  case 'marker'
    % Check if the cfg is correct for this method
    cfg = ft_checkconfig(cfg, 'dataset2files', 'yes');
    ft_checkconfig(cfg, 'required', 'trialdef');
    cfg.trialfun = ft_getopt(cfg, 'trialfun', 'ft_trialfun_general');
    trialdef = cfg.trialdef;
    trialdef.prestim = cfg.prestim;
    trialdef.poststim = cfg.poststim;
    cfg.trialdef = ft_checkconfig(trialdef, 'required', {'eventvalue', 'eventtype'});

    % Get the trialfun
    cfg.trialfun = ft_getuserfun(cfg.trialfun, 'trialfun');

    % Evaluate the trialfun
    fprintf('evaluating trialfunction ''%s''\n', func2str(cfg.trialfun));
    trl   = feval(cfg.trialfun, cfg);

    % Prepare the found events for output
    artifact = trl(:,1:2);
    cfg.artfctdef.tms.artifact = artifact;
    fprintf('found %d events\n', size(artifact,1));

  otherwise
    % This should be redundant as ft_checkconfig does not allow other methods than the supported ones.
    ft_error('unsupported method');
end

% FIXME - not removing this causes problems when passing to ft_preprocessing
cfg = rmfield(cfg, 'method');


text:
function success = test_most(verbose, exit_on_fail)
%TEST_MOST  Run all MOST tests.
%   TEST_MOST
%   TEST_MOST(VERBOSE)
%   TEST_MOST(VERBOSE, EXIT_ON_FAIL)
%   SUCCESS = TEST_MOST(...)
%
%   Runs all of the MOST tests. If VERBOSE is true (false by default),
%   it prints the details of the individual tests. If EXIT_ON_FAIL is true
%   (false by default), it will exit MATLAB or Octave with a status of 1
%   unless T_RUN_TESTS returns ALL_OK.
%
%   See also T_RUN_TESTS.

%   MOST
%   Copyright (c) 2004-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MOST.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/most for more info.

if nargin < 2
    exit_on_fail = 0;
    if nargin < 1
        verbose = 0;
    end
end

tests = {};

%% MOST tests
have_c3sopf = exist('c3sopf', 'file') == 2;
tests{end+1} = 't_most_3b_1_1_0';
tests{end+1} = 't_most_3b_3_1_0';
if have_c3sopf
    tests{end+1} = 't_most_3b_1_1_2';
    tests{end+1} = 't_most_3b_3_1_2';
end
tests{end+1} = 't_most_30b_1_1_0';
tests{end+1} = 't_most_30b_3_1_0';
if have_c3sopf
    tests{end+1} = 't_most_30b_1_1_17';
    tests{end+1} = 't_most_30b_3_1_17';
end
tests{end+1} = 't_most_fixed_res';
tests{end+1} = 't_most_30b_1_1_0_uc';
if have_c3sopf
    tests{end+1} = 't_most_sp';
    tests{end+1} = 't_most_spuc';
end
tests{end+1} = 't_most_mpopf';
tests{end+1} = 't_most_uc';
tests{end+1} = 't_most_suc';
tests{end+1} = 't_most_tlmp';
tests{end+1} = 't_most_w_ds';

%% run the tests
all_ok = t_run_tests( tests, verbose );

%% handle success/failure
if exit_on_fail && ~all_ok
    exit(1);
end
if nargout
    success = all_ok;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function []=plot_dramapB(vResults)
% --------------------------------------------
% Plot Mc map on topography
%
% Incoming variables:
% vResults: Necessarily includes vResults.mPolygon, vResults.mValueGrid, vResults.fSpacingHorizontal
%
% Author: J. Woessner, woessner@seismo.ifg.ethz.ch
% last update: 27.02.03

% get tbase.bin values
[tmap,tmapleg] = tbase(3);

% Construct a map graticule mesh for surface object display
[lat,lon] = meshgrat(tmap,tmapleg);

% Get coordinates from calculated Mc values, set up matrices
vX=min(vResults.mPolygon(:,1)):vResults.fSpacingHorizontal:max(vResults.mPolygon(:,1));
vY=min(vResults.mPolygon(:,2)):vResults.fSpacingHorizontal:max(vResults.mPolygon(:,2));
[X , Y]  = meshgrid(vX,vY);

% % Set values of Mc to NaN
 vSel = find(vResults.mValueGrid(:,9)> 2.2);
 vResults.mValueGrid(vSel,9) = NaN;

% Reshape
mBvalues=reshape(vResults.mValueGrid(:,9),length(vY),length(vX));

% Interpolate Values
mBGrid = interp2(X,Y,mBvalues,lon,lat);

% Create figure
exfig=figure_exists('b-value worldmap',0);
if exfig==0
    ui=figure_w_normalized_uicontrolunits( 'Name','b-value worldmap',...
        'NumberTitle','off');
    figure_w_normalized_uicontrolunits(ui)
end

clf
hold on;
axis off;
axesm('MapProjection','robinson')

% Magic coloring
mi = min(min(mBGrid));
l =  isnan(mBGrid);
mBGrid(l) = mi-10;
ll = tmap < 0 & mBGrid < 0;
mBGrid(ll) = mBGrid(ll)*0 + 10;
hMap=meshm(mBGrid,tmapleg,size(tmap),tmap);

daspectm('m',20);
tightmap
view([0 90])
camlight; lighting phong
set(gca,'projection','perspective');

j = jet(60000);
j = j(30000:60000,:);
%j = winter(80000);
%j = flipud(j);
j = [ [ 0.9 0.9 0.9 ] ; j; [ 0.5 0.5 0.5] ];

%caxis([ min(min(mBvalues)) max(max(mBvalues)) ]);
%Harvard 0-70km
%caxis([0.6 2.2])
% ISC 0-70km
caxis([ 0.5 2]);
colormap(j); brighten(0.1);

axis off; set(gcf,'color','k')

setm(gca,'ffacecolor','k')
setm(gca,'fedgecolor','w','flinewidth',1);

setm(gca,'mlabellocation',60)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',25)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','w','FontSize',10,'Fontweight','bold','Fontname','times','Labelunits','degrees')
setm(gca,'Labelformat','none')
setm(gca,'Grid','on')
setm(gca,'GLineStyle','-')

h5 = colorbar;
% Harvard 0-70km
% set(h5,'position',[0.82 0.35 0.01 0.3],'TickDir','out','Ycolor','w','Xcolor','w',...
%     'Fontweight','bold','FontSize',10,'YTick',[0.6 1 1.4  1.8 2.2] );
% ISC 0-70km
set(h5,'position',[0.82 0.35 0.01 0.3],'TickDir','out','Ycolor','w','Xcolor','w',...
    'Fontweight','bold','FontSize',10,'Ytick',[0.5 1 1.5 2]);

set(gcf,'Inverthardcopy','off');

author: CelsoReyes


text:
% script for compiled EEGLAB version about button

tmptmptext = strvcat('                                  EEGLAB 1997-2020 (compiled version)', ...
                  '       ', ...
                  '       ', ...
                  'The EEGLAB electrophysiological data analysis software is an open source research', ...
                  'project developed at the Swartz Center for Computational Neuroscience, Institute', ...
                  'for Neural Computation at the University of California San Diego in collaboration', ...
                  'with The Salk Institute (La Jolla) and the CERCO laboratory (Toulouse) where part', ...
                  'of the core code was developed.', ...
                  ' ', ...
                  'EEGLAB incorporates several third party suites of functions, Fieldtrip for source ', ...
                  'localization, and BIOSIG for binary file import. A list of plugins is also included', ...
                  'in this compiled release. See http://sccn.ucsd.edu/eeglab/plugins.html for more', ...
                  'information.', ...
                  ' ', ...
                  'The EEGLAB software, its source code, and the source code of all additional third', ...
                  'party functions it includes is protected by the BSD open source license.', ...
                  ' ');

warndlg2(tmptmptext);
clear tmptmptext;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function serverStarted = mrmCheckServer(host)
%Test whether the server has been started.
%
%     serverStarted = mrmCheckServer(host)
%
% There should be a way to test without creating a new window.  But I am
% not sure how.
%
% (c) Stanford Vista Team, 2008

if ieNotDefined('host'), host = 'localhost'; end

% Try opening a window
stat = mrMesh(host, 999, 'refresh');

% If the server is not started, the window won't open and we get a return
% of -1000. Otherwise, the window was properly opened.
if stat ~= -1000
    % Window was opened properly by server.  So close it and return.
    serverStarted = 1;
    mrMesh(host,999,'close');
    return;
else
    % Tell 'em that the call didn't work, so the server isn't running.
    serverStarted = 0;
end

return;


text:
function varargout = process_test_parametric1( varargin )
% PROCESS_TEST_PARAMETRIC1: Parametric one-sample tests (test vs zero).

% @=============================================================================
% 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, Dimitrios Pantazis, 2015-2016

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Parametric test against zero';
    sProcess.Category    = 'Stat1';
    sProcess.SubGroup    = 'Test';
    sProcess.Index       = 701;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Statistics';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',  'results',  'timefreq',  'matrix'};
    sProcess.OutputTypes = {'pdata', 'presults', 'ptimefreq', 'pmatrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 2;

    % === GENERIC EXTRACT OPTIONS
    % Label
    sProcess.options.extract_title.Comment    = 'Select data to test:';
    sProcess.options.extract_title.Type       = 'label';
    sProcess.options.extract_title.InputTypes = {'data', 'results', 'timefreq', 'matrix'};
    % Options
    sProcess = process_extract_values('DefineExtractOptions', sProcess);
    % DISABLE ABSOLUTE VALUE
    sProcess.options.isabs.Value = 0;
    sProcess.options.isnorm.Value = 0;
    sProcess.options.isabs.Hidden = 1;
    sProcess.options.isnorm.Hidden = 1;

    % === OUTPUT COMMENT
    sProcess.options.Comment.Comment = 'Comment (empty=default): ';
    sProcess.options.Comment.Type    = 'text';
    sProcess.options.Comment.Value   = '';

    % === TEST: title
    sProcess.options.test_title.Comment    = '
Test statistic:';
    sProcess.options.test_title.Type       = 'label';
    % === TEST: type
    sProcess.options.test_type.Comment = {['One-sample Student''s t-test    X~N(m,s)
' ...
                                           't = mean(X) ./ std(X) .* sqrt(n)      df=n-1'], ...
                                          ['One-sample Chi2 test     Zi~N(0,1), i=1..n
' ...
                                           'Q = sum(|Zi|^2)     Q~Chi2(n)'], ...
                                          ['One-sample Chi2 test (unconstrained sources)     Zix,Ziy,Ziz~N(0,1), i=1..n
' ...
                                           'Q = sum(|Zix|^2 + |Ziy|^2 + |Ziz|^2)     Q~Chi2(3*n)']; ...
                                          'ttest_onesample', 'chi2_onesample', 'chi2_onesample_unconstr'};
    sProcess.options.test_type.Type    = 'radio_label';
    sProcess.options.test_type.Value   = 'ttest_onesample';
    % === TAIL FOR THE TEST STATISTIC
    sProcess.options.tail.Comment  = {'One-tailed (-)', 'Two-tailed', 'One-tailed (+)', ''; ...
                                      'one-', 'two', 'one+', ''};
    sProcess.options.tail.Type     = 'radio_linelabel';
    sProcess.options.tail.Value    = 'two';
end

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

%% ===== RUN =====
function sOutput = Run(sProcess, sInputs) %#ok<DEFNU>
    sOutput = process_test_parametric2('Run', sProcess, sInputs, []);
end

author: brainstorm-tools


text:
function dmapErrors = evalDepthMaps(class,jobID,metric)
    globals;
    if(nargin<3)
        metric = 'zmae';
    end
    dmapDir = jobDirs(class,jobID,'dmap');
    evalFileName = jobDirs(class,jobID,'evalDepth');
    statesDir = jobDirs(class, jobID, 'state');
    gtDmapDir = fullfile(cachedir,class,'gtDepthMap');
    if(~exist(dmapDir,'dir'))
        error('Depth maps not found. Generate them first!');
    end

    if(~exist(gtDmapDir,'dir'))
        error('Ground truth depth maps not found. Generate them first!');
    end

    switch metric
        case 'corr'
            disp('Using Depth Correlation');
        case 'zmae'
            disp('Using Z-MAE');
        case 'absrel'
            disp('Using absolute relative error');
    end

    fnames = getFileNamesFromDirectory(dmapDir,'types',{'.mat'});
    fnames = removeFlipNames(fnames);

    disp('Depth map evaluation');

    dmapErrors = zeros(length(fnames),1);
    p =  TimedProgressBar( length(fnames), 30, ...
    'Depth Map Errors: Remaining ', ', Completed ', 'Depth Map evaluation Time: ' );

    parfor i=1:length(fnames)
        GtFile = fullfile(gtDmapDir,fnames{i});
        dmapFile = fullfile(dmapDir,fnames{i});
        state = load(fullfile(statesDir,fnames{i}),'state');
        state = state.state;
        if(~exist(GtFile,'file'))
            warning('Gt depth map not found: %s\n',GtFile);
            dmapErrors(i) = nan;
            continue;
        end
        gtdmap = load(GtFile,'dmap');
        gtdmap = gtdmap.dmap;
        dmap = load(dmapFile,'dmap');dmap = dmap.dmap;
        switch metric
            case 'corr'
                dmapErrors(i) = dmapMetricCorr(dmap,gtdmap,state.gtmask);
            case 'zmae'
                dmapErrors(i) = dmapMetricZMAE(dmap,gtdmap,state.gtmask);
            case 'absrel'
                dmapErrors(i) = dmapMetricRel(dmap,gtdmap,state.gtmask);
        end

        % Visualization
        if 0
            gtdmap(isinf(gtdmap)) = nan;
            dmap(isinf(dmap)) = nan;
            gtdmap(~state.mask) = nan;
            dmap(~state.mask) = nan;
            subplot(1,2,1);
            out = visualizeDEM(dmap);
            imagesc(imcrop(out,state.bbox));axis equal off;
            title(sprintf('Opt Depth Map - %.4f',dmapErrors(i)));
            subplot(1,2,2)
            out = visualizeDEM(gtdmap);
            imagesc(imcrop(out,state.bbox)); axis equal off;
            title('Gt Depth Map');
            pause;
            clf;
        end
        p.progress;
    end
    p.stop;
    save(evalFileName,'dmapErrors','fnames');
end

repo: CategoryShapes


text:
classdef generalization_DS < handle

%  This datasource object (DS) allows one to train a classifier on a specific
%  set of labels, and then test the classifier on a different set of
%  labels - which enables one to evaluate how similar neural representations
%  are across different but related conditions (i.e., does training on one set of
%  conditions generalization to a different but related set of conditions?). This datasource
%  is a subclass of the handle class (i.e., it has a persistent state) and
%  contains a basic_DS where it gets most of its functionality.
%
%  The constructor for this datasource has the same arguments
%  as basic_DS, plus two additional arguments 'the_training_label_names'
%  and 'the_test_label_names' i.e., the constructor has the form:
%
%     ds = generalization_DS(binned_data_name, specific_binned_label_name, num_cv_splits, the_training_label_names, the_test_label_names, load_data_as_spike_counts)
%
%      the_training_label_names and the_test_label_names are cell arrays that
%        specify which labels should belong to which class, with the first element
%        of these cells arrays specifying the training/test labels for the first class
%        the second element of the cell array specifies which labels belong to
%        the second class, etc..  For example, suppose one was interested in testing
%        position invariance, and had done an experiment in which data was recorded
%        while 7 different objects were shown at three different locations.  If the
%        labels for the 7 objects at the first location had labels 'obj1_loc1', 'obj2_loc1', ..., 'obj7_loc1',
%        at the second location were 'obj1_loc2', 'obj2_loc2', ..., 'obj7_loc2',
%        and at the third location were 'obj1_loc3', 'obj2_loc3', ..., 'obj7_loc3',
%        then one could do a test of position invariance by setting the_training_label_names{1} = {'obj1_loc1},
%        setting the_training_label_names{2} = {'obj2_loc1'}, ...,  the_training_label_names{7} = {'obj7_loc1},
%        and setting the the_test_label_names{1} = {'obj1_loc2', 'obj1_loc3'},
%        the_test_label_names{2} = {'obj2_loc2', 'obj2_loc3'}, ..., the_test_label_names{7} = {'obj7_loc2', 'obj7_loc3'}.
%        The object is able to test such generalization from training on one set of labels and
%        testing on a different set of labels by remapping the training label numbers to the
%        index number in the_training_label_names cell array, and remapping the
%        test label numbers with the the index number into the the_test_label_names cell array.
%
%  This object has all of the same properites as the basic_DS object (except that label_names_to_use which
%    has been replaced by the the_training_label_names, and the_test_label_names properties).  There
%    is also an additional property that can be set for this object which is:
%
%       use_unique_data_in_each_CV_split  (default value is 0).
%
%       When this argument is set to 0, the get_data method returns the normal leave one split
%        out training and test data sets (i.e., the training set consists of
%        (num_cv_splits - 1) splits of the data and the test set consists of 1 split of the data).
%        The data in the training still comes from different splits as the data in the
%        test set, thus one can have some of the same labels in the both
%        the_training_label_names and in the_test_label_names (in fact, if ones sets
%        the_test_label_names = the_training_label_names, then the get_data
%        method will be the same as the basic_DS get_data method).  However,
%        if use_unique_data_in_each_CV_split = 1, then each training
%        and test set will consist data from only split, and thus each cross-validation
%        run is essentially like running an independent decoding experiment.
%        In this case the_training_label_names and the_test_label_names must not contain
%        any of the same labels (otherwise, they would be copies of the same data
%        which would violate the fact that the training and the test set must not have any of the same data).

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

%     This code is part of the Neural Decoding Toolbox.
%     Copyright (C) 2011 by Ethan Meyers (emeyers@mit.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 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/>.

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

properties

    the_basic_DS = [];      % the basic_DS that will be used to give this object most of its functionality

    the_training_label_names = [];   % a cell array specifying which label names (or numbers) should be used for training the classifier
                                       %   i.e., the_training_label_names{1} = {class_1_training_names}; the_training_label_names{2} = {class_2_training_names}; etc.
    the_test_label_names = [];       % a cell array specifying which label names (or numbers) should be used for testing the classifier
                                       %   i.e., the_test_label_names{1} = {class_1_test_names}; the_test_label_names{2} = {class_2_test_names}; etc.

    use_unique_data_in_each_CV_split = 0;  % if this is set to 1 then each CV splits has unique data, i.e.,
                                                         %  each CV split has the amount of data that is typically only in the test set,
                                                         %  and every CV training set does not consist of (num_cv_splits -1) * num_labels data points
                                                         %  but instead consists of length(cell2mat(the_test_label_names)) training points.

    % some properties of basic_DS that will also be available in generalization_DS by setting the basic_DS properties

    num_times_to_repeat_each_label_per_cv_split = 1;     %  how many of each unique label should be in each CV block

    sample_sites_with_replacement = 0;  %  specify whether to sample neurons with replacement - if they are sampled with replacement, then some features will be repeated within a single point data vector
    num_resample_sites = -1;            %  how many sites should be used for each resample iteration - must be less than length(the_data)

    create_simultaneously_recorded_populations = 0;   % to use pseudo-populations or simultaneous populations (2 => that the training set is pseudo and test is simultaneous)

    sites_to_use = -1;                   %  a list of indices of which sites (neurons) to use in the the_data cell array
    sites_to_exclude = [];               %  a list of features that should explicitly be excluded
    time_periods_to_get_data_from = [];  %  a cell array containing vectors that specify which time bins to use from the_data

    % randomly shuffles the labels prior to the get_data method being called - which is useful for creating one point in a null distribution to check if decoding results are above what is expected by change.
    randomly_shuffle_labels_before_running = 0;

end

properties (GetAccess = 'public', SetAccess = 'private')

      initialized = 0;

      the_training_label_numbers = [];   % numbers that the_training_label_names were mapping on to
      the_test_label_numbers = [];       % numbers that the_test_label_names were mapping on to
end

methods

    function ds = generalization_DS(binned_data_name, specific_binned_label_name, num_cv_splits, the_training_label_names, the_test_label_names, load_data_as_spike_counts)

       if nargin < 6
           load_data_as_spike_counts = 0;
       end

       ds.the_basic_DS = basic_DS(binned_data_name, specific_binned_label_name, num_cv_splits, load_data_as_spike_counts);   % set properties using parent constructor

       ds.the_training_label_names = the_training_label_names;   % set properties that are unique to this object
       ds.the_test_label_names = the_test_label_names;

       % get all the unique labels
       if length(the_training_label_names) ~= length(the_test_label_names)
            error('The cell arrays the_training_label_names and the_test_label_names must have the same number of cells (with the names/numbers in cell i containing which labels belong to class i)');
       end

    end

    function the_properties = get_DS_properties(ds)

        the_properties = ds.the_basic_DS.get_DS_properties;

        the_properties.the_training_label_names = ds.the_training_label_names;
        the_properties.the_test_label_names  = ds.the_test_label_names;

        the_properties.use_unique_data_in_each_CV_split = ds.use_unique_data_in_each_CV_split;

    end

    function  [XTr_all_time_cv YTr_all XTe_all_time_cv YTe_all] = get_data(ds)
    %function  [XTr_all_time_cv YTr_all XTe_all_time_cv YTe_all ADDITIONAL_DATASOURCE_INFO] = get_data(ds)

        % initialize things here....
        if ds.initialized == 0

            % creating separate variables for this since calling a field of an object in Matlab is slow
            the_training_label_names = ds.the_training_label_names;
            the_test_label_names = ds.the_test_label_names;

            ds.the_basic_DS.num_times_to_repeat_each_label_per_cv_split = ds.num_times_to_repeat_each_label_per_cv_split;

            ds.the_basic_DS.sample_sites_with_replacement = ds.sample_sites_with_replacement ;
            ds.the_basic_DS.num_resample_sites = ds.num_resample_sites;

            ds.the_basic_DS.create_simultaneously_recorded_populations = ds.create_simultaneously_recorded_populations;

            ds.the_basic_DS.sites_to_use = ds.sites_to_use;
            ds.the_basic_DS.sites_to_exclude = ds.sites_to_exclude;
            ds.the_basic_DS.time_periods_to_get_data_from = ds.time_periods_to_get_data_from;

            ds.the_basic_DS.randomly_shuffle_labels_before_running = ds.randomly_shuffle_labels_before_running;

            % if the_labels are strings, convert these names into label numbers...
            if iscell(ds.the_basic_DS.the_labels{1})

                % put all the training and test names into one cell array called all_training_and_test_names
                cTrainingNames = 0;
                cTestNames = 0;
                for iClass = 1:length(the_training_label_names)
                    for iTrain = 1:size(the_training_label_names{iClass}, 2)
                        cTrainingNames = cTrainingNames + 1;
                        all_training_names{cTrainingNames} = the_training_label_names{iClass}{iTrain};
                    end

                    for iTest = 1:size(the_test_label_names{iClass}, 2)
                        cTestNames = cTestNames + 1;
                        all_test_names{cTestNames} = the_test_label_names{iClass}{iTest};
                    end
                end

               all_training_and_test_names = union(all_training_names, all_test_names);

               % sanity checka that the same label is not in multiple training classes (and the the same label is not in multiple test classes)
               if length(unique(all_training_names)) ~= length(all_training_names)
                    warning('some of the same strings in the_training_label_names are in multiple classes');
               end
               if length(unique(all_test_names)) ~= length(all_test_names)  % same sanity check for the_test_label_names
                    warning('some of the same strings in the_test_label_names are in multiple classes');
               end

               % convert the_labels into numebrs
               ignore_case_of_strings = 0;  % for now, always respect the case of the strings used in the labels

                % by passing all_training_and_test_names as the 3rd argument convert_label_strings_into_numbers this keeps the mapping from the_test_label_names to the_test_label_numbers correct
                % an erorr in convert_label_strings_into_numbers will be thrown if some of the training or test label names are strings that are not in binned labels the_labels
                the_labels_as_numbers = convert_label_strings_into_numbers(ds.the_basic_DS.the_labels, ignore_case_of_strings, all_training_and_test_names);

                ds.the_basic_DS.the_labels = the_labels_as_numbers;

                % remap the_training_label_names and the_training_label_names into numbers
                cTrainingNames = 0;
                cTestNames = 0;
                for iClass = 1:length(the_training_label_names)
                    for iTrain = 1:size(the_training_label_names{iClass}, 2)
                        cTrainingNames = cTrainingNames + 1;
                        the_training_label_numbers{iClass}(iTrain) = find(ismember(all_training_and_test_names, the_training_label_names{iClass}{iTrain}));
                    end

                    for iTest = 1:size(the_test_label_names{iClass}, 2)
                        cTestNames = cTestNames + 1;
                        the_test_label_numbers{iClass}(iTest) = find(ismember(all_training_and_test_names, the_test_label_names{iClass}{iTest}));
                    end
                end

               %  sanity check to make sure that one is not training with label l in class k and then has test label l in class j
               some_of_the_same_labels_are_in_multiple_classes = 0;
               for iClass = 1:length(the_training_label_names)
                    temp_test_label_numbers = the_test_label_numbers;
                    temp_test_label_numbers{iClass} = [];
                    temp_test_label_numbers = cell2mat(temp_test_label_numbers);
                    if ~isempty(intersect(temp_test_label_numbers, the_test_label_numbers{iClass}))
                        some_of_the_same_labels_are_in_multiple_classes = 1;
                    end
               end

               if  some_of_the_same_labels_are_in_multiple_classes == 1
                    warning('some labels that are in training class k, are in a different test class j')
               end

               ds.the_training_label_numbers = the_training_label_numbers;
               ds.the_test_label_numbers = the_test_label_numbers;

            else  % if numbers for labels have been specified instead of names, just use those numbers

               ds.the_training_label_numbers = ds.the_training_label_names;
               ds.the_test_label_numbers = ds.the_test_label_names;

            end

            % only use labels that are listed in the_training_label_numbers and the_test_label_numbers
            the_training_nums = cell2mat(ds.the_training_label_numbers);
            the_test_nums = cell2mat(ds.the_test_label_numbers);
            label_numbers_to_use = unique([the_training_nums(:); the_test_nums(:)]);
            ds.the_basic_DS.label_names_to_use = label_numbers_to_use;

            ds.initialized = 1;
        end

         % creating separate variables for this since calling a field of an object in Matlab is slow
         the_training_label_numbers = ds.the_training_label_numbers;
         the_test_label_numbers = ds.the_test_label_numbers;

        %[XTr_all_time_cv YTr_all XTe_all_time_cv YTe_all] = get_data@basic_DS(ds);  % old version where this object was a subclass of basic_DS
         [XTr_all_time_cv YTr_all XTe_all_time_cv YTe_all] = ds.the_basic_DS.get_data;  % this might be much slower than inheriting from basic_DS :(

        if ds.use_unique_data_in_each_CV_split == 1

            % if running each training and test set separately, there can not be any overlap between the labels listed in the training and test sets
            %(otherwise there could be some of the same data in the training and test sets which is completely forbidden!!!!!
            all_unique_training_labels = unique(cell2mat(the_training_label_numbers));
            all_unique_test_labels = unique(cell2mat(the_test_label_numbers));

            if length(intersect(all_unique_training_labels,  all_unique_test_labels)) > 0
                error('if running each split separately (i.e., use_unique_data_in_each_CV_split == 1), then none of the same labels can be in the training and test sets, otherwise there will be some of the same data in the training and test sets')
            end

            % remap labels (only using old 'test' data because I want each CV split to be independent)
            remapped_YTr_all = NaN .* ones(size(YTr_all));  % only getting labels (and data) from test set b/c want each CV split to contain unique data
            remapped_YTe_all = NaN .* ones(size(YTe_al));

            for iGroup = 1:length(the_training_label_numbers)
                remapped_YTr_all(ismember(YTr_all_cv, the_training_label_numbers{iGroup})) = iGroup;
                remapped_YTe_all(ismember(YTe_all, the_test_label_numbers{iGroup})) = iGroup;
            end

           % remove data from trials in which the labels are not appropriate for the training/test sets
           train_inds  = ~isnan(remapped_YTr_all);
           test_inds  = ~isnan(remapped_YTe_all);

           YTr_all = remapped_YTr_all(train_inds);   % remove NaNs from remapped labels
           YTe_all = remapped_YTe_all(test_inds);

            for iCV = 1:length(XTr_all_time_cv{1})
               for iTime = 1:length(XTr_all_time_cv)
                   XTr_all_time_cv{iTime}{iCV} = XTr_all_time_cv{iTime}{iCV}(:, train_inds);   % only getting data from test set b/c want each CV split to be independent
                   XTe_all_time_cv{iTime}{iCV} = XTe_all_time_cv{iTime}{iCV}(:, test_inds);
               end
            end

        else

            % remap labels
            remapped_YTr_all = NaN .* ones(size(YTr_all));
            remapped_YTe_all = NaN .* ones(size(YTe_all));

            for iGroup = 1:length(the_training_label_numbers)
                remapped_YTr_all(ismember(YTr_all, the_training_label_numbers{iGroup})) = iGroup;
                remapped_YTe_all(ismember(YTe_all, the_test_label_numbers{iGroup})) = iGroup;
            end

           % remove data from trials in which the labels are not appropriate for the training/test sets
           training_data_inds_to_remove = isnan(remapped_YTr_all);
           test_data_inds_to_remove = isnan(remapped_YTe_all);

           for iTime = 1:length(XTr_all_time_cv)
               for iCV = 1:length(XTr_all_time_cv{1})
                    XTr_all_time_cv{iTime}{iCV}(:, training_data_inds_to_remove) = [];
                    XTe_all_time_cv{iTime}{iCV}(:, test_data_inds_to_remove) = [];
               end
           end

           remapped_YTr_all(training_data_inds_to_remove) = [];
           remapped_YTe_all(test_data_inds_to_remove) = [];

           YTr_all = remapped_YTr_all;
           YTe_all = remapped_YTe_all;

        end

    end  % end get_data method

end   % end methods

end % end class

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function showsolution3(node,elem,u,expr,varargin)
%% SHOWSOLUTION3 plots the solution u on a tetrahedron mesh in 3-D.
%
%    showsolution3(node,elem,u) displays the functoin u on a topological
%    2-dimensional mesh given by node and elem matrices. The function u
%    could be piecewise constant or piecewise linear.
%
%    showsolution3(node,elem,u,expr) displays the function u on the
%    boundary of parts of the mesh specificed by the expression. For
%    example, showsoluiton3(node,elem,'z==0') will show the function on the
%    cross section of z=0.
%
%    showsolution3(node,elem,u,viewangle) changes the display angle. The
%    deault view angle on planar meshes is view(2) and view(3) for surface
%    meshes.
%
%    showsolution3(node,elem,u,expr,'param','value','param','value'...) allows
%    additional patch param/value pairs to be used when displaying the
%    mesh.
%
%   Example:
%     f = inline('x.^2 + y.^2 + z.^2');
%     node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
%     elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7];
%     for k=1:4
%       [node,elem] = uniformbisect3(node,elem);
%     end
%     u = f(node(:,1),node(:,2),node(:,3));
%     subplot(1,2,1);
%     showsolution3(node,elem,u);
%     subplot(1,2,2);
%     showsolution3(node,elem,u,'~(x>0 & y>0)',[139,16],'EdgeColor','k');
%
%   See also showmesh, showsolution3.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if (nargin >= 4) && (any(expr))
    x = node(:,1);  y = node(:,2);  z = node(:,3); %#ok<*NASGU>
    incl = find(eval(expr));
    elem = elem(any(ismember(elem,incl),2),:);
end
[bdNode, bdFace] = findboundary3(elem); %#ok<ASGLU>
if isempty(varargin)
    showsolution(node,bdFace,u,'EdgeColor','k');
elseif (nargin == 5) && isnumeric(varargin{1})
    showsolution(node,bdFace,u,'EdgeColor','k');
    view(varargin{1});
else
    showsolution(node,bdFace,u,varargin{1:end});
end

repo: ifem


text:
function C = train_HLDA(XTr, YTr, nSegments, varargin)
% TRAIN_HLDA - Hierarchical linear discriminant analysis
%
%Synopsis:
%   C = train_HLDA(XTr, YTr, nSegments, varargin)
%
%Arguments:
%   XTR: DOUBLE [TxNxM] - Data matrix, with T temporal features, N
%                           N spatial features, and M
%                           training points/examples.
%   YTR: INT [CxM]      - Class membership labels of points in X_TR. C by M
%                           matrix of training labels, with C representing
%                           the number of classes and M the number of
%                           training examples/points.
%                           YTR(i,j)==1 if the point j belongs to class i.
%   nSegments: INT      - the number of non-overlapping segments the interval
%                           should be separated in
%   OPT: PROPLIST       - Structure or property/value list of optional
%                           properties. Options are also passed to clsutil_shrinkage.
%     'Regression'  - BOOL (default 0): If true, the top level classifier is
%                           a logistic regression classifier.
%     'nChannels'   - INT (default 0): Used when using 'crossvalidation' in order to
%                           reconstruct the 2D feature matrix to its original 3D
%                           [T x Ch x N] shape.
%Returns:
%   C: STRUCT           - Structure containing the trained classifiers for
%                           each segment and the final top level
%                           classifier. The structure C includes the fields:
%       'seg': STRUCT []      - contains a trained LDA classifier for each segment
%       'final': STRUCT       - the final top-level LDA classifier
%       'nChannels': INT      - (optional) the number of channels in the training data
%Description:
%   train_HLDA trains a hierarchical LDA classifier given training data,
%   labels and a number of segments. Either LDA (default) or logistic regression
%   is used as a top-level classifier.
%
%   References:Gerson, A.D., Parra, L.C., Sajda, P.: Cortically coupled
%   computer vision for rapid image search. IEEE Transactions on Neural
%   Systems and Rehabilitation Engineering 14, 174–179 (2006).

%
%
%Examples:
%   train_HLDA(XTr, YTr, nSegments)
%   train_HLDA(XTr, YTr, nSegments, 'Regression', 1, 'nChannels', 64)
%
%See also:
%   APPLY_HLDA

props= {'Regression'      0                             'BOOL'
        'nChannels'       0                             'INT'
       };

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

% validate argument types
misc_checkType(XTr, 'DOUBLE');

if opt.nChannels ~= 0
  misc_checkType(XTr, 'DOUBLE[- -]');
end

misc_checkType(YTr, 'DOUBLE[2 -]');
misc_checkType(nSegments, 'INT');

dims = size(XTr);

% make data matrix 3D if it has been tranformed to 2D for crossvalidation
if (length(size(XTr)) == 2) && opt.nChannels > 0
  XTr = reshape(XTr, [], opt.nChannels, dims(2));
end

dims = size(XTr);

%boundary indices between segments
seg_idx = round(linspace(0, dims(1), nSegments+1));

seg_scores = zeros(nSegments, dims(end));

for i = 1:nSegments
    if length(size(XTr)) > 2
      seg = XTr(seg_idx(i) + 1 : seg_idx(i + 1), :, :); % i-th segment of X
      seg = reshape(seg, [], dims(end));
    else
      seg = XTr(seg_idx(i) + 1 : seg_idx(i + 1), :); % i-th segment of X
    end
    %train LDA classifier for this segment
    seg_LDA = train_RLDAshrink(seg, YTr, 'Scaling', 1);
    C.seg(i) = seg_LDA;
    seg_scores(i,:) = apply_separatingHyperplane(seg_LDA, seg); %get classifier scores for segment
end

if opt.Regression
    %logistic regression as top-level classifier
    C.final.B = mnrfit(seg_scores', YTr');
else
    %LDA as top-level classifier
    C.final = train_RLDAshrink(seg_scores, YTr, 'Scaling', 1);
end

% add number of channels to classifier so they can be used in the apply function
if opt.nChannels > 0
  C.nChannels = opt.nChannels;
end

author: bbci


text:
% DATA = HADSST2_PARSEDATA()

% Last modified 2011-01-11
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function data = hadsst3_load_data()

filename = '/share/climate/data/HadSST3/HadSST3_mostlikely.nc';
%filename = '/share/climate/data/UK_Met_Office/HadSST2/HadSST2_nobs.nc';

longitude = nc_varget(filename, 'longitude');
latitude = nc_varget(filename, 'latitude');
[LON,LAT] = meshgrid(longitude, latitude);
coordinates = [LON(:)';LAT(:)'];

time = nc_varget(filename, 'time');
time = time(:)' + datenum([1850 0 0]);

observations = nc_varget(filename, 'sst');
%observations(observations==raw.data.FillValue) = NaN;
observations = reshape(observations, [length(time), size(coordinates, 2)])';

data = struct('observations', double(observations), ...
              'coordinates', double(coordinates), ...
              'longitude', double(longitude), ...
              'latitude', double(latitude), ...
              'time', double(time));

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab