text:
function [P]=tesSmooth_LAP(TES,V,IND_V,cPar)

%% CONTROL PARAMETERS

if isfield(cPar,'LambdaSmooth')
    LambdaSmooth=cPar.LambdaSmooth;
else
    LambdaSmooth=0.5; %DEFAULT
end

if isfield(cPar,'n')
    nMax=cPar.n;
else
    nMax=1; %DEFAULT
end

if isfield(cPar,'RigidConstraints')
    indRigid=cPar.RigidConstraints;
else
    indRigid=[]; %DEFAULT
end

if isfield(cPar,'Tolerance')
    SSQD_Tol=cPar.Tolerance;
else
    SSQD_Tol=[]; %DEFAULT
end

%%

if ~isempty(SSQD_Tol)
    SSQD_old=[];
    SSQD_ratio=0;
end

if isempty(IND_V)
    [~,IND_V]=patchIND(TES,V,2);
end
logicValid=IND_V>0;
indNoneValid=find(sum(logicValid,2)==0);

%%

nDims=size(V,2); %Number of dimensions

P=V;
PP=V;
Q=V;
for qIter=1:nMax

    %% SIMPLE LAPLACIAN SMOOTHENING

    %Loop for all dimensions
    for qDim=1:1:nDims
        Xp=NaN(size(IND_V,1),size(IND_V,2));
        Xp(logicValid)=P(IND_V(logicValid),qDim);
        Xp=gnanmean(Xp,2);
        Xp(indNoneValid)=V(indNoneValid,qDim);
        PP(:,qDim)=Xp;
    end
    P=P+LambdaSmooth.*(PP-P);

    %%

    %Put back constrained points
    if ~isempty(indRigid)
       P(indRigid,:)=V(indRigid,:);
    end

    if ~isempty(SSQD_Tol)
        %Compute sum of squared differences with respect to previous iteration
        SSQD_new=gnansum((P(:)-Q(:)).^2);
        if ~isempty(SSQD_old)
            SSQD_ratio=SSQD_new./SSQD_old;
        end

        %Store current metrics
        Q=P;
        SSQD_old=SSQD_new;
        if abs(1-SSQD_ratio)<=SSQD_Tol
           break %STOP SMOOTHING LOOP IF TOLERANCE IS REACHED
        end
    end

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function gopt = illustration_options(gdefault, varargin)
%ILLUSTRATION_OPTIONS Options for illustrations of EQ partitions
%
%Syntax
% gopt = illustration_options(gdefault,options);
%
%Description
% GOPT = ILLUSTRATION_OPTIONS(GDEFAULT,options) collects illustration options,
% specified as name, value pairs, and places these into the structure GOPT.
% The structure GDEFAULT is used to define default option values.
%
% The structures gdefault and gopt may contain the following fields:
% fontsize:      numeric
% long_title:    boolean
% stereo:        boolean
% show_points:   boolean
% show_sphere:   boolean
% show_surfaces: boolean
%
% The following illustration options are available.
%
% 'fontsize':    Font size used in titles.
%      number    Assigns number to field gopt.fontsize.
%
% 'title':       Length of titles.
%     'long':    Long titles.
%                Sets gopt.show_title to true.
%                Sets gopt.long_title to true.
%     'short':   Short titles.
%                Sets gopt.show_title to true.
%                Sets gopt.long_title to false.
%     'none':    No titles.
%                Sets gopt.show_title to false.
%                Sets gopt.long_title to false.
%     'show':    Show default titles.
%                Sets gopt.show_title to true.
%     'hide':    Same as 'none'.
%
% 'proj':        Projection from the sphere to the plane R^2 or the space R^3.
%     'stereo':  Stereographic projection from the sphere to the whole space.
%                Sets gopt.stereo to true.
%     'eqarea':  Equal area projection from the sphere to the disk or ball.
%                Sets gopt.stereo to false.
%
% 'points':      Show or hide center points of regions.
%     'show':    Show center points of regions.
%                Sets gopt.show_points to true.
%     'hide':    Hide center points of regions.
%                Sets gopt.show_points to false.
%
% 'sphere':      Show or hide the sphere S^2.
%     'show':    Show sphere.
%                Sets gopt.show_sphere to true.
%     'hide':    Hide sphere.
%                Sets gopt.show_sphere to false.
%
% 'surf':        Show or hide surfaces of regions of a partition of S^3.
%     'show':    Show surfaces of regions.
%                Sets gopt.show_surfaces to true.
%     'hide':    Hide surfaces of regions.
%                Sets gopt.show_surfaces to false.
%
%Examples
% > gdefault.fontsize=14;
% > gopt=illustration_options(gdefault,'proj','stereo')
% gopt =
%     fontsize: 14
%       stereo: 1
%
% > gopt=illustration_options(gdefault,'proj','stereo','fontsize',12)
% gopt =
%    fontsize: 12
%      stereo: 1

% Copyright 2004-2005 Paul Leopardi for the University of New South Wales.
% $Revision 1.10 $ $Date 2005-06-01 $
% Documentation files renamed
% $Revision 1.00 $ $Date 2005-02-12 $
%
% For licensing, see COPYING.
% For references, see AUTHORS.
% For revision history, see CHANGELOG.

gopt = gdefault;
nargs = length(varargin);
nopts = floor(nargs/2);
opt_args = {varargin{1:2:2*nopts-1}};
for k=1:nopts
    if ~ischar([opt_args{k}])
        fprintf('Option names must be character strings\n');
        option_error(varargin{:});
    end
end
opt_vals = {varargin{2:2:2*nopts}};

option_name = 'fontsize';
pos = strmatch(option_name,opt_args,'exact');
if ~isempty(pos)
    if (length(pos) == 1)
        gopt.fontsize = opt_vals{pos};
    else
        duplicate_error(option_name,varargin{:});
    end
end

option_name = 'title';
pos = strmatch(option_name,opt_args,'exact');
if ~isempty(pos)
    if (length(pos) == 1)
        value = opt_vals{pos};
    else
        duplicate_error(option_name,varargin{:});
    end
    switch value
    case 'long'
        gopt.show_title = true;
        gopt.long_title = true;
    case 'short'
        gopt.show_title = true;
        gopt.long_title = false;
    case 'none'
        gopt.show_title = false;
        gopt.long_title = false;
    case 'hide'
        gopt.show_title = false;
    case 'hide'
        gopt.show_title = false;
        gopt.long_title = false;
    case 'show'
        gopt.show_title = true;
    otherwise
        value_error(value,varargin{:});
    end
end

option_name = 'proj';
pos = strmatch(option_name,opt_args);
if ~isempty(pos)
    if (length(pos) == 1)
        value = opt_vals{pos};
    else
        duplicate_error(option_name,varargin{:});
    end
    switch value
    case 'stereo'
        gopt.stereo = true;
    case 'eqarea'
        gopt.stereo = false;
    otherwise
        value_error(value,varargin{:});
    end
end

option_name = 'points';
pos = strmatch(option_name,opt_args,'exact');
if ~isempty(pos)
    if (length(pos) == 1)
        value = opt_vals{pos};
    else
        duplicate_error(option_name,varargin{:});
    end
    switch value
    case 'show'
        gopt.show_points = true;
    case 'hide'
        gopt.show_points = false;
    otherwise
        value_error(value,varargin{:});
    end
end

option_name = 'surf';
pos = strmatch(option_name,opt_args);
if ~isempty(pos)
    if (length(pos) == 1)
        value = opt_vals{pos};
    else
        duplicate_error(option_name,varargin{:});
    end
    switch value
    case 'show'
        gopt.show_surfaces = true;
    case 'hide'
        gopt.show_surfaces = false;
    otherwise
        value_error(value,varargin{:});
    end
end

option_name = 'sphere';
pos = strmatch(option_name,opt_args);
if ~isempty(pos)
    if (length(pos) == 1)
        value = opt_vals{pos};
    else
        duplicate_error(option_name,varargin{:});
    end
    switch value
    case 'show'
        gopt.show_sphere = true;
    case 'hide'
        gopt.show_sphere = false;
    otherwise
        value_error(value,varargin{:});
    end
end
%
% end function

function duplicate_error(option_name,varargin)
fprintf('Duplicate option %s\n',option_name);
option_error(varargin{:});
%
% end function

function value_error(value,varargin)
fprintf('Invalid option value ');
disp(value);
option_error(varargin{:});
%
% end function

function option_error(varargin)
fprintf('Error in options:\n');
disp(varargin);
error('Please check "help illustration_options" for options');
%
% end function

repo: gspbox


text:
% ITK_PSTRANSFORM: Spatial transformation (i.e. warp) of a set of points,
% defined from a known landmark correspondence
%
% This MEX-function provides a Matlab interface to run ITK's Kernel
% and B-spline transforms:
%
%   itk::ElasticBodySplineKernelTransform
%   itk::ElasticBodyReciprocalSplineKernelTransform
%   itk::ThinPlateSplineKernelTransform
%   itk::ThinPlateR2LogRSplineKernelTransform
%   itk::VolumeSplineKernelTransform
%   itk::BSplineScatteredDataPointSetToImageFilter
%
%
% YI = itk_pstransform(TRANSFORM, X, Y, XI)
%
%   X, Y are 2-column (2D) or 3-column (3D) matrices. Each row has
%   the coordinates of a point. The warp is defined so that
%   X(i,:)->Y(i,:).
%
%   XI is a matrix with the same number of columns as X, Y. Each row
%   has the coordinates of a point to be warped.
%
%   YI has the same dimensions as XI. YI contains the coordinates of
%   the warped points.
%
%   TRANSFORM is a string that allows to select the type of warp (no
%   defaults):
%
% YI = itk_pstransform('elastic', X, Y, XI)
% YI = itk_pstransform('elasticr', X, Y, XI)
% YI = itk_pstransform('tps', X, Y, XI)
% YI = itk_pstransform('tpsr2', X, Y, XI)
% YI = itk_pstransform('volume', X, Y, XI)
%
%   'elastic':  itk::ElasticBodySplineKernelTransform
%   'elasticr': itk::ElasticBodyReciprocalSplineKernelTransform
%   'tps':      itk::ThinPlateSplineKernelTransform
%   'tpsr2':    itk::ThinPlateR2LogRSplineKernelTransform
%   'volume':   itk::VolumeSplineKernelTransform
%
%   Note that 'tpsr2' produces the same result as our Matlab implementation
%   pts_tps_map(), as it implements the classic kernel proposed by
%   Bookstein, r^2 ln(r^2).
%
% YI = itk_pstransform('bspline', X, Y, XI, ORDER, LEVELS)
%
%   'bspline':  itk::BSplineScatteredDataPointSetToImageFilter
%
%   By Nicholas J. Tustison, James C. Gee in the Insight Journal
%   paper: http://hdl.handle.net/1926/140
%
%   ORDER is an integer with the B-spline order. By default, ORDER=3,
%   and the B-spline is cubic.
%
%   LEVELS is an integer with the number of multi-resolution levels
%   in the algorithm. A higher number of levels will make the spline
%   more flexible and match the landmarks better. By default, LEVELS=5.
%
% See also: pts_tps_map, pts_tps_weights.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2011 University of Oxford
% Version: 0.1.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

error('MEX file not found')

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [cnn_feature1, cnn_feature2] = FeatureMapNormalization(cnn_feature)

r = size(cnn_feature,1);
c = size(cnn_feature,2);
f = size(cnn_feature,3);
t = size(cnn_feature,4);
cnn_feature1 = permute(cnn_feature,[1,2,4,3]);
cnn_feature1 = reshape(cnn_feature1,r*c*t,[]);
cnn_feature1 = bsxfun(@rdivide,cnn_feature1,max(cnn_feature1,[],1)+eps);
cnn_feature1 = reshape(cnn_feature1,r,c,t,f);
cnn_feature1 = permute(cnn_feature1,[1,2,4,3]);

cnn_feature2 = bsxfun(@rdivide,cnn_feature,max(cnn_feature,[],3)+eps);

end

author: wanglimin

repo: TDD

save_path: github-repos/MATLAB/wanglimin-TDD


text:
% data,labels are the training data and labels.
% info has the information about these examples (expi,flies,t)
% newdata,newlabels are the data and labels for examples not part of training data.
% newinfo has information about these examples.
% expdirs have the expdirs, trx has the loaded tracks for that experiment, and moviefilename is 'movie.ufmf';

%% Bouts.
bouts = struct('ndx',[],'label',[],'timestamp',[]);
for expNdx = 1:obj.nexps
  for flyNdx = 1:obj.nflies_per_exp(expNdx)
    curLabels = obj.GetLabels(expNdx,flyNdx);
    for boutNum = 1:numel(curLabels.t0s)
      idx =  obj.FlyNdx(expNdx,flyNdx) & ...
        obj.windowdata.t >= curLabels.t0s(boutNum) & ...
        obj.windowdata.t < curLabels.t1s(boutNum);
      bouts.ndx(end+1,:) = obj.FlyNdx(expNdx,flyNdx) & ...
        obj.windowdata.t >= curLabels.t0s(boutNum) & ...
        obj.windowdata.t < curLabels.t1s(boutNum);
      bouts.label(end+1) = find(strcmp(obj.labelnames,curLabels.names{boutNum}));
      bouts.timestamp(end+1) = curLabels.timestamp(boutNum);
    end

  end
end

zz = bouts.timestamp;

moviefilename = 'movie.ufmf';

clparams = struct('iter',100,'iter_updates',10,...
      'numSample',2500,'numBins',30,'CVfolds',7,...
      'baseClassifierTypes',{'Decision Stumps'},'baseClassifierSelected',1);

allbagModels = {}; celldistmat = {};

%% To show the images.. load the trx.
for ndx = 1:obj.nexps
  trxfile = fullfile(obj.expdirs{ndx},obj.GetFileName('trx'));
  trx{ndx} = load_tracks(trxfile);
end
expdirs = obj.expdirs;
%%

for outqndx = 1:numel(Q.traintimes_chaseactive)-1
%%
  qq(1) = Q.traintimes_chaseactive(outqndx);
  if outqndx+1>numel(Q.traintimes_chaseactive)
    qq(2) = Q.traintimes_chaseactive(end);
  else
    qq(2) = Q.traintimes_chaseactive(outqndx+1);
  end

%%

oldbouts = find(zz<qq(1));
oldidx = false(size(bouts.ndx(1,:)));
oldidxall = false(size(bouts.ndx(1,:)));
for ndx = oldbouts(:)'
  nn = find(bouts.ndx(ndx,:),1);
  xx = find(bouts.ndx(ndx,:),1,'last');
  ii = round( (nn+xx)/2);
  oldidx(ii) = true;
  oldidxall = oldidxall | bouts.ndx(ndx,:);
end

oldndxintoall = find(oldidx(oldidxall));
data = obj.windowdata.X(oldidxall,:);
labels = obj.windowdata.labelidx_cur(oldidxall);
info = struct;
info.expi = obj.windowdata.exp(oldidxall);
info.flies = obj.windowdata.flies(oldidxall);
info.t = obj.windowdata.t(oldidxall);

newbouts = find(zz>=qq(1) & zz<qq(2) );
newidx = false(size(bouts.ndx(1,:)));
newidxall = false(size(bouts.ndx(1,:)));
for ndx = newbouts(:)'
  nn = find(bouts.ndx(ndx,:),1);
  xx = find(bouts.ndx(ndx,:),1,'last');
  ii = round( (nn+xx)/2);
  newidx(ii) = true;
  newidxall = newidxall | bouts.ndx(ndx,:);
end

newdata = obj.windowdata.X(newidxall,:);
newlabels = obj.windowdata.labelidx_cur(newidxall);
newinfo = struct;
newinfo.expi = obj.windowdata.exp(newidxall);
newinfo.flies = obj.windowdata.flies(newidxall);
newinfo.t = obj.windowdata.t(newidxall);
fprintf('Training,pos:%d,neg:%d Test:pos%d,neg:%d\n',nnz(labels==1),...
  nnz(labels==2),nnz(newlabels==1),nnz(newlabels==2));

%%
binVals = findThresholds(data,clparams);
bins = findThresholdBins(data,binVals);

%%
%{
[bagModels distmat] = doBagging(data,labels,[],binVals,bins,clparams);
allbagModels{outqndx} = bagModels;
celldistmat{outqndx} = distmat;
%}

%%
bagModels = allbagModels{outqndx};
distmat = celldistmat{outqndx};

%%
newdistmat = findBagCoords(newdata,bagModels);

alldistmat = [distmat; newdistmat];
alllabels = [labels(:)' (newlabels(:)'+2)];

numex = size(alldistmat,1);
xdistmat = zeros(numex);
for ndx = 1:numex
    tt = repmat(alldistmat(ndx,:),numex,1);
    xdistmat(:,ndx) = nanmean(abs(tt-alldistmat),2);
end

wt = getWeights(sign([labels(:)' newlabels(:)']-1.5));
wt = sqrt(wt*wt');
% wt = ones(size(xdistmat));
Z = mdscale(xdistmat,2,'Start','random','Weights',wt);
allnlabels = [labels(:)' newlabels(:)'];
R(1) = sum(Z(1:numel(labels),1).*sign(labels-1.5));
R(2) = sum(Z(1:numel(labels),2).*sign(labels-1.5));
Z = Z(:,1:2);
theta = -atan2(R(2),R(1));
rotmat = [cos(theta) sin(theta); -sin(theta) cos(theta)];
Zr = Z*rotmat;
Zold = Z;
Z = Zr;
Z(:,1) = (Z(:,1)-min(Z(:,1)))/(max(Z(:,1))-min(Z(:,1)));
Z(:,2) = (Z(:,2)-min(Z(:,2)))/(max(Z(:,2))-min(Z(:,2)));

%%
%{
end
save embedBagsChase allbagModels celldistmat;
%}
%% Scatter plot the embedded points.

wt = getWeights(sign(labels-1.5));
[~,bb] = loglossboostLearnRandomFeatures(data,sign(labels-1.5),100,wt,binVals,bins,clparams);
allpred = myBoostClassify([data;newdata],bb);
tr = prctile(abs(allpred),80);
allnlabels = [labels(:)' newlabels(:)'];

%%
C = zeros(numel(alllabels),3);
% S = max(10,1./(1+exp(allpred.*(sign(allnlabels'-1.5))))*100);
S = [];
for ndx = 1:numel(alllabels)
  switch alllabels(ndx)
    case 1
      C(ndx,:) = [1 0 0];
    case 2
      C(ndx,:) = [0 0 1];

    case 3
      C(ndx,:) = [1 0.6 0];
    case 4
      C(ndx,:) = [0 0.6 1];
  end

%  S(ndx) = 1/(min(1,abs(allpred(ndx)/tr)))*10;

%   if sign(allpred(ndx))== sign(allnlabels(ndx)-1.5)
%     S(ndx) = 20;
%   else
%     S(ndx) = 100;
%   end
end
% Scale by mfac so that the image patches are more separate

%%
%{
mfac = 6000;
figure; f = scatter(Z(:,1)*mfac,Z(:,2)*mfac,S,C,'filled');

%% Connect frames from same bouts

hold on;
curZ = Z(1:numel(labels),:);
for ndx = 1:numel(bouts.timestamp)

  if bouts.timestamp(ndx)>qq(1); continue; end;
  curidx = bouts.ndx(ndx,:);
  lndx = find(curidx(oldidxall));
  if bouts.label(ndx)==1
    col = [1 0 0];
  else
    col = [0 0 1];
  end
  plot(curZ(lndx,1)*mfac,curZ(lndx,2)*mfac,'Color',col);
end

curZ = Z(numel(labels)+1:end,:);
for ndx = 1:numel(bouts.timestamp)
  if bouts.timestamp(ndx)>qq(2) || bouts.timestamp(ndx)<qq(1); continue; end;
  curidx = bouts.ndx(ndx,:);
  lndx = find(curidx(newidxall));
  if bouts.label(ndx)==1
    col = [1 0.6 0];
  else
    col = [0 0.6 1];
  end
  plot(curZ(lndx,1)*mfac,curZ(lndx,2)*mfac,'Color',col);
end

%% To show the images.. load the trx.

for ndx = 1:obj.nexps
  trxfile = fullfile(obj.expdirs{ndx},obj.GetFileName('trx'));
  trx{ndx} = load_tracks(trxfile);
end

%% Load images and tracks for training data.
im = {}; traj = {};
for ndx = 1:numel(info.expi)
  expi = info.expi(ndx);
  flies = info.flies(ndx);
  t = info.t(ndx);
  mname = fullfile(expdirs{expi},moviefilename);
  [im{ndx},traj{ndx}] = vid2im(mname,flies,t,trx{expi});
end

%% Load images and tracks for test data.
newim = {}; newtraj = {};
for ndx = 1:numel(newinfo.expi)
  expi = newinfo.expi(ndx);
  flies = newinfo.flies(ndx);
  t = newinfo.t(ndx);
  mname = fullfile(expdirs{expi},moviefilename);
  [newim{ndx},newtraj{ndx}] = vid2im(mname,flies,t,trx{expi});
end

%% Show the most negative images on out of training bouts

hold on;

curZ = Z(numel(labels)+1:end,:);
for ndx = 1:numel(bouts.timestamp)
  if bouts.timestamp(ndx)>qq(2) || bouts.timestamp(ndx)<qq(1); continue; end;
  curidx = bouts.ndx(ndx,:);
  lndx = find(curidx(newidxall));
  curScores = myBoostClassify(newdata(lndx,:),bb);
  [~,minndx] = min(curScores*sign(bouts.label(ndx)-1.5));
  curlndx = lndx(minndx);
  ii = uint8(repmat(newim{curlndx},[1 1 3]));
  newX = curZ(curlndx,1)*mfac; newY = curZ(curlndx,2)*mfac;
  image(newX,newY,ii);
  plot(newtraj{curlndx}.Y+newX,newtraj{curlndx}.X+newY)
end
axis equal;

%% Show the most negative images on  training bouts

hold on;

curZ = Z(1:numel(labels),:);
for ndx = 1:numel(bouts.timestamp)
  if bouts.timestamp(ndx)>qq(1), continue; end;
  curidx = bouts.ndx(ndx,:);
  lndx = find(curidx(oldidxall));
  curScores = myBoostClassify(data(lndx,:),bb);
  [~,minndx] = min(curScores*sign(bouts.label(ndx)-1.5));
  curlndx = lndx(minndx);
  ii = uint8(repmat(im{curlndx},[1 1 3]));
  newX = curZ(curlndx,1)*mfac; newY = curZ(curlndx,2)*mfac;
  image(newX,newY,ii);
  plot(traj{curlndx}.Y+newX,traj{curlndx}.X+newY)
end
axis equal;

%% Show the training images.

hold on;
curZ = Z(1:size(distmat,1),:);
for ndx = 1:numel(im)
  ii = uint8(repmat(im{ndx},[1 1 3]));
  newX = curZ(ndx,1)*mfac; newY = curZ(ndx,2)*mfac;
  image(newX,newY,ii);
  plot(traj{ndx}.Y+newX,traj{ndx}.X+newY)
end
axis equal;

%% Show the "out of training" images

hold on;
curZ = Z( (size(distmat,1)+1):end,:);
for ndx = 1:numel(newim)
  ii = uint8(repmat(newim{ndx},[1 1 3]));
  newX = curZ(ndx,1)*mfac; newY = curZ(ndx,2)*mfac;
  image(newX,newY,ii);
  plot(newtraj{ndx}.Y+newX,newtraj{ndx}.X+newY)
end
axis equal;
%
%%
mfac = 1200; sfac = 30;
clear Zm;
Zm(:,1) = Z(:,1)-min(Z(:,1));
Zm(:,2) = Z(:,2)-min(Z(:,2));

zszx = max(Zm(:,1))*mfac/sfac;
zszy = max(Zm(:,2))*mfac/sfac;
bimg = ones(round(zszy)+1,round(zszx)+1,3);

counts = zeros(round(zszy)+1,round(zszx)+1);
for ndx= 1:numel(labels)
  curx = round( Zm(ndx,1) *mfac/sfac)+1;
  cury = round( Zm(ndx,2) *mfac/sfac)+1;
  bimg(cury,curx,:) = reshape(C(ndx,:),[1 1 3]);
end

bbimg = imfilter(bimg,fspecial('gaussian',2*3*8,8),'symmetric');
hsvbimg = rgb2hsv(bbimg);
satimg = hsvbimg(:,:,2);
satimg = satimg-mean(satimg(:))+0.1;
satimg(satimg>1) = 1;
hsvbimg(:,:,2) = satimg;
bbimg = hsv2rgb(hsvbimg);
fbimg = imresize(bbimg,sfac);

fbimg(fbimg(:)>1) = 1;

figure; image(fbimg);
%}

%%

mfac = 400; sfac = 3;

clear Zm;
Zm(:,1) = Z(:,1)-min(Z(:,1))+0.15;
Zm(:,2) = Z(:,2)-min(Z(:,2))+0.15;

zszx = (max(Zm(:,1))+0.15)*mfac/sfac;
zszy = (max(Zm(:,2))+0.15)*mfac/sfac;

numpos = nnz(labels==1);
numneg = nnz(labels==2);
ratio = numneg/numpos;
rad = sqrt(ratio);
totArea = 4*rad.^2*numpos + numneg;
bfac = sqrt(8000/totArea);

bimgy = round(zszy)+round(bfac+rad)+sfac;
bimgx = round(zszx)+round(bfac+rad)+sfac;
bimg = ones(bimgy,bimgx,3);

for ndx= 1:numel(labels)
  curx = round( Zm(ndx,1) *mfac/sfac)+1;
  cury = round( Zm(ndx,2) *mfac/sfac)+1;

  bimg(cury,curx,:) = reshape(C(ndx,:),[1 1 3]);
  if labels(ndx) ==1
    bboxy = cury+round(-rad*bfac:rad*bfac);
    bboxy(bboxy<1) = [];
    bboxx = curx+round(-rad*bfac:rad*bfac);
    bboxx(bboxx<1) = [];
    bimg(bboxy,bboxx,1) = C(ndx,1);
    bimg(bboxy,bboxx,2) = C(ndx,2);
    bimg(bboxy,bboxx,3) = C(ndx,3);
  else
    bboxy = cury+round(-bfac:bfac);
    bboxy(bboxy<1) = [];
    bboxx = curx+round(-bfac:bfac);
    bboxx(bboxx<1) = [];
    bimg(bboxy,bboxx,1) = C(ndx,1);
    bimg(bboxy,bboxx,2) = C(ndx,2);
    bimg(bboxy,bboxx,3) = C(ndx,3);

  end

end
aoccupied = bimg(:,:,2)~=1;
fprintf('Area occupied for %d:%d\n ',outqndx,nnz(aoccupied));
bbimg = imfilter(bimg,fspecial('gaussian',2*3*32,32),'symmetric');
fbimg = imresize(bbimg,sfac);

fbimg(fbimg(:)>1) = 1;
%%

overall = figure;
image(fbimg); hold on;

overall1 = figure;
image(fbimg); hold on;

bname = sprintf('/groups/branson/home/kabram/paperFigs/classifier_chase_images_%d/',outqndx);
if ~exist(bname,'dir'), mkdir(bname); end
imwrite(fbimg,fullfile(bname,'background.png'));

%
% scatter(Zm(numel(labels)+1:end,1)*mfac,Zm(numel(labels)+1:end,2)*mfac,...
%   S(numel(labels)+1:end),C(numel(labels)+1:end,:),'filled');
% hold on;

occupied = false(size(fbimg));
rem = 55;

curZ = Zm(numel(labels)+1:end,:);
px = []; py = [];
bx_p = []; by_p = [];
bx_n = []; by_n = [];
bx_np = []; by_np = [];
bx_nn = []; by_nn = [];
phandle = [];
count = 0;
for ndx = 1:numel(bouts.timestamp)
  if bouts.timestamp(ndx)>qq(2) || bouts.timestamp(ndx)<qq(1); continue; end;
  curidx = bouts.ndx(ndx,:);
  lndx = find(curidx(newidxall));
  curScores = myBoostClassify(newdata(lndx,:),bb);
  [~,minndx] = min(curScores*sign(bouts.label(ndx)-1.5));
  curlndx = lndx(minndx);
  expi = newinfo.expi(curlndx);
  flies = newinfo.flies(curlndx);
  t = newinfo.t(curlndx);
  mname = fullfile(expdirs{expi},moviefilename);
  [im,traj] = vid2im(mname,flies,t,trx{expi});

  ii = uint8(repmat(flipud(im),[1 1 3]));
  ii(:,[1:rem end-rem-1:end],:)=[];
  halfsz = round(size(ii,1)/2);
  bottomout = size(ii,1)-halfsz-30+1;
  ii(halfsz+30:end,:,:) = [];

  alpha = (ii(:,:,2)<170);
  newX = curZ(curlndx,1)*mfac-size(ii,1)/2;
  newY = curZ(curlndx,2)*mfac-size(ii,2)/2;

  testX = floor(newX:newX+size(ii,1)); testX(testX<1) = []; testX(testX>size(occupied,1))=[];
  testY = floor(newY:newY+size(ii,1)); testY(testY<1) = []; testY(testY>size(occupied,2))=[];
  oo = occupied(testX,testY);
  if(nnz(oo(:))/numel(oo)>0.4), continue; end;
  count = count+1;

  if bouts.label(ndx) == 1
    bx_np = [bx_np nan newX newX newX+size(ii,2) newX+size(ii,2) newX];
    by_np = [by_np nan newY newY+size(ii,1) newY+size(ii,1) newY newY];
  else
    bx_nn = [bx_nn nan newX newX newX+size(ii,2) newX+size(ii,2) newX];
    by_nn = [by_nn nan newY newY+size(ii,1) newY+size(ii,1) newY newY];
  end

  hszx = size(ii,1) + bottomout;

  occupied(testX,testY) = true;

  figure(overall1);
  if bouts.label(ndx) == 1
    plot(newX+size(ii,1)/2,newY+size(ii,2)/2,'Marker','*','MarkerEdgeColor',[1 0.6 0],'MarkerSize',3);
  else
    plot(newX+size(ii,2)/2,newY+size(ii,2)/2,'Marker','*','MarkerEdgeColor',[0 0.6 1],'MarkerSize',3);
  end

  imwrite(ii,sprintf('%s/individual_images_test_%d.png',bname,count));

  imagetrxtemp = figure; imagebothtemp = figure; hold on;
  figure(imagebothtemp);
  image(1,1,ii,'AlphaData',alpha);
  plot(-traj.Y+hszx-rem,-traj.X+hszx,'LineStyle','none','Marker','.','MarkerSize',10,'MarkerEdgeColor','k');
  axis off; axis image;
  fname = sprintf('%s/individual_both_test_%d',bname,count);
  SaveFigLotsOfWays(imagebothtemp,fname);

  figure(imagetrxtemp); axis off; axis image;
  plot(-traj.Y+hszx-rem,-traj.X+hszx,'LineStyle','none','Marker','.','MarkerSize',10,'MarkerEdgeColor','k');
  fname = sprintf('%s/individual_trx_test_%d',bname,count);
  axis off; axis image;
  SaveFigLotsOfWays(imagetrxtemp,fname);

  close(imagebothtemp); close(imagetrxtemp);

  figure(overall);
  image(newX,newY,ii,'AlphaData',alpha);
  hszx = size(ii,1) + bottomout;
  px = [px nan -traj.Y+newX+hszx-rem];
  py = [py nan -traj.X+newY+hszx];

end

pxtest = px;
pytest = py;
px = []; py = [];

curZ = Zm(1:numel(labels),:);
count = 0;
for ndx = 1:numel(bouts.timestamp)
  if bouts.timestamp(ndx)>qq(1), continue; end;
  curidx = bouts.ndx(ndx,:);
  lndx = find(curidx(oldidxall));
  curlndx = round(median(lndx));
  expi = info.expi(curlndx);
  flies = info.flies(curlndx);
  t = info.t(curlndx);
  mname = fullfile(expdirs{expi},moviefilename);
  [im,traj] = vid2im(mname,flies,t,trx{expi});

  ii = uint8(repmat(flipud(im),[1 1 3]));
  ii(:,[1:rem end-rem-1:end],:)=[];
  halfsz = round(size(ii,1)/2);
  bottomout = size(ii,1)-halfsz-30+1;
  ii(halfsz+30:end,:,:) = [];
  alpha = (ii(:,:,2)<170);
  newX = curZ(curlndx,1)*mfac - size(ii,1)/2;
  newY = curZ(curlndx,2)*mfac - size(ii,2)/2;

  testX = floor(newX:newX+size(ii,1)); testX(testX<1) = []; testX(testX>size(occupied,1))=[];
  testY = floor(newY:newY+size(ii,1)); testY(testY<1) = []; testY(testY>size(occupied,2))=[];
  oo = occupied(testX,testY);
  if(nnz(oo(:))/numel(oo)>0.4), continue; end;
  count = count+1;

  if bouts.label(ndx) == 1
    bx_p = [bx_p nan newX newX newX+size(ii,2) newX+size(ii,2) newX];
    by_p = [by_p nan newY newY+size(ii,1) newY+size(ii,1) newY newY];
  else
    bx_n = [bx_n nan newX newX newX+size(ii,2) newX+size(ii,2) newX];
    by_n = [by_n nan newY newY+size(ii,1) newY+size(ii,1) newY newY];
  end

  occupied(testX,testY) = true;

  imwrite(ii,sprintf('%s/individual_images_train_%d.png',bname,count));
  imagetrxtemp = figure; imagebothtemp = figure; hold on;
  figure(imagebothtemp);axis off; axis image;
  image(1,1,ii,'AlphaData',alpha);
  plot(-traj.Y+hszx-rem,-traj.X+hszx,'LineStyle','none','Marker','.','MarkerSize',6,'MarkerEdgeColor','k');
  fname = sprintf('%s/individual_both_train_%d',bname,count);
  axis off; axis image;
  SaveFigLotsOfWays(imagebothtemp,fname);

  figure(imagetrxtemp);axis off; axis image;
  plot(-traj.Y+hszx-rem,-traj.X+hszx,'LineStyle','none','Marker','.','MarkerSize',6,'MarkerEdgeColor','k');
  fname = sprintf('%s/individual_trx_train_%d',bname,count);
  axis off; axis image;
  SaveFigLotsOfWays(imagetrxtemp,fname);

  close(imagebothtemp); close(imagetrxtemp);

  figure(overall);
  image(newX,newY,ii,'AlphaData',alpha);
  figure(overall1);
  image(newX,newY,ii,'AlphaData',alpha);

  hszx = size(ii,1) + bottomout;
  px = [px nan -traj.Y+newX+hszx-rem];
  py = [py nan -traj.X+newY+hszx];
end

figure(overall);
phandle = [phandle plot(bx_np,by_np,'-','Color',[1 0.6 0])];
phandle = [phandle plot(bx_nn,by_nn,'-','Color',[0 0.6 1])];
phandle = [phandle plot(bx_p,by_p,'-','Color',[1 0 0])];
phandle = [phandle plot(bx_n,by_n,'-','Color',[0 0 1])];

figure(overall1);
phandle = [phandle plot(bx_p,by_p,'-','Color',[1 0 0])];
phandle = [phandle plot(bx_n,by_n,'-','Color',[0 0 1])];

for pndx = 1:numel(phandle);
  uistack(phandle(pndx),'bottom');
  uistack(phandle(pndx),'up');
end

figure(overall);
ppx = [pxtest px]; ppy = [pytest py];
plot(ppx,ppy,'LineStyle','none','Marker','.','MarkerSize',6,'MarkerEdgeColor','k');
% plot(ppx+0.5,ppy,'LineStyle','none','Marker','.','MarkerSize',1.5,'MarkerEdgeColor','k');
% plot(px,py+0.5,'LineStyle','none','Marker','.','MarkerSize',1.5,'MarkerEdgeColor','k');
% plot(px+0.5,py+0.5,'LineStyle','none','Marker','.','MarkerSize',1.5,'MarkerEdgeColor','k');
axis equal; axis tight; axis off;
saveas(overall,[bname '/overall.fig']);
SaveFigLotsOfWays(overall,[bname '/overall']);
close(overall);

figure(overall1);
plot(px,py,'LineStyle','none','Marker','.','MarkerSize',6,'MarkerEdgeColor','k');
% plot(px+0.5,py,'LineStyle','none','Marker','.','MarkerSize',1.5,'MarkerEdgeColor','k');
% plot(px,py+0.5,'LineStyle','none','Marker','.','MarkerSize',1.5,'MarkerEdgeColor','k');
% plot(px+0.5,py+0.5,'LineStyle','none','Marker','.','MarkerSize',1.5,'MarkerEdgeColor','k');
axis equal; axis tight; axis off;
saveas(overall1,[bname '/overall1.fig']);
SaveFigLotsOfWays(overall1,[bname '/overall1']);
close(overall1);

% title(sprintf('Train:%d',outqndx));
% fname = sprintf('/groups/branson/home/kabram/paperFigs/classifier_chase_new_plotted_%d',outqndx);
% SaveFigLotsOfWays(hfig,fname);
% close(hfig);
end

author: kristinbranson

repo: JAABA


text:
function c = acorr(x,maxlag)
%ACORR Estimate autocorrelation function of time series
%
%   C = ACORR(X,MAXLAG) returns normalized autocorrelation
%   sequences for each column of X using
%   C(:,i)=XCORR(X(:,i)-MEAN(X(:,i)),MAXLAG,'coeff'), but returns only
%   lags 1:MAXLAG. Default MAXLAG = M-1;
%
%   See also
%     XCORR
%
% Copyright (C) 2000 Aki Vehtari

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

if nargin < 1
  error('Not enough input arguments.');
end
if nargin < 2
  maxlag=length(x)-1;
end
[m,n]=size(x);
c=zeros(maxlag,n);
for i1=1:n
  ct=xcorr(x(:,i1)-mean(x(:,i1)),maxlag,'coeff');
  ct=ct(maxlag+2:end);
  c(:,i1)=ct;
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
classdef NNOP < Algorithm
    %NNOP Neural Network with Ordered Partitions (NNOP). This model
    % considers the OrderedPartitions coding scheme for the labels and a
    % rule for decisions based on the first node whose output is higher
    % than a predefined threshold (T=0.5, in our experiments). The
    % model has one hidden layer with hiddenN neurons and one outputlayer
    % with as many neurons as the number of classes minus one. The learning
    % is based on iRProp+ algorithm and the implementation provided by
    % Roberto Calandra in his toolbox Rprop Toolbox for {MATLAB}:
    % http://www.ias.informatik.tu-darmstadt.de/Research/RpropToolbox
    % The model is adjusted by minimizing mean squared error. A regularization
    % parameter "lambda" is included based on L2, and the number of
    % iterations is specified by the "iter" parameter.
    %   NNPOM methods:
    %      fitpredict               - runs the corresponding algorithm,
    %                                   fitting the model and testing it
    %                                   in a dataset.
    %      fit                        - Fits a model from training data
    %      predict                    - Performs label prediction
    %
    %   NNPOM properties:
    %      epsilonInit                - Range for initializing the weights.
    %      parameters.hiddenN         - Number of hidden neurons of the
    %                                   model.
    %      parameters.iter            - Number of iterations for iRProp+
    %                                   algorithm.
    %      parameters.lambda          - Regularization parameter.
    %
    %   References:
    %     [1] J. Cheng, Z. Wang, and G. Pollastri, "A neural network
    %         approach to ordinal regression," in Proc. IEEE Int. Joint
    %         Conf. Neural Netw. (IEEE World Congr. Comput. Intell.), 2008,
    %         pp. 1279-1284.
    %     [2] P.A. Gutiérrez, M. Pérez-Ortiz, J. Sánchez-Monedero,
    %         F. Fernández-Navarro and C. Hervás-Martínez
    %         Ordinal regression methods: survey and experimental study
    %         IEEE Transactions on Knowledge and Data Engineering, Vol. 28.
    %         Issue 1, 2016
    %         http://dx.doi.org/10.1109/TKDE.2015.2457911
    %
    %   This file is part of ORCA: https://github.com/ayrna/orca
    %   Original authors: Pedro Antonio Gutiérrez, María Pérez Ortiz, Javier Sánchez Monedero
    %   Citation: If you use this code, please cite the associated paper http://www.uco.es/grupos/ayrna/orreview
    %   Copyright:
    %       This software is released under the The GNU General Public License v3.0 licence
    %       available at http://www.gnu.org/licenses/gpl-3.0.html

    properties
        description = 'Neural Network with Ordered Partitions';
        % Weights range
        epsilonInit = 0.5;
        parameters = struct('iter', 500,'hiddenN', 50,'lambda', 0.01);
    end

    methods

        function obj = NNOP(varargin)
            %NNOP constructs an object of the class NNOP and sets its default
            %   characteristics
            %   obj = NNOP('epsilonInit', 0.5) sets initialization of
            %   epsilon to 0.5
            obj.parseArgs(varargin);
        end

        function obj = set.epsilonInit(obj,e)
            if strcmp(class(obj.epsilonInit), class(e))
                obj.epsilonInit= e;
            else
                error('epsilonInit type is ''%s'' and ''%s'' was provided', class(obj.epsilonInit), class(e))
            end
        end

        function [projectedTrain, predictedTrain] = privfit( obj, train, parameters)
            %PRIVFIT trains the model for the NNOP method with TRAIN data and
            %vector of parameters PARAM.

            % Aux variables
            X = train.patterns;
            y = train.targets;
            input_layer_size  = size(X,2);
            hidden_layer_size = parameters.hiddenN;
            num_labels = numel(unique(y));
            m = size(X,1);

            % Recode y to Y using ordered partitions
            Y = repmat(y,1,num_labels) <= repmat((1:num_labels),m,1);

            % Hidden layer weigths (with bias)
            initial_Theta1 = obj.randInitializeWeights(input_layer_size+1, hidden_layer_size);
            % Output layer weigths (without bias, the biases will be the
            %                       Thresholds)
            initial_Theta2 = obj.randInitializeWeights(hidden_layer_size+1, num_labels-1);

            % Pack parameters
            initial_nn_params = [initial_Theta1(:) ; initial_Theta2(:)];

            % Set regularization parameter
            lambda = parameters.lambda;

            % Create "short hand" for the cost function to be minimized
            costFunction = @(p) obj.nnOPCostFunction(p, ...
                input_layer_size, ...
                hidden_layer_size, ...
                num_labels, X, Y, lambda);

            % RProp options
            p.verbosity = 0;                    % Increase indent
            p.MaxIter   = parameters.iter;     	% Maximum number of iterations
            p.d_Obj     = -1;                   % Objective cost
            p.method    = 'IRprop+';            % Use IRprop- algorithm
            p.display   = 0;

            % Running RProp
            [nn_params,cost,exitflag,stats1] = rprop(costFunction,initial_nn_params,p);

            %             options = optimoptions('fminunc','Algorithm','quasi-newton','SpecifyObjectiveGradient',true,'Diagnostics','on','Display','iter-detailed','UseParallel',true,'MaxIter', 1000,'CheckGradients',true);
            %             [nn_params, cost, exitflag, output] = fminunc(costFunction, initial_nn_params, options);

            % Unpack the parameters
            [Theta1, Theta2] = obj.unpackParameters(nn_params,input_layer_size,hidden_layer_size,num_labels);
            model.Theta1=Theta1;
            model.Theta2=Theta2;
            model.num_labels=num_labels;
            model.m = m;
            model.parameters = parameters;
            obj.model = model;
            [projectedTrain, predictedTrain] = obj.predict(train.patterns);
        end

        function [projected, predicted]= predict(obj,test)
            %PREDICT predicts labels of TEST patterns labels. The object needs to be fitted to the data first.
            m = size(test,1);
            a1 = [ones(m, 1) test];
            z2 = [ones(m, 1) a1*obj.model.Theta1'];
            a2 =  1.0 ./ (1.0 + exp(-z2));
            projected=a2*obj.model.Theta2';
            projected=1.0 ./ (1.0 + exp(-projected));

            a3 = ([projected ones(m,1)] > 0.5).*repmat(1:obj.model.num_labels,m,1);
            a3(a3==0)=obj.model.num_labels+1;

            predicted = min(a3,[],2);
        end

    end

    methods(Access = private)

        function [Theta1, Theta2] = unpackParameters(obj,nn_params,input_layer_size,hidden_layer_size,num_labels)
            % UNPACKPARAMETERS obtains Theta1 and Theta2
            % back from the whole array nn_params
            nTheta1 = hidden_layer_size * (input_layer_size + 1);
            Theta1 = reshape(nn_params(1:nTheta1), ...
                hidden_layer_size, (input_layer_size + 1));
            Theta2 = reshape(nn_params((1+nTheta1):end), ...
                num_labels-1, (hidden_layer_size+1));
        end

        function W = randInitializeWeights(obj, L_in, L_out)
            %RANDINITIALIZEWEIGHTS randomly initializes the weights of a layer with L_in
            %incoming connections and L_out outgoing connections
            W = rand(L_out, L_in)*2*obj.epsilonInit - obj.epsilonInit;
        end

        function [J,grad] = nnOPCostFunction(obj, nn_params, ...
                input_layer_size, ...
                hidden_layer_size, ...
                num_labels, ...
                X, Y, lambda)
            %NNPOMCOSTFUNCTION implements the cost function and obtains the
            %corresponding derivatives.

            % Unroll all the parameters
            [Theta1, Theta2] = unpackParameters(obj,...
                nn_params,input_layer_size,hidden_layer_size,num_labels);

            % Setup some useful variables
            m = size(X, 1);

            % Neural Network model
            a1 = [ones(m, 1) X];
            z2 = a1*Theta1';
            a2 =  [ones(m, 1) (1.0 ./ (1.0 + exp(-z2)))];
            z3=a2*Theta2';
            h =  [1.0 ./ (1.0 + exp(-z3)) ones(m, 1)];

            % Final output
            out = h;

            % calculte penalty (regularización L2)
            p = sum(sum(Theta1(:, 2:end).^2, 2))+sum(sum(Theta2(:, 2:end).^2, 2));

            % MSE
            J = sum(sum((out-Y).^2, 2))/(2*m) + lambda*p/(2*m);
            % Cross entropy
            %J = sum(-log(out(Y==1)), 1)/m + lambda*p/(2*m);
            if nargout > 1
                % Cross entropy
                %out(out<0.00001)=0.00001;
                %errorDer = zeros(size(Y));
                %errorDer(Y~=0) = (-Y(Y~=0)./out(Y~=0));

                % MSE
                errorDer=(out-Y);

                % Calculate sigmas
                sigma3 = errorDer.*h.*(1-h);
                sigma3 = sigma3(:,1:(end-1));
                sigma2 = (sigma3*Theta2).*a2.*(1-a2);
                sigma2 = sigma2(:, 2:end);

                % Accumulate gradients
                delta_1 = (sigma2'*a1);
                delta_2 = (sigma3'*a2);

                % calculate regularized gradient
                p1 = (lambda/m)*[zeros(size(Theta1, 1), 1) Theta1(:, 2:end)];
                p2 = (lambda/m)*[zeros(size(Theta2, 1), 1) Theta2(:, 2:end)];
                Theta1_grad = delta_1./m + p1;
                Theta2_grad = delta_2./m + p2;

                % Unroll gradients
                grad = [Theta1_grad(:) ; Theta2_grad(:)];
            end

        end
    end
end

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
function model = weakTrain(X, Y, opts)
% weak random learner
% can currently train:
% 1. decision stump: look along random dimension of data, choose threshold
% that maximizes information gain in class labels
% 2. 2D linear decision learner: same as decision stump but in 2D. I know,
% in general this could all be folded into single linear stump, but I make
% distinction for historical, cultural, and efficiency reasons.
% 3. Conic section learner: second order learning in 2D. i.e. x*y is a
% feature in addition to x, y and offset (as in 2.)
% 4. Distance learner. Picks a data point in train set and a threshold. The
% label is computed based on distance to the data point

classifierID= 1; % by default use decision stumps only
numSplits= 30;
classifierCommitFirst= true;

if nargin < 3, opts = struct; end
if isfield(opts, 'classifierID'), classifierID = opts.classifierID; end
if isfield(opts, 'numSplits'), numSplits = opts.numSplits; end
if isfield(opts, 'classifierCommitFirst'), classifierCommitFirst = opts.classifierCommitFirst; end

if classifierCommitFirst
    % commit to a weak learner first, then optimize its parameters only. In
    % this variation, different weak learners don't compete for a node.
    if length(classifierID)>1
        classifierID= classifierID(randi(length(classifierID)));
    end
end

u= unique(Y);
[N, D]= size(X);

if N == 0
    % edge case. No data reached this leaf. Don't do anything...
    model.classifierID= 0;
    return;
end

bestgain= -100;
model = struct;
% Go over all applicable classifiers and generate candidate weak models
for classf = classifierID

    modelCandidate= struct;
    maxgain= -1;

    if classf == 1
        % Decision stump

        % proceed to pick optimal splitting value t, based on Information Gain
        for q= 1:numSplits

            if mod(q-1,5)==0
                r= randi(D);
                col= X(:, r);
                tmin= min(col);
                tmax= max(col);
            end

            t= rand(1)*(tmax-tmin)+tmin;
            dec = col < t;
            Igain = evalDecision(Y, dec, u);

            if Igain>maxgain
                maxgain = Igain;
                modelCandidate.r= r;
                modelCandidate.t= t;
            end
        end

    elseif classf == 2
        % Linear classifier using 2 dimensions

        % Repeat some number of times:
        % pick two dimensions, pick 3 random parameters, and see what happens
        for q= 1:numSplits

            r1= randi(D);
            r2= randi(D);
            w= randn(3, 1);

            dec = [X(:, [r1 r2]), ones(N, 1)]*w < 0;
            Igain = evalDecision(Y, dec, u);

            if Igain>maxgain
                maxgain = Igain;
                modelCandidate.r1= r1;
                modelCandidate.r2= r2;
                modelCandidate.w= w;
            end
        end

    elseif classf == 3
        % Conic section weak learner in 2D (not too good presently, what is the
        % best way to randomly suggest good parameters?

        % Pick random parameters and see what happens
        for q= 1:numSplits

            if mod(q-1,5)==0
                r1= randi(D);
                r2= randi(D);
                w= randn(6, 1);
                phi= [X(:, r1).*X(:, r2), X(:,r1).^2, X(:,r2).^2, X(:, r1), X(:, r2), ones(N, 1)];
                mv= phi*w;
            end

            t1= randn(1);
            t2= randn(1);
            if rand(1)<0.5, t1=-inf; end
            dec= mv<t2 & mv>t1;
            Igain = evalDecision(Y, dec, u);

            if Igain>maxgain
                maxgain = Igain;
                modelCandidate.r1= r1;
                modelCandidate.r2= r2;
                modelCandidate.w= w;
                modelCandidate.t1= t1;
                modelCandidate.t2= t2;
            end
        end

    elseif classf==4
        % RBF weak learner: Picks an example and bases decision on distance
        % threshold

        % Pick random parameters and see what happens
        for q= 1:numSplits

            % this is expensive, lets only recompute every once in a while...
            if mod(q-1,5)==0
                x= X(randi(size(X, 1)), :);
                dsts= pdist2(X, x);
                maxdsts= max(dsts);
                mindsts= min(dsts);
            end

            t= rand(1)*(maxdsts - mindsts)+ mindsts;
            dec= dsts < t;
            Igain = evalDecision(Y, dec, u);

            if Igain>maxgain
                maxgain = Igain;
                modelCandidate.x= x;
                modelCandidate.t= t;
            end
        end

    else
        fprintf('Error in weak train! Classifier with ID = %d does not exist.\n', classf);
    end

    % see if this particular classifier has the best information gain so
    % far, and if so, save it as the best choice for this node
    if maxgain >= bestgain
        bestgain = maxgain;
        model= modelCandidate;
        model.classifierID= classf;
    end

end

end

function Igain= evalDecision(Y, dec, u)
% gives Information Gain provided a boolean decision array for what goes
% left or right. u is unique vector of class labels at this node

    YL= Y(dec);
    YR= Y(~dec);
    H= classEntropy(Y, u);
    HL= classEntropy(YL, u);
    HR= classEntropy(YR, u);
    Igain= H - length(YL)/length(Y)*HL - length(YR)/length(Y)*HR;

end

% Helper function for class entropy used with Decision Stump
function H= classEntropy(y, u)

    cdist= histc(y, u) + 1;
    cdist= cdist/sum(cdist);
    cdist= cdist .* log(cdist);
    H= -sum(cdist);

end

author: karpathy

repo: Random-Forest-Matlab


text:
function [p] = dot(tt1,tt2)
%Frobenius dot product of two TT-matrices
%   [PR]=DOT(TT1,TT2) Frobenius dot product of two TT-matrices
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
p=dot(tt1.tt,tt2.tt);
return
end

author: oseledets

repo: TT-Toolbox


text:
function [stat,pval,statAll,pvalAll]=grangercause(y,constant,lags,het,uncorr,inference)
% Granger causality testing with a variance-covariance matrix estimated under a variety of
% assumptions on the covariance of the errors:
%   * Conditionally Homoskedastic and Uncorrelated
%   * Conditionally Homoskedastic but Correlated
%   * Heteroskedastic but Conditionally Uncorrelated
%   * Heteroskedastic and Correlated
%
% USAGE:
%   [STAT] = grangercause(Y,CONSTANT,LAGS)
%   [STAT,PVAL,STATALL,PVALALL] = grangercause(Y,CONSTANT,LAGS,HET,UNCORR,INFERENCE)
%
% INPUTS:
%   Y             - A T by K matrix of data
%   CONSTANT      - Scalar variable: 1 to include a constant, 0 to exclude
%   LAGS          - Non-negative integer vector representing the VAR orders to include in the model.
%   HET           - [OPTIONAL] A scalar integer indicating the type of covariance estimator
%                      0 - Homoskedastic
%                      1 - Heteroskedastic [DEFAULT]
%   UNCORR        - [OPTIONAL] A scalar integer indicating the assumed structure of the error
%                     covariance matrix
%                      0 - Correlated errors  [DEFAULT]
%                      1 - Uncorrelated errors
%   INFERENCE     - [OPTIONAL] Inference method
%                      1 - Likelihood ratio [DEFAULT]
%                      2 - LM test
%                      3 - Wald test
%
% OUTPUTS:
%   STAT          - K by K matrix of Granger causality statistics computed using the specified
%                     covariance estimator and inference method. STAT(i,j) corresponds to a test that
%                     y(i) is not caused by y(j)
%   PVAL          - K by K matrix of p-values corresponding to STAT
%   STATALL       - K by 1 vector of Granger causality statistics computed using the specified
%                     covariance estimator and inference method. STATALL(i) corresponds to a test that
%                     y(i) is not caused by any y(j), j neq i
%   PVALALL       - K by 1 vector of p-values corresponding to STATALL
%
% COMMENTS:
%   Granger causality tests based on a VAR including any lags.
%
%   y(:,t)' = CONST + P(1) * y(:,y-1) + P(2)*y(:,y-2) + ... + P(1)*y(:,t-K)'
%
%   where P(j) are K by K parameter matrices and CONST is a K by 1 parameter matrix (if CONSTANT==1)
%
% EXAMPLE:
%   Conduct GC testing in a VAR(1) with a constant
%        parameters = grangercause(y,1,1)
%   Conduct GC testing in a VAR(3) with no constant
%        parameters = grangercause(y,0,[1:3])
%   Conduct GC testing in a VAR that includes lags 1 and 3 with a constant
%        parameters = grangercause(y,1,[1 3])
%
% See also VECTORAR, VECTORARVCV

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3.0    Date: 1/1/2007

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<3 || nargin>6
    error('3 to 6 inputs required')
end
if nargin==3
    het=1;
    uncorr=0;
    inference=1;
elseif nargin==4
    uncorr=0;
    inference=1;
elseif nargin==5
    inference=1;
end
% Check Y
if ndims(y)~=2
    error('Y must be T by K')
end
% Check constant
if ~isscalar(constant)
    error('CONSTANT must be either 0 or 1')
end
if ~ismember(constant,[0 1])
    error('CONSTANT must be either 0 or 1')
end
% Check lags
if ndims(lags)~=2
    error('LAGS must be a vector of positive integers containing lags to include')
end
if size(lags,1)>size(lags,2)
    lags = lags';
end
if ~all(lags>0)
    error('LAGS must be a vector of positive integers containing lags to include')
end
if ~all(floor(lags)==lags)
    error('LAGS must be a vector of positive integers containing lags to include')
end
if length(lags)~=length(unique(lags))
    error('LAGS must be a vector of unique elements')
end
lags=sort(lags);
% Check het
if ~isscalar(het)
    error('HET must be either 0 or 1')
end
if ~ismember(het,[0 1])
    error('HET must be either 0 or 1')
end
% Check uncorr
if ~isscalar(uncorr)
    error('UNCORR must be either 0 or 1')
end
if ~ismember(uncorr,[0 1])
    error('UNCORR must be either 0 or 1')
end
% Inference mathod
if ~isscalar(inference)
    error('INFERENCE must be either 0 or 1')
end
if ~ismember(inference,[1 2 3])
    error('INFERENCE must be either 0 or 1')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Inference Methods
% 1. LR MLE : Relatively easy, need to iteratively depete the correct
% columns and reestimate.  Use a LR test
% 2. Wald: Just compute the VCV, then figure out how to select the correct
%       rows, and you are done
% 3. LM: Have to re-estimate under the null and use those errors to compute
%       the scores.  Test is based on scores.
% 4. LR Robust: Just like LM above but using errors estimated under the
%       alternative

% LR MLE
% First generate the unrestricted errors, then iterate across series and
% drop regressors

T=size(y,1);
K=size(y,2);
P=length(lags);
m=max(lags);
ylags=cell(K,1);
for k=1:K
    ytemp = [ones(m,1)*mean(y(:,k));
        y(:,k)];
    [nothing,ylags{k}]=newlagmatrix(ytemp,m,0);
end

X=zeros(T,K+constant);
index=1;
if constant
    X(:,1)=ones(T,1);
    index=index+1;
end
p=length(lags);
for i=1:p
    for k=1:K
        X(:,index)=ylags{k}(:,lags(i));
        index=index+1;
    end
end

% Each column contains the parameters for a single y
paramvec = X\y;
errors = y-X*paramvec;
s2=errors'*errors/T;
e=cell(K,K);
for i=1:K
    for j=1:K
        % i will be the LHS variable, j will be the lags
        tempX=X;
        allcols = 1:size(X,2);
        % Drop columns j, K+j, 2K+j,...P(K-1)*j (+constant if present)
        if constant
            drop = (j:K:size(X,2)-1)+constant;
        else
            drop = (j:K:size(X,2));
        end
        remain = setdiff(allcols,drop);
        beta = tempX(:,remain)\y(:,i);
        % These errors are estiamted under the null and can be used for LR and LM testing
        e{i,j} = y(:,i)-tempX(:,remain)*beta;
    end
end
% Construct the regressions needed to test the "all" hypothesis
eall = cell(K,1);
for i=1:K
    tempX = X;
    % Keep columns i, K+i, 2K+i, P(K-1)*i (+constant if present)
    if constant
        remain = [constant (i:K:(size(X,2)-1))+constant];
    else
        remain = i:K:size(X,2);
    end
    beta = tempX(:,remain)\y(:,i);
    eall{i} = y(:,i)-tempX(:,remain)*beta;
end

Np = p*K+constant;
stat = zeros(K);
statAll = zeros(K,1);
% Homoskedastic Likelihood Ratio
if inference==1 && het==0
    for i=1:K
        for j=1:K
            e2=errors;
            e2(:,i)=e{i,j};
            sR=e2'*e2/T;
            stat(i,j)=(T-P*K^2+P)*(log(det(sR))-log(det(s2)));
        end
    end
    for i=1:K
        e2=errors;
        e2(:,i) = eall{i};
        sR=e2'*e2/T;
        statAll(i)=(T-P*K^2+(K-1)*P)*(log(det(sR))-log(det(s2)));
    end
    % Heteroskedasticity robust LR;  key here is score covariance is computed
    % using errors estimated under the alternative
elseif inference==1 && het==1
    for i=1:K
        for j=1:K
            e2=errors;
            e2(:,i)=e{i,j};
            X2=repmat(X,1,K);
            e2=reshape(repmat(e2,Np,1),T,K*Np);
            s=X2.*e2;
            sbar=mean(s);
            % Estimate S using the errors, not e2
            S=vectorarscorecov(errors,X,het,uncorr,K,T,Np);
            stat(i,j)=(T-P*K^2+P)*sbar*S^(-1)*sbar';
        end
    end
    for i=1:K
        e2=errors;
        e2(:,i)=eall{i};
        X2=repmat(X,1,K);
        e2=reshape(repmat(e2,Np,1),T,K*Np);
        s=X2.*e2;
        sbar=mean(s);
        % Estimate S using the errors, not e2
        S=vectorarscorecov(errors,X,het,uncorr,K,T,Np);
        statAll(i)=(T-P*K^2+(K-1)*P)*sbar*S^(-1)*sbar';
    end
    % All LM tests since function will handle;  the key here is that the errors
    % are computed under the null
elseif inference==2
    for i=1:K
        for j=1:K
            e2=errors;
            e2(:,i)=e{i,j};
            X2=repmat(X,1,K);
            e2=reshape(repmat(e2,Np,1),T,K*Np);
            s=X2.*e2;
            sbar=mean(s);
            e2=errors;
            e2(:,i)=e{i,j};
            S=vectorarscorecov(e2,X,het,uncorr,K,T,Np);
            stat(i,j)=(T-P*K^2+P)*sbar*S^(-1)*sbar';
        end
    end
    for i=1:K
        e2=errors;
        e2(:,i)=eall{i};
        X2=repmat(X,1,K);
        e2=reshape(repmat(e2,Np,1),T,K*Np);
        s=X2.*e2;
        sbar=mean(s);
        e2=errors;
        e2(:,i)=eall{i};
        S=vectorarscorecov(e2,X,het,uncorr,K,T,Np);
        statAll(i)=(T-P*K^2+(K-1)*P)*sbar*S^(-1)*sbar';
    end
    % Wald tests
elseif inference==3
    XpXi = ((X'*X)/T)^(-1);
    Ainv = kron(eye(K),XpXi);
    B = vectorarscorecov(errors,X,het,uncorr,K,T,Np);
    V = Ainv*B*Ainv;
    %Vinv=V^(-1);
    % Now I have to cleverly select the parameters to set to 0 and then
    % compute the wald tests
    % Easier to iterate over the j's in the inside loop
    for i=1:K
        for j=1:K
            temp=zeros(size(paramvec))';
            % Select the columns of the parameters to test
            if constant
                pl = (j:K:size(X,2)-1)+constant;
            else
                pl = (j:K:size(X,2));
            end

            temp(i,pl)=1;
            temp=temp';
            temp=temp(:);
            pl = find(temp);
            p = paramvec;
            p = p(:);
            stat(i,j)=(T-P*K^2+P)*(p(pl)'*V(pl,pl)^(-1)*p(pl));
        end
    end
    for i=1:K
        temp=zeros(size(paramvec))';
        % Select the columns of the parameters to test
        if constant
            remain = [constant (i:K:size(X,2)-1)+constant];
        else
            remain = (i:K:size(X,2)-1);
        end
        pl = setdiff(1:size(X,2),remain);
        temp(i,pl)=1;
        temp=temp';
        temp=temp(:);
        pl = find(temp);
        p = paramvec;
        p = p(:);
        statAll(i)=(T-P*K^2+(K-1)*P)*(p(pl)'*V(pl,pl)^(-1)*p(pl));
    end
end
% All stats have the same dist
pval=1-chi2cdf(stat,length(lags));
pvalAll = 1-chi2cdf(statAll,(K-1)*length(lags));

function S=vectorarscorecov(errors,X,het,uncorr,K,T,Np)
s2=errors'*errors/T;
XpX=X'*X/T;
if ~het && uncorr
    S=kron(diag(diag(s2)),XpX);
elseif ~het && ~uncorr
    S=kron(s2,XpX);
elseif het && uncorr
    X2=repmat(X,1,K);
    e2=reshape(repmat(errors,Np,1),T,K*Np);
    s=X2.*e2;
    s=s-repmat(mean(s),T,1);
    S=zeros(Np*K);
    for i=1:K
        sel=(i-1)*Np+1:i*Np;
        temp = s(:,sel);
        S(sel,sel)=temp'*temp/T;
    end
elseif het && ~uncorr
    X2=repmat(X,1,K);
    e2=reshape(repmat(errors,Np,1),T,K*Np);
    s=X2.*e2;
    s=s-repmat(mean(s),T,1);
    S=s'*s/T;
end

author: bashtage

repo: mfe-toolbox


text:
function [fileInfo] = getFileInfo(fileStr, varargin)
	% Gets file information for subject based on the file path, returns a structure with various information.
	% Biafra Ahanonu
	% started: 2013.11.04 [12:38:42]
	% inputs
		% fileStr - character string
	% options
		% assayList
	% outputs
		%

	[fileInfo] = ciapkg.io.getFileInfo(fileStr,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [M]=gen_MassMatrix(q,Pcii,Icii,mcii)
%% This function is used to calculate the mass matrix of the KUKA iiwa 7 R 800 manipulator

% Arreguments:
%--------------------
% q: is 1x7 vector, joint nagles vector of the manipulator
% Pcii: is 3X7 matrix while each column represents the local coordinates
% of the center of mass of each link.
% Icii: is (3x3x7) matrix, each 3x3 matrix of which represnets the
% associated link inertial tensor represented in its local inertial frame
% mcii: is (1x7) vector, each element of which specifies a mass of one of
% the links

% Return value:
%--------------------
% M: 7x7 matrix, the mass matrix of the KUKA iiwa 7 R 800.

% Copyright: Mohammad SAFEEA

T=kukaGetKenimaticModelAccelerated(q);
[M]=GetInertiaMatrix5(T,Pcii,Icii,mcii);

repo: KST-Kuka-Sunrise-Toolbox


text:
function score = rgb_fitness ( m, n, k, a, b )

%*****************************************************************************80
%
%% RGB_FITNESS scores image B as a match for image A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Nick Berry,
%    A "Practical" Use for Genetic Programming,
%    http://www.datagenetics.com/blog.html
%
%  Parameters:
%
%    Input, integer M, N, K, the dimensions of the images.
%
%    Input, uint8 A(M,N,K), the original image.
%
%    Input, uint8 B(M,N,K), the approximate image.
%
%    Output, int SCORE, a score for the matching.
%
  score = sum ( sum ( sum ( abs ( double ( a ) - double ( b ) ) ) ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% import_STL_txt
% Below is a demonstration of the features of the |import_STL_txt| function

%% Syntax
% |[stlStruct] = import_STL_txt(fileName);|

%% Description
% Use |import_STL_txt| to import .txt type (as apposed to binary) STL
% files. The function supports multi-solid STL.

%% Examples

clear; close all; clc;

%%
% Plot settings
faceAlpha=0.5;
fontSize=25;

%% Import STL file as patch data

%Set main folder
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
pathName=fullfile(defaultFolder,'data','STL');
fileName=fullfile(pathName,'femur.stl');
[stlStruct] = import_STL_txt(fileName);

F=stlStruct.solidFaces{1};
V=stlStruct.solidVertices{1};

%%
% Merging nodes example
[F,V]=mergeVertices(F,V);

%%
% Plotting the model

cFigure;
title('Imported patch data from STL','fontSize',fontSize);
gpatch(F,V,'g');
axisGeom;
camlight('headlight');
lighting phong; axis off;
drawnow;

%% Importing a multi-solid STL file as patch data

%Set main folder
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
pathName=fullfile(defaultFolder,'data','STL');
fileName=fullfile(pathName,'stanford_bunny_multi.stl');
[stlStruct] = import_STL_txt(fileName);

%%
% Plotting the models
pColors=autumn(numel(stlStruct.solidNames));

cFigure;
title('Imported patch data from multi-solid STL','fontSize',fontSize);
for q=1:1:numel(stlStruct.solidNames)
    F=stlStruct.solidFaces{q};
    V=stlStruct.solidVertices{q};
    gpatch(F,V,pColors(q,:),'k',faceAlpha);
end
axisGeom;
camlight('headlight');
lighting phong; axis off;
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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


text:
function node_boundary = triangulation_order6_boundary_node ( node_num, ...
  triangle_num, triangle_node )

%*****************************************************************************80
%
%% TRIANGULATION_ORDER6_BOUNDARY_NODE indicates which nodes are on the boundary.
%
%  Discussion:
%
%    This routine is given an order 6 triangulation, an abstract list of sets
%    of six nodes.  The vertices are listed clockwise, then the
%    midside nodes.
%
%    It is assumed that each edge of the triangulation is either
%    * an INTERIOR edge, which is listed twice, once with positive
%      orientation and once with negative orientation, or;
%    * a BOUNDARY edge, which will occur only once.
%
%    This routine should work even if the region has holes - as long
%    as the boundary of the hole comprises more than 3 edges!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, integer TRIANGLE_NODE(6,TRIANGLE_NUM), the nodes that make up the
%    triangles.
%
%    Output, logical NODE_BOUNDARY(NODE_NUM), is TRUE if the node
%    is on a boundary edge.
%
  m = 3;
  n = 3 * triangle_num;
%
%  Set up the edge array.  The midside node is listed last, as
%  it is not needed for the sorting process.
%
  edge(1,               1:  triangle_num) = triangle_node(1,1:triangle_num);
  edge(2,               1:  triangle_num) = triangle_node(4,1:triangle_num);
  edge(3,               1:  triangle_num) = triangle_node(2,1:triangle_num);

  edge(1,  triangle_num+1:2*triangle_num) = triangle_node(2,1:triangle_num);
  edge(2,  triangle_num+1:2*triangle_num) = triangle_node(5,1:triangle_num);
  edge(3,  triangle_num+1:2*triangle_num) = triangle_node(3,1:triangle_num);

  edge(1,2*triangle_num+1:3*triangle_num) = triangle_node(3,1:triangle_num);
  edge(2,2*triangle_num+1:3*triangle_num) = triangle_node(6,1:triangle_num);
  edge(3,2*triangle_num+1:3*triangle_num) = triangle_node(1,1:triangle_num);
%
%  In each column, force the smaller of the two vertices to appear first.
%
  e1(1:n) = min ( edge(1:2:3,1:n) );
  e2(1:n) = max ( edge(1:2:3,1:n) );

  edge(1,1:n) = e1(1:n);
  edge(3,1:n) = e2(1:n);
%
%  Ascending sort the column array.
%
  edge = ( sortrows ( edge' ) )';
%
%  Records which appear twice are internal edges and can be ignored.
%
  node_boundary(1:node_num) = 0;

  i = 0;

  while ( i < 3 * triangle_num )

    i = i + 1;

    if ( i == 3 * triangle_num )
      node_boundary(edge(1:m,i)) = 1;
    elseif ( all ( edge(1:m,i) == edge(1:m,i+1) ) )
      i = i + 1;
    else
      node_boundary(edge(1:m,i)) = 1;
    end

  end

  return
end


text:
function [grad] = l2rowscaledg(x,y,outderv, alpha)

normeps = 1e-5;
if (~exist('outderv','var')||isempty(outderv))
    error('Requires outderv of previous layer to compute gradient!');
end

epssumsq = sum(x.^2,2) + normeps;

l2rows = sqrt(epssumsq)*alpha;

if (~exist('y','var')||isempty(y))
     y = bsxfun(@rdivide,x,l2rows);
end

grad = bsxfun(@rdivide, outderv, l2rows) - ...
       bsxfun(@times, y, sum(outderv.*x, 2) ./ epssumsq);

repo: REST


text:
function c = comment(d)

%tstoolbox/@description/comment
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,1);

c = cellstr(d.comment);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function g = sigmoid(z)
%SIGMOID Compute sigmoid functoon
%   J = SIGMOID(z) computes the sigmoid of z.

% You need to return the following variables correctly
g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the sigmoid of each value of z (z can be a matrix,
%               vector or scalar).

%g = 1 ./ (1 .+ e .^ -z);
g = 1.0 ./ (1.0 + exp(-z));%exp(n)-->e的n次方

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

end

author: scruel

repo: Notes-ML-AndrewNg


text:
function [xn,dxdf,dxdc,dxdk,dxdalpha] = normalize2(x_kk,fc,cc,kc,alpha_c),

%normalize
%
%[xn] = normalize(x_kk,fc,cc,kc,alpha_c)
%
%Computes the normalized coordinates xn given the pixel coordinates x_kk
%and the intrinsic camera parameters fc, cc and kc.
%
%INPUT: x_kk: Feature locations on the images
%       fc: Camera focal length
%       cc: Principal point coordinates
%       kc: Distortion coefficients
%       alpha_c: Skew coefficient
%
%OUTPUT: xn: Normalized feature locations on the image plane (a 2XN matrix)
%
%Important functions called within that program:

k1 = kc(1);
k2 = kc(2);
k3 = kc(5);
p1 = kc(3);
p2 = kc(4);

N = size(x_kk,2);

% First: Subtract principal point, and divide by the focal length:
x_distort = [(x_kk(1,:) - cc(1))/fc(1);(x_kk(2,:) - cc(2))/fc(2)];

v1 = - x_distort(1,:) / fc(1);
v2 = - x_distort(2,:) / fc(1);

dx_distortdfc = zeros(2*N,2);
dx_distortdfc(1:2:end,1) = v1';
dx_distortdfc(2:2:end,2) = v2';

v1 = - x_distort(1,:) / fc(1);
v2 = - x_distort(2,:) / fc(1);

dx_distortdcc = zeros(2*N,2);
dx_distortdcc(1:2:end,1) = -(1/fc(1)) * ones(N,1);
dx_distortdcc(2:2:end,2) = -(1/fc(2)) * ones(N,1);

% Second: undo skew
x_distort(1,:) = x_distort(1,:) - alpha_c * x_distort(2,:);

dx_distort2dfc = [ dx_distortdfc(:,1)-alpha_c *dx_distortdfc(:,2)   dx_distortdfc(:,2)];
dx_distort2dcc = [ dx_distortdcc(:,1)-alpha_c *dx_distortdcc(:,2)   dx_distortdcc(:,2)];

dx_distort2dalpha_c = zeros(2*N,1);
dx_distort2dalpha_c(1:2:end) = -x_distort(2,:)';

x = x_distort; 				% initial guess

for kk=1:20,

    r_2 = sum(x.^2);
    k_radial =  1 + k1 * r_2 + k2 * r_2.^2 + k3 * r_2.^3;
    delta_x = [2*p1*x(1,:).*x(2,:) + p2*(r_2 + 2*x(1,:).^2); p1 * (r_2 + 2*x(2,:).^2)+2*p2*x(1,:).*x(2,:)];
    x = (x_distort - delta_x)./(ones(2,1)*k_radial);

end;

xn = x;

dxdk = zeros(2*N,5); % Approximation (no time)
dxdf = dx_distort2dfc;
dxdc = dx_distort2dcc;
dxdalpha = dx_distort2dalpha_c;

repo: ekfmonoslam