text:
% DEMTOYUNSUPERVISED A script to run unsupervised deep GP on toy hierarchical data.
%
% DESC A script to run unsupervised deep GP on toy hierarchical data. The script provides
% the option of parametrising the model and initialisation in many many
% different ways... but the core of the demo (define a deep GP and train
% it) is actually not that big, if you decide to use the default options.
%
% COPYRIGHT: Andreas C. Damianou, 2013
%
% SEE ALSO: demToyRegression.m
%
% DEEPGP

experimentNo = 1;
toyType = 'hgplvmSampleTr1';
baseKern='rbfardjit'; % The mapping kernel between the layers
Q = {6,4}; % Dimensionality of the latent space in each layer
initSNR = {100, 50}; % Initial Signal to Noise ration per layer
% How to initialise X when multiple output modalities are present. See
% hsvargplvm_init for details
initial_X = 'separately';
% How to initialise X (more specifically, means of q(X)) for each layer.
% The selected method (potentially
% different per layer) will be applied in a sequential fashion, eg if we
% use PCA we obtain X_1 from PCA on Y, then X_2 from PCA on X_1 etc. Deep
% GPs are completely different than a stacked method, since X's will be
% integrated out (the initial X's above are actually the initial means of the var.
% distribution) and everythin will be optimised jointly.
% Here we opt for a Bayesian GPLVM that  gives the initial X.
% See hsvargplvm_init for other options (eg pca).
initX = 'vargplvm';
%- options for the BayesianGPLVM used to initialise the variational means
stackedInitIters = 200;
stackedInitVardistIters = 100;
stackedInitSNR = 100;
initVardistIters = 100;
demToyHsvargplvm1; % Run the actual demo

%% --- Plot true data
subplot(3,2,1)
myPlot(Z{3},'X2',[],[],{3,8},0)
subplot(3,2,3)
myPlot(Z{1},'XA',[],[],{3,8},0)
subplot(3,2,4)
myPlot(Z{2},'XB',[],[],{3,8},0)
subplot(3,2,5)
plot(Ytr{1},'x-'); title('YA');
subplot(3,2,6)
plot(Ytr{2},'x-'); title('YB');

%% -- Plot spaces discovered by deep GPs (two most dominant dimensions for
%% top layer and similarly for each of the two modalities of layer 1)
figure
hsvargplvmShowScales(model);

s2 = sort(vargplvmRetainedScales(model.layer{2}.comp{1}));
sA =  sort(vargplvmRetainedScales(model.layer{1}.comp{1}));
sB =  sort(vargplvmRetainedScales(model.layer{1}.comp{2}));

figure
subplot(2,2,1)
myPlot(model.layer{2}.vardist.means(:,s2(1:2)),'deepGP_X2',[],[],{3,8},0)
subplot(2,2,3)
myPlot(model.layer{1}.vardist.means(:,sA(1:2)),'deepGP_XA',[],[],{3,8},0)
subplot(2,2,4)
myPlot(model.layer{1}.vardist.means(:,sB(1:2)),'deepGP_XB',[],[],{3,8},0)

%% --- Compare with stacked Bayesian GP-LVM % TODO
%[XA, s, WA, modelA] = vargplvmEmbed(Ytr{1}, 5, initXOptions{1}{:});
%[XB, s, WB, modelB] = vargplvmEmbed(Ytr{2}, 5, initXOptions{1}{:});
%[X2, s, W2, model2]  = vargplvmEmbed([XA XB], 5, initXOptions{2}{:});

%% --- Compare with stacked PCA and isomap

figure
pcaXA = ppcaEmbed(Ytr{1}, 2);
pcaXB = ppcaEmbed(Ytr{2},2);
pcaX2 = ppcaEmbed([pcaXA pcaXB],2);
subplot(2,2,1)
myPlot(pcaX2,'pcaX2',[],[],{3,8},0)
subplot(2,2,3)
myPlot(pcaXA,'pcaXA',[],[],{3,8},0)
subplot(2,2,4)
myPlot(pcaXB,'pcaXB',[],[],{3,8},0)

figure
isomapXA = isomap2Embed(Ytr{1}, 2);
isomapXB = isomap2Embed(Ytr{2},2);
isomapX2 = isomap2Embed([isomapXA isomapXB],2);
subplot(2,2,1)
myPlot(isomapX2,'isomapX2',[],[],{3,8},0)
subplot(2,2,3)
myPlot(isomapXA,'isomapXA',[],[],{3,8},0)
subplot(2,2,4)
myPlot(isomapXB,'isomapXB',[],[],{3,8},0)

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
classdef svm < dml.method
% SVM support vector machine.
%
% DESCRIPTION
% Linear support vector machine classifier
%
%   EXAMPLE
%   X = rand(10,20); Y = [1 1 1 1 1 2 2 2 2 2]';
%   m = dml.svm
%   m = m.train(X,Y);
%   Z = m.test(X);
%
%   DEVELOPER
%   Jason Farquhar (j.farquhar@donders.ru.nl)

  properties

    X % training data

    primal % weights in primal form

    dual % weights in dual form

    C % regularization parameter

    kernel = 'linear'; % type of kernel

    Ktrain % precomputed kernel for training data
    Ktest  % precomputed kernel for test data

    distance = false; % return distance from decision boundary instead of class label
    native   = false; % uses native Bioinformatics toolbox SVM implementation if true
    issqrtmK = false;
  end

  methods

    function obj = svm(varargin)

      obj = obj@dml.method(varargin{:});

    end

    function obj = train(obj,X,Y)

      opts = {'verb' -1};

      % handle multiple datasets
      if iscell(X)
        obj = dml.ndata('method',obj);
        obj = obj.train(X,Y);
        return;
      end

      if obj.native
        obj.Ktrain = fitcsvm(X,Y);
        return
      end

      if obj.restart || isempty(obj.dual)

        obj.X      = X;
        obj.Ktrain = compKernel(X,X,obj.kernel,obj.Ktrain);
        obj.Ktest  = [];

        if isempty(obj.C)
          if obj.issqrtmK
            % the actual kernel will be K*K'
            diagK = sum(obj.Ktrain.^2,2);
            meanK = mean(reshape(obj.Ktrain*obj.Ktrain',[],1));
          else
            diagK = diag(obj.Ktrain);
            meanK = mean(obj.Ktrain(:));
          end
          obj.C = .1*(mean(diagK)-meanK);
          if obj.verbose, fprintf('using default C=%.2f\n',obj.C); end
        end

        if obj.issqrtmK
          opts = cat(2, opts, {'issqrtmK' true});
        end
        obj.dual = l2svm_cg(obj.Ktrain,2*(Y-1)-1,obj.C, opts);

      else

        % facilitate warm restarts
        obj.dual = l2svm_cg(obj.Ktrain,2*(Y-1)-1,obj.C, opts); %,'alphab',obj.dual);

      end

      obj.primal = 0;
      for j = 1:size(X,1)
        obj.primal = obj.primal + obj.dual(j)*X(j,:);
      end

    end

    function Y = test(obj,X)

      if obj.native
        C = svmclassify(obj.Ktrain,X);
        Y = zeros(size(C,1),2);
        Y(C + (0:2:(2*(size(C,1)-1)))')=1;
        return
      end

      if isempty(obj.Ktest)
        obj.Ktest = compKernel(X,obj.X,'linear');
      end

      probs = obj.Ktest * obj.dual(1:end-1) + obj.dual(end);

      Y = zeros(numel(probs),2);
      if obj.distance

        Y(probs<0,1) = abs(probs(probs<0));
        Y(probs>0,2) = abs(probs(probs>0));

      else
        % post is just the sign and does not have a probabilistic interpretation
        Y(:,1) = (probs < 0);
        Y(:,2) = (probs > 0);
      end

    end

    function m = model(obj)
    % returns
    %
    % m.primal primal form SVM parameters

      m.primal = obj.primal;

    end

  end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function x = mg_ns_iter(mgdata,x0,f,level,npre,npost,sweeps)
%mg_ns_iter    performs one MG iteration (Navier-Stokes)
%   x = mg_ns_iter(mgdata,x0,f,level,npre,npost,sweeps)
%   input
%          mgdata       structure containing matrices, grid transfer
%                       operators and smoothing operators
%          x0           initial iterate
%          f            right-hand side
%          level        grid level
%          npre         number of presmoothing steps
%          npost        number of postsmoothing steps
%          sweeps       type of sweeping strategy used for Gauss-Seidel
%                       smoothing
%   output
%          x            result of one MG iteration
%
%   IFISS function: HCE; 18 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

A = mgdata(level).matrix;
P = mgdata(level).prolong;
for i=2:level, smoother(i) = mgdata(i).smoother; end

if level==2,
   n = length(A);
   [U,Sig,V] = svd(full(A));
   s = diag(Sig);
   if abs(s(n))<1.d-14, si = [1./s(1:n-1);0];
   else                 si = 1./s(1:n);
   end
   Sigi = diag(si);
   x = V*(Sigi*(U'*f));
else

   % presmooth
   x = mg_pre(A,x0, f,npre,smoother,level,sweeps);
   % Restrict residual
   r = f - A*x;
   rc = P'*r;
   % coarse grid correction
   cc = mg_ns_iter(mgdata,zeros(size(rc)),rc,level-1,npre,npost,sweeps);
   % add this line for W-cycle
%  cc = mg_ns_iter(mgdata,cc,rc,level-1,npre,npost,sweeps);
   x = x + P*cc;
   % postsmooth
   x = mg_post(A,x,f,npost,smoother,level,sweeps);
end

author: johannesgerer

repo: jburkardt-m


text:
function [nWritten] = write_plx(filename, ch, freq, npw, n, ts, wave, units)
% write_plx: Write waveform data to a .plx file
%
% [nWritten] = write_plx(filename, ch, freq, npw, n, ts, wave, units)
%
% INPUT:
%   filename - if empty string, will use File Save dialog
%   ch - 1-based channel number that will be written to the .plx file
%   freq - the timestamp frequency in Hz
%   npw - number of points in each waveform
%   n - number of waveforms
%   ts - array of n timestamps (in seconds)
%   wave - array of waveforms [npw, n], raw a/d values
%   units - array of n unit assignments for each waveform, 0 = unsorted
% OUTPUT:
%   nWritten - number of waveforms successfully written to the output file
%
% This function will write a .plx file that consists of a single channel
% of spike data only. It is not intended to solve the general problem of
% writing out .plx files from Matlab. This code is intended to provide a
% starting point that can be customized as needed.
%
% Some of the limitations of this function include:
% - There is no facility to write multiple channels, slow a/d data, or event data.
% - Lots of assumptions are made (no comment, non-trodal, prethresh = npw/4,
% the header date is today, etc)
% - The count arrays in the header (TSCounts, WFCounts, EVCounts) are not
% filled (this is fixable with PlexUtil)
% - the gain for the channel is set to 1, which will result in the
% reconstructed voltages for the spikes to be different.
% - unit templates are not written out, they can be reconstructed with OFS
% - long timestamps are not supported
%
% Note that the PlexUtil program from Plexon can merge the output files
% from this routine together to form multi-channel files.
%
%

if(length(filename) == 0)
   [fname, pathname] = uigetfile('*.plx', 'Select a plx file');
	filename = strcat(pathname, fname);
end

pad256(1:256) = uint8(0);

% create the file and write the file header

plx_id = fopen(filename, 'w');
fwrite(plx_id, 1480936528, 'integer*4');    % 'PLEX' magic code
fwrite(plx_id, 101, 'integer*4');           % the version no.
fwrite(plx_id, pad256(1:128), 'char');      % placeholder for comment
fwrite(plx_id, freq, 'integer*4');          % timestamp frequency
fwrite(plx_id, 1, 'integer*4');             % no. of DSP channels
fwrite(plx_id, 0, 'integer*4');             % no. of event channels
fwrite(plx_id, 0, 'integer*4');             % no. of A/D (slow-wave) channels
fwrite(plx_id, npw, 'integer*4');           % no. points per waveform
fwrite(plx_id, npw/4, 'integer*4');         % (fake) no. pre-threshold points
[YR, MO, DA, HR, MI, SC] = datevec(now);    % current date & time
fwrite(plx_id, YR, 'integer*4');            % year
fwrite(plx_id, MO, 'integer*4');            % month
fwrite(plx_id, DA, 'integer*4');            % day
fwrite(plx_id, HR, 'integer*4');            % hour
fwrite(plx_id, MI, 'integer*4');            % minute
fwrite(plx_id, SC, 'integer*4');            % second
fwrite(plx_id, 0, 'integer*4');             % fast read (reserved)
fwrite(plx_id, freq, 'integer*4');          % waveform frequency
fwrite(plx_id, max(ts)*freq, 'double');     % last timestamp
fwrite(plx_id, pad256(1:56), 'char');       % should make 256 bytes

% now the count arrays (with counts of zero)

for i = 1:40
  fwrite(plx_id, pad256(1:130), 'char');    % first 20 are TSCounts, next 20 are WFCounts
end
for i = 1:8
  fwrite(plx_id, pad256(1:256), 'char');    % all of these make up EVCounts
end

% now the single PL_ChanHeader

% assume simple channel names
[DSPname, errmsg] = sprintf('DSP%d',ch);
[SIGname, errmsg] = sprintf('SIG%d',ch);

fwrite(plx_id, DSPname, 'char');
fwrite(plx_id, pad256(1:32-length(DSPname)));
fwrite(plx_id, SIGname, 'char');
fwrite(plx_id, pad256(1:32-length(SIGname)));
fwrite(plx_id, ch, 'integer*4');            % DSP channel number
fwrite(plx_id, 0, 'integer*4');             % waveform rate limit (not used)
fwrite(plx_id, ch, 'integer*4');            % SIG associated channel number
fwrite(plx_id, ch, 'integer*4');            % SIG reference  channel number
fwrite(plx_id, 1, 'integer*4');             % dummy for gain
fwrite(plx_id, 0, 'integer*4');             % filter off
fwrite(plx_id, -12, 'integer*4');           % (fake) detection threshold value
fwrite(plx_id, 0, 'integer*4');             % sorting method (dummy)
fwrite(plx_id, 0, 'integer*4');             % number of sorted units
for i = 1:10
  fwrite(plx_id, pad256(1:64), 'char');     % filler for templates (5 * 64 * short)
end
fwrite(plx_id, pad256(1:20), 'char');       % template fit (5 * int)
fwrite(plx_id, npw, 'integer*4');           % sort width (template only)
fwrite(plx_id, pad256(1:80), 'char');       % boxes (5 * 2 * 4 * short)
fwrite(plx_id, 0, 'integer*4');             % beginning of sorting window
fwrite(plx_id, pad256(1:128), 'char');      % comment
fwrite(plx_id, pad256(1:44), 'char');       % padding

% now the spike waveforms, each preceded by a PL_DataBlockHeader

for ispike = 1:n

  fwrite(plx_id, 1, 'integer*2');           % type: 1 = spike
  fwrite(plx_id, 0, 'integer*2');           % upper byte of 5-byte timestamp
  fwrite(plx_id, ts(ispike)*freq, 'integer*4');  % lower 4 bytes
  fwrite(plx_id, ch, 'integer*2');          % channel number
  fwrite(plx_id, units(ispike), 'integer*2');  % unit no. (0 = unsorted)
  fwrite(plx_id, 1, 'integer*2');           % no. of waveforms = 1
  fwrite(plx_id, npw, 'integer*2');         % no. of samples per waveform

  fwrite(plx_id, wave(ispike, 1:npw), 'integer*2');

end

nWritten = n;

fclose(plx_id);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% minCEntropy clustering: partitional clustering using the minimum conditional Entropy objective
%(C) Nguyen Xuan Vinh, 2010. Contact:  vinh.nguyenx@gmail.com, vinh.nguyen@monash.edu
%Input:
%   a: data, rows for objects, cols for features
%   K: number of desired clusters
%   sigma_factor: default kernel with sigma_0, specify sigma_factor to obtain
%               a new kernel width,  sigma=sigma_0/sigma_factor
%   n_run: optional, number of runs, default: 1
%   init_mem: optional, initial clustering
%Output:
%   max_mem: best clustering over n_run runs
%   max_obj: max objective value
%   S: kernel similarity matrix
%Reference:
%   [1] N. X. Vinh, Epps, J., "minCEntropy: a Novel Information Theoretic Approach for the Generation of
%       Alternative Clusterings,"  in IEEE Int. Conf. on Data Mining (ICDM) 2010.
%Example:
%
% X = [randn(100,2)+5;randn(100,2)+[5*ones(100,1) -5*ones(100,1)];  randn(100,2)-5;randn(100,2)+[-5*ones(100,1) 5*ones(100,1)]];
% [mem]=minCEntropy(X,2,1,10);  %% run minCEntropy+ 10 times,
%                               %% with K=2, sigma=sigma_0
% figure;scatter(X(:,1),X(:,2),30,mem);title('minCEntropy clustering');

function [max_mem,max_obj,all_mem,S]=minCEntropy(a,K,sigma_factor,n_run,SS,init_mem)

if nargin<3 sigma_factor=1;end;
if nargin<4 n_run   = 1;end;
if nargin<6 hasC=-1;end;

[n dim]=size(a);

if nargin<5
    SE=sqdistance(a');
    sigma0=sum(sum(sqrt(SE)))/n^2/2; %1/2 average pairwise distance
    sigma0=real(sigma0);
    sigma=sigma0/sigma_factor;
    sig2=4*sigma^2;
    S=exp(-SE/sig2);
else
    S=SS; %preprovided kernel matrix
end
Spc=zeros(n,K);            %point->cluster similarity
G=zeros(1,K);              %cluster quality
Nj=zeros(1,K);             %cluster size

max_obj=-inf;%best objective
max_mem=[];
all_mem=zeros(n_run,n);

for run=1:n_run
    change_count=0;

    if hasC==-1
        %initialization with Kmeans, using only 1 iteration
        mem=kmeans(a,K,'Maxiter',0,'EmptyAction','singleton');
        %random initialization
        %mem=round(rand(1,n)*(K-1))+1;
    else
        n_run=1;
        mem=init_mem;
    end;

    setup();

    obj=sum(G./Nj);
    isContinue=1;
while isContinue
isContinue=0;
for i=1:n

   cur_clus=mem(i);
   %check point->cluster similarity
   max_inc=-inf; % maximum objective increase
   for new_clus=1:K
      if new_clus==cur_clus continue;end;
      cond2=G(cur_clus)/(Nj(cur_clus)-1)/Nj(cur_clus)-G(new_clus)/(Nj(new_clus)+1)/Nj(new_clus)-2*Spc(i,cur_clus)/(Nj(cur_clus)-1)+2*Spc(i,new_clus)/(Nj(new_clus)+1);
      if cond2>max_inc
        max_inc=cond2;
        max_clus=new_clus;
      end
   end

    if(max_inc>0) %make change
         new_clus=max_clus;
         change_count=change_count+1;
         isContinue=1;

         %update tables
         for t=1:n
            if t==i continue;end;
            Spc(t,cur_clus)=Spc(t,cur_clus)-S(t,i);
            Spc(t,new_clus)=Spc(t,new_clus)+S(t,i);
         end

         G(cur_clus)=G(cur_clus)-2*Spc(i,cur_clus);
         G(new_clus)=G(new_clus)+2*Spc(i,new_clus);

         %update membership
         mem(i)=new_clus;
         Nj(cur_clus)=Nj(cur_clus)-1;
         Nj(new_clus)=Nj(new_clus)+1;

         CC(i,cur_clus)=false;
         CC(i,new_clus)=true;

         cur_clus=new_clus;

         change_count=change_count+1;

    end %make change

end%for i=1:n  one round through the data set
end%while point can still move

obj=sum(G./Nj);
fprintf('Sigma factor: %d, changes: %d, quality: %f\n',sigma_factor,change_count,obj);
if max_obj<obj
   max_obj=obj;
   max_mem=mem;
end
all_mem(run,:)=mem;
end%for run

fprintf('>>>>>>>>Finished clustering. best quality: %f\n',max_obj);
mem=max_mem;

%----------end main function------------------
    function setup()
        for i=1:n
           for j=1:K
                Spc(i,j)=sum(S(i,mem==j));  % point -> cluster similarity
           end
        end

        for j=1:K
            G(j)=sum(Spc(mem==j,j));
            Nj(j)=sum(mem==j);
        end
    end

end%main function

author: Sable

repo: mcbench-benchmarks


text:
function ct = ChannelTag(scnl)
   %ChannelTag Converter from scnlobject to ChannelTag
  ct = reshape([scnl.tag],size(scnl));
end

author: geoscience-community-codes

repo: GISMO


text:
% The COBRA Toolbox: testPruneModel.m
%
% Purpose:
%     - test the pruneModel function
%
% Author:
%     - Loic Marx, December 2018

% save the current path
currentDir = pwd;

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

% load reference data
load('refData_pruneModel.mat')

% define input
model = getDistributedModel('ecoli_core_model.mat');

% retrieve the default biomass reaction
biomassRxn = checkObjective(model);
biomassRxn = biomassRxn{1};

minGrowth = 1; % this value is for testing only

% generate data
[modelUpdated, modelPruned, Ex_Rxns] = pruneModel(model, minGrowth, biomassRxn);

% if biomassRxn is not define
[modelUpdated_2inputs, modelPruned_2inputs, Ex_Rxns_2inputs] = pruneModel(model, minGrowth);

% comparison between refData and generated data
assert(isequal(modelUpdated_ref, modelUpdated))
assert(isequal(modelPruned_ref, modelPruned))
assert(isequal(Ex_Rxns_ref, Ex_Rxns))

% comparison between refData and data generated without the biomassRxn input
assert(isequal(modelUpdated_ref, modelUpdated_2inputs))
assert(isequal(modelPruned_ref, modelPruned_2inputs))
assert(isequal(Ex_Rxns_ref, Ex_Rxns_2inputs))

author: opencobra

repo: cobratoolbox


text:
function laguerre_general_values_test ( )

%*****************************************************************************80
%
%% LAGUERRE_GENERAL_VALUES_TEST demonstrates the use of LAGUERRE_GENERAL_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGUERRE_GENERAL_VALUES_TEST:\n' );
  fprintf ( 1, '  LAGUERRE_GENERAL_VALUES stores values of\n' );
  fprintf ( 1, '  the generalized Laguerre function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N     A    X             L(N,A)(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, a, x, fx ] = laguerre_general_values ( n_data );

    if ( n_data == 0 )
      break
    end

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

  end

  return
end


text:
function cleanClass = cleanSingletonWhitePoints(class);
%
%Author:  RFD, BW
%Date:
%Purpose:
%   Remove singleton white matter points from each plane
%

% Mask will find points that are singletons
mask = -1*ones(3,3);
mask(2,2) = 1;

%
im = class.data;
imSize = size(im);
cleanIm = im;

%

changes = [];
for ii=1:imSize(1);
    curSlice = squeeze(cleanIm(ii,:,:));
    tst = conv2(curSlice,mask,'same');
    changes = [changes, sum(sum(tst > 0))];
    curSlice(tst > 0) = class.type.unknown;
    cleanIm(ii,:,:) = curSlice;
end
changes, sum(changes)

changes = [];
for ii=1:imSize(2);
    curSlice = squeeze(cleanIm(:,ii,:));
    tst = conv2(curSlice,mask,'same');
    changes = [changes, sum(sum(tst > 0))];
    curSlice(tst > 0) = class.type.unknown;
    cleanIm(:,ii,:) = curSlice;
end
changes, sum(changes)

changes = [];
for ii=1:imSize(3);
    curSlice = squeeze(cleanIm(:,:,ii));
    tst = conv2(curSlice,mask,'same');
    changes = [changes, sum(sum(tst > 0))];
    curSlice(tst > 0) = class.type.unknown;
    cleanIm(:,:,ii) = curSlice;
end
changes, sum(changes)

% Assign results to returned structure
cleanClass = class;
cleanClass.data = cleanIm;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = sin(varargin)
%SIN (overloaded)

switch class(varargin{1})

    case 'double'
        error('Overloaded SDPVAR/SIN CALLED WITH DOUBLE. Report error')

    case 'sdpvar'
        varargout{1} = InstantiateElementWiseUnitary(mfilename,varargin{:});
        %varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        % General operator
        operator = struct('convexity','none',...
            'monotonicity','none',...
            'definiteness','none',...
            'model','callback');

        operator.bounds     = @bounds;
        operator.convexhull = @convexhull;
        operator.derivative = @(x)(cos(x));
        operator.range = [-1 1];
        operator.domain = [-inf inf];

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

    otherwise
        error('SDPVAR/SIN called with CHAR argument?');
end

function [L,U] = bounds(xL,xU)
if xU-xL >= 2*pi
    L = -1;
    U = 1;
else
    n = floor(( (xL + xU)/2/(2*pi)));
    xL = xL - n*2*pi;
    xU = xU - n*2*pi;
    yL = sin(xL);
    yU = sin(xU);
    L = min([yL yU]);
    U = max([yL yU]);
    if (xL<pi/2 & xU>pi/2) |  (xL<-3*pi/2 & xU>-3*pi/2)
        U = 1;
    end
    if (xL < 3*pi/2 & xU > 3*pi/2) | (xL < -pi/2 & xU > -pi/2)
        L = -1;
    end
end

function [Ax, Ay, b] = convexhull(xL,xU)
if sin(xL)>=0 & sin(xU)>=0 & xU-xL<pi
    xM = (xL+xU)/2;
    fL = sin(xL);
    fM = sin(xM);
    fU = sin(xU);
    dfL = cos(xL);
    dfM = cos(xM);
    dfU = cos(xU);
    [Ax,Ay,b] = convexhullConcave(xL,xM,xU,fL,fM,fU,dfL,dfM,dfU);
elseif sin(xL)<=0 & sin(xU)<=0 & xU-xL<pi
    fL = sin(xL);
    fU = sin(xU);
    dfL = cos(xL);
    dfU = cos(xU);
    [Ax,Ay,b] = convexhullConvex(xL,xU,fL,fU,dfL,dfU);
else
    [Ax,Ay,b] = convexhullGeneral(xL,xU,@sin);
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function p_ = ViewImpliedVol(X,p)

[J,K]=size(X);

% constrain probabilities to sum to one...
Aeq = ones(1,J);
beq=1;

% ...constrain the expectation...
V = [X(:,12) - X(:,11)];
m=mean(V);
s=std(V);

A = V';
b = m-s;

% ...compute posterior probabilities
p_ = EntropyProg(p,A,b,Aeq ,beq);


text:
function KeggSpeciespKa = getKeggpKas(target_cids, target_inchi, n_pkas)
if nargin < 3
    n_pkas = 20;
end

if ismac
    cxcalc_cmd = '/Applications/ChemAxon/JChem/bin/cxcalc';
    babel_cmd = '/usr/local/bin/babel';
else
    cxcalc_cmd = 'cxcalc';
    babel_cmd = 'babel';
end

[success, ~] = system(cxcalc_cmd);
if success ~= 0
    error('Please make sure the command line program "babel" is installed and in the path');
end

KeggSpeciespKa = [];

for i = 1:length(target_cids)
    cid = target_cids(i);
    inchi = target_inchi{i};
    if isempty(inchi)
        continue
    end

    if ispc
        [success, smiles] = system(['echo ' inchi ' | ' babel_cmd ' -iinchi -osmi']);
    else
        [success, smiles] = system(['echo "' inchi '" | ' babel_cmd ' -iinchi -osmi']);
    end

    if success == 0
       smiles = strtok(smiles);
       structure = smiles;
    else
        structure = inchi;
    end

    fprintf('Using cxcalc on C%05d: %s\n', cid, structure);

    if ispc
        cmd = [cxcalc_cmd ' "' structure '" pka -a ' num2str(n_pkas) ' -b ' num2str(n_pkas) ' majorms -M true --pH 7'];
    else
        cmd = ['echo "' structure '" | ' cxcalc_cmd ' pka -a ' num2str(n_pkas) ' -b ' num2str(n_pkas) ' majorms -M true --pH 7'];
    end
    [success, cxcalc_stdout] = system(cmd);

    if ~isempty(strfind(cxcalc_stdout,'Inconsistent molecular structure'))
       success = -1;
    end

    if success == 0
        %fprintf(cxcalc_stdout);
        pkalist = regexp(cxcalc_stdout,'\n','split');
        titles = regexp(pkalist{1,1}, '\t', 'split');
        vals = regexp(pkalist{1,2}, '\t', 'split');
        if all(cellfun(@isempty,vals(2:end)))
           vals = cell(1,2*n_pkas + 3);
        end

        inds = zeros(2*n_pkas, 1);
        for i = 1:n_pkas
            inds(2*i-1, 1) = find(strcmp(titles, ['apKa' num2str(i)]));
            inds(2*i, 1) = find(strcmp(titles, ['bpKa' num2str(i)]));
        end

        pkalist = vals(1, inds);
        pkalist = regexprep(pkalist, ',', '\.');
        pkalist = str2double(pkalist);
        pkalist = sort(pkalist,'descend');
        pkalist = pkalist(pkalist>=0 & pkalist<=14);

        % find the nH and charge of the major macrospecies
        ind = find(strcmp(titles, 'major-ms'));
        majorms_smiles = vals{1, ind};
        if isempty(majorms_smiles)
            majorms_smiles = smiles;
        end

        if ispc
            cmd = ['echo ' majorms_smiles ' | babel -ismiles -oinchi ---errorlevel 0 -xFT/noiso'];
        else
            cmd = ['echo "' majorms_smiles '" | babel -ismiles -oinchi ---errorlevel 0 -xFT/noiso'];
        end
        [success, babel_stdout] = system(cmd);
        if success == 0 && ~isempty(babel_stdout) && strcmp('InChI=',babel_stdout(1:6))
            majorms_nstd_inchi = strtok(babel_stdout);
            [~, nH, charge] = getFormulaAndChargeFromInChI(majorms_nstd_inchi);
        else
            nH = 0;
            charge = 0;
        end

        idx = length(KeggSpeciespKa) + 1;
        if ~isempty(pkalist)
            pkas = zeros(length(pkalist)+1,length(pkalist)+1);
            pkas(2:end,1:end-1) = diag(pkalist);
            pkas = pkas + pkas';
            KeggSpeciespKa(idx).pKas = pkas;

            mmsbool = false(size(pkas,1),1);
            if any(pkalist <= 7)
                mmsbool(find(pkalist <= 7,1)) = true;
            else
                mmsbool(end) = true;
            end
            KeggSpeciespKa(idx).majorMSpH7 = mmsbool;
            zs = 1:size(pkas,1);
            zs = zs - find(mmsbool);
            zs = zs + charge;
            KeggSpeciespKa(idx).zs = zs;

            nHs = 1:size(pkas,1);
            nHs = nHs - find(mmsbool);
            nHs = nHs + nH;
            KeggSpeciespKa(idx).nHs = nHs;
            KeggSpeciespKa(idx).cid = cid;
        else
            KeggSpeciespKa(idx).pKas = [];
            KeggSpeciespKa(idx).majorMSpH7 = true;
            KeggSpeciespKa(idx).zs = charge;
            KeggSpeciespKa(idx).nHs = nH;
            KeggSpeciespKa(idx).cid = cid;
        end
    end
end
KeggSpeciespKa = KeggSpeciespKa';


text:
% Runs a step of non-local Bayes
%
% USAGE: [deno, aggw] = nlbayes_step(nisy, bsic, sigma, prms)
%
%  -> nisy    : noisy image
%  -> bsic    : basic estimate (can be empty for the first step)
%  -> sigma   : noise std. dev.
%  -> prms    : struc with prms (wx, px, np, r)
%
%  <- deno    : denoised image
%  <- aggw    : aggregation weights
function [deno, aggw] = nlbayes_step(nisy, bsic, sigma, prms)

% image size and channels
w = size(nisy,2);
h = size(nisy,1);
chnls = size(nisy,3);

% in the absence of the basic estimate, define basic as the noisy image
if isempty(bsic),
	bsic = nisy;
	step2 = false;
else
	step2 = true
end

% aggregation
aggw = zeros(size(nisy)); % aggregation weights
aggp = zeros(size(nisy)); % mask of already processed patches (for speed-up trick)
aggu = zeros(size(nisy)); % aggregated image
deno = zeros(size(nisy)); % denoised image (deno = aggu ./ aggw);

% step sizes
stepx = floor(prms.px/2);
stepy = floor(prms.px/2);
ii = 0; % iteration counter

% patch dimensionality
pdim = prms.px * prms.px * chnls;

% use dct basis
% U = kron(eye(chnls), dct_basis(prms.px, prms.px));
% use pca basis
U = [];

% use an aggregation window
if prms.pw,
	wx = chebwin(prms.px);
else
	wx = ones(prms.px,1);
end
wwx = repmat(reshape(wx*wx',[prms.px*prms.px 1]),[chnls prms.np]);

% main loop
for pay = [1:stepy:h - prms.px+1,h - prms.px+1],
for pax = [1:stepx:w - prms.px+1,w - prms.px+1], ii = ii + 1;

	% acceleration: skip iteration if patch has been been already denoised
	if aggp(pay,pax,1),
		continue;
	end

	% -------------------------------------------------- compute patch group

	% patches in search region
	srch_region = bsic(max(1,pay - prms.wx):min(h,pay + prms.wx + prms.px - 1),...
	                   max(1,pax - prms.wx):min(w,pax + prms.wx + prms.px - 1),:);

	srch_patches = im2col_ch(srch_region, [prms.px prms.px]);
	refe_patch = bsic(pay:pay+prms.px-1,pax:pax+prms.px-1,:);

	[distances, idx] = sort(L2_distance(refe_patch(:), srch_patches));

	% coordinates of the np nearest neighbors to the ref patch
	idx = idx(1:prms.np)';

	srch_h = size(srch_region,1) - prms.px + 1;
	patches.coords = [max(1,pax - prms.wx) + floor((idx-1)/srch_h),...
	                  max(1,pay - prms.wx) +   mod( idx-1 ,srch_h)];

	% ---------------------------------------------- extract similar patches
	if step2,
		patches.bsic = srch_patches(:,idx);

		% extract noisy patches
		srch_region = nisy(max(1,pay - prms.wx):min(h,pay + prms.wx + prms.px - 1),...
		                   max(1,pax - prms.wx):min(w,pax + prms.wx + prms.px - 1),:);
		patches.nisy = im2col_ch(srch_region, [prms.px prms.px]);
		patches.nisy = patches.nisy(:,idx);
	else
		patches.nisy = srch_patches(:,idx);
	end

	% ----------------------------------------------- compute bayes estimate
	if step2,
		aa = reshape(patches.bsic,[pdim prms.np]);
		[dd,gg] = compute_bayes_estimate(patches.nisy,aa,sigma,prms.r,'pos',U);
	else
		[dd,gg] = compute_bayes_estimate(patches.nisy,[],sigma,prms.r,'pos',U);
	end

	% ------------------------------------------- aggregate patches on image
	gg = ones(size(dd,1),1)*gg;
	aggu = aggregate_patches(aggu, wwx.*gg.*dd, [prms.px prms.px], patches.coords);
	aggw = aggregate_patches(aggw, wwx.*gg    , [prms.px prms.px], patches.coords);

	aggp(patches.coords(:,2) + (patches.coords(:,1)-1)*h) = 1+...
	           aggp(patches.coords(:,2) + (patches.coords(:,1)-1)*h);

	% draw
	if (mod(ii,50) == 1)
		nonzero = find(aggw ~= 0);
		deno(nonzero) = min(255, max(0, aggu(nonzero) ./ aggw(nonzero)));

		% draw a red box indicating limits of search region
		deno(max(1,pay-prms.wx):min(h,pay+prms.wx),max(1,pax-prms.wx),:) = 0;
		deno(max(1,pay-prms.wx):min(h,pay+prms.wx),min(w,pax+prms.wx),:) = 0;
		deno(max(1,pay-prms.wx),max(1,pax-prms.wx):min(w,pax+prms.wx),:) = 0;
		deno(min(h,pay+prms.wx),max(1,pax-prms.wx):min(w,pax+prms.wx),:) = 0;
		deno(max(1,pay-prms.wx):min(h,pay+prms.wx),max(1,pax-prms.wx),1) = 255;
		deno(max(1,pay-prms.wx):min(h,pay+prms.wx),min(w,pax+prms.wx),1) = 255;
		deno(max(1,pay-prms.wx),max(1,pax-prms.wx):min(w,pax+prms.wx),1) = 255;
		deno(min(h,pay+prms.wx),max(1,pax-prms.wx):min(w,pax+prms.wx),1) = 255;

		imagesc(max(min([deno;255-aggw],255),0)/255,[0 1]);
		axis equal, axis off,
		colormap gray
		drawnow
		pause(.01)

		% reset red pixels to 0
		deno(max(1,pay-prms.wx):min(h,pay+prms.wx),max(1,pax-prms.wx),1) = 0;
		deno(max(1,pay-prms.wx):min(h,pay+prms.wx),min(w,pax+prms.wx),1) = 0;
		deno(max(1,pay-prms.wx),max(1,pax-prms.wx):min(w,pax+prms.wx),1) = 0;
		deno(min(h,pay+prms.wx),max(1,pax-prms.wx):min(w,pax+prms.wx),1) = 0;
	end

end
end

% compute denoised image
nonzero = find(aggw ~= 0);
deno(nonzero) = min(255, max(0, aggu(nonzero) ./ aggw(nonzero)));

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ a, b ] = assemble_fem ( x_num, x, element_num, element_node, ...
  quad_num, t, k_fun, rhs_fun )

%*****************************************************************************80
%
%% ASSEMBLE_FEM assembles the finite element stiffness matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X_NUM, the number of nodes.
%
%    Input, real X(X_NUM), the coordinates of nodes.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(2,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Input, integer QUAD_NUM, the number of quadrature points used in assembly.
%
%    Input, real T, the current time.
%
%    Input, real K_FUN(), a function to evaluate the heat conductivity.
%
%    Input, real RHS_FUN(), a function to evaluate the right hand side.
%
%    Output, sparse real A(X_NUM,X_NUM), the finite element stiffness matrix.
%
%    Output, real B(X_NUM), the right hand side.
%
%  Local parameters:
%
%    Local, real BI, DBIDX, the value of some basis function
%    and its first derivative at a quadrature point.
%
%    Local, real BJ, DBJDX, the value of another basis
%    function and its first derivative at a quadrature point.
%

%
%  Initialize the arrays.
%
  b = zeros ( x_num, 1 );
  a = sparse ( x_num, x_num );
%
%  Get the quadrature weights and nodes.
%
  [ reference_w, reference_q ] = quadrature_set ( quad_num );
%
%  Consider each ELEMENT.
%
  for element = 1 : element_num

    element_x(1:2,1) = x(element_node(1:2,element));

    element_q(1:quad_num,1) = reference_to_physical ( element, element_node, ...
      x, quad_num, reference_q );

    element_area = element_x(2) - element_x(1);

    element_w(1:quad_num,1) = ( element_area / 2.0 ) * reference_w(1:quad_num);
%
%  Consider the QUAD-th quadrature point in the element.
%
    k_value = k_fun ( quad_num, element_q, t );
    rhs_value = rhs_fun ( quad_num, element_q, t );

    for quad = 1 : quad_num
%
%  Consider the TEST-th test function.
%
%  We generate an integral for every node associated with an unknown.
%
      for i = 1 : 2

        test = element_node(i,element);

        [ bi, dbidx ] = basis_function ( test, element, x, 1, element_q(quad) );

        b(test) = b(test) + element_w(quad) * rhs_value(quad) * bi;
%
%  Consider the BASIS-th basis function, which is used to form the
%  value of the solution function.
%
        for j = 1 : 2

          basis = element_node(j,element);

          [ bj, dbjdx ] = basis_function ( basis, element, x, 1, element_q(quad) );

          a(test,basis) = a(test,basis) + element_w(quad) * ( ...
            + k_value(quad) * dbidx * dbjdx );

        end

      end

    end

  end

  return
end


text:
% MATGEOM Geometric Computing Toolbox
% Version 1.0 26-07-2017 .
%
%   MatGeom Provides low-level functions for geometric computing. It is
%   possible to create, display, compute intersections... of various
%   geometrical primitives, in 2D and 3D.
%
%   The library is organized into several modules:
%   geom2d              - General function in euclidean plane
%   polygons2d          - Functions operating on point lists
%   graphs              - Manipulation of geometric graphs
%   polynomialCurves2d  - Representation of smooth polynomial curves
%   geom3d              - General function in 3D euclidean space
%   meshes3d            - Manipulation of 3D surfacic meshes
%
%   Type 'help(MODULENAME)' for further info.
%
%   To install the library, with all sub-directories, run the script
%   'setupMatGeom'.
%
%   More information on the project homepage:
%   https://github.com/mattools/matGeom
%
%   Online documentation:
%   https://github.com/mattools/matGeom
%
% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2011-03-21,    using Matlab 7.9.0.529 (R2009b)
% Project homepage: http://github.com/mattools/matGeom
% http://www.pfl-cepia.inra.fr/index.php?page=geom3d
% Copyright 2011 INRA - Cepia Software Platform.

help(mfilename);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [mmps2] = mps22mmps2(mps2)
% Convert acceleration from meters per square-second to millimeters per
% second squared.
% Chad A. Greene 2012
mmps2 = mps2*1e+3;


text:
function test_nint_test01 ( )

%*****************************************************************************80
%
%% TEST01 just prints out the title information.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
  clear

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  GET_PROBLEM_NUM returns the number of problems.\n' );
  fprintf ( 1, '  P00_NAME(#) returns the name for problem #.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We use these two routines to print a directory\n' );
  fprintf ( 1, '  of all the problems.\n' );

  problem_num = get_problem_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of problems available is %d\n', problem_num );
  fprintf ( 1, '\n' );

  for problem = 1 : problem_num

    name = p00_name ( problem );
    fprintf ( 1, '  %8d  "%s";\n', problem, name );

  end

  return
end


text:
function value = year_length_months_ethiopian ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_MONTHS_ETHIOPIAN returns the number of months in an Ethiopian year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of months in the year.
%
  value = 13;

  return
end


text:
function r8but_print ( n, mu, a, title )

%*****************************************************************************80
%
%% R8BUT_PRINT prints a R8BUT matrix.
%
%  Discussion:
%
%    The R8BUT storage format is used for a banded upper triangular matrix.
%    The matrix is assumed to be zero above the MU-th superdiagonal.
%    The matrix is stored in an MU+1 by N array.
%    Columns are preserved.
%
%    The diagonal is stored in row MU+1 of the array.
%    The first superdiagonal in row MU, columns 2 through N.
%    The second superdiagonal in row MU-1, columns 3 through N.
%    The MU-th superdiagonal in row 1, columns MU+1 through N.
%
%  Example:
%
%    N = 5, MU = 2
%
%    A11 A12 A13   0   0
%      0 A22 A23 A24   0
%      0   0 A33 A34 A35
%      0   0   0 A44 A45
%      0   0   0   0 A55
%                --- ---
%                    ---
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer MU, the upper bandwidth.
%
%    Input, real A(MU+1,N), the R8BUT matrix.
%
%    Input, string TITLE, a title to be printed.
%
  r8but_print_some ( n, mu, a, 1, 1, n, n, title );

  return
end


text:
function H = ldivide(F, G)
%.\   Pointwise left divide for two CHEBFUN3V objects.
%
% See also CHEBFUN3V/RDIVIDE.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

if ( ( isempty(F) ) || ( isempty(G) ) )
    H = chebfun3v;
    return
end

if ( isa(F, 'double') || isa(F, 'chebfun3') ) % CHEBFUN3 .\ CHEBFUN3V
    nG = G.nComponents;
    Gc = G.components;
    H = G;
    for j = 1:nG
        H.components{j} = ldivide(F, Gc{j});
    end

elseif ( isa(G, 'double') || isa(G, 'chebfun3'))  % CHEBFUN3V .\ CHEBFUN3
    nF = F.nComponents;
    Fc = F.components;
    H = F;
    for j = 1:nF
        H.components{j} = ldivide(Fc{j}, G);
    end

else                                     % CHEBFUN3V .\ CHEBFUN3V
    nF = F.nComponents;
    nG = G.nComponents;
    if ( nF ~= nG )
        error('CHEBFUN:CHEBFUN3V:ldivide:dimensionMismatch', ...
            'CHEBFUN3V do not have the same number of components.')
    end
    H = F;
    for j = 1:nF
        H.components{j} = ldivide(F.components{j}, G.components{j});
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function mtrFigureGASMatchedPathDist(machineDir, imageDir)

subjVec = {'tony_nov05','thor_nov05','sil_nov05'};
%subjVec = {'thor_nov05'};
threshVec = [1000];
ctFile = 'paths500k_5k.dat';
sttFile = 'stt_clip.dat';
distMatricesFile = 'distMatricesGAS.mat';
bSavePlot = 0;

for ss = 1:length(subjVec)
    subjDir = fullfile(machineDir,subjVec{ss});
    figure;

    fgDir = fullfile(subjDir, 'conTrack');
    disp(['cd ' fgDir]);
    cd(fgDir);
    [meanD, maxD] = mtrMatchPathways(fullfile(subjDir,'dt6.mat'), fullfile(subjDir,'fibers',sttFile), ctFile);

    % Save out the figures to the image dir
    subplot(2,1,1); hist((min(meanD,[],2)));
    xlabel('Mean distance to closest conTrack path (mm)');
    subplot(2,1,2); hist((min(maxD,[],2)));
    xlabel('Max distance to closest conTrack path (mm)');
    save(distMatricesFile,'meanD','maxD');

    if(bSavePlot)
        set(gcf,'Position',[504 687 436 259]);
        figFilename = fullfile(subjDir,imageDir,['distHist.png'])
        set(gcf,'PaperPositionMode','auto');
        print('-dpng', figFilename);
    end
end


text:
function [snr per] = get_per(in_vec)
% get_per.m - get PER data
global numPkts PER_snr PER_exp;

snr = PER_snr';
per = PER_exp';


text:
function test_bug1600

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_checkdata ft_datatype_source

% The problem: ft_checkdata(volume, 'datatype', 'source') does not seem to
% convert the inside back to vectorial representation

% verify the problem
volume         = [];
volume.dim     = [3 4 5];
volume.anatomy = zeros(3,4,5);
volume.anatomy(2,2:3,2:4) = rand(1,2,3);
volume.inside  = volume.anatomy>0;
volume.transform = eye(4);

source = ft_checkdata(volume, 'datatype', 'source');

assert(all(size(source.inside)==[60 1]));


text:
function [rs, rs2, rs3] = dimensions(s, bins)

%tstoolbox/@signal/dimensions
%   Syntax:
%     * [bc,in,co] = dimensions(s, bins)
%
%   Input arguments:
%     * s - data points (row vectors)
%     * bins - maximal number of partition per axis, default is 100
%
%   Output arguments:
%     * bc - scaling of boxes with partititon sizes (log[2]-log[2])
%     * in - scaling of information with partititon sizes (log[2]-log[2])
%     * co - scaling of correlation with partititon sizes (log[2]-log[2])
%
%   Compute boxcounting, information and correlation dimension of a
%   time-delay reconstructed timeseries s for dimensions from 1 to D,
%   where D is the dimension of the input vectors using boxcounting
%   approach.
%
%   Scale data to be within 0 and 1. Give a sortiment of (integer)
%   partitionsizes with almost exponential behaviour.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,2);

if ndim(s) ~= 2
    error('Signal must contain vector data');
end

if nargin<2
    bins = 100;
end

points = data(s);
[N,dim] = size(points);

% scale data to be within 0 and 1
points = points - min(min(points));
points = points / max(max(points));

% give a sortiment of (integer) partitionsizes with almost exponential behaviour
par = [2 3 4 5 6 7 8 10 12 14 16 20 23 27 32 39 46 54 64 77 91 108 128 153 182 216 256 ...
      305 363 431 512 609 725 862 1024 1218 1449 1723 2048 2436 2897 3445 4096 4871 ...
	  5793 6889 8192 9742 11586 13778 16384 19484 23171 27555 32768 38968 46341 55109 65536];

partitions = par(find(par<=bins));	% use no sizes greater than bins

[c,d,e] = boxcount(points, partitions);

c = [zeros(1,dim) ; c];             % add zeros for partition size 1
d = [zeros(1,dim) ; d];             % add zeros for partition size 1
e = [zeros(1,dim) ; e];             % add zeros for partition size 1
partitions = [1 ; partitions(:)];

a1 = achse(-log2(partitions));     		% create axis with arbitrary spacing
a1 = setname(a1, 'ld r');

a2 = setname(achse(unit, 1, 1), 'Embedding dimension');

rs = signal(core(c), s);
rs = setaxis(rs, 1, a1);
rs = setaxis(rs, 2, a2);
rs = setplothint(rs, 'multigraph');
rs = addhistory(rs,  ['Computed boxcounting dimension']);
rs = addcommandlines(rs, 's = boxdim(s', bins);
rs = setyname(rs, 'ld N(r)');
rs = setlabel(rs, 'Scaling of D0');

rs2 = signal(core(d), s);
rs2 = setaxis(rs2, 1, a1);
rs2 = setaxis(rs2, 2, a2);
rs2 = setplothint(rs2, 'multigraph');
rs2 = addhistory(rs2,  ['Computed information dimension']);
rs2 = addcommandlines(rs2, 's = infodim(s', bins);
rs2 = setyname(rs2, 'I(r)');
rs2 = setlabel(rs2, 'Scaling of D1');

rs3 = signal(core(e), s);
rs3 = setaxis(rs3, 1, a1);
rs3 = setaxis(rs3, 2, a2);
rs3 = setplothint(rs3, 'multigraph');
rs3 = addhistory(rs3,  ['Computed correlation dimension']);
rs3 = addcommandlines(rs3, 's = corrdim(s', bins);
rs3 = setyname(rs3, 'ld C(r)');
rs3 = setlabel(rs3, 'Scaling of D2');

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function wrapper_rgbd_to_ucm(imset, modelfile, ucm_dir, contours_cues_dir, C)
  sc = [2 1 0.5];
  imlist = getImageSet(imset);

  model = load(modelfile); %fullfile_ext(contours_model_dir, 'forest', modelname, 'mat'));
  model = model.model;

  exists_or_mkdir(fullfile(ucm_dir, 'multi'));
  exists_or_mkdir(fullfile(ucm_dir, 'scale_2.0'));
  exists_or_mkdir(fullfile(ucm_dir, 'scale_1.0'));
  exists_or_mkdir(fullfile(ucm_dir, 'edges_1.0'));
  exists_or_mkdir(fullfile(ucm_dir, 'scale_0.5'));
  exists_or_mkdir(fullfile(ucm_dir, 'multi-png'));

  for i = 1:length(imlist),
    id = imlist{i};
    I = getImage(id, 'images');
    D = getImage(id, 'depth');
    cacheFile = fullfile_ext(contours_cues_dir, id, 'mat');
    [E, Es, O] = detectEdge(I, D, [], C, model, sc, [], cacheFile);
    [ucm2 ucms] = contours_to_ucm(I, sc, Es, O);
    parsave(fullfile_ext(ucm_dir, 'multi', id, 'mat'), 'ucm2', ucm2);
    parsave(fullfile_ext(ucm_dir, 'edges_1.0', id, 'mat'), 'E', E{2}, 'Es', Es{2}, 'O', O{2});
    parsave(fullfile_ext(ucm_dir, 'scale_2.0', id, 'mat'), 'ucm2', ucms(:,:,1));
    parsave(fullfile_ext(ucm_dir, 'scale_1.0', id, 'mat'), 'ucm2', ucms(:,:,2));
    parsave(fullfile_ext(ucm_dir, 'scale_0.5', id, 'mat'), 'ucm2', ucms(:,:,3));
    imwrite(im2uint8(ucm2(3:2:end, 3:2:end)), fullfile_ext(ucm_dir, 'multi-png', id, 'png'));
  end
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function [pb] = pbHuman(pres,iid,r)
if nargin<3, r=[0.25 0.99]; end
segs = readSegs(pres,iid);
pb = zeros(size(segs{1}));
for i = 1:numel(segs),
  bmap = seg2bmap(segs{i});
  pb = pb + bmap;
end
pb = pb / numel(segs);
pb = (pb~=0) .* (pb * (r(2)-r(1)) + r(1));
pb = min(pb,r(2));

author: Cloud-CV

repo: object-proposals


text:
function draw_spheres(time, period, S, figure_handle)
% DRAW_SPHERES - Plot the spherical obstacles

persistent scatter_handle

if mod(time, period) == 0 && S.is_active_spheres
    [~,n_obs] = size(S.S);
    s = repmat(400*600,n_obs,1);
    X = S.S(1,:)';
    Y = S.S(2,:)';
    Z = S.S(3,:)';
    colors = repmat(0.5,3, n_obs);

    if time == 0
        scatter_handle = scatter3(X, Y, Z, s, colors', '.');
        hold on

        axis square;
        view(32,47);
        xlabel('x position [m]');
        ylabel('y position [m]');
        zlabel('z position [m]');

%     else
%         set(scatter_handle, 'Xdata', Y, ...
%             'Ydata', X, ...
%             'Zdata', Z, ...
%             'Marker', '.', ...
%             'SizeData', s, ...
%             'CData', colors');

    end
end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function [ y2, m2, d2, f2 ] = yjf_to_ymdf_english ( y1, j1, f1 )

%*****************************************************************************80
%
%% YJF_TO_YMDF_ENGLISH converts an English date from YJF to YMDF format.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, J1, real F1, the YJF date.
%
%    Output, integer Y2, M2, D2, real F2,
%    the YMDF date.
%

%
%  Copy the input.
%
  y2 = y1;
  j2 = j1;
  f2 = f1;
%
%  Check the input.
%
  [ y2, j2, ierror ] = yj_check_english ( y2, j2 );

  if ( ierror ~= 0 )
    y2 = 0;
    m2 = 0;
    d2 = 0;
    f2 = 0.0;
    return
  end
%
%  Convert the input.
%
  m2 = 1;
  d2 = j2;

  [ y2, m2, d2 ] = day_borrow_english ( y2, m2, d2 );

  [ y2, m2, d2 ] = day_carry_english ( y2, m2, d2 );

  return
end


text:
function annotation = addDetector2LM(boundingBox, boxScores, boxClass, annotation)
%
% annotation = addDetector2LM(boundingBox, boxScores, boxClass, annotation)
%
% Input
%    boundingBox: matrix [Nobjects x 4]. Each row is boundingBox(i,:) = [xmin xmax ymin ymax]
%    boxScores  : vector [Nobjects]
%    boxClass   : string with the object name. I assume all bounding boxes
%                 belong to a single class
%    annotation : the object will be added to the annotation struct
%                 (optional)
%
% Output
%    annotation

nDetections = size(boundingBox, 1);

% check if there are other objects
nObjects = 0;

if nargin == 4
    if isfield(annotation, 'object')
        nObjects = length(annotation.object);
    end
else
    annotation = [];
end

% Translate bounding boxes into LabelMe format:
if nDetections>0
    k = nObjects;
    for i = 1:nDetections;%1:nDetections
        k = k+1;
        annotation.object(k).name = boxClass;
        annotation.object(k).confidence = boxScores(i);
        %X = [boundingBox(i,1) boundingBox(i,2) boundingBox(i,2) boundingBox(i,1) boundingBox(i,1)];
        %Y = [boundingBox(i,3) boundingBox(i,3) boundingBox(i,4) boundingBox(i,4) boundingBox(i,3)];

        X = [boundingBox(i,1) boundingBox(i,2) boundingBox(i,2) boundingBox(i,1)];
        Y = [boundingBox(i,3) boundingBox(i,3) boundingBox(i,4) boundingBox(i,4)];

        annotation.object(k).polygon = setLMpolygon(X,Y);
        annotation.object(k).deleted = 0;
        annotation.object(k).date = '';
        annotation.object(k).id = i;
    end
end

author: CSAILVision

repo: LabelMeToolbox


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

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

% Last Modified by GUIDE v2.5 06-Feb-2008 23:04:44

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

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

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

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

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes HELP2 wait for user response (see UIRESUME)
% uiwait(handles.figure1);

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

% Get default command line output from handles structure
varargout{1} = handles.output;


text:
%% Template Matlab script to create an BIDS compatible electrodes.tsv file
% This example lists all required and optional fields.
% When adding additional metadata please use CamelCase
%
% DHermes, 2017
% modified Jaap van der Aar 30.11.18

% Writing json files relies on the JSONio library
% https://github.com/bids-standard/bids-matlab
%
% Make sure it is in the matab/octave path
try
    bids.bids_matlab_version;
catch
    warning('%s\n%s\n%s\n%s', ...
            'Writing the JSON file seems to have failed.', ...
            'Make sure that the following library is in the matlab/octave path:', ...
            'https://github.com/bids-standard/bids-matlab');
end

%%
clear;

this_dir = fileparts(mfilename('fullpath'));
root_dir = fullfile(this_dir, '..', filesep, '..');

project = 'templates';

sub_label = '01';
ses_label = '01';
task_label = 'LongExample';
run_label = '01';

name_spec.modality = 'ieeg';
name_spec.suffix = 'channels';
name_spec.ext = '.tsv';
name_spec.entities = struct('sub', sub_label, ...
                            'ses', ses_label, ...
                            'task', task_label, ...
                            'run', run_label);

% using the 'use_schema', true
% ensures that the entities will be in the correct order
bids_file = bids.File(name_spec, 'use_schema', true);

% Contrust the fullpath version of the filename
channels_tsv_name = fullfile(root_dir, project, bids_file.bids_path, bids_file.filename);

%% make a participants table and save

%% required columns
name = {'n/a'}; % Label of the channel, only contains letters and numbers. The label must
% correspond to _electrodes.tsv name and all ieeg type channels are required to have \
% a position. The reference channel name MUST be provided in the reference column

type = {'n/a'}; % Type of channel, see below for adequate keywords in this field

units = {'n/a'}; % Physical unit of the value represented in this channel, for example: V for Volt,
% specified according to the SI unit symbol and possibly prefix symbol (for example: mV, ?V),
% see the BIDS spec (section 15 Appendix V: Units) for guidelines for Units and Prefixes.

low_cutoff = [0]; % Frequencies used for the low pass filter applied to the
% channel in Hz. If no low pass filter was applied, use n/a. Note that
% anti-alias is a low pass filter, specify its frequencies here if applicable.

high_cutoff = [0]; % Frequencies used for the high pass filter applied to
% the channel in Hz. If no high pass filter applied, use n/a.

%% recommended columns:

reference = {'n/a'}; % Specification of the reference (for example: "mastoid", "ElectrodeName01",
% "intracranial", "CAR", "other", "n/a"). If the channel is not an electrode channel
% (for example: a microphone channel) use `n/a`.

group = {'n/a'}; % Which group of channels (grid/strip/probe) this channel belongs to.
% One group has one wire and noise can be shared. This can be a name or number.
% Note that any groups specified in `_electrodes.tsv` must match those present here.

%% optional columns

sampling_frequency = [0]; % Sampling rate of the channel in Hz.

description = {'n/a'}; % Brief free-text description of the channel, or other information of
% interest (for example: position (for example: 'left lateral temporal surface', 'unipolar/bipolar', etc.)).

notch = [0]; % Frequencies used for the notch filter applied to the channel,
% in Hz. If no notch filter applied, use n/a.

status = {'bad'}; % Data quality observed on the channel (good/bad). A channel is considered bad
% if its data quality is compromised by excessive noise. Description of noise type SHOULD be
% provided in [status_description].

status_description = {'n/a'}; % Freeform text description of noise or artifact affecting data
% quality on the channel. It is meant to explain why the channel was declared bad in [status].

%% write
t = table(name, type, units, low_cutoff, high_cutoff, reference, ...
          group, sampling_frequency, description, notch, status, status_description);

writetable(t, channels_tsv_name, 'FileType', 'text', 'Delimiter', '\t');

author: bids-standard

repo: bids-starter-kit


text:
% SERVER represents a handle to a server instance.
%
% To initialize a server, use:
%
% s = Server(options);
%
% See help for Server.Server and Server.rpc for usage.
classdef Server < handle

properties
  slaveCount
  serv_sock
  slave_processes
  slave_sock
  options
  hostname
  rootpath
  savedObjects
end

methods

function server = Server(options)
% SERVER.SERVER is the constructor for a server instance.
%
% This initializes the server instance using the provided options
% (documented below).  Server() will block until all of the slaves have
% connected.  After it returns, use Server.rpc() to call remote functions.
%
%'options' is a struct containing configuration information and
% must include at least one of options.slavelist or options.slavecount.
% You can launch the slave instances of MATLAB by hand, or use the
% slaveRun.sh script.  The function returns a handle to the
% server state.  The server object is then used to perform operations on the slaves.
%
% OPTIONS:
%
% options.slavecount
%
% options.port (default=10000)
% Port on which server will listen for slave connections.  If the server gets orphaned or you
% didn't destroy the server properly with serverDestroy(), then this may complain that a port
% is in use.  Restarting the server MATLAB instance will fix the problem.
%
% options.timeout (default=5000)
% Time (in milliseconds) for timeout of server socket.  This affects how frequently you
% can Ctl-C the server while it's waiting for slaves.  The clock resets each time a slave
% connects.
%
% options.retries (default=forever)
% Number of times to retry the accept() call to wait for a client.  By default, the server
% waits forever (i.e., until the user Ctl-C's the function or the clients all connect).  For
% batch jobs, don't forget to set this to something reasonable!
%
% options.slavedir (default = `pwd`)
% Working directory on the slave.  The slaves will all 'cd' to this directory after connecting.
%
% options.startslaves
% If nonzero, serverInit will attempt to start the slaves on its own.  This is presently
% a bit unstable but will hopefully make batch jobs easier in the future.
%
% options.slavepath (default = setup_paths())
% When options.startslaves is set, specifies a single path string passed to addpath() before
% the slaves are started.  This is usually the path to the matlabserver/ directory (which includes
% the slaveRun.m file needed by the slaves.)
%
% options.slavelist
% When options.startslaves is set, contains a cell array of host names of all the slaves to start.

  import java.io.*;
  import java.net.ServerSocket;
  import java.net.Socket;

  rootpath=setup_paths();

  if (~isfield(options, 'emulate'))
    options.emulate = 0;
  end
  if (~isfield(options, 'startslaves'))
    options.startslaves = 0;
  end
  if (~isfield(options, 'matlabpool'))
    options.matlabpool = 0;
  end
  if (~isfield(options, 'timeout'))
    options.timeout = 5000; % 5 second timeout
  end
  if (~isfield(options, 'retries'))
    options.retries = 100000; % try for a really long time
  end
  if (~isfield(options, 'port'))
    options.port = 10000;
  end
  if (~isfield(options, 'slavedir'))
    options.slavedir=pwd;
  end
  if (~isfield(options, 'slavepath'))
    options.slavepath=rootpath;
  end
  if (~isfield(options, 'slavecount'))
    if (~isfield(options, 'slavelist'))
      error('Must provide options.slavelist cell array or options.slavecount.');
    elseif (options.emulate)
      error('Must provide options.slavecount when "emulate" is used.');
      options.slavecount = length(options.slavelist);
    end
  else
    if (isfield(options, 'slavelist'))
      if (options.slavecount ~= length(options.slavelist))
        error('options.slavelist doesn''t match options.slavecount.');
      end
    end
  end
  if (options.startslaves)
    if (~isfield(options, 'slavelist'))
      error('Must provide options.slavelist cell array if using options.startslaves.');
    end
  end
  slaveCount = options.slavecount;
  server.slaveCount=slaveCount;
  server.options = options;
  server.rootpath = rootpath;

  if (options.emulate)
    % do nothing.

  else % emulate == 0

  % get local hostname
  addr=java.net.InetAddress.getLocalHost();
  hostname=char(addr.getHostName());
  server.hostname=hostname;

  % Open server socket
  server.serv_sock = ServerSocket(options.port, slaveCount);
  server.serv_sock.setReuseAddress(true);
  if (options.timeout > 0)
    server.serv_sock.setSoTimeout(options.timeout);
  end

  try

    % start up slaves
    if (options.startslaves)
      RT = java.lang.Runtime.getRuntime();
      for i = 1:slaveCount
        if (strcmp(options.slavelist{i}, 'localhost') == 1)
          cmd=sprintf('matlab -nodesktop -nosplash -r addpath(''%s''),slaveRun(''%s'',%d),quit;', options.slavepath, hostname, options.port);
          server.slave_process{i}=RT.exec(cmd, '', File(options.slavedir));
        else
          cmd=sprintf('ssh %s "matlab -nodesktop -nosplash -r addpath(''%s''),slaveRun(''%s'',%d),quit;"', options.slavelist{i}, options.slavepath, hostname, options.port);
          server.slave_process{i}=RT.exec(cmd);
        end
      end
    end

    % wait for slaves to connect; NOTE: slave_sock{i} and slave_process{i}
    % may not correspond.
    for i=1:slaveCount
      for r=1:options.retries
        try
          % accept connection from slave
          server.slave_sock{i} = server.serv_sock.accept();
          break ;
        catch
          % timeout
          if (r < options.retries)
            fprintf(1,'Waiting for slave to connect... %d/%d\n',i, server.slaveCount);
          else
            display('Timed out waiting for slaves.');
            rethrow(lasterror);
          end
        end
      end
    end

    % all are connected; close listener.
    server.serv_sock.close();

    try
      % exchange info with clients
      display('Handshaking...');
      slaveinfo.workdir = options.slavedir;
      server.rpc('handshake', {slaveinfo});

      % start matlabpool sequentially...
      if (options.matlabpool ~= 0)
        for i=1:server.slaveCount
          fprintf(1, 'Starting matlabpool on slave %d...', i);
          server.sendRequest(i, 0, 'matlabpool');
          dummy=server.getReply(i);
        end
      end
    catch
      display('Handshake failed.');
      rethrow(lasterror);
    end

  catch
    delete(server);
    server=[];
    rethrow(lasterror);
  end

  end % if options.emulate == 0
end

function clearReplies(server)
  if (isfield(server, 'slave_sock'))
    for i=1:length(server.slave_sock)
      server.slave_sock{i}.setSoTimeout(200);
      try
        server.getReply(i);
      end
      server.slave_sock{i}.setSoTimeout(0);
    end
  end
end

function delete(server)
% SERVER.DELETE destroys a server instance.
%
% This closes the listening socket, all slave sockets, and terminates any
% launched slave processes.

  fprintf(1, 'Destroying server on %s:%d.\n', ...
          server.hostname, server.options.port);

  if (isfield(server, 'serv_sock'))
    server.serv_sock.close();
  end

  if (isfield(server, 'slave_sock'))
    for i=1:length(server.slave_sock)
      server.slave_sock{i}.close();
    end
  end

  if (isfield(server, 'slave_process'))
    display('Shutting down slaves...');
    for i=1:length(server.slave_process)
      server.slave_process{i}.destroy();
    end
  end

end

end % end methods

end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function varargout = transpose(varargin)
% .'   SPHEREFUN transpose.
%    F.' is the non-conjugate transpose of a F.
%
% See also CTRANSPOSE.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

[varargout{1:nargout}] = transpose@separableApprox(varargin{:});

end


text:
function g = mlpderiv(net, x)
%MLPDERIV Evaluate derivatives of network outputs with respect to weights.
%
%	Description
%	G = MLPDERIV(NET, X) takes a network data structure NET and a matrix
%	of input vectors X and returns a three-index matrix G whose I, J, K
%	element contains the derivative of network output K with respect to
%	weight or bias parameter J for input pattern I. The ordering of the
%	weight and bias parameters is defined by MLPUNPAK.
%
%	See also
%	MLP, MLPPAK, MLPGRAD, MLPBKP
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Check arguments for consistency
errstring = consist(net, 'mlp', x);
if ~isempty(errstring);
  error(errstring);
end

[y, z] = mlpfwd(net, x);

ndata = size(x, 1);

if isfield(net, 'mask')
  nwts = size(find(net.mask), 1);
  temp = zeros(1, net.nwts);
else
  nwts = net.nwts;
end

g = zeros(ndata, nwts, net.nout);
for k = 1 : net.nout
  delta = zeros(1, net.nout);
  delta(1, k) = 1;
  for n = 1 : ndata
    if isfield(net, 'mask')
      temp = mlpbkp(net, x(n, :), z(n, :), delta);
      g(n, :, k) = temp(logical(net.mask));
    else
      g(n, :, k) = mlpbkp(net, x(n, :), z(n, :),...
	delta);
    end
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [ y, m ] = month_carry_gregorian ( y, m )

%*****************************************************************************80
%
%% MONTH_CARRY_GREGORIAN carries a year of months on the Gregorian calendar.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, integer M, the year and month.
%
%    Output, integer Y, integer M, the year and month.
%    On output, M is no greater than 12.
%
  while ( 1 )

    months = year_length_months_gregorian ( y );

    if ( m <= months )
      break
    end

    m = m - months;
    y = y + 1;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%Create A figure
figure
mainAxes = axes;
plot(mainAxes ,1:5,1:5,'r',1:5,1:0.5:3,'b');
xlabel 'This is an X Axis Label'
ylabel 'This is a Y Axis Label'
title('A Generic Title');
legend('Red Line','Blue Line');
set(gca,'Layer','Top')

%Break The Axes
h = breakxaxis([2 3]);

%Un-Break The X Axes
%unbreakxaxes(h)

author: Sable

repo: mcbench-benchmarks


text:
function re=rev(x);
% REV
% REV(x) reverses the elements of x
d=size(x); li=d(1,1);
ind=li:-1:1;
re=x(ind,:);


text:
classdef nme_line3p < mp.nm_element & mp.form_acp

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

%     properties
%     end

    methods
        function name = name(obj)
            name = 'line3p';
        end

        function np = np(obj)
            np = 6;     %% this is a 6 port element
        end

        function obj = build_params(obj, nm, dm)
            build_params@mp.nm_element(obj, nm, dm);    %% call parent

            dme = obj.data_model_element(dm);
            bus_dme = dm.elements.(dme.cxn_type);
            nk = 3*obj.nk;

            base_kv = bus_dme.tab.base_kv(bus_dme.on(dme.fbus(dme.on))) / sqrt(3);
            base_z = 1000 / dm.base_kva * base_kv .^ 2;

            sf = (base_z ./ dme.len) * ones(1, 6);      %% scale factor
            y_od = dme.ys(dme.lc, :) .* sf;             %% - of off-diagonal
            y_d  = y_od + dme.yc(dme.lc, :)/2 .* sf;    %% diagonal

            Y_od = obj.vec2symmat_stacked(y_od);        %% - of off-diagonal
            Y_d  = obj.vec2symmat_stacked(y_d);         %% diagonal

            obj.Y = [ Y_d   -Y_od;
                     -Y_od   Y_d ];
        end

        function M = vec2symmat_stacked(obj, vv)
            %% analogous to making a symmetric matrix from a vector of 6 values
            %% M = [v(1) v(2) v(3);
            %%      v(2) v(4) v(5);
            %%      v(3) v(5) v(6) ];
            %% except v(k) is now a diagonal matrix build from vector vv(:, k)
            n = size(vv, 1);
            v = cell(6, 1);
            for k = 1:6
                v{k} = spdiags(vv(:, k), 0, n, n);
            end
            M = [v{1} v{2} v{3};
                 v{2} v{4} v{5};
                 v{3} v{5} v{6} ];
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function y = spdiag(x)
y = sparse(1:numel(x),1:numel(x),x);
end

author: jacobeisenstein

repo: SAGE


text:
function out = interp3wrapper(x,y,z,v,xi,yi,zi,method,exterpval,limitelem)
%
% out = interp3wrapper(x,y,z,v,xi,yi,zi,method,exterpval,limitelem)
%
dim = size(xi);
t = prod(dim);

if ndims(x) == 3
	% We don't need x y z to be a whole 3D matrix
	x = squeeze(x(1,:,1));
	y = squeeze(y(:,1,1));
	z = squeeze(z(1,1,:));
end

if( ~exist('limitelem','var') || isempty(limitelem) )
	limitelem = 200*200*50;
end

if t <= limitelem || ndims(xi) == 2
	if exist('exterpval','var') && ~isempty(exterpval)
		out = interp3(x,y,z,v,xi,yi,zi,method,exterpval);
	else
		out = interp3(x,y,z,v,xi,yi,zi,method);
	end
	return;
end

N = ceil( t / limitelem );
spacing = max(floor(dim(3)/N),2);
N = ceil(dim(3)/spacing);

%fprintf('Interp3wrapper ');
for k = 1:N
	%fprintf('.');
	zmin = (k-1)*spacing+1;
	zmax = min(dim(3),k*spacing);

	zmin2 = min(min(min(zi(:,:,zmin:zmax))));
	zmax2 = max(max(max(zi(:,:,zmin:zmax))));
	zminidx = find(z<zmin2,1,'last');
	if isempty(zminidx)
		zminidx = 1;
	end
	zmaxidx = find(z>zmax2,1,'first');
	if isempty(zmaxidx)
		zmaxidx = length(z);
	end
	if zmaxidx == zminidx
		if zmaxidx >= length(z)
			zminidx = zminidx-1;
		else
			zmaxidx = zmaxidx+1;
		end
	end

	xmin2 = min(min(min(xi(:,:,zmin:zmax))));
	xmax2 = max(max(max(xi(:,:,zmin:zmax))));
	xminidx = find(x<xmin2,1,'last');
	if isempty(xminidx)
		xminidx = 1;
	end
	xmaxidx = find(x>xmax2,1,'first');
	if isempty(xmaxidx)
		xmaxidx = length(x);
	end

	ymin2 = min(min(min(yi(:,:,zmin:zmax))));
	ymax2 = max(max(max(yi(:,:,zmin:zmax))));
	yminidx = find(y<ymin2,1,'last');
	if isempty(yminidx)
		yminidx = 1;
	end
	ymaxidx = find(y>ymax2,1,'first');
	if isempty(ymaxidx)
		ymaxidx = length(y);
	end

	if exist('exterpval','var') && ~isempty(exterpval)
		out(:,:,zmin:zmax) = interp3(x(xminidx:xmaxidx),y(yminidx:ymaxidx),z(zminidx:zmaxidx),v(yminidx:ymaxidx,xminidx:xmaxidx,zminidx:zmaxidx),xi(:,:,zmin:zmax),yi(:,:,zmin:zmax),zi(:,:,zmin:zmax),method,exterpval);
	else
		out(:,:,zmin:zmax) = interp3(x(xminidx:xmaxidx),y(yminidx:ymaxidx),z(zminidx:zmaxidx),v(yminidx:ymaxidx,xminidx:xmaxidx,zminidx:zmaxidx),xi(:,:,zmin:zmax),yi(:,:,zmin:zmax),zi(:,:,zmin:zmax),method);
	end
end
%fprintf('\n');

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function point_num = sparse_grid_cc_size ( dim_num, level_max )

%*****************************************************************************80
%
%% SPARSE_GRID_CC_SIZE sizes a sparse grid using Clenshaw Curtis rules.
%
%  Discussion:
%
%    The grid is defined as the sum of the product rules whose LEVEL
%    satisfies:
%
%      0 <= LEVEL <= LEVEL_MAX.
%
%    This calculation is much faster than a previous method.  It simply
%    computes the number of new points that are added at each level in the
%    1D rule, and then counts the new points at a given DIM_NUM dimensional
%    level vector as the product of the new points added in each dimension.
%
%    This approach will work for nested families, and may be extensible
%    to other families, and to mixed rules.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 December 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Output, integer POINT_NUM, the total number of unique
%    points in the grids.
%

%
%  Special case.
%
  if ( level_max < 0 )
    point_num = 0;
    return
  end

  if ( level_max == 0 )
    point_num = 1;
    return
  end
%
%  Construct the vector that counts the new points in the 1D rule.
%
  new_1d = zeros ( level_max+1, 1 );

  new_1d(0+1) = 1;
  new_1d(1+1) = 2;

  j = 1;
  for l = 2 : level_max
    j = j * 2;
    new_1d(l+1) = j;
  end
%
%  Count the number of points by counting the number of new points
%  associated with each level vector.
%
  level_1d = zeros ( dim_num, 1 );

  point_num = 0;

  for level = 0 : level_max

    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );

      point_num = point_num + prod ( new_1d(level_1d(1:dim_num)+1) );

      if ( ~more )
        break
      end

    end

  end

  return
end


text:
function Draw_MPC_point_stabilization_v1 (t,xx,xx1,u_cl,xs,N,rob_diam)

set(0,'DefaultAxesFontName', 'Times New Roman')
set(0,'DefaultAxesFontSize', 12)

line_width = 1.5;
fontsize_labels = 14;

%--------------------------------------------------------------------------
%-----------------------Simulate robots -----------------------------------
%--------------------------------------------------------------------------
x_r_1 = [];
y_r_1 = [];

r = rob_diam/2;  % obstacle radius
ang=0:0.005:2*pi;
xp=r*cos(ang);
yp=r*sin(ang);

figure(500)
% Animate the robot motion
%figure;%('Position',[200 200 1280 720]);
set(gcf,'PaperPositionMode','auto')
set(gcf, 'Color', 'w');
set(gcf,'Units','normalized','OuterPosition',[0 0 0.55 1]);

tic
for k = 1:size(xx,2)
    h_t = 0.14; w_t=0.09; % triangle parameters

    x1 = xs(1); y1 = xs(2); th1 = xs(3);
    x1_tri = [ x1+h_t*cos(th1), x1+(w_t/2)*cos((pi/2)-th1), x1-(w_t/2)*cos((pi/2)-th1)];%,x1+(h_t/3)*cos(th1)];
    y1_tri = [ y1+h_t*sin(th1), y1-(w_t/2)*sin((pi/2)-th1), y1+(w_t/2)*sin((pi/2)-th1)];%,y1+(h_t/3)*sin(th1)];
    fill(x1_tri, y1_tri, 'g'); % plot reference state
    hold on;
    x1 = xx(1,k,1); y1 = xx(2,k,1); th1 = xx(3,k,1);
    x_r_1 = [x_r_1 x1];
    y_r_1 = [y_r_1 y1];
    x1_tri = [ x1+h_t*cos(th1), x1+(w_t/2)*cos((pi/2)-th1), x1-(w_t/2)*cos((pi/2)-th1)];%,x1+(h_t/3)*cos(th1)];
    y1_tri = [ y1+h_t*sin(th1), y1-(w_t/2)*sin((pi/2)-th1), y1+(w_t/2)*sin((pi/2)-th1)];%,y1+(h_t/3)*sin(th1)];

    plot(x_r_1,y_r_1,'-r','linewidth',line_width);hold on % plot exhibited trajectory
    if k < size(xx,2) % plot prediction
        plot(xx1(1:N,1,k),xx1(1:N,2,k),'r--*')
    end

    fill(x1_tri, y1_tri, 'r'); % plot robot position
    plot(x1+xp,y1+yp,'--r'); % plot robot circle

    hold off
    %figure(500)
    ylabel('$y$-position (m)','interpreter','latex','FontSize',fontsize_labels)
    xlabel('$x$-position (m)','interpreter','latex','FontSize',fontsize_labels)
    axis([-0.2 1.8 -0.2 1.8])
    pause(0.0)
    box on;
    grid on
    %aviobj = addframe(aviobj,gcf);
    drawnow
    % for video generation
    F(k) = getframe(gcf); % to get the current frame
end
toc
close(gcf)
%viobj = close(aviobj)
%video = VideoWriter('exp.avi','Uncompressed AVI');

% video = VideoWriter('exp.avi','Motion JPEG AVI');
% video.FrameRate = 5;  % (frames per second) this number depends on the sampling time and the number of frames you have
% open(video)
% writeVideo(video,F)
% close (video)

figure
subplot(211)
stairs(t,u_cl(:,1),'k','linewidth',1.5); axis([0 t(end) -0.35 0.75])
ylabel('v (rad/s)')
grid on
subplot(212)
stairs(t,u_cl(:,2),'r','linewidth',1.5); axis([0 t(end) -0.85 0.85])
xlabel('time (seconds)')
ylabel('\omega (rad/s)')
grid on

author: MMehrez


text:
%--- help for create_variable_names ---
%
%  INTERNAL FUNCTION
%
%

author: jmaih


text:
function plotTrajectories(data, IDstart, IDend)
%
% Authors:       Alexander Wischnewski
%
% Description:
%   used to visualize the timely behavior of new target trajectories
%
% Inputs:
%   IDstart:    Trajectory ID start (see debug_slow.debug_slow_tartraj_TrajCnt)
%   IDend:      Trajectory ID end (see debug_slow.debug_slow_tartraj_TrajCnt)

    % colors
    C = {'k','b','r','g','y',[.5 .6 .7],[.8 .2 .6],[.3 .1 .2],[.1 .5 .9]}; % Cell array of colros.
    figure;
    xy_fig = gca;
    hold on; grid on;
    xlabel('x East in m');
    ylabel('y North in m');
    axis equal;
    figure;
    kappa_fig = gca;
    hold on; grid on;
    xlabel('Global s-coordinate in m');
    ylabel('Kappa in radpm');
    figure;
    v_fig = gca;
    hold on; grid on;
    xlabel('Global s-coordinate in m');
    ylabel('Velocity in mps');
    for i = IDstart:1:IDend
        % find idx with trajectory
        idx = find(data.debug_slow.debug_slow_tartraj_TrajCnt.Data==i, 1, 'first');
        disp(idx);
        % map color index to available colors
        color_idx = mod(i-IDstart, length(C)-1) + 1;
        plot(xy_fig, data.debug_slow.debug_slow_tartraj_x_m.Data(idx, :), ...
            data.debug_slow.debug_slow_tartraj_y_m.Data(idx, :), 'color', C{color_idx});
        scatter(xy_fig, data.debug_slow.debug_slow_x_real_m.Data(idx), ...
            data.debug_slow.debug_slow_y_real_m.Data(idx), 30, 'filled', ...
            'MarkerEdgeColor', C{color_idx}, 'MarkerFaceColor', C{color_idx});
        plot(kappa_fig, data.debug_slow.debug_slow_tartraj_s_glob_m.Data(idx, :), ...
            data.debug_slow.debug_slow_tartraj_kappa_radpm.Data(idx, :));
        plot(v_fig, data.debug_slow.debug_slow_tartraj_s_glob_m.Data(idx, :), ...
            data.debug_slow.debug_slow_tartraj_v_mps.Data(idx, :));
    end
end

author: TUMFTM


text:
function hole_point = p06_hole_point ( hole_index, m )

%*****************************************************************************80
%
%% P06_HOLE_POINT returns a point inside a given hole in problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer HOLE_INDEX, the index of the hole.
%
%    Input, integer M, the spatial dimension.
%
%    Output, real HOLE_POINT(M), a point in the hole
%
  hole_point(1:2) = [ 0.0, 0.0 ];

  return
end


text:
function varargout = trigonometric(X)
% Trigonometric function
%
%   TRIGONOMETRIC([x1, x2, ..., xn]) returns the value of the
%   trigonometric fnuction at the specified points. All [xi] may be
%   vectors.
%
%   The global minimum is
%
%               f(x1, x2, ..., xn) = f(0, 0, ..., 0) = 0.

% Author: Rody P.S. Oldenhuis
% Delft University of Technology
% E-mail: oldenhuis@dds.nl
% Last edited 28/Feb/2009

    % if no input is given, return dimensions, bounds and minimum
    if (nargin == 0)
        varargout{1} = inf;  % # dims
        varargout{2} = -inf; % LB
        varargout{3} = +inf; % LB
        varargout{4} = 0; % solution
        varargout{5} = 0; % function value at solution

    % otherwise, output function value
    else

        % create basevectors
        n = size(X, 2);
        i = repmat(1:n, size(X, 1), 1);

        % compute cosine only once
        cosX = cos(X);

        % compute intermediate sum
        sumcosX = repmat(sum(cosX, 2), 1, n);

        % output rowsum
        varargout{1} = sum( (n + i.*(1 - cosX) - sin(X) - sumcosX).^2, 2);

    end

end


text:
function [v,r] = dquat2line(dq)

% DQUAT2LINE Transforms a line position dual quaternion into
%            its vector representation (orientation + position).
%
%     [V,R] = DQUAT2LINE(DQ) transforms the double quaternion
%         representation of a line position (in the euclidean space) into a
%         vector V which represents the line orientation (unitary vector)
%         and a vector R which is the coordinates of the nearest point from
%         the origin belonging to the line (it could be another r, it is
%         NOT unique).
%         DQ is either a vector of size 8 or an array of size 8*N (each
%         column represents a line position dual quaternion) where N is the
%         number of lines. V and R have the same size. They are either a
%         vector of size 3 or an array of size 3*N depending on the input
%         format.
%
% See also DQUAT2POS, DQUAT2VEL, DQUAT2LINEVEL, LINE2DQUAT

if nargout < 2 || nargout > 2
    error('DualQuaternion:dquat2line:wrongNumberOfOutputs',...
        'You specified %d outputs. It should be 2 (see documentation).',nargout);
end

sdq = size(dq);
if sdq == [1 8], dq = dq'; sdq = size(dq); end

% wrong format
if sdq(1) ~= 8
    error('DualQuaternion:dquat2line:wrongsize',...
        '%d rows in the DQ array. It should be 8.',sdq(1));
end

% check line position dual quaternion
tol = 1e-6;
[maxval,imax] = max(max(abs(dq([1 5],:))));
if maxval > tol
   warning('DualQuaternion:dquat2line:wrongFormat',...
        'At least one dual quaternion is not a line position dual quaternion (tol = %.1e).\n Indices of max values: %d \n Max value = %.2e',...
        tol,imax,maxval);
end

% extraction of the line position parameters, v and r
v = dq(2:4,:);
w = dq(6:8,:);

% check if the norm of v is unitary (it should be)
normv2 = sum(v.^2);
[maxval2,imax2] = max(normv2);
[minval2,imin2] = min(normv2);
i2 = union(imax2,imin2);
if maxval2 > 1+tol || minval2 < 1-tol
       warning('DualQuaternion:dquat2line:wrongFormatNotUnitary',...
        'At least one dual quaternion has not a unitary orientation ( is it a line position dual quaternion?) (tol = %.1e).\n Indices of max values: %d \n Min and Max value = [%.2e %.2e]',...
        tol,i2,minval2,maxval2);
end
% r is the position vector of the point belonging to the line such that
% this point is the nearest point from the origin (it is orthogonal to the
% line orientation)
normv2 = repmat(normv2,3,1);
r = cross(v,w)./normv2;


text:
% eeg_mergechan() - merge channel structure while preserving channel
%                   order
%
% >> mergelocs = eeg_mergechan(locs1, locs2);
%
% Inputs:
%     locs1     - EEGLAB channel location structure
%     locs2     - second EEGLAB channel location structure
%
% Output:
%     mergelocs - merged channel location structure
%
% Author: Arnaud Delorme, August 2006

% Copyright (C) Arnaud Delorme, CERCO, 2006, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% union of two channel location structure
% without loosing the order information
% ---------------------------------------
function alllocs = myunion(locs1, locs2)

   labs1 = { locs1.labels };
   labs2 = { locs2.labels };

   count1 = 1;
   count2 = 1;
   count3 = 1;
   alllocs = locs1; alllocs(:) = [];
   while count1 <= length(locs1) | count2 <= length(locs2)

       if count1 > length(locs1)
           alllocs(count3) = locs2(count2);
           count2 = count2 + 1;
           count3 = count3 + 1;
       elseif count2 > length(locs2)
           alllocs(count3) = locs1(count1);
           count1 = count1 + 1;
           count3 = count3 + 1;
       elseif strcmpi(labs1(count1), labs2(count2))
           alllocs(count3) = locs1(count1);
           count1 = count1 + 1;
           count2 = count2 + 1;
           count3 = count3 + 1;
       elseif isempty(strmatch(labs1(count1), labs2))
           alllocs(count3) = locs1(count1);
           count1 = count1 + 1;
           count3 = count3 + 1;
       else
           alllocs(count3) = locs2(count2);
           count2 = count2 + 1;
           count3 = count3 + 1;
       end;

   end;

author: PatternRecognition

repo: OpenBMI


text:
function pars = build_parameters_structure_v4(R, opts, varargin)

%BUILD_PARAMETERS_STRUCTURE Build parameters structure used for function ibTB.

%   Copyright (C) 2009 Cesare Magri
%   Version 4.0.0

% -------
% LICENSE
% -------
% This software is distributed free under the condition that:
%
% 1. it shall not be incorporated in software that is subsequently sold;
%
% 2. the authorship of the software shall be acknowledged and the following
%    article shall be properly cited in any publication that uses results
%    generated by the software:
%
%      Magri C, Whittingstall K, Singh V, Logothetis NK, Panzeri S: A
%      toolbox for the fast information analysis of multiple-site LFP, EEG
%      and spike train recordings. BMC Neuroscience 2009 10(1):81;
%
% 3.  this notice shall remain in place in each source file.

if isempty(varargin)
    msg = 'No output option specified';
    error('buildParsStruct:noOutputOptSpecified', msg);
end

% R -----------------------------------------------------------------------
if ndims(R)==3
    [pars.Nc, size2ofR, pars.Ns] = size(R);
else
    msg = 'Response matrix must be a 3D matrix.';
    error('buildParsStruct:Rnot3D', msg);
end;

% NT ----------------------------------------------------------------------
% If the number of trials per stimulus is constant, NT can be provided
% as a scalar. In this case, the array version is built internally.
if isscalar(opts.nt)
    pars.Nt = ones(pars.Ns,1) * opts.nt;
    maxNt = pars.Nt;
else
    if length(opts.nt)~=pars.Ns
        msg = 'size(R,3) must match length(opts.nt). Try transposing nt.';
        error('buildParsStruct:lengthNt_vs_size3R', msg);
    end

    % R and NT compatibility:
    pars.Nt = opts.nt(:);
    maxNt = max(pars.Nt);
end

% R and NT compatibility:
if maxNt~=size2ofR
    msg = 'max(nt) must be equal to size(R,2).';
    error('buildParsStruct:maxNt_vs_size2R', msg);
end

% METHOD ------------------------------------------------------------------
pars.method = opts.method;
switch lower(opts.method)
    case {'dr'}
        pars.methodFunc = @direct_method_v5b;
        pars.methodNum = 1;
    case {'gs'}
        pars.methodFunc = @gaussian_method_v7_1_0;
        pars.methodNum = 2;
    otherwise
        msg = ['Undefined method ' pars.method];
        error('buildParsStruct:methodNotFound', msg);
end;

% BIAS --------------------------------------------------------------------
switch lower(opts.bias)
    case 'naive'
        pars.biasCorrNum = 0;
    case 'qe'
        pars.biasCorrNum = 1;
    case 'pt'
        pars.biasCorrNum = 2;
    case 'gsb'
        pars.biasCorrNum = 3;
    otherwise
        msg = ['Bias correction option ''' opts.bias ''' not found'];
        error('buildParsStruct:biasCorrNotFound', msg);
end

% BTSP (optional) ---------------------------------------------------------
pars.numberOfSpecifiedOptions = 0;
if isfield(opts, 'btsp')
    pars.numberOfSpecifiedOptions = pars.numberOfSpecifiedOptions + 1;

    if round(opts.btsp)~=opts.btsp
        msg = 'opts.btsp must be an integer.';
        error('buildParsStruct:btspNotInteger', msg);
    else
        pars.btsp = opts.btsp;
    end
else
    % No bootstrap:
    pars.btsp = 0;
end

% OUTPUT LIST -------------------------------------------------------------
% Checking which output-options have been selected and keeping track of the
% position they have in VARARGIN: this will allow do provide the outputs in
% the correct order.
pars.Noutput = length(varargin); % number of ouputs

varargin = lower(varargin);
pars.HR    = strcmpi(varargin, 'hr'      );
pars.HRS   = strcmpi(varargin, 'hrs'     );
pars.HlR   = strcmpi(varargin, 'hlr'     );
pars.HlRS  = strcmpi(varargin, 'hirs'    );
pars.HiR   = strcmpi(varargin, 'hir'     );
pars.ChiR  = strcmpi(varargin, 'chir'    );
pars.HshR  = strcmpi(varargin, 'hshr'    );
pars.HshRS = strcmpi(varargin, 'hshrs'   );
% See note in ENTROPY regarding this quantity:
pars.HiRS  = strcmpi(varargin, 'hirsdef' );

if pars.Nc>1
    pars.doHR    = any(pars.HR   );
    pars.doHRS   = any(pars.HRS  );
    pars.doHlR   = any(pars.HlR  );
    pars.doHlRS  = any(pars.HlRS );
    pars.doHiR   = any(pars.HiR  );
    pars.doHiRS  = any(pars.HiRS );
    pars.doChiR  = any(pars.ChiR );
    pars.doHshR  = any(pars.HshR );
    pars.doHshRS = any(pars.HshRS);
else
    pars.doHR    = any(pars.HR | pars.HlR | pars.HiR | pars.ChiR | pars.HshR);
    pars.doHRS   = any(pars.HRS | pars.HlRS | pars.HiRS | pars.HshRS);
    pars.doHlR   = false;
    pars.doHlRS  = false;
    pars.doHiR   = false;
    pars.doHiRS  = false;
    pars.doChiR  = false;
    pars.doHshR  = false;
    pars.doHshRS = false;
end

% ADDCHECKS (optionals) ---------------------------------------------------
if isfield(opts, 'verbose') && opts.verbose
    pars.numberOfSpecifiedOptions = pars.numberOfSpecifiedOptions + 1;
    pars.lengthVarargin = length(varargin);

    pars.addChecks = true;
else
    pars.addChecks = false;
end

if any(pars.addChecks)
	additional_checks_v4(R, pars, opts);
end;

author: fieldtrip

repo: fieldtrip