text:
function [xi,w]=firstOrderTriangleCubPoints()
%%FIRSTORDERTRIANGLECUBPOINTS Obtain first-order cubature points for
%   integration over a triangle in 2D. The points and weights are for the
%   triangle with vertices (1,0), (0,1), (0,0), but can be transformed to
%   any triangle using transformSimplexTriPoints.
%
%INPUTS: None
%
%OUTPUTS: xi A 2XnumCubPoints set of points for the standard triangle.
%          w A 1XnumCubPoints set of cubature weights. This sums to the
%            volume of the triangle (1/2).
%
%This function implements the points given in [1] (3 points).
%
%EXAMPLE:
%Given the vertices of the simplex, we compare a first-order moment
%computed using these cubature points to one computed using
%monomialIntSimplex. The results are the same within typical finite
%precision limits.
% [xi,w]=firstOrderTriangleCubPoints();
% alpha=[1;0];
% theMoment=findMomentFromSamp(alpha,xi,w)
% intVal=monomialIntSimplex(alpha)
%
%REFERENCES:
%[1] F. D. Witherden and P. E. Vincent, "On the identification of symmetric
%    quadrature rules for finite element methods," Computer and Mathematics
%    with Applications, vol. 69, no. 10, pp. 1232-1241, May 2015.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

M=[-0.33333333333333333333333333333333333333, -0.33333333333333333333333333333333333333, 2];
w=M(:,3);
xi=M(:,1:2)';
%Transform the points to the standard triangle.
v1=[-1,-1, 1;
    -1, 1,-1];
v2=[1,0,0;
    0,1,0];
[A,d]=affineTransBetweenTriangles(v1,v2);
xi=bsxfun(@plus,A*xi,d);
w=w/4;
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 e = strainRate(L)
% strain rate

if all(L.isSymmetric(:))
  e = 0.5*max(eig(L));
else

  for n = 1:length(L)

    e(n) = 0.5 * svds(L.M(:,:,n),1);

  end

end

author: mtex-toolbox

repo: mtex

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


text:
%kasi2kdf 'imports a file from the SpectraCube system of ASI '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros asi2kdf.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% MultiChoice: t 'File Type', default: 2: 'Integer toggle selection'
%    Choices are:
%   1: 'Spectral raw-data file (.raw)'
%   2: 'Interferogram raw-data (.sdi) type 1'
%   3: 'Interferogram binary exprot (.sdm)'
%   4: 'Interferogram raw-data (.sdi) type 2'
% Integer: px 'PosX', default: 326: 'X-size position in headerfile'
% Integer: py 'PosY', default: 327: 'Y-size position in headerfile'
% Integer: pz 'PosChan', default: 325: 'Channel-size position in headerfile'
% Integer: hs 'Headersize', default: 334: 'Size of the header in Unsigned SHorts'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kasi2kdf(i, {'i','';'t',2;'px',326;'py',327;'pz',325;'hs',334;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% asi2kdf - imports a file from the SpectraCube system of ASI
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996-2003, Rainer Heintzmann,  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [X1, Y1] = Direction2Angular(D, r, c)
%
%        [X1, Y1] = Direction2Angular(D, r, c)
%
%
%        Input:
%           -D: 3D directions of the img format
%           -r: height of the angular image
%           -c: width of the angular image
%        Output:
%           -X1: X coordinates in the Angular format
%           -Y1: Y coordinates in the Angular format
%
%     Copyright (C) 2011-15  Francesco Banterle
%
%     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/>.
%

%Coordinates generation
R = acos(-D(:,:,3)) ./ (pi * 2 * sqrt(D(:,:,1).^2 + D(:,:,2).^2));

X1 = (0.5 + R .* D(:,:,1)) * c;
Y1 = (0.5 - R .* D(:,:,2)) * r;

X1 = RemoveSpecials(X1);
Y1 = RemoveSpecials(Y1);

end

author: banterle


text:
function mColormap = gui_Colormap_HSVCut(nSize)

% If size is not specified, set it to 256
if nargin < 1
  nSize = 256;
end

nTmpSize = floor(nSize/0.85);
mColormap = hsv(nTmpSize);
mColormap = mColormap(1:nSize,:);

author: CelsoReyes


text:
function dag = test_dag(A,varargin)
% TEST_DAG Tests if a graph is directed and acyclic
%
% dag = test_dag(A) returns 1 if A is the adjacency matrix for a dag and
% returns 0 otherwise.  This function just calls topological order and
% checks the output to determine if the graph is a dag.  It is provided
% merely for expository purposes.
%
% Example:
%   n = 10; A = sparse(1:n-1, 2:n, 1, n, n); % construct a simple dag
%   test_dag(A)
%   A(10,1) = 1; % complete the cycle
%   test_dag(A)

% David Gleich
% Copyright, Stanford University, 2007-2008

%% History
%%

dag = 1;
if isempty(topological_order(A,varargin{:}))
    dag = 0;
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function check = exponential_check ( a, b )

%*****************************************************************************80
%
%% EXPONENTIAL_CHECK checks the parameters of the Exponential CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameter of the PDF.
%    0.0 < B.
%
  if ( b <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'EXPONENTIAL_CHECK - Fatal error!\n' );
    fprintf ( 1, '  B <= 0.0\n' );
    check = 0;
    return
  end

  check = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function bands = extract_subbands(coeffs,pos)
pos=round(pos);
for i=1:length(pos)
    if i>1
        bands{i} = coeffs(pos(i):pos(i-1)-1);
    else
        bands{1} = coeffs(pos(1):end);
    end
end

author: alexandrebarachant


text:
function varargout = sl2dmatcov(type, data, matsiz, n, meanmat, PL, PR, w)
%SL2DMATCOV Computes the 2D matrix-covariances
%
% $ Syntax $
%   - CL = sl2dmatcov('CL', data, matsiz, n, meanmat, PL, PR, w)
%   - CR = sl2dmatcov('CR', data, matsiz, n, meanmat, PL, PR, w)
%   - [CL, CR] = sl2dmatcov('Both', data, matsiz, n, meanmat, PL, PR, w)
%
% $ Arguments $
%   - data:     the stack of matrices or the cell array of filenames of the
%               array files storing the matrices.
%   - matsiz:   the size of each matrix
%   - n:        the number of samples
%   - PL:       the left-projection matrix, default = []
%   - PR:       the right-projection matrix, default = []
%   - w:        the weights of the matrix samples, default = []
%
% $ Remarks $
%   - The function computes the 2D matrix-covariances according to
%     the following formulas:
%       Y_i = PL^T * (X - meanX)* PR
%       CL = ( sum_{i=1}^n w(i) * Y_i * Y_i^T ) / ( sum_{i=1}^n w(i) )
%       CR = ( sum_{i=1}^n w(i) * Y_i^T * Y_i ) / ( sum_{i=1}^n w(i) )
%     Following special cases are considered:
%       If w is empty, then we take w(i) = 1 for all i
%       If PL is empty, then we take PL as an identity matrix
%       If PR is empty, then we take PR as an identity matrix
%       If meanmat is 0, then we take meanmat as a zero matrix
%
% $ History $
%   - Created by Dahua Lin, on Jul 31st, 2006
%

%% parse and verify input arguments

if nargin < 5
    raise_lackinput('sl2dmatcov', 5);
end
if nargin < 6
    PL = [];
end
if nargin < 7
    PR = [];
end
if nargin < 8
    w = [];
end

matsiz = matsiz(:)';
if length(matsiz) ~= 2
    error('sltoolbox:invalidarg', ...
        'matsiz shoudl be a 2-elem vector');
end
if ~isequal(size(meanmat), matsiz)
    error('sltoolbox:invalidarg', ...
        'The size of mean matrix is not as specified');
end

d1 = matsiz(1);
d2 = matsiz(2);
if ~isempty(PL) && size(PL, 1) ~= d1
    error('sltoolbox:invalidarg', ...
        'The size of PL is inconsistent with the matrix dimension');
end
if ~isempty(PR) && size(PR, 1) ~= d2
    error('sltoolbox:invalidarg', ...
        'The size of PR is inconsistent with the matrix dimension');
end

if ~isempty(w)
    if length(w) ~= n
        error('The weights length is inconsistent with the number of samples');
    end
    tw = sum(w);
else
    tw = n;
end

%% Main body

if isnumeric(data)

    if ~isequal(size(data), [matsiz, n])
        error('sltoolbox:invalidarg', ...
            'The size of data is inconsistent with specified');
    end

    Y = compute_Y(data, meanmat, PL, PR);

    switch type
        case 'CL'
            CL = compute_SL(Y, w);
            CL = CL / tw;
            varargout = {CL};
        case 'CR'
            CR = compute_SR(Y, w);
            CR = CR / tw;
            varargout = {CR};
        case 'Both'
            [CL, CR] = compute_SLSR(Y, w);
            CL = CL / tw;
            CR = CR / tw;
            varargout = {CL, CR};
        otherwise
            error('sltoolbox:invalidarg', ...
                'invalid type: %s', type);
    end

elseif iscell(data)

    nfiles = length(data);
    cf = 0;

    switch type
        case 'CL'
            for i = 1 : nfiles
                curarr = slreadarray(data{i});
                curn = size(curarr, 3);
                Y = compute_Y(curarr, meanmat, PL, PR);
                curCL = compute_SL(Y, w);
                if i == 1
                    CL = curCL;
                else
                    CL = CL + curCL;
                end
                clear Y curCL;
                cf = cf + curn;
            end
            CL = CL / tw;
            varargout = {CL};

        case 'CR'
            for i = 1 : nfiles
                curarr = slreadarray(data{i});
                curn = size(curarr, 3);
                Y = compute_Y(curarr, meanmat, PL, PR);
                curCR = compute_SR(Y, w);
                if i == 1
                    CR = curCR;
                else
                    CR = CR + curCR;
                end
                clear Y curCR;
                cf = cf + curn;
            end
            CR = CR / tw;
            varargout = {CR};

        case 'Both'
            for i = 1 : nfiles
                curarr = slreadarray(data{i});
                curn = size(curarr, 3);
                Y = compute_Y(curarr, meanmat, PL, PR);
                [curCL, curCR] = compute_SLSR(Y, w);
                if i == 1
                    CL = curCL;
                    CR = curCR;
                else
                    CL = CL + curCL;
                    CR = CR + curCR;
                end
                clear Y curCL curCR;
                cf = cf + curn;
            end
            CL = CL / tw;
            CR = CR / tw;
            varargout = {CL, CR};

        otherwise
            error('sltoolbox:invalidarg', ...
                'invalid type: %s', type);
    end

    if cf ~= n
        error('sltoolbox:sizmismatch', ...
            'The total number of samples is not n');
    end

else
    error('sltoolbox:invalidarg', ...
        'data should be a numeric array or a cell array of filenames');
end

%% Core routine

function Y = compute_Y(X, M, PL, PR)

n = size(X, 3);
if isempty(PL)
    if isempty(PR)
        if isequal(M, 0)
            Y = X;
        else
            Y = zeros(size(X, 1), size(X, 2), n);
            for i = 1 : n
                Y(:,:,i) = X(:,:,i) - M;
            end
        end
    else
        Y = zeros(size(X, 1), size(PR, 2), n);
        if isequal(M, 0)
            for i = 1 : n
                Y(:,:,i) = X(:,:,i) * PR;
            end
        else
            for i = 1 : n
                Y(:,:,i) = (X(:,:,i) - M) * PR;
            end
        end
    end
else
    PLT = PL';
    if isempty(PR)
        Y = zeros(size(PL, 2), size(X, 2), n);
        if isequal(M, 0)
            for i = 1 : n
                Y(:,:,i) = PLT * X(:,:,i);
            end
        else
            for i = 1 : n
                Y(:,:,i) = PLT * (X(:,:,i) - M);
            end
        end
    else
        Y = zeros(size(PL, 2), size(PR, 2), n);
        if isequal(M, 0)
            for i = 1 : n
                Y(:,:,i) = PLT * X(:,:,i) * PR;
            end
        else
            for i = 1 : n
                Y(:,:,i) = PLT * (X(:,:,i) - M) * PR;
            end
        end
    end
end

function SL = compute_SL(Y, w)

SL = zeros(size(Y, 1));
n = size(Y, 3);
if isempty(w)
    for i = 1 : n
        curY = Y(:, :, i);
        SL = SL + curY * curY';
    end
else
    for i = 1 : n
        curY = Y(:, :, i);
        SL = SL + w(i) * curY * curY';
    end
end

function SR = compute_SR(Y, w)

SR = zeros(size(Y, 2));
n = size(Y, 3);
if isempty(w)
    for i = 1 : n
        curY = Y(:, :, i);
        SR = SR + curY' * curY;
    end
else
    for i = 1 : n
        curY = Y(:, :, i);
        SR = SR + w(i) * curY' * curY;
    end
end

function [SL, SR] = compute_SLSR(Y, w)

SL = zeros(size(Y, 1));
SR = zeros(size(Y, 2));
n = size(Y, 3);
if isempty(w)
    for i = 1 : n
        curY = Y(:, :, i);
        SL = SL + curY * curY';
        SR = SR + curY' * curY;
    end
else
    for i = 1 : n
        curY = Y(:, :, i);
        SL = SL + w(i) * curY * curY';
        SR = SR + w(i) * curY' * curY;
    end
end

author: lmthang


text:
%信号解扩

function out = despread(data, code)

% ****************************************************************
%   data    : 输入数据序列
%   code    : 解扩使用的扩频码序列
%   out     : 解扩后的输出数据序列
% ****************************************************************

switch nargin                           %如果输入参数个数不对，提示错误
case { 0 , 1 }
    error('缺少输入参数');
end

[hn,vn] = size(data);
[hc,vc] = size(code);

out    = zeros(hc,vn/vc);

for ii=1:hc
    xx=reshape(data(ii,:),vc,vn/vc);
    out(ii,:)= code(ii,:)*xx/vc;
end

%******************************** end of file ********************************

repo: OFDM


text:
%%%%%%%%%%Programmed by: Chi-Hang Kwan%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%Completion Date: December 13, 2012%%%%%%%%%%%%%%%%%%%

function posc_new = AI(window,pos,vel,posc,tstep)
tb = get(window,'UserData');

posf = pos + vel*tstep; %the position of the puck at the next frame
maxspeed = 1; %maximum moving speed of the computer's mallet

if posf(2) < 0
    pos_desire = [posf(1) tb.fieldh/4];
else
    if (posc(2)-posf(2)) > tb.r_mallet && abs(posf(1)-posc(1)) < tb.r_mallet %sweet spot to strike
        pos_desire = posf;
        if (posc(2)-posf(2))<(tb.r_mallet + tb.r_puck + 0.06)
            maxspeed = 2; %accelerate to hit puck when it is within 6cm of impact
        end
    else
        pos_desire(2) = posf(2) + tb.r_mallet+tb.r_puck + 0.03;
        if (posf(2)-posc(2))>=(tb.r_mallet-tb.r_puck)
            pos_desire(1) = posf(1) + sign(posc(1)-posf(1))*(tb.r_mallet+tb.r_puck+0.01);
        else
            pos_desire(1) = posf(1);
        end
    end
end

%%%%%%%%%%To ensure the mallet stays within its own half%%%%%%%%%%%%%
if abs(pos_desire(1))>tb.fieldw/2-tb.r_mallet
    pos_desire(1)= sign(pos_desire(1))*(tb.fieldw/2-tb.r_mallet);
end
if pos_desire(2) > tb.fieldh/2-tb.r_mallet
    pos_desire(2)= tb.fieldh/2-tb.r_mallet;
elseif pos_desire(2) < tb.r_mallet
    pos_desire(2) = tb.r_mallet;
end

%%%%%%%%%%To limit the speed of travel of the mallet%%%%%%%%%%%%%
delta = pos_desire - posc;
deltam = sqrt(delta(1)^2 + delta(2)^2);
if deltam == 0
    posc_new = posc;
else
    posc_new = posc + min(deltam,tstep*maxspeed)/deltam*delta;
end

author: Sable

repo: mcbench-benchmarks


text:
%MATHCONVENTIONS specifying mathematical conventions for IMAEL m-files
%
% Angles are counted in degrees, unless specified otherwise.
% Many old functions use radians, so it is advised to always handle angle
% data with care.
%

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function r8vec_mask_print ( n, a, mask_num, mask, title )

%*****************************************************************************80
%
%% R8VEC_MASK_PRINT prints a masked R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of components of the vector.
%
%    Input, real A(N), the vector to be printed.
%
%    Input, integer MASK_NUM, the number of masked elements.
%
%    Input, integer MASK(MASK_NUM), the indices of the vector to be printed.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Masked vector printout:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  for i = 1 : mask_num
    fprintf ( 1, '  %6d: %6d  %14g\n', i, mask(i), a(mask(i)) );
  end

  return
end


text:
function linplus_test1566 ( )

%*****************************************************************************80
%
%% TEST1566 tests R8BUT_MXV, R8BUT_PRINT, R8BUT_RANDOM, R8BUT_SL, R8BUT_VXM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  mu = 3;
  n = 10;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST1566\n' );
  fprintf ( 1, '  For a band matrix in upper triangular storage,\n' );
  fprintf ( 1, '  R8BUT_RANDOM sets a random value;\n' );
  fprintf ( 1, '  R8BUT_SL solves systems;\n' );
  fprintf ( 1, '  R8BUT_MXV computes matrix-vector products;\n' );
  fprintf ( 1, '  R8BUT_VXM computes vector-matrix products;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N =     %d\n', n );
  fprintf ( 1, '  Upper bandwidth MU = %d\n', mu );

  [ a, seed ] = r8but_random ( n, mu, seed );

  r8but_print ( n, mu, a, '  The R8BUT matrix:' );

  for job = 0 : 1
%
%  Set the desired solution.
%
    x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
    if ( job == 0 )
      b = r8but_mxv ( n, mu, a, x );
    else
      b = r8but_vxm ( n, mu, a, x );
    end

    r8vec_print ( n, b, '  The right hand side:' );
%
%  Solve the linear system.
%
    x = r8but_sl ( n, mu, a, b, job );

    if ( job == 0 )
      r8vec_print ( n, x, '  Solution to the untransposed system:' );
    else
      r8vec_print ( n, x, '  Solution to the transposed system:' );
    end

  end

  return
end


text:
% mattocell() - convert matrix to cell array
%
% Usage: >> C = mattocell( M );
%
% Author: Arnaud Delorme, CNL / Salk Institute, Jan 25 2002
%
% Note: this function overload the nnet toolbox function mattocell,
% but does not have all its capacities. You can delete the current
% function if you have the toolbox.

% Copyright (C) Jan 25 2002 Arnaud Delorme, CNL / Salk Institute
%
% 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 C = mattocell( M, varargin );

if nargin < 1
	help mattocell;
	return;
end;
if isempty(M)
	C = [];
	return;
end;
C = cell(size(M));
for i=1:size(M,1)
    for j=1:size(M,2)
        C{i,j} = M(i,j);
    end;
end;

author: PatternRecognition

repo: OpenBMI


text:
classdef InheritStateEnum < matlab.mixin.SetGet
    %InheritStateEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        InheritFromLastState('Inherit From Last State')
        InheritFromSpecifiedEvent('Inherit From Specified Event')
    end

    properties
        name
    end

    methods
        function obj = InheritStateEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('InheritStateEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('InheritStateEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('InheritStateEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function gX = matern32KernDiagGradX(kern, X)

% MATERN32KERNDIAGGRADX Gradient of MATERN32 kernel's diagonal with respect to X.
% FORMAT
% DESC computes the gradient of the diagonal of the matern kernel with nu=3/2 kernel matrix with
% respect to the elements of the design matrix given in X.
% ARG kern : the kernel structure for which gradients are being computed.
% ARG X : the input data in the form of a design matrix.
% RETURN gX : the gradients of the diagonal with respect to each element
% of X. The returned matrix has the same dimensions as X.
%
% SEEALSO : matern32KernParamInit, kernDiagGradX, matern32kernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2006

% KERN

gX = zeros(size(X));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function spec = fwDs_spec(hatRxx, f, d, tauGrid, c)

[nbin,nFrames] = size(hatRxx(:,:,1,1));
ngrid = length(tauGrid);
R11 = hatRxx(:,:,1,1);
R12 = hatRxx(:,:,1,2);
R22 = hatRxx(:,:,2,2);
TR = real(R11 + R22);
SINC = sinc(2*f*d/c);

SNR = zeros(nbin,nFrames,ngrid);
for pkInd=1:ngrid,
    EXP = repmat(exp(-2*1i*pi*tauGrid(pkInd)*f),1,nFrames);
    SNR(:,:,pkInd) = repmat(-(1+SINC)/2,1,nFrames) + repmat((1-SINC)/2,1,nFrames).*(TR + 2*real(R12.*EXP))./(TR - 2*real(R12.*EXP));
end
spec = SNR;

end

author: WenzheLiu-Speech


text:
%
% [idx,netsim,dpsim,expref]=apclusterSparse(s,p)
%
% APCLUSTER uses affinity propagation (Frey and Dueck, Science,
% 2007) to identify data clusters, using a set of real-valued
% pair-wise data point similarities as input. Each cluster is
% represented by a data point called a cluster center, and the
% method searches for clusters so as to maximize a fitness
% function called net similarity. The method is iterative and
% stops after maxits iterations (default of 500 - see below for
% how to change this value) or when the cluster centers stay
% constant for convits iterations (default of 50). The command
% apcluster(s,p,'plot') can be used to plot the net similarity
% during operation of the algorithm.
%
% For N data points, there may be as many as N^2-N pair-wise
% similarities (note that the similarity of data point i to k
% need not be equal to the similarity of data point k to i).
% These may be passed to APCLUSTER in an NxN matrix s, where
% s(i,k) is the similarity of point i to point k. In fact, only
% a smaller number of relevant similarities are needed for
% APCLUSTER to work. If only M similarity values are known,
% where M < N^2-N, they can be passed to APCLUSTER in an Mx3
% matrix s, where each row of s contains a pair of data point
% indices and a corresponding similarity value: s(j,3) is the
% similarity of data point s(j,1) to data point s(j,2).
%
% APCLUSTER automatically determines the number of clusters,
% based on the input p, which is an Nx1 matrix of real numbers
% called preferences. p(i) indicates the preference that data
% point i be chosen as a cluster center. A good choice is to
% set all preference values to the median of the similarity
% values. The number of identified clusters can be increased or
% decreased  by changing this value accordingly. If p is a
% scalar, APCLUSTER assumes all preferences are equal to p.
%
% The fitness function (net similarity) used to search for
% solutions equals the sum of the preferences of the the data
% centers plus the sum of the similarities of the other data
% points to their data centers.
%
% The identified cluster centers and the assignments of other
% data points to these centers are returned in idx. idx(j) is
% the index of the data point that is the cluster center for
% data point j. If idx(j) equals j, then point j is itself a
% cluster center. The sum of the similarities of the data
% points to their cluster centers is returned in dpsim, the
% sum of the preferences of the identified cluster centers is
% returned in expref and the net similarity (sum of the data
% point similarities and preferences) is returned in netsim.
%
% EXAMPLE
%
% N=100; x=rand(N,2); % Create N, 2-D data points
% M=N*N-N; s=zeros(M,3); % Make ALL N^2-N similarities
% j=1;
% for i=1:N
%   for k=[1:i-1,i+1:N]
%     s(j,1)=i; s(j,2)=k; s(j,3)=-sum((x(i,:)-x(k,:)).^2);
%     j=j+1;
%   end;
% end;
% p=median(s(:,3)); % Set preference to median similarity
% [idx,netsim,dpsim,expref]=apclusterSparse(s,p,'plot');
% fprintf('Number of clusters: %d\n',length(unique(idx)));
% fprintf('Fitness (net similarity): %f\n',netsim);
% figure; % Make a figures showing the data and the clusters
% for i=unique(idx)'
%   ii=find(idx==i); h=plot(x(ii,1),x(ii,2),'o'); hold on;
%   col=rand(1,3); set(h,'Color',col,'MarkerFaceColor',col);
%   xi1=x(i,1)*ones(size(ii)); xi2=x(i,2)*ones(size(ii));
%   line([x(ii,1),xi1]',[x(ii,2),xi2]','Color',col);
% end;
% axis equal tight;
%
% PARAMETERS
%
% [idx,netsim,dpsim,expref]=apclusterSparse(s,p,'NAME',VALUE,...)
%
% The following parameters can be set by providing name-value
% pairs, eg, apcluster(s,p,'maxits',1000):
%
%   Parameter    Value
%   'maxits'     Any positive integer. This specifies the
%                maximum number of iterations performed by
%                affinity propagation. Default: 1000.
%   'convits'    Any positive integer. APCLUSTER decides that
%                the algorithm has converged if the estimated
%                cluster centers stay fixed for convits
%                iterations. Increase this value to apply a
%                more stringent convergence test. Default: 100.
%   'dampfact'   A real number that is less than 1 and
%                greater than or equal to 0.5. This sets the
%                damping level of the message-passing method,
%                where values close to 1 correspond to heavy
%                damping which may be needed if oscillations
%                occur. Default: .9
%   'plot'       No value needed. This creates a figure that
%                plots the net similarity after each iteration
%                of the method. If the net similarity fails to
%                converge, consider increasing the values of
%                dampfact and maxits.
%   'details'    No value needed. This causes idx, netsim,
%                dpsim and expref to be stored after each
%                iteration.
%   'nonoise'    No value needed. Degenerate input similarities
%                (eg, where the similarity of i to k equals the
%                similarity of k to i) can prevent convergence.
%                To avoid this, APCLUSTER adds a small amount
%                of noise to the input similarities. This flag
%                turns off the addition of noise.
%
% Copyright (c) Brendan J. Frey and Delbert Dueck (2006). This
% software may be freely used and distributed for
% non-commercial purposes.

function [idx,netsim,dpsim,expref]=apclusterSparse(s,p,varargin);

% Handle arguments to function
if nargin<2 error('Too few input arguments');
else
    maxits=1000; convits=100; lam=0.9; plt=0; details=0; nonoise=0;
    i=1;
    while i<=length(varargin)
        if strcmp(varargin{i},'plot')
            plt=1; i=i+1;
        elseif strcmp(varargin{i},'details')
            details=1; i=i+1;
        elseif strcmp(varargin{i},'nonoise')
            nonoise=1; i=i+1;
        elseif strcmp(varargin{i},'maxits')
            maxits=varargin{i+1};
            i=i+2;
            if maxits<=0 error('maxits must be a positive integer'); end;
        elseif strcmp(varargin{i},'convits')
            convits=varargin{i+1};
            i=i+2;
            if convits<=0 error('convits must be a positive integer'); end;
        elseif strcmp(varargin{i},'dampfact')
            lam=varargin{i+1};
            i=i+2;
            if (lam<0.5)||(lam>=1)
                error('dampfact must be >= 0.5 and < 1');
            end;
        else i=i+1;
        end;
    end;
end;
if lam>0.9
    fprintf('\n*** Warning: Large damping factor in use. Turn on plotting\n');
    fprintf('    to monitor the net similarity. The algorithm will\n');
    fprintf('    change decisions slowly, so consider using a larger value\n');
    fprintf('    of convits.\n\n');
end;

% Check that standard arguments are consistent in size
if length(size(s))~=2 error('s should be a 2D matrix');
elseif length(size(p))>2 error('p should be a vector or a scalar');
elseif size(s,2)==3
    tmp=max(max(s(:,1)),max(s(:,2)));
    if length(p)==1 N=tmp; else N=length(p); end;
    if tmp>N
        error('data point index exceeds number of data points');
    elseif min(min(s(:,1)),min(s(:,2)))<=0
        error('data point indices must be >= 1');
    end;
elseif size(s,1)==size(s,2)
    N=size(s,1);
    if (length(p)~=N)&&(length(p)~=1)
        error('p should be scalar or a vector of size N');
    end;
else error('s must have 3 columns or be square'); end;

% Make vector of preferences
if length(p)==1 p=p*ones(N,1); end;

% Append self-similarities (preferences) to s-matrix
tmps=[repmat([1:N]',[1,2]),p]; s=[s;tmps];
M=size(s,1);

% In case user did not remove degeneracies from the input similarities,
% avoid degenerate solutions by adding a small amount of noise to the
% input similarities
if ~nonoise
    rns=randn('state'); randn('state',0);
    s(:,3)=s(:,3)+(eps*s(:,3)+realmin*100).*rand(M,1);
    randn('state',rns);
end;

% Construct indices of neighbors
ind1e=zeros(N,1);
for j=1:M k=s(j,1); ind1e(k)=ind1e(k)+1; end;
ind1e=cumsum(ind1e); ind1s=[1;ind1e(1:end-1)+1]; ind1=zeros(M,1);
for j=1:M k=s(j,1); ind1(ind1s(k))=j; ind1s(k)=ind1s(k)+1; end;
ind1s=[1;ind1e(1:end-1)+1];
ind2e=zeros(N,1);
for j=1:M k=s(j,2); ind2e(k)=ind2e(k)+1; end;
ind2e=cumsum(ind2e); ind2s=[1;ind2e(1:end-1)+1]; ind2=zeros(M,1);
for j=1:M k=s(j,2); ind2(ind2s(k))=j; ind2s(k)=ind2s(k)+1; end;
ind2s=[1;ind2e(1:end-1)+1];

% Allocate space for messages, etc
A=zeros(M,1); R=zeros(M,1); t=1;
if plt netsim=zeros(1,maxits+1); end;
if details
    idx=zeros(N,maxits+1);
    netsim=zeros(1,maxits+1);
    dpsim=zeros(1,maxits+1);
    expref=zeros(1,maxits+1);
end;

% Execute parallel affinity propagation updates
e=zeros(N,convits); dn=0; i=0;
while ~dn
    i=i+1;

    % Compute responsibilities
    for j=1:N
        ss=s(ind1(ind1s(j):ind1e(j)),3);
        as=A(ind1(ind1s(j):ind1e(j)))+ss;
        [Y,I]=max(as); as(I)=-realmax; [Y2,I2]=max(as);
        r=ss-Y; r(I)=ss(I)-Y2;
        R(ind1(ind1s(j):ind1e(j)))=(1-lam)*r+ ...
            lam*R(ind1(ind1s(j):ind1e(j)));
    end;

    % Compute availabilities
    for j=1:N
        rp=R(ind2(ind2s(j):ind2e(j)));
        rp(1:end-1)=max(rp(1:end-1),0);
        a=sum(rp)-rp;
        a(1:end-1)=min(a(1:end-1),0);
        A(ind2(ind2s(j):ind2e(j)))=(1-lam)*a+ ...
            lam*A(ind2(ind2s(j):ind2e(j)));
    end;

    % Check for convergence
    E=((A(M-N+1:M)+R(M-N+1:M))>0); e(:,mod(i-1,convits)+1)=E; K=sum(E);
    if i>=convits || i>=maxits
        se=sum(e,2);
        unconverged=(sum((se==convits)+(se==0))~=N);
        if (~unconverged&&(K>0))||(i==maxits) dn=1; end;
    end;

    % Handle plotting and storage of details, if requested
    if plt||details
        if K==0
            tmpnetsim=nan; tmpdpsim=nan; tmpexpref=nan; tmpidx=nan;
        else
            tmpidx=zeros(N,1); tmpdpsim=0;
            tmpidx(find(E))=find(E); tmpexpref=sum(p(find(E)));
            discon=0;
            for j=find(E==0)'
                ss=s(ind1(ind1s(j):ind1e(j)),3);
                ii=s(ind1(ind1s(j):ind1e(j)),2);
                ee=find(E(ii));
                if length(ee)==0 discon=1;
                else
                    [smx imx]=max(ss(ee));
                    tmpidx(j)=ii(ee(imx));
                    tmpdpsim=tmpdpsim+smx;
                end;
            end;
            if discon
                tmpnetsim=nan; tmpdpsim=nan; tmpexpref=nan; tmpidx=nan;
            else tmpnetsim=tmpdpsim+tmpexpref;
            end;
        end;
    end;
    if details
        netsim(i)=tmpnetsim; dpsim(i)=tmpdpsim; expref(i)=tmpexpref;
        idx(:,i)=tmpidx;
    end;
    if plt
        netsim(i)=tmpnetsim;
        figure(234);
        tmp=1:i; tmpi=find(~isnan(netsim(1:i)));
        plot(tmp(tmpi),netsim(tmpi),'r-');
        xlabel('# Iterations');
        ylabel('Net similarity of quantized intermediate solution');
        drawnow;
    end;
end;
% Identify  exemplars
E=((A(M-N+1:M)+R(M-N+1:M))>0); K=sum(E);
if K>0
    tmpidx=zeros(N,1); tmpidx(find(E))=find(E); % Identify clusters
    for j=find(E==0)'
        ss=s(ind1(ind1s(j):ind1e(j)),3);
        ii=s(ind1(ind1s(j):ind1e(j)),2);
        ee=find(E(ii));
        [smx imx]=max(ss(ee));
        tmpidx(j)=ii(ee(imx));
    end;
    EE=zeros(N,1);
    for j=find(E)'
        jj=find(tmpidx==j); mx=-Inf;
        ns=zeros(N,1); msk=zeros(N,1);
        for m=jj'
            mm=s(ind1(ind1s(m):ind1e(m)),2);
            msk(mm)=msk(mm)+1;
            ns(mm)=ns(mm)+s(ind1(ind1s(m):ind1e(m)),3);
        end;
        ii=jj(find(msk(jj)==length(jj)));
        [smx imx]=max(ns(ii));
        EE(ii(imx))=1;
    end;
    E=EE;
    tmpidx=zeros(N,1); tmpdpsim=0;
    tmpidx(find(E))=find(E); tmpexpref=sum(p(find(E)));
    for j=find(E==0)'
        ss=s(ind1(ind1s(j):ind1e(j)),3);
        ii=s(ind1(ind1s(j):ind1e(j)),2);
        ee=find(E(ii));
        [smx imx]=max(ss(ee));
        tmpidx(j)=ii(ee(imx));
        tmpdpsim=tmpdpsim+smx;
    end;
    tmpnetsim=tmpdpsim+tmpexpref;
else
    tmpidx=nan*ones(N,1); tmpnetsim=nan; tmpexpref=nan;
end;
if details
    netsim(i+1)=tmpnetsim; netsim=netsim(1:i+1);
    dpsim(i+1)=tmpnetsim-tmpexpref; dpsim=dpsim(1:i+1);
    expref(i+1)=tmpexpref; expref=expref(1:i+1);
    idx(:,i+1)=tmpidx; idx=idx(:,1:i+1);
else
    netsim=tmpnetsim; dpsim=tmpnetsim-tmpexpref;
    expref=tmpexpref; idx=tmpidx;
end;
if plt||details
    fprintf('\nNumber of identified clusters: %d\n',K);
    fprintf('Fitness (net similarity): %f\n',tmpnetsim);
    fprintf('  Similarities of data points to exemplars: %f\n',dpsim(end));
    fprintf('  Preferences of selected exemplars: %f\n',tmpexpref);
    fprintf('Number of iterations: %d\n\n',i);
end;
if unconverged
    fprintf('\n*** Warning: Algorithm did not converge. The similarities\n');
    fprintf('    may contain degeneracies - add noise to the similarities\n');
    fprintf('    to remove degeneracies. To monitor the net similarity,\n');
    fprintf('    activate plotting. Also, consider increasing maxits and\n');
    fprintf('    if necessary dampfact.\n\n');
end;

repo: cluster


text:
function m = gpComputeM(model)

% GPCOMPUTEM Compute the matrix m given the model.
% FORMAT
% DESC computes the matrix m (the scaled, bias and mean function
% removed matrix of the targets), given the model.
% ARG model : the model for which the values are to be computed.
% RETURN m : the scaled, bias and mean function removed values.
%
% SEEALSO : gpCreate, gpComputeAlpha, gpUpdateAD
%
% COPYRIGHT : Neil D. Lawrence, 2006

% GP

% Remove mean function value from m (if mean function present).
if isfield(model, 'meanFunction') && ~isempty(model.meanFunction)
  m = model.y - modelOut(model.meanFunction, model.X);
else
  m = model.y;
end

% Remove bias and apply scale.
for i = 1:model.d
  m(:, i) = m(:, i) - model.bias(i);
  if model.scale(i)
    m(:, i) = m(:, i)/model.scale(i);
  end
end


text:
% Matlab Toolbox for Dimensionality Reduction
% Version 0.8 18-APR-2012
%
%
%Main functions
%----------------------
% compute_mapping    Performs dimension reduction using the specified technique
% out_of_sample      Performs out-of-sample extension using trained DR technique
% out_of_sample_est  Performs out-of-sample extension by approximation
% reconstruct_data   Reconstruct data from low-dimensional representation
% intrinsic_dim      Estimates intrinsic dimensionality of the data
%
%Helper functions
%----------------------
% drgui              Shows a GUI that facilitates access to all functions
% generate_data      Generates artificial data manifolds
% prewhiten          Performs whitening of a data set
%
%Installation functions
%----------------------
% mexall             Compiles all C++ code in the toolbox
% test_toolbox       Tests all functionalities of the toolbox
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology


text:
function [t,conse,allt,glm] = designsim(noisetype,dmodel,HRF,ISI,TR,noise_var,c,beta,S,varargin)
% function [t,modelse,allt,glm] = designsim(noisetype,dmodel,HRF,ISI,TR,noise_var,c,betas,S,xcfunction [if using 'myxc'],num2avgover)
%  [t,modelse,allt,glm] = designsim('myxc',M.modelatTR,HRF,ISI,TR,noise_var,c,beta,fullS,xc)
%
% noisetype:
%   '1overf' : use Luis Hernandez' 1/f model with your specified noise variance
%   'myxc'   : use your own autocorrelation function and noise variance
%
% model = matrix of regressors (cols) sampled at frequency of ISI
% noise_var = variance of 1/f noise from your scanner
% c = a contrast across the regressors
% beta = a list of betas to test
%   each beta should be a row vector, and you can have multiple rows
%
% S: smoothing matrix, or 0 for no smoothing.
% _______________________________________________________________________
% Output is a distribution of t-scores in a vector and the maximum correlation between predictors
% t-scores are in a column vector, one for each set of betas.
%
% Special mode: if 4 output arguments specified, 'plot mode'.  Produces graph and extended output.
% outputs:
%	t 		mean of t-scores for specified number of noise vectors tested, numnoise in script
%	allt	vector of all t values generated for each noise vector
%
% specifying a 4th output puts designsim.m into 'verbose' mode.
%
% Last modified   5/24/01 Tor Wager		add noise first, then smooth and filter, etc.

% if no stimuli for one condition, gives error - create extra col of zeros if this happens.
while size(beta,2) > size(dmodel,2)
 	warning('      Designsim.m : dmodel has too few columns or beta is too long! No stim. in one cond.?')
	whos dmodel
     	dmodel = [dmodel zeros(size(dmodel,1),1)];
end

if isempty(c)
	if nargout > 3,disp('	...designsim.m: No contrasts found.  Using first predictor for one sample T-test.'),end
	c = zeros(1,size(dmodel,2));
	c(1) = 1;
end

c = c';

done = 0;
%if nargout > 3,numnoise = 1;,else numnoise = 1;,end
% if average over noise models is given, return the average - otherwise return each individual t and se score.
if nargin > 10, numnoise = varargin{2};,else numnoise = 1;,end

% make the noise vectors
% =====================================================================================
switch noisetype
case '1overf'
   noise = make1overf(size(dmodel,1), 1/ISI) * sqrt(noise_var);
case 'myxc'

    for i = 1:numnoise
        noise(:,i) = noisevector(size(dmodel,1),varargin{1},noise_var)';
    end
otherwise error('unsupported noise type.')
end

   		% Create the response data by weighting all the regressors by a beta parameter
   		% adding all the regressors together
   		% and adding noise to the result
 beta = beta';
 if nargout > 3,
    	%figure;,
    	%disp('testing with only one noise model, returning 1st beta vector only in glm.')
    	%beta = beta(:,1); noise = noise(:,1);
	c = c(:,1);	% test first contrast here - this should be changed back for normal use. this is for rundesignsimMAP
 else c = c(:,1); % tests only first contrast if you specify many iterations.
 end

% make the data vectors and fit
% =====================================================================================
for i = 1:size(beta,2)

	if nargout > 3
        	subplot(4,size(beta,2),3*i-2);hold on
        	plot(dmodel * beta(:,i),'r');title(['ideal model, betas = ' num2str(beta(:,i)')]);drawnow
        	clear glm
	end

	% make data
	% --------------------------------------------------
	for j = 1:size(noise,2)

		try
			myfit = dmodel * beta(:,i);
		catch
			disp(['	designsim.m: dmodel and beta are not the same size - wrong beta size?'])
			whos dmodel
			beta
			error('exiting.')
		end
		myfit = (myfit - mean(myfit));		% mean center so betas are not 'correlated' with intercept.
		data(:,j) = myfit + noise(:,j);

		% testing stuff -
		%disp(['designsim.m: beta = '])
		%beta(:,i)
		%figure;subplot(3,1,1);plot(noise(:,j))
		%subplot(3,1,2);plot(dmodel * beta(:,i))
		%subplot(3,1,3);plot(data(:,j));
	end

	% smoothing
	% --------------------------------------------------
	if S,
		dmodel = S * dmodel;
		data = S * data;
		glm.smoothed = 'yes';
	end

	% fit model to data and save t and se
	% --------------------------------------------------
	for j = 1:size(noise,2)
        	[myt,myconse] = my_glm(dmodel,data(:,j),c);		% fit model
        	allt(i,j) = myt(1);
		allse(i,j) = myconse(1);

		% make glm structure if necessary
		% ---------------------------------------------
        	if nargout > 3
            		glm.y = data;
            		glm.X = [dmodel ones(size(dmodel,1),1)];
            		glm.betas = (glm.X \ glm.y)';
            		glm.fit = glm.X * glm.betas';
            		glm.e = glm.y - glm.fit;
            		glm.xtxi = inv(glm.X' * glm.X);
            		glm.df = (size(glm.y,1) - size(glm.X,2));
            		glm.evar = glm.e' * glm.e / glm.df;
            		glm.se = sqrt(diag(glm.xtxi .* glm.evar))';
            		glm.t = (glm.betas ./ glm.se);
            		glm.p = tdist(glm.t,glm.df);
            		glm.c = [c;zeros(1,size(c,2))];
            		glm.con_se = diag(sqrt(glm.evar .* (glm.c' * glm.xtxi * glm.c)))';
            		glm.con_t = glm.betas * glm.c ./ glm.con_se;
            		glm.con_p = tdist(glm.con_t,glm.df);

	    		if j == 1							% plot the first noise model.
        			subplot(4,size(beta,2),3*i-1);hold on
        			plot(data(:,j));hold on;plot(glm.fit,'r--'); title(['Data with noise added and fit']);drawnow
        			subplot(4,size(beta,2),3*i);hold on
        			plot(data(:,j) - glm.fit); title('Residuals');drawnow
 	    		end
        	end

	end
	% loop thru noise models.
 end	% loop thru betas.

t = mean(allt,2);
conse = mean(allse,2);

return

% rmat = abs(corrcoef(dmodel));%rmat(rmat == 1) = 0;		% zero diagonals%r = max(max(rmat));


text:
function varargout = subsref(pf,s)
% overide polefigure() and polefigure{}

if any(strcmp(s(1).type,{'()','{}'}))

  if strcmp(s(1).type,'{}') || iscell(s(1).subs{1}) || ...
      (ischar(s(1).subs{1}) && ~strcmp(s(1).subs{1},':')) || isa(s(1).subs{1},'Miller')

    pf = pf.select(s(1).subs);

  else

    id = false(size(pf.intensities));

    if isa(s(1).subs{1},'vector3d')
      id(find(pf.r,s(1).subs{1}))=true;
    else
      id = subsasgn(id,s(1),true);
    end

    cs = cumsum([0,cellfun('prodofsize',pf.allI)]);

    for k = 1:length(pf.allI)

      idi = id(cs(k)+1:cs(k+1));

      pf.allI{k} = pf.allI{k}(idi);
      pf.allR{k} = pf.allR{k}(idi);
    end
  end

  % is there something more to do?
  if numel(s)>1
    s = s(2:end);
  else
    varargout{1} = pf;
    return
  end
end

% maybe reference to a dynamic property
try %#ok<TRYNC>
  [varargout{1:nargout}] = subsref@dynProp(pf,s);
  return
end

% maybe reference to a dynamic option
try %#ok<TRYNC>
  [varargout{1:nargout}] = subsref@dynOption(pf,s);
  return
end

end


text:
% Add multiple edges in an edge list
%
% INPUTS: original (non-compact) edge list
% OUTPUTS: final compact edge list (no row repetitions)
%
% Example: [1 2 2; 2 2 1; 4 5 1] -> [1 2 3; 4 5 1]
% GB: last updated, Sep 25 2012

function elc=addEdgeWeights(el)

el2=[el(:,1), el(:,2)]; % make the edge list searchable w/o the weights
visited=[];             % mark visited edges

elc=[];
for e=1:size(el,1)
    if sum(ismember(visited,el2(e,:),'rows'))==0  % if not visited yet
        ind=ismember(el2,el2(e,:),'rows');
        ind=find(ind==1);     % these are all the ocurrences of el(e,:)
        elc=[elc; el(e,1), el(e,2), sum(el(ind,3))];
        visited=[visited; el2(e,:)];
    end
end

author: aeolianine

repo: octave-networks-toolbox


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcMagnify(Temp,Event,handles)
handles = guidata(handles.MU_matrix_display);

MSize = inputdlg('Please specify the size of magnified region (pixels).','Specify Size',1,{'30'});
if isempty(MSize)
    warndlg('Magnifier was cancelled.');
    return;
end

try
    handles.MSize = str2double(MSize{1});
    A=handles.BMatrix(1:handles.MSize,1:handles.MSize); % quick test for size
    handles.MSize = max(2,handles.MSize - mod(handles.MSize,2)); % even value for size, minimum 2
    figure('Name', 'Instant Magnifier');
    imagesc(zeros(handles.MSize,handles.MSize));
catch me
    errordlg('Input size is invalid.');
    return;
end
handles.MHandle = get(gca,'Children');
colormap(handles.V.Color_map);
axis off;

guidata(handles.MU_matrix_display, handles);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [c]=ref_dgt(f,g,a,M)
%REF_DGT  Reference Discrete Gabor transform.
%   Usage:  c=ref_dgt(f,g,a,M);
%
%   Linear algebra version of the algorithm. Create big matrix
%   containing all the basis functions and multiply with the transpose.

g = double(g);
f = double(f);
% Calculate the parameters that was not specified.
L=size(g,1);

b=L/M;
N=L/a;
W=size(f,2);
R=size(g,2);

% Create 2x2 grid matrix..
V=[a,0;
   0,b];

% Create lattice and Gabor matrix.
lat=ref_lattice(V,L);
G=ref_gaboratoms(g,lat);

% Apply matrix to f.
c=G'*f;

% reshape to correct output format.

c=reshape(c,M,N,R*W);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function m = copymatrix(varargin)

% function m = copymatrix(m,sub,val)
%
% <m> is a matrix.  can be [] when <sub> is a logical, in which case
%   we assume <m> is zeros(size(<sub>)).
% <sub> is
%   (1) some sort of index (e.g. vector of indices, logical mask)
%   (2) a function that accepts <m> and outputs an index
% <val> is something that can be assigned to the <sub> indices
%
% return a copy of <m> that has <val> shoved into <sub>.
% this function is useful for making modifications of a matrix on-the-fly.
%
% example:
% imagesc(copymatrix(randn(10,10),rand(10,10)>.5,0));
% isequal(copymatrix([1 2 3],@(x) x > 1,1),[1 1 1])
%
% OR
%
% function m = copymatrix(m,sub,dim,val)
%
% <m> is a matrix
% <sub> is a vector of indices or a vector of logicals
% <dim> is a dimension of <m>
% <val> is something that can be assigned to the <sub> indices of <m>,
%   assuming ':' for all other dimensions
%
% return a copy of <m> that has <val> shoved into <sub>.
% this function is useful for making modifications of a matrix on-the-fly.
%
% example:
% copymatrix([1 2 3; 4 5 6],2,1,[1 2 3])

if nargin==3
  m = varargin{1};
  sub = varargin{2};
  val = varargin{3};
  if isa(sub,'logical') && isempty(m)
    m = zeros(size(sub));
  end
  if isa(sub,'function_handle')
    m(feval(sub,m)) = val;
  else
    m(sub) = val;
  end
else
  m = varargin{1};
  sub = varargin{2};
  dim = varargin{3};
  val = varargin{4};
  ix = repmat({':'},[1 max(ndims(m),dim)]);
  ix{dim} = sub;
  m(ix{:}) = val;
end

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function file_list=get_file_list(directory)
%file_list=get_file_list(directory)
%returns a filelist cell{}, if no directory is given, '.' is used.
if nargin<1
  directory='.';
end

f_struct=dir(directory);
file_list={};
j=1;
for i=1:length(f_struct)
  if ~strcmp(f_struct(i).name, '.') & ~strcmp(f_struct(i).name, ...
					      '..') & ~strcmp(f_struct(i).name, '.svn')
    file_list{j}=f_struct(i).name;
    j=j+1;
  end %if
end %for i
end

author: thomaskuestner


%**************************************************************************
% Q Learning applied to Cart-Pole balancing problem.
% The environment of the learning system is a black box, from which it has
% several lines and a reinforcement line. Its task is to learn to give
% responses which maximize the scalar signals on its reinforcement line.
% The Q-value update is the following, if the system takes action a from
%   state s at time t, and arrives at state ss with feedback r at time t+1:
%
%   Q(t+1, s, a) = Q(t, s, a)
%                  + alpha (r + gamma max_{b}Q(t,ss, b) - Q(t, s, a))

% get_box:  Given the current state, returns a number from 1 to 162
%           designating the region of the state space encompassing the current
%           state.
%           Returns a value of -1 if a failure state is encountered.

% cart_pole: The cart and pole dynamics; given the Force and
%            current state, estimates next state

%Code written by: Savinay Nagendra
%email id:        nagsavi17@gmail.com
%**************************************************************************
clc;
clear all;
close all;
% Initialization

NUM_BOXES = 163;
ALPHA = 0.4;             % Learning rate parameter
GAMMA = 0.999;           % Discount factor for future reinf
Q = zeros(NUM_BOXES,2);  % State-Action Values
action = [30 -30];
MAX_FAILURES = 1000;
MAX_STEPS = 150000;
epsilon = 0;
steps = 0;
failures = 0;
thetaPlot = 0;
xPlot = 0;
%Pendulum state initialization
theta = 0;
thetaDot = 0;
x = 0;
xDot = 0;
box = getBox4(theta,thetaDot,x,xDot);

while(steps<=MAX_STEPS && failures<+MAX_FAILURES)
    steps = steps + 1;

    % choose either explore or exploit
    if(rand>epsilon)       % exploit
        [~,actionMax] = max(Q(box,:));
        currentAction = action(actionMax);
    else                   % explore
        currentAction = datasample(action,1);
    end

    actionIndex = find(action == currentAction); % index of chosen action
    %Apply action to the simulated cart pole
    [thetaNext,thetaDotNext,thetaacc,xNext,xDotNext] = cart_pole2(currentAction,theta,thetaDot,x,xDot);
    %Get box of state space containing the resulting state
    thetaPlot(end + 1) = thetaNext;
    xPlot(end + 1) = xNext;
    newBox = getBox4(thetaNext,thetaDotNext,xNext,xDotNext);
    theta = thetaNext;
    thetaDot = thetaDotNext;
    x = xNext;
    xDot = xDotNext;
    if(newBox==163)
        r = -1;
        Q(newBox,:) = 0;
        figure(2);
        plot((1:length(thetaPlot)),thetaPlot,'-ob');
        figure(3);
        plot((1:length(xPlot)),xPlot,'-og');

        thetaPlot = 0;
        xPlot = 0;
        theta = 0;
        thetaDot = 0;
        x = 0;
        xDot = 0;
        newBox = getBox4(theta,thetaDot,x,xDot);
        failures = failures + 1;
        fprintf('Trial %d was %d steps. \n',failures,steps);
        figure(1);
        plot(failures,steps,'-or');
        hold on;
        steps = 0;
    else
        r = 0;
    end
    Q(box,actionIndex) = Q(box,actionIndex) + ALPHA*(r + GAMMA*max(Q(newBox,:)) - Q(box,actionIndex));
    box = newBox;
end
if(failures == MAX_FAILURES)
    fprintf('Pole not balanced. Stopping after %d failures.',failures);
else
    fprintf('Pole balanced successfully for at least %d steps\n', steps);
    figure(1);
    plot(failures+1,steps,'-or');
    hold on;
    figure(2);
    plot((1:length(thetaPlot)),thetaPlot,'-ob');
    figure(3);
    plot((1:length(xPlot)),xPlot,'-og');
    figure(4);
    plot((1:301),thetaPlot(1:301),'-ob');
    hold on;
    figure(5);
    plot((1:301),xPlot(1:301),'-og');
    hold on;
end


text:
function value = diaedg ( x0, y0, x1, y1, x2, y2, x3, y3 )

%*****************************************************************************80
%
%% DIAEDG chooses a diagonal edge.
%
%  Discussion:
%
%    The routine determines whether 0--2 or 1--3 is the diagonal edge
%    that should be chosen, based on the circumcircle criterion, where
%    (X0,Y0), (X1,Y1), (X2,Y2), (X3,Y3) are the vertices of a simple
%    quadrilateral in counterclockwise order.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2005
%
%  Author:
%
%    Original FORTRAN77 version by Barry Joe,
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Barry Joe,
%    GEOMPACK - a software package for the generation of meshes
%    using geometric algorithms,
%    Advances in Engineering Software,
%    Volume 13, pages 325-331, 1991.
%
%  Parameters:
%
%    Input, real X0, Y0, X1, Y1, X2, Y2, X3, Y3, the
%    coordinates of the vertices of a quadrilateral, given in
%    counter clockwise order.
%
%    Output, integer VALUE, chooses a diagonal:
%    +1, if diagonal edge 02 is chosen;
%    -1, if diagonal edge 13 is chosen;
%     0, if the four vertices are cocircular.
%
  tol = 100.0 * eps;

  dx10 = x1 - x0;
  dy10 = y1 - y0;
  dx12 = x1 - x2;
  dy12 = y1 - y2;
  dx30 = x3 - x0;
  dy30 = y3 - y0;
  dx32 = x3 - x2;
  dy32 = y3 - y2;

  tola = tol * max ( abs ( dx10 ), ...
               max ( abs ( dy10 ), ...
               max ( abs ( dx30 ), abs ( dy30 ) )));

  tolb = tol * max ( abs ( dx12 ), ...
               max ( abs ( dy12 ), ...
               max ( abs ( dx32 ), abs ( dy32 ) )));

  ca = dx10 * dx30 + dy10 * dy30;
  cb = dx12 * dx32 + dy12 * dy32;

  if ( tola < ca & tolb < cb )

    value = -1;

  elseif ( ca < -tola & cb < -tolb )

    value = 1;

  else

    tola = max ( tola, tolb );
    s = ( dx10 * dy30 - dx30 * dy10 ) * cb ...
      + ( dx32 * dy12 - dx12 * dy32 ) * ca;

    if ( tola < s )
      value = -1;
    elseif ( s < -tola )
      value = 1;
    else
      value = 0;
    end

  end

  return
end


text:
function color = isColorStack(img)
%ISCOLORSTACK  Check if a 3D stack is color or gray-scale
%
%   COL = isColorStack(IMG)
%   Returns TRUE if the 3D image has a color dimension, or FASLE if the
%   image is gray scale image.
%
%   This function was mainly created for simplifying the listings of other
%   functions in the same package.
%
%   Example
%   img = analyze75read(analyze75info('brainMRI.hdr'));
%   isColorStack(img)
%   ans =
%       0
%
%
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-07-02,    using Matlab 7.9.0.529 (R2009b)
% http://www.pfl-cepia.inra.fr/index.php?page=slicer
% Copyright 2010 INRA - Cepia Software Platform.

dim = size(img);
color = length(dim)>3;


text:
function Res = PlotCube(hAxes,Cube)
% plot the cube Cube in axes hAxes
axes(hAxes); % current axes for draw
cla
hold on
for k=1:27,
  CurrNodes = Cube.Nodes(:,:,k);
  CurrColor = Cube.Color(:,:,k);
  patch('Vertices',CurrNodes,'Faces',Cube.Order,...
    'FaceVertexCData',CurrColor,'FaceColor','flat',...
    'ButtonDownFcn','Cube=RotateWithMouse(hAxes,Cube);');
end
hold off
return


text:
function dUpsilon = lfmvpGradientUpsilonVector(gamma, sigma2, t, upsilon)

% LFMVPGRADIENTUPSILONVECTOR Gradient upsilon vector vel. pos.
% FORMAT
% DESC computes the gradient of a portion of the LFM kernel.
% ARG gamma : Gamma value for system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG upsilon : precomputation of the upsilon matrix.
% ARG mode : operation mode, according to the derivative (mode 0,
% derivative wrt t1, mode 1 derivative wrt t2)
% RETURN upsilon : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio A. Alvarez, 2010
%
% SEEALSO : lfmvpComputeUpsilonMatrix.m

% KERN

if nargin<4
    upsilon = lfmComputeUpsilonVector(gamma, sigma2, t);
end

dUpsilon = -upsilon - gamma*lfmGradientUpsilonVector(gamma, sigma2, t);


text:
function [out] = VBA_spm_eeg_displayECD(Pos,Orient,Var,Names,options)
% Plot dipole positions onto the SPM canonical mesh
% FORMAT [out] = spm_eeg_displayECD(Pos,Orient,Var,Names,options)
%
% IN (admissible choices):
%   - Pos: a 3xndip matrix containing the positions of the dipoles in
%   the canonical frame of reference
%   - Orient: the same with dipole orientations
%   - Var: the same with position variance
%   - Names: the same with dipole names
%   - options: an optional structure containing
%       .hfig: the handle of the display figure
%       .tag: the tag to be associated with the created UI objects
%       .add: binary variable ({0}, 1: just add dipole in the figure .hfig)
%
% OUT:
%   - out: a structure containing the handles of the object in the figure
%   (including the mesh, the dipoles, the transparency slider, etc...)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Jean Daunizeau
% $Id: spm_eeg_displayECD.m 5737 2013-11-10 20:23:49Z karl $

% checks and defaults
%--------------------------------------------------------------------------
hfig       = [];
ParentAxes = [];
query      = [];
handles    = [];
tag        = '';
try, options; catch, options = [];      end
try, hfig        = options.hfig;        end
try, tag         = options.tag;         end
try, ParentAxes  = options.ParentAxes;  end
try, query       = options.query;       end
try, handles     = options.handles;     end
try
    figure(hfig);
catch
    hfig  = VBA_spm_figure('GetWin','Graphics');
    VBA_spm_figure('Clear',hfig);
    ParentAxes = axes('parent',hfig);
end
try
    markersize = options.markersize;
catch
    markersize = 20;
end
try
    meshsurf = options.meshsurf;
catch
    meshsurf = fullfile(spm('Dir'),'canonical','cortex_5124.surf.gii');
end

if isscalar(Var), Var = Pos*0 + Var^2;   end
try, Pos{1};    catch, Pos = {Pos};      end
try, Orient{1}; catch, Orient = {Orient};end
try, Var{1};    catch, Var = {Var};      end

ndip = size(Pos{1},2);
if ~exist('Names','var') || isempty(Names)
    for i=1:ndip
        Names{i} = num2str(i);
    end
end

col = ['b','g','r','c','m','y','k','w'];
tmp = ceil(ndip./numel(col));
col = repmat(col,1,tmp);
pa  = get(ParentAxes,'position');

if ndip > 0

    if isempty(query)
        opt.hfig = hfig;
        opt.ParentAxes = ParentAxes;
        opt.visible = 'off';
        pos2 = [pa(1),pa(2)+0.25*pa(4),0.03,0.5*pa(4)];
        out  = VBA_spm_eeg_render(meshsurf,opt);
        handles.mesh = out.handles.p;
        handles.BUTTONS.transp = out.handles.transp;
        handles.hfig = out.handles.fi;
        handles.ParentAxes  = out.handles.ParentAxes;
        set(handles.mesh,...
            'facealpha',0.1,...
            'visible','on')
        set(handles.BUTTONS.transp,...
            'value',0.1,...
            'position',pos2,...
            'visible','on')
    end

    set(ParentAxes,'nextplot','add')
    for j=1:length(Pos)
        for i =1:ndip
            try
                set(handles.hp(j,i),...
                    'xdata',Pos{j}(1,i),...
                    'ydata',Pos{j}(2,i),...
                    'zdata',Pos{j}(3,i));
            catch
                handles.hp(j,i) = plot3(handles.ParentAxes,...
                    Pos{j}(1,i),Pos{j}(2,i),Pos{j}(3,i),...
                    [col(i),'.'],...
                    'markerSize',markersize,...
                    'visible','off');
            end
            try
                no = sqrt(sum(Orient{j}(:,i).^2));
                if no > 0
                    Oi = 1e1.*Orient{j}(:,i)./no;
                else
                    Oi = 1e-5*ones(3,1);
                end
                try
                    set(handles.hq(j,i),...
                        'xdata',Pos{j}(1,i),...
                        'ydata',Pos{j}(2,i),...
                        'zdata',Pos{j}(3,i),...
                        'udata',Oi(1),...
                        'vdata',Oi(2),...
                        'wdata',Oi(3))
                catch
                    handles.hq(j,i) = quiver3(handles.ParentAxes,...
                        Pos{j}(1,i),Pos{j}(2,i),Pos{j}(3,i),...
                        Oi(1),Oi(2),Oi(3),col(i),...
                        'lineWidth',2,'visible','off');
                end
                if isequal(query,'add')
                    set(handles.hq(j,i),...
                        'LineStyle','--',...
                        'lineWidth',1)
                end
            end
            [x,y,z]= ellipsoid(Pos{j}(1,i),Pos{j}(2,i),Pos{j}(3,i),...
                1.*sqrt(Var{j}(1,i)),1.*sqrt(Var{j}(2,i)),1.*sqrt(Var{j}(1,i)),20);
            try
                set(handles.hs(j,i),...
                    'xdata',x,...
                    'ydata',y,...
                    'zdata',z);
            catch
                handles.hs(j,i) = surf(handles.ParentAxes,...
                    x,y,z,...
                    'edgecolor','none',...
                    'facecolor',col(i),...
                    'facealpha',0.2,...
                    'visible','off');
            end
            try
                set(handles.ht(j,i),...
                    'position',Pos{j}(:,i));
            catch
                handles.ht(j,i) = text(...
                    Pos{j}(1,i),Pos{j}(2,i),Pos{j}(3,i),...
                    Names{i},...
                    'Parent',handles.ParentAxes,...
                    'visible','off');
            end
        end
    end

    if length(Pos) > 1

        try, set(handles.hp(end,:),'visible','on'); end
        try, set(handles.hq(end,:),'visible','on'); end
        try, set(handles.hs(end,:),'visible','on'); end
        try, set(handles.ht(end,:),'visible','on'); end

        handles.uic(1) = uicontrol(handles.fi,...
            'units','normalized',...
            'position',[0.45,0.5,0.2,0.03],...
            'style','radio','string','Show priors',...
            'callback',@doChange1,...
            'BackgroundColor',[1 1 1],...
            'tooltipstring','Display prior locations',...
            'userdata',handles,'value',0,...
            'BusyAction','cancel',...
            'Interruptible','off',...
            'tag','plotEEG');

        handles.uic(2) = uicontrol(handles.fi,...
            'units','normalized',...
            'position',[0.45,0.53,0.2,0.03],...
            'style','radio','string','Show posteriors',...
            'callback',@doChange2,...
            'BackgroundColor',[1 1 1],...
            'tooltipstring','Display posterior locations',...
            'userdata',handles,'value',1,...
            'BusyAction','cancel',...
            'Interruptible','off',...
            'tag','plotEEG');

    else

        try, set(handles.hp(1,:),'visible','on'); end
        try, set(handles.hq(1,:),'visible','on'); end
        try, set(handles.hs(1,:),'visible','on'); end
        try, set(handles.ht(1,:),'visible','on'); end

    end

end

try
    clear out
    out.handles = handles;
catch
    out = [];
end

%==========================================================================
function doChange1(i1,i2)
val = get(i1,'value');
handles = get(i1,'userdata');
if ~val
    try, set(handles.hp(1,:),'visible','off'); end
    try, set(handles.hq(1,:),'visible','off'); end
    try, set(handles.hs(1,:),'visible','off'); end
    try, set(handles.ht(1,:),'visible','off'); end
else
    try, set(handles.hp(1,:),'visible','on');  end
    try, set(handles.hq(1,:),'visible','on');  end
    try, set(handles.hs(1,:),'visible','on');  end
    try, set(handles.ht(1,:),'visible','on');  end
end

%==========================================================================
function doChange2(i1,i2)
val = get(i1,'value');
handles = get(i1,'userdata');
if ~val
    try, set(handles.hp(2,:),'visible','off'); end
    try, set(handles.hq(2,:),'visible','off'); end
    try, set(handles.hs(2,:),'visible','off'); end
    try, set(handles.ht(2,:),'visible','off'); end
else
    try, set(handles.hp(2,:),'visible','on');  end
    try, set(handles.hq(2,:),'visible','on');  end
    try, set(handles.hs(2,:),'visible','on');  end
    try, set(handles.ht(2,:),'visible','on');  end
end

author: MBB-team

repo: VBA-toolbox


text:
% SIGNALSTAT  -  Computes and plots statistical characteristics of a signal,
%                  including the data histogram, a fitted normal distribution,
%                  a normal distribution fitted on trimmed data, a boxplot, and
%                  the QQ-diagram. The estimates value are printed in a panel and
%                  can be read as output. Optionally, a topographic map (see TOPOPLOT)
%                  can be plotted.
%                  The boxplot and the Kolmogorov-Smirnov test require the
%                  MATLAB Statistics Toolbox.
%
% Usage:
%   >>  signalstat( data )
%   >>  signalstat( data, plotlab, dlabel, percent );
%   >>  [M,SD,sk,k,med,zlow,zhi,tM,tSD,tndx,ksh] = ...
%                 signalstat( data, plotlab, dlabel, percent, dlabel2, map, chan_locs );
%
% Inputs:
%   data        - data vector
%
% Optional inputs:
%   plotlab     - 1: default->plot  |  0: ->no plot
%   dlabel      - A label for the data ([]: default->'Potential [V]')
%   percent     - percentage of data to exclude for trimmed mean & SD ([]:default->5)
%                 Excluded is 'percent'/2 high % and 'percent'/2 low %
%   dlabel2     - A title label for the statistics table
%   map         - Data vector to be displayed as topographic map. If a single integer,
%                 only the corresponding electrode location is displayed
%   chan_locs   - name of an EEG electrode position file (See >> topoplot example for format).
%                 Can also be a structure (see >> help pop_editset)
%
% Outputs:
%   M,SD        - mean and standard deviation
%   sk,k        - skewness and excess kurtosis
%   med         - median
%   zlow,zhi    - low and high 'percent/2'-Percentile ('percent/2'/100-Quantile)
%   tM,tSD      - trimmed mean and SD, removing data<zlow and data>zhigh
%   tndx        - index of the data retained after trimming
%   ksh         - output flag of the Kolmogorov-Smirnov test at level p=0.05
%                 0: data could be normally distributed; 1: data are not normally distributed
%                 -1: test could not be executed
%
% Author: Luca Finelli, CNL / Salk Institute - SCCN, 2 August 2002
%
% See also:
%   POP_SIGNALSTAT, QQDIAGRAM, EEGLAB

% Copyright (C) 2002 Luca Finelli, Salk/SCCN, La Jolla, CA

% Note:
% QQDIAGRAM IS EQUIVALENT TO PERCENTILE/PERCENTILE PLOT
% X = EEG.data(5,:); % data
% Y = randn(1, 1000); % gaussan random distribution
% figure; qqdiagram(X, Y,  2);
% figure; plot(prctile(X,2), prctile(Y,2));

% 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.

function [M,SD,sk,k,med,zlow,zhi,tM,tSD,tndx,ksh] = signalstat( data, plotlab, dlabel, percent, dlabel2, map, chan_locs);

M=[]; SD=[]; sk=[]; k=[]; med=[]; zlow=[]; zhi=[]; tM=[]; tSD=[]; tndx=[]; ksh=[];

istats=1;
%hs = help('stats');
% toolbx = ver;
% if isempty(hs) || all([~any(strcmpi({toolbx.Name},'statistics toolbox')), ~any(strcmpi({toolbx.Name},'statistics and machine learning toolbox'))])
%     disp('SIGNALSTAT note: the boxplot (not shown) requires the MATLAB Statistics Toolbox or Statistics and Machine Learning Toolbox');
%     istats=0;
% end

if (nargin<8 && nargin>5) && min(size(map))~=1
		error('signalstat(): the map input must be a vector')
end

if nargin<7 && nargin>5
	disp('signalstat(): no location file for the topographic map')
	help signalstat;
	return
end

if nargin < 6
	map = [];
end

if nargin < 5
	dlabel2 = '';
end

if nargin>3
	if isempty(percent)
		percent=5;
	end
	if any(percent > 100) || any(percent < 0)
		error('signalstat(): percent must be between 0 and 100');
	end
end

if nargin < 4
	percent = 5;
end

if (nargin < 3 || isempty(dlabel))
	dlabel='Potential [V]';
end

if nargin < 2
  plotlab=1;
end;

if ~isnumeric(plotlab)
	error('signalstat(): plotlab must be numeric');
end

if plotlab ~= 0 && plotlab ~= 1
		error('signalstat(): plotlab must be 0 or 1');
end
if nargin < 1
	help signalstat;
	return;
end;

if ndims(data)>2
	error('signalstat(): data must be a vector (1-dim signal)')
end

if ~isreal(data)
	error('signalstat(): data cannot be complex')
end

fprintf('signalstat(): computing statistics...\n');

% Statistical characteristics
%----------------------------
pnts=length(data);   % number of data points
rg=max(data)-min(data);

M=mean(data);        % mean
med=median(data);    % median

vr=var(data);        % variance (N-1 normalized)
SD=std(data);        % standard deviation

if istats
	sk=skewness(data,0); % skewness (third central moment divided by
                         % the cube of the standard deviation)
    k=kurtosis(data,0)-3;  % kurtosis (fourth central  moment divided by
                         % fourth power of the standard deviation)
else
	sk=NaN;
	k=kurt(data)-3;
end

% Checks on skewness and kurtosis
%--------------------------------
sklab='Distribution is symmetric';
if sk>0.01
	sklab='Distribution is right-skewed';
elseif sk < -0.01
	sklab='Distribution is left-skewed';
end

klab='';
if k>0.01
	klab='Distribution is super-Gaussian'; % i.e. kurtosis bigger then Gaussian
elseif k < -0.01
	klab='Distribution is sub-Gaussian';
end

% Estimates without the highest and lowest 'percent'/2 % of data
%---------------------------------------------------------------
pc=percent/100;

zlow = quantile(data,(pc / 2));   % low  quantile
zhi  = quantile(data,1 - pc / 2); % high quantile
tndx = find((data >= zlow & data <= zhi & ~isnan(data)));

tM=mean(data(tndx)); % mean with excluded pc/2*100% of highest and lowest values
tSD=std(data(tndx)); % trimmed SD

% Selected central tendency estimator
%------------------------------------
cte=M;

% Normal fit
%-----------
if istats
	alpha=0.05;          % 1-alpha confidence interval
	[muhat,sigmahat,muci,sigmaci] = normfit(data,alpha);
end

nbins=max(50,round(pnts/100));
[nel,binpos]=hist(data,nbins);
dx=binpos(2)-binpos(1);               % bin width

datafit=normpdf(binpos,cte,SD);       % estimated pdf
datafit=datafit*pnts*dx;

tdatafit=normpdf(binpos,cte,tSD);     % estimated pdf with trimmed SD
tdatafit=tdatafit*pnts*dx;

%datarnd=normrnd(cte,sigmahat,1,pnts); % synthetic data

if istats
	% Goodness-of-fit hypothesis test
	%--------------------------------
	kstail = 0; % 0 = 2-sided test

	CDF=normcdf(data,cte,sigmahat);        % estimated cdf

	[ksh,ksp,ksstat,kscv] = kstest(data,[data', CDF'],alpha,kstail); % Kolmogorov-Smirnov test

	kstestlab='Kolmogorov-Smirnov test: verified Gaussian';
	kscol=[0.2 1 0.2];
	if ksh
		kstestlab='Kolmogorov-Smirnov test: not Gaussian';
		kscol=[.7 .3 .3];
	end
end

% Graphics
%-------------------------------------

if plotlab
  figure
  try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end
  COLOR = [0.56 .66 .9];
  set(gcf,'NumberTitle','off','Name','Signal statistics -- signalstat()')
  fwidth=800;  % figure size in pixels
  fheight=600;

  funits=get(0,'Units');
  set(0,'Units','Pixel')
  scnsize=get(0,'ScreenSize');
  fpos=[round((scnsize(3)-fwidth)/2),round((scnsize(4)-fheight)/2),fwidth,fheight];
  set(0,'Units',funits)

  set(gcf,'Position',fpos)

  % Plotting the histogram
  %------------------------
  subplot(2,2,1)

  hist(data,nbins);  % 'XLim',[-125 125]
  uapos = get(gca,'Position');
  xlim  = get(gca,'XLim');
  set(gca,'FontSize',14)
  xlabel(dlabel)
  title('Data Histogram and Fitted Normal PDF')

  % HM=pnts*normpdf(M,muhat,sigmahat)/2; % FWHM height
  % plot([cte-sd cte+sd],[HM HM],'r--','LineWidth',2)

  % Overplotting a normal distribution
  %-----------------------------------
  hold on
  h1=plot(binpos,datafit,'c','LineWidth',2);
  set(gca,'XLim',xlim)

  % Overplotting a normal distribution from trimmed SD
  %---------------------------------------------------
  h2=plot(binpos,tdatafit,'y');
  ymin=get(gca,'YLim');
  plot([zlow zlow],[0 ymin(2)/20],'y','LineWidth',2) % low  percentile
  plot([zhi  zhi], [0 ymin(2)/20],'y','LineWidth',2) % high percentile
  set(gca,'XLim',xlim)

  % Overplotting a mean and zero line
  %----------------------------------
 % xmin=get(gca,'XLim');
 % ymin=get(gca,'YLim');
  plot([0 0],ymin,'k')
  h3=plot([cte cte],ymin,'r--','LineWidth',2);
  set(gca,'Color',COLOR,'XMinorTick','on','XLim',xlim)

  if istats
	  set(gca,'XTick',[])
  elseif ~istats && strcmp(dlabel,'Potential [V]')
	  set(gca,'XTick',[-125, -75, -25, 0, 25, 75,  125],...
			  'XTickLabel',['-125' ; ' -75' ; ' -25' ; '  0 ' ; ' 25 ' ; ' 75 ' ; ' 125'])
  end

  if strcmp(dlabel,'Potential [V]')
	  set(gca,'XLim',[-125 125])
  end

  set(gca,'FontSize',10)
  H=[h1 h2 h3];
  legend(H,'Gaussian fit','Trimmed G.fit','Mean')
  legend boxoff

  zoom off

  % Boxplot
  %--------------------
  if istats
	  subplot(2,2,3)

	  boxplot(data,1,'+',0,1.5)
	  lapos=get(gca,'Position');
	  set(gca,'Position', [uapos(1) uapos(2)-uapos(4)/3 uapos(3) uapos(4)/3])
	  nlapos=get(gca,'Position');
	  hold on
	  ymin2=get(gca,'YLim');
	  plot([0 0],[0 ymin2(2)],'k')
	  plot([cte cte],[0 ymin(2)],'r--','LineWidth',2)
	  set(gca,'FontSize',14,'XMinorTick','on')
      set(gca,'XLim',xlim)

      if strcmp(dlabel,'Potential [V]')
		  set(gca,'XTick',[-125 -75 -25 0 25 75  125],...
				  'XTickLabel',['-125' ; ' -75' ; ' -25' ; '  0 ' ; ' 25 ' ; ' 75 ' ; ' 125'],...
				  'XLim',[-125 125])
	  end

	  xlabel(dlabel)
	  ylabel('')
	  zoom off
  end

  % QQ plot
  %--------
  subplot(2,2,2)

  qqdiagram(data)
  apos=get(gca,'Position');
  set(gca,'Position', [1-uapos(1)-uapos(3) uapos(2)-uapos(4)/3 (uapos(3)) (uapos(4)+uapos(4)/3)])
  set(gca,'XTick',[-4 -2 0 2 4])
  xmin=get(gca,'XLim');
  ymin=get(gca,'YLim');
  hold on
  plot([xmin(1) xmin(1)+diff(xmin)/20],[zlow zlow],'y-','LineWidth',2)
  plot([xmin(1) xmin(1)+diff(xmin)/20],[zhi  zhi] ,'y-','LineWidth',2)
  set(gca,'XLim',xmin);
  %plot([0 0],ymin,'k--')
  set(gca,'FontSize',14)
  xlabel('Standard Normal Quantiles [Std.Dev.]')
  if strcmp(dlabel,'Potential [V]')
	  ylabel('Ordered Observations [V]')
  elseif strcmp(dlabel,'Component Activity')
	  ylabel('Ordered Observations [rel. V]')
  else
	  ylabel('Ordered Observations')
  end

  title('QQ Plot (Data vs Standard Normal)')
  set(gca,'Color',COLOR)

  % TOPO plot
  %---------
  if (~isempty(map))
	  sbplot(7,9,6)
	  % th=axes('Position',[]);
	  % subplot('Position',[.10 .86 .20 .14]);
	  fprintf('signalstat(): plotting a topographic map...\n');
	  if length(map) == 1
		  topoplot(map,chan_locs,'electrodes','off', ...
				   'style', 'blank', 'emarkersize1chan', 10);
	  else
		  topoplot(map,chan_locs,'electrodes','off');
	  end
	  axis('square')
  end

  % Color schemes
  %--------------
  nero    = [0 0 0];
  bordeau = [0.7 0.3 0.3];
  rosso   = [1 0 0];
  roschi1 = [1 .3 0.4];
  giallo1 = [1 .9 0];
  arancio = [1 .5 .3];
  verchi1 = [0.2 1 0.2];
  verchi2 = [0.7 1 0.7];
  verchi3 = [0.4 1 0.4];
  verscu1 = [0.1 0.7 0.2];
  bluchi1 = [0.4 0.9 1];
  bluchi2 = [.2 .5 .7];
  grichia = [.95 .95 .95];

  % Data axis
  %------------
  bgcolor=COLOR;

  dah = axes('Position',[uapos(1) .1 1-2*uapos(1) .25]);
  set(dah,'Box','on','Color',bgcolor,'XTick',[],'YTick',[],'FontName','Courier','FontSize',12,'FontWeight','demi')
  title(dlabel2,'FontWeight','bold','FontSize',14,'FontName','Arial');
  text(0.05,0.9,['Mean:          ' num2str(M,3)]  ,'Color',bordeau,'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.05,0.8,['Trimmed mean:  ' num2str(tM,3)] ,'Color',giallo1,'FontName','Courier','FontSize',12,'FontWeight','demi')

  text(0.05,0.5,['Standard dev.: ' num2str(SD,4)] ,'Color',bordeau,'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.05,0.4,['Trimmed st.d.: ' num2str(tSD,4)],'Color',giallo1,'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.05,0.3,['Variance:      ' num2str(vr,4)] ,'Color',giallo1,'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.05,0.2,['Range:         ' num2str(rg,4)] ,'Color',giallo1,'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.05,0.1,['Data points:   ' num2str(pnts)] ,'Color',bordeau,'FontName','Courier','FontSize',12,'FontWeight','demi')

  text(0.4,0.9,[num2str(percent/2/100,'%1.3f') '-quantile: ' num2str(zlow,3)] ,'Color',verchi2,...
	   'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.4,0.8,['0.5  -quantile: ',num2str(med,3),'  (median)'],'Color',verchi1,'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.4,0.7,[num2str((100-percent/2)/100,'%1.3f') '-quantile:  ' num2str(zhi,3)] ,'Color',verchi2,...
	   'FontName','Courier','FontSize',12,'FontWeight','demi')

  text(0.4,0.3,['Excess kurtosis: ' num2str(k, 3) ' (near 0 if Gaussian)'] ,'Color',verchi1,...
				'FontName','Courier','FontSize',12,'FontWeight','demi')
  text(0.4,0.2,klab,'Color',verchi2,'FontName','Courier','FontSize',12,'FontWeight','demi')

  if istats
	  text(0.4,0.5,['Skewness: ' num2str(sk,3) ' (near 0 if Gaussian)'] ,'Color',verchi1,...
		   			'FontName','Courier','FontSize',12,'FontWeight','demi')
	  text(0.4,0.4,sklab,'Color',verchi2,'FontName','Courier','FontSize',12,'FontWeight','demi')
      text(0.4,0.1,kstestlab,'Color',kscol,'FontName','Courier','FontSize',12,'FontWeight','demi')
  end
  axcopy;
end

%--------------------------------------------------
% clone of the normpdf function of the stat toolbox
function fitvals = normpdf(myvals,mymean,mystd)
if nargin < 3,
    mystd = 1;
end
if nargin < 2;
    mymean = 0;
end
if length(mymean) < length(myvals)
	tmpmean = mymean;
	mymean = zeros(size(myvals));
	mymean(:) = tmpmean;
end
if length(mystd) < length(myvals)
	tmpmean = mystd;
	mystd = zeros(size(myvals));
	mystd(:) = tmpmean;
end
mymean(1:10);
mystd(1:10);

fitvals = zeros(size(myvals));
tmp = find(mystd > 0);
if any(tmp)
    myvalsn = (myvals(tmp) - mymean(tmp)) ./ mystd(tmp);
    fitvals(tmp) = exp(-0.5 * myvalsn .^2) ./ (sqrt(2*pi) .* mystd(tmp));
end
tmp1 = find(mystd <= 0);
if any(tmp1)
    tmp2   = NaN;
    fitvals(tmp1) = tmp2(ones(size(tmp1)));
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text: %LDIVIDE Mapping overload

author: marianux

repo: ecg-kit

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


text:
function [m] = um2m(um)
% Convert length from micrometers (or microns) to meters.
% Chad A. Greene 2012
m = um*1e-6;


text:
function [Population,FrontNo,CrowdDis] = DominationSelection(Global,Population)
% The dominant relationship and crowding based environmental selection

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Shufen Qin
% E-mail: shufen.qin@stu.tyust.edu.cn

    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(Population.objs,Global.N);
    Next = false(1,length(FrontNo));
    Next(FrontNo<MaxFNo) = true;

    %% Calculate the crowding distance of each solution
    CrowdDis = CrowdingDistance(Population.objs,FrontNo);

    %% Select the solutions in the last front based on their crowding distances
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(CrowdDis(Last),'descend');
    Next(Last(Rank(1:Global.N-sum(Next)))) = true;

    %% Population for next generation
    FrontNo    = FrontNo(Next);
    CrowdDis   = CrowdDis(Next);
    Population = Population(Next);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function desc = sift( patches )
psz = [size(patches, 1), size(patches, 2)];
frm = [(psz(1) ./ 2 + 0.5) * ones(2, 1) ; psz(1) ./ 2; 0];
desc = [];
for pi = 1:size(patches, 3)
  I = single(patches(:, :, pi));
  [Ix, Iy] = vl_grad(I) ;
  mod      = sqrt(Ix.^2 + Iy.^2) ;
  ang      = atan2(Iy, Ix) ;
  grd      = shiftdim(cat(3, mod, ang), 2) ;
  d        = vl_siftdescriptor(grd, frm, 'magnif', 0.5) ;
  if isempty(desc)
    desc = zeros(numel(d), size(patches, 3), 'single');
  end
  desc(:, pi) = d;
end

end

author: hpatches

repo: hpatches-benchmark


text:
function [completed_data, predicted_label, energy] = rec_completion_test(model, test_data, mask, show, param)
% Given incomplete 3D shape(tsdf), recognition and completion test for
% multi-class model.

% Input the trained model, and some fixed masks, this function perform
% completion and recognition simultaneously. Returns the recognition
% accuracy for the incomplete data(test_data) and the completed_data as
% well as the free energy and predicted_label associated with that
% completion.

if ~exist('param','var');
    param.batch_size = 32;
    param.epochs = 50;
    param.gibbs_iter = 1;
    param.earlyStop = false;
end

addpath voxelization;
addpath 3D;
global kConv_backward  kConv_backward_c kConv_forward2 kConv_forward_c;

if ~isfield(model.layers{2},'uw')
    model = merge_model(model);
end

n = size(test_data,1);
num_layer = length(model.layers);
batch_size = 32;
batch_num = ceil(n / batch_size);

predicted_label = zeros(n, model.classes);
completed_data = initialize_missing(test_data, mask);

ncount = 0;
running_mean = zeros(1, model.classes);

for epoch = 1 : param.epochs
    for b = 1 : batch_num
        batch_end = min(n, b * batch_size);
        batch_data = completed_data((b-1) * batch_size+1 : batch_end, :,:,:,:);
        batch_label = predicted_label((b-1) * batch_size+1 : batch_end, :,:,:,:);
        batch_mask = mask((b-1) * batch_size+1 : batch_end,:,:,:);
        this_size = size(batch_data, 1);
        % propagate/inference bottum up using recognition weight.
        for l = 2 : num_layer - 1
            if l == 2
                hidden_presigmoid = myConvolve2(kConv_forward2, batch_data, model.layers{l}.uw, model.layers{l}.stride, 'forward');
                hidden_presigmoid = bsxfun(@plus, hidden_presigmoid, permute(model.layers{l}.c, [2,3,4,5,1]));
                batch_hidden_prob = sigmoid(hidden_presigmoid);
            elseif strcmp(model.layers{l}.type, 'convolution')
                hidden_presigmoid = myConvolve(kConv_forward_c, batch_data, model.layers{l}.uw, model.layers{l}.stride, 'forward');
                hidden_presigmoid = bsxfun(@plus, hidden_presigmoid, permute(model.layers{l}.c, [2,3,4,5,1]));
                batch_hidden_prob = sigmoid(hidden_presigmoid);
            else
                batch_data = reshape(batch_data, this_size, []);
                hidden_presigmoid = bsxfun(@plus, ...
                    batch_data * model.layers{l}.uw, model.layers{l}.c);
                batch_hidden_prob = 1 ./ ( 1 + exp(- hidden_presigmoid) );
            end
            batch_data = batch_hidden_prob;
        end

        batch_data = reshape(batch_data, this_size, []);
        % calculate the free energy for each label hypothesis
        %for c = 1 : model.classes
        %    try_label = zeros(batch_size, model.classes);
        %    try_label(:,c) = 1;
        %    batch_label(:,c) = free_energy(model, [try_label, batch_data], num_layer);
        %end

        hn_1 = batch_data;
        hn_1 = single(hn_1 > rand(size(hn_1)));

        temp_w = model.layers{num_layer}.w;
        temp_w(1:model.classes,:) = temp_w(1:model.classes,:) * model.duplicate;
        for i = 1 : param.gibbs_iter
            % alternating gibbs
            % prop up
            hn = bsxfun(@plus, [batch_label, hn_1] * temp_w, model.layers{num_layer}.c);
            hn = 1 ./ (1 + exp(-hn));
            hn = single(hn > rand(size(hn)));

            % prop down
            hn_1 = bsxfun(@plus, hn * model.layers{num_layer}.w', model.layers{num_layer}.b);
            batch_label = exp(bsxfun(@minus, hn_1(:,1:model.classes), max(hn_1(:,1:model.classes), [], 2)));
            batch_label = bsxfun(@rdivide, batch_label, sum(batch_label, 2));
            hn_1 = 1 ./ ( 1 + exp(-hn_1(:,model.classes+1:end)));
        end

        batch_data = reshape(hn_1, [this_size, model.layers{num_layer-1}.layerSize]);

        for l = num_layer - 1 : -1 : 2
            if l == 2
                batch_data = reshape(batch_data, [this_size, model.layers{l}.layerSize]);
                presigmoid = myConvolve(kConv_backward, batch_data, model.layers{l}.dw, model.layers{l}.stride, 'backward');
                presigmoid = bsxfun(@plus, presigmoid, permute(model.layers{l}.b, [5,1,2,3,4]));
                batch_data = 1 ./ ( 1 + exp(-presigmoid));
            elseif strcmp(model.layers{l}.type, 'convolution')
                batch_data = reshape(batch_data, [this_size, model.layers{l}.layerSize]);
                presigmoid = myConvolve(kConv_backward_c, batch_data, model.layers{l}.dw, model.layers{l}.stride, 'backward');
                presigmoid = bsxfun(@plus, presigmoid, permute(model.layers{l}.b, [5,1,2,3,4]));
                batch_data = 1 ./ ( 1 + exp(-presigmoid));
            else
                batch_data = reshape(batch_data, [this_size, model.layers{l}.layerSize]);
                presigmoid = bsxfun(@plus, ...
					batch_data * model.layers{l}.dw', model.layers{l}.b);
                batch_data = 1 ./ ( 1 + exp(-presigmoid) );
            end
        end
        % clamp the real data
        this_data = test_data((b-1) * batch_size + 1 : batch_end,:,:,:);
        batch_data(~batch_mask) = this_data(~batch_mask);

        completed_data((b-1) * batch_size + 1 : batch_end,:,:,:) = batch_data;
        predicted_label((b-1) * batch_size + 1 : batch_end,:) = batch_label;
    end

    if all(mean(predicted_label,1) == running_mean)
        ncount = ncount + 1;
    else
        running_mean = mean(predicted_label,1);
        ncount = 1;
    end
    % early stop
    if ((~param.earlyStop && epoch >= 100) || (param.earlyStop)) && ncount > 30 || epoch == param.epochs
        break;
    end
end

% calculate the free_energy
energy = zeros(n,1);
for b = 1 : batch_num
    batch_end = min(n, b * batch_size);
    batch_data = completed_data((b-1) * batch_size+1 : batch_end, :,:,:,:);
    batch_label = predicted_label((b-1) * batch_size+1 : batch_end, :,:,:,:);
    this_size = size(batch_data, 1);
    % propagate/inference bottum up using recognition weight.
    for l = 2 : num_layer - 1
        if l == 2
            hidden_presigmoid = myConvolve2(kConv_forward2, batch_data, model.layers{l}.uw, model.layers{l}.stride, 'forward');
            hidden_presigmoid = bsxfun(@plus, hidden_presigmoid, permute(model.layers{l}.c, [2,3,4,5,1]));
            batch_hidden_prob = sigmoid(hidden_presigmoid);
        elseif strcmp(model.layers{l}.type, 'convolution')
            hidden_presigmoid = myConvolve(kConv_forward_c, batch_data, model.layers{l}.uw, model.layers{l}.stride, 'forward');
            hidden_presigmoid = bsxfun(@plus, hidden_presigmoid, permute(model.layers{l}.c, [2,3,4,5,1]));
            batch_hidden_prob = sigmoid(hidden_presigmoid);
        else
            batch_data = reshape(batch_data, this_size, []);
            hidden_presigmoid = bsxfun(@plus, ...
                batch_data * model.layers{l}.uw, model.layers{l}.c);
            batch_hidden_prob = 1 ./ ( 1 + exp(- hidden_presigmoid) );
        end
        batch_data = batch_hidden_prob;
    end
    energy((b-1) * batch_size+1 : batch_end) = free_energy(model, [batch_label, batch_data], num_layer);
end

% show the result
if show
    for i = 1 : n
        the_sample = test_data(i,:,:,:,:);

        figure;
        subplot(1,2,1);
        plot3D(squeeze(completed_data(i,:,:,:,:)) > 0.1);
        p = patch(isosurface(squeeze(completed_data(i,:,:,:)),0.1));
        set(p,'FaceColor','red','EdgeColor','none');
        daspect([1,1,1])
        view(3); axis tight
        camlight
        lighting gouraud

        subplot(1,2,2);
        plot3D(squeeze(the_sample) > 0);
        p = patch(isosurface(squeeze(the_sample),0.1));
        set(p,'FaceColor','red','EdgeColor','none');
        daspect([1,1,1])
        view(3); axis tight
        camlight
        lighting gouraud
        pause;
        close(gcf);
    end
end

function data = initialize_missing(data, mask)
    n = size(data,1);
    bin  = 9;
    div = floor(n / bin); rev = mod(n,bin);
    rev_start = 5 - floor(rev / 2); rev_end = 5 + floor((rev-1) / 2);
    if rev == 0
        for i = 1 : bin
            this_data = data((i-1)*div+1:i*div,:,:,:);
            this_mask = mask((i-1)*div+1:i*div,:,:,:);
            this_data(this_mask) = single(rand(size(this_data(this_mask))) > (0.1 * i ));
            data((i-1)*div+1:i*div,:,:,:) = this_data;
        end
    else
        for i = 1 : rev_start-1
            this_data = data((i-1)*div+1:i*div,:,:,:);
            this_mask = mask((i-1)*div+1:i*div,:,:,:);
            this_data(this_mask) = single(rand(size(this_data(this_mask))) > (0.1 * i ));
            data((i-1)*div+1:i*div,:,:,:) = this_data;
        end
        np = (rev_start - 1) * div;
        for i = rev_start : rev_end
            this_data = data(np + (i-rev_start)*(div+1)+1 : np + (i-rev_start+1) * (div+1),:,:,:);
            this_mask = mask(np + (i-rev_start)*(div+1)+1 : np + (i-rev_start+1) * (div+1),:,:,:);
            this_data(this_mask) = single(rand(size(this_data(this_mask))) > 0.1 * i);
            data(np + (i-rev_start)*(div+1)+1 : np + (i-rev_start+1) * (div+1),:,:,:) = this_data;
        end
        np = (rev_start - 1) * div + rev * (div+1);
        for i = rev_end+1 : bin
            this_data = data(np + (i-rev_end-1)*div+1 : np + (i-rev_end) * div,:,:,:);
            this_mask = mask(np + (i-rev_end-1)*div+1 : np + (i-rev_end) * div,:,:,:);
            this_data(this_mask) = single(rand(size(this_data(this_mask))) > 0.1 * i);
            data(np + (i-rev_end-1)*div+1 : np + (i-rev_end) * div,:,:,:) = this_data;
        end
    end

function [y] = sigmoid(x)
	y = 1 ./ (1 + exp(-x));

author: zhirongw


text:
% Multiple plots using plotPub

clear all;
addpath('../lib');

%% lets plot 3 cycles of 50Hz AC voltage
f = 50;
Vm = 10;
phi = pi/4;

% generate the signal
t = [0:0.0001:3/f];
th = 2*pi*f*t;
v1 = Vm*sin(th);
v2 = Vm*sin(th - phi);
v3 = Vm*sin(th - phi*2);

figure;
plot(t*1E3, v1);
hold on;
plot(t*1E3, v2);
plot(t*1E3, v3);
hold off;

%% settings
opt = [];
opt.XLabel = 'Time, t (ms)'; % xlabel
opt.YLabel = 'Voltage, V (V)'; %ylabel
opt.YTick = [-10, 0, 10];
opt.YLim = [-11, 11];

% Save? comment the following line if you do not want to save
opt.FileName = 'plotMultiple.tiff';

% create the plot
setPlotProp(opt);

author: masumhabib

repo: PlotPub

save_path: github-repos/MATLAB/masumhabib-PlotPub