text:
function status=sdephaseplot3(t,y,flag,w)
%SDEPHASEPLOT3  3-D phase space SDE output function.
%   When the function SDEPHASEPLOT3 is passed to an SDE solver as the
%   'OutputFUN' property, i.e., OPTIONS = SDESET('OutputFUN',@SDEPHASEPLOT3),
%   the solver calls SDEPHASEPLOT3(T,Y,'') after every timestep. The
%   SDEPHASEPLOT3 function plots the first three components of the solution, Y,
%   as they are computed, adapting the axis limits of the plot dynamically. To
%   plot particular components of the solution, specify their indices via the
%   'OutputYSelect' SDE options property. The first specified component is
%   plotted with respect to the X-axis, the second component with respect to the
%   Y-axis, and the third component with respect to the Z-axis.
%
%   At the start of integration, the solver calls SDEPHASEPLOT3(TSPAN,Y0,'init')
%   to initialize the output function. After each integration step to new time,
%   T, and solution vector, Y, the solver calls STATUS = SDEPHASEPLOT3(T,Y,'').
%   If the solver's 'Refine' property is greater than one (see SDESET), T is a
%   column vector of new output times and Y is an array of corresponding column
%   vectors. The STATUS return value is 1 if the figure window and plot axis are
%   still open and 0 otherwise. When the integration is complete, the solver
%   calls SDEPHASEPLOT3([],[],'done').
%
%   If fewer than three components of the solution, Y, are specified, one or two
%   components of the integrated Wiener increments, W, can be plotted versus Y
%   solution components or three integrated Wiener increment components can be
%   plotted. Set the 'OutputWSelect' SDE options property to 'yes' or to a
%   vector of indices to pass the integrated Wiener increments to SDEPHASEPLOT3
%   as a fourth argument, SDEPHASEPLOT3(T,Y,'',W). If the 'OutputWSelect'
%   property is set to 'yes' or a non-empty vector of indices the solver calls
%   SDEPHASEPLOT3(TSPAN,Y0,'init',W0) to initialize the output function at the
%   start of integration and SDEPHASEPLOT3([],[],'done',[]) when the integration
%   is complete.
%
%   See also:
%       SDEPLOT, SDEPHASEPLOT2, SDEIMGPLOT, SDEPRINT, SDESET, SDEGET, SDE_EULER,
%       SDE_MILSTEIN, SDE_BM, SDE_GBM, SDE_OU, ODEPHAS3

%   SDEPHASEPLOT3 is based on an updating of Matlab's ODEPHAS3,
%   version 1.27.4.10

%   Andrew D. Horchler, horchler @ gmail . com, 5-11-13
%   Revision: 1.2, 5-13-13

persistent FIG_HANDLE AX_HANDLE;
status = 1;         % Figure widow still open and and plot axis still present
isW = (nargin > 3); % Have integrated Wiener increments been passed
if nargin < 3
    flag = '';
end

switch flag
    case 'init'
        % Initialize persitent handle variables
        FIG_HANDLE = figure(gcf);
        AX_HANDLE = gca;

        % Set units to pixels to get width of axis, set back to default
        units = get(AX_HANDLE,'Units');
        set(AX_HANDLE,'Units','Pixels');
        pos = get(AX_HANDLE,'OuterPosition');
        set(AX_HANDLE,'Units',units);

        % Number of time samples to expect
        LEN_TSPAN = length(t);

        % Use figure axis width and TSPAN length determine redraw chunk
        chunk = min(ceil(LEN_TSPAN/pos(3)),LEN_TSPAN);

        % Number of solution variables, Y (cannot change)
        N = length(y);

        % Initialize UserData and Y
        ud = [];
        ud.y(3,chunk) = 0;
        ud.i = 1;

        if N >= 3
            ud.y(:,1) = y(1:3);
        elseif isW
            % Number of integrated Wiener increment variables, W (cannot change)
            D = length(w);

            if D+N >= 3
                if N == 2
                    ud.y(:,1) = [y(1:2);w(1)];
                elseif N == 1
                    ud.y(:,1) = [y(1);w(1:2)];
                else
                    ud.y(:,1) = w(1:3);
                end
            else
                error('SDETools:sdephaseplot3:TooFewInputsW',...
            	     ['Output function requires at least three solution '...
                      'components from Y or W.']);
            end
        else
            error('SDETools:sdephaseplot3:TooFewInputs',...
            	 ['Output function requires at least three solution '...
                  'components from Y.']);
        end

        % Plot initial data
        ud.lines = plot3(ud.y(1),ud.y(2),ud.y(3),'-');

        % Plot start and end markers, turn on grid
        nextplot = get(AX_HANDLE,'NextPlot');
        set(AX_HANDLE,'NextPlot','add');
        plot3(ud.y(1),ud.y(2),ud.y(3),'g.');
        ud.marker = plot3(ud.y(1),ud.y(2),ud.y(3),'r.');
        set(AX_HANDLE,'NextPlot',nextplot,'DrawMode','fast');
        grid on;

        % Store UserData and draw
        set(FIG_HANDLE,'UserData',ud);
        drawnow;
    case ''
        if isempty(FIG_HANDLE) || isempty(AX_HANDLE)
            if isW
                error('SDETools:sdephaseplot3:NotCalledWithInitW',...
                     ['Output function has not been initialized. Use syntax '...
                      'OutputFUN(TSPAN,Y0,''init'',W0).']);
            else
                error('SDETools:sdephaseplot3:NotCalledWithInit',...
                     ['Output function has not been initialized. Use syntax '...
                      'OutputFUN(TSPAN,Y0,''init'').']);

            end
        end

        % If figure is open
        if ishghandle(FIG_HANDLE) && ishghandle(AX_HANDLE)
            % Get UserData
            ud = get(FIG_HANDLE,'UserData');
            lt = length(t);
            N = size(y,1);
            lY = size(ud.y,2);

            % Update UserData
            oldi = ud.i;
            newi = oldi+lt;
            if newi > lY
                % Set line data
                XYZData = get(ud.lines,{'XData','YData','ZData'});
                set(ud.lines,{'XData','YData','ZData'},...
                             {[XYZData{1} ud.y(1,:)],...
                              [XYZData{2} ud.y(2,:)],...
                              [XYZData{3} ud.y(3,:)]});

                % Set marker point
                set(ud.marker,{'XData','YData','ZData'},...
                              {ud.y(1,end),ud.y(2,end),ud.y(3,end)});

                % Reset UserData
                ud.y(:,lY) = 0;
                oldi = 0;
                newi = lt;
            end

            % Append new data to UserData
            if N >= 3
                ud.y(:,oldi+1:newi) = y(1:3,:);
            elseif N == 2
                ud.y(:,oldi+1:newi) = [y(1:2,:);w(1,:)];
            elseif N == 1
                ud.y(:,oldi+1:newi) = [y(1,:);w(1:2,:)];
            else
                ud.y(:,oldi+1:newi) = w(1:3,:);
            end
            ud.i = newi;

            % Store updated UserData and draw if redraw chunk was full
            set(FIG_HANDLE,'UserData',ud);
            if oldi == 0
                drawnow;
            end
        else
            status = 0;
        end
    case 'done'
        % Rename and delete persistent handles
        hf = FIG_HANDLE;
        FIG_HANDLE = [];
        ha = AX_HANDLE;
        AX_HANDLE = [];

        % If figure is open
        if ~isempty(hf) && ishghandle(hf) && ~isempty(ha) && ishghandle(ha)
            % Get non-zero UserData
            ud = get(hf,'UserData');
            ud.y = ud.y(:,1:ud.i);

            % Set any remaining line data
            XYZData = get(ud.lines,{'XData','YData','ZData'});
            set(ud.lines,{'XData','YData','ZData'},...
                         {[XYZData{1} ud.y(1,:)],...
                          [XYZData{2} ud.y(2,:)],...
                          [XYZData{3} ud.y(3,:)]});

            % Set final marker point
            set(ud.marker,{'XData','YData','ZData'},...
                          {ud.y(1,end),ud.y(2,end),ud.y(3,end)});

            % Delete UserData
            set(hf,'UserData',[]);

            % Refresh or draw
            if ishold(ha)
                drawnow;
            else
                refresh;
            end
        else
            status = 0;
        end
    otherwise
        error('SDETools:sdephaseplot3:InvalidFlag',...
              'Invalid status flag passed to output function.');
end

author: horchler

repo: SDETools

save_path: github-repos/MATLAB/horchler-SDETools


text:
function alpha = formatAngle(alpha)
%FORMATANGLE  Ensure an angle value is comprised between 0 and 2*PI
%   ALPHA2 = formatAngle(ALPHA)
%   ALPHA2 is the same as ALPHA modulo 2*PI and is positive.
%
%   Example:
%   formatAngle(5*pi)
%   ans =
%       3.1416
%
%   See also
%   vectorAngle, lineAngle
%
% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2008-03-10,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% HISTORY
% 2010-03-31 deprecate and replace by function 'normalizeAngle'

% deprecation warning
warning('geom2d:deprecated', ...
    '''formatAngle'' is deprecated, use ''normalizeAngle'' instead');

alpha = mod(alpha, 2*pi);
alpha(alpha<0) = 2*pi + alpha(alpha<0);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function job = clusterVariantGraph(job,varargin)
% the MCL algorithm generalized to hyper graphs
%
% Syntax
%
%   job.calcVariantGraph
%   job.clusterVariantGraph
%   job.calcParentFromVote('minProb',0.5)
%
% Input
%  job - @parentGrainReconstructor
%
% Output
%  job.votes - computed parent votes
%
% Options
%  inflationPower - inflation power of the MCL algorithm (default=1.05)
%  numIter        - number of iterations (default=10)
%  cutOff         - default 0.0001
%  withDiagonal   -
%  keepGraph      - do not kill graph after clustering
%  includeSimilar - similarly oriented variants share probability
%  includeTwins   - potential twins share probability
%  minCluster     - minimum cluster size (default 4)
%
% References
%
% * <https://arxiv.org/abs/2201.02103 The variant graph approach to
% improved parent grain reconstruction>, arXiv, 2022
%

% ensure we have a variant graph
if isempty(job.graph), job.calcVariantGraph(varargin{:}); end

% extract parameters
p = get_option(varargin,'inflationPower', 1.05);
numIter = get_option(varargin,'numIter', 10);
cutOff = get_option(varargin,'cutOff',0.0001);

% some constants
A = job.graph;
nVG = size(A,1); % size of the variant graph
job.graph = [];
p2cV = variants(job.p2c,'parent');
numV = length(p2cV);

% indices to the variant graph
isCP = job.isChild | job.isParent; % only child and parent phases go into the variant graph

% rep2VG - how often it repeats to fit in the variant graph [numV numV 1 numV 0 1 numV ...]
rep2VG = numV * job.isChild + job.isParent;

% VG2ind - variant graph index -> grain index
VG2ind = 1:length(job.grains);
VG2ind = repelem(VG2ind(isCP), rep2VG(isCP));

% indeces of the pseudo diagonal
if check_option(varargin,'noPseudoDiagonal')

  % position of variant 1 of each child grain in the Matrix
  indV1 = [0;cumsum(rep2VG)];
  indV1(end) = [];
  indV1 = indV1(job.isChild);

  % a numV x numV matrix without diagonal
  [i,j] = meshgrid(1:numV,1:numV);
  i(eye(numV)==1) = [];
  j(eye(numV)==1) = [];

  pseudoDiag = sparse(indV1 + i(:).',indV1 + j(:).',true,size(A,1),size(A,2));
end

normalize

% this is an adaptation of the MCL algorithm
for iter = 1:numIter

  % expansion
  A = A * A;

  % remove pseudo diagonal entries, i.e., A(j+k,j+l)
  if exist('pseudoDiag','var'), A(pseudoDiag) = 0; end

  % inflation by Hadamard potentiation
  A = A.^p;

  % prune elements of A that are below minval
  %A = (A > cutOff) .* A;
  A = spfun(@(x) (x > cutOff) .* x,A);

  if check_option(varargin,'sym')
    A = sqrt(A.' .* A);
    normalize
  end

  % column re-normalisation
  normalize

  if check_option(varargin,'test2')
    A = sqrt(diag(diag(A)) * A);
    normalize
  end

  if check_option(varargin,'verbose'), disp(nnz(A)); end
end

if check_option(varargin,'test1')
  A = diag(diag(A)) * A;
  normalize
end

if check_option(varargin,'sym')
  A = sqrt(A.' .* A);
  normalize
end

% create a table of probabilities for the different parentIds of each child
% grains
if p>1
  s = full(sum(A,2)); % sum columns
else
  s = full(sum(A,2) - diag(A)); % sum columns
end

% ensure minimum cluster size
minCluster = get_option(varargin,'minCluster',4);
s = s .* (full(sum(A>0,2))>=minCluster);

% store in numGrains x numVariant matrix
isVGChild = job.isChild(VG2ind); % isChild in variant graph
pIdP = nan(length(job.grains),numV);
pIdP(job.isChild,:) = reshape(s(isVGChild),numV,[]).';

% some post processing
if check_option(varargin,'includeTwins')

  % define theoretical twinning of the packet
  tp2cV = p2cV(:) .* orientation.byAxisAngle(round(p2cV \ Miller(0,1,1,p2cV.SS),'maxHKL',1),60*degree);

  M = 0.5*triu(angle_outer(p2cV,tp2cV,'noSym2') < 5 * degree);

  M(~sum(M),~sum(M)) = eye(nnz(~sum(M)));
  pIdP = pIdP * M;
end

if check_option(varargin,'includeSimilar')
  %M = triu(angle_outer(job.p2c.variants('parent'),job.p2c.variants('parent'),'noSym2')<5*degree);
  %M(:,sum(M)==1) = 0;

  M = angle_outer(p2cV,p2cV,'noSym2')<5*degree;
  M = 0.9*M + 0.1*eye(size(M));

  pIdP = pIdP * M;

elseif check_option(varargin,'mergeSimilar')

  M = angle_outer(p2cV,p2cV,'noSym2')<5*degree;
  pIdP = pIdP * M;

end

% sort the table and store the highest probabilities in the job class
[pIdP,pId] = sort(pIdP,2,'descend');
job.votes = table(pId,pIdP,'VariableNames',{'parentId','prob'});

if check_option(varargin,'keepGraph')
  job.graph = A;
end

  function normalize2

    % column re-normalisation
    % sum over all targets
    s = full(sum(A));

    % sum over all variants
    s = accumarray(VG2ind.',s);
    s = repelem(s,rep2VG(isPC));

    % create sparse diagonal matrix for
    dinv = spdiags(1./s,0,nVG,nVG);
    A = A * dinv;

  end
  function normalize

    % column re-normalisation
    % sum over all targets
    s = full(sum(A,2)).';

    % sum over all variants
    s = accumarray(VG2ind.',s);
    s = repelem(s(isCP),rep2VG(isCP));

    % create sparse diagonal matrix for
    dinv = spdiags(1./s,0,nVG,nVG);
    A = dinv * A;

  end
end

author: mtex-toolbox

repo: mtex

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


text:
classdef SVOREX < Algorithm
    %SVOREX Support Vector for Ordinal Regression (Explicit constraints)
    %   This class derives from the Algorithm Class and implements the
    %   SVOREX method. This class uses SVOREX implementation by
    %   W. Chu et al (http://www.gatsby.ucl.ac.uk/~chuwei/svor.htm)
    %
    %   SVOREX methods:
    %      fitpredict               - runs the corresponding algorithm,
    %                                   fitting the model and testing it in a dataset.
    %      fit                        - Fits a model from training data
    %      predict                    - Performs label prediction
    %
    %   References:
    %     [1] P.A. Gutiérrez, M. Pérez-Ortiz, J. Sánchez-Monedero,
    %         F. Fernández-Navarro and C. Hervás-Martínez
    %         Ordinal regression methods: survey and experimental study
    %         IEEE Transactions on Knowledge and Data Engineering, Vol. 28. Issue 1
    %         2016
    %         http://dx.doi.org/10.1109/TKDE.2015.2457911
    %     [2] W. Chu and S. S. Keerthi, Support Vector Ordinal Regression,
    %         Neural Computation, vol. 19, no. 3, pp. 792–815, 2007.
    %         http://10.1162/neco.2007.19.3.792
    %
    %   This file is part of ORCA: https://github.com/ayrna/orca
    %   Original authors: Pedro Antonio Gutiérrez, María Pérez Ortiz, Javier Sánchez Monedero
    %   Citation: If you use this code, please cite the associated paper http://www.uco.es/grupos/ayrna/orreview
    %   Copyright:
    %       This software is released under the The GNU General Public License v3.0 licence
    %       available at http://www.gnu.org/licenses/gpl-3.0.html

    properties
        description = 'Support Vector for Ordinal Regression (Explicit constraints)';
        parameters = struct('C', 0.1, 'k', 0.1);
    end
    properties (Access = private)
        algorithmMexPath = fullfile(fileparts(which('Algorithm.m')),'SVOREX');
    end

    methods
        function obj = SVOREX(varargin)
            %SVOREX constructs an object of the class SVOREX and sets its default
            %   characteristics
            %   OBJ = SVOREX(KERNEL) builds SVOREX with RBF as kernel function
            obj.parseArgs(varargin);
        end

        function [projectedTrain, predictedTrain] = privfit(obj,train,parameters)
            %PRIVFIT trains the model for the SVOREX method with TRAIN data and
            %vector of parameters PARAMETERS.
            if isempty(strfind(path,obj.algorithmMexPath))
                addpath(obj.algorithmMexPath);
            end
            [alpha, thresholds, projectedTrain] = svorex([train.patterns train.targets],parameters.k,parameters.C,0,0,0);
            predictedTrain = obj.assignLabels(projectedTrain, thresholds);
            model.projection = alpha;
            model.thresholds = thresholds;
            model.parameters = parameters;
            model.train = train.patterns;
            obj.model = model;
            projectedTrain = projectedTrain';
            if ~isempty(strfind(path,obj.algorithmMexPath))
                rmpath(obj.algorithmMexPath);
            end
        end

        function [projected, predicted] = privpredict(obj, test)
            %PREDICT predicts labels of TEST patterns labels. The object needs to be fitted to the data first.
            kernelMatrix = computeKernelMatrix(obj.model.train',test','rbf',obj.model.parameters.k);
            projected = obj.model.projection*kernelMatrix;

            predicted = SVOREX.assignLabels(projected, obj.model.thresholds);
            projected = projected';
        end
    end

    methods (Static = true)
        function predicted = assignLabels(projected, thresholds)
            numClasses = size(thresholds,2)+1;
            %TEST assign the labels from projections and thresholds
            project2 = repmat(projected, numClasses-1,1);
            project2 = project2 - thresholds'*ones(1,size(project2,2));

            % Asignation of the class
            % f(x) = max {Wx-bk<0} or Wx - b_(K-1) > 0
            wx=project2;

            % The procedure for that is the following:
            % We assign the values > 0 to NaN
            wx(wx(:,:)>0)=NaN;

            % Then, we choose the biggest one.
            [maximum,predicted]=max(wx,[],1);

            % If a max is equal to NaN is because Wx-bk for all k is >0, so this
            % pattern belongs to the last class.
            predicted(isnan(maximum(:,:)))=numClasses;

            predicted = predicted';
        end
    end
end

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
classdef EventWithJointLimits < Simulink.IntEnumType
    %EVENTWITHJOINTLIMITS Class for selecting a None, Warning or EmergencyStop signal
    %
    %   Class for selecting a None, Warning or EmergencyStop signal triggered by the joints hitting
    %   the limits.

    enumeration
        None(0)
        Warning(1)
        EmergencyStop(2)
    end
end

author: robotology

repo: whole-body-controllers


text:
for i=1:2000
    input(i,:)=10*rand(1,2)-5;
    output(i)=input(i,1)^2+input(i,2)^2;
end
output=output';

save data input output

author: HuangCongQing


text:
function sim = gaussianKernel(x1, x2, sigma)
%RBFKERNEL returns a radial basis function kernel between x1 and x2
%   sim = gaussianKernel(x1, x2) returns a gaussian kernel between x1 and x2
%   and returns the value in sim

% Ensure that x1 and x2 are column vectors
x1 = x1(:); x2 = x2(:);

% You need to return the following variables correctly.
sim = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the similarity between x1
%               and x2 computed using a Gaussian kernel with bandwidth
%               sigma
%
%

sim = exp(-sum((x1 - x2).^2)/(2 * sigma^2));

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

end

author: rmarquis

repo: coursera-machinelearning


text:
function [Hdraw,HvarsDraw,phi_Hdraw,Adraw,h0]=sampleH_Ben(yData,Psi,A,B,phi_H,HvarsOld,priorValues,dataValues,n,h0,HH)

%yData=YData
%Psi=PsiDraw_prop
%A=Adraw
%B=Bdraw
%phi_H=phi_Hdraw'
%HvarsOld=HvarsDraw
%% Initialize
offset_c=priorValues.offset_c; %constant for log transformation

[T,M]=size(yData);
p = size(B,1)/M; %lags of the B vector that has dimensions MxM*p (no constant), each column is one regression

% obtain prior data
b0=priorValues.mean_ln_h0;  %log mean of the initial state (variance scaling parameters) as residuals from an AR(4) in the training sample
a0=eye(n)/priorValues.var_ln_h0(1,1);    %variance of the initial state for the elements of lambda

priorPhi_H=priorValues.phi_h;            %centering parameter for the inverse gamma of the variance of the innovations governing the random walk for lambda
priorD_H  =priorValues.d_h;              %scaling parameter for the inverse gamma distribution for the variance of the innovations governing the random walk for lambda

S_h = priorPhi_H*ones(n,1);
nu_h = 3*ones(n,1);
%% Prepare Data

Y_Psi=yData-Psi;  %subtract the local mean
Y_Psi(1:p,:)=yData(1:p,:)-ones(p,1)*mean(yData(1:p,:)); %also generate initial conditions for the construction of the lagmatrix
X_Psi = lagmatrix(Y_Psi,1:p); %create RHS of the VAR part
X_Psi = X_Psi(p+1:end,:);     %remove the first p rows of RHS
Y_Psi=Y_Psi(p+1:end,:);       %and do so for LHS

E=Y_Psi-X_Psi*B;              %VAR residuals
Escaled=E*A';                 %transform such that VAR residuals have variance Lambda

Hvars=zeros(T,M);             %sampled states (diagonal elements of lambda)

%% Sample (States,Vars,phi)

for i=1:M
    % prepare for each i
    residsTemp=Escaled(:,i);
    yStar=log(residsTemp.^2+offset_c); %transform scaled residuals
    phi=phi_H(i); %previous draw for the variance of this particular variable
    lnSigma2=log(HvarsOld(p+1:T,i)); % note log => h = ln sigma2

    % sample offset mixture
    [yStarAdj, Ht] = bear.statesMix(yStar,lnSigma2);

   % sample h conditional on the state
    iSig_s = sparse(1:T-p,1:T-p,1./Ht); %state dependent variance
    Kh = HH/phi + iSig_s;
    h_hat = Kh\(h0(i,1)/phi*HH*ones(T-p,1) + iSig_s*(yStarAdj));
    logVarsDraw_Hi = h_hat + chol(Kh,'lower')'\randn(T-p,1);
    Hvars(:,i)=exp([zeros(1,p) logVarsDraw_Hi']');

    % sample phi
%     [phiDraw_Hi]=samplePhi(logVarsDraw_Hi,priorD_H,priorPhi_H);
%     phi_Hdraw(i)=phiDraw_Hi;

end

   % sample phi
    diff = (log(Hvars(p+1:end,:)) - [h0'; log(Hvars(p+1:end-1,:))]).^2;
    phi_Hdraw = 1./gamrnd(nu_h + size(diff,1)/2, 1./(S_h + sum(diff)'/2))';

% sample the initial value
    Kh0 = a0 + sparse(1:n,1:n,1./phi_Hdraw');
    h0_hat = Kh0\(a0*b0 + log(Hvars(p+1,:)')./phi_Hdraw');
    h0 = h0_hat + chol(Kh0,'lower')'\randn(n,1);

%% Sample A
[Adraw]=bear.sampleA_H(yData,Psi,B,Hvars,T,priorValues,dataValues);
AdrawInv=Adraw\eye(M);

%% Construct H
Hdraw=zeros(M,M,T);

for t=1:T
    Hdraw(:,:,t)=AdrawInv*diag(Hvars(t,:))*AdrawInv';
end

HvarsDraw=Hvars;
end

author: european-central-bank

repo: BEAR-toolbox


text:
function printAPResults( classes, results )
if ~isnumeric(results)
    aps = [results(:).ap]';
else
    aps = results;
end

for i = 1:numel(classes)
    class_string = classes{i}(1:min(5,length(classes{i})));
    fprintf('& %5s ', class_string)
end
fprintf('& %5s \\\\ \n', 'mean')
for i = 1:numel(classes)
    fprintf('& %2.3f ', aps(i))
end
fprintf('& %2.3f \\\\ \n', mean(aps))

end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function PLOT_Map(States,Params, Trajectory)

%This function plots an empty version of the map.

x = States(1,:);
y = States(2,:);
th = States(3,:);
phi = States(4,:);

Lt = Params.Dyn.Lt;
Lc = Params.Dyn.Lc;

Npts = Params.Traj.Npts;

%[dStates, A, B, Positions] = Derive_EoM()

Pa = [...   %Position of the back of the trailer
    x;
    y];

Pb = [...    %Position of the back of the cab
 x - Lt*sin(th);
 y + Lt*cos(th)];

Pc = [...   %Position of the front of the cab
 x - Lc*(cos(phi).*sin(th) + cos(th).*sin(phi)) - Lt*sin(th);
 y + Lc*(cos(phi).*cos(th) - sin(phi).*sin(th)) + Lt*cos(th)];

if strcmp(Params.Sim.Direction,'reverse')
    Dir = 2;
else
    Dir = 1;
end
Road = zeros(2,Npts);
for i=1:Npts
    Road(:,i) = Trajectory.States{Dir,i}(1:2);
end

figure(1)
clf; hold on;
plot(Pa(1,1),Pa(2,1),'r','LineWidth',5)
plot(Pb(1,1),Pb(2,1),'b','LineWidth',3)
plot(Pc(1,1),Pc(2,1),'g','LineWidth',3)
legend('Back of Trailer','Back of Cab', 'Front of Cab')
plot(Road(1,:),Road(2,:),'k','LineWidth',Params.Sim.RoadWidth)
title('Map of the Road and Traces')
axis equal

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function Population = SubcomponentOptimizer(Problem,Population,Neighbour,indices)
% Subcomponent optimizer

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

    for i = 1 : length(Population)
        if rand < 0.9
            P = Neighbour(i,randperm(size(Neighbour,2),2));
        else
            P = randperm(length(Population),2);
        end
        OffDec          = Population(i).dec;
        NewDec          = OperatorDE(Problem,OffDec,Population(P(1)).dec,Population(P(2)).dec,{1,0.5,length(OffDec)/length(indices)/2,20});
        OffDec(indices) = NewDec(indices);
        Offspring       = Problem.Evaluation(OffDec);
        if sum(Offspring.obj) < sum(Population(i).obj)
            Population(i) = Offspring;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function model = updateFieldOrderForType(model,type,newOrder)
% Reorder the fields associated with the provided type to the new order.
% USAGE:
%     matchingFields = updateFieldOrderForType(model, type, newOrder)
%
% INPUTS:
%
%    model:              the model to update
%    type:               the Type of field to update one of
%                        ('rxns','mets','comps','genes')
%    newOrder:           The new Order. must be of the same size as the
%                        requested model.(type) field.
%
% OUTPUT:
%
%    model:              The model with all field associated with type
%                        reordered
%
% .. Authors:
%                   - Thomas Pfau June 2017

parser = inputParser();
parser.addRequired('model',@isstruct);
parser.addRequired('type',@(x) any(ismember(x,{'rxns','mets','comps','genes','proteins'})));
parser.addRequired('newOrder',@(x) isnumeric(x) && all(ismember(newOrder,1:numel(model.(type)))));

parser.parse(model,type,newOrder);

[fields,dimensions] = getModelFieldsForType(model,type);

for i = 1:numel(fields)
    %For now, we will only handle 2 dimensional arrays
    if dimensions(i) == 1
        model.(fields{i})(:,:) = model.(fields{i})(newOrder,:);
    elseif dimensions(i) == 2
        model.(fields{i})(:,:) = model.(fields{i})(:,newOrder);
    end
end

if strcmp(type,'genes')
    %the gene positions changed, so we have to rebuild the rules field.

    if isfield(model,'rules')
        for i = 1:numel(model.genes)
            if i ~= newOrder(i)
                %replace by new with an indicator that this is new.
                model.rules = strrep(model.rules,['x(' num2str(newOrder(i)) ')'],['x(' num2str(i) '$)']);
            end
        end
        %remove the indicator.
        model.rules = strrep(model.rules,'$','');
    end
end

author: opencobra

repo: cobratoolbox


text:
function options = insert_param_when_absent(options, fieldname, default_value)
if ~isfield(options, fieldname)
    options.(fieldname) = default_value;
end

end

author: zhouyuanzxcv

repo: Hyperspectral


text:
classdef GemLabelButton < GemVirtualPanel
    % GemLabelButton GEM class for a button image with label text below
    %
    %
    %
    %     Licence
    %     -------
    %     Part of GEM. https://github.com/tomdoel/gem
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = protected)
        Button
        Text
        Tag
    end

    properties
        DefaultButtonWidth = 50
        DefaultButtonHeight = 50
        TextHeight = 25
        VerticalSpacing = 2
        ButtonHorizontalSpacing = 0
        LabelFontSize = 9
    end

    properties (Dependent)
        ButtonWidth
        ButtonHeight
    end

    events
        ButtonClicked
    end

    methods
        function obj = GemLabelButton(parent, text, tooltip, tag, icon)
            obj = obj@GemVirtualPanel(parent);
            obj.Button = GemButton(parent, text, tooltip, tag, @obj.ButtonClickedCallback);
            obj.Button.ButtonWidth = obj.DefaultButtonWidth;
            obj.Button.ButtonHeight = obj.DefaultButtonHeight;
            obj.Button.UnSelectedColour = [50, 50, 50];
            obj.Button.ShowTextOnButton = false;
            obj.Button.AddAndResizeImageWithHighlightMask(icon, [0, 0, 0]);
            obj.Button.AutoUpdateStatus = true;

            obj.Button.UnSelectedColour = uint8(255*(obj.StyleSheet.BackgroundColour));
            obj.Button.SelectedColour = uint8(255*(obj.StyleSheet.IconSelectedColour));
            obj.Button.HighlightColour = uint8(255*(obj.StyleSheet.IconHighlightColour));
            obj.Button.HighlightSelectedColour = uint8(255*(obj.StyleSheet.IconHighlightSelectedColour));

            obj.AddChild(obj.Button);
            obj.Text = GemText(parent, text, tooltip, tag);
            obj.Text.HorizontalAlignment = 'center';
            obj.Text.FontSize = obj.LabelFontSize;
            obj.Tag = tag;
            obj.AddChild(obj.Text);
        end

        function Resize(obj, new_position)
            Resize@GemVirtualPanel(obj, new_position);

            button_x_pos = new_position(1) + obj.ButtonHorizontalSpacing;
            button_width = obj.Button.ButtonWidth;
            text_width = button_width + 2*obj.ButtonHorizontalSpacing;
            button_y_pos = new_position(2) + obj.TextHeight + obj.VerticalSpacing;
            button_height = obj.Button.ButtonHeight;
            obj.Button.Resize([button_x_pos, button_y_pos, button_width, button_height]);
            obj.Text.Resize([new_position(1), new_position(2), text_width, obj.TextHeight]);
        end

        function set.ButtonWidth(obj, width)
            obj.Button.ButtonWidth = width;
        end

        function width = get.ButtonWidth(obj)
            width = obj.Button.ButtonWidth;
        end

        function set.ButtonHeight(obj, height)
            obj.Button.ButtonHeight = height;
        end

        function height = get.ButtonHeight(obj)
            height = obj.Button.ButtonHeight;
        end

        function width = GetWidth(obj)
            width = obj.Button.ButtonWidth + 2*obj.ButtonHorizontalSpacing;
        end

        function height = GetRequestedHeight(obj, width)
            height = obj.Button.ButtonHeight + obj.TextHeight + obj.VerticalSpacing;
        end

        function Enable(obj)
            Enable@GemVirtualPanel(obj);
            obj.Button.Enable;
            obj.Text.Enable;
        end

        function Disable(obj)
            Disable@GemVirtualPanel(obj);
            obj.Button.Disable;
            obj.Text.Disable;
        end

        function Select(obj, is_selected)
            obj.Button.Select(is_selected);
        end
    end

    methods (Access = protected)
        function ButtonClickedCallback(obj, tag)
            notify(obj, 'ButtonClicked', CoreEventData(tag));
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function C = estimateC(W,R,B,C,lam,tol)

% INPUT:
% W: 2-by-P matrix
% R: 2-by-3 matrix
% B: 3K-by-P matrx
% C: 1-by-K vector

% This function estimates c by minimizing
% f(C) + lam*r(C), where
% f(C) = 0.5 * norm(W-R*S,'fro')^2 where S = \sum_i C_i*B_i
% r(C) = \|C\|_1
% It implements proximal gradient + nesterov

P = size(W,2);
K = size(B,1)/3;

C = C'; % transpose to be a column vector
C0 = C; % C0: the previous estimate
t = 1; % auxiliary variable for nesterov
t0 = 1; % auxiliary variable for nesterov
fvalue = inf;

% next we work on the linear system y = X*C
y = W(:); % vectorized W
X = zeros(2*P,K); % each column is a rotated Bk
for k = 1:K
    RBk = R*B(3*k-2:3*k,:);
    X(:,k) = RBk(:);
end

if lam == 0
    C = X\y;
    C = C';
    return
end

% mu is set as the 2-norm of the Hessian of f(C)
mu = norm(X'*X);

for iter = 1:1000

    % Z is an auxiliary variable in nesterov method
    Z = C + (t0-1)/t*(C-C0);

    % gradient descent
    Z = Z + X'*(y-X*Z)/mu;

    % nonegative thresholding
%     Z = Z - lam/mu;
%     Z = max(Z,0);

    % soft thresholding
    Z = sign(Z).*max(abs(Z)-lam/mu,0);

    % update C
    C0 = C;
    C = Z;

    % update t
    t0 = t;
    t = (1+sqrt(1+4*t^2))/2;

    % function value
    fvalue0 = fvalue;
    fvalue = 0.5 * norm(y-X*C,'fro')^2 + lam*sum(abs(C));
    if fvalue > fvalue0
        t0 = 1; % APG with restart
        t = 1;
    end

    % check convergence
    RelChg = norm(C-C0,'fro')/(norm(C0,'fro')+eps) ;
%     fprintf('Iter %d: FunVal = %f, RelChg = %f\n',iter,fvalue,RelChg);
    if RelChg < tol
        break
    end

end

C = C'; % transpose back to be a row vector

author: geopavlakos


text:
%This function develop by Renoald Tang
%University Teknologi Malaysia, Photogrammetry and Laser scanning group
%for academic purpose
%Email:renoald@live.com
%This function read Off file
function readOFF(file)
[np,nf]=value(file);
fprintf('Number_of_point:%d\n',np);
fprintf('Number_of_face:%d\n',nf);
ext(file,np);

%return number of face and point
function [np,nf]=value(file)
C=0;
fid=fopen(file);
while 1
    tline = fgetl(fid);
    C=C+1;
    if C==2

      num=textscan(tline,'%d %d %d');
        break;
    end
    if ~ischar(tline),   break,   end

end
fclose(fid);
np=num{1,1};
nf=num{1,2};
function ext(file,n)
C=0;
K=0;
fid=fopen(file);
fid1=fopen('point','w');
fid2=fopen('face','w');
while 1
    tline=fgetl(fid);
    C=C+1;
     if ~ischar(tline) , break , end
    if C > 2
       K=K+1;
       if K <=n
           fprintf(fid1,'%s\n',tline);
       else
           fprintf(fid2,'%s\n',tline);
       end

    end

end
fclose(fid);
fclose(fid1);
fclose(fid2);%return and extract face and point value
%disp(K)

author: Sable

repo: mcbench-benchmarks


text:
classdef GradientVariationExperiment < handle

    properties (Access = public)
        backgroundMesh
        boundaryMesh
        levelSet
        regularizedPerimeter
        domainLength
    end

    properties (Access = private)
        iMesh
        inputFile
        levelSetParams
    end

    methods (Access = public)

        function obj = GradientVariationExperiment(cParams)
            obj.init(cParams)
            obj.createBackgroundAndBoundaryMesh();
            obj.computeDomainLength();
            obj.createLevelSet();
            obj.createRegularizedPerimeters();
        end

    end

    methods (Access = protected)

        function init(obj,cParams)
            obj.inputFile      = cParams.inputFile;
            obj.iMesh          = cParams.iMesh;
            obj.levelSetParams = cParams.levelSetParams;
        end

        function createBackgroundAndBoundaryMesh(obj)
            s.inputFile = obj.inputFile;
            s.isBackgroundMeshRectangularBox = true;
            mCreator = BackgroundAndBoundaryMeshCreatorFromInputFile(s);
            obj.backgroundMesh = mCreator.backgroundMesh;
            obj.boundaryMesh = mCreator.boundaryMesh;
        end

        function createLevelSet(obj)
            s = obj.levelSetParams;
            s.coord      = obj.backgroundMesh.coord;
            s.ndim       = obj.backgroundMesh.ndim;
            lsCreator = LevelSetCreator.create(s);
            obj.levelSet = lsCreator.getValue();
        end

        function computeDomainLength(obj)
            x = obj.backgroundMesh.coord;
            d = max(x(:,1)) - min(x(:,1));
            obj.domainLength = d;
        end

        function createRegularizedPerimeters(obj)
            s.inputFile        = obj.inputFile;
            s.backgroundMesh   = obj.backgroundMesh;
            s.scale            = 'MACRO';
            s.designVariable   = obj.levelSet;
            s.outputFigureName = ['SmoothedCircleMesh',num2str(obj.iMesh)];
            s.plotting         = false;
            s.printing         = false;
            s.capturingImage   = false;
            s.isRobinTermAdded = true;
            s.perimeterType    = 'perimeterInterior';
            rPerimeter = RegularizedPerimeterComputer(s);
            rPerimeter.compute();
            obj.regularizedPerimeter = rPerimeter;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%% Example 4.5: Solution of the Ornstein–Uhlenbeck process
%
% Copyright:
%   2018 - Simo Särkkä and Arno Solin
%
% License:
%   This software is provided under the MIT License. See the accompanying
%   LICENSE file for details.

%% Simulate trajectories from the OU process

  if exist('rng') % Octave doesn't have rng
      rng(10,'twister')
  else
      randn('state',1)
  end

  lambda = 0.5;
  q = 1;
  dt = 0.01;
  T = (0:dt:1);
  x0 = 4;
  M = exp(-lambda*T)*x0;
  P = q/(2*lambda)*(1 - exp(-2*lambda*T));

  XX = zeros(50,length(T));
  for n=1:size(XX,1)
    x = x0;
    for k=1:length(T)
      XX(n,k) = x;
      x = x - lambda * x * dt + sqrt(dt)*randn;
    end
  end

  figure(1); clf; hold on

    % Shade the 95% quantiles
    fill([T fliplr(T)],[M+1.96*sqrt(P) fliplr(M-1.96*sqrt(P))],1, ...
      'FaceColor',[.9 .9 .9],'EdgeColor',[.9 .9 .9])

    % Plot realizations
    h1 = plot(T,XX,'-','Color',[.5 .5 .5],'LineWidth',0.5);

    % Plot mean and quantiles
    h2 = plot(T,M,'k-','LineWidth',1);
    h34 = plot(T,M+1.96*sqrt(P),'--k', ...
               T,M-1.96*sqrt(P),'--k','LineWidth',0.7);

    %set(h(1:3),'Linewidth',2)

    legend([h2(1) h34(1) h1(1)],'Mean','95\% quantiles','Realizations')
    xlabel('Time, $t$'); ylabel('$x(t)$')

    ylim([0 5])
    box on
    set(gca,'Layer','Top')

author: AaltoML

repo: SDE

save_path: github-repos/MATLAB/AaltoML-SDE


text:
function [out] = melt_3(p1,p2,T,S1,S2,St,dt,varargin)
%melt_3

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Glacier melt provided no snow is stored on the ice layer
% Constraints:  f <= S1/dt
% @(Inputs):    p1   - degree-day factor [mm/oC/d]
%               p2   - temperature threshold for snowmelt [oC]
%               T    - current temperature [oC]
%               S1   - current storage in glacier [mm]
%               S2   - current storage in snowpack [mm]
%               St   - storage in S2 threshold below which glacier melt occurs [mm]
%               dt   - time step size [d]
%               varargin(1) - smoothing variable r (default 0.01)
%               varargin(2) - smoothing variable e (default 5.00)

if size(varargin,2) == 0
    out = min(max(p1*(T-p2),0),S1/dt).*smoothThreshold_storage_logistic(S2,St);
elseif size(varargin,2) == 1
    out = min(max(p1*(T-p2),0),S1/dt).*smoothThreshold_storage_logistic(S2,St,varargin(1));
elseif size(varargin,2) == 2
    out = min(max(p1*(T-p2),0),S1/dt).*smoothThreshold_storage_logistic(S2,St,varargin(1),varargin(2));
end

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function [data] = tapas_linear_prepare_data(data, model, inference)
%%
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

data = data;

end

author: translationalneuromodeling

repo: tapas


text:
function Y = identityPhi(X)
    Y = X;

repo: REST


text:
function [ishappy, cutoff] = looseCheck(f, values, data, pref) %#ok<INUSD,STOUT>
%LOOSECHECK   Loose happiness check for CHEBTECH construction.
%   [ISHAPPY, CUTOFF] = LOOSECHECK(F, VALUES, DATA, PREF). This function
%   is not yet implemented.

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

error('CHEBFUN:CHEBTECH:looseCheck:notDone', ...
    'This file is not implemented yet.')

% TODO: Implement this file.

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% This code demonstrates how the boundary matching code works,
% and how the benchmark uses it.  This is not meant to show
% how to run the benchmark.  See the README file for that.

% setup
present = 'color';
iid = 101085;
nthresh = 10;

% read the image
im = rgb2gray(double(imread(imgFilename(iid)))/255);
figure(1); clf;
imshow(im);

% create a pb image
%pb = pbNitzberg(im);
pb = pbGM(im);
figure(2); clf;
imagesc(pb,[0 1]);
axis image; axis off; truesize;

% read segmentations
segs = readSegs(present,iid);

% match the first seg and a thresholded pb
bmap1 = double(seg2bmap(segs{1}));
bmap2 = double(pb > 0.5);
[match1,match2,cost,oc] = correspondPixels(bmap1,bmap2,0.01,1000);
h=figure(3); clf;
plotMatch(h,bmap1,bmap2,match1,match2);
title('Seg1 vs. Pb>0.5','Color',[1 1 1]);

% compare the pb and segs
[thresh,cntR,sumR,cntP,sumP] = boundaryPR(pb,segs,nthresh);

% precision/recall plot
r = cntR./(sumR+(sumR==0));
p = cntP./(sumP+(sumP==0));
f = 2.*r.*p./(r+p+((r+p)==0));
figure(4); clf;
plot(r,p,'-o');
axis equal; axis([0 1 0 1]);
xlabel('Recall'); ylabel('Precision');

% find best F-measure (should interpolate)
[t,idx] = max(f(:));
title(sprintf('F=%.2g at (R,P)=(%.2g,%.2g) t=%.2g',...
              f(idx),r(idx),p(idx),thresh(idx)));

author: Cloud-CV

repo: object-proposals


text:
%REMCLASS Remove small classes
%
%   B = REMCLASS(A,N)
%
% INPUT
%   A   Dataset
%   N   Integer, maximum class size to be removed (optional; default 0)
%
% OUTPUT
%   B   Dataset
%
% DESCRIPTION
% Classes having N objects or less are removed. The corresponding objects
% are made unlabeled. Use SELDAT to remove unlabeled objects.
%
% SEE ALSO
% SELDAT, GENDAT

function [b, L] = my_remclass(a,n)
	prtrace(mfilename);
	if nargin < 2, n = 0; end

	N = classsizes(a);
	J = find(N <= n);
	L = findnlab(a,J);
	b = setnlab(a,0,L);
	b = setlablist(b);

return

author: marianux

repo: ecg-kit

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


text:
function [trainSets,testSets] = getSplits(outcome,nSplit,testSplit)

indNeg = find(outcome == 0); nNeg = numel(indNeg);
indPos = find(outcome == 1); nPos = numel(indPos);
nPosTest = round(testSplit*nPos); % Number of positive instances in the test sets
nNegTest = round(testSplit*nNeg); % Number of negative instances in the test sets
nInst = numel(outcome);
nTest = nPosTest + nNegTest;
nTrain = nInst - nTest;

trainSets = zeros(nTrain,nSplit);
testSets = zeros(nTest,nSplit);
for s = 1:nSplit
    indPosTest = indPos(datasample(1:nPos,nPosTest,'Replace',false)');
    indNegTest = indNeg(datasample(1:nNeg,nNegTest,'Replace',false)');
    indTest = [indPosTest;indNegTest]; indPerm = randperm(nTest)'; indTest = indTest(indPerm);
    indTrain = (1:nInst)'; indTrain(indTest) = []; indPerm = randperm(nTrain)'; indTrain = indTrain(indPerm);
    trainSets(:,s) = indTrain; testSets(:,s) = indTest;
end

end

author: mvallieres

repo: radiomics


text:
%
% Usage:   B =mexInvSym(A);
%
% Name: mexInvSym
%
% Description: returns the inverse of a symmetric matrix A
%
% Inputs: A:  double n x n matrix
%
% Output: B: double n x n matrix
%
% Author: Julien Mairal, 2009

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function [ problem ] = barbara_problem( )
%sBARBARATEST Prepares the test s of learning dictionary from
% Barbara image
s.image = spx.data.standard_images.barbara_gray_512x512;
% Block size
s.blkSize = 8;
% Number of atoms in dictionary
s.D = 121;
% Signal dimension
s.N = s.blkSize * s.blkSize;
% Let us create all the patches from image
s.patches = im2col(s.image, [s.blkSize, s.blkSize], 'sliding');
% Select 10 % of patches as training patches
s.training_patches = s.patches(:, 1:10:end);
% Let us construct a 1D-DCT matrix of size 8x11.
dr = sqrt(s.D);
DCT=zeros(s.blkSize,dr);
for k=0:1:dr-1,
    is = 0:1:s.blkSize-1;
    % Create one column
    col=cos(is'*k*pi/dr);
    if k>0
        col=col-mean(col);
    end;
    DCT(:,k+1)=col/norm(col);
end;
% We create our initial dictionary as a 64x121 size
s.initial_dictionary = kron(DCT, DCT);
s.DCT = DCT;
% We return the setup
problem = s;
end

author: indigits

repo: sparse-plex


text:
Iref=(im2double(imread('test.jpg')));
tau=[0.2,0.3,0.1,0.2,200,100];
% tau=[0,0,0,0,50,100];
It=Iref;
for i=1:3
It(:,:,i)=warpImg(Iref(:,:,i),tau);
end
It=1.4*abs(It+0.0001)+0.1;
imshow([Iref,It])

tic
[ImTrans,tau,ImTrans_c] = align_c(It,Iref, zeros(6,1),2);
toc;
imshow([It,Iref,ImTrans,ImTrans_c,10*abs(ImTrans_c-Iref)])

author: csjcai

repo: RealSR

save_path: github-repos/MATLAB/csjcai-RealSR


text:
classdef RWMOP34 < PROBLEM
% <multi> <real> <constrained>
% Synchronous pptimal pulse-width modulation of 11-level inverters

%------------------------------- Reference --------------------------------
% A. Kumar, G. Wu, M. Ali, Q. Luo, R. Mallipeddi, P. Suganthan, and S. Das,
% A benchmark-suite of real-world constrained multi-objective optimization
% problems and some baseline results, Swarm and Evolutionary Computation,
% 2021, 67: 100961.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Initialization
        function Setting(obj)
            obj.M        = 2;
            obj.D        = 30;
            obj.lower    = zeros(1,30);
            obj.upper    = 90*ones(1,30);
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            x = varargin{1};
            m = 0.3333;
            s = [1,-1,1,1,1,-1,-1,-1,1,1,1,1,-1,-1,1,-1,-1,-1,1,1,1,1,-1,1,1,-1,-1,1,-1,-1];
            k = [5,7,11,13,17,19,23,25,29,31,35,37,41,43,47,49,53,55,59,61,65,67,71,73,77,79,83,85,91,95,97];
            % Objective function
            for i = 1 : size(x,1)
                su = 0;
                for j = 1 : 31
                    su2 = 0;
                    for l = 1 : size(x,2)
                        su2 = su2 + s(l).*cos(k(j).*x(i,l)*pi/180);
                    end
                    su = su + su2.^2./k(j).^4;
                end
                f(i,1) = (su).^0.5./(sum(1./k.^4)).^0.5;
            end
            f(:,2) = (sum(s.*cos(x*pi/180),2)-m).^2;
            % Constraints
            for i = 1 : size(x,2)-1
                g(:,i) = x(:,i)-x(:,i+1)+1e-6;
            end
            Population = SOLUTION(varargin{1},f,g,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            R = [6.1063821e-01   1.1108847e-01];
        end
    end
end

author: BIMK

repo: PlatEMO


text:
function msns_l_i
%
%  Generates the data used in 'msns_l'. Data are stored in global
%  variables.
%
%  Calling sequence:
%
%    msns_l_i
%
%  Remarks:
%
%    Call 'msns_m_i' before calling this routine.
%
%
%  LYAPACK 1.0 (Thilo Penzl, October 1999)

global LP_N LP_NU

if ~length(LP_N)
  error('This routine needs global data which must be generated by calling ''msns_m_i'' first.');
end

[LP_NU,t] = chol(-LP_N);        % Note the minus!

if t~=0
  error('N is not (numerically) negative definite!');
end


text:
% pop_rejcont() - reject continuous portions of data based on spectrum
%                 thresholding. First, contiguous data epochs are extracted
%                 and a standard spectrum thresholding algorithm is
%                 applied. Regions of contiguous epochs larger than a
%                 specified size are then labeled as artifactual.
%
% Usage:
%   >> pop_rejcont( INEEG ) % pop-up interative window mode
%   >> [OUTEEG, selectedregions] = pop_rejcont( INEEG, 'key', 'val');
%
% Inputs:
%   INEEG      - input dataset
%
% Optional inputs:
%  'elecrange'   - [integer array] electrode indices {Default: all electrodes}
%  'epochlength' - [float] epoch length in seconds {Default: 0.5 s}
%  'overlap'     - [float] epoch overlap in seconds {Default: 0.25 s}
%  'freqlimit'   - [min max] frequency range too consider for thresholding
%                  Default is [35 128] Hz.
%  'threshold'   - [float] frequency upper threshold in dB {Default: 10}
%  'contiguous'  - [integer] number of contiguous epochs necessary to
%                  label a region as artifactual {Default: 4 }
%  'addlength'   - [float] once a region of contiguous epochs has been labeled
%                  as artifact, additional trailing neighboring regions on
%                  each side may also be added {Default: 0.25 s}
%  'eegplot'     - ['on'|'off'] plot rejected portions of data in a eegplot
%                  window. Default is 'off'.
%  'onlyreturnselection'  - ['on'|'off'] this option when set to 'on' only
%                  return the selected regions and does not remove them
%                  from the datasets. This allow to perform quick
%                  optimization of the rejected portions of data.
%  'precompstruct' - [struct] structure containing precomputed spectrum (see
%                  Outputs) to be used instead of computing the spectrum.
%  'verbose'       - ['on'|'off'] display information. Default is 'off'.
%  'taper'         - ['none'|'hamming'] taper to use before FFT. Default is
%                    'none' for backward compatibility but 'hamming' is
%                    recommended.
%
% Outputs:
%   OUTEEG          - output dataset with updated joint probability array
%   selectedregions - frames indices of rejected electrodes. Array of n x 2
%                     n being the number of regions and 2 for the beginning
%                     and end of each region.
%   precompstruct   - structure containing precomputed data. This structure
%                     contains the spectrum, the frequencies and the EEGLAB
%                     dataset used as input with epochs extracted.
%
% Author: Arnaud Delorme, CERCO, UPS/CNRS, 2009-
%
% Example:
% EEG = pop_rejcont(EEG, 'elecrange',[1:32] ,'freqlimit',[20 40] ,'threshold',...
%    10,'epochlength',0.5,'contiguous',4,'addlength',0.25, 'taper', 'hamming');
%
% See also: eegthresh()

% Copyright (C) 2009 Arnaud Delorme, CERCO, UPS/CNRS
%
% 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 [EEG selectedregions precompstruct com ] = pop_rejcont(EEG, varargin);

com = '';
if nargin < 1
    help pop_rejcont;
    return;
end;

if nargin < 2
    firstelec = 'EXG1';                % first non EEG channel
    % take all scalp electrodes
    % -------------------------
    if ~isempty(EEG.chanlocs)
        tmpchanlocs = EEG.chanlocs;
        indelec = strmatch( firstelec, { tmpchanlocs.labels });

        if isempty(indelec), elecrange = 1:EEG.nbchan;
        else                 elecrange = 1:(indelec-1);
        end;
    else
        elecrange = 1:EEG.nbchan;
    end;
    elecrange = deblank(vararg2str(elecrange));
    %elecrange = elecrange(2:end-1);

%     promptstr = { 'Channel range' ...
%                   'Frequency range (Hz)' ...
%                   'Frequency threshold in dB' ...
%                   'Epoch segment length (s)' ...
%                   'Minimum number of contiguous epochs' ...
%                   'Add trails before and after regions (s)' ...
%                   };
%     initstr = { elecrange '20 40' '10' '0.5' '4' '0.25' };
%     result = inputdlg2(promptstr, 'Reject portions of continuous data - pop_rejcont', 1, initstr);
    uilist = { { 'style' 'text' 'string' 'Channel range' } ...
               { 'style' 'edit' 'string' elecrange } ...
               { 'style' 'text' 'string' 'Frequency range (Hz)' } ...
               { 'style' 'edit' 'string' '20 40' } ...
               { 'style' 'text' 'string' 'Frequency threshold in dB' } ...
               { 'style' 'edit' 'string' '10' } ...
               { 'style' 'text' 'string' 'Epoch segment length (s)' } ...
               { 'style' 'edit' 'string' '0.5' } ...
               { 'style' 'text' 'string' 'Minimum number of contiguous epochs' } ...
               { 'style' 'edit' 'string' '4' } ...
               { 'style' 'text' 'string' 'Add trails before and after regions (s)' } ...
               { 'style' 'edit' 'string' '0.25' } ...
               { 'style' 'text' 'string' 'Use hanning window before computing FFT' } ...
               { 'style' 'checkbox' 'string' '' 'value' 1 } ...
                  };
    geom = { [2 1] [2 1] [2 1] [2 1] [2 1] [2 1] [2 1] };
    result = inputgui('uilist', uilist, 'geometry', geom);
    if length( result ) == 0 return; end;

    options = { 'elecrange'     str2num(result{1}) ...
                'freqlimit'     str2num(result{2}) ...
                'threshold'     str2double(result{3}) ...
                'epochlength'   str2double(result{4}) ...
                'contiguous'    str2double(result{5}) ...
                'addlength'     str2double(result{6}) ...
                'taper'         fastif(result{7}, 'hamming', 'none') };
else
    options = varargin;
end;

opt = finputcheck(options, { 'threshold'     { 'real';'cell' }  []    10;
                             'freqlimit'     { 'real';'cell' }  []    [35 128];
                             'elecrange'     'real'   []    [1:EEG.nbchan];
                             'rejectori'     'real'   []    [];
                             'contiguous'    'real'   []    4;
                             'addlength'     'real'   []    0.25;
                             'precompstruct' 'struct' []    struct([]);
                             'eegplot'       'string' { 'on';'off' } 'off';
                             'onlyreturnselection' 'string' { 'on';'off' } 'off';
                             'verbose'       'string' { 'on';'off' } 'on';
                             'taper'         'string' { 'none' 'hamming' } 'none';
                             'overlap'       'real'   []    0.25;
                             'epochlength'   'real'   []    0.5 }, 'pop_rejcont');
if isstr(opt), error(opt); end;
if ~iscell(opt.threshold) && length(opt.threshold) == 2 && ...
    iscell(opt.freqlimit) && length(opt.freqlimit) == 2
    opt.threshold = { opt.threshold(1) opt.threshold(2) };
end;
if ~iscell(opt.threshold), opt.threshold = { opt.threshold }; end;
if ~iscell(opt.freqlimit), opt.freqlimit = { opt.freqlimit }; end;

%EEG.event = [];
grouplen  = opt.contiguous/2*opt.epochlength*EEG.srate+1; % maximum number of points for grouping regions
color     = [ 0 0.9 0];             % color of rejection window

NEWEEG = EEG;
if isempty(opt.precompstruct)
    % compute power spectrum
    % ----------------------
    % average reference
    % NEWEEG.data(opt.elecrange,:) = NEWEEG.data(opt.elecrange,:)-repmat(mean(NEWEEG.data(opt.elecrange,:),1), [length(opt.elecrange) 1]);

    % only keep boundary events
    % -------------------------
    tmpevent = NEWEEG.event;
    if ~isempty(tmpevent)
        if isnumeric( tmpevent(1).type )
            NEWEEG.event = [];
        else
            boundEvent = strmatch('boundary', { tmpevent.type }, 'exact');
            NEWEEG.event = NEWEEG.event(boundEvent);
        end;
    end;

    [TMPNEWEEG] = eeg_regepochs(NEWEEG, opt.overlap, [0 opt.epochlength], NaN);
    %[TMPNEWEEG indices] = pop_rejspec(TMPNEWEEG, 1, [1:64], -100, 15, 30, 45, 0, 0);
    %rejepoch = find(indices);

    tmpdata = TMPNEWEEG.data;
    if strcmpi(opt.taper, 'hamming'),
        tmpdata = bsxfun(@times, tmpdata, hamming(size(TMPNEWEEG.data,2))');
    end;
    tmp   = fft(tmpdata, [], 2);
    freqs = linspace(0, TMPNEWEEG.srate/2, size(tmp,2)/2);
    freqspectrum = freqs(2:end); % remove DC (match the output of PSD)
    tmp   = tmp(:,2:size(tmp,2)/2,:);
    warning('off', 'MATLAB:log:logOfZero');
    tmpspec  = 10*log10(abs(tmp).^2);
    warning('on', 'MATLAB:log:logOfZero');
    tmpspec  = tmpspec - repmat( mean(tmpspec,3), [1 1 TMPNEWEEG.trials]);
    specdata = tmpspec;

    % compute mean spectrum
    % ---------------------
    meanspectrum = nan_mean(specdata(opt.elecrange, :, :), 1);
    precompstruct.spec  = meanspectrum;
    precompstruct.freqs = freqspectrum;
    precompstruct.EEG   = TMPNEWEEG;
else
    meanspectrum  = opt.precompstruct.spec;
    freqspectrum  = opt.precompstruct.freqs;
    TMPNEWEEG     = opt.precompstruct.EEG;
    precompstruct = opt.precompstruct;
end;

% apply threshold to average of all electrodes
% --------------------------------------------
rejepoch = [];
for iReject = 1:length(opt.threshold)
    threshold = opt.threshold{iReject};
    freqLim   = opt.freqlimit{iReject};
    if length(threshold) == 1, threshold = [ -100 threshold ]; end;
    [I1 tmpRejEpoch NS Erej] = eegthresh( meanspectrum, size(meanspectrum,2), 1, threshold(1), threshold(2), [freqspectrum(1) freqspectrum(end)], freqLim(1), freqLim(2));
    rejepoch = union_bc(rejepoch, tmpRejEpoch);
    if strcmpi(opt.verbose, 'on')
        fprintf('%d regions selected for rejection, threshold %3.2f-%3.2f dB, frequency limits %3.1f-%3.1f\n', length(tmpRejEpoch), threshold(1), threshold(2), freqLim(1), freqLim(2));
    end;
end;

% build the winrej array for eegplot
% ----------------------------------
winrej   = [];
if ~isempty(find(cellfun(@isempty, { TMPNEWEEG.event.epoch }) == 1))
    error('Some events are not associated with any epoch');
end;
tmpevent = TMPNEWEEG.event;
allepoch = [ tmpevent.epoch ];
if ~isempty(rejepoch)
    for index = 1:length(rejepoch)
        eventepoch = find( rejepoch(index) == allepoch );
        if strcmpi(TMPNEWEEG.event(eventepoch(1)).type, 'X')
            urevent = TMPNEWEEG.event(eventepoch(1)).urevent;
            lat     = TMPNEWEEG.urevent(urevent).latency;
            winrej  = [ winrej; lat lat+opt.epochlength*TMPNEWEEG.srate-1 color ]; %Erej(:,index)'];
        else
            error('Wrong type for epoch');
        end;
    end;
    winrej(:,6:6+length(opt.elecrange)-1) = 0;
end;

% remove isolated regions and merge others
% ----------------------------------------
merged = 0;
isolated = 0;
for index = size(winrej,1):-1:1
    if size(winrej,1) >= index && winrej(index,2) - winrej(index,1) > grouplen, winrej(index,:) = []; isolated = isolated + 1;
    elseif index == 1 && size(winrej,1) > 1 && winrej(index+1,1) - winrej(index,2) > grouplen, winrej(index,:) = []; isolated = isolated + 1;
    elseif index == size(winrej,1) && size(winrej,1) > 1 && winrej(index,1) - winrej(index-1,2) > grouplen, winrej(index,:) = []; isolated = isolated + 1;
    elseif index > 1 && size(winrej,1) > 1 && index < size(winrej,1) && winrej(index+1,1) - winrej(index,2) > grouplen && ...
            winrej(index,1) - winrej(index-1,2) > grouplen
        winrej(index,:) = [];
        isolated = isolated + 1;
    elseif index < size(winrej,1) && size(winrej,1) > 1 && winrej(index+1,1) - winrej(index,2) <= grouplen
        winrej(index,2) = winrej(index+1,2);
        winrej(index+1,:) = [];
        merged = merged + 1;
    end;
end;
if strcmpi(opt.verbose, 'on')
    fprintf('%d regions merged\n', merged);
    fprintf('%d regions removed\n', isolated);
end;

% add time before and after each region
% -------------------------------------
for index = 1:size(winrej,1)
    winrej(index,1) = max(1,         winrej(index,1)-opt.addlength*EEG.srate);
    winrej(index,2) = min(EEG.pnts,  winrej(index,2)+opt.addlength*EEG.srate);
end;

% plot result
% -----------
if ~isempty(winrej)
    selectedregions = winrej(:,1:2);
    if strcmpi(opt.onlyreturnselection, 'off')
        % merge with initial regions
        if ~isempty(opt.rejectori)
            winrej(:,3) = 1; % color
            for iRow = 1:size(opt.rejectori,1)
                winrej(end+1,1:2) = opt.rejectori(iRow,:);
                winrej(end  ,4)   = 1; % color
                winrej(end  ,5)   = 1; % color
            end;
        end;

        command = '[EEG LASTCOM] = pop_select(EEG, ''nopoint'', TMPREJ(:,1:2)); eegh(LASTCOM); [ALLEEG EEG CURRENTSET LASTCOM] = pop_newset(ALLEEG, EEG, CURRENTSET, ''study'', ~isempty(STUDY)+0); eegh(LASTCOM); eeglab redraw';
        if nargin < 2 || strcmpi(opt.eegplot, 'on')
            eegplot(NEWEEG.data(opt.elecrange,:), 'srate', NEWEEG.srate, 'winrej', winrej, 'command', command, 'events', EEG.event, 'winlength', 50);
            disp('Green is overlap');
            disp('Light blue is ORIGINAL rejection');
            disp('Yellow is AUTOMATIC rejection');
        else
            NEWEEG = pop_select(EEG, 'nopoint', round(selectedregions));
        end;
        EEG = NEWEEG;
    else
        EEG = [];
    end;
else
    selectedregions = [];
    if strcmpi(opt.verbose, 'on')
        disp('No region removed');
    end;
end;
if nargout > 3
    com = sprintf('EEG = pop_rejcont(EEG, %s);', vararg2str(options));
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function pass = test_sin(pref)
% This tests the sine function in Chebfun3.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e2 * pref.cheb3Prefs.chebfun3eps;

% This just constructs a function involving a sine:
f = chebfun3(@(x,y,z) sin(x+y+z));
pass(1) = abs(f(.1,.2,.3) - sin(.6)) < tol;

% This tests the Chebfun3 sine function:
f2 = sin(f);
pass(2) = abs(f2(.2,.2,.2) - sin(sin(.6))) < tol;

% Test computing the sine from a Chebfun3 object in different ways:
f2 = chebfun3(@(x,y,z) sin(f(x,y,z)));
pass(3) = abs(f2(.2,.2,.2) - sin(sin(.6))) < tol;

f2 = chebfun3(@(x,y,z) sin(f(x,y,z)), 'fiberDim', 1);
pass(4) = abs(f2(.2,.2,.2) - sin(sin(.6))) < tol;

f2 = chebfun3(@(x,y,z) sin(f(x,y,z)), 'fiberDim', 2);
pass(5) = abs(f2(.2,.2,.2) - sin(sin(.6))) < tol;

f2 = chebfun3(@(x,y,z) sin(f(x,y,z)), 'fiberDim', 3);
pass(6) = abs(f2(.2,.2,.2) - sin(sin(.6))) < tol;

% This varies the domain:
d = [1 3 1 3 1 3];
x = chebfun3(@(x,y,z) x, d);
y = chebfun3(@(x,y,z) y, d);
z = chebfun3(@(x,y,z) z, d);
f3 = sin(2*x+3*y+z);
pass(7) = abs(f3(1, 2, 3) - sin(11)) < tol;

% Here we check something in trig mode
f4 = chebfun3(@(x,y,z) sin(2*x + 3*y), [-pi pi -pi pi -pi pi], 'trig');
pass(8) = abs(sin(f4(1,1,1)) - sin(sin(5))) < tol;

ep = 1e-8;
tol2 = 1e2 * ep;
f5 = chebfun3(@(x,y,z) exp(x.*y.*z) , 'eps', ep);
pass(9) = abs(sin(f5(.5,.5,.5)) - sin(exp(0.125))) < tol2;

end


text:
function f = perform_mesh_smoothing(face,vertex,f,options)

% perform_mesh_smoothing - smooth a function defined on a mesh by averaging
%
%   f = perform_mesh_smoothing(face,vertex,f,options);
%
%   Smooth a function f on a width of options.niter_averaging vertices.
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;
naver = getoptions(options, 'niter_averaging', 1);
type = getoptions(options, 'averaging_type', 'combinatorial');

if nargin<3
    f = [];
end
if isempty(f)
    f = vertex;
end
if size(f,1)<size(f,2)
    f = f';
end
[vertex,face] = check_face_vertex(vertex,face);

if size(f,2)>1
    for i=1:size(f,2)
        f(:,i) = perform_mesh_smoothing(face,vertex,f(:,i),options);
    end
    return;
end

n = max(face(:));

% compute normalized averaging matrix
if strcmp(type, 'combinatorial')
    %add diagonal
    W = triangulation2adjacency(face) + speye(n);
    D = spdiags(full(sum(W,2).^(-1)),0,n,n);
    W = D*W;
else
    options.normalize=1;
    W = compute_mesh_weight(vertex,face,type,options);
end

% do averaging to smooth the field
for k=1:naver
    f = W*f;
end

author: gpeyre

repo: matlab-toolboxes


text:
function varargout = flipdim(varargin)
%FLIPDIM   Flip/reverse a SPHEREFUN in a chosen direction.
%   G = FLIPDIM(F, DIM) returns a SPHEREFUN G with the same domain as F but
%   reversed in a direction, i.e., G(x,y)=F(x, c+d-y). If DIM = 2 (default) then
%   G(x,y) = F(x, c+d-y).  Otherwise DIM = 1 and G(x,y) = F(a+b-x, y). The
%   domain of F is [a, b, c, d].
%
% See also SPHEREFUN/FLIPLR, SPHEREFUN/FLIPUD.

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

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

end


text:
function R = GenerateRefPoints(Q,diff,alpha,N)
% Generate the reference points according to the combined population

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

    Q    = Q(NDSort(Q.objs,1)==1);
    R    = [];
    subN = min(length(Q),ceil(alpha*N));
    CrowdDis = CrowdingDistanceInEachObj(Q.objs);
    [~,rank] = sort(CrowdDis,1,'descend');
    for m = 1 : length(Q(1).obj)
        Rm      = Q(rank(1:subN,m)).objs;
        Rm(:,m) = Rm(:,m) - diff(m);
        R       = [R;Rm];
    end
    R = R(NDSort(R,1)==1,:);
    if size(R,1) > N
        CrowdDis = CrowdingDistanceInEachObj(R);
        [~,rank] = sort(sum(CrowdDis,2),'descend');
        R        = R(rank(1:N),:);
    end
end

function CrowdDis = CrowdingDistanceInEachObj(PopObj)
% Calculate the crowding distance of each solution in each objective

    [N,M]    = size(PopObj);
    CrowdDis = zeros(N,M);
    Fmax     = max(PopObj,[],1);
    Fmin     = min(PopObj,[],1);
    for i = 1 : M
        [~,rank] = sortrows(PopObj(:,i));
        CrowdDis(rank(1),i)   = inf;
        CrowdDis(rank(end),i) = inf;
        for j = 2 : N-1
            CrowdDis(rank(j),i) = (PopObj(rank(j+1),i)-PopObj(rank(j-1),i))/(Fmax(i)-Fmin(i));
        end
    end
end


text:
classdef cluster

methods(Static)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Following functions are meant for assisting in
% setting up experiments.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function labels = labels_from_cluster_sizes(cluster_sizes)
        % total number of points
        S = sum(cluster_sizes);
        % Number of points in array
        K = numel(cluster_sizes);
        labels = zeros(1, S);
        i = 0;
        for k=1:K
            Sk  = cluster_sizes(k);
            for s=1:Sk
                i = i+1;
                labels(i) = k;
            end
        end
        return;
    end

    function [cluster_sizes, labels] = cluster_sizes_from_labels(labels)
        tbl = tabulate(labels);
        labels = tbl(:, 1)';
        cluster_sizes =  tbl(:, 2)';
    end

    function [start_indices, end_indices] = start_end_indices(cluster_sizes)
        % Returns start and end indices from cluster sizes
        start_indices = cumsum(cluster_sizes) + 1;
        start_indices = [1 start_indices(1:end-1)];
        end_indices = start_indices + cluster_sizes -1;
    end

    function result = clustering_error(estimated_labels, true_labels, num_clusters)
        % finds out an appropriate mapping between
        % true labels and estimated labels and
        % calculates the corresponding clustering error

        % make sure that both are row vectors
        if ~isrow(estimated_labels)
            estimated_labels = estimated_labels';
        end
        if ~isrow(true_labels)
            true_labels = true_labels';
        end
        if size(estimated_labels) ~= size(true_labels)
            error('Both estimated and true label arrays should have same size.');
        end
        % Let us consider all possible mappings
        possible_mappings = perms(1:num_clusters);
        % number of possible mappings
        num_mappings = size(possible_mappings,1 );
        % number of labels
        num_labels = size(estimated_labels, 2);
        % number of misclustered points for each mapping
        missed_points = zeros(num_mappings, 1);
        for j=1:num_mappings
            %cur_mapping = possible_mappings(j, :)
            mapped_labels = possible_mappings(j, true_labels);
            missed_points(j) = sum(estimated_labels ~= mapped_labels);
        end
        % missed_points
        % we choose the mapping with minimum mismatch
        [miss, index] = min(missed_points, [], 1);
        % final clustering error value
        result.num_labels = num_labels;
        result.num_missed_points = miss;
        result.error  = miss / num_labels;
        result.error_perc = result.error * 100;
        % corresponding mapping
        result.mapping = possible_mappings(index, :);
        result.mapped_labels = result.mapping(true_labels);
        result.misses = estimated_labels ~= result.mapped_labels;
    end

    function result = clustering_error_hungarian_mapping(estimated_labels, true_labels, num_clusters)
        % make sure that both are row vectors
        if ~isrow(estimated_labels)
            estimated_labels = estimated_labels';
        end
        if ~isrow(true_labels)
            true_labels = true_labels';
        end
        if size(estimated_labels) ~= size(true_labels)
            error('Both estimated and true label arrays should have same size.');
        end
        % number of labels
        num_labels = size(estimated_labels, 2);
        mapped_labels = bestMap(estimated_labels, true_labels)';
        miss = sum(estimated_labels ~= mapped_labels);
        % final clustering error value
        result.num_labels = num_labels;
        result.num_missed_points = miss;
        result.error  = miss / num_labels;
        result.error_perc = result.error * 100;
        result.mapped_labels = mapped_labels;
        result.misses = estimated_labels ~= result.mapped_labels;
    end

end

end


text:
function imo = cnn_get_im_flow_batch2(images, varargin)

opts.subTractFlow = 'off';
opts.nFramesPerVid = 1;
opts.numAugments = 1;
opts.frameSample = 'uniformly';
opts.flowDir = '';
opts.imageDir = '';
opts.temporalStride = 0;

opts.imageSize = [227, 227] ;
opts.border = [29, 29] ;
opts.averageImage = [] ;
opts.rgbVariance = [] ;
opts.augmentation = 'croponly' ;
opts.interpolation = 'bilinear' ;
opts.numAugments = 1 ;
opts.numThreads = 0 ;
opts.prefetch = false ;
opts.keepAspect = true;
opts.flowScales = [];
opts.cheapResize = 0;
opts.nFrameStack = 10;
opts.frameList = NaN;
opts.nFrames = [];
opts.subMedian = false;
opts.stretchAspect = 4/3 ;
opts.stretchScale = 1.2 ;
opts.fetchGPU = true ;
[opts, varargin] = vl_argparse(opts, varargin);

flowDir = opts.flowDir;
imgDir = opts.imageDir;
% prefetch is used to load images in a separate thread
prefetch = opts.prefetch & isempty(opts.frameList);
fetchOpts= {'numThreads', opts.numThreads};
if opts.fetchGPU
  fetchOpts{end+1} = 'Gpu' ;
end

switch opts.augmentation
  case 'croponly'
    tfs = [.5 ; .5 ; 0 ];
  case 'f5'
    tfs = [...
      .5 0 0 1 1 .5 0 0 1 1 ;
      .5 0 1 0 1 .5 0 1 0 1 ;
       0 0 0 0 0  1 1 1 1 1] ;
  case 'f25'
    [tx,ty] = meshgrid(linspace(0,1,5)) ;
    tfs = [tx(:)' ; ty(:)' ; zeros(1,numel(tx))] ;
    tfs_ = tfs ;
    tfs_(3,:) = 1 ;
    tfs = [tfs,tfs_] ;
  case 'f25noCtr'
    [tx1,ty1] = meshgrid(linspace(.75,1,20)) ;
    [tx2,ty2] = meshgrid(linspace(0,.25,20)) ;
    tx = [tx1 tx2];     ty = [ty1 ty2];
    tfs = [tx(:)' ; ty(:)' ; zeros(1,numel(tx))] ;
    tfs_ = tfs ;
    tfs_(3,:) = 1 ;
    tfs = [tfs,tfs_] ;
end

nStack = opts.imageSize(3);

if iscell(opts.frameList)
  im = vl_imreadjpeg(opts.frameList{1}, fetchOpts{:}) ;
  sampled_frame_nr = opts.frameList{2};
else
  sampleFrameLeftRight = floor(nStack/4); % divide by 4 because of left,right,u,v
  frameOffsets = [-sampleFrameLeftRight:sampleFrameLeftRight-1]';

  frames = cell(numel(images), nStack, opts.nFramesPerVid);
  frames_rgb = cell(numel(images), 1, opts.nFramesPerVid);

  sampled_frame_nr = cell(numel(images),1);

  for i=1:numel(images)
    vid_name = images{i};
    nFrames = opts.nFrames(i);

    if  strcmp(opts.frameSample, 'uniformly')
      sampleRate = max(floor((nFrames-nStack/2)/opts.nFramesPerVid),1);
      frameSamples = nStack/4+1:sampleRate:nFrames - nStack/4 ;
      opts.temporalStride = sampleRate;
      frameSamples = vl_colsubset(nStack/4+1:nFrames-nStack/4, opts.nFramesPerVid, 'uniform') ;
      opts.temporalStride =  frameSamples(2) - frameSamples(1);
    elseif strcmp(opts.frameSample, 'temporalStride')
      frameSamples = nStack/4+1:opts.temporalStride:nFrames-nStack/4 ;
      if length(frameSamples) < opts.nFrameStack,
          frameSamples = round(linspace(nStack/4+1, nFrames - nStack/4, opts.nFramesPerVid)) ;
          opts.temporalStride = frameSamples(2) - frameSamples(1);
      end
    elseif strcmp(opts.frameSample, 'random')
      frameSamples = randperm(nFrames-nStack/2)+nStack/4;
    elseif strcmp(opts.frameSample, 'temporalStrideRandom')
      frameSamples = nStack/4 +1:opts.temporalStride:nFrames - nStack/4 ;
      if length(frameSamples) < opts.nFrameStack,
          frameSamples = round(linspace(nStack/4+1, nFrames - nStack/4, opts.nFrameStack)) ;
          opts.temporalStride = frameSamples(2) - frameSamples(1);
      end
    end

    if length(frameSamples) < opts.nFramesPerVid,
        if length(frameSamples) > opts.nFrameStack
          frameSamples = frameSamples(1:length(frameSamples)-mod(length(frameSamples),opts.nFrameStack));
        end
        diff =  opts.nFramesPerVid - length(frameSamples);
        addFrames = 0;
        while diff > 0
          last_frame = min(frameSamples(end), max(nFrames - nStack/4 - opts.nFrameStack,nStack/4 ));
          if mod(addFrames,2) % add to the front
            addSamples = nStack/4+1:opts.temporalStride:nFrames - nStack/4;
            addSamples = addSamples(1: length(addSamples) - mod(length(addSamples),opts.nFrameStack));
            if length(addSamples) > diff, addSamples = addSamples(1:diff); end
          else % add to the back
            addSamples = fliplr(nFrames - nStack/4 : -opts.temporalStride: nStack/4+1);
            addSamples = addSamples(mod(length(addSamples),opts.nFrameStack)+1:length(addSamples));
            if length(addSamples) > diff, addSamples = addSamples(end-diff+1:end); end
          end

          if addFrames > 20
            addSamples = round(linspace(nStack/4+1, nFrames - nStack/4, opts.nFrameStack)) ;
          end
          frameSamples = [frameSamples addSamples];
          diff = opts.nFramesPerVid - length(frameSamples);
          opts.temporalStride = max(ceil(opts.temporalStride-1), 1);
          addFrames = addFrames+1;

        end
    end
    if length(frameSamples) > opts.nFramesPerVid
      if strcmp(opts.frameSample, 'temporalStride')
        s = fix((length(frameSamples)-opts.nFramesPerVid)/2);
      else % random
        s = randi(length(frameSamples)-opts.nFramesPerVid);
      end
        frameSamples = frameSamples(s+1:s+opts.nFramesPerVid);
    end

    for k = 1:opts.nFramesPerVid
        frames_rgb{i,1,k} = [vid_name 'frame' sprintf('%06d.jpg', frameSamples(k))] ;
    end

      frameSamples =  repmat(frameSamples,nStack/2,1) +  repmat(frameOffsets,1,size(frameSamples,2));
      for k = 1:opts.nFramesPerVid
        for j = 1:nStack/2
            frames{i,(j-1)*2+1, k} = ['u' filesep vid_name 'frame' sprintf('%06d.jpg', frameSamples(j,k)) ] ;
            frames{i,(j-1)*2+2, k} = ['v' frames{i,(j-1)*2+1, k}(2:end)];
        end
      end

      sampled_frame_nr{i} = frameSamples;
  end

    if iscell(opts.imageDir)
          imgDir = opts.imageDir{i};
          flowDir = opts.flowDir{i};
    end

  frames_rgb = strcat([imgDir filesep], frames_rgb);
  if ~isempty(flowDir)
    frames = strcat([flowDir filesep], frames);
    frames = cat(2, frames, frames_rgb);
  else
    frames = frames_rgb;
  end
  if opts.numThreads > 0
    if prefetch
      vl_imreadjpeg(frames, fetchOpts{:}, 'prefetch') ;
      imo = {frames sampled_frame_nr}  ;
      return ;
    end
    im = vl_imreadjpeg(frames, fetchOpts{:} ) ;
  end

end

if strcmp(opts.augmentation, 'none')

  szw = cellfun(@(x) size(x,2),im);
  szh = cellfun(@(x) size(x,1),im);

  h_min = min(szh(:));
  w_min =  min(szw(:));
  sz = [h_min w_min] ;

  sz = max(opts.imageSize(1:2), sz);
  sz = min(2*opts.imageSize(1:2), sz);

  scal = ([h_min w_min] ./ sz);

  imo =  zeros(sz(1), sz(2), opts.imageSize(3)+3, ...
            numel(images), 2 * opts.nFramesPerVid, 'single') ;
  if opts.fetchGPU
    imo = gpuArray(imo);
  end

  for i=1:numel(images)
    si = 1 ;
    for k = 1:opts.nFramesPerVid

      if numel(unique(szw)) > 1 || numel(unique(szh)) > 1
        for l=1:size(im,2)
            im{i,l,k} = im{i,l,k}(1:h_min,1:w_min,:);
        end
      end
        imt = cat(3, im{i,:,k}) ;

        if any(scal ~= 1)
          imo(:, :, :, i, si) = imresize(cat(3, im{i,:,k}),sz) ;
        else
          imo(:, :, :, i, si) = imt ;
        end
        imt = [];
        imo(:, :, :, i, si+1) =  imo(:, end:-1:1, :, i, si);
        imo(:, :, 1:2:nStack, i, si+1) = -imo(:, :, 1:2:nStack, i, si+1) + 255; %invert u if we flip

        si = si + 2 ;

    end
  end

  if opts.subMedian
    median_flow = median(imo(:,:,1:nStack),1);
    median_flow = median(median_flow,2);
    imo(:,:,1:nStack,:,:) = bsxfun(@minus, imo(:,:,1:nStack,:,:), median_flow ) ;
    imo(:,:,1:nStack,:,:) = bsxfun(@plus, imo(:,:,1:nStack,:,:), 128 ) ;
  end

  if ~isempty(opts.averageImage)
    opts.averageImage = mean(mean(opts.averageImage,1),2) ;
    imo = bsxfun(@minus, imo,opts.averageImage) ;
  end
  return;
end

% augment now
if exist('tfs', 'var')
  [~,transformations] = sort(rand(size(tfs,2), numel(images)*opts.nFramesPerVid), 1) ;
end

imo = ( zeros(opts.imageSize(1), opts.imageSize(2), opts.imageSize(3)+3, ...
            numel(images), opts.numAugments * opts.nFramesPerVid, 'single') ) ;

if opts.fetchGPU
  imo = gpuArray(imo);
end

for i=1:numel(images)
  si = 1 ;

  szw = cellfun(@(x) size(x,2),im);
  szh = cellfun(@(x) size(x,1),im);

  h_min = min(szh(:));
  w_min =  min(szw(:));

  if  strcmp( opts.augmentation, 'multiScaleRegular')
    reg_szs = [256, 224, 192, 168] ;
    sz(1) = reg_szs(randi(4)); sz(2) = reg_szs(randi(4));
  elseif strcmp( opts.augmentation, 'stretch')
    aspect = exp((2*rand-1) * log(opts.stretchAspect)) ;
    scale = exp((2*rand-1) * log(opts.stretchScale)) ;
    tw = opts.imageSize(2) * sqrt(aspect) * scale ;
    th = opts.imageSize(1) / sqrt(aspect) * scale ;
    reduce = min([w_min / tw, h_min / th, 1]) ;
    sz = round(reduce * [th ; tw]) ;
  else
    sz = round(min(opts.imageSize(1:2)' .* (.75+0.5*rand(2,1)), [h_min; w_min])) ; % 0.75 +- 0.5, not keep aspect
  end

  for k = 1:opts.nFramesPerVid

      if numel(unique(szw)) > 1 || numel(unique(szh)) > 1
        for l=1:size(im,2)
          im{i,l,k} = im{i,l,k}(1:h_min,1:w_min,:);
        end
      end

      imt = cat(3, im{i,:,k}) ;
      if opts.subMedian
          median_flow = median(imt(:,:,1:nStack),1);
          median_flow = median(median_flow,2);
          imt(:,:,1:nStack) = bsxfun(@minus, imt(:,:,1:nStack), median_flow ) ;
          imt(:,:,1:nStack) = bsxfun(@plus, imt(:,:,1:nStack), 128 ) ;
      end
%       imt = gpuArray(imt);

    w = size(imt,2) ;
    h = size(imt,1) ;
    if ~strcmp(opts.augmentation, 'uniform')
      if ~isempty(opts.rgbVariance) % colour jittering only in training case
        offset = zeros(size(imt));
        offset = bsxfun(@minus, offset, reshape(opts.rgbVariance * randn(opts.imageSize(3),1), 1,1,opts.imageSize(3))) ;
        imt = bsxfun(@minus, imt, offset) ;
      end

      for ai = 1:opts.numAugments
        switch opts.augmentation
          case 'stretch'
            dx = randi(w - sz(2) + 1 ) ;
            dy = randi(h - sz(1) + 1 ) ;
            flip = rand > 0.5 ;
          case 'multiScaleRegular'
            dy = [0 h-sz(1) 0 h-sz(1)  floor((h-sz(1)+1)/2)] + 1; % 4 corners & centre
            dx = [0 w-sz(2) w-sz(2) 0 floor((w-sz(2)+1)/2)] + 1;
            corner = randi(5);
            dx = dx(corner); dy = dy(corner); % pick one corner of the image
            flip = rand > 0.5 ;
          case 'f25noCtr'
            tf = tfs(:, transformations(mod(i+ai-1, numel(transformations)) + 1)) ;
            dx = floor((w - sz(2)) * tf(2)) + 1 ;
            dy = floor((h - sz(1)) * tf(1)) + 1 ;
            flip = tf(3) ;
          otherwise
            sz = opts.imageSize(1:2) ;
            tf = tfs(:, transformations(mod(ai-1, numel(transformations)) + 1)) ;
            dx = floor((w - sz(2)) * tf(2)) + 1 ;
            dy = floor((h - sz(1)) * tf(1)) + 1 ;
            flip = tf(3) ;
        end

        if opts.cheapResize
          sx = round(linspace(dx, sz(2)+dx-1, opts.imageSize(2))) ;
          sy = round(linspace(dy, sz(1)+dy-1, opts.imageSize(1))) ;
        else
          factor = [opts.imageSize(1)/sz(1) ...
              opts.imageSize(2)/sz(2)];

          if any(abs(factor - 1) > 0.0001)
            imt =   imresize(imt(dy:sz(1)+dy-1,dx:sz(2)+dx-1,:), [opts.imageSize(1:2)]);
          end

          sx = 1:opts.imageSize(2); sy = 1:opts.imageSize(1);
        end

        if flip
          sx = fliplr(sx) ;
          imo(:,:,:,i,si) = imt(sy,sx,:) ;
          imo(:,:,1:2:nStack,i,si) = -imt(sy,sx,1:2:nStack) + 255; %invert u if we flip
        else
          imo(:,:,:,i,si) = imt(sy,sx,:) ;
        end

        si = si + 1 ;
      end
    else

      w = size(imt,2) ; h = size(imt,1) ;

      indices_y = [0 h-opts.imageSize(1)] + 1;
      indices_x = [0 w-opts.imageSize(2)] + 1;
      center_y = floor(indices_y(2) / 2)+1;
      center_x = floor(indices_x(2) / 2)+1;

      if opts.numAugments == 6,  indices_y = center_y;
      elseif opts.numAugments == 2,  indices_x = [];   indices_y = [];
      elseif opts.numAugments ~= 10, error('only 6 or 10 uniform crops allowed');  end
        for y = indices_y
        for x = indices_x
          imo(:, :, :, i, si) = ...
              imt(y:y+opts.imageSize(1)-1, x:x+opts.imageSize(2)-1, :);

          imo(:, :, :, i, si+1) = imo(:, end:-1:1, :, i, si);
          imo(:, :, 1:2:nStack, i, si+1) = -imo(:, end:-1:1, 1:2:nStack, i, si) + 255; %invert u if we flip

          si = si + 2 ;
        end
        end
        imo(:,:,:, i,si) = imt(center_y:center_y+opts.imageSize(1)-1,center_x:center_x+opts.imageSize(2)-1,:);

        imo(:,:,:, i,si+1) = imo(:, end:-1:1, :, i, si);
        imo(:,:,1:2:nStack, i,si+1) = -imo(:, end:-1:1, 1:2:nStack, i, si) + 255; %invert u if we flip

        si = si + 2;
    end
  end
end

if ~isempty(opts.averageImage)
  imo = bsxfun(@minus, imo, opts.averageImage) ;
end

end

author: feichtenhofer

repo: st-resnet


text:
function [S] = cell2nestedstruct(C)
%
%  S = CELL2NESTEDSTRUCT(C) converts an MxN cell array C into a
%     nested 1x1 structure array S
%
%     Examples:
%
%       C = {'a','b',1;  'qqq','',cell(6,2);   'aa','',[1 3]};
%       [S] = cell2nestedstruct(C)
%
%       S =
%           a: [1x1 struct]
%           qqq: {6x2 cell}
%           aa: [1 3]
%
%
%       C = {'water'    'a'     ''     ''     ''    111;...
%            'water'    'b'     ''     ''     ''    222;...
%            'melon'     ''     ''     ''     ''    [];...
%            'a'        'b'    'c'    'd'    'e'    'aaa';...
%            'q'        'r'     ''     ''     ''    cell(6,2)};
%       [S] = cell2nestedstruct(C);
%       [C2] = nestedstruct2cell(S);
%       isequal(C,C2)
%
%       ans =
%
%           1
%
%
%     Notes:
%       1. The fields in the first N-1 rows of C must be be left-justified.
%       2. The data must be in the right-most column.
%       3. The structure is produced row-by-row using the general syntax:
%                   S.row1.row2. ... rowN-1 = rowN;
%           which means that a sigle row, for example:
%                   {'water'    'a'     ''     ''     ''    111}
%           is equivalent to entering:
%                   S.water.a = 111;
%
%
%
%     See also CELL2STRUCT, STRUCT2CELL
%
%     See also NESTEDSTRUCT2CELL (available at Matlab
%     Central File Exchange)
%
%Author:
%
% Todd C Pataky (0todd0@gmail.com)  26-May-2006

%Programming notes:
%1. Currently working on implementating of MxN structural array, but
%   blank fields, for example, make the N-D case difficult to deal with.
%   Suggestions appreciated.

%%%%%%%%%%%%%%%%%%%%
%DATA ORGANIZATION CHECKS
%%%%%%%%%%%%%%%%%%%%
if ~iscell(C) | ndims(C)~=2 | size(C,2)<2
    error('Must input an MxN cell array (N>1).')
end
if ~iscellstr(C(:,1:end-1))
    error('First %.0f columns of cell array must contain strings.',size(C,2)-1)
end
%%%%%%%%%%%%%%%%%%%%
%check that empty cells are to the right
%%%%%%%%%%%%%%%%%%%%
c = C(:,1:end-1);
L = ones(size(c));
ind = strmatch({' '},c);
L(ind) = 0;   %logical array indicating non-blank entries
for k=1:size(L,1)
    f = find(L(k,:));
    if ~isequal(f,1:length(f))
        error('All non-blank cells in row %.0f must be left-justified',k)
    end
end
%%%%%%%%%%%%%%%%%%%%
%check that there are no duplicate rows
%%%%%%%%%%%%%%%%%%%%
for k = 1:size(C,1)-1
    c = C(k,1:end-1);
    for kk = (k+1):size(C,1);
        cc = C(kk,1:end-1);
        if isequal(c,cc)
            error('Rows %.0f and %.0f are identical branches',k,kk);
        end
    end
end
%%%%%%%%%%%%%%%%%%%%
%check that branches go to the same depth (i.e. check for branch overwriting)
%%%%%%%%%%%%%%%%%%%%
hereCheckOverWrite(C(:,1:end-1));    %see subfunction below

%%%%%%%%%%%%%%%%%%%%
%ASSEMBLE STRUCTURE ROW-BY-ROW
%%%%%%%%%%%%%%%%%%%%
for k=1:size(C,1)
    a = C(k,:);
    ind = strmatch('',a(1:end-1),'exact');
    a(ind) = [];

    s = a(1:end-1);
    data = a{end};

    s = [s; repmat({'.'},size(s))];
    s = s(:)';
    s = cat(2,s{1:end-1});

    eval(['S.',s,'= data;']);
end

%%%%%%%%%%%%%%%%%%%%
%SUBFUNCTIONS
%%%%%%%%%%%%%%%%%%%%
function hereCheckOverWrite(C)
%%%%%%%%%%%%%%%%%%%%
if size(C,2)==1
    return
end
%%%%%%%%%%%%%%%%%%%%
%if a zero stem has the same root as any other branch, error
%%%%%%%%%%%%%%%%%%%%
row = strmatch({' '},C(:,end));
for k=1:length(row)
    root0 = C(row(k),1:end-1);
    for kk=1:size(C,1)
        if kk~=row(k)
            root1 = C(kk,1:end-1);
            if isequal(root0,root1)
                error('Overwriting error. A stem from row %.0f is empty, but rows %.0f and %.0f have the same roots.',...
                    row(k) , min(row(k),kk) , max(row(k),kk) );
            end
        end
    end
end
hereCheckOverWrite(C(:,1:end-1));   %recurse

repo: mcbench-benchmarks


text:
function [inputMovie] = addDroppedFramesToMovie(inputMovie,droppedFrames,varargin)
	% Fixes movie by adding dropped frames with the mean of the movie (to reduce impact on cell extraction algorithms).
	% Biafra Ahanonu
	% started: 2016.10.04 [20:31:15]
	% inputs
		% inputMovie - matrix dims are [X Y t] - where t = number of time points
		% path to inscopix file or list of dropped frames
	% outputs
		% inputMovie with dropped frames added back in.

	[inputMovie] = ciapkg.video.addDroppedFramesToMovie(inputMovie,droppedFrames,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function iq = acquireIqVector(mxa_ip)

% MATLAB/MXA example 7
% Getting IQ data using the MXA driver and plot display

% Version: 1.0
% Date: Sep 11, 2006
% 2006 Agilent Technologies, Inc.

% TCPIP parameters
mxa_port = 5025;

% MXA Interface creation and connection opening
mxa_if = tcpip(mxa_ip,mxa_port);
mxa = icdevice('agilent_mxa.mdd', mxa_if);
connect(mxa,'object')

% Get IQ data
iq = invoke(mxa,'getIQ');

% Close the MXA connection and clean up
disconnect(mxa);
delete(mxa);
clear mxa;

author: Sable


text:
%%% void show_video(struct)
%
function show_video(video)
  for i = 1 : video.nrFramesTotal
    frame = video.frames(i).cdata;
    imshow(frame,[],'InitialMagnification','fit');
    disp(i);
    pause(0.01);
  end
end

author: andrewssobral

repo: lrslibrary


text:
%
% Usage:   [V [val_regularizer]]=mexProximalGraph(U,graph,param);
%
% Name: mexProximalGraph
%
% Description: mexProximalGraph computes a proximal operator. Depending
%         on the value of param.regul, it computes
%
%         Given an input matrix U=[u^1,\ldots,u^n], and a set of groups G,
%         it computes a matrix V=[v^1,\ldots,v^n] such that
%
%         if param.regul='graph'
%         for every column u of U, it computes a column v of V solving
%             argmin 0.5||u-v||_2^2 + lambda\sum_{g \in G} \eta_g||v_g||_inf
%
%         if param.regul='graph+ridge'
%         for every column u of U, it computes a column v of V solving
%             argmin 0.5||u-v||_2^2 + lambda\sum_{g \in G} \eta_g||v_g||_inf + lambda_2||v||_2^2
%
%
%         if param.regul='multi-task-graph'
%            V=argmin 0.5||U-V||_F^2 + lambda \sum_{i=1}^n\sum_{g \in G} \eta_g||v^i_g||_inf + ...
%                                                lambda_2 \sum_{g \in G} \eta_g max_{j in g}||V_j||_{inf}
%
%         it can also be used with any regularization addressed by mexProximalFlat
%
%         for all these regularizations, it is possible to enforce non-negativity constraints
%         with the option param.pos, and to prevent the last row of U to be regularized, with
%         the option param.intercept
%
% Inputs: U:  double p x n matrix   (input signals)
%               m is the signal size
%         graph: struct
%               with three fields, eta_g, groups, and groups_var
%
%               The first fields sets the weights for every group
%                  graph.eta_g            double N vector
%
%               The next field sets inclusion relations between groups
%               (but not between groups and variables):
%                  graph.groups           sparse (double or boolean) N x N matrix
%                  the (i,j) entry is non-zero if and only if i is different than j and
%                  gi is included in gj.
%
%               The next field sets inclusion relations between groups and variables
%                  graph.groups_var       sparse (double or boolean) p x N matrix
%                  the (i,j) entry is non-zero if and only if the variable i is included
%                  in gj, but not in any children of gj.
%
%               examples are given in test_ProximalGraph.m
%
%         param: struct
%               param.lambda  (regularization parameter)
%               param.regul (choice of regularization, see above)
%               param.lambda2  (optional, regularization parameter)
%               param.lambda3  (optional, regularization parameter)
%               param.verbose (optional, verbosity level, false by default)
%               param.intercept (optional, last row of U is not regularized,
%                 false by default)
%               param.pos (optional, adds positivity constraints on the
%                 coefficients, false by default)
%               param.numThreads (optional, number of threads for exploiting
%                 multi-core / multi-cpus. By default, it takes the value -1,
%                 which automatically selects all the available CPUs/cores).
%
% Output: V: double p x n matrix (output coefficients)
%         val_regularizer: double 1 x n vector (value of the regularization
%         term at the optimum).
%
% Author: Julien Mairal, 2010


text:
function v = dquat2pos(dq)

% DQUAT2POS  Transforms a double quaternion point position representation
%            into a vector representation.
%
%     V = DQUAT2POS(DQ) transforms the double quaternion representation of
%        a point position (in the euclidean space) into a vector v, which
%        represents the point coordinates. DQ is either a vector of size 8
%        or an array of size 8*N (each column represents a point position
%        dual quaternion) where N is the number of point locations. V is a
%        vector of size 3 or an array of size 3*N depending on the input
%        format.
%
% See also POS2DQUAT, DQUAT2VEL, DQUAT2LINE, DQUAT2LINEVEL

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

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

% check point position dual quaternion
tol = 1e-6;
[maxval1,imax1] = max(abs(dq(1,:)-1));
[maxval2,imax2] = max(max(abs(dq(2:5,:))));
imax = union(imax1,imax2);
if maxval1 > tol || maxval2 > tol
    warning('DualQuaternion:dquat2pos:wrongFormat',...
        'At least one dual quaternion is not a point position dual quaternion (tol = %.1e).\n Indices of max values: %d',...
        tol,imax);
end

% extraction of the point position coordinates
v = dq(6:8,:);


text:
function [ grid_num, point_num ] = cc_grids_minmax_size ( dim_num, q_min, ...
  q_max )

%*****************************************************************************80
%
%% CC_GRIDS_MINMAX_SIZE returns sizes for the CC_GRIDS_MINMAX.
%
%  Discussion:
%
%    This routine can be used to determine the necessary size to be
%    allocated to arrays GRID_ORDER and GRID_POINT in a call to
%    CC_GRIDS_MINMAX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 October 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer Q_MIN, Q_MAX, the minimum and maximum values of
%    Q, the sum of the orders in each spatial coordinate.
%
%    Output, integer GRID_NUM, the number of Clenshaw Curtis
%    grids whose Q value is between Q_MIN and Q_MAX.
%
%    Output, integer POINT_NUM, the total number of points in the grids.
%

%
%  Determine the total number of points that will be generated
%  by "going through the motions".
%
  point_num = 0;
  grid_num = 0;

  for q = q_min : q_max

    more = 0;
    order_1d = [];

    while ( 1 )

      [ order_1d, more ] = compnz_next ( q, dim_num, order_1d, more );

      order_nd = prod ( order_1d(1:dim_num) );

      point_num = point_num + order_nd;

      grid_num = grid_num + 1;

      if ( ~more )
        break
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function inspect_pull1970

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_databrowser

cfg = [];
cfg.dataset = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.ds');
cfg = ft_definetrial(cfg);
cfg.trl = cfg.trl(1:20,:);
cfg.trl(:,3) = cfg.trl(:,1)-1;
cfg.demean = 'yes';
cfg.channel = 'MEG';
data = ft_preprocessing(cfg);

%%

cfg = [];
cfg.plotevents = 'no';
cfg.viewmode = 'vertical';
cfg.artfctdef.other.artifact = zeros(0,2);
cfg.artfctdef.blink.artifact = zeros(0,2);
cfg.artfctdef.movement.artifact = zeros(0,2);
ft_databrowser(cfg, data)

%%
% all viewmodes can work with component data

cfg = [];
cfg.method = 'pca';
comp = ft_componentanalysis(cfg, data);

%%

cfg = [];
cfg.plotevents = 'no';
cfg.viewmode = 'vertical';
ft_databrowser(cfg, comp)

%%

cfg = [];
cfg.plotevents = 'no';
cfg.viewmode = 'butterfly';
ft_databrowser(cfg, comp)

%%

cfg = [];
cfg.plotevents = 'no';
cfg.viewmode = 'component';
ft_databrowser(cfg, comp)

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip