text:
% EEG_EPOCHFORMAT - Convert the epoch information of a dataset from struct
%                     to array or vice versa.
%
% Usage: >> [epochsout fields] = eeg_epochformat( epochs, 'format', fields, events );
%
% Input:
%   epochs - epoch numerical or cell array or epoch structure
%   format - ['struct'|'array'] convert epoch array to structure and epoch
%            structure to array.
%   fields - [optional] cell array of strings containing the names of
%            the epoch struct fields. If this field is empty, it uses the
%            following list for the names of the fields { 'var1' 'var2' ... }.
%            For structure conversion, this field helps export a given
%            event type. If this field is left empty, the time locking
%            event for each epoch is exported.
%   events - numerical array of event indices associated with each epoch.
%            For array conversion, this field is ignored.
%
% Outputs:
%   epochsout - output epoch array or structure
%   fields    - output cell array with the name of the fields
%
% Epoch format:
%   struct - Epoch information is organised as an array of structs
%   array  - Epoch information is organised as an 2-d array of numbers,
%            each column representing a user-defined variable (the
%            order of the variable is a function of its order in the
%            struct format).
%
% Note: 1) The epoch structure is defined only for epoched data.
%       2) The epoch 'struct' format is more comprehensible.
%          For instance, to see all the properties of epoch i,
%          type >> EEG.epoch(i)
%          Unfortunately, structures are awkward for expert users to deal
%          with from the command line (Ex: To get an array of 'var1' values,
%           >> celltomat({EEG.epoch(:).var1})')
%          In array format, assuming 'var1' is the first variable
%          declared, the same information is obtained by
%           >> EEG.epoch(:,1)
%       3) This function automatically updates the 'epochfields'
%          cell array depending on the format.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 12 Feb 2002
%
% See also: EEGLAB

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) Arnaud Delorme, CNL / Salk Institute, 12 Feb 2002, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% $Log: eeg_epochformat.m,v $
% Revision 1.4  2005/05/24 16:57:09  arno
% cell2mat
%
% Revision 1.3  2003/07/20 19:32:20  scott
% typos
%
% Revision 1.2  2002/04/21 01:10:35  scott
% *** empty log message ***
%
% Revision 1.1  2002/04/05 17:32:13  jorn
% Initial revision
%

% 03/13/02 added field arrays options -ad

function [ epoch, fields, epocheventout] = eeg_epochformat( epoch, format, fields, epochevent);

if nargin < 2
   help eeg_epochformat;
   return;
end;

if nargin < 3
    fields = {};
end
epocheventout = [];

switch format
case 'struct'
   if ~isempty(epoch) && ~isstruct(epoch)

      fields = getnewfields( fields, size(epoch,2) - length(fields));

      % generate the structure
      % ----------------------
      command = 'epoch = struct(';
      for index = 1:length(fields)
          if iscell(epoch)
              command = [ command '''' fields{index} ''', epoch(:,' num2str(index) ')'',' ];
          else
              command = [ command '''' fields{index} ''', mattocell( epoch(:,' num2str(index) ')'',' ...
                            '[1], ones(1,size(epoch,1))),' ];
          end;
      end
      eval( [command(1:end-1) ');' ] );

      if exist('epochevent') == 1
         for index = 1:size(epoch,2)
             if iscell(epochevent)
                 epoch(index).event = epochevent{index};
             else
                 epoch(index).event = epochevent(index);
             end
         end
      end
   end

case 'array'
    if isstruct(epoch)

        % note that the STUDY std_maketrialinfo also gets the epoch info for the
        % time locking event

        selectedType = fields;
        if iscell(fields) && ~isempty(fields), selectedType = fields{1}; end
  	    fields = fieldnames( epoch );

        eval( [ 'values = { epoch.' fields{1} ' };' ]);

        if any(cellfun(@length, values) > 1)
            if ~isfield(epoch, 'eventlatency')
                error('eventlatency field not present in data epochs');
            end

            if isempty(selectedType)
                % find indices of time locking events
                for index = 1:length(epoch)
                    epochlat = [ epoch(index).eventlatency{:} ];
                    tmpevent = find( abs(epochlat) < 0.02 );
                    if isempty(tmpevent)
                        error('time locking event missing, cannot convert to array');
                    end
                    epochSubIndex(index) = tmpevent;
                end
            else
                % find indices of specific event type (if several take the
                % first one
                for index = 1:length(epoch)
                    epochtype = epoch(index).eventtype;
                    tmpeventind = strmatch( selectedType, epochtype );
                    if length(tmpeventind) > 1
                        fprintf('Warning: epoch %d has several events of "type" %s, taking the fist one\n', index, selectedType);
                    end
                    if isempty(tmpeventind)
                         epochSubIndex(index) = NaN;
                    else epochSubIndex(index) = tmpeventind(1);
                    end
                end
            end
        else
            epochSubIndex = ones(1, length(epoch));
        end

        % copy values to array
        tmp = cell( length(epoch), length( fields ));
        for index = 1:length( fields )
            for trial = 1:length(epoch)
                tmpval = getfield(epoch, {trial}, fields{index});
                if isnan(epochSubIndex(trial))
                     tmp(trial, index) = { NaN };
                elseif iscell(tmpval)
                     tmp(trial, index) = tmpval(epochSubIndex(trial));
                elseif ~ischar(tmpval)
                     tmp(trial, index) = { tmpval(epochSubIndex(trial)) };
                else tmp(trial, index) = { tmpval };
                end
            end
        end
        epoch = tmp;
    end
    otherwise, error('unrecognised format');
end

return;

% create new field names
% ----------------------
function epochfield = getnewfields( epochfield, nbfields )
   count = 1;
   if nbfields > 0
        while nbfields > 0
            if isempty( strmatch([ 'var' int2str(count) ], epochfield ) )
               epochfield =  { epochfield{:} [ 'var' int2str(count) ] };
               nbfields = nbfields-1;
            else    count = count+1;
            end
        end;
   else
        epochfield = epochfield(1:end+nbfields);
   end;
return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [v,t]=histndim(x,b,mode)
%HISTNDIM - generates and/or plots an n-dimensional histogram
%
%  Inputs:  X(m,d)   is the input data: each row is one d-dimensiona data point
%           B(3,d)   specifies the histogram bins.
%                         B(1,:) gives the number of bins in each dmension [default 10]
%                         B(2,:) gives the minimum of the first bin in each dimension [default min(X)]
%                         B(3,:) gives the maximum of the last bin in each dimension [default max(X)]
%                    If B has only one column, the same values are use for al dimensions
%                    If B(1,i)=0 then that dimension will be ignored (and excluded from V)
%           MODE     is a character string containing a combination of the following:
%                        'z' for zero base in the 2D plot [default base = min(V)]
%                        'p' to scale V as probabilities [default actual counts]
%                        'h' to plot a histogam even if output arguments are present
%
% Outputs:  V        d-dimensional array containing the histogram counts
%           T        d-element cell array. d{i} contains the bin boundary values for
%                    the i'th dimension. The length of d{i} is one more than the number of bins
%                    in that dimension.
%
%                    Note that if any of B(1,:) are zero then the number of dimensions in V and elements
%                    of T will be correspondingly reduced.
%
% Example: histndim(randn(100000,2),[20 -3 3]','pz');

%	   Copyright (C) Mike Brookes 2004
%      Version: $Id: histndim.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n,d]=size(x);
if nargin<3
    mode=' ';
    if(nargin<2)
        b=repmat(10,1,d);
    end
end

if size(b,2)==1
    b=repmat(b,1,d);
end
if size(b,1)<3
    mi=min(x,[],1);
    ma=max(x,[],1);
    w=(ma-mi)./(b(1,:)-0.001);  % nudge slightly to make sure al points included
    b(3,:)=ma+0.0005*w;
    b(2,:)=mi-0.0005*w;
end

acd=find(b(1,:)>0);
sv=b(1,acd);
nbt=prod(sv);
t=cell(length(acd),1);

% loop through each dimension
k=1;        % indexing factor
ok=repmat(1>0,n,1);
ix=repmat(nbt-sum(cumprod(sv)),n,1);
for i=1:length(acd)
    j=acd(i);
    bw=b(1,j)/(b(3,j)-b(2,j));
    bi=ceil((x(:,j)-b(2,j))*bw);
    ok=ok & (bi>0) & (bi<=b(1,j));
    ix(ok)=ix(ok)+k*bi(ok);
    k=k*b(1,j);
    t{i}=b(2,j)+(0:b(1,j))/bw;
end
v=full(sparse(ix(ok),1,1,nbt,1));
if length(sv)>1
    v=reshape(v,sv);
end
if any(mode=='p')
    v=v/n;
end

if ~nargout | any(mode=='h')
    svg=find(sv>1);
    if length(svg)==1
        j=acd(svg);
        bar(b(2,j)+(0.5:sv(svg)-0.5)*(b(3,j)-b(2,j))/b(1,j),v(:));
    elseif length(svg)==2
        j=acd(svg(1));
        k=acd(svg(2));
        bj=b(1,j);
        bk=b(1,k);
        %     imagesc(b(2:3,k),b(2:3,j),reshape(v,b(1,j),b(1,k)));
        vda=kron(reshape(v,bj,bk),[1 1;1 1]);
        if any(mode=='z')
            ba=0;
        else
            ba=min(vda(:));
        end
        vda=[repmat(ba,1,2*bk+2);repmat(ba,2*bj,1) vda repmat(ba,2*bj,1);repmat(ba,1,2*bk+2)];
        jda=kron(t{svg(1)},[1 1]);
        jda=jda-(jda(3)-jda(2))*0.01*[-0.5 (-1).^(1:2*bj) 0.5]; % nudge slightly to avoid MATLAB plotting bug
        kda=kron(t{svg(2)},[1 1]);
        kda=kda-(kda(3)-kda(2))*0.01*[-0.5 (-1).^(1:2*bk) 0.5];
        surf(jda,kda,vda');
        ylabel(sprintf('Axis %d',k));
        xlabel(sprintf('Axis %d',j));
        colorbar;
    else
        fprintf(2,'Error in %s: Cannot plot 3-D histogram\n',mfilename);
    end
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [fieldfeatures]=segfieldsfeat(imsegs,nsegments,smap,xspfield,spfind,spdata,features)
sinds=[1:nsegments];
xspfield(isnan(xspfield))=0;

if exist('spfind','var')
    spfind1{1}=spfind{1};%floor
    spfind1{2}=[spfind{2};spfind{3};spfind{4};spfind{5}];% walls+ceiling
    clear spfind;
    spfind=spfind1;
end

%get model for fields
if exist('spfind','var')
    bgmean=[];
    bgmedian=[];
    for f =1:numel(spfind)
        if numel(spfind1{f})==0
            bgmean(f,1:21)=0;
            bgmedian(f,1:21)=0;
            continue;
        end
        sparea = imsegs.npixels(spfind{f});
        npix = sum(sparea);
        spnorm = sparea / npix;

        bgmean(f,1:3) = sum( spdata(spfind{f}, 1:3).*repmat(spnorm, [1 3]), 1);
        bgmean(f,4:6) = rgb2hsv(bgmean(f, 1:3));
        bgmean(f,7:21)=sum(spdata(spfind{f}, 15:29) .* repmat(spnorm, [1 15]), 1);

        featinds = [1:3];
        for i=1:length(featinds)
            [vv,ii] = sort(spdata(spfind{f},featinds(i)));
            temp = cumsum(spnorm(ii));
            tempfeat1(i) = vv(min(find(temp>0.5)));
        end
        bgmedian(f,1:3) = tempfeat1(:)';
        bgmedian(f,4:6) = rgb2hsv(bgmedian(f, 1:3));

        featinds = [15:29];
        for i=1:length(featinds)
            [vv,ii] = sort(spdata(spfind{f},featinds(i)));
            temp = cumsum(spnorm(ii));
            tempfeat2(i) = vv(min(find(temp>0.5)));
        end
        bgmedian(f,7:21) = tempfeat2(:)';

    end
end

for k = 1:nsegments

    spind = find(smap==sinds(k));
    sparea = imsegs.npixels(spind);
    npix = sum(sparea);

    oarea=sum(xspfield(spind,:),1);
    tfeatures = oarea(:)';
    tfeatures=tfeatures/(npix+(npix==0));
    tempfieldfeatures=tfeatures;
    tempfieldfeatures(tfeatures==0)=eps;

    entfieldfeat=-1*sum(tfeatures.*log2(tempfieldfeatures));
    tfeatures=[tfeatures entfieldfeat ];
    fieldfeatures(k,:)=tfeatures;

    nf=numel(fieldfeatures);
    if exist('spfind','var')

        for tempf=1:size(bgmedian,1)
            fieldfeatures(k,nf+(1:5)) = [abs(bgmedian(tempf,1:4)-features(k,1:4)) ...
                norm((features(k,15:29)- bgmedian(tempf,7:21)),2)] ;%r g b h and mean tex res diff
            nf=nf+5;
        end
    end

end

end

author: zouchuhang

repo: LayoutNet


text:
function [x, funVal, ValueL, res]=overlapping_LeastR(A, y, z, opts)
%
%%
% Function overlapping_LeastR
%      Least Squares Loss with the
%           overlapping group Lasso
%
%% Problem
%
%  min  1/2 || A x - y||^2 + z_1 \|x\|_1 + z_2 * sum_i w_i ||x_{G_i}||
%
%  G_i's are nodes
%
%    we have L1 for each element
%    and the L2 for the overlapping group
%
%  The overlapping group information is contained in
%
%   opts.G- a row vector containing the indices of all the overlapping
%           groups G_1, G_2, ..., G_groupNum
%
%   opts.ind- a 3 x groupNum matrix
%           opts.ind(1,i): the starting index of G_i in opts.G
%           opts.ind(2,i): the ending index of G_i in opts.G
%           opts.ind(3,i): the weight for the i-th group
%
% For better illustration, we consider the following example of four groups:
%   G_1={1,2,3}, G_2={2,4}, G_3={3,5}, G_4={1,5}.
% Let us assume the weight for each group is 123.
%
%   opts.G=[1,2,3,2,4,3,5,1,5];
%   opts.ind=[ [1, 3, 123]', [4,5,123]',[6,7,123]',[8,9,123]' ];
%
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size mx1)
%  z -        The regularization parameter (z=[z_1,z_2] >=0)
%  opts-      Optimal inputs (default value: opts=[])
%
%% Output parameters:
%
%  x-         Solution
%  funVal-    Function value during iterations
%
%% Copyright (C) 2010-2011 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on April 21, 2010.
%
%% Related papers
%
% [1]  Jun Liu and Jieping Ye, Fast Overlapping Group Lasso,
%      arXiv:1009.0306v1, 2010
%
%% Related functions:
%
%  sll_opts, initFactor, pathSolutionLeast
%
%%

%% Verify and initialize the parameters
%%
if (nargin <4)
    error('\n Inputs: A, y, z, and opts (.ind) should be specified!\n');
end

[m,n]=size(A);

if (length(y) ~=m)
    error('\n Check the length of y!\n');
end

if (z(1)<0 || z(2)<0)
    error('\n z should be nonnegative!\n');
end

lambda1=z(1);
lambda2=z(2);

opts=sll_opts(opts); % run sll_opts to set default values (flags)

% restart the program for better efficiency
%  this is a newly added function
if (~isfield(opts,'rStartNum'))
    opts.rStartNum=opts.maxIter;
else
    if (opts.rStartNum<=0)
        opts.rStartNum=opts.maxIter;
    end
end

%% Detailed initialization
%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Group & Others

% Initialize maxIter2
if (~isfield(opts,'maxIter2'))
    maxIter2=1000;
else
    maxIter2=opts.maxIter2;
end
% the maximal number of iteration for the projection

% Initialize tol2
if (~isfield(opts,'tol2'))
    tol2=1e-8;
else
    tol2=opts.tol2;
end
% the duality gap of the projection

% Initialize flag2
if (~isfield(opts,'flag2'))
    flag2=2;
else
    flag2=opts.flag2;
end

% Initialize G
if (~isfield(opts,'G'))
    error('\n In overlapping_LeastR, the field .G should be specified');
else
    G=opts.G-1;
    % we substract 1, as in C, the index begins with 0
end

% Initialize w
if (~isfield(opts,'ind'))
    error('\n In overlapping_LeastR, the field .ind should be specified');
else
    w=opts.ind;

    if (size(w,1)~=3)
        error('\n w is a 3 x groupNum matrix');
    end

    w(1:2,:)=w(1:2,:)-1;
    % we substract 1, as in C, the index begins with 0
end

groupNum=size(w,2);
% the number of groups

Y=zeros(length(G),1);
% the starting point for the projection

%% Starting point initialization

% compute AT y
if (opts.nFlag==0)
    ATy =A'*y;
elseif (opts.nFlag==1)
    ATy= A'*y - sum(y) * mu';  ATy=ATy./nu;
else
    invNu=y./nu;              ATy=A'*invNu-sum(invNu)*mu';
end

% process the regularization parameter
if (opts.rFlag~=0)
    % z here is the scaling factor lying in [0,1]
%     if (lambda1<0 || lambda1>1 || lambda2<0 || lambda2>1)
%         error('\n opts.rFlag=1, and z should be in [0,1]');
%     end

    % compute lambda1_max
    temp=abs(ATy);
    lambda1_max=max(temp);

    lambda1=lambda1*lambda1_max;

    % compute lambda2_max(lambda_1)

    % lambda_2_max to be added later

    lambda2_max=1;

    lambda2=lambda2*lambda2_max;
end

% initialize a starting point
if opts.init==2
    x=zeros(n,1);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if (length(x)~=n)
            error('\n Check the input .x0');
        end
    else
        x=ATy;  % if .x0 is not specified, we use ratio*ATy,
        % where ratio is a positive value
    end
end

% compute A x
if (opts.nFlag==0)
    Ax=A* x;
elseif (opts.nFlag==1)
    invNu=x./nu; mu_invNu=mu * invNu;
    Ax=A*invNu -repmat(mu_invNu, m, 1);
else
    Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
end

if (opts.init==0)
    % ------  This function is not available
    %
    % If .init=0, we set x=ratio*x by "initFactor"
    % Please refer to the function initFactor for detail
    %

    % Here, we only support starting from zero, due to the complex
    % structure

    %x=zeros(n,1);
end

%% The main program
% The Armijo Goldstein line search schemes + accelearted gradient descent

bFlag=0; % this flag tests whether the gradient step only changes a little

if (opts.mFlag==0 && opts.lFlag==0)
    L=1;
    % We assume that the maximum eigenvalue of A'A is over 1

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,1);

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

        % compute AT As
        if (opts.nFlag==0)
            ATAs=A'*As;
        elseif (opts.nFlag==1)
            ATAs=A'*As - sum(As) * mu';  ATAs=ATAs./nu;
        else
            invNu=As./nu;                ATAs=A'*invNu-sum(invNu)*mu';
        end

        % obtain the gradient g
        g=ATAs-ATy;

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;

        firstFlag=1;
        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the L1/Lq-norm regularized projection
            v=s-g/L;

            % projection
            [x,gap,penalty2]=overlapping(v,  n, groupNum, lambda1/L, lambda2/L,...
                w, G, Y, maxIter2, flag2, tol2);

            if (nargout ==4)
                % record the number of iterations

                if (firstFlag)
                    res.projStep(iterStep)=penalty2(4);
                else
                    res.projStep(iterStep)=...
                        max(res.projStep(iterStep),penalty2(4) );
                end
            end
            firstFlag=0;

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute A x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./nu; mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
            end

            Av=Ax -As;
            r_sum=v'*v; l_sum=Av'*Av;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is ||Av||_2^2 <= L * ||v||_2^2
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                %fprintf('\n L=%5.6f',L);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        xxp=x-xp;   Axy=Ax-y;

        ValueL(iterStep)=L;

        % function value = loss + regularizatioin
        funVal(iterStep)=Axy'* Axy/2 + lambda1 * sum(abs(x)) + lambda2 *penalty2 (1);

        if (nargout ==4)
            % record pp and gg;
            res.pp(iterStep)=penalty2 (2);
            res.qq(iterStep)=penalty2 (3);

            % record gap
            res.gap(iterStep)=gap;

            % record the number of zero group in the solution
            res.zg(iterStep)=penalty2(5);
        end

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end

        % restart the program every opts.rStartNum
        if (~mod(iterStep, opts.rStartNum))
            alphap=0; alpha=1;
            xp=x; Axp=Ax; xxp=zeros(n,1); L =L/2;
        end
    end
else
    error('\n The function does not support opts.mFlag neq 0 & opts.lFlag neq 0!');
end

author: thomaskuestner


text:
function writeCbToSBML(model,fileName)
%writeCbToSBML write a COBRA model to SBML
%
% writeCbToSBML(model,fileName)
%
%INPUTS
% model         COBRA model structure
% fileName      Name of xml file output
%
        sbmlModel = convertCobraToSBML(model);
        OutputSBML(sbmlModel, fileName);

author: opencobra

repo: cobratoolbox


text:
function value = i4vec_norm_l0 ( n, a )

%*****************************************************************************80
%
%% I4VEC_NORM_L0 returns the l0 "norm" of an I4VEC.
%
%  Discussion:
%
%    An I4VEC is a vector of I4's.
%
%    The l0 "norm" simply counts the number of nonzero entries in the vector.
%    It is not a true norm, but has some similarities to one.  It is useful
%    in the study of compressive sensing.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, integer A(N), the vector.
%
%    Output, integer I4VEC_NORM_L0, the value of the norm.
%
  value = sum ( a(1:n) ~= 0 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = identity ( m, n )

%*****************************************************************************80
%
%% IDENTITY returns the IDENTITY matrix.
%
%  Formula:
%
%    if ( I = J )
%      A(I,J) = 1
%    else
%      A(I,J) = 0
%
%  Example:
%
%    M = 4, N = 5
%
%    1 0 0 0 0
%    0 1 0 0 0
%    0 0 1 0 0
%    0 0 0 1 0
%
%  Rectangular properties:
%
%    A is integral: int ( A ) = A.
%
%    A is a zero/one matrix.
%
%  Square Properties:
%
%    A is nonsingular.
%
%    A is involutional: A * A = I.
%
%    A is diagonal.
%
%    Because A is diagonal, it has property A.
%
%    A is symmetric: A' = A.
%
%    A is a circulant matrix: each row is shifted once to get the next row.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    LAMBDA(1:N) = 1
%
%    The matrix of eigenvectors of A is A.
%
%    det ( A ) = 1.
%
%    A is unimodular.
%
%    For any vector v, A*v = v.
%
%    For any matrix B, A*B = B*A=B.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    A is centrosymmetric: A(I,J) = A(N+1-I,N+1-J).
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of
%    the matrix.
%
%    Output, real A(M,N), the matrix.
%
  a = zeros ( m, n );

  for i = 1 : m
    for j = 1 : n
      if ( i == j )
        a(i,j) = 1.0;
      else
        a(i,j) = 0.0;
      end
    end
  end

  return
end


text:
function [area,centroids] = calcVoronoiArea(v,varargin)
% compute the spherical area of the Voronoi decomposition
%
% Input
%  v - @vector3d
%
% Output
%  area - area of the corresponding Voronoi cells
%  centroids - centroid of the voronoi cell
%
% Options
% incomplete -

v = reshape(v,[],1);
N = length(v);

% in case of antipodal symmetry - add antipodal points
antipodal = v.antipodal || check_option(varargin, 'antipodal');
if antipodal
  v.antipodal = false;
  [v,~,IC] = unique([v;-v],'noSymmetry');
end

[V,C] = calcVoronoi(v);

nd = ~cellfun('isempty',C);
v = v.subSet(nd);

last = cumsum(cellfun('prodofsize',C(nd)));

left = [C{nd}];
shift = 2:last(end)+1;             % that is the shift
shift(last) = [0;last(1:end-1)]+1; % and the last gets the first
right = left(shift);

center = cumsum([1 diff(shift)>1]);

va = v.subSet(center);
vb = V.subSet(left);
vc = V.subSet(right);    % next vertex around

% calculate the area for each triangle around generator (va)
A = real(sphericalTriangleArea(va,vb,vc));

if nargout>1
  [x,y,z]= double(A.*(va+vb+vc));
  x = full(sparse(center,1,x,length(S2G),1));
  y = full(sparse(center,1,y,length(S2G),1));
  z = full(sparse(center,1,z,length(S2G),1));
  centroids = vector3d(x,y,z);
end

% accumulate areas of spherical triangles around generator
A = full(sparse(center,1,A,length(v),1));

area = zeros(size(nd));
area(nd) = A(1:nnz(nd));

if antipodal
  idx = ( accumarray(IC,ones(size(IC))) == 2 ); % find all double occurences
  area(idx) = area(idx)/2; % halve their weight
  area = area(IC); % go back to original order
  area = area(1:N); % only the original nodes
end

author: mtex-toolbox

repo: mtex

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


text:
function res4=writeRes4(res4File,res4,MAX_COILS);
%  Write the new .res4 file.  Use ieee-be (big endian) format
%  Character-string output is done using function writeCTFstring which
%  checks that strings are the correct length for the .res4 file format.

%  Function calls: - writeCTFstring (included in this listing).

fid_res4=fopen(res4File,'w','ieee-be');
if fid_res4<0
  fprintf('writeCTFds (writeRes4): Could not open file %s\n',res4File);
  return
end

fwrite(fid_res4,[res4.header(1:7),char(0)],'uint8');   % 8-byte header

%  meg41GeneralResRec
res4.appName=writeCTFstring(res4.appName,-256,fid_res4);
res4.dataOrigin=writeCTFstring(res4.dataOrigin,-256,fid_res4);
res4.dataDescription=writeCTFstring(res4.dataDescription,-256,fid_res4);
fwrite(fid_res4,res4.no_trials_avgd,'int16');
res4.data_time=writeCTFstring(res4.data_time,255,fid_res4);
res4.data_date=writeCTFstring(res4.data_date,255,fid_res4);

% new_general_setup_rec_ext part of meg41GeneralResRec

fwrite(fid_res4,res4.no_samples,'int32');        % 4
fwrite(fid_res4,[res4.no_channels 0],'int16');   % 2*2
fwrite(fid_res4,res4.sample_rate,'double');      % 8
fwrite(fid_res4,res4.epoch_time,'double');       % 8
fwrite(fid_res4,[res4.no_trials 0],'int16');     % 2*2
fwrite(fid_res4,res4.preTrigPts,'int32');        % 4
fwrite(fid_res4,res4.no_trials_done,'int16');    % 2
fwrite(fid_res4,res4.no_trials_display,'int16'); % 2
fwrite(fid_res4,res4.save_trials,'int32');       % 4 CTFBoolean

% meg41TriggerData part of new_general_setup_rec_ext   10 bytes total
fwrite(fid_res4,res4.primaryTrigger,'int32');      % 4
fwrite(fid_res4,res4.triggerPolarityMask,'int32'); % 4
fwrite(fid_res4,[0 0],'uint8');                    % 2 bytes
% end of meg41TriggerData part of new_general_setup_rec_ext

fwrite(fid_res4,[0 0],'uint8');                    % 2 bytes padding
fwrite(fid_res4,[res4.trigger_mode 0],'int16');    % 2*2
fwrite(fid_res4,res4.accept_reject_Flag,'int32');  % 4 CTFBoolean
fwrite(fid_res4,[res4.run_time_display 0],'int16'); % 2*2

fwrite(fid_res4,res4.zero_Head_Flag,'int32');      % 4 CTFBoolean
fwrite(fid_res4,res4.artifact_mode,'int32');       % 4 CTFBoolean
%  end of new_general_setup_rec_ext part of meg41GeneralResRec

% meg4FileSetup part of meg41GeneralResRec
res4.nf_run_name=writeCTFstring(res4.nf_run_name,32,fid_res4);
res4.nf_run_title=writeCTFstring(res4.nf_run_title,-256,fid_res4);
res4.nf_instruments=writeCTFstring(res4.nf_instruments,32,fid_res4);
res4.nf_collect_descriptor=writeCTFstring(res4.nf_collect_descriptor,32,fid_res4);
res4.nf_subject_id=writeCTFstring(res4.nf_subject_id,-32,fid_res4);
res4.nf_operator=writeCTFstring(res4.nf_operator,32,fid_res4);
res4.nf_sensorFileName=writeCTFstring(res4.nf_sensorFileName,56,fid_res4);
res4.rdlen=length(res4.run_description);      % Run_description may have been changed.
fwrite(fid_res4,[0 res4.rdlen 0],'int32');    %  12-byte structure of the .res4 file.
res4.run_description=writeCTFstring(res4.run_description,res4.rdlen,fid_res4);
%  end of meg4FileSetup part of meg41GeneralResRec

%  filter descriptions

fwrite(fid_res4,res4.num_filters,'int16');
for kfilt=1:res4.num_filters
  fwrite(fid_res4,res4.filters(kfilt).freq,'double');
  fwrite(fid_res4,res4.filters(kfilt).fClass,'int32');
  fwrite(fid_res4,res4.filters(kfilt).fType,'int32');
  fwrite(fid_res4,res4.filters(kfilt).numParam,'int16');
  for kparm=1:res4.filters(kfilt).numParam
    fwrite(fid_res4,res4.filters(kfilt).Param(kparm),'double');
  end
end

%  Write channel names.   Must have size(res4.chanNames)=[nChan 32]
[no_chanNames len_chanNames]=size(res4.chanNames);
for kchan=1:res4.no_channels
  res4.chanNames(kchan,:)=writeCTFstring(res4.chanNames(kchan,:),len_chanNames,fid_res4);
end

%  Write sensor resource table
for kchan=1:res4.no_channels
  fwrite(fid_res4,res4.senres(kchan).sensorTypeIndex,'int16');
  fwrite(fid_res4,res4.senres(kchan).originalRunNum,'int16');
  fwrite(fid_res4,res4.senres(kchan).coilShape,'int32');
  fwrite(fid_res4,res4.senres(kchan).properGain,'double');
  fwrite(fid_res4,res4.senres(kchan).qGain,'double');
  fwrite(fid_res4,res4.senres(kchan).ioGain,'double');
  fwrite(fid_res4,res4.senres(kchan).ioOffset,'double');
  fwrite(fid_res4,res4.senres(kchan).numCoils,'int16');
  numCoils=res4.senres(kchan).numCoils;
  fwrite(fid_res4,res4.senres(kchan).grad_order_no,'int16');
  fwrite(fid_res4,0,'int32');  % Padding to 8-byte boundary

  % coilTbl
  for qx=1:numCoils
    fwrite(fid_res4,[res4.senres(kchan).pos0(:,qx)' 0],'double');
    fwrite(fid_res4,[res4.senres(kchan).ori0(:,qx)' 0],'double');
    fwrite(fid_res4,[res4.senres(kchan).numturns(qx) 0 0 0],'int16');
    fwrite(fid_res4,res4.senres(kchan).area(qx),'double');
  end
  if numCoils<MAX_COILS
    fwrite(fid_res4,zeros(10*(MAX_COILS-numCoils),1),'double');
  end

  %HdcoilTbl
  for qx=1:numCoils
    fwrite(fid_res4,[res4.senres(kchan).pos(:,qx)' 0],'double');
    fwrite(fid_res4,[res4.senres(kchan).ori(:,qx)' 0],'double');
    fwrite(fid_res4,[res4.senres(kchan).numturns(qx) 0 0 0],'int16');
    fwrite(fid_res4,res4.senres(kchan).area(qx),'double');
  end
  if numCoils<MAX_COILS
    fwrite(fid_res4,zeros(10*(MAX_COILS-numCoils),1),'double');
  end
end
%  End writing sensor resource table

%  Write the table of balance coefficients.
if res4.numcoef<=0
  fwrite(fid_res4,res4.numcoef,'int16');  % Number of coefficient records
elseif res4.numcoef>0
  scrx_out=[];
  for kx=1:res4.numcoef
    sName=strtok(char(res4.scrr(kx).sensorName),['- ',char(0)]);
    if ~isempty(strmatch(sName,res4.chanNames))
      scrx_out=[scrx_out kx];
    end
  end
  %  Remove the extra coefficient records
  res4.scrr=res4.scrr(scrx_out);
  res4.numcoef=size(res4.scrr,2);
  fwrite(fid_res4,res4.numcoef,'int16');  % Number of coefficient records
  %  Convert res4.scrr to double before writing to output file.  In MATLAB 5.3.1,
  %  when the 'ieee-be' option is applied, fwrite cannot write anything except
  %  doubles and character strings, even if fwrite does allow you to specify short
  %  integer formats in the output file.
  for nx=1:res4.numcoef
    fwrite(fid_res4,double(res4.scrr(nx).sensorName),'uint8');
    fwrite(fid_res4,[double(res4.scrr(nx).coefType) 0 0 0 0],'uint8');
    fwrite(fid_res4,double(res4.scrr(nx).numcoefs),'int16');
    fwrite(fid_res4,double(res4.scrr(nx).sensor),'uint8');
    fwrite(fid_res4,res4.scrr(nx).coefs,'double');
  end
end
fclose(fid_res4);
return

% *************************************************************************************
%*************** Function writeCTFstring ************************************************
function strng=writeCTFstring(instrng,strlen,fid)

%  Writes a character string to output unit fid.  Append nulls to get the correct length.
%  instrng : Character string.  size(instrng)=[nLine nPerLine].  strng is reformulated as a
%            long string of size [1 nChar].  Multiple lines are allowed so the user can
%            easily append text.  If necessary, characters are removed from
%            instrng(1:nLine-1,:) so all of strng(nLine,:) can be accomodated.
%    strlen: Number of characters to write.  strlen<0 means remove leading characters.
%        If abs(strlen)>length(strng) pad with nulls (char(0))
%        If 0<strlen<length(strng), truncate strng and terminate with a null.
%        If -length(string)<strlen<0, remove leading characters and terminate with a null.

%  Form a single long string
nLine=size(instrng,1);
strng=deblank(instrng(1,:));

if nLine>1
  %  Concatenate lines 1:nLine-1
  for k=2:nLine-1
    if length(strng)>0
      if ~strcmp(strng(length(strng)),'.') & ~strcmp(strng(length(strng)),',')
        strng=[strng '.'];   % Force a period at the end of the string.
      end
    end
    strng=[strng '  ' deblank(instrng(k,:))];
  end

  if length(strng)>0
    if ~strcmp(strng(length(strng)),'.')  % Force a period at the end of the string.
      strng=[strng '.'];
    end
  end
  %  Add all of the last line.
  nChar=length(strng);
  nLast=length(deblank(instrng(nLine,:)));
  strng=[strng(1:min(nChar,abs(strlen)-nLast-4)) '  ' deblank(instrng(nLine,:))];
end

if length(strng)<abs(strlen)
  strng=[strng char(zeros(1,abs(strlen)-length(strng)))];
elseif length(strng)>=strlen & strlen>0
  strng=[strng(1:strlen-1) char(0)];
else
   strng=[strng(nLast+[strlen+2:0]) char(0)];
end

fwrite(fid,strng,'char');
return
%  ************* End of function writeCTFstring ***************************************
% ************************************************************************************

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% read_eeglabheader() - import EEGLAB dataset files
%
% Usage:
%   >> header = read_eeglabheader(filename);
%
% Inputs:
%   filename - [string] file name
%
% Outputs:
%   header   - FILEIO toolbox type structure
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, 2008-

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD, 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

function header = read_eeglabheader(filename)

if nargin < 1
  help read_eeglabheader;
  return;
end;

if ~isstruct(filename)
  load('-mat', filename);
else
  EEG = filename;
end;

header.Fs          = EEG.srate;
header.nChans      = EEG.nbchan;
header.nSamples    = EEG.pnts;
header.nSamplesPre = -EEG.xmin*EEG.srate;
header.nTrials     = EEG.trials;
try
  header.label       = { EEG.chanlocs.labels }';
catch
  warning('creating default channel names');
  for i=1:header.nChans
    header.label{i} = sprintf('chan%03d', i);
  end
end
ind = 1;
for i = 1:length( EEG.chanlocs )
    if isfield(EEG.chanlocs(i), 'X') && ~isempty(EEG.chanlocs(i).X)
        header.elec.label{ind, 1} = EEG.chanlocs(i).labels;
        % this channel has a position
        header.elec.pnt(ind,1) = EEG.chanlocs(i).X;
        header.elec.pnt(ind,2) = EEG.chanlocs(i).Y;
        header.elec.pnt(ind,3) = EEG.chanlocs(i).Z;
        ind = ind+1;
    end;
end;

% remove data
% -----------
%if isfield(EEG, 'datfile')
%    if ~isempty(EEG.datfile)
%        EEG.data = EEG.datfile;
%    end;
%else
%    EEG.data = 'in set file';
%end;
EEG.icaact = [];

header.orig = EEG;

author: PatternRecognition

repo: OpenBMI


text:
% isnumeric() - returns 1
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Nov. 2008

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD
%
% 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

function r = isnumeric(a)

    r = 1;


text:
function [isvalid msg]  = verify_tracks(filename,targettype)
% function [isvalid msg] = verify_tracks(filename,targettype)
% Check integrity of tracks file (see published doc as a guide).
%
% filename: full path to trx file
% targettype: JLabelData targetType
%
% isvalid: true if trx file passes all checks, false otherwise
% msg: if isvalid is false, string containing reason(s) for failure.

TRX_FLDS_REQ = {'x' 'y' 'theta' 'a' 'b' 'nframes' 'firstframe' 'endframe' ...
  'off' 'id' 'x_mm' 'y_mm' 'theta_mm' 'a_mm' 'b_mm' 'sex' 'dt' 'fps'};
TRX_FLDS_OPT = {'timestamps'};
TRX_FLDS_LAR = {'area' 'xcontour' 'ycontour' 'xspine' 'yspine' 'area_mm' ...
  'xspine_mm' 'yspine_mm'};

try
  mat = load(filename','-mat');
catch %#ok<CTCH>
  isvalid = false;
  msg = sprintf('Can''t load MAT file ''%s''.',filename);
  return;
end

if ~isfield(mat,'trx')
  isvalid = false;
  msg = 'MAT file does not contain variable ''trx''.';
  return;
end
trx = mat.trx;

% check trx fields
tf = isfield(trx,TRX_FLDS_REQ);
if ~all(tf)
  isvalid = false;
  msg = sprintf('trx struct array missing required fields: %s',...
    civilizedStringFromCellArrayOfStrings(TRX_FLDS_REQ(~tf)));
  return;
end

unkFlds = setdiff(fieldnames(trx),[TRX_FLDS_REQ TRX_FLDS_OPT TRX_FLDS_LAR]);
if ~isempty(unkFlds)
  warning('verify_tracks:unknownFields','Unknown fields in trx struct array: %s',...
    civilizedStringFromCellArrayOfStrings(unkFlds));
end

Ntrx = numel(trx);
issues = cell(0,1);
NFRAMES_ROW_VECS = {'x' 'y' 'theta' 'a' 'b' 'x_mm' 'y_mm' 'theta_mm' 'a_mm' 'b_mm'};
VALID_SEX_VALS = {'M' 'F' '?'};
% Check each trx element. This is not a 100% comprehensive check.
for iTrx = 1:Ntrx

  trxEl = trx(iTrx);
  nf = trxEl.nframes;

  % sizes of row vectors
  for fld = NFRAMES_ROW_VECS, fld = fld{1}; %#ok<*FXSET>
    val = trxEl.(fld);
    issues = lclChk(issues,iTrx,isrow(val) && numel(val)==nf,...
      sprintf('trx field ''%s'' has unexpected size',fld));
  end

  % firstframe/endframe/nframes/off
  issues = lclChk(issues,iTrx,(trxEl.endframe-trxEl.firstframe+1)==trxEl.nframes,...
    'inconsistent values for nframes/firstframe/endframe');
  issues = lclChk(issues,iTrx,trxEl.off==1-trxEl.firstframe,...
    'inconsistent values for firstframe/off');

  val = trxEl.sex;
  issues = lclChk(issues,iTrx,...
    (isscalar(val) || isvector(val) && iscellstr(val) && numel(val)==nf) ...
      && all(ismember(val,VALID_SEX_VALS)),'invalid value for field ''sex''');

  val = trxEl.dt;
  issues = lclChk(issues,iTrx,isrow(val) && numel(val)==nf-1,...
      'field ''dt'' has unexpected size');
end

isvalid = isempty(issues);
msg = sprintf('%s\n',issues{:});
msg = msg(1:end-1);

function issues = lclChk(issues,iTrx,tf,str)
if ~tf
  issues{end+1,1} = sprintf('Trx %d: %s',iTrx,str);
end

author: kristinbranson

repo: JAABA


text:
function ex1bvp
%EX1BVP  Example 1 of the BVP tutorial.
%   This is the example for MUSN in U. Ascher, R. Mattheij, and R. Russell,
%   Numerical Solution of Boundary Value Problems for Ordinary Differential
%   Equations, SIAM, Philadelphia, PA, 1995.  MUSN is a multiple shooting
%   code for nonlinear BVPs.  The problem is
%
%      u' =  0.5*u*(w - u)/v
%      v' = -0.5*(w - u)
%      w' = (0.9 - 1000*(w - y) - 0.5*w*(w - u))/z
%      z' =  0.5*(w - u)
%      y' = -100*(y - w)
%
%   The interval is [0 1] and the boundary conditions are
%
%      u(0) = v(0) = w(0) = 1,  z(0) = -10,  w(1) = y(1)
%
%   The example uses a guess for the solution coded here in EX1INIT.
%   The results of a run of the FORTRAN code MUSN are here compared to
%   the curves produced by BVP4C.

% Copyright 2004, The MathWorks, Inc.

solinit = bvpinit(linspace(0,1,5),@ex1init);
options = bvpset('Stats','on','RelTol',1e-5);

sol = bvp4c(@ex1ode,@ex1bc,solinit,options);

% The solution at the mesh points
x = sol.x;
y = sol.y;

% Solution obtained using MUSN:
amrx = [ 0. .1 .2 .3 .4 .5 .6 .7 .8 .9 1.]';
amry = [1.00000e+00   1.00000e+00   1.00000e+00  -1.00000e+01   9.67963e-01
        1.00701e+00   9.93036e-01   1.27014e+00  -9.99304e+00   1.24622e+00
        1.02560e+00   9.75042e-01   1.47051e+00  -9.97504e+00   1.45280e+00
        1.05313e+00   9.49550e-01   1.61931e+00  -9.94955e+00   1.60610e+00
        1.08796e+00   9.19155e-01   1.73140e+00  -9.91915e+00   1.72137e+00
        1.12900e+00   8.85737e-01   1.81775e+00  -9.88574e+00   1.80994e+00
        1.17554e+00   8.50676e-01   1.88576e+00  -9.85068e+00   1.87957e+00
        1.22696e+00   8.15025e-01   1.93990e+00  -9.81503e+00   1.93498e+00
        1.28262e+00   7.79653e-01   1.98190e+00  -9.77965e+00   1.97819e+00
        1.34161e+00   7.45374e-01   2.01050e+00  -9.74537e+00   2.00827e+00
        1.40232e+00   7.13102e-01   2.02032e+00  -9.71310e+00   2.02032e+00];

% Shift up the fourth component for the plot.
amry(:,4) = amry(:,4) + 10;
y(4,:) = y(4,:) + 10;

figure
plot(x,y',amrx,amry,'*')
axis([0 1 -0.5 2.5])
title('Example problem for MUSN')
ylabel('bvp4c and MUSN (*) solutions')
xlabel('x')

% --------------------------------------------------------------------------

function dydx = ex1ode(x,y)
%EX1ODE  ODE function for Example 1 of the BVP tutorial.
%   The components of y correspond to the original variables
%   as  y(1) = u, y(2) = v, y(3) = w, y(4) = z, y(5) = y.
dydx =  [ 0.5*y(1)*(y(3) - y(1))/y(2)
         -0.5*(y(3) - y(1))
         (0.9 - 1000*(y(3) - y(5)) - 0.5*y(3)*(y(3) - y(1)))/y(4)
          0.5*(y(3) - y(1))
          100*(y(3) - y(5)) ];

%-------------------------------------------------------------------------

function res = ex1bc(ya,yb)
%EX1BC  Boundary conditions for Example 1 of the BVP tutorial.
%   RES = EX1BC(YA,YB) returns a column vector RES of the
%   residual in the boundary conditions resulting from the
%   approximations YA and YB to the solution at the ends of
%   the interval [a b]. The BVP is solved when RES = 0.
%   The components of y correspond to the original variables
%   as  y(1) = u, y(2) = v, y(3) = w, y(4) = z, y(5) = y.
res = [ ya(1) - 1
        ya(2) - 1
        ya(3) - 1
        ya(4) + 10
        yb(3) - yb(5)];

%-------------------------------------------------------------------------

function v = ex1init(x)
%EX1INIT  Guess for the solution of Example 1 of the BVP tutorial.
v = [        1
             1
     -4.5*x^2+8.91*x+1
            -10
     -4.5*x^2+9*x+0.91 ];

author: Sable

repo: mcbench-benchmarks


text:
function runIntersect = mintersect(varargin)
%MINTERSECT Multiple set intersection.
%   MINTERSECT(A,B,C,...) when A,B,C... are vectors returns the values
%   common to all A,B,C... The result will be sorted.  A,B,C... can be cell
%   arrays of strings.
%
%   MINTERSECT repeatedly evaluates INTERSECT on successive pairs of sets,
%   which may not be very efficient.  For a large number of sets, this should
%   probably be reimplemented using some kind of tree algorithm.
%
%   MINTERSECT(A,B,'rows') when A,B,C... are matrices with the same
%   number of columns returns the rows common to all A,B,C...
%
%   See also INTERSECT

flag = 0;
if isempty(varargin),
    error('No inputs specified.')
else
    if isequal(varargin{end},'rows'),
        flag = 'rows';
        setArray = varargin(1:end-1);
    else
        setArray = varargin;
    end
end

runIntersect = setArray{1};
for i = 2:length(setArray),

    if isequal(flag,'rows'),
        runIntersect = intersect(runIntersect,setArray{i},'rows');
    elseif flag == 0,
        runIntersect = intersect(runIntersect,setArray{i});
    else
        error('Flag not set.')
    end

    if isempty(runIntersect),
        return
    end

end


text:
function cond = condition_sample1 ( n, a, m )

%*****************************************************************************80
%
%% CONDITION_SAMPLE1 estimates the L1 condition number of a matrix.
%
%  Discussion:
%
%    A naive sampling method is used.
%
%    Only "forward" sampling is used, that is, we only look at results
%    of the form y=A*x.
%
%    Presumably, solving systems A*y=x would give us a better idea of
%    the inverse matrix.
%
%    Moreover, a power sequence y1 = A*x, y2 = A*y1, ... and the same for
%    the inverse might work better too.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the matrix.
%
%    Input, real A(N,N), the matrix.
%
%    Input, integer M, the number of samples to use.
%
%    Output, real COND, an estimate of the L1 condition number.
%
  a_norm = 0.0;
  ainv_norm = 0.0;
  seed = 123456789;

  for i = 1 : m

    [ x, seed ] = r8vec_uniform_unit ( n, seed );
    x_norm = norm ( x, 1 );
    ax = a * x;
    ax_norm = norm ( ax, 1 );

    if ( ax_norm == 0.0 )
      cond = 0.0;
      return
    end

    a_norm    = max ( a_norm,    ax_norm / x_norm  );
    ainv_norm = max ( ainv_norm, x_norm  / ax_norm );

  end

  cond = a_norm * ainv_norm;

  return
end


text:
function CRB = crb_general_det_1d(design, wavelength, doas, P_est, noise_var, snapshot_count)
%CRB_GENERAL_DET_1D CRB for general 1D arrays based on the deterministic
%(conditional) model, in radians.
%Inputs:
%   design - Array design.
%   wavelength - Wavelength.
%   doas - DOA vector in radians.
%   P_est - Estimated source covariance matrix from deterministic source
%           signals.
%   noise_var - Noise power.
%   snapshot_count - (Optional) number of snapshots. Default is one.
%Reference:
%   [1] P. Stoica and A. Nehorai, "Performance study of conditional and
%       unconditional direction-of-arrival estimation," IEEE Transactions
%       on Acoustics, Speech and Signal Processing, vol. 38, no. 10,
%       pp. 1783-1795, Oct. 1990.
if design.dim ~= 1
    error('1D array expected.');
end
if nargin <= 5
    snapshot_count = 1;
end
[A, D] = steering_matrix(design, wavelength, doas);
[m, k] = size(A);
H = D'*(eye(m) - A/(A'*A)*A')*D;
CRB = real(H .* P_est.');
CRB = eye(k) / CRB * (noise_var / snapshot_count / 2);
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


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

% FT_ARTIFACT_NAN identifies artifacts that are indicated in the data as NaN (not a
% number) values.
%
% Use as
%   [cfg, artifact] = ft_artifact_nan(cfg, data)
% where the input data is a structure as obtained from FT_REJECTARTIFACT with the
% option cfg.artfctdef.reject='nan', or from FT_REJECTVISUAL with cfg.keeptrial='nan'
% or cfg.keepchannel='nan'.
%
% The configuration can contain
%   cfg.artfctdef.nan.channel = Nx1 cell-array with selection of channels, see FT_CHANNELSELECTION for details
%
% 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, you can use
%   cfg.inputfile   =  ...
% to read the input data 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) 2017, 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 debug
ft_preamble loadvar    data
ft_preamble provenance data

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

% set the default options
cfg.feedback              = ft_getopt(cfg, 'feedback', 'text');
cfg.artfctdef             = ft_getopt(cfg, 'artfctdef');
cfg.artfctdef.nan         = ft_getopt(cfg.artfctdef, 'nan');
cfg.artfctdef.nan.channel = ft_getopt(cfg.artfctdef.nan, 'channel', {'all'});

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

% read the header, or get it from the input data
if ~hasdata
  cfg = ft_checkconfig(cfg, 'dataset2files', 'yes');
  cfg = ft_checkconfig(cfg, 'required', {'headerfile', 'datafile'});
  hdr = ft_read_header(cfg.headerfile, 'headerformat', cfg.headerformat);
else
  data = ft_checkdata(data, 'datatype', 'raw', 'hassampleinfo', 'yes');
  cfg  = ft_checkconfig(cfg, 'forbidden', {'dataset', 'headerfile', 'datafile'});
  hdr  = ft_fetch_header(data);
end

% set default cfg.continuous
if ~isfield(cfg, 'continuous')
  if hdr.nTrials==1
    cfg.continuous = 'yes';
  else
    cfg.continuous = 'no';
  end
end

% get the specification of the data segments that should be scanned for artifacts
if ~isfield(cfg, 'trl') && hasdata
  trl = data.sampleinfo;
  for k = 1:numel(data.trial)
    trl(k,3) = time2offset(data.time{k}, data.fsample);
  end
elseif isfield(cfg, 'trl') && ischar(cfg.trl)
  trl = loadvar(cfg.trl, 'trl');
elseif isfield(cfg, 'trl') && isnumeric(cfg.trl)
  trl = cfg.trl;
else
  ft_error('cannot determine which segments of data to scan for artifacts');
end

% get the remaining settings
artfctdef     = cfg.artfctdef.nan;
artfctdef.trl = trl;
numtrl        = size(trl,1);
channel       = ft_channelselection(artfctdef.channel, hdr.label);
chanindx      = match_str(hdr.label, channel);
nchan         = length(chanindx);
artifact      = zeros(0,2);

ft_progress('init', cfg.feedback, ['searching for artifacts in ' num2str(nchan) ' channels']);
for trlop=1:numtrl
  ft_progress(trlop/numtrl, 'searching in trial %d from %d\n', trlop, numtrl);
  if hasdata
    dat = ft_fetch_data(data,        'header', hdr, 'begsample', trl(trlop,1), 'endsample', trl(trlop,2), 'chanindx', chanindx, 'checkboundary', strcmp(cfg.continuous, 'no'));
  else
    dat = ft_read_data(cfg.datafile, 'header', hdr, 'begsample', trl(trlop,1), 'endsample', trl(trlop,2), 'chanindx', chanindx, 'checkboundary', strcmp(cfg.continuous, 'no'), 'dataformat', cfg.dataformat);
  end

  tmp = any(isnan(dat),1);
  if any(tmp)
    % there can be multiple segments with nans
    begsample = find(diff([0 tmp])>0);
    endsample = find(diff([tmp 0])<0);
    for j=1:numel(begsample)
      artifact(end+1,:) = [begsample(j) endsample(j)] + data.sampleinfo(trlop,1) - 1;
    end
  end
end % for trlop
ft_progress('close');

% remember the details that were used here
cfg.artfctdef.nan          = artfctdef;
cfg.artfctdef.nan.artifact = artifact;

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   data
ft_postamble savevar


text:
function labels = getClassLabels()

labels={'ped', ...			% 1
	'person_on_vhcl', ...	% 2
	'car', ...				% 3
	'bicycle', ...			% 4
	'mbike', ...			% 5
	'non_mot_vhcl', ...		% 6
	'static_person', ...	% 7
	'distractor', ...		% 8
	'occluder', ...			% 9
	'occluder_on_grnd', ...		%10
	'occluder_full', ...		% 11
	'reflection', ...		% 12
	'crowd' ...			% 13
	};

end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
% loaddat() - loading neuroscan format data file into matlab.
%
% Usage:
%   >> [typeeeg, rt, response, n] = loaddat( filename );
%
% Inputs:
%   filename     - input Neuroscan .dat file
%
% Outputs:
% 	typeeeg	 - type of the sweeps
%	rt	 - reaction time of the subject
%	response - response of the subject
%	n	 - number of sweeps
%
% See also: loadeeg()

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, 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

function [tmptypeeeg, tmprt, tmpresponseeeg, n] = loaddat( FILENAME)

if nargin<1
	fprintf('Not enought arguments\n');
	help loaddat
	return;
end;

BOOL='int16';
ULONG='int32';
FLOAT='float32';

fid=fopen(FILENAME,'r','ieee-le');
if fid<0
	fprintf(2,['Error LOADEEG: File ' FILENAME ' not found\n']);
	return;
end;

% skip the first 20 lines
% -----------------------
for index=1:20	fgetl(fid); end;

% read the matrix
% ---------------
tmpMat 		    = fscanf(fid, '%f', [5, inf]);
tmptypeeeg  	= tmpMat(3,:);
tmpresponseeeg  = tmpMat(4,:);
tmprt       	= tmpMat(5,:);
n        	    = size( tmpMat, 2);

fclose(fid);
return;


text:
function [A,B,C]=triAngles4Sides(a,b,c)
%%TRIANGLES4SIDES Given the sides of a triangle, find the angles of the
%         triangle. This just applies the law of cosines.
%
%INPUTS: a, b, c The scalar lengths of the sides of the triangle in any
%                order.
%
%OUTPUTS: A, B, C The angles of the triangle in radians. A is the angle
%                 opposite side a, B the angle opposite side b and C the
%                 angle opposite side C.
%
%The law of cosines is given in [1], among many other places.
%
%EXAMPLE:
%This recovers the angles of a 30-60-90 traingle.
% a=1;
% b=sqrt(3);
% c=2;
% [A,B,C]=triAngles4Sides(a,b,c);
% ADeg=rad2deg(A)
% BDeg=rad2deg(B)
% CDeg=rad2deg(C)
%
%REFERENCES:
%[1] Weisstein, Eric W. "Law of Cosines." From MathWorld--A Wolfram Web
%    Resource. https://mathworld.wolfram.com/LawofCosines.html
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

a2=a.*a;
b2=b.*b;
c2=c.*c;

A=acos((-a2+b2+c2)./(2*b.*c));
B=acos((a2-b2+c2)./(2*a.*c));
C=pi-A-B;
%We could have used:
%C=acos((a2+b2-c2)./(2*a.*b));
%But we know that the angles in a triangle in Euclidean geometry sum to 180
%degrees.

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [elem2dof,dofSign,edge] = dofBDM1(elem)
%% DOFBDM1 dof structure for BDM1 element
%
% [elem2dof,dofSign,edge] = dofBDM1(elem) constructs data structure for
% the BDM face element in 3-D. elem is the connectivity matrix for a 3-D
% triangulation. elem2dof is the elementwise pointer from elem to dof
% indices. dofSign records the consistency of the local and global edge
% orientation. edge is the edge matrix.
%
% Added by Ming Wang.
%
%  See also dof3BDM1, dofRT0, dof3RT0
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

[elem2edge,dofSign,edge] = dofRT0(elem);
NT = size(elem,1);  NE = size(edge,1);
totaledge = uint32([elem(:,[2 3]); elem(:,[3 1]); elem(:,[1 2])]);
[tempvar,i] = sort(totaledge,2); %#ok<*ASGLU>
[tempvar,j]= sort(i,2);
elem2dof = repmat(elem2edge(:),1,2) + uint32((j-1)*NE);
elem2dof = reshape(elem2dof(:),NT,6);
dofSign = repmat(dofSign,1,2);

repo: ifem


text:
function res = imrecerode(img, h, varargin)
%IMRECERODE Perform a morphological reconstruction by erosion
%
%   usage:
%   RES = imrecerode(IMG, H)
%   performs morphological reconstruction by erosion of img+h over img.
%
%   RES = imrecerode(IMG, H, CONN)
%   specifies connectivity, by default 8 or 26, depending on image
%   dimension
%

% Process input arguments
if length(size(img))==3
    conn = 26;
else
    conn = 8;
end

if ~isempty(varargin)
    conn=varargin{1};
end

% apply reconstruction
marker = imcomplement(img);
res = imcomplement(imreconstruct(marker, imadd(marker, h), conn));

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
classdef AnalysisMagnitudeHist < HistAnalysisWindow
    methods
        function obj=AnalysisMagnitudeHist(ax)
            obj@HistAnalysisWindow(ax);
        end

        function prepare_axes(obj)
            % prepare the cumulative event rate axes
            if isempty(obj.ax.Tag)
                obj.ax.Tag = 'dvCumrate';
            end
            obj.ax.Title.String='Cumulative Rate';
            obj.ax.XLabel.String='Time';
            obj.ax.YLabel.String='Cumulative Events';
        end

        function [x,y]=calculate(obj,catalog)
            % calculate the cumulative number of events (y) over time (x)
            x=sort(catalog.Date);
            y=[]; % histogram?
        end

    end
end

author: CelsoReyes


text:
function chooseBestModels_HN(pathExperimentsBinary,pathExperimentsTime,fSetNames,nameOutcomes,metric,maxOrder)
% -------------------------------------------------------------------------
% function chooseBestModels_HN(pathWORK,nExp,fSetName,nameOutcomes,maxOrder)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function requires user input to choose the set types and model orders
% providing the best parsimonious models for all outcomes analyzed in the
% HN study. See ref.[1] for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathExperiments: Full path to the directory containing all experiments.
%                     --> Ex: '/myProject/WORKSPACE/CV-BASED_RESULTS'
% 2. nExp: Numerical value specifying the number of experiments to analyze.
%          --> Ex: 10
% 3. fSetName: Cell of strings specifying the name of the type of feature
%              set analyzed.
%              --> Ex: {'PET','CT','SEPARATE','FUSED'}
% 4. nameOutcomes: Cell of strings specifying the outcome names to analyze.
%                  --> Ex: {'Failure','Locoregional','Distant','Death'}
% 5. maxOrder: Integer specifying the maximal multivariable model order.
%              --> Ex: 10
% -------------------------------------------------------------------------
% OUTPUTS: Final prediction models saved in a folder named 'FINAL_MODELS/nameOutcome/fSetName'
%          in the given experiment folder.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% 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/>.
% -------------------------------------------------------------------------

startpath = pwd;
cd(pathExperimentsBinary)
nOutcomes = length(nameOutcomes);
nType = length(fSetNames);

fprintf('\n')
cd(pathExperimentsTime), pathExperimentTime = pathExperimentsTime;
cd(pathExperimentTime), mkdir('FINAL_MODELS'), cd('FINAL_MODELS'), pathFinalModelsTime = pwd;
cd(pathExperimentsBinary), pathExperimentBinary = pathExperimentsBinary;
cd(pathExperimentBinary), mkdir('FINAL_MODELS'), cd('FINAL_MODELS'), pathFinalModelsBinary = pwd;
for o = 1:nOutcomes
    cd(pathFinalModelsTime), mkdir(nameOutcomes{o}), cd(fullfile(pathFinalModelsTime,nameOutcomes{o})), pathOutcomeTime = pwd;
    cd(pathFinalModelsBinary), mkdir(nameOutcomes{o}), cd(fullfile(pathFinalModelsBinary,nameOutcomes{o})), pathOutcomeBinary = pwd;
    plotPredictionResults_HN(fullfile(pathExperimentBinary,'RESULTS'),nameOutcomes{o},fSetNames,metric,maxOrder)
    fprintf('\n======== DISPLAYING PREDICTION RESULTS FOR "%s" OUTCOME, "%s" METRIC  ========\n',nameOutcomes{o},metric)
    for f = 1:nType
        cd(pathOutcomeTime), mkdir(fSetNames{f}), cd(fullfile(pathOutcomeTime,fSetNames{f})), pathFinalTime = pwd;
        cd(pathOutcomeBinary), mkdir(fSetNames{f}), cd(fullfile(pathOutcomeBinary,fSetNames{f})), pathFinalBinary = pwd;
        while 1
            order = input(['Which model order of the ',fSetNames{f},' feature set provides the best parsimonious model? \n' ...
                           '--> Type a number between 1 to ',num2str(maxOrder),' and press ENTER \n' ...
                           'ANSWER: ']);
            fprintf('\n')
            if isnumeric(order) && order <= maxOrder && order >= 1
                break
            end
        end
        cd(fullfile(pathExperimentBinary,'RESULTS'))
        results = load(['RESULTS_',fSetNames{f},'_',nameOutcomes{o}]); results = struct2cell(results); results = results{1};
        finalModel = results.(['Order',num2str(order)]);
        finalModel.Order = order;
        finalModel.outcome = nameOutcomes{o};
        cd(pathFinalBinary), save('finalModel','finalModel')
        cd(pathFinalTime), save('finalModel','finalModel')
    end
    close all
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function []=anim8_DIC3DPP_faceMeasureDirection(DIC3DPPresults,faceMeasureString,RBMlogic,varargin)
%% function for plotting 3D-DIC results of face measures as color + direction as arrows in STEP4.
% plotting 3D surfaces from camera pairs, animation changing
% with time, and the faces colored according to faceMeasureString
% this function is called in plotMultiDICPairResults
%
% Options:
% anim8_DIC3DPP_faceMeasureDirection(DIC_3Dallpairs_results,faceMeasureString)
% anim8_DIC3DPP_faceMeasureDirection(DIC_3Dallpairs_results,faceMeasureString,optStruct)
%
% Inputs:
% * DIC3DAllPairsResults
% * faceMeasureString: can be any of the following:
%   'Epc1','Epc2','epc1','epc2','Lamda1','Lamda2', or a combination of two  in a cell array to plot side by side
% * optStruct: optional structure for plotting options which may include any of the following fields:
%   - smoothLogic: logical variable for smoothing (true)/not smoothing (false) the face measure
%   - FaceAlpha: transparacy of the faces (scalar between 0 and 1, where zero is transparent and 1 is opaque)
%   - colorBarLimits: a 2x1 scalar vector for the colobar limits. if not set, it's automatic
%   - dataLimits: a 2x1 scalar vector for the data limits of the face measure. if a face measure is outside these limits, it is set to NaN. if not set no face is set to NaN
%   - colorMap
%   - zDirection: 1 for z up and -1 for z down
%   - lineColor: line color for the mesh. can be for example 'b','k','none',etc...
%   - supTitleString=faceMeasureString;

%% Assign plot options
Narg=numel(varargin);
switch Narg
    case 1
        optStruct=varargin{1};
    case 0
        optStruct=struct;
    otherwise
        ('wrong number of input arguments');
end

% complete the struct fields
if ~isfield(optStruct,'smoothLogic')
    optStruct.smoothLogic=0;
end
if ~isfield(optStruct,'FaceAlpha')
    optStruct.FaceAlpha=1;
end
if ~isfield(optStruct,'dataLimits')
    optStruct.dataLimits=[-inf inf];
end
if ~isfield(optStruct,'zDirection') % 1 or -1
    optStruct.zDirection=1;
end
if ~isfield(optStruct,'lineColor') % 'none' or 'k'
    optStruct.lineColor='none';
end
if ~isfield(optStruct,'maxCorrCoeff')
    optStruct.maxCorrCoeff=[];
end
if ~isfield(optStruct,'quiverScaleFactor')
    optStruct.quiverScaleFactor=20;
end
%%
nFrames=numel(DIC3DPPresults.Points3D);

[xl,yl,zl]=axesLimits(DIC3DPPresults.Points3D);
meanEdgeLength=nanmean(patchEdgeLengths(DIC3DPPresults.Faces,DIC3DPPresults.Points3D{1}));

%% Assign the right face measure into FC

if iscell(faceMeasureString)
    faceMeasureCell=faceMeasureString;
    nStrains=numel(faceMeasureString);
    switch nStrains
        case 1
            FC=cell(nFrames,nStrains);
            D=cell(nFrames,nStrains);
            Ds=cell(nFrames,nStrains);
            Vc=cell(nFrames,nStrains);
        case 2
            FC=cell(nFrames,nStrains);
            D=cell(nFrames,nStrains);
            Ds=cell(nFrames,nStrains);
            Vc=cell(nFrames,nStrains);
        otherwise
            error('wrong number of cells in faceMeasureString cell array (must be 1 or 2)');
    end
elseif ischar(faceMeasureString)
    nStrains=1;
    faceMeasureCell=cell(1);
    faceMeasureCell{1}=faceMeasureString;
    FC=cell(nFrames,1);
    D=cell(nFrames,1);
    Ds=cell(nFrames,1);
    Vc=cell(nFrames,1);
else
    error('wrong face measure (second input variable)');
end

directionStringCell=cell(nStrains,1);
for is=1:nStrains
    switch faceMeasureCell{is}
        case 'Epc1'
            directionStringCell{is}='Epc1vecCur';
            optStruct.supTitleString{is}='1st principal Lagrangian strain';
        case 'Epc2'
            directionStringCell{is}='Epc2vecCur';
            optStruct.supTitleString{is}='2nd principal Lagrangian strain';
        case 'epc1'
            directionStringCell{is}='epc1vec';
            optStruct.supTitleString{is}='1st principal Eulerian strain';
        case 'epc2'
            directionStringCell{is}='epc2vec';
            optStruct.supTitleString{is}='2nd principal Eulerian strain';
        case 'Lamda1'
            directionStringCell{is}='Epc1vecCur';
            optStruct.supTitleString{is}='1st principal stretch';
        case 'Lamda2'
            directionStringCell{is}='Epc2vecCur';
            optStruct.supTitleString{is}='2nd principal stretch';
        otherwise
            error('unexpected face measure string. plots not created');
    end
    for it=1:nFrames
        if RBMlogic
            FC{it,is}=DIC3DPPresults.Deform_ARBM.(faceMeasureCell{is}){it}; % face color (strain)
            D{it,is}=DIC3DPPresults.Deform_ARBM.(directionStringCell{is}){it}; % direction (unit vector)

            switch faceMeasureCell{is}
                case {'Epc1','Epc2','epc1','epc2'}
                    Ds{it,is}=optStruct.quiverScaleFactor*FC{it,is}.*D{it,is}; % direction with magnitude (scaled vector)
                    DsLengths=sqrt((sum(Ds{it,is}.^2,2)));
                    LogicTooLong=DsLengths>meanEdgeLength;
                    Ds{it,is}(LogicTooLong,:)=meanEdgeLength*Ds{it,is}(LogicTooLong,:)./DsLengths(LogicTooLong);
                case {'Lamda1','Lamda2'}
                    Ds{it,is}=optStruct.quiverScaleFactor*(FC{it,is}-1).*D{it,is}; % direction with magnitude (scaled vector)
                    DsLengths=sqrt((sum(Ds{it,is}.^2,2)));
                    LogicTooLong=DsLengths>meanEdgeLength;
                    Ds{it,is}(LogicTooLong,:)=meanEdgeLength*Ds{it,is}(LogicTooLong,:)./DsLengths(LogicTooLong);
            end
            Vc{it,is}=DIC3DPPresults.FaceCentroids_ARBM{it}-.5*Ds{it,is};
        else
            FC{it,is}=DIC3DPPresults.Deform.(faceMeasureCell{is}){it}; % face color (strain)
            D{it,is}=DIC3DPPresults.Deform.(directionStringCell{is}){it}; % direction (unit vector)

            switch faceMeasureCell{is}
                case {'Epc1','Epc2','epc1','epc2'}
                    Ds{it,is}=optStruct.quiverScaleFactor*FC{it,is}.*D{it,is}; % direction with magnitude (scaled vector)
                    DsLengths=sqrt((sum(Ds{it,is}.^2,2)));
                    LogicTooLong=DsLengths>meanEdgeLength;
                    Ds{it,is}(LogicTooLong,:)=meanEdgeLength*Ds{it,is}(LogicTooLong,:)./DsLengths(LogicTooLong);
                case {'Lamda1','Lamda2'}
                    Ds{it,is}=optStruct.quiverScaleFactor*(FC{it,is}-1).*D{it,is}; % direction with magnitude (scaled vector)
                    DsLengths=sqrt((sum(Ds{it,is}.^2,2)));
                    LogicTooLong=DsLengths>meanEdgeLength;
                    Ds{it,is}(LogicTooLong,:)=meanEdgeLength*Ds{it,is}(LogicTooLong,:)./DsLengths(LogicTooLong);
            end
            Vc{it,is}=DIC3DPPresults.FaceCentroids{it}-.5*Ds{it,is};
        end

        if ~isempty(optStruct.maxCorrCoeff)
            corrNow=DIC3DPPresults.FaceCorrComb{it};
            FC{it,is}(corrNow>optStruct.maxCorrCoeff,:)=NaN;
            D{it,is}(corrNow>optStruct.maxCorrCoeff,:)=NaN;
            Vc{it,is}(corrNow>optStruct.maxCorrCoeff,:)=NaN;
        end

    end
end
FCmat = cell2mat(FC);
if ~isfield(optStruct,'colorBarLimits')
    switch faceMeasureCell{is}
        case {'Epc1','Epc2','epc1','epc2'}
            Emax=max(abs(FCmat(:)));
            optStruct.colorBarLimits=[-Emax Emax];
        case {'Lamda1','Lamda2'}
            Lmax=max(abs(FCmat(:)-1));
            if Lmax>1
                optStruct.colorBarLimits=[0 2];
            else
                optStruct.colorBarLimits=[1-prctile(abs(FCmat(:)-1),100) 1+prctile(abs(FCmat(:)-1),100)];
            end
    end
end
colorBarLogic=1;
if ~isfield(optStruct,'colorMap')
    optStruct.colorMap=.8*coldwarm;
end

%% Plot

animStruct=struct;

hf=cFigure;
hf.Units='normalized'; hf.OuterPosition=[.05 .05 .9 .9]; hf.Units='pixels';

for is=1:nStrains
    subplot(1,nStrains,is);

    axisGeom;
    ax=gca;
    ax.CameraUpVector=[0 0 optStruct.zDirection];
    colormap(optStruct.colorMap);
    if colorBarLogic
        colorbar;
        caxis(optStruct.colorBarLimits);
    end

    title(optStruct.supTitleString{is});
%     axis off
    % camlight headlight

    it=1;
    Fnow=DIC3DPPresults.Faces;
    if RBMlogic
        Pnow=DIC3DPPresults.Points3D_ARBM{it};
    else
        Pnow=DIC3DPPresults.Points3D{it};
    end

    Vnow=Vc{it,is};
    FCnow=FC{it,is};
    Dnow=Ds{it,is};
    if optStruct.smoothLogic
        smoothPar.lambda=0.5;
        smoothPar.n=2;
        [FCnow]=patchSmoothFaceMeasure(Fnow,Pnow,FCnow,smoothPar);
    end
    FCnow(FCnow<optStruct.dataLimits(1))=NaN;
    FCnow(FCnow>optStruct.dataLimits(2))=NaN;

    hp(is)=gpatch(Fnow,Pnow,FCnow,optStruct.lineColor,optStruct.FaceAlpha); hold on
    hq(is)=quiver3(Vnow(:,1),Vnow(:,2),Vnow(:,3),Dnow(:,1),Dnow(:,2),Dnow(:,3),0,'Color',.2*[1 1 1],'ShowArrowHead','off','AutoScale','off'); hold on;
%
    h_ax=gca;
    h_ax.XLim = xl; h_ax.YLim = yl; h_ax.ZLim = zl;

end

%% fill in the animstruct

animStruct.Time=1:nFrames;
animStruct.Handles=cell(1,nFrames);
animStruct.Props=cell(1,nFrames);
animStruct.Set=cell(1,nFrames);

for it=1:nFrames
    animStruct.Handles{it}=[];
    animStruct.Props{it}=cell(1,8*nStrains);
    animStruct.Set{it}=cell(1,8*nStrains);

    for is=1:nStrains
        if RBMlogic
            Pnow=DIC3DPPresults.Points3D_ARBM{it};
        else
            Pnow=DIC3DPPresults.Points3D{it};
        end
        Fnow=DIC3DPPresults.Faces;
        FCnow=FC{it,is};
        if optStruct.smoothLogic
            [FCnow]=patchSmoothFaceMeasure(Fnow,Pnow,FCnow,smoothPar);
        end
        FCnow(FCnow<optStruct.dataLimits(1))=NaN;
        FCnow(FCnow>optStruct.dataLimits(2))=NaN;

        Vnow=Vc{it,is};
        Dnow=Ds{it,is};

        animStruct.Handles{it}=[animStruct.Handles{it} hp(is) hp(is) hq(is) hq(is) hq(is) hq(is) hq(is) hq(is)]; %Handles of objects to animate (add one every pair)

        animStruct.Props{it}{1+8*(is-1)}='CData';
        animStruct.Props{it}{2+8*(is-1)}='Vertices'; %Properties of objects to animate
        animStruct.Props{it}{3+8*(is-1)}='XData'; %Properties of objects to animate
        animStruct.Props{it}{4+8*(is-1)}='YData'; %Properties of objects to animate
        animStruct.Props{it}{5+8*(is-1)}='ZData'; %Properties of objects to animate
        animStruct.Props{it}{6+8*(is-1)}='UData'; %Properties of objects to animate
        animStruct.Props{it}{7+8*(is-1)}='VData'; %Properties of objects to animate
        animStruct.Props{it}{8+8*(is-1)}='WData'; %Properties of objects to animate

        animStruct.Set{it}{1+8*(is-1)}=FCnow;
        animStruct.Set{it}{2+8*(is-1)}=Pnow; %Property values for to set in order to animate
        animStruct.Set{it}{3+8*(is-1)}=Vnow(:,1); %Property values for to set in order to animate
        animStruct.Set{it}{4+8*(is-1)}=Vnow(:,2); %Property values for to set in order to animate
        animStruct.Set{it}{5+8*(is-1)}=Vnow(:,3); %Property values for to set in order to animate
        animStruct.Set{it}{6+8*(is-1)}=Dnow(:,1); %Property values for to set in order to animate
        animStruct.Set{it}{7+8*(is-1)}=Dnow(:,2); %Property values for to set in order to animate
        animStruct.Set{it}{8+8*(is-1)}=Dnow(:,3); %Property values for to set in order to animate

        h_ax.XLim = xl; h_ax.YLim = yl; h_ax.ZLim = zl;

    end
end

anim8(hf,animStruct);

addColorbarLimitsButton(hf);
addColormapButton(hf);
addEdgeColorButton(hf);
addFaceAlphaButton(hf);
addLightButton(hf);
addAmbientStrengthButton(hf);
addDiffuseStrengthButton(hf);
addSpecularStrengthButton(hf);
addQuiverFactorButton(hf);
addFaceLightingButton(hf);

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
% %
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function U = create_guess(varargin)
%CREATE_GUESS Creates initial guess for CP or Tucker fitting.
%
%  U = CREATE_GUESS('Param',value,...) creates an initial guess at the
%  factor matrices for a CP or Tucker decomposition. The factors can be
%  generated randomly, random orthogonal, etc. If the tensor is provided,
%  it can be alternatively generated via the HO-SVD.
%
%   --- Parameters ---
%
%   'Factor_Generator' - Method to be used to generate the factor matrices.
%   Options:
%      - 'rand' (uniform on [0,1])
%      - 'randn' (standard normal distribution)
%      - 'orthogonal'
%      - 'stochastic' (uniform on [0,1] with column sums rescaled to 1)
%      - 'nvecs' (HOSVD solution)
%      - 'pertubation' of the true solution
%   Alternatively, pass in a function that accepts two arguments (the size
%   of the matrix) and generates the desired factor. Default: 'rand'
%
%   'Size' - Size of the tensor. Required to be specified unless 'Data' or
%   'Soln' is given. Default: []
%
%   'Num_Factors' - Number of factors (can be either a single value for CP
%   or a vector for Tucker). Required to be specified unless 'Soln' is
%   given. Default: []
%
%   'Data' - The actual tensor to be fit. Required if 'nvecs' is the
%   selected Factor Generator. The 'Size' parameter is ignored if this
%   is specified. Default: []
%
%   'Soln' - The actual solution to the problem. Required if 'pertubation'
%   is the selected Factor Generator. The 'Size' and 'Num_Factors'
%   parameters are ignored if this is specified. Default: []
%
%   'Pertubation' - Size of the pertubation is the 'pertubation' option is
%   selected under 'Factor_Generator'. The pertubation is of the form U+p*N
%   where U is the original factor matrix, N is a noise matrix with entries
%   selected for a standard normal distribution, and p is the pertubation
%   parameter times ||U||/||N||. Default: 0.10
%
%   'Skip' - Specifies mode to skip in initial guess generation (this is
%   useful for ALS). Default: 0 (no skipping)
%
%   'State' - State of the random number generator. This can be used
%   to reproduce results.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Random set-up
defaultStream = RandStream.getDefaultStream;

%% Parse inputs
p = inputParser;
p.addParamValue('Factor_Generator', 'rand', @(x) isa(x,'function_handle') || ...
    ismember(lower(x),{'rand','randn','orthogonal','stochastic','nvecs','pertubation'}));
p.addParamValue('Size', [], @(x) isempty(x) || all(x));
p.addParamValue('Num_Factors', [], @(x) isempty(x) || all(x));
p.addParamValue('Soln', [], @(x) isempty(x) || isa(x,'ktensor') || isa(x,'ttensor'));
p.addParamValue('Data', [], @(x) isempty(x) || isa(x,'tensor') || isa(x,'sptensor'));
p.addParamValue('Pertubation', 0.10, @(x) x >= 0 & x < 1);
p.addParamValue('Skip', 0);
p.addParamValue('State', defaultStream.State, @(x) true);
p.parse(varargin{:});
params = p.Results;

%% Initialize random number generator with specified state.
defaultStream.State = params.State;

%% Determine problem size
if ~isempty(params.Soln)
    sz = size(params.Soln);
elseif ~isempty(params.Data)
    sz = size(params.Data);
else
    sz = params.Size;
end
if isempty(sz)
    error('Size must be specified');
end
nd = length(sz);
modes = setdiff(1:nd,params.Skip);

%% Determine number of factors
if ~isempty(params.Soln)
    nf = zeros(nd,1);
    for n = 1:nd
        nf(n) = size(params.Soln.U{n},2);
    end
else
    nf = params.Num_Factors;
    if length(nf) == 1
        nf = nf * ones(nd,1);
    end
end

%% Create factor matrices
U = cell(nd,1);
if isa(params.Factor_Generator,'function_handle')
    for n = modes
        U{n} = params.Factor_Generator(sz(n), nf(n));
    end
    return;
end

switch(params.Factor_Generator)
    case 'rand'
        for n = modes
            U{n} = rand(sz(n), nf(n));
        end
    case 'randn'
        for n = modes
            U{n} = randn(sz(n), nf(n));
        end
    case 'orthogonal'
        for n = modes
            X = tt_RandOrthMat(sz(n));
            U{n} = X(:,1:nf(n));
        end
    case 'stochastic'
        for n = modes
            X = rand(sz(n), nf(n));
            S = sum(X,1);
            U{n} = X * diag(1./S);
        end
    case 'nvecs'
        if isempty(params.Data)
            error('Data required for nvecs initialization');
        end
        for n = modes
            U{n} = nvecs(params.Data,n,nf(n));
        end
    case 'pertubation'
        if isempty(params.Soln)
            error('Soln required for pertubation initialization');
        end
        for n = modes
            X = params.Soln{n};
            N = rand(size(X));
            p = params.Pertubation * norm(X,'fro') / norm(N,'fro');
            U{n} = X + p * N;
        end
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function polygons = splitPolygons(polygon)
%SPLITPOLYGONS Convert a NaN separated polygon list to a cell array of polygons.
%
%   POLYGONS = splitPolygons(POLYGON);
%   POLYGON is a N-by-2 array of points, possibly with pairs of NaN values.
%   The functions separates each component separated by NaN values, and
%   returns a cell array of polygons.
%
%   See also:
%   polygons2d
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2007-10-12,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

if iscell(polygon)
    % case of a cell array
    polygons = polygon;

elseif sum(isnan(polygon(:))) == 0
    % single polygon -> no break
    polygons = {polygon};

else
    % find indices of NaN couples
    inds = find(sum(isnan(polygon), 2) > 0);

    % number of polygons
    N = length(inds) + 1;
    polygons = cell(N, 1);

    % iterate over NaN-separated regions to create new polygon
    inds = [0; inds; size(polygon, 1)+1];
    for i = 1:N
        polygons{i} = polygon((inds(i)+1):(inds(i+1)-1), :);
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [xlat, alt] = geodet2(decl, rmag)

% geocentric to geodetic coordinates
% exact solution (Borkowski, 1989)

% input

%  decl = geocentric declination (radians)
%  rmag = geocentric distance (kilometers)

% output

%  xlat = geodetic latitude (radians)
%  alt  = geodetic altitude (kilometers)

% Orbital Mechanics with Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global req flat

fr = 1 / flat;

% determine x and z components of the geocentric distance

rx = rmag * cos(decl);

rz = rmag * sin(decl);

% compute geodetic latitude and altitude

if (rz == 0)
   % special case - equatorial

   xlat = 0;
   alt = rmag - req;

elseif (abs(decl) == 0.5 * pi)
   % special case - polar

   xlat = decl;
   alt = rmag - (req - req/fr);

else
   % general case

   b = sign(rz) * (req - req/fr);

   e = ((rz + b) * b/req - req)/rx;

   f = ((rz - b) * b/req + req)/rx;

   p = (e * f + 1) * 4/3;

   q = (e * e - f * f) * 2;

   d = p * p * p + q * q;

   if (d >= 0.0)
      s = sqrt(d) + q;
      s = sign(s) * (exp(log(abs(s))/3));
      v = p/s - s;
      v = -(q + q + v * v * v)/(3 * p);
   else
      v = 2 * sqrt(-p) * cos(acos(q/p/sqrt(-p))/3);
   end

   g = 0.5 * (e + sqrt(e * e + v));

   t = sqrt(g * g + (f - v * g)/(g + g - e)) - g;

   xlat = atan((1 - t * t) * req/(2 * b * t));

   alt = (rx - req * t) * cos(xlat) + (rz - b) * sin(xlat);
end


text:
%  Generating a Givens' (2x2) matrix T such that
%            T*v = [d; 0]
%
%   Input:  v -- a vector of dimension 2
%
%  Output:  T -- Givens' matrix
%           d -- the 2-norm of v
%
%  syntax  >> [T,d] = GivensMatrix(v)

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function out = apply_HLDA(C, X)
% APPLY_HLDA - Hierarchical linear discriminant analysis
%
%Synopsis:
%   out = apply_HLDA(C, X)
%
%Arguments:
%   C: STRUCT           - a hierarchical LDA classifier structure that
%                           contains the individual segment classifiers and
%                           the top-level classifier. Must include the
%                           fields 'seg' and 'final'.
%   X: DOUBLE [TxNxM] - Data matrix, with T temporal features, N
%                           N spatial features, and M
%                           training points/examples.
%Returns:
%   out: FLOAT[]        - an array containing the classifier score for each
%                           sample
%Description:
%   APPLY_HLDA applies a hierarchical LDA classifier given data and a
%   trained HLDA 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:
%   apply_HLDA(C, X)
%
%See also:
% train_HLDA

% validate argument types
misc_checkType(C, 'STRUCT(seg final)');
misc_checkType(X, 'DOUBLE');
if isfield(C,'nChannels')
  misc_checkType(X, 'DOUBLE[- -]');
end

nSegments = length(C.seg);

dims = size(X);

% make data matrix 3D if it has been tranformed to 2D for crossvalidation
if (length(size(X)) == 2) && isfield(C, 'nChannels')
  X = reshape(X, [], C.nChannels, dims(2));
end

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

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

for i = 1:nSegments
    seg = X(seg_idx(i) + 1 : seg_idx(i + 1), :, :); % i-th segment of X
    %apply LDA classifier for this segment
    seg = reshape(seg, [], dims(end));
    seg_scores(i,:) = apply_separatingHyperplane(C.seg(i), seg); %get classifier scores for segment
end

if isfield(C.final, 'B')
    %logistic regression as top-level classifier
%    exps = C.final.B(1) + C.final.B(2:end)'*seg_scores;
%    out = 1./(1 + exp(-exps));

    pihat = mnrval(C.final.B, seg_scores');

    pihat = pihat';
    out = pihat(1,:);
    %transform scores to interval [-1 1]
    out = -2*(out - 0.5);
else
    %LDA as top-level classifier
    out = apply_separatingHyperplane(C.final, seg_scores);
end

author: bbci


text:
function  write2hdf5_test_sequence
% write image data (and simple feature data) from all sequences to hdf5 file
% positive and negative example mixed

is_sp = 1; % is_sp = 1: use simple feature
sp_dim = 7;

data_dir = '/scratch/BS/pool2/tang/multicut/pairwise_example/';
% save dir
if ~is_sp
    save_dir = '/BS/multicut_tracking/work/Dataset/multicut/pairwise_example/tud_crossing_sequences/mixed_hdf5/';
else
    save_dir = '/BS/multicut_tracking/work/Dataset/multicut/pairwise_example/tud_crossing_sequences/mixed_sp_hdf5/';
end
if ~exist(save_dir)
    mkdir(save_dir);
end
% data set
data_set = {'tud-crossing'};
frame_rate = [25];

% load images for each sequence
num_sequeces = length(data_set);
positive_img_dir = cell(num_sequeces,1);
negative_img_dir = cell(num_sequeces,1);

positive_img_list = cell(num_sequeces,1);
negative_img_list = cell(num_sequeces,1);

rep_positive_dir = cell(num_sequeces,1);
rep_negative_dir = cell(num_sequeces,1);

rep_positive_label = cell(num_sequeces,1);
rep_negative_label = cell(num_sequeces,1);

sp_positive = cell(num_sequeces,1);
sp_positive_name = cell(num_sequeces,1);

sp_negative = cell(num_sequeces,1);
sp_negative_name = cell(num_sequeces,1);
for s = 1:num_sequeces
    positive_img_dir{s} = [data_dir data_set{s} '/positive/'];
    negative_img_dir{s} = [data_dir data_set{s} '/negative/'];

    fprintf('Listing all the images for sequence %d...\n',s);

    % positive examples
    tmp = dir([positive_img_dir{s} '*.jpg']);
    positive_img_list{s} = {tmp.name}';

    rep_positive_dir{s} = repmat(s, length(positive_img_list{s}),1);
    rep_positive_label{s} = ones(length(positive_img_list{s}),1);
    if is_sp ==1
        load([positive_img_dir{s} 'simple_feature.mat']);
        sp_positive{s} = [simple_feature simple_feature];
        sp_positive_name{s} = [simple_feature_name simple_feature_name];
    end

    % negative examples
    tmp = dir([negative_img_dir{s} '*.jpg']);
    negative_img_list{s} = {tmp.name}';

    rep_negative_dir{s} = repmat(s, length(negative_img_list{s}),1);
    rep_negative_label{s} = zeros(length(negative_img_list{s}),1);
    if is_sp ==1
        load([negative_img_dir{s} 'simple_feature.mat']);
        sp_negative{s} = [simple_feature simple_feature];
        sp_negative_name{s} = [simple_feature_name simple_feature_name];
    end
end

positive_img_list_all = cat(1,positive_img_list {:});
rep_positive_dir = cat(1,rep_positive_dir {:});
rep_positive_label = cat(1,rep_positive_label {:});

negative_img_list_all = cat(1,negative_img_list {:});
rep_negative_dir = cat(1,rep_negative_dir {:});
rep_negative_label = cat(1,rep_negative_label {:});

% image name list
img_list = [positive_img_list_all; negative_img_list_all];
% image dir list
rep_dir = [rep_positive_dir; rep_negative_dir];
% image label list
rep_label = [rep_positive_label; rep_negative_label];
% simple feature list
if is_sp == 1
    sp_positive = cat(2,sp_positive {:});
    sp_negative = cat(2,sp_negative {:});
    sp = [sp_positive';sp_negative'];
else
    sp = [];
end

num_img = length(positive_img_list_all) + length(negative_img_list_all);

% randomize
rand_list = randi([1,num_img],num_img,1);rand_list(end:num_img) =1; % make sure rand_list length
img_list = img_list(rand_list);
rep_dir = rep_dir(rand_list);
rep_label = rep_label(rand_list);
if is_sp ==1
    sp = sp(rand_list);
else
    sp = [];
end

% chunk size
chunksz = 1000;
% number of files
num_round = (num_img)/chunksz;
for k = 1:num_round
    fprintf('Saving hdf5 file %d out of %d ...\n', k, num_round);
    % must be single format
    batchdata = single(zeros(227,227,3,chunksz));

    % if no sp diminsion, all the labels are zeros.
    batchlabs = single(zeros(1+sp_dim,chunksz));

    % the end idx of last round
    last_end_img_idx = (k-1)*chunksz;
    % each example
    for i = 1:chunksz
        idx = last_end_img_idx + i;

        % current dataset label
        cur_dataset_label = rep_dir(idx);

        % current image name
        cur_image_name = img_list{idx};

        % current example's label
        cur_class_label = rep_label(idx);

        % current example's directory
        if cur_class_label
            cur_img_dir = positive_img_dir{cur_dataset_label};
        else
            cur_img_dir = negative_img_dir{cur_dataset_label};
        end
        % read example
        img = imread([cur_img_dir cur_image_name]);
        % unit to single
        img = imresize(img,[227,227]);
        % change to single format
        img =  single(img);
        % RGB to BGR
        img = img(:,:,[3 2 1]);
        % make row first
        img = permute(img, [2 1 3]);
        %
        % load to batchdata and label
        batchdata(:,:,:,i) = img;
        batchlabs(1,i) = cur_class_label;
        %
        if is_sp == 1
            cur_sp = sp{idx};
            cur_frame_rate = frame_rate(cur_dataset_label);
            % frame rate normalize ratio
            fr_norm = 25/cur_frame_rate;
            feature = compute_simple_feature(cur_sp,fr_norm);
            feature = single(feature);
            batchlabs(2:end,i) = feature;
        end
    end
    % save
    filename = [save_dir sprintf('%06d',last_end_img_idx + 1) '.h5'];
    store2hdf5(filename, batchdata, batchlabs, 1);
end
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
classdef DesignVarMonitor_LevelSet_3D < DesignVarMonitor_LevelSet

    properties (Access = protected)
        unfittedType = 'BOUNDARY'
        meshIncludeBoxContour = true
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan