text:
function [Data_PCA,F,angles] = locPCA_affnity_fast(data,index_mat,params)
n_neighbors = size(index_mat,1)%params.n_neighbors;
%n_neighbors = params.n_neighbors;
int_dim = params.int_dim;
N = size( data, 2 );
D = size( data, 1 );
Data_eig=cell(1,N);
% retain = params.m*params.n_neighbors;
% n_neighbors=retain
%n_neighbors=retain;
%n_neighbors = params.n_neighbors;
%[index_mat]=K_nearest_neighbors(data,n_neighbors);
%F=sparse(zeros(N,N));
F=spalloc(N,N,300*N);
%F=(zeros(N,N));
[m,N] = size(data);  % m is the dimensionality of the input sample points.

eig_dim1=zeros(m,N);
if length(n_neighbors)==1
    K = repmat(n_neighbors,[1,N]);
end;

parfor i=1:N

    % Compute the d largest right singular eigenvectors of the centered matrix
    Ii = index_mat(:,i); ki = K(i);
    Xi = data(:,Ii)-repmat(mean(data(:,Ii),2),[1,ki]);

    %construct the matrix Y
    Y = Xi' / sqrt(n_neighbors-1);
    %[U,Sing,V] = svd( Xi );
    [U,Sing,V] = svd( Y,'econ');
    %   Zi = U(:,1:d)*U(1:d,:)*Xi;
    %  [signals,PC,V] = pca2(Xi) ;
    % Data_PCA(1:D,i) = data(1:D,i);
    %Data_eig{i}=V(:,1:d);
    Data_eig{i}= V(:,1:int_dim);
    eig_dim1(:,i)=(Data_eig{i}(1:m));
    %Data_eig{i}=V(:,:1:d);

end

angles=spalloc(N,N,300*N);
%h = waitbar(0,'Local PCA...');
if params.int_dim>1
for i=1:N
    waitbar(i/N);
    normal_i=Data_eig{i}';
    temp_nei_i=index_mat(:,i);
    for j=1:length(temp_nei_i)
        idxNeig = temp_nei_i(j);
        normal_j=(Data_eig{idxNeig})';
        [theta] = max(subspaceangle(normal_i',normal_j'));
        angles(idxNeig,i) = theta;
        angles(i,idxNeig) = theta;
       % theta = subspace(normal_i',normal_j');
        F(idxNeig,i)=cos(theta).^params.powerCos;
        % F(j,i)=abs(acosd(theta));
        %F(j,i) = sind(F(j,i));
        %F(j,i)=exp((-F(j,i).^2)./1.0);
        F(i,idxNeig)=F(idxNeig,i);
    end
end
elseif params.int_dim==1
%  % s = sqrt( sum( ( data_noise(:,I)-data_noise(:,J) ).^2, 1) );
% for i=1:N
%     waitbar(i/N);
%     normal_i   =   eig_dim1(:,i)';
%     temp_nei_i =   index_mat(:,i);
%     normal_j=( eig_dim1(:, temp_nei_i));
%     F(i, temp_nei_i)=normal_i*normal_j;
%     F(temp_nei_i,i)=F(i,temp_nei_i);
% end
% [I]=find(F);
% s1 =  (acosd(F(I)));
% s2=  180- (acosd(F(I)));
% s=min(s1,s2);
% angles(I)=s;
% [IDX_vote,S]=knnsearch(data',data','k',retain+1,'distance','euclidean');
%  IDX_vote=IDX_vote(:,2:end);
%  S=S(:,2:end);
%  IDX_vote=IDX_vote';
%  S=S';
%  S = S(1:n_neighbors,:);
 IDX2 = index_mat;%(1:n_neighbors,:);
 a1=1:N;
 a2=repmat(a1,n_neighbors);
 a2=a2(1:n_neighbors,1:N);
 a3=a2(:);
 I=a3;
 J=IDX2(:);
 r1=eig_dim1(:,I);
 r2=eig_dim1(:,J);
 %A1=zeros(D,N); A2=zeros(D,N);
 a3=r1(1,:);
 a4=r2(1,:);
 a5=r1(2,:);
 a6=r2(2,:);
 A=r1.*r2;
 A=sum(A);
 %A=(a3.*a4)+(a5.*a6);
 s=A(:);

 w_temp=abs(s);
 waff = sparse(I,J,w_temp,N,N);
 W_k=waff;
 W_k= (W_k + W_k')./2;
 F=sparse(W_k);
% for i=1:N
%     waitbar(i/N);
%     normal_i   =   eig_dim1(:,i)';
%     temp_nei_i =   index_mat(:,i);
%     normal_j=( eig_dim1(:, temp_nei_i));
%     F(i, temp_nei_i)=normal_i*normal_j;
%     F(temp_nei_i,i)=F(i,temp_nei_i);
% end
% [I]=find(F);
% s1 =  (acosd(F(I)));
% s2=  180- (acosd(F(I)));
% s=min(s1,s2);
% angles(I)=s;
% end
end

%close(h);
Data_PCA = Data_eig;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function flag=istfraff(method);
% flag=istfr2(method) returns true is method is an affine
% time frequency representation.
%	See also istfr1, istfr2.

%	F. Auger, may 98
%	Copyright (c) CNRS - France 1998.
%
%  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

method=upper(method);
if strcmp(method,'TFRASPW' ) | strcmp(method,'TFRSCALO') | ...
   strcmp(method,'TFRDFLA' ) | strcmp(method,'TFRSPAW' ) | ...
   strcmp(method,'TFRUNTER') | strcmp(method,'TFRBERT' ) | ...
   strcmp(method,'TFRSPBK' ),
 flag=1;
else
 flag=0;
end;

author: HeLiangHIT


text:
%% This file is the ODE file used to simulate the Michaelis-Menten kinetics.
% Date: 04/24/2019
% Coded By: K

% Here Vmax is rhe maximum rate of reaction time
% Km is the concentration of half-maximal reaction rate
%

%% ODE functions
function dx=MMK_ODE(t,x,jx,Vmax,Km)
dx=jx-(Vmax*x(1,:)./(Km+x(1,:)));

author: dynamicslab

repo: SINDy-PI


text:
function [corrTime] = check_t(time)

% SYNTAX:
%   [corrTime] = check_t(time);
%
% INPUT:
%   time = GPS time
%
% OUTPUT:
%   corrTime = corrected GPS time
%
% DESCRIPTION:
%   Function accounting for beginning or end of week crossover. From the
%   Interface Specification document revision E (IS-GPS-200E), page 93.

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

%----------------------------------------------------------------------------------------------
%                           goGPS v0.4.3
%
% Copyright (C) Kai Borre
% Kai Borre 04-01-96
%
% Adapted by Mirko Reguzzoni, Eugenio Realini, 2009
%----------------------------------------------------------------------------------------------

half_week = 302400;     % seconds

corrTime = time;
corrTime(time > half_week) =time(time > half_week) - 2*half_week;
corrTime(time < - half_week) =time(time < - half_week) + 2*half_week;
% if time > half_week
%     corrTime = time - 2*half_week;
% elseif time < -half_week
%     corrTime = time + 2*half_week;
% end

author: goGPS-Project


text:
function nDefined = notDefined(varString)
%
% nDefined= notDefined( varString )
%
% Determines if a variable is defined in the workspace.
% A variable is defined if (a) it exists and (b) it is not empty.
%
% nDefined: 1 (true) if the variable is not defined in the calling workspace
%           0 (false) if the variable is defined in the calling workspace
%
%  This routine replaces calls of the form:
%    if ~exist('varname','var') || isempty(xxx) with the call
%    if notDefined('varname')
%

if (~ischar(varString)), error('Varible name must be a string'); end

nDefined = 1;

str = sprintf('''%s''',varString);
cmd1 = ['~exist(' str ',''var'') == 1'];
cmd2 = ['isempty(',varString ') == 1'];

% If either of these conditions holds, Defined == true
nDefined = evalin('caller',cmd1);     % Check that the variable exists in the caller space
if nDefined, return;                  % If it does not, return with a status of 0
else
    nDefined = evalin('caller',cmd2); % Check if the variable is empty in the caller space
    if nDefined return;
    else nDefined = 0;

    end
end

return;

author: dustinstansbury

repo: medal


text:
function value = r8mat_nonzeros ( m, n, a )

%*****************************************************************************80
%
%% R8MAT_NONZEROS counts the nonzeros in an R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 August 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), the matrix.
%
%    Output, integer VALUE, the number of nonzeros.
%
  value = length ( find ( a ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% DEMCMU35GPLVMVARGPLVM3 Run variational GPLVM with dynamics on CMU35 data.
%
% COPYRIGHT :  Andreas C. Damianou, Michalis K. Titsias, 2011
%
% SEEALSO : demCmu35vargplvmReconstructTaylor.m
% VARGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

% Define constants (in a manner that allows other scripts to parametrize
% this one).
if ~exist('experimentNo') ,  experimentNo = 404;      end
if ~exist('itNo')         ,  itNo = [200 1000 1300];              end     % Default: 2000
if ~exist('indPoints')    ,  indPoints = 100;          end     % Default:
if ~exist('latentDim')    ,  latentDim = 9;          end
% Set to 1 to use dynamics or to 0 to use the standard var-GPLVM
if ~exist('dynUsed')      ,  dynUsed = 1;             end

if ~exist('fixedBetaIters'), fixedBetaIters = 50;      end     % DEFAULT:
% Set to 1 to tie the inducing points with the latent vars. X
if ~exist('fixInd')        ,     fixInd = 0;                             end
if ~exist('baseKern')   ,     baseKern = {'rbfard2', 'white'}; end
if ~exist('dynamicKern')   ,     dynamicKern = {'rbf', 'bias', 'white'}; end
if ~exist('reconstrIters') ,     reconstrIters = 2000;                   end
if ~exist('learnVariance'),     learnVariance =0;   end
%if ~exist('vardistCovarsMult'),  vardistCovarsMult=1.3;                  end
if ~exist('initX'),     initX ='ppca';   end
if ~exist('doReconstr'),     doReconstr=1;   end
if ~exist('optimiser', 'var'), optimiser = 'scg'; end
if ~exist('testReoptimise', 'var'), testReoptimise = 0; end
if ~exist('parallelVardist', 'var'), parallelVardist = false; end

% Get the sequence numbers.
[Y, lbls] = lvmLoadData('cmu35WalkJog');
seq = cumsum(sum(lbls)) - [1:31];

dataSetName = 'cmu35gplvm';

% load data
[Y, lbls, Ytest, lblstest] = lvmLoadData(dataSetName);

fprintf(1,'\n#----------------------------------------------------\n');
fprintf(1,'# Dataset: %s\n',dataSetName);
fprintf(1,'# ExperimentNo: %d\n', experimentNo);
fprintf(1,'# Inducing points: %d\n',indPoints);
fprintf(1,'# Latent dimensions: %d\n',latentDim);
fprintf(1,'# Iterations (with/without fixed Beta): %d / %s\n',fixedBetaIters, num2str(itNo));
%fprintf(1,'# VardistCovarsMult: %d\n',vardistCovarsMult);
fprintf(1,'# Reconstruction iterations: %d\n', reconstrIters);
fprintf(1,'# Dataset size used (train) : %d \n', size(Y,1));
fprintf(1,'# Tie Inducing points: %d\n',fixInd);
fprintf(1,'# InitX: %s\n',initX);
fprintf(1,'# Base kern: ');
disp(baseKern);
fprintf(1,'# Dynamics used: %d\n', dynUsed);
if dynUsed
    fprintf(1,'# Dynamics kern: ');
    disp(dynamicKern);
end

fprintf(1,'#----------------------------------------------------\n');

% %%% Uncomment this code to get a subset of the sequences
% seqFrom=2;
% seqEnd=4;
%
% if seqFrom ~= 1
%     Yfrom = seq(seqFrom-1)+1;
% else
%     Yfrom = 1;
% end
% Yend=seq(seqEnd);
% Y=Y(Yfrom:Yend,:);
% seq=seq(seqFrom:seqEnd);
% seq=seq - ones(1,length(seq)).*(Yfrom-1);

% Fix times:
prevSeq = 1;
timeStampsTraining = [];
dt=0.05;
for i=1:length(seq)
    t = ([0:(seq(i)-prevSeq)].*dt)';
    prevSeq = seq(i)+1;
    timeStampsTraining = [timeStampsTraining ;t];
end;
timeStampsTest = ([0:size(Ytest,1)-1].*dt)';

% Set up model
options = vargplvmOptions('dtcvar');
%options.kern = {'rbfard2', 'bias', 'white'};
options.kern = baseKern;
options.numActive = indPoints; % DEfault: 100

options.optimiser = optimiser;

d = size(Y, 2);
fprintf(1,'# Creating the model...\n');
if fixInd
    options.fixInducing=1;
    options.fixIndices=1:indPoints;
end
model = vargplvmCreate(latentDim, d, Y, options);
model = vargplvmParamInit(model, model.m, model.X);
model.beta=1/(0.01*var(model.m(:)));
model.vardist.covars = 0.5*ones(size(model.vardist.covars)) + 0.001*randn(size(model.vardist.covars));

model.reconstrIters = reconstrIters;

%-------- Add dynamics to the model -----
if dynUsed
    optionsDyn.type = 'vargpTime';
    optionsDyn.t=timeStampsTraining;
    optionsDyn.inverseWidth=30;
    %   optionsDyn.vardistCovars = vardistCovarsMult;
    optionsDyn.seq = seq;
    optionsDyn.learnVariance = learnVariance;
    optionsDyn.initX = initX;
    optionsDyn.regularizeMeans = 0;

    % Dynamic kernel:
    kern = kernCreate(t, dynamicKern); % Default: {'rbf','white','bias'}

    if strcmp(kern.comp{2}.type, 'white')
        kern.comp{2}.variance = 1e-2; % Usual values: 1e-1, 1e-3
    end

    if strcmp(kern.comp{2}.type, 'whitefixed')
        if ~exist('whiteVar')
            whiteVar = 1e-6;
        end
        kern.comp{2}.variance = whiteVar;
        fprintf(1,'# fixedwhite variance: %d\n',whiteVar);
    end

    if strcmp(kern.comp{1}.type, 'rbfperiodic')
        if exist('periodicPeriod')
            kern.comp{1}.period = periodicPeriod;
        end
        fprintf(1,'# periodic period: %d\n',kern.comp{1}.period);
    end

    % The following is related to the expected number of
    % zero-crossings.(larger inv.width numerator, rougher func)
    if ~strcmp(kern.comp{1}.type,'ou')
        kern.comp{1}.inverseWidth = optionsDyn.inverseWidth./(((max(t)-min(t))).^2);
        kern.comp{1}.variance = 1;
    end
    optionsDyn.kern = kern;

    if exist('vardistCovarsMult')
        optionsDyn.vardistCovars = vardistCovarsMult;
    end

    % Fill in with default values whatever is not already set
    optionsDyn = vargplvmOptionsDyn(optionsDyn);
    model = vargplvmAddDynamics(model, 'vargpTime', optionsDyn, optionsDyn.t, 0, 0,optionsDyn.seq);

    fprintf(1,'# Further calibration of the initial parameters...\n');
    model = vargplvmInitDynamics(model,optionsDyn);
end

if parallelVardist
    model.vardist.parallel = true;
end

modelInit = model;

% do not learn beta for few iterations for intitilization
if fixedBetaIters > 0
    model.learnBeta = 0; model.initVardist = 1; model.learnSigmaf = 0;
    display = 1;
    fprintf(1,'# Intitiliazing the model (fixed beta) %d iterations...\n',fixedBetaIters);
    model = vargplvmOptimise(model, display, fixedBetaIters);
    model.fixedBetaIters = fixedBetaIters;
    disp('# Saving model after optimising beta...')
    try
        vargplvmWriteResult(model, model.type, dataSetName, experimentNo);
    catch e
        warning('Error while saving model: ')
        fprintf('%s\n',e.message)
    end
end

% Optimise the model.
display = 1;
model.learnBeta = 1; model.initVardist = 0; model.learnSigmaf = 1;
model.iters=0;

for i=1:length(itNo)
    iters = itNo(i); % Default: 1000
    fprintf(1,'# Optimising the model for %d iterations (session %d)...\n',iters,i);
    model = vargplvmOptimise(model, display, iters);
    fprintf('1/b=%.4d var(model.m)=%.4d\n', 1/model.beta, var(model.m(:)));
    model.iters = model.iters + iters;
    model.date = date;
    % Save the results.
    fprintf(1,'# Saving model after doing %d iterations',iters)
    try
        vargplvmWriteResult(model, model.type, dataSetName, experimentNo);
    catch e
        warning('Error while saving model: ')
        fprintf('%s\n',e.message)
    end
end

%bar(model.kern.comp{1}.inputScales)
%prefix = 'scales';
%saveAllOpenFigures(['Results/CMU/NEW/' num2str(experimentNo) '/'], prefix,1)

if testReoptimise
    model.dynamics.reoptimise=1;
end

% Reconstruction can also be done separately calling demCmu35vargplvmReconstructTaylor
if doReconstr
    %---- RECONSTRUCTION ---%
    save 'TEMPExperimentNo.mat' 'experimentNo'
    clear
    load 'TEMPExperimentNo'
    delete TEMPExperimentNo.mat
    fprintf('# Taylor Reconstruction for expNo:%d\n', experimentNo);
    demCmu35vargplvmReconstructTaylor
    %prefix = [num2str(experimentNo) '/'];
    %saveAllOpenFigures('Results/CMU/NEW/', prefix,1)
end

%% ----------- VISUALISATION
%{
clear;ca
%load('/home/andreas/SoftwareNotDrpBox/oldMatFiles/vargplvm/2/demCmu35gplvmVargplvm33.mat');
load('/home/andreas/Dropbox/_PhD/Software/github/private/vargplvm/matlab/private/priv/matFiles/Nips11/demCmu35gplvmVargplvm33.mat');
baseDir = datasetsDirectory;
dirSep = filesep;
[Y1, lbls, Ytest, lblstest] = lvmLoadData('cmu35WalkJog', [], baseDir);
skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35.asf']);
[tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35_01.amc'], skel);

[Y, lbls] = lvmLoadData('cmu35WalkJog');
seq = cumsum(sum(lbls)) - [1:31];
dataSetName = 'cmu35gplvm';
% load data
[Y, lbls, Ytest, lblstest] = lvmLoadData(dataSetName);
channels = demCmu35VargplvmLoadChannels(Y, skel);

modelTmp = model;
modelTmp.dynamics = [];
lvmVisualise(modelTmp, [], 'skelVargplvmVisualise', 'skelVargplvmModify', channels, skel)
figure; vargplvmShowScales(model)
%}
%------------------------

author: SheffieldML

repo: vargplvm


text:
function s=get_name(a)

s=[get_name(a.algorithm) ' (' get_name(a.child) ')'];  %% print name of algorithm
eval_name                   %% print values of hyperparameters

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function val = viewGetSession(vw,param,varargin)
% Get data from various view structures
%
% This function is wrapped by viewGet. It should not be called by anything
% else other than viewGet.
%
% This function retrieves information from the view that relates to a
% specific component of the application.
%
% We assume that input comes to us already fixed and does not need to be
% formatted again.

if notDefined('vw'), vw = getCurView; end
if notDefined('param'), error('No parameter defined'); end

mrGlobals;
val = [];

switch param

    case 'homedir'
        % Return full path to directory.
        %   homedir = viewGet(vw, 'Home Directory');
        val = HOMEDIR;
    case 'sessionname'
        % Retrun name of session, such as 'BW090616-8Bars-14deg'.
        %   sessionName = viewGet(vw, 'session name');
        val = mrSESSION.sessionCode;
    case 'subject'
        % Return name of subject, such as 'Wandell'
        %   subject = viewGet(vw, 'subject')
        val = mrSESSION.subject;
    case 'name'
        % Return name of view, such as 'INPLANE{1}'
        %    name = viewGet(vw, 'view name');
        val = vw.name;
    case 'annotation'
        % Return description of currently selected scan (string, such as
        % '14 Deg 8 Bars with blanks')
        %       annotation = viewGet(vw, 'annotation');
        if length(varargin) >= 1,   scan = varargin{1};
        else                        scan = viewGet(vw, 'CurScan'); end
        dt = viewGet(vw, 'DTStruct');
        val = dt.scanParams(scan).annotation;
    case 'annotations'
        % Return a cell array with descriptions of each scan in current
        % dataType
        %   annotations = viewGet(vw, 'annotations');
        dt  = viewGet(vw, 'DT Struct');
        val = {dt.scanParams.annotation};
    case 'viewtype'
        % Return the view type ('Gray', 'Volume', 'Inplane', 'Flat', etc)
        %   viewType = viewGet(vw, 'View Type');
        val = vw.viewType;
    case 'subdir'
        % Return the sub directory name (not the full path) with data for
        % current view
        %   subdir = viewGet(vw, 'sub directory');
        val = vw.subdir;
    case 'viewdir'
        % Return the complete path combination of homedir and subdir
        % Will then try to create them if not already created
        val = fullfile(viewGet(vw,'Home Directory'),viewGet(vw,'Sub Directory'));
        if ~exist(val,'dir')
            fprintf('Trying to make %s...',str);
            try
                [~, message] = mkdir(val);
            catch
                fprintf('Whoops, didn''t succed. Maybe a permissions problem?');
                fprintf('\n Message: %s',message);
            end
            fprintf('\n');
        end
    case 'curscan'
        % Return the currently selected scan number
        %   curscan = viewget(vw, 'Current Scan');
        if checkfields(vw,'curScan')
            val = vw.curScan;
        else
            if checkfields(vw,'ui','scan','sliderHandle')
                val = round(get(vw.ui.scan.sliderHandle,'value'));
            else
                % Sometimes there is no window interface (it is hidden).  Then, we have
                % to find another way to determine the current scan.  Here, we ask the
                % user.  It would be possible to store this information in the VIEW
                % structure.  But we don't.  Ugh.
                val = 1; ieReadNumber('Enter scan number');
            end
        end

    case 'curslice'
        % Return the current slice number. This is the actual slice number
        % if we are in the Inplane view. It is the plane number (sag, cor,
        % or axi) of the currently selected plane if we are in the Volume
        % view. And it is 1 or 2 in the Flat view (for left or right).
        %   curslice = viewGet(vw, 'Current Slice');
        if isequal(vw.name,'hidden')
            % no UI or slider -- use tSeries slice
            curSlice = vw.tSeriesSlice;
            if isnan(curSlice) | isempty(curSlice), val = 1; end
            return
        end
        switch vw.viewType
            case 'Inplane'
                val = vw.tSeriesSlice; % err on the side of not needing a UI
                if isnan(val) && checkfields(vw, 'ui', 'slice')
                    val = get(vw.ui.slice.sliderHandle,'val');
                end
            case {'Volume','Gray','generalGray'}
                sliceOri=getCurSliceOri(vw);
                val=str2double(get(vw.ui.sliceNumFields(sliceOri),'String'));
            case 'Flat'
                if isfield(vw,'numLevels') % test for levels view
                    % flat-levels view (older, but still supported)
                    val = getFlatLevelSlices(vw);
                    val = val(1);
                else
                    % regular flat view: slice is hemisphere, slice 3 means both
                    val = findSelectedButton(vw.ui.sliceButtons);
                end
        end
    case 'nscans'
        % Return the number of scans in the currently selected dataTYPE
        %   nscans = viewGet(vw, 'Number of Scans');
        if length(varargin) < 1, dataType = vw.curDataType;
        else dataType = varargin{1}; end
        if ischar(dataType)
            dataType = existDataType(dataType);
        end
        if dataType==0
            error('Invalid data type specified: %i', dataType);
        end
        val = length(dataTYPES(dataType).scanParams); %TODO: Use dtGet
    case 'nslices'
        % Return the number of slices in the current view struct
        %   nslices = viewGet(vw, 'Number of Slices');
        switch vw.viewType
            case 'Inplane'
                if ~checkfields(vw, 'anat'), vw = loadAnat(vw); end
                val = niftiGet(viewGet(vw,'Anatomy Nifti'),'num slices');
            case {'Volume' 'Gray'}
                val = 1;
            case 'Flat'
                if isfield(vw,'numLevels') % acr levels view
                    val = 2 + sum(vw.numLevels);
                else
                    val = 2;
                end
        end
    case 'montageslices'
        % Return the current subset of slices that are visible in the GUI.
        %   montageSlices = viewGet(vw, 'Montage Slices');
        %
        % only for some view types: inplane montage, flat level
        % Well, this seems not to be true since it returns a value in the
        % Volume view.
        val = viewGet(vw, 'Current Slice');
        if viewGet(vw, 'ishidden') || ~ismember(vw.refreshFn,...
                {'refreshMontageView' 'refreshFlatLevelView'})
            warning('vista:viewError', 'MontageSlices View does not apply to this view type');
            return  % return current slice
        end
        nSlices = get(vw.ui.montageSize.sliderHandle, 'Value');
        val = val:val+nSlices-1;
        val = val(val <= viewGet(vw, 'numSlices'));

    case 'dtname'
        % Return the name of the currently selected dataTYPE
        %   dtName = viewGet(vw, 'Data TYPE Name', [dtnum]);
        if ~isempty(varargin), dtnum = varargin{1};
        else dtnum =  viewGet(vw, 'curdt') ; end
        val = dataTYPES(dtnum).name;
    case 'curdt'
        % Return the number of the currently selected dataTYPE
        %   dtNum = viewGet(vw, 'Current Data TYPE');
        if isfield(vw, 'curDataType')
            val = vw.curDataType;
        else
            val = 0;
        end
    case 'dtstruct'
        % Return the currently selected dataTYPE struct
        %   dtStruct = viewGet(vw, 'DT struct');
        curdt = viewGet(vw, 'Current Data TYPE');
        val   = dataTYPES(curdt);

    case 'size'

        switch viewGet(vw,'View Type')
            case 'Inplane'
                val = viewGet(vw,'anatsize');
            case {'Volume','Gray','generalGray'}
                if isfield(vw, 'anat')
                    % In the GUI this field is populated so we get it
                    if ~isempty(vw.anat), val = viewGet(vw,'Anat Size'); end
                end
                % If it is a hidden volume view then this field is not
                % populated so we need to load the anatomy image in order
                % to calculate it's size
                if ~exist('val','var') || isempty(val)
                    pth = getVAnatomyPath; % assigns it if it's not set
                    % Here we load the volume anatomy - Note that we need
                    % to apply this transform between the dimensions of the
                    % raw image and the mrvista coordinates
                    val = size(readVolAnat(pth));

                    % This was the original call but we note that the
                    % function nifti2mrVistaAnat applies a transform to the
                    % data that makes the data dimensions different from
                    % those in the header.
                    % [~, val] = readVolAnatHeader(pth);
                end
            case 'Flat'
                val = [vw.ui.imSize,2];
        end

    otherwise
        error('Unknown viewGet parameter');

end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%GetNextItem - Get next item from a Batch object (iterator mode).
%
% This is a helper class to easily read and parse batch files.
%
%  USAGE
%
%    [b,item] = GetNextItem(b);
%
%    b              batch object

% Copyright (C) 2007 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

function [b,item] = GetNextItem(b)

% Check number of parameters
if nargin ~= 1,
	error('Incorrect number of parameters (type ''help GetNextItem'' for details).');
end

b.currentItem = b.currentItem+1;
b.currentField = 0;
if b.currentItem > size(b.field,1),
	item = [];
else
	item = b.currentItem;
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%% clear the workspace and select data
% clear; clc; close all;

%% choose data
neuron = Sources2D();
nam = get_fullname('./data_1p.tif');          % this demo data is very small, here we just use it as an example
nam = neuron.select_data(nam);  %if nam is [], then select data interactively

%% parameters
% -------------------------    COMPUTATION    -------------------------  %
pars_envs = struct('memory_size_to_use', 8, ...   % GB, memory space you allow to use in MATLAB
    'memory_size_per_patch', 0.6, ...   % GB, space for loading data within one patch
    'patch_dims', [64, 64]);  %GB, patch size

% -------------------------      SPATIAL      -------------------------  %
gSig = 3;           % pixel, gaussian width of a gaussian kernel for filtering the data. 0 means no filtering
gSiz = 13;          % pixel, neuron diameter
ssub = 1;           % spatial downsampling factor
with_dendrites = false;   % with dendrites or not
if with_dendrites
    % determine the search locations by dilating the current neuron shapes
    updateA_search_method = 'dilate';  %#ok<UNRCH>
    updateA_bSiz = 5;
    updateA_dist = neuron.options.dist;
else
    % determine the search locations by selecting a round area
    updateA_search_method = 'ellipse'; %#ok<UNRCH>
    updateA_dist = 5;
    updateA_bSiz = neuron.options.dist;
end
spatial_constraints = struct('connected', true, 'circular', false);  % you can include following constraints: 'circular'
spatial_algorithm = 'hals_thresh';

% -------------------------      TEMPORAL     -------------------------  %
Fs = 10;             % frame rate
tsub = 1;           % temporal downsampling factor
deconv_flag = true;     % run deconvolution or not
deconv_options = struct('type', 'ar1', ... % model of the calcium traces. {'ar1', 'ar2'}
    'method', 'foopsi', ... % method for running deconvolution {'foopsi', 'constrained', 'thresholded'}
    'smin', -5, ...         % minimum spike size. When the value is negative, the actual threshold is abs(smin)*noise level
    'optimize_pars', true, ...  % optimize AR coefficients
    'optimize_b', true, ...% optimize the baseline);
    'max_tau', 100);    % maximum decay time (unit: frame);

nk = 3;             % detrending the slow fluctuation. usually 1 is fine (no detrending)
% when changed, try some integers smaller than total_frame/(Fs*30)
detrend_method = 'spline';  % compute the local minimum as an estimation of trend.

% -------------------------     BACKGROUND    -------------------------  %
bg_model = 'ring';  % model of the background {'ring', 'svd'(default), 'nmf'}
nb = 1;             % number of background sources for each patch (only be used in SVD and NMF model)
ring_radius = 18;  % when the ring model used, it is the radius of the ring used in the background model.
%otherwise, it's just the width of the overlapping area
num_neighbors = []; % number of neighbors for each neuron
bg_ssub = 2;        % downsample background for a faster speed

% -------------------------      MERGING      -------------------------  %
show_merge = false;  % if true, manually verify the merging step
merge_thr = 0.65;     % thresholds for merging neurons; [spatial overlap ratio, temporal correlation of calcium traces, spike correlation]
method_dist = 'max';   % method for computing neuron distances {'mean', 'max'}
dmin = 5;       % minimum distances between two neurons. it is used together with merge_thr
dmin_only = 2;  % merge neurons if their distances are smaller than dmin_only.
merge_thr_spatial = [0.8, 0.4, -inf];  % merge components with highly correlated spatial shapes (corr=0.8) and small temporal correlations (corr=0.1)

% -------------------------  INITIALIZATION   -------------------------  %
K = [];             % maximum number of neurons per patch. when K=[], take as many as possible.
min_corr = 0.8;     % minimum local correlation for a seeding pixel
min_pnr = 8;       % minimum peak-to-noise ratio for a seeding pixel
min_pixel = gSig^2;      % minimum number of nonzero pixels for each neuron
bd = 0;             % number of rows/columns to be ignored in the boundary (mainly for motion corrected data)
frame_range = [];   % when [], uses all frames
save_initialization = false;    % save the initialization procedure as a video.
use_parallel = true;    % use parallel computation for parallel computing
show_init = true;   % show initialization results
choose_params = true; % manually choose parameters
center_psf = true;  % set the value as true when the background fluctuation is large (usually 1p data)
% set the value as false when the background fluctuation is small (2p)

% -------------------------  Residual   -------------------------  %
min_corr_res = 0.7;
min_pnr_res = 6;
seed_method_res = 'auto';  % method for initializing neurons from the residual
update_sn = true;

% ----------------------  WITH MANUAL INTERVENTION  --------------------  %
with_manual_intervention = false;

% -------------------------  FINAL RESULTS   -------------------------  %
save_demixed = true;    % save the demixed file or not
kt = 3;                 % frame intervals

% -------------------------    UPDATE ALL    -------------------------  %
neuron.updateParams('gSig', gSig, ...       % -------- spatial --------
    'gSiz', gSiz, ...
    'ring_radius', ring_radius, ...
    'ssub', ssub, ...
    'search_method', updateA_search_method, ...
    'bSiz', updateA_bSiz, ...
    'dist', updateA_bSiz, ...
    'spatial_constraints', spatial_constraints, ...
    'spatial_algorithm', spatial_algorithm, ...
    'tsub', tsub, ...                       % -------- temporal --------
    'deconv_flag', deconv_flag, ...
    'deconv_options', deconv_options, ...
    'nk', nk, ...
    'detrend_method', detrend_method, ...
    'background_model', bg_model, ...       % -------- background --------
    'nb', nb, ...
    'ring_radius', ring_radius, ...
    'num_neighbors', num_neighbors, ...
    'bg_ssub', bg_ssub, ...
    'merge_thr', merge_thr, ...             % -------- merging ---------
    'dmin', dmin, ...
    'method_dist', method_dist, ...
    'min_corr', min_corr, ...               % ----- initialization -----
    'min_pnr', min_pnr, ...
    'min_pixel', min_pixel, ...
    'bd', bd, ...
    'center_psf', center_psf);
neuron.Fs = Fs;

%% distribute data and be ready to run source extraction
neuron.getReady(pars_envs);

%% initialize neurons from the video data within a selected temporal range
if choose_params
    % change parameters for optimized initialization
    [gSig, gSiz, ring_radius, min_corr, min_pnr] = neuron.set_parameters();
end

[center, Cn, PNR] = neuron.initComponents_parallel(K, frame_range, save_initialization, use_parallel);
neuron.compactSpatial();
if show_init
    figure();
    ax_init= axes();
    imagesc(Cn, [0, 1]); colormap gray;
    hold on;
    plot(center(:, 2), center(:, 1), '.r', 'markersize', 10);
end

%% estimate the background components
neuron.update_background_parallel(use_parallel);
neuron_init = neuron.copy();

%%  merge neurons and update spatial/temporal components
neuron.merge_neurons_dist_corr(show_merge);
neuron.merge_high_corr(show_merge, merge_thr_spatial);

%% update spatial components

%% pick neurons from the residual
[center_res, Cn_res, PNR_res] =neuron.initComponents_residual_parallel([], save_initialization, use_parallel, min_corr_res, min_pnr_res, seed_method_res);
if show_init
    axes(ax_init);
    plot(center_res(:, 2), center_res(:, 1), '.g', 'markersize', 10);
end
neuron_init_res = neuron.copy();

%% udpate spatial&temporal components, delete false positives and merge neurons
% update spatial
if update_sn
    neuron.update_spatial_parallel(use_parallel, true);
    udpate_sn = false;
else
    neuron.update_spatial_parallel(use_parallel);
end
% merge neurons based on correlations
neuron.merge_high_corr(show_merge, merge_thr_spatial);

for m=1:2
    % update temporal
    neuron.update_temporal_parallel(use_parallel);

    % delete bad neurons
    neuron.remove_false_positives();

    % merge neurons based on temporal correlation + distances
    neuron.merge_neurons_dist_corr(show_merge);
end

%% add a manual intervention and run the whole procedure for a second time
neuron.options.spatial_algorithm = 'nnls';
if with_manual_intervention
    show_merge = true;
    neuron.orderROIs('snr');   % order neurons in different ways {'snr', 'decay_time', 'mean', 'circularity'}
    neuron.viewNeurons([], neuron.C_raw);

    % merge closeby neurons
    neuron.merge_close_neighbors(true, dmin_only);

    % delete neurons
    tags = neuron.tag_neurons_parallel();  % find neurons with fewer nonzero pixels than min_pixel and silent calcium transients
    ids = find(tags>0);
    if ~isempty(ids)
        neuron.viewNeurons(ids, neuron.C_raw);
    end
end
%% run more iterations
neuron.update_background_parallel(use_parallel);
neuron.update_spatial_parallel(use_parallel);
neuron.update_temporal_parallel(use_parallel);

K = size(neuron.A,2);
tags = neuron.tag_neurons_parallel();  % find neurons with fewer nonzero pixels than min_pixel and silent calcium transients
neuron.remove_false_positives();
neuron.merge_neurons_dist_corr(show_merge);
neuron.merge_high_corr(show_merge, merge_thr_spatial);

if K~=size(neuron.A,2)
    neuron.update_spatial_parallel(use_parallel);
    neuron.update_temporal_parallel(use_parallel);
    neuron.remove_false_positives();
end

%% save the workspace for future analysis
neuron.orderROIs('snr');
cnmfe_path = neuron.save_workspace();

%% show neuron contours
Coor = neuron.show_contours(0.6);

%% create a video for displaying the
amp_ac = 140;
range_ac = 5+[0, amp_ac];
multi_factor = 10;
range_Y = 1300+[0, amp_ac*multi_factor];

avi_filename = neuron.show_demixed_video(save_demixed, kt, [], amp_ac, range_ac, range_Y, multi_factor);

%% save neurons shapes
neuron.save_neurons();

author: zhoupc


text:
function gX = lfmKernGradX(kern, x, x2)

% LFMKERNGRADX Gradient of LFM kernel with respect to a point x.
% FORMAT
% DESC computes the gradient of the single input motif
% kernel with respect to the input positions.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x : locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData x numInputs x numData. Where numData is
% the number of data points and numInputs is the number of input
% dimensions in X.
%
% FORMAT
% DESC computes the gradident of the single input motif
% kernel with respect to the input positions where both the row
% positions and column positions are provided separately.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x1 : row locations against which gradients are being computed.
% ARG x2 : column locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData2 x numInputs x numData1. Where numData1 is
% the number of data points in X1, numData2 is the number of data
% points in X2 and numInputs is the number of input
% dimensions in X.
%
% SEEALSO lfmKernParamInit, kernGradX, lfmKernDiagGradX
%
% COPYRIGHT : Neil D. Lawrence, 2006

% KERN

error('lfmKernGradX not yet implemented.')

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% Init Script
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     08/11/2016  Royi Avital
%   *   First release.
%

%% General Parameters
close('all');
clear();
clc();

% set(0, 'DefaultFigureWindowStyle', 'docked');
% defaultLoosInset = get(0, 'DefaultAxesLooseInset');
% set(0, 'DefaultAxesLooseInset', [0.05, 0.05, 0.05, 0.05]);

figPosSmall     = [100, 100, 400, 300];
figPosMedium    = [100, 100, 800, 600];
figPosLarge     = [100, 100, 960, 720];
figPosDefault   = [100, 100, 560, 420];

fontSizeTitle   = 14;
fontSizeAxis    = 12;
fontSizeString  = 12;

lineWidthThin   = 1;
lineWidthNormal = 3;
lineWidthThick  = 4;

markerSizeSmall     = 4;
markerSizeNormal    = 8;
markerSizeLarge     = 10;

% https://www.mathworks.com/help/matlab/graphics_transition/why-are-plot-lines-different-colors.html
% https://www.mathworks.com/matlabcentral/answers/160332
mColorOrder = get(groot, 'DefaultAxesColorOrder');

randomNumberStream  = RandStream('mlfg6331_64', 'NormalTransform', 'Ziggurat');
subStreamNumber     = round(sum(clock()));
subStreamNumber     = 1620;
set(randomNumberStream, 'Substream', subStreamNumber);
RandStream.setGlobalStream(randomNumberStream);

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

author: RoyiAvital

repo: StackExchangeCodes


text:
% [INPUT]
% data = A float t-by-2 matrix (-Inf,Inf) representing the model input.
% type = A string representing the type of metric to calculate:
%   - 'GG' for Gonzalo-Granger Component Metric;
%   - 'H' for Hasbrouck Information Metric.
% lag_max = An integer [2,t-2] representing the maximum lag order to be evaluated (optional, default=10).
% lag_sel = A string representing the lag order selection criteria (optional, default='AIC'):
%   - 'AIC' for Akaike's Information Criterion;
%   - 'BIC' for Bayesian Information Criterion;
%   - 'FPE' for Final Prediction Error;
%   - 'HQIC' for Hannan-Quinn Information Criterion.
%
% [OUTPUT]
% m1 = A float [0,1] representing the first value of the metric.
% m2 = A float [0,1] representing the second value of the metric.
% lag = An integer [1,lag_max] representing the selected lag order.

function [m1,m2,lag] = price_discovery(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('data',@(x)validateattributes(x,{'double'},{'real' 'finite' '2d' 'nonempty' 'size' [NaN 2]}));
        ip.addRequired('type',@(x)any(validatestring(x,{'GG' 'H'})));
        ip.addOptional('lag_max',10,@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 1 'scalar'}));
        ip.addOptional('lag_sel','AIC',@(x)any(validatestring(x,{'AIC' 'BIC' 'FPE' 'HQIC'})));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    [data,lag_max] = validate_input(ipr.data,ipr.lag_max);
    type = ipr.type;
    lag_sel = ipr.lag_sel;

    nargoutchk(2,3);

    [m1,m2,lag] = price_discovery_internal(data,type,lag_max,lag_sel);

end

function [m1,m2,lag] = price_discovery_internal(data,type,lag_max,lag_sel)

    t = size(data,1);
    b = t - 2;

    if (lag_max == 1)
        lag = 1;
    else
        lag = select_lag_order(data,lag_max,lag_sel);
    end

    lag_seq = 1:(lag + 1);

    [~,~,e] = regress(data(:,2),[ones(size(data,1),1) data(:,1)]);
    ect = e((lag + 2):end);

    x1_diff = diff(data(:,1));
    x1 = cell2mat(arrayfun(@(i)x1_diff(i:(b - lag + i)),lag_seq,'UniformOutput',false));

    x2_diff = diff(data(:,2));
    x2 = cell2mat(arrayfun(@(i)x2_diff(i:(b - lag + i)),lag_seq,'UniformOutput',false));

    [b1,~,e1] = regress(x1(:,1),[ones(size(ect,1),1) ect x1(:,2:end-1) x2(:,2:end-1)]);
    [b2,~,e2] = regress(x2(:,1),[ones(size(ect,1),1) ect x1(:,2:end-1) x2(:,2:end-1)]);

    d1 = b1(2);
    d2 = b2(2);

    if (strcmp(type,'GG'))
        v = abs(d1) / (abs(d1) + abs(d2));
    else
        s1 = std(e1);
        s2 = std(e2);
        rho = corr(e1,e2);

        vp = ((-d2 * s1) + (d1 * rho * s2)) ^ 2;
        v = vp / (vp + ((d1 * s2 * sqrt(1 - rho^2)) ^ 2));
    end

    m1 = min(max(v,0),1);
    m2 = 1 - m1;

end

function lag = select_lag_order(data,lag_max,lag_sel)

    t = size(data,1);
    lag_sam = lag_max + 1;
    k = t - lag_sam + 1;

    data_lag = zeros(k,lag_sam * 2);

    for i = 1:2
        sv = ((lag_sam - 1) * 2) + i;

        data_i = data(:,i);
        offs = repmat(1:k,1,lag_sam) + repelem(lag_sam:-1:1,1,k) - 1;

        data_lag(:,i:2:sv) = reshape(data_i(offs),[k lag_sam]);
    end

    data_lag = data_lag(:,3:end);
    s = size(data_lag,1);

    data_y = data(lag_sam:end,:);

    ni = 2:2:(2 * lag_max);
    rhs = [ones(k,1) (lag_sam:k+lag_sam-1).'];
    crit = zeros(lag_max,1);

    for i = 1:lag_max
        data_x = [data_lag(:,1:ni(i)) rhs];

        r = zeros(k,2);

        parfor j = 1:2
            [~,~,e] = regress(data_y(:,j),data_x);
            r(:,j) = e;
        end

        cp = zeros(2);

        for cp_i = 1:2
            for cp_j = 1:2
                cp(cp_i,cp_j) = r(:,cp_i).' * r(:,cp_j);
            end
        end

        sigmad = det(cp / k);
        d = (i * 4) + 4;

        switch (lag_sel)
            case 'AIC'
                crit(i) = log(sigmad) + ((2 / s) * d);
            case 'BIC'
                crit(i) = log(sigmad) + ((log(s) / s) * d);
            case 'FPE'
                ns = size(data_x,2);
                crit(i) = ((s + ns) / (s - ns))^2 * sigmad;
            otherwise
                crit(i) = log(sigmad) + (2 * (log(log(s)) / s) * d);
        end
    end

    [~,lag] = min(crit);

end

function [data,lag_max] = validate_input(data,lag_max)

    t = size(data,1);
    b = t - 2;

    if (t < 5)
        error('The value of ''data'' is invalid. Expected input to be a matrix with at least 5 rows.');
    end

    if (lag_max > b)
        error(['The value of ''lag_max'' is invalid. Expected input to be less than or equal to ' num2str(b) '.']);
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
%cropRectanglesMex crops multiple bounding boxes from the initial image and resizes them to the standard output size.
% The operation is performed on a GPU using NVIDIA Performance Primitives (NPP) library
% cropRectanglesMex was created to prepare batches for training CNNs using MatConvNet (http://www.vlfeat.org/matconvnet/).
%
% Usage:
% crops = cropRectanglesMex( im, boundingBoxes, outputSize);
%
% Inputs:
% im  - the image to crop from, should be a 3 channel image (dimension order: height, width, channels) of type single.
%        Normalization (e.g. [0,1] or [0, 255]) is not important. The image should be stored in RAM (not GPU).
% boundingBoxes - bounding boxes to crop, double[ numBoundingBoxes x 4 ], each line corresponds to one bounding box.
%       The bounding box format is y1, x1, y2, x2, where the origin is in the top-left corner.
%       Pixels are indexed starting from 1 (e.g. [1 1 2 2] corresponds to the box containing the 4 top-left pixels of the image).
%       Bounding boxes can be partially outside of the image. The default value for filling such areas is 0 in all the channels.
% outputSize - the target size of the resized crops, double[2 x 1]. outputSize(1) - the height, outputSize(2) - the width.
%
% Outputs:
% crops - the cropped and resized patches, gpuArray, single[ outputSize(1), outputSize(2), numChannels = 3, numBoundingBoxes ]
%
% The function can be compiled using build_cropResizeMex.m.
% example_cropRectanglesMex.m provides the example of usage

% Anton Osokin, firstname.lastname@gmail.com, March 2015

author: aosokin


text:
% ======================================================== %
% Files of the Matlab programs included in the book:       %
% Xin-She Yang, Nature-Inspired Metaheuristic Algorithms,  %
% Second Edition, Luniver Press, (2010).   www.luniver.com %
% ======================================================== %

% -------------------------------------------------------- %
% Bat-inspired algorithm for continuous optimization (demo)%
% Programmed by Xin-She Yang @Cambridge University 2010    %
% -------------------------------------------------------- %
% Usage: bat_algorithm([20 1000 0.5 0.5]);                 %

% -------------------------------------------------------------------
% This is a simple demo version only implemented the basic          %
% idea of the bat algorithm without fine-tuning the parameters,     %
% Then, though this demo works very well, it is expected that       %
% this demo is much less efficient than the work reported in        %
% the following papers:                                             %
% (Citation details):                                               %
% 1) Yang X.-S., A new metaheuristic bat-inspired algorithm,        %
%    in: Nature Inspired Cooperative Strategies for Optimization    %
%    (NISCO 2010) (Eds. J. R. Gonzalez et al.), Studies in          %
%    Computational Intelligence, Springer, vol. 284, 65-74 (2010).  %
% 2) Yang X.-S., Nature-Inspired Metaheuristic Algorithms,          %
%    Second Edition, Luniver Presss, Frome, UK. (2010).             %
% 3) Yang X.-S. and Gandomi A. H., Bat algorithm: A novel           %
%    approach for global engineering optimization,                  %
%    Engineering Computations, Vol. 29, No. 5, pp. 464-483 (2012).  %
% -------------------------------------------------------------------

% Main programs starts here
function [best,fmin,N_iter]=bat_algorithm(para)
% Display help
 help bat_algorithm.m

% Default parameters
if nargin<1,  para=[20 1000 0.5 0.5];  end
n=para(1);      % Population size, typically 10 to 40
N_gen=para(2);  % Number of generations
A=para(3);      % Loudness  (constant or decreasing)
r=para(4);      % Pulse rate (constant or decreasing)
% This frequency range determines the scalings
% You should change these values if necessary
Qmin=0;         % Frequency minimum
Qmax=2;         % Frequency maximum
% Iteration parameters
N_iter=0;       % Total number of function evaluations
% Dimension of the search variables
d=10;           % Number of dimensions
% Lower limit/bounds/ a vector
Lb=-2*ones(1,d);
% Upper limit/bounds/ a vector
Ub=2*ones(1,d);
% Initializing arrays
Q=zeros(n,1);   % Frequency
v=zeros(n,d);   % Velocities
% Initialize the population/solutions
for i=1:n,
  Sol(i,:)=Lb+(Ub-Lb).*rand(1,d);
  Fitness(i)=Fun(Sol(i,:));
end
% Find the initial best solution
[fmin,I]=min(Fitness);
best=Sol(I,:);

% ======================================================  %
% Note: As this is a demo, here we did not implement the  %
% reduction of loudness and increase of emission rates.   %
% Interested readers can do some parametric studies       %
% and also implementation various changes of A and r etc  %
% ======================================================  %

% Start the iterations -- Bat Algorithm (essential part)  %
for t=1:N_gen,
% Loop over all bats/solutions
        for i=1:n,
          Q(i)=Qmin+(Qmin-Qmax)*rand;
          v(i,:)=v(i,:)+(Sol(i,:)-best)*Q(i);
          S(i,:)=Sol(i,:)+v(i,:);
          % Apply simple bounds/limits
          Sol(i,:)=simplebounds(Sol(i,:),Lb,Ub);
          % Pulse rate
          if rand>r
          % The factor 0.001 limits the step sizes of random walks
              S(i,:)=best+0.001*randn(1,d);
          end

     % Evaluate new solutions
           Fnew=Fun(S(i,:));
     % Update if the solution improves, or not too loud
           if (Fnew<=Fitness(i)) & (rand<A) ,
                Sol(i,:)=S(i,:);
                Fitness(i)=Fnew;
           end

          % Update the current best solution
          if Fnew<=fmin,
                best=S(i,:);
                fmin=Fnew;
          end
        end
        N_iter=N_iter+n;
end
% Output/display
disp(['Number of evaluations: ',num2str(N_iter)]);
disp(['Best =',num2str(best),' fmin=',num2str(fmin)]);

% Application of simple limits/bounds
function s=simplebounds(s,Lb,Ub)
  % Apply the lower bound vector
  ns_tmp=s;
  I=ns_tmp<Lb;
  ns_tmp(I)=Lb(I);

  % Apply the upper bound vector
  J=ns_tmp>Ub;
  ns_tmp(J)=Ub(J);
  % Update this new move
  s=ns_tmp;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Objective function: your own objective function can be written here
% Note: When you use your own function, please remember to
%       change limits/bounds Lb and Ub (see lines 52 to 55)
%       and the number of dimension d (see line 51).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z=Fun(u)
% Sphere function with fmin=0 at (0,0,...,0)
z=sum(u.^2);
%%%%% ============ end ====================================

author: Sable

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal least squares.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/ols_phase_transition_gaussian_dict_gaussian_data.mat';
N = 1024;
pta = spx.pursuit.PhaseTransitionAnalysis(N);
% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.single.OrthogonalLeastSquares(Phi, K).solve(y).z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);

author: indigits

repo: sparse-plex


text:
function [MHz] = kHz2MHz(kHz)
% Convert frequency from kilohertz to megahertz.
% Chad A. Greene 2012
MHz = kHz*1e-3;


text:
function op = proj_linfl2( q )

%PROJ_LINFL2   Projection of each row onto the scaled l2 norm ball.
%    OP = PROJ_LINFL2( Q ) returns an operator implementing the
%    indicator function for the set of l2 norm ball of size q,
%    { X | for all rows i, norm( X(i,:),2) <= q }. Q is optional; if omitted,
%    Q=1 is assumed. But if Q is supplied, it must be a positive
%    real scalar.
% Dual: prox_l1l2.m
% See also: prox_l1, prox_linf, proj_l1

if nargin == 0,
	q = 1;
elseif ~isnumeric( q ) || ~isreal( q ) || numel( q ) ~= 1 || q <= 0,
	error( 'Argument must be positive.' );
end

% In r2007a and later, we can use bsxfun instead of the spdiags trick
if exist('OCTAVE_VERSION','builtin')
    vr = '2000';
else
    vr=version('-release');
end
if str2num(vr(1:4)) >= 2007
    op = @(varargin)proj_linfl2_q_bsxfun( q, varargin{:} );
else
    % the default, using spdiags
    op = @(varargin)proj_linfl2_q( q, varargin{:} );
end

function [ v, x ] = proj_linfl2_q( q, x, t )
v = 0;
switch nargin,
	case 2,
		if nargout == 2,
			error( 'This function is not differentiable.' );
		elseif norm( x(:), Inf ) > q,
			v = Inf;
		end
	case 3,
        % Compute the norms of the rows
        m = size(x,1);
        nrms = sqrt( sum( abs(x).^2 , 2 ) );
        % Scale the rows using left diagonal multiplication
        x = spdiags( min(1,q./nrms), 0, m, m )*x;
	otherwise,
		error( 'Not enough arguments.' );
end

function [ v, x ] = proj_linfl2_q_bsxfun( q, x, t )
v = 0;
switch nargin,
	case 2,
		if nargout == 2,
			error( 'This function is not differentiable.' );
		elseif norm( x(:), Inf ) > q,
			v = Inf;
		end
	case 3,
        nrms = sqrt( sum( abs(x).^2 , 2 ) );
        bsxfun( @times, x, min(1,q./nrms) );
	otherwise,
		error( 'Not enough arguments.' );
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
% Local Regression and Likelihood, Figure 7.4.
%
% Censored Local Likelihood
%
% Author: Catherine Loader
%
% NEED: this is given problems for unknown reasons. something
% to do with censored observations and the ibeta routine.
% if i replace this by
% res[ZLIK] = th*y-y*log(p)
% res[ZDLL] = y*p
% res[ZDDLL] = y*p*(1-p)
% works fine. Memory leak somewhere??
% also works if i ignore censoring.

load border;
fit = locfit(day,runs,'cens',no,'family','geom','alpha',0.7);
figure('Name','fig7_4: censored local likelihood;' );
lfplot(fit);
xlabel('Date');
ylabel('Runs');


text:
function change = delta(ct, mt, y, b)

% delta 函数是非均匀突变使用的非均匀分布。
% 此函数根据当前发电量、最大发电量和可能的偏差量返回变化。
%
% ct - current generation
% mt - maximum generation
% y  - maximum amount of change, i.e. distance from parameter value to bounds
% b  - shape parameter

%%
r = ct / mt;
if(r > 1)
  r = 0.99;
end
change = y * (rand * (1 - r)) ^ b;

repo: Matlab-Machine

path: github-repos/MATLAB/Time9Y-Matlab-Machine/Matlab-Machine-6590d27c49bb2c9415613f132fc2b50fea982c60/009_Optimize time series prediction of BP neural network based on genetic algorithm/goat/delta.m


text:
%plugin_movepath()- Given a path to a plugin folder, this function will
%                     put the plugin at the bottom of the path.
%
% Usage:
%         plugin_movepath('x','begin');           % Put plugin 'x' at the top of the path
%         plugin_movepath('x','end');             % Put plugin 'x' at the bottom of the path
%         plugin_movepath('x','begin','warns',1); % Put plugin 'x' at the top of the path and show warnings
%
%  Inputs:
%         foldername   - [string] Plugin name or part of it
%         pluginpos    - {'begin','end'} Position to move the plugin in the path.
%                      To the top ('begin') or to the bottom ('end').
% Optional inputs:
%        warns        -[0,1] Allow isplay [1] or do not display [0] warnings.
%                      Warnings are restored at the end of the process. Default[0]
%  Outputs:
%        oldpath  - Original MATLAB path before entering this function
%        newpath  - MATLAB path after being modified by this function
%
% Author: Ramon Martinez-Cancino and Arnaud Delorme, SCCN, 2016
%
% Copyright (C) 2016  Ramon Martinez-Cancino,INC, SCCN
%
% 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[oldpath, newpath] = plugin_movepath(foldername,pluginpos, varargin)
oldpath = []; newpath = [];

try
    options = varargin;
    if ~isempty( varargin ),
        for i = 1:2:numel(options)
            g.(options{i}) = options{i+1};
        end
    else g = []; end;
catch
    disp('plugin_movepath() error: calling convention {''key'', value, ... } error'); return;
end;

try g.warns;   catch, g.warns  = 0;  end; % NO warnings by default

% Checking entries
if sum(strcmp(pluginpos,{'begin','end'})) ~= 1
    fprintf(2,'plugin_movepath error: Invalid function argument\n');
    return;
end

% Look in the plugin list for the plugin name provided
global PLUGINLIST
hitindx = find(~cellfun(@isempty,strfind(lower({PLUGINLIST.plugin}),lower(foldername))));
if isempty(hitindx)
    fprintf(2,'plugin_movepath error: Unidentified plugin folder\n')
    return;
else
    eeglabfolder = fileparts(which('eeglab.m'));
    pluginfolder = fullfile(eeglabfolder,'plugins',PLUGINLIST(hitindx).foldername );
end

% Backing up old path
if ismatlab
    oldpath = matlabpath;
else
    oldpath = path;
end;

% Retreiving path
comp = computer;
if strcmpi(comp(1:2), 'PC')
    newpathtest = [ pluginfolder ';' ];
else
    newpathtest = [ pluginfolder ':' ];
end;
ind = strfind(oldpath, newpathtest);

% Checking out if the work is already done
if strcmp(pluginpos,'begin') && ind == 1, return; end;
if strcmp(pluginpos,'end')   && strcmp(oldpath(end-length(pluginfolder):end),pluginfolder), return; end;

% Shooting down warnings
if ~g.warns
    tmpwarn = warning;
    warning off;
end

% Remove folder and subfolders from path
rmpath(genpath(pluginfolder));

% Add folder to the path again in the requested position
 if strcmp(pluginpos,'begin')
     addpath(genpath(pluginfolder),'-begin');
     fprintf(1,['EEGLAB warning: to avoid name conflict ' PLUGINLIST(hitindx).foldername ' functions relocated at the end of the path\n']);
 elseif strcmp(pluginpos,'end')
     addpath(genpath(pluginfolder),'-end');
     fprintf(1,['EEGLAB warning: to avoid name conflict ' PLUGINLIST(hitindx).foldername ' functions relocated at the end of the path\n']);
 end

 % Restoring warnings
 if ~g.warns
     warning(tmpwarn);
 end

 newpath = path; % Retreiving new path

 % required here because path not added yet
% to the admin folder
function res = ismatlab;

v = version;
if v(1) > '4'
    res = 1;
else
    res = 0;
end;


text:
function W = mywigner(Ex)
%MYWIGNER: Calculates the Wigner distribution from a column vector
%
%	W  = mywigner(Ex)
%
%	W  = output Wigner distribution
%	Ex = Input electric field (MUST be a column vector
%
%	Notes:
%		W = Int(-inf..inf){E(x+y)E(x-y)exp[2ixy]}
%
%		E(x+y) & E(x-y) are calculated via a FFT (fast Fourier transform) using the
%		shift theorem. The integration is performed via a FFT. Thus it is important
%		for the data to satisfy the sampling theorem:
%		dy = 2*pi/X			X = span of all x-values	dy = y resolution
%		dx = 2*pi/Y			Y = span of all y-values	dx = x resolution
%		The data must be completely contained within the range x(0)..x(N-1) &
%		y(0)..y(N-1) (i.e. the function must fall to zero within this range).
%
%	v1.0
%
%	Currently waiting for update:
%		Remove the fft/ifft by performing this inside the last function calls
%		Allow an arbitrary output resolution
%		Allow an input vector for x (and possibly y).

if (size(Ex, 2)-1)
    error('E(x) must be a column vector');
end

N = length(Ex);														%   Get length of vector
x = ifftshift(((0:N-1)'-N/2)*2*pi/(N-1));							%   Generate linear vector
X = (0:N-1)-N/2;
EX1 = ifft( (fft(Ex)*ones(1,N)).*exp( i*x*X/2 ));					%   +ve shift
EX2 = ifft( (fft(Ex)*ones(1,N)).*exp( -i*x*X/2 ));					%   -ve shift
W = real(fftshift(fft(fftshift(EX1.*conj(EX2), 2), [], 2), 2));		%   Wigner function


text:
function [ MRI_signal_noisy ] = add_Rician_noise_to_simulated_dMRI_data( MRI_signal, SNR )
%ADD_RICIAN_NOISE_TO_SIMULATED_DMRI_DATA adds Rician noise to simulated
%MRI_signals.
%
% inputs: noisefree MRI_signal and SNR
% outputs: noisy MRI_signal
%
% noisy signal, S = |S* + N1 + iN2|
% where S* = noise-free signal
% N1 and N2 are random numbers, drawn from normal distribution, with mean=0, and std = 1/SNR.

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

no_of_directions =  size(MRI_signal);
N1 = (1/SNR).*randn(no_of_directions);
N2 = (1/SNR).*randn(no_of_directions);

MRI_signal_noisy = abs(MRI_signal + N1 + 1i.*N2);
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%% function
clear; clc; close all;
T = 1000;
s = double(rand(1, T)>0.98);
sig = 0.04;

% example
tau_d = 5;
tau_r = 1;
nMax = 100;
pars = [tau_d, tau_r];
kernel = create_kernel('exp2', pars, nMax);

t = 1:kernel.nMax;
gt = kernel.fhandle(kernel.pars, t);
c1 = conv(s, gt);
c1 = c1(1:T);
y1 = c1 + randn(1, T) * sig;

%% use the true convolution kernel
kernel0 = kernel;
kernel0.pars = [10, 0.1];
kernel0.bound_pars = false;
figure('position', [1,1,1500, 200]);
plot(y1);
hold on;
plot(c1, 'r', 'linewidth', 2);
plot(-s*0.1, 'r', 'linewidth', 2);
tic;
[chat, shat, kernel_fit, iters] = deconvCa(y1, kernel0, 2, true, false);
toc;
plot(chat,'-.g','linewidth', 2);
plot(-shat*0.1, '-.g', 'linewidth', 2); %, '-.');
legend('data', 'ground truth: c','ground truth: s', 'OASIS:c', 'OASIS: s'); % 'gound truth: s', 'OASIS: c', 'OASIS: s');

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function a = r8bb_set ( n1, n2, ml, mu, a, i, j, value )

%*****************************************************************************80
%
%% R8BB_SET sets an entry of a R8BB matrix.
%
%  Discussion:
%
%    The R8BB storage format is for a border banded matrix.  Such a
%    matrix has the logical form:
%
%      A1 | A2
%      ---+---
%      A3 | A4
%
%    with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
%    are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
%    respectively.
%
%    A should be defined as a vector.  The user must then store
%    the entries of the four blocks of the matrix into the vector A.
%    Each block is stored by columns.
%
%    A1, the banded portion of the matrix, is stored in
%    the first (2*ML+MU+1)*N1 entries of A, using standard LINPACK
%    general band format.  The reason for the factor of 2 in front of
%    ML is to allocate space that may be required if pivoting occurs.
%
%    The following formulas should be used to determine how to store
%    the entry corresponding to row I and column J in the original matrix:
%
%    Entries of A1:
%
%      1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
%
%      Store the I, J entry into location
%      (I-J+ML+MU+1)+(J-1)*(2*ML+MU+1).
%
%    Entries of A2:
%
%      1 <= I <= N1, N1+1 <= J <= N1+N2.
%
%      Store the I, J entry into location
%      (2*ML+MU+1)*N1+(J-N1-1)*N1+I.
%
%    Entries of A3:
%
%      N1+1 <= I <= N1+N2, 1 <= J <= n1.
%
%      Store the I, J entry into location
%      (2*ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
%
%    Entries of A4:
%
%      N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
%
%      Store the I, J entry into location
%      (2*ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
%      (same formula used for A3).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N1, N2, the order of the banded and dense blocks.
%    N1 and N2 must be nonnegative, and at least one must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative, and no greater than N1-1.
%
%    Input, real A((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the R8BB matrix.
%
%    Input, integer I, J, the row and column of the entry to be set.
%
%    Input, real VALUE, the value to be assigned to the (I,J) entry.
%
%    Output, real A((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the updated R8BB matrix.
%
  if ( i <= 0 | n1+n2 < i )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8BB_SET - Fatal error!\n' );
    fprintf ( 1, 'R8BB_SET - Illegal value of row index I = %d\n', i );
    error ( 'R8BB_SET - Fatal error!' );
  end

  if ( j <= 0 | n1+n2 < j )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8BB_SET - Fatal error!\n' );
    fprintf ( 1, 'R8BB_SET - Illegal value of column index J = %d\n', j );
    error ( 'R8BB_SET - Fatal error!' );
  end
%
%  The A1 block of the matrix.
%
%  Check for out of band problems.
%
%  Normally, we would check the condition MU < (J-I), but the storage
%  format requires extra entries be set aside in case of pivoting, which
%  means that the condition becomes MU+ML < (J-I).
%
  if ( i <= n1 & j <= n1 )
    if ( mu+ml < (j-i) | ml < (i-j) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8BB_SET - Warning!\n' );
      fprintf ( 1, 'R8BB_SET - Unable to set entry A(%d,%d).\n', i, j );
      error ( 'R8BB_SET - Warning!' );
    else
      ij = (i-j+ml+mu+1)+(j-1)*(2*ml+mu+1);
    end
%
%  The A2 block of the matrix.
%
  elseif ( i <= n1 & n1 < j )
    ij = (2*ml+mu+1)*n1+(j-n1-1)*n1+i;
%
%  The A3 and A4 blocks of the matrix.
%
  elseif ( n1 < i )
    ij = (2*ml+mu+1)*n1+n2*n1+(j-1)*n2+(i-n1);
  end

  a(ij) = value;

  return
end


text:
%  Save NIFTI dataset. Support both *.nii and *.hdr/*.img file extension.
%  If file extension is not provided, *.hdr/*.img will be used as default.
%
%  Usage: save_nii(nii, filename, [old_RGB])
%
%  nii.hdr - struct with NIFTI header fields (from load_nii.m or make_nii.m)
%
%  nii.img - 3D (or 4D) matrix of NIFTI data.
%
%  filename - NIFTI file name.
%
%  old_RGB    - an optional boolean variable to handle special RGB data
%       sequence [R1 R2 ... G1 G2 ... B1 B2 ...] that is used only by
%       AnalyzeDirect (Analyze Software). Since both NIfTI and Analyze
%       file format use RGB triple [R1 G1 B1 R2 G2 B2 ...] sequentially
%       for each voxel, this variable is set to FALSE by default. If you
%       would like the saved image only to be opened by AnalyzeDirect
%       Software, set old_RGB to TRUE (or 1). It will be set to 0, if it
%       is default or empty.
%
%  Tip: to change the data type, set nii.hdr.dime.datatype,
%	and nii.hdr.dime.bitpix to:
%
%     0 None                     (Unknown bit per voxel) % DT_NONE, DT_UNKNOWN
%     1 Binary                         (ubit1, bitpix=1) % DT_BINARY
%     2 Unsigned char         (uchar or uint8, bitpix=8) % DT_UINT8, NIFTI_TYPE_UINT8
%     4 Signed short                  (int16, bitpix=16) % DT_INT16, NIFTI_TYPE_INT16
%     8 Signed integer                (int32, bitpix=32) % DT_INT32, NIFTI_TYPE_INT32
%    16 Floating point    (single or float32, bitpix=32) % DT_FLOAT32, NIFTI_TYPE_FLOAT32
%    32 Complex, 2 float32      (Use float32, bitpix=64) % DT_COMPLEX64, NIFTI_TYPE_COMPLEX64
%    64 Double precision  (double or float64, bitpix=64) % DT_FLOAT64, NIFTI_TYPE_FLOAT64
%   128 uint RGB                  (Use uint8, bitpix=24) % DT_RGB24, NIFTI_TYPE_RGB24
%   256 Signed char            (schar or int8, bitpix=8) % DT_INT8, NIFTI_TYPE_INT8
%   511 Single RGB              (Use float32, bitpix=96) % DT_RGB96, NIFTI_TYPE_RGB96
%   512 Unsigned short               (uint16, bitpix=16) % DT_UNINT16, NIFTI_TYPE_UNINT16
%   768 Unsigned integer             (uint32, bitpix=32) % DT_UNINT32, NIFTI_TYPE_UNINT32
%  1024 Signed long long              (int64, bitpix=64) % DT_INT64, NIFTI_TYPE_INT64
%  1280 Unsigned long long           (uint64, bitpix=64) % DT_UINT64, NIFTI_TYPE_UINT64
%  1536 Long double, float128  (Unsupported, bitpix=128) % DT_FLOAT128, NIFTI_TYPE_FLOAT128
%  1792 Complex128, 2 float64  (Use float64, bitpix=128) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
%  2048 Complex256, 2 float128 (Unsupported, bitpix=256) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
%
%  Part of this file is copied and modified from:
%  http://www.mathworks.com/matlabcentral/fileexchange/1878-mri-analyze-tools
%
%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%  - "old_RGB" related codes in "save_nii.m" are added by Mike Harms (2006.06.28)
%
function save_nii(nii, fileprefix, old_RGB)

   if ~exist('nii','var') || isempty(nii) || ~isfield(nii,'hdr') || ...
	~isfield(nii,'img') || ~exist('fileprefix','var') || isempty(fileprefix)

      error('Usage: save_nii(nii, filename, [old_RGB])');
   end

   if isfield(nii,'untouch') && nii.untouch == 1
      error('Usage: please use ''save_untouch_nii.m'' for the untouched structure.');
   end

   if ~exist('old_RGB','var') || isempty(old_RGB)
      old_RGB = 0;
   end

   v = version;

   %  Check file extension. If .gz, unpack it into temp folder
   %
   if length(fileprefix) > 2 && strcmp(fileprefix(end-2:end), '.gz')

      if ~strcmp(fileprefix(end-6:end), '.img.gz') && ...
	 ~strcmp(fileprefix(end-6:end), '.hdr.gz') && ...
	 ~strcmp(fileprefix(end-6:end), '.nii.gz')

         error('Please check filename.');
      end

      if str2num(v(1:3)) < 7.1 || ~usejava('jvm')
         error('Please use MATLAB 7.1 (with java) and above, or run gunzip outside MATLAB.');
      else
         gzFile = 1;
         fileprefix = fileprefix(1:end-3);
      end
   end

   filetype = 1;

   %  Note: fileprefix is actually the filename you want to save
   %
   if ~isempty(findstr('.nii',fileprefix)) && ...
                        strcmp(fileprefix(end-3:end), '.nii')
      filetype = 2;
      fileprefix(end-3:end)='';
   end

   if ~isempty(findstr('.hdr',fileprefix)) && ...
                        strcmp(fileprefix(end-3:end), '.hdr')
      fileprefix(end-3:end)='';
   end

   if ~isempty(findstr('.img',fileprefix)) && ...
                        strcmp(fileprefix(end-3:end), '.img')
      fileprefix(end-3:end)='';
   end

   write_nii(nii, filetype, fileprefix, old_RGB);

   %  gzip output file if requested
   %
   if exist('gzFile', 'var')
      if filetype == 1
         gzip([fileprefix, '.img']);
         delete([fileprefix, '.img']);
         gzip([fileprefix, '.hdr']);
         delete([fileprefix, '.hdr']);
      elseif filetype == 2
         gzip([fileprefix, '.nii']);
         delete([fileprefix, '.nii']);
      end;
   end;

   if filetype == 1

      %  So earlier versions of SPM can also open it with correct originator
      %
      M=[[diag(nii.hdr.dime.pixdim(2:4)) -[nii.hdr.hist.originator(1:3).*nii.hdr.dime.pixdim(2:4)]'];[0 0 0 1]];
      save([fileprefix '.mat'], 'M');
   end

   return					% save_nii

%-----------------------------------------------------------------------------------
function write_nii(nii, filetype, fileprefix, old_RGB)

   hdr = nii.hdr;

   if isfield(nii,'ext') && ~isempty(nii.ext)
      ext = nii.ext;
      [ext, esize_total] = verify_nii_ext(ext);
   else
      ext = [];
   end

   switch double(hdr.dime.datatype),
   case   1,
      hdr.dime.bitpix = int16(1 ); precision = 'ubit1';
   case   2,
      hdr.dime.bitpix = int16(8 ); precision = 'uint8';
   case   4,
      hdr.dime.bitpix = int16(16); precision = 'int16';
   case   8,
      hdr.dime.bitpix = int16(32); precision = 'int32';
   case  16,
      hdr.dime.bitpix = int16(32); precision = 'float32';
   case  32,
      hdr.dime.bitpix = int16(64); precision = 'float32';
   case  64,
      hdr.dime.bitpix = int16(64); precision = 'float64';
   case 128,
      hdr.dime.bitpix = int16(24); precision = 'uint8';
   case 256
      hdr.dime.bitpix = int16(8 ); precision = 'int8';
   case 511,
      hdr.dime.bitpix = int16(96); precision = 'float32';
   case 512
      hdr.dime.bitpix = int16(16); precision = 'uint16';
   case 768
      hdr.dime.bitpix = int16(32); precision = 'uint32';
   case 1024
      hdr.dime.bitpix = int16(64); precision = 'int64';
   case 1280
      hdr.dime.bitpix = int16(64); precision = 'uint64';
   case 1792,
      hdr.dime.bitpix = int16(128); precision = 'float64';
   otherwise
      error('This datatype is not supported');
   end

   hdr.dime.glmax = round(double(max(nii.img(:))));
   hdr.dime.glmin = round(double(min(nii.img(:))));

   if filetype == 2
      fid = fopen(sprintf('%s.nii',fileprefix),'w');

      if fid < 0,
         msg = sprintf('Cannot open file %s.nii.',fileprefix);
         error(msg);
      end

      hdr.dime.vox_offset = 352;

      if ~isempty(ext)
         hdr.dime.vox_offset = hdr.dime.vox_offset + esize_total;
      end

      hdr.hist.magic = 'n+1';
      save_nii_hdr(hdr, fid);

      if ~isempty(ext)
         save_nii_ext(ext, fid);
      end
   else
      fid = fopen(sprintf('%s.hdr',fileprefix),'w');

      if fid < 0,
         msg = sprintf('Cannot open file %s.hdr.',fileprefix);
         error(msg);
      end

      hdr.dime.vox_offset = 0;
      hdr.hist.magic = 'ni1';
      save_nii_hdr(hdr, fid);

      if ~isempty(ext)
         save_nii_ext(ext, fid);
      end

      fclose(fid);
      fid = fopen(sprintf('%s.img',fileprefix),'w');
   end

   ScanDim = double(hdr.dime.dim(5));		% t
   SliceDim = double(hdr.dime.dim(4));		% z
   RowDim   = double(hdr.dime.dim(3));		% y
   PixelDim = double(hdr.dime.dim(2));		% x
   SliceSz  = double(hdr.dime.pixdim(4));
   RowSz    = double(hdr.dime.pixdim(3));
   PixelSz  = double(hdr.dime.pixdim(2));

   x = 1:PixelDim;

   if filetype == 2 && isempty(ext)
      skip_bytes = double(hdr.dime.vox_offset) - 348;
   else
      skip_bytes = 0;
   end

   if double(hdr.dime.datatype) == 128

      %  RGB planes are expected to be in the 4th dimension of nii.img
      %
      if(size(nii.img,4)~=3)
         error(['The NII structure does not appear to have 3 RGB color planes in the 4th dimension']);
      end

      if old_RGB
         nii.img = permute(nii.img, [1 2 4 3 5 6 7 8]);
      else
         nii.img = permute(nii.img, [4 1 2 3 5 6 7 8]);
      end
   end

   if double(hdr.dime.datatype) == 511

      %  RGB planes are expected to be in the 4th dimension of nii.img
      %
      if(size(nii.img,4)~=3)
         error(['The NII structure does not appear to have 3 RGB color planes in the 4th dimension']);
      end

      if old_RGB
         nii.img = permute(nii.img, [1 2 4 3 5 6 7 8]);
      else
         nii.img = permute(nii.img, [4 1 2 3 5 6 7 8]);
      end
   end

   %  For complex float32 or complex float64, voxel values
   %  include [real, imag]
   %
   if hdr.dime.datatype == 32 || hdr.dime.datatype == 1792
      real_img = real(nii.img(:))';
      nii.img = imag(nii.img(:))';
      nii.img = [real_img; nii.img];
   end

   if skip_bytes
      fwrite(fid, zeros(1,skip_bytes), 'uint8');
   end

   fwrite(fid, nii.img, precision);
%   fwrite(fid, nii.img, precision, skip_bytes);        % error using skip
   fclose(fid);

   return;					% write_nii

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function fx = p10_fx ( x )

%*****************************************************************************80
%
%% P10_FX evaluates the Repeller.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(*), the point at which F is to be evaluated.
%
%    Output, real FX(*), the value of the function at X.
%
  fx = 20.0 * x ./ ( 100.0 * x .* x + 1.0 );

  return
end


text:
% test for image rotation

M = load_image('lena', 256);

theta = pi/50;
for i=1:100
    M = perform_image_rotation(M,theta);
    imageplot(M); drawnow;
end

author: gpeyre

repo: matlab-toolboxes