text:
function M = legendreMass(k)
% Compute the mass matrix using legendre polinomials of order 'k' as a base
% for our expansion.

M = diag(2./(2*(0:k)+1));


text:
function roi = roiEdit(roi);
% Edit / View roi properties dialog.
%
%  roi = roiEdit(roi);
%
% ras 09/06/2005: wrote it (based on editROIFields).

%%%%%%%%%%
% PARAMS %
%%%%%%%%%%
% convert color into [R G B] vector if needed
colorList = {'yellow' 'magenta' 'cyan' 'red' 'green' 'blue' ...
            'white' 'kblack' 'user'};
colorChar = char(colorList)';
colorChar = colorChar(1,:);
if ischar(roi.color)
    colorNum = findstr(roi.color,colorChar);
    switch (colorNum)
        case 'yellow', colorRGB=[1 1 0];
        case 'magenta', colorRGB=[1 0 1];
        case 'cyan', colorRGB=[0 1 1];
        case 'red', colorRGB=[1 0 0];
        case 'green', colorRGB=[0 1 0];
        case 'blue', colorRGB=[0 0 1];
        case 'white', colorRGB=[1 1 1];
        case 'k', colorRGB=[0 0 0];
        otherwise
            colorRGB=[0 0 0];
    end
else
    colorNum=8; colorRGB=roi.color;
end

% options for ROI type
typeList = {'volume' 'area' 'line' 'point'};
typeNum = cellfind(typeList, lower(roi.type));

% options for ROI fill mode
fillList = {'perimeter' 'filled' 'patches'};
fillNum = cellfind(fillList, lower(roi.fillMode));

%%%%%%%%%%%%%%%%%
% Create Dialog %
%%%%%%%%%%%%%%%%%
dlg(1).string = 'ROI name:';
dlg(1).fieldName = 'name';
dlg(1).style = 'edit';
dlg(1).value = roi.name;

dlg(end+1).string = 'Color:';
dlg(end).fieldName = 'color';
dlg(end).list = colorList;
dlg(end).choice = colorNum;
dlg(end).style = 'popupmenu';
if ischar(roi.color)
    dlg(end).value = colorList{colorNum};
else
    dlg(end).value = 'user';
end

dlg(end+1).string = 'Size:';
dlg(end).fieldName = 'size';
nVoxels = size(roi.coords, 2);
volume = nVoxels * prod(roi.voxelSize);
dlg(end).value = sprintf('%i voxels, %3.1f mm^3', nVoxels, volume);
dlg(end).style = 'text';

dlg(end+1).string = 'Fill ROI When Drawing?';
dlg(end).fieldName = 'fillMode';
dlg(end).list = fillList;
dlg(end).value = fillNum;
dlg(end).style = 'popupmenu';

dlg(end+1).string = 'Save ROI As:';
dlg(end).fieldName = 'viewType';
dlg(end).list = {'Inplane' 'Volume'};
dlg(end).value = roi.viewType;
dlg(end).style = 'popup';

dlg(end+1).string = 'Created:';
dlg(end).fieldName = 'created';
dlg(end).value = roi.created;
dlg(end).style = 'text';

dlg(end+1).string = 'Last Modified:';
dlg(end).fieldName = 'modified';
dlg(end).list = colorList;
dlg(end).value = roi.modified;
dlg(end).style = 'text';

dlg(end+1).string = 'Comments:';
dlg(end).fieldName = 'comments';
dlg(end).value = roi.comments;
dlg(end).style = 'edit4';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% put up dialog, get response %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
resp = generalDialog(dlg, 'Edit ROI Fields', [.24 .5 .25 .25]);
% resp.color = cellfind(colorList, resp.color);

% If user selects 'OK', change the parameters.  Otherwise the
% user isn't happy with these settings so bail out.
if isempty(resp), return; end

resp = rmfield(resp, 'size');

roi = mergeStructures(roi, resp);
if (strcmp(roi.color,'user'))
    roi.color = uisetcolor(colorRGB);
else
    roi.color = roi.color(1);
end

roi.modified = datestr(now);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function dpixc = makeavg(dpixc_ind, blinkmat, offset, avg)
% dpixc = makeavg(dpixc_ind, blinkmat, offset, avg)
% avg = 0 -> noise free image
if size(blinkmat,2) == 1 %one slice only
    dpixc_nonoise = dpixc_ind'*blinkmat + offset;
else
    dpixc_nonoise = array2im(dpixc_ind'*blinkmat + offset);
end
dpixc_dip = newim(dpixc_nonoise);
for ii=1:avg
    dpixc_dip = dpixc_dip + noise(dpixc_nonoise,'poisson');
end
if avg==0
    dpixc_dip = dpixc_nonoise;
end
navg=max(1,avg);
dpixc = double(dpixc_dip/navg);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [e,E_q] = q2e(q)

% Q2E  Quaternion to Euler angles conversion.
%   Q2E(Q) returns an Euler angles vector [roll;pitch;yaw] corresponding to
%   the orientation quaternion Q.
%
%   [E,Jq] = Q2E(Q) returns also the Jacobian matrix.
%
%   See also QUATERNION, EULERANGLES, R2Q, E2Q, Q2V.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.
%   Revised 2018 Joan Sola - handle full gimbal range (2pi, pi, 2pi)

a  = q(1); % a: real part
b  = q(2);
c  = q(3);
d  = q(4);

r32 =  2*c*d + 2*a*b;
r33 =  a^2   - b^2   - c^2 + d^2;
r31 =  2*b*d - 2*a*c;
r21 =  2*b*c + 2*a*d;
r11 =  a^2   + b^2   - c^2 - d^2;

n   = sqrt(r11^2+r21^2);

e   = [ atan2( r32,r33)
        atan2(-r31, n)
        atan2( r21,r11) ];

if nargout >1

    % partials of R wrt q
    dr11dq  = [  2*a,  2*b, -2*c, -2*d];
    dr21dq  = [  2*d,  2*c,  2*b,  2*a];
    dr31dq  = [ -2*c,  2*d, -2*a,  2*b];
    dr32dq  = [  2*b,  2*a,  2*d,  2*c];
    dr33dq  = [  2*a, -2*b, -2*c,  2*d];

    % partials of n wrt R elements
    dndr11 = r11/n;
    dndr21 = r21/n;

    % partials of e wrt n
    de2dn   = r31/(n^2 + r31^2);

    % partials of e wrt R elements
    de1dr32 =  r33/(r33^2 + r32^2);
    de1dr33 = -r32/(r33^2 + r32^2);
    de2dr31 = -n /(n^2 + r31^2);
    de3dr11 = -r21/(r11^2 + r21^2);
    de3dr21 =  r11/(r11^2 + r21^2);

    % chain rule: rows of Jacobian
    de1dq  = de1dr33*dr33dq + de1dr32*dr32dq;
    de2dq  = de2dr31*dr31dq + de2dn*(dndr11*dr11dq+dndr21*dr21dq);
    de3dq  = de3dr11*dr11dq + de3dr21*dr21dq;

    % assemble Jacobian
    E_q    = [de1dq;de2dq;de3dq];
end

return

%% jac

syms a b c d real
q=[a;b;c;d];
[e,E_q] = q2e(q);
simplify(E_q-jacobian(e,q))

%% test
for i = 1 : 1000
    e = [2/3;1/3;2/3] .* randn(3,1);
    eo = q2e(e2q(e));
    if any(e-eo > 1e-10) % Should never enter this IF
        e
        eo
        disp('###########')
    end
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [f] = spm_mountaincar_fun(P,G)
% [Cross-entropy] objective function for mountain car problem
% FORMAT [f] = spm_mountaincar_fun(P,G)
%
% P = spm_vec(P)
% P.a - 0th order coefficients of force
% P.b - 1st order coefficients of force
% P.c - 2nd order coefficients of force
% P.d - action efficacy
%
% G   - world model; including
%    G.fq : function fq(x) returning desired equilibrium density at x
%    G.X  : matrix of locations in x
%
% f   - KL divergence between actual and desired equilibrium densities
% x   - cell of grid point support
%
% see:
% Gaussian Processes in Reinforcement Learning
% Carl Edward Rasmussen and Malte Kuss
% Max Planck Institute for Biological Cybernetics
% Spemannstraße 38, 72076 T¨ubingen, Germany
% {carl,malte.kuss}@tuebingen.mpg.de
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_mountaincar_fun.m 7679 2019-10-24 15:54:07Z spm $

% place paramters in model
%--------------------------------------------------------------------------
G(1).pE = P;

% get equilibrium density
%--------------------------------------------------------------------------
[M0,q0] = spm_fp(G);

% desired equilibrium density
%--------------------------------------------------------------------------
q   = feval(G(1).fq,G(1).X);

% KL divergence or cross entropy
%--------------------------------------------------------------------------
q   = q/sum(q(:));
D   = q(:)./q0(:);
i   = find(D > exp(-16));
f   = q(i)'*log(D(i));

author: spm


text:
function [noiseVariance, debiasedNoiseVariance, IQEstimate, noiseEstimateOomen] = realized_noise_estimate(price, time, timeType, options)
% Estimation of the optimal bandwidth to use when estimating the quadratic variation using a Realized Kernel
%
% USAGE:
%   [NOISEVARIANCE,DEBIASEDNOISEVARIANCE,IQESIQESTIMATETIMATE,NOISEVARIANCEOOMEN]
%                        = realized_kernel_select_bandwidth(PRICE,TIME,TIMETYPE,OPTIONS)
%
% INPUTS:
%   PRICE       - m by 1 vector of prices
%   TIME        - m by 1 vector of times where TIME(i) corresponds to PRICE(i)
%   TIMETYPE    - String describing the way times are measured
%               'wall'    24-hour clock of the form HHMMSS, e.g. 101543 or 153217
%               'seconds' Time measured in seconds past midnight on the first day
%               'unit'  Unit normalized date format, e.g. .1, .234, .9. Unit normalized times are
%                 more general than the other types and can be applied to data from more than one
%                 calendar day
%   OPTIONS     - Realized kernel options structure.  See help realized_kernel_options
%
% OUTPUTS:
%   NOISEVARIANCE         - Bandi-Russel noise variance estimate
%   DEBIASEDNOISEVARIANCE - Debiased Bandi-Russel noise variance estimate using the adjustment of BNHLS
%   IQESTIMATE            - An estimate of the lower bound of the IQ based on low-freuqency returns
%   NOISEVARIANCEOOMEN    - Estimate using Oomen (2006) alternative AC(1) estimator
%
% COMMENTS:
%   This is a helper function for REALIZED_KERNEL.  See Barndorf-Nielsen, Hansen, Lunde and Shephard
%   (2008) for details about the optimal selection of bandwidth
%
%  See also REALIZED_KERNEL, REALIZED_PRICE_FILTER, REALIZED_KERNEL_WEIGHTS
%  REALIZED_KERNEL_SELECT_LAG_LENGTH, REALIZED_VARIANCE

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin~=4
    error('Seven inputs required.')
end
if size(price,2)>size(price,1)
    price=price';
end
if size(price,2)>1
    error('PRICE must be a m by 1 vector.')
end
if size(time,2)>size(time,1)
    time=time';
end
if any(diff(time)<0)
    error('TIME must be sorted and increasing')
end
if size(time,2)>1 || length(time)~=length(price)
    error('TIME must be a m by 1 vector.')
end
time = double(time);
timeType=lower(timeType);
if ~ismember(timeType,{'wall','seconds','matlab'})
    error('TIMETYPE must be one of ''wall'', ''seconds'' or ''Matlab''.');
end

% List of flat top kernels
flatTopKernelList = {'bartlett','twoscale','2ndorder','epanechnikov',...
    'cubic','multiscale','5thorder','6thorder','7thorder','8thorder','parzen',...
    'th1','th2','th5','th16'};

% List of non flat top kernels
nonFlatTopKernelList = {'nonflatparzen','qs','fejer','thinf','bnhls'};

% Combined kernel list
kernelList = [flatTopKernelList nonFlatTopKernelList];

if ~isfield(options,'medFrequencyKernel') || ~ismember(options.medFrequencyKernel,kernelList)
    error('KERNEL must be a field of OPTIONS and one of the listed types.')
end

medFrequencySamplingType = options.medFrequencySamplingType;
medFrequencySamplingInterval = options.medFrequencySamplingInterval;
medFrequencyKernel = options.medFrequencyKernel;
medFrequencyBandwidth = options.medFrequencyBandwidth;
if ~isscalar(medFrequencySamplingInterval) || floor(medFrequencySamplingInterval)~=medFrequencySamplingInterval || medFrequencySamplingInterval<=0
    error('MEDIUMFREQUENCYTIME must be a postive integer scalar value.');
end

IQEstimationSamplingType = options.IQEstimationSamplingType;
IQEstimationSamplingInterval = options.IQEstimationSamplingInterval;
if ~isscalar(IQEstimationSamplingInterval) || floor(IQEstimationSamplingInterval)~=IQEstimationSamplingInterval || IQEstimationSamplingInterval<=0
    error('LOWFREQUENCYTIME must be a postive integer scalar value.');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% To compute the "optimal" bandwidth an estimate of the noise and the QV using low frequency data
lowFrequencyRealizedVariance = realized_variance(price,time, timeType , IQEstimationSamplingType, IQEstimationSamplingInterval);

noiseVarianceSamplingType = options.noiseVarianceSamplingType;
noiseVarianceSamplingInterval = options.noiseVarianceSamplingInterval;

% Compute the noiseVariance using the debiased version of Bandi-Russell, using all prices
noiseVariance = realized_variance(price, time, timeType , noiseVarianceSamplingType , noiseVarianceSamplingInterval);
% Compute the effective 'n' to use in the Bandi-Russel estimator
noiseFilteredPrice = realized_price_filter(price, time, timeType , noiseVarianceSamplingType , noiseVarianceSamplingInterval);
if options.useAdjustedNoiseCount
    % Only count the number of non-zero returns
    n = sum(diff(noiseFilteredPrice)~=0);
else
    n = length(noiseFilteredPrice) - 1;
end
noiseVariance = noiseVariance/(2*n);

noiseReturns = diff(log(noiseFilteredPrice));
n = length(noiseReturns);
noiseEstimateOomen = -1/(n-1) * noiseReturns(1:end-1)'*noiseReturns(2:end);

% Require a realized kernel estimate to adjust the variance
medFrequencyOptions = realized_options('kernel');
medFrequencyOptions.kernel = medFrequencyKernel;
medFrequencyOptions.bandwidth = medFrequencyBandwidth;
medFrequencyOptions.endTreatment = 'stagger';

medFreuqencyRealizedKernel = realized_kernel(price, time, timeType, medFrequencySamplingType, medFrequencySamplingInterval, medFrequencyOptions);
medFrequencyRealizedVariance = realized_variance(price,time, timeType , medFrequencySamplingType, medFrequencySamplingInterval);
% Use the BNHLS corrected noise estimate
debiasedNoiseVariance = exp( log(noiseVariance) - medFreuqencyRealizedKernel/medFrequencyRealizedVariance);
IQEstimate = lowFrequencyRealizedVariance^2;

author: bashtage

repo: mfe-toolbox


text:
function ax = axis_pct(pct)
% AXIS_PCT       Set reasonable axis limits.
%
% AXIS_PCT(pct) sets axis limits to extend pct% beyond limits of plotted
% objects.  Default is 5%.
% Works for linear or log scale.
% Unfortunately, the axes won't change when new points are plotted.

% Written by Tom Minka

if nargin < 1
  pct = 0.05;
end
ax = [Inf -Inf Inf -Inf Inf -Inf];

% find bounding box of plotted objects
children = get(gca,'children');
for child = children'
  if strcmp(get(child,'type'),'text')
    xyz = get(child,'position');
    % need to determine bounding box of the text
    if strcmp(get(gca,'xscale'), 'log')
      xyz(1) = log10(xyz(1));
    end
    if strcmp(get(gca,'yscale'), 'log')
      xyz(2) = log10(xyz(2));
    end
    if strcmp(get(gca,'zscale'), 'log')
      xyz(3) = log10(xyz(3));
    end
    c([1 2]) = xyz(1);
    c([3 4]) = xyz(2);
    c([5 6]) = xyz(3);
  else
    x = get(child,'xdata');
    x = x(finite(x));
    if strcmp(get(gca,'xscale'), 'log')
      x = x(x > 0);
      x = log10(x);
    end
    if isempty(x)
      c([1 2]) = 0;
    else
      c([1 2]) = [min(x(:)) max(x(:))];
    end
    y = get(child,'ydata');
    y = y(finite(y));
    if strcmp(get(gca,'yscale'), 'log')
      y = y(y > 0);
      y = log10(y);
    end
    if isempty(y)
      c([3 4]) = 0;
    else
      c([3 4]) = [min(y(:)) max(y(:))];
    end
    try
      z = get(child,'zdata');
      z = z(finite(z));
      if isempty(z)
	c([5 6]) = 0;
      else
	if strcmp(get(gca,'zscale'), 'log')
	  z = z(z > 0);
	  z = log10(z);
	end
	c([5 6]) = [min(z(:)) max(z(:))];
      end
    end
  end
  ax([1 3 5]) = min(ax([1 3 5]), c([1 3 5]));
  ax([2 4 6]) = max(ax([2 4 6]), c([2 4 6]));
end
dx = ax(2)-ax(1);
if dx == 0
  dx = 1;
end
dy = ax(4)-ax(3);
if dy == 0
  dy = 1;
end
dz = ax(6)-ax(5);
if dz == 0
  dz = 1;
end
ax = ax + [-dx dx -dy dy -dz dz]*pct;
if strcmp(get(gca,'xscale'), 'log')
  ax([1 2]) = 10.^(ax([1 2]));
end
if strcmp(get(gca,'yscale'), 'log')
  ax([3 4]) = 10.^(ax([3 4]));
end
if strcmp(get(gca,'zscale'), 'log')
  ax([5 6]) = 10.^(ax([5 6]));
end
% clip for 2D
ax = ax(1:length(axis));
if ~isempty(children)
  axis(ax);
end
if nargout < 1
  clear ax
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function [H,univariate] = dcc_fit_variance(data,p,o,q,gjrType,startingVals)
% Fits TARCH models for use in DCC and related estimators
%
% USAGE:
%   [H,UNIVARIATE] = dcc_fit_variance(DATA,P,O,Q,GJRTYPE)
%
% INPUTS:
%   DATA    - A column of mean zero data
%   P       - K by 1 vector of positive, scalar integers representing the number of symmetric innovations
%   O       - K by 1 vector of non-negative scalar integers representing the number of asymmetric innovations (0
%                    for symmetric processes)
%   Q       - K by 1 vector of non-negative, scalar integers representing the number of lags of conditional
%                    variance (0 for ARCH)
%   GJRTYPE - K by 1 vector of model types:
%                    1 - Model evolves in absolute values
%                    2 - Model evolves in squares [DEFAULT]
%   STARTINGVALS - [OPTIONAL] K+sum(P)+sum(O)+sum(Q) vector of starting values
%
% OUTPUTS:
%   H          - T by K matrix of conditional variances
%   UNIVARIATE - A cell array of structures used to reconstruct the variance
%
% COMMENTS:
%
%  See also TARCH, DCC, CCC_MVGARCH

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 17/4/2012

if size(startingVals,2)>size(startingVals,1)
    startingVals = startingVals';
end

[T,k] = size(data);
H = zeros(T,k);
univariate = cell(k,1);
univariteOptions = optimset('fminunc');
univariteOptions.Display = 'none';
univariteOptions.LargeScale = 'off';
offset = 0;

for i=1:k
    if ~isempty(startingVals)
        count = 1+p(i)+o(i)+q(i);
        volStartingVals = startingVals(offset + (1:count));
        offset = offset + count;
    else
        volStartingVals = [];
    end
    [parameters, ~, ht, ~, ~, scores, diagnostics] = tarch(data(:,i),p(i),o(i),q(i), [], gjrType(i), volStartingVals, univariteOptions);
    % Store output for later use
    univariate{i}.p = p(i);
    univariate{i}.o = o(i);
    univariate{i}.q = q(i);
    univariate{i}.fdata = diagnostics.fdata;
    univariate{i}.fIdata = diagnostics.fIdata;
    univariate{i}.back_cast = diagnostics.back_cast;
    univariate{i}.m = diagnostics.m;
    univariate{i}.T = diagnostics.T;
    univariate{i}.tarch_type = gjrType(i);
    univariate{i}.parameters = parameters;
    univariate{i}.ht = ht;
    univariate{i}.A = diagnostics.A;
    univariate{i}.scores = scores;
    H(:,i) = ht;
end


text:
function [x, cost, info, options] = neldermead(problem, x, options)
% Nelder Mead optimization algorithm for derivative-free minimization.
%
% function [x, cost, info, options] = neldermead(problem)
% function [x, cost, info, options] = neldermead(problem, x0)
% function [x, cost, info, options] = neldermead(problem, x0, options)
% function [x, cost, info, options] = neldermead(problem, [], options)
%
% Apply a Nelder-Mead minimization algorithm to the problem defined in
% the problem structure, starting with the population x0 if it is provided
% (otherwise, a random population on the manifold is generated). A
% population is a cell containing points on the manifold. The number of
% elements in the cell must be dim+1, where dim is the dimension of the
% manifold: problem.M.dim().
%
% To specify options whilst not specifying an initial guess, give x0 as []
% (the empty matrix).
%
% This algorithm is a plain adaptation of the Euclidean Nelder-Mead method
% to the Riemannian setting. It comes with no convergence guarantees and
% there is room for improvement. In particular, we compute centroids as
% Karcher means, which seems overly expensive: cheaper forms of
% average-like quantities might work better.
% This solver is useful nonetheless for problems for which no derivatives
% are available, and it may constitute a starting point for the development
% of other Riemannian derivative-free methods.
%
% None of the options are mandatory. See in code for details.
%
% Requires problem.M.pairmean(x, y) to be defined (computes the average
% between two points, x and y).
%
% If options.statsfun is defined, it will receive a cell of points x (the
% current simplex being considered at that iteration), and, if required,
% one store structure corresponding to the best point, x{1}. The points are
% ordered by increasing cost: f(x{1}) <= f(x{2}) <= ... <= f(x{dim+1}),
% where dim = problem.M.dim().
%
% Based on http://www.optimization-online.org/DB_FILE/2007/08/1742.pdf.
%
% See also: manopt/solvers/pso/pso

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%
%   Apr.  4, 2015 (NB):
%       Working with the new StoreDB class system.
%       Clarified interactions with statsfun and store.
%
%   Nov. 11, 2016 (NB):
%       If options.verbosity is < 2, prints minimal output.
%
%   Sep.  6, 2018 (NB):
%       Using retraction instead of exponential.

    % Verify that the problem description is sufficient for the solver.
    if ~canGetCost(problem)
        warning('manopt:getCost', ...
                'No cost provided. The algorithm will likely abort.');
    end

    % Dimension of the manifold
    dim = problem.M.dim();

    % Set local defaults here
    localdefaults.storedepth = 0;                     % no need for caching
    localdefaults.maxiter = max(2000, 4*dim);

    localdefaults.reflection = 1;
    localdefaults.expansion = 2;
    localdefaults.contraction = .5;
    % forced to .5 to enable using pairmean functions in manifolds.
    % localdefaults.shrinkage = .5;

    % Merge global and local defaults, then merge w/ user options, if any.
    localdefaults = mergeOptions(getGlobalDefaults(), localdefaults);
    if ~exist('options', 'var') || isempty(options)
        options = struct();
    end
    options = mergeOptions(localdefaults, options);

    % Start timing for initialization.
    timetic = tic();

    % If no initial simplex x is given by the user, generate one at random.
    if ~exist('x', 'var') || isempty(x)
        x = cell(dim+1, 1);
        for i = 1 : dim+1
            x{i} = problem.M.rand();
        end
    end

    % Create a store database and a key for each point.
    storedb = StoreDB(options.storedepth);
    key = cell(size(x));
    for i = 1 : dim+1;
        key{i} = storedb.getNewKey();
    end

    % Compute objective-related quantities for x, and setup a
    % function evaluations counter.
    costs = zeros(dim+1, 1);
    for i = 1 : dim+1
        costs(i) = getCost(problem, x{i}, storedb, key{i});
    end
    costevals = dim+1;

    % Sort simplex points by cost.
    [costs, order] = sort(costs);
    x = x(order);
    key = key(order);

    % Iteration counter.
    % At any point, iter is the number of fully executed iterations so far.
    iter = 0;

    % Save stats in a struct array info, and preallocate.
    % savestats will be called twice for the initial iterate (number 0),
    % which is unfortunate, but not problematic.
    stats = savestats();
    info(1) = stats;
    info(min(10000, options.maxiter+1)).iter = [];

    % Start iterating until stopping criterion triggers.
    while true

        % Make sure we don't use to much memory for the store database.
        storedb.purge();

        stats = savestats();
        info(iter+1) = stats; %#ok<AGROW>
        iter = iter + 1;

        % Start timing this iteration.
        timetic = tic();

        % Sort simplex points by cost.
        [costs, order] = sort(costs);
        x = x(order);
        key = key(order);

        % Log / display iteration information here.
        if options.verbosity >= 2
            fprintf('Cost evals: %7d\tBest cost: %+.4e\t', ...
                    costevals, costs(1));
        end

        % Run standard stopping criterion checks.
        [stop, reason] = stoppingcriterion(problem, x, options, info, iter);

        if stop
            if options.verbosity >= 1
                fprintf([reason '\n']);
            end
            break;
        end

        % Compute a centroid for the dim best points.
        xbar = centroid(problem.M, x(1:end-1));

        % Compute the direction for moving along the axis xbar - worst x.
        vec = problem.M.log(xbar, x{end});

        % Reflection step
        xr = problem.M.retr(xbar, vec, -options.reflection);
        keyr = storedb.getNewKey();
        costr = getCost(problem, xr, storedb, keyr);
        costevals = costevals + 1;

        % If the reflected point is honorable, drop the worst point,
        % replace it by the reflected point and start new iteration.
        if costr >= costs(1) && costr < costs(end-1)
            if options.verbosity >= 2
                fprintf('Reflection\n');
            end
            costs(end) = costr;
            x{end} = xr;
            key{end} = keyr;
            continue;
        end

        % If the reflected point is better than the best point, expand.
        if costr < costs(1)
            xe = problem.M.retr(xbar, vec, -options.expansion);
            keye = storedb.getNewKey();
            coste = getCost(problem, xe, storedb, keye);
            costevals = costevals + 1;
            if coste < costr
                if options.verbosity >= 2
                    fprintf('Expansion\n');
                end
                costs(end) = coste;
                x{end} = xe;
                key{end} = keye;
                continue;
            else
                if options.verbosity >= 2
                    fprintf('Reflection (failed expansion)\n');
                end
                costs(end) = costr;
                x{end} = xr;
                key{end} = keyr;
                continue;
            end
        end

        % If the reflected point is worse than the second to worst point,
        % contract.
        if costr >= costs(end-1)
            if costr < costs(end)
                % do an outside contraction
                xoc = problem.M.retr(xbar, vec, -options.contraction);
                keyoc = storedb.getNewKey();
                costoc = getCost(problem, xoc, storedb, keyoc);
                costevals = costevals + 1;
                if costoc <= costr
                    if options.verbosity >= 2
                        fprintf('Outside contraction\n');
                    end
                    costs(end) = costoc;
                    x{end} = xoc;
                    key{end} = keyoc;
                    continue;
                end
            else
                % do an inside contraction
                xic = problem.M.retr(xbar, vec, options.contraction);
                keyic = storedb.getNewKey();
                costic = getCost(problem, xic, storedb, keyic);
                costevals = costevals + 1;
                if costic <= costs(end)
                    if options.verbosity >= 2
                        fprintf('Inside contraction\n');
                    end
                    costs(end) = costic;
                    x{end} = xic;
                    key{end} = keyic;
                    continue;
                end
            end
        end

        % If we get here, shrink the simplex around x{1}.
        if options.verbosity >= 2
            fprintf('Shrinkage\n');
        end
        for i = 2 : dim+1
            x{i} = problem.M.pairmean(x{1}, x{i});
            key{i} = storedb.getNewKey();
            costs(i) = getCost(problem, x{i}, storedb, key{i});
        end
        costevals = costevals + dim;

    end

    info = info(1:iter);

    % Iteration done: return only the best point found.
    cost = costs(1);
    x = x{1};
    key = key{1};

    % Routine in charge of collecting the current iteration stats.
    function stats = savestats()
        stats.iter = iter;
        stats.cost = costs(1);
        stats.costevals = costevals;
        if iter == 0
            stats.time = toc(timetic);
        else
            stats.time = info(iter).time + toc(timetic);
        end
        % The statsfun can only possibly receive one store structure. We
        % pass the key to the best point, so that the best point's store
        % will be passed. But the whole cell x of points is passed through.
        stats = applyStatsfun(problem, x, storedb, key{1}, options, stats);
    end

end

author: NicolasBoumal

repo: manopt


text:
width = 1024;
height = 600;
mf = spx.graphics.Figures(width, height);

mf.new_figure('Figure 1');

mf.new_figure('Figure 2');

mf.new_figure('Figure 3');

mf.new_figure('Figure 4');

author: indigits

repo: sparse-plex


text:
function polar_info_esti = FastSCL_decoder(llr, L, info_bits, lambda_offset, llr_layer_vec, bit_layer_vec, psi_vec, node_type_matrix, H_crc)
%2018.1.7.14:16 Yu Y. R.
%LLR-based SCL deocoder
%Systematic polar code used
%4 nodes are identified and decoded, i.e., rate0, rate1, rep, spc,
%No BLER loss, while using fast algorithms.
%const
N = length(llr);
m = log2(N);
%memory declared
P = zeros(N - 1, L); %llr is public-used, so N - 1 is enough.
C = zeros(2 * N - 1, 2 * L);
PM = zeros(L, 1);
activepath = zeros(L, 1);
lazy_copy = zeros(m, L);
%initialize the 1st SC decoder
activepath(1) = 1;
lazy_copy(:, 1) = 1;
%decoding starts
%default: in the case of path clone in REP and leaf node,
%origianl always corresponds to bit 0s, while the new path bit 1s.
for i_node = 1 : size(node_type_matrix, 1)
    current_index = node_type_matrix(i_node, 1);
    llr_layer = llr_layer_vec(current_index);
    M = node_type_matrix(i_node, 2);%number of leaf nodes in this constituent node
    reduced_layer = log2(M); %For LLR recursion, this reduced layer denotes where to stop; For Bit recursion, this denotes where to start.
    psi = psi_vec(i_node);%This psi is used for bits recursion
    psi_mod_2 = mod(psi, 2);%To decide whether the bit recuision should continue. 1 : continue; 0 : stop.
    for l_index = 1 : L
        if activepath(l_index) == 0
            continue;
        end
        switch current_index
            case 1
                index_1 = lambda_offset(m);
                for beta = 0 : index_1 - 1
                    P(beta + index_1, l_index) = sign(llr(beta + 1)) *  sign(llr(beta + index_1 + 1)) * min(abs(llr(beta + 1)), abs(llr(beta + index_1 + 1)));
                end
                for i_layer = m - 2 : -1 : reduced_layer
                    index_1 = lambda_offset(i_layer + 1);
                    index_2 = lambda_offset(i_layer + 2);
                    for beta = index_1 : index_2 - 1
                        P(beta, l_index) = sign(P(beta + index_1, l_index)) *  sign(P(beta + index_2, l_index)) * ...
                            min(abs(P(beta + index_1, l_index)), abs(P(beta + index_2, l_index)));
                    end
                end
            case N/2 + 1
                index_1 = lambda_offset(m);
                for beta = 0 : index_1 - 1
                    x_tmp = C(beta + index_1, 2 * l_index - 1);
                    P(beta + index_1, l_index) = (1 - 2 * x_tmp) * llr(beta + 1) + llr(beta + index_1 + 1);
                end
                for i_layer = m - 2 : -1 : reduced_layer
                    index_1 = lambda_offset(i_layer + 1);
                    index_2 = lambda_offset(i_layer + 2);
                    for beta = index_1 : index_2 - 1
                        P(beta, l_index) = sign(P(beta + index_1, l_index)) *  sign(P(beta + index_2, l_index)) * ...
                            min(abs(P(beta + index_1, l_index)), abs(P(beta + index_2, l_index)));
                    end
                end
            otherwise
                index_1 = lambda_offset(llr_layer + 1);
                index_2 = lambda_offset(llr_layer + 2);
                for beta = index_1 : index_2 - 1
                    %lazy copy
                    P(beta, l_index) = (1 - 2 * C(beta, 2 * l_index - 1)) * P(beta + index_1, lazy_copy(llr_layer + 2, l_index)) + P(beta + index_2, lazy_copy(llr_layer + 2, l_index));
                end
                for i_layer = llr_layer - 1 : -1 : reduced_layer
                    index_1 = lambda_offset(i_layer + 1);
                    index_2 = lambda_offset(i_layer + 2);
                    for beta = index_1 : index_2 - 1
                        P(beta, l_index) = sign(P(beta + index_1, l_index)) *  sign(P(beta + index_2, l_index)) * ...
                            min(abs(P(beta + index_1, l_index)), abs(P(beta + index_2, l_index)));
                    end
                end
        end
    end%LLR calculations. You may fold thiis code block.
    index_vec = M : 2 * M - 1;
    switch node_type_matrix(i_node, 3)
        case -1%****RATE 0*****
            x = zeros(M, 1);
            for l_index = 1 : L
                if activepath(l_index) == 0
                    continue;
                end
                C(index_vec, psi_mod_2 + 2 * l_index - 1) = x;
                PM(l_index) = PM(l_index) + sum(log(1 + exp(-P(index_vec, l_index))));
                %Do not forget updating path metric in rate 0 nodes
                %No copy in Rate 0 node
            end
        case 1%*****RATE 1*****
            %input LLR: P(index_vec, :)
            [PM, activepath, selected_subcodeword, lazy_copy] = Rate1(M, P(index_vec, :), activepath, PM, L, lazy_copy);
            %output: estimated sub polar codeword: 'selected_subcodeword'.
            for l_index = 1 : L
                if activepath(l_index) == 0
                    continue
                end
                C(index_vec, 2 * l_index - 1 + psi_mod_2) = selected_subcodeword(:, l_index);
            end
        case 2%*****REP*****
           [PM, activepath, lazy_copy, selected_subcodeword] = Rep(PM, P(index_vec, :), activepath, L, M, lazy_copy);
           for l_index = 1 : L
               if activepath(l_index) == 0
                   continue
               end
               C(index_vec, 2 * l_index - 1 + psi_mod_2) = selected_subcodeword(:, l_index);
           end
        case 3%*******SPC*******
            [PM, activepath, selected_subcodeword, lazy_copy] = SPC(index_vec, P(index_vec, :), activepath, PM, L, lazy_copy);
            for l_index = 1 : L
                if activepath(l_index) == 0
                    continue
                end
                C(index_vec, 2 * l_index - 1 + psi_mod_2) = selected_subcodeword(:, l_index);
            end
    end
    %Internal Bits Update for each active path
    bit_layer = bit_layer_vec(current_index + M - 1);
    for l_index = 1 : L
        if activepath(l_index) == 0
            continue;
        end
        if psi_mod_2 == 1%right child of its mother
            for i_layer = reduced_layer : bit_layer - 1
                index_1 = lambda_offset(i_layer + 1);
                index_2 = lambda_offset(i_layer + 2);
                for beta = index_1 : index_2 - 1
                    C(beta + index_1, 2 * l_index) = mod(C(beta, 2 * lazy_copy(i_layer + 1, l_index) - 1) + C(beta, 2 * l_index), 2);
                    C(beta + index_2, 2 * l_index) = C(beta, 2 * l_index);
                end
            end
            index_1 = lambda_offset(bit_layer + 1);
            index_2 = lambda_offset(bit_layer + 2);
            for beta = index_1 : index_2 - 1
                C(beta + index_1, 2 * l_index - 1) = mod(C(beta, 2 * lazy_copy(bit_layer + 1, l_index) - 1) + C(beta, 2 * l_index), 2);
                C(beta + index_2, 2 * l_index - 1) = C(beta, 2 * l_index);
            end
        end
    end
    %lazy copy
    if i_node < size(node_type_matrix, 1)
        for i_layer = 1 : llr_layer_vec(current_index + M) + 1
            for l_index = 1 : L
                lazy_copy(i_layer, l_index) = l_index;
            end
        end
    end
end
%select the best path
[~, path_ordered] = sort(PM);
for l_index = 1 : L
    path_num = path_ordered(l_index);
    x = C(end - N + 1 : end, 2 * path_num - 1);
    info_with_crc = x(info_bits);
    err = sum(mod(H_crc * info_with_crc, 2));
    if err == 0
        polar_info_esti = info_with_crc;
        break;
    else
        if l_index == L
            x = C(end - N + 1 : end, 2 * path_ordered(1) - 1);
            polar_info_esti = x(info_bits);
        end
    end
end
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
function history(dat)
% Display history for image_vector object

u = '_________________________________________________';
disp(u)
disp('History for image data object')
disp(u)

if ischar(dat.dat_descrip), fprintf('Data description:\n%s\n', dat.dat_descrip); end

disp(char(dat.history))

disp(u)

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function write_sicd_imsubhdr(obj, n)
%WRITE_SICD_IMSUBHDR Writes the NITF image segment sub-header for SICD file
%
% It's not expected that the typical user would call this function, rather
% it is called as part of a bigger NITF/SICD writer.
%
% Inputs:
%       obj:       This is the 'hidden' object-specific handle.  It is
%                  similar to 'this' in Java.
%
% Outputs:
%       lengthWritten: The number of bytes written to disk
%
% References:
%       MIL-STD-2500C, Department of Defense Interface Standard
%                      National Imagery Transmission Format version 2.1
%                      01 May 2006
%
% Written by: Tom Krauss and Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Always writing version 2.1 NITF

fwriten(obj.FID, 'IM', 2);
if obj.NumIS==1
    IID1 = 'SICD000';
else
    IID1 = ['SICD' num2str(n,'%03d')];
end
fwriten(obj.FID, IID1, 10); % IID1

if isfield(obj.sicdmeta,'ImageCreation')&&isfield(obj.sicdmeta.ImageCreation,'DateTime')
    fdt = datestr(obj.sicdmeta.ImageCreation.DateTime,'yyyymmddHHMMSS');
else
    fdt = '';
end
fwriten(obj.FID, fdt, 14); % IDATIM

fwriten(obj.FID, '', 17); % TGTID

if isfield(obj.sicdmeta,'CollectionInfo')&&isfield(obj.sicdmeta.CollectionInfo,'CoreName')
    IID2 = obj.sicdmeta.CollectionInfo.CoreName;
else
    IID2 = 'Unknown';
end
fwriten(obj.FID, ['SICD: ' IID2], 80); % IID2

obj.write_sicd_security_tags();  % Security tags go here
fwriten(obj.FID, '0', 1);  % ENCRYPT
if isfield(obj.sicdmeta,'CollectionInfo')&&isfield(obj.sicdmeta.CollectionInfo,'CollectorName')
    ISORCE = obj.sicdmeta.CollectionInfo.CollectorName;
else
    ISORCE = '';
end
fwriten(obj.FID, ['SICD: ' ISORCE], 42); % ISORCE
fwriten(obj.FID, uint64(obj.NumRowsIS(n)), 8); % NROWS
fwriten(obj.FID, uint64(obj.sicdmeta.ImageData.NumCols), 8); % NCOLS

switch obj.sicdmeta.ImageData.PixelType
    case 'RE32F_IM32F'
        PVTYPE = 'R';
        ABPP = 32;
        ISUBCAT1 = 'I';
        ISUBCAT2 = 'Q';
    case 'RE16I_IM16I'
        PVTYPE = 'SI';
        ABPP = 16;
        ISUBCAT1 = 'I';
        ISUBCAT2 = 'Q';
    case 'AMP8I_PHS8I'
        PVTYPE = 'INT';
        ABPP = 8;
        ISUBCAT1 = 'M';
        ISUBCAT2 = 'P';
end
fwriten(obj.FID, PVTYPE, 3); % PVTYPE (here signed integer)
fwriten(obj.FID, 'NODISPLY',8); % IREP
fwriten(obj.FID, 'SAR', 8); % ICAT
fwriten(obj.FID, uint64(ABPP),2);  % ABPP
fwriten(obj.FID, 'R', 1);  % PJUST
fwriten(obj.FID, 'G', 1);  % ICORDS
if isfield(obj.sicdmeta,'GeoData')&&isfield(obj.sicdmeta.GeoData,'ImageCorners')
    geo_format = {'num_units',3,'include_symbols',false};
    % TODO: Adjust to handle multiple image segments
    geolo = [latlonstr(obj.sicdmeta.GeoData.ImageCorners.ICP.FRFC.Lat,'lat',geo_format{:}) ...
        latlonstr(obj.sicdmeta.GeoData.ImageCorners.ICP.FRFC.Lon,'lon',geo_format{:}) ...
        latlonstr(obj.sicdmeta.GeoData.ImageCorners.ICP.FRLC.Lat,'lat',geo_format{:}) ...
        latlonstr(obj.sicdmeta.GeoData.ImageCorners.ICP.FRLC.Lon,'lon',geo_format{:}) ...
        latlonstr(obj.sicdmeta.GeoData.ImageCorners.ICP.LRLC.Lat,'lat',geo_format{:}) ...
        latlonstr(obj.sicdmeta.GeoData.ImageCorners.ICP.LRLC.Lon,'lon',geo_format{:}) ...
        latlonstr(obj.sicdmeta.GeoData.ImageCorners.ICP.LRFC.Lat,'lat',geo_format{:}) ...
        latlonstr(obj.sicdmeta.GeoData.ImageCorners.ICP.LRFC.Lon,'lon',geo_format{:})];
else
    geolo = '';
end
fwriten(obj.FID, geolo, 60); % IGEOLO

fwriten(obj.FID, uint64(0), 1);  % NICOM (# image comments)
fwriten(obj.FID, 'NC', 2);  % IC (No image compression)
fwriten(obj.FID, uint64(2), 1);  % NBANDS

% 'Band' 1: in phase
fwriten(obj.FID, '',  2); % IREPBAND1
fwriten(obj.FID, ISUBCAT1, 6); % ISUBCAT1
fwriten(obj.FID, 'N', 1); % IFC1
fwriten(obj.FID, '',  3); % IMFLT1
fwriten(obj.FID, '0', 1); % NLUTS1

% 'Band' 2: quadrature
fwriten(obj.FID, '',  2); % IREPBAND2
fwriten(obj.FID, ISUBCAT2, 6); % ISUBCAT2
fwriten(obj.FID, 'N', 1); % IFC2
fwriten(obj.FID, '',  3); % IMFLT2
fwriten(obj.FID, '0', 1); % NLUTS2

fwriten(obj.FID, '0', 1); % ISYNC, always 0
fwriten(obj.FID, 'P', 1);  % IMODE (band interleaved by pixel)
fwriten(obj.FID, uint64(1), 4);  % NBPR - blocks per row
fwriten(obj.FID, uint64(1), 4);  % NBPC - blocks per column
if obj.sicdmeta.ImageData.NumCols>8192
    NPPBH = 0; % (zero means "use NCOLS")
else
    NPPBH = obj.sicdmeta.ImageData.NumCols;
end
fwriten(obj.FID, uint64(NPPBH), 4);  % NPPBH - Number of Pixels Per Block Horizontal
if obj.NumRowsIS(n)>8192
    NPPBV = 0; % (zero means "use NROWS")
else
    NPPBV = obj.NumRowsIS(n);
end
fwriten(obj.FID, uint64(NPPBV), 4);  % NPPBV - Number of Pixels Per Block Vertical
fwriten(obj.FID, uint64(ABPP),2);  % NBPP
fwriten(obj.FID, uint64(n), 3);  % IDLVL
fwriten(obj.FID, uint64(n-1), 3);  % IALVL
if n==1
    fwriten(obj.FID, uint64(0), 5); % ILOC (row) - Image location offsets
else
    fwriten(obj.FID, uint64(obj.NumRowsIS(n-1)), 5); % ILOC (row) - Image location offsets
end
fwriten(obj.FID, uint64(0), 5); % ILOC (column) - Image location offsets
fwriten(obj.FID, '1.0 ', 4);     % IMAG - image magnification
fwriten(obj.FID, uint64(0), 5);  % UDIDL - Generally not used for SICD
fwriten(obj.FID, uint64(0), 5);  % IXSHDL - Generally not used for SICD

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function test04 ( dim_num, n, z )

%*****************************************************************************80
%
%% TEST04 tests GAMMA_MEASURE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04\n' );
  fprintf ( 1, '  GAMMA_MEASURE computes the GAMMA measure of quality.\n' );
  fprintf ( 1, '  The mesh ratio               Gamma = %14f\n', ...
    gamma_measure ( dim_num, n, z ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [FileName, FileType, isAnatomy] = file_fullpath( FileName )
% FILE_FULLPATH: Return the full filename from a relative filename
%
% USAGE:  [FileName, FileType, isAnatomy] = file_fullpath( FileName )
%         [FileName, FileType, isAnatomy] = file_fullpath( FileNames )

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2011-2022

% Empty input
if isempty(FileName)
    disp('FILE_FULLPATH> Input filename is empty');
    return;
% List of files: recursive calls
elseif iscell(FileName)
    for i = 1:length(FileName)
        FileName{i} = file_fullpath( FileName{i} );
    end
    return
end

% Get protocol folders
ProtocolInfo = bst_get('ProtocolInfo');
% Get file type
FileType = file_gettype(FileName);
isAnatomy = 0;
tmpFile = [];
% Add protocol path
switch lower(FileType)
    case {'brainstormsubject', 'subject', 'subjectimage', 'anatomy', 'scalp', 'outerskull', 'innerskull', 'cortex', 'fibers', 'fem', 'other', 'tess'}
        if ~file_exist(FileName)
            tmpFile = bst_fullfile(ProtocolInfo.SUBJECTS, FileName);
        end
        isAnatomy = 1;
    case {'brainstormstudy', 'study', 'studysubject', 'condition', 'rawcondition', 'channel', 'headmodel', 'data', 'rawdata', 'results', 'kernel', 'pdata', 'presults', 'noisecov', 'ndatacov', 'dipoles', 'timefreq', 'spectrum', 'ptimefreq', 'pspectrum', 'matrix', 'pmatrix', 'proj', 'image', 'video', 'videolink', 'spike', 'dirspikes'}
        if ~file_exist(FileName)
            tmpFile = bst_fullfile(ProtocolInfo.STUDIES, FileName);
        end
    case 'link'
        tmpFile = file_resolve_link(FileName);
    otherwise
        error('Unsupported file type.');
end

% If something was changed
if ~isempty(tmpFile)
    % Check for file existence
    if ~isempty(dir(tmpFile))
        FileName = tmpFile;
    else
        warning(['File not found: ', tmpFile]);
    end
end
end

author: brainstorm-tools


text:
function [ state ] = movePTP_ConditionalTorque_HomeJointSpace( t_Kuka , relVel,joints_indices,max_torque,min_torque)
%% This function is used for performing point to point motion
%  in joint space to the home position of the robot, for the KUKA iiwa 7 R 800.
% this motion is interruptible when one or more of the joints-torques exceed the predefined limits

%% Syntax:
% [ state ] = movePTP_ConditionalTorque_HomeJointSpace( t_Kuka , relVel,joints_indices,max_torque,min_torque)

%% About:
% This function is used to move the robot to home position, the motion is
% performed in joint space.

%% Arreguments:
% t_Kuka: is the TCP/IP connection
% relVel: is a double, from zero to one, specifies the over-ride relative
% velocity.
% joints_indices: is a vector of the indices of the joints where the
% torques limits are to be imposed, the joints are indexed starting from one.
% max_torque: is a vector of the maximum torque limits for the joints
% specified in the (joints_indices) vector.
% min_torque: is a vecotr of the minimum torque limits for the joints
% specified in the (joints_indices) vector.

%% Return value:
% state: a number equals to one if the motion is completed sccessfully or a
% zero if the motion was interrupted due to external contact with the
% robot.
% if there is an error, the return value is minus one

% Copy right, Mohammad SAFEEA, 9th of April 2018

if((size(relVel,1)==1)&&(size(relVel,2)==1))
else
    disp('Error, relative velocity is a double and shall not be an array');
    state=-1;
    return;
end
% Check if the inputs "joints_indices,max_torque,min_torque" are vectors
if(isVec(joints_indices)==0)
    disp('Error, joints_indices shall be a vector');
    state=-1;
    return;
end
if(isVec(max_torque)==0)
    disp('Error, max_torque shall be a vector');
    state=-1;
    return;
end
if(isVec(min_torque)==0)
    disp('Error, min_torque shall be a vector');
    state=-1;
    return;
end
joints_indices=colVec(joints_indices);
max_torque=colVec(max_torque);
min_torque=colVec(min_torque);

% Check if size of vectrs are equal
j_num=size(joints_indices,1);
n=size(min_torque,1);
if(j_num==n)
else
    disp('Error, sizes of vectors "joints_indices" and "min_torque" shall be equal');
    state=-1;
    return;
end
n=size(max_torque,1);
if(j_num==n)
else
    disp('Error, sizes of vectors "joints_indices" and "max_torque" shall be equal');
    state=-1;
    return;
end

theCommand=['jRelVel_',num2str(relVel),'_']; % set over ride.
fprintf(t_Kuka, theCommand);
message=fgets(t_Kuka);

jPos=[];
for i=1:7
    jPos{i}=0;
end

sendJointsPositions( t_Kuka ,jPos); % send destination joint positions.

theCommand='doPTP!inJS';
for i=1:j_num
    datemp=joints_indices(i)-1;
    theCommand=[theCommand,'_',num2str(datemp)];
    datemp=min_torque(i);
    theCommand=[theCommand,'_',num2str(datemp)];
    datemp=max_torque(i);
    theCommand=[theCommand,'_',num2str(datemp)];
end

fprintf(t_Kuka, theCommand); % start the point to point motion.
message=fgets(t_Kuka);
if checkAcknowledgment(message)
else
    disp('Error, could not perform the motion due to the acknowledgment message was not received');
    disp('Check connection between PC and Robot');
    state=-1;
    return;
end

%% Check the return from the server
message='';

while true
    message=fgets(t);
    daSize=max(size(message));
    if(daSize>2)
        if(strfind(message,'done')==1)
            state=1;
            return;
        end
        if(strfind(message,'interrupted')==1)
            state=0;
            return;
        end
    end
end

end

repo: KST-Kuka-Sunrise-Toolbox


text:
function check = half_normal_check ( a, b )

%*****************************************************************************80
%
%% HALF_NORMAL_CHECK checks the parameters of the Half Normal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, logical CHECK, is true if the parameters are legal.
%
  if ( b <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HALF_NORMAL_CHECK - Fatal error!\n' );
    fprintf ( 1, '  B <= 0.\n' );
    error ( 'HALF_NORMAL_CHECK - Fatal error!' );
  end

  check = 1;

  return
end


text:
function [H,Hnorm,inv_Hnorm] = compute_homography(m,M);

%compute_homography
%
%[H,Hnorm,inv_Hnorm] = compute_homography(m,M)
%
%Computes the planar homography between the point coordinates on the plane (M) and the image
%point coordinates (m).
%
%INPUT: m: homogeneous coordinates in the image plane (3xN matrix)
%       M: homogeneous coordinates in the plane in 3D (3xN matrix)
%
%OUTPUT: H: Homography matrix (3x3 homogeneous matrix)
%        Hnorm: Normalization matrix used on the points before homography computation
%               (useful for numerical stability is points in pixel coordinates)
%        inv_Hnorm: The inverse of Hnorm
%
%Definition: m ~ H*M where "~" means equal up to a non zero scalar factor.
%
%Method: First computes an initial guess for the homography through quasi-linear method.
%        Then, if the total number of points is larger than 4, optimize the solution by minimizing
%        the reprojection error (in the least squares sense).
%
%
%Important functions called within that program:
%
%comp_distortion_oulu: Undistorts pixel coordinates.
%
%compute_homography.m: Computes the planar homography between points on the grid in 3D, and the image plane.
%
%project_points.m: Computes the 2D image projections of a set of 3D points, and also returns te Jacobian
%                  matrix (derivative with respect to the intrinsic and extrinsic parameters).
%                  This function is called within the minimization loop.

Np = size(m,2);

if size(m,1)<3,
   m = [m;ones(1,Np)];
end;

if size(M,1)<3,
   M = [M;ones(1,Np)];
end;

m = m ./ (ones(3,1)*m(3,:));
M = M ./ (ones(3,1)*M(3,:));

% Prenormalization of point coordinates (very important):
% (Affine normalization)

ax = m(1,:);
ay = m(2,:);

mxx = mean(ax);
myy = mean(ay);
ax = ax - mxx;
ay = ay - myy;

scxx = mean(abs(ax));
scyy = mean(abs(ay));

Hnorm = [1/scxx 0 -mxx/scxx;0 1/scyy -myy/scyy;0 0 1];
inv_Hnorm = [scxx 0 mxx ; 0 scyy myy; 0 0 1];

mn = Hnorm*m;

% Compute the homography between m and mn:

% Build the matrix:

L = zeros(2*Np,9);

L(1:2:2*Np,1:3) = M';
L(2:2:2*Np,4:6) = M';
L(1:2:2*Np,7:9) = -((ones(3,1)*mn(1,:)).* M)';
L(2:2:2*Np,7:9) = -((ones(3,1)*mn(2,:)).* M)';

if Np > 4,
	L = L'*L;
end;

[U,S,V] = svd(L);

hh = V(:,9);
hh = hh/hh(9);

Hrem = reshape(hh,3,3)';
%Hrem = Hrem / Hrem(3,3);

% Final homography:

H = inv_Hnorm*Hrem;

if 0,
   m2 = H*M;
   m2 = [m2(1,:)./m2(3,:) ; m2(2,:)./m2(3,:)];
   merr = m(1:2,:) - m2;
end;

%keyboard;

%%% Homography refinement if there are more than 4 points:

if Np > 4,

   % Final refinement:
   hhv = reshape(H',9,1);
   hhv = hhv(1:8);

   for iter=1:10,

		mrep = H * M;

		J = zeros(2*Np,8);

		MMM = (M ./ (ones(3,1)*mrep(3,:)));

		J(1:2:2*Np,1:3) = -MMM';
		J(2:2:2*Np,4:6) = -MMM';

		mrep = mrep ./ (ones(3,1)*mrep(3,:));

		m_err = m(1:2,:) - mrep(1:2,:);
		m_err = m_err(:);

		MMM2 = (ones(3,1)*mrep(1,:)) .* MMM;
		MMM3 = (ones(3,1)*mrep(2,:)) .* MMM;

		J(1:2:2*Np,7:8) = MMM2(1:2,:)';
		J(2:2:2*Np,7:8) = MMM3(1:2,:)';

		MMM = (M ./ (ones(3,1)*mrep(3,:)))';

		hh_innov  = inv(J'*J)*J'*m_err;

		hhv_up = hhv - hh_innov;

		H_up = reshape([hhv_up;1],3,3)';

		%norm(m_err)
		%norm(hh_innov)

		hhv = hhv_up;
      H = H_up;

   end;

end;

if 0,
   m2 = H*M;
   m2 = [m2(1,:)./m2(3,:) ; m2(2,:)./m2(3,:)];
   merr = m(1:2,:) - m2;
end;

return;

%test of Jacobian

mrep = H*M;
mrep = mrep ./ (ones(3,1)*mrep(3,:));

m_err = mrep(1:2,:) - m(1:2,:);
figure(8);
plot(m_err(1,:),m_err(2,:),'r+');
std(m_err')

repo: ekfmonoslam


text:
function cdf = erlang_cdf ( x, a, b, c )

%*****************************************************************************80
%
%% ERLANG_CDF evaluates the Erlang CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%
%    Input, real A, B, integer C, the parameters of the PDF.
%    0.0 < B.
%    0 < C.
%
%    Output, real CDF, the value of the CDF.
%
  if ( x < a )

    cdf = 0.0;

  else

    x2 = ( x - a ) / b;
    p2 = c;

    cdf = gamma_inc ( p2, x2 );

  end

  return
end


text:
function rmPlotMap (v,  saveFlag, method, fieldname, prf_size, fieldRange)
% rmPlotMap (v, fieldname, saveFlag, prf_size, fieldRange)
%
%   Purpose:
%       Visualize the responses within an ROI in stimulus-referred space.
%       The stimulus referred space is created via a retinotopic model,
%       which gives an x, y, and sigma value to each voxel. The responses
%       come from the current scan - it does not need to be the scan used
%       to make the pRF model. The model of each voxel's RF is then scaled
%       to the current scan's  map. The models from all the maps are then
%       added together and plotted as figure in stimulus space.
%
%   Note: works only in Gray view
%
%
%
%   INPUT
%       v:           mrVista view structure
%       fieldname:   'ph', 'co', 'amp', or 'map': the parameter used to
%                           scale the RF plot
%       saveFlag:    boolean (if true, then save the plot as a jpg)
%       prf_size:    boolean (if true, use prf sigma to generate RFs,
%                           if false assume fixed size for all pRFs)
%       fieldRange:  deg (size of visual field to plot)
%
%
%
%   6/19/2008: JW wrote it, adpating from KA's rmPlotCoverage
%   7/3/2008:  Divided the image at each point by the total pRF coverage of
%              that point in space. This puts the image in  % signal units.
%

% Check the arguments and set defaults
if ieNotDefined('v'), v = getCurView; end
if ieNotDefined('fieldname'), fieldname = v.ui.displayMode; end
if ieNotDefined('saveFlag'), saveFlag = true;  end
if ieNotDefined('prf_size'), prf_size = 'from model';  end
if ieNotDefined('fieldRange'), fieldRange = 15;  end % deg
if ieNotDefined('nSamples'), nSamples = 75;  end
if ieNotDefined('normResponses'), normResponses = true;  end
if ieNotDefined('method'), method = 0; end %0 = sum

if isequal(saveFlag, 'dialog')
	% get parameters from a dialog
	[prf_size method normResponses fieldRange nSamples saveFlag ok] = rmPlotMapParams;
    if ~ok, return; end
end

% Get pRF model
try
    rmModel   = viewGet(v,'rmSelectedModel');
    rmfname   = viewGet(v,'rmfile');
    [foo,rmfname,foo,foo] = fileparts(rmfname);
catch
    error('Need to load retModel into curent view');
end

% Get ROI
try
    ROIcoords = getCurROIcoords(v);
    ROIname = v.ROIs(v.selectedROI).name;
catch
    error('Need to select ROI in GUI')
end

% Get scan
curScan = getCurScan(v);

% If cothresh is set in GUI, use it to restrict ROI
try
    co  = getCurDataROI(v,'co',curScan,ROIcoords);
    cothresh = viewGet(v, 'cothresh');
    ROIcoords = ROIcoords(:, co > cothresh);
end

nVoxels = size(ROIcoords,2);

% Get the map
curData  = getCurDataROI(v,fieldname,curScan,ROIcoords);

% Get X, Y, and sigma for each voxel
[tmp1 tmp2 indices] = intersectCols(ROIcoords, v.coords);
clear tmp1 tmp2;
subSize = rmModel.sigma.major(indices);
subX = rmModel.x0(indices);
subY = rmModel.y0(indices);
% if prf_size is false, then make all pRFs the same size, instead of using
%   the fitted sigmas.
if strcmp(prf_size, 'equal for every voxel')
   subSize=ones(size(subSize))*0.5;
end

% Set up stimulus-referred visual field
x = linspace(-fieldRange,fieldRange,nSamples);
[X,Y] = meshgrid(x,x);
mask = makecircle(size(X,1));

% Build pRF for each voxel:
%all_models = rfGaussian2d(single(X(:)),single(Y(:)),...
%    single(subSize),single(subSize),single(0), single(subX),single(subY));

% Plot the total coverage of the RF.
RFcov = zeros(nSamples^2,1);
for ii = 1:nVoxels
    thisModel = rfGaussian2d(single(X(:)),single(Y(:)),...
    single(subSize(ii)),single(subSize(ii)),single(0), single(subX(ii)),single(subY(ii)));
    if method == 0,
        RFcov = RFcov + thisModel;
    else
        RFcov = max(RFcov, thisModel);
    end
end

RFcov = reshape(RFcov,[1 1].*sqrt(numel(RFcov)));

h = figure;
subplot(2,1,1)
imagesc (X(1,:),Y(:,1),RFcov .* mask);
colorbar
axis equal tight;
title([ROIname, ', Visual field coverage (sum)']);

% Plot the current map in stimulus-referred space
%   (This is the main point of the function)
RF = zeros(nSamples^2,1);
for ii = 1:nVoxels
    thisModel = rfGaussian2d(single(X(:)),single(Y(:)),...
        single(subSize(ii)),single(subSize(ii)),single(0), single(subX(ii)),single(subY(ii)));
    RF = RF + thisModel * curData(ii);
end
% convert RF model from 1D to 2D
RF = reshape(RF,[1 1].*sqrt(numel(RF)));

% divide the movie by the coverage map to normalize (approximatley) to % signal
if normResponses == true, RF = RF ./ RFcov; end

%set color range
imMax = max(max(RF(:)), -min(RF(:)));
imMin = -imMax;

% plotting
figure(h);
headerStr = [ROIname,', ', fieldname, ', scan', num2str(curScan) ];
subplot(2,1,2)
imagesc(X(1,:),Y(:,1),RF .* mask, [imMin imMax])
colorbar
axis equal tight;
title(headerStr);
xlabel(rmfname);

% saving
if saveFlag,
    fname = [headerStr,  '.jpg'];
    saveas(gcf, fname)
end

return

% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function [prf_size method normResponses fieldRange nSamples saveFlag ok] = rmPlotMapParams;
%% dialog to get parameters for rmPlotParam.
dlg(1).fieldName = 'prf_size';
dlg(end).style = 'listbox';
dlg(end).list = {'from model', 'equal for every voxel'};
dlg(end).string = 'pRF sigma';
dlg(end).value = 1;

dlg(end+1).fieldName = 'fieldRange';
dlg(end).style = 'number';
dlg(end).string = 'Visual Field Range (deg)?';
dlg(end).value = '20';

dlg(end+1).fieldName = 'nSamples';
dlg(end).style = 'number';
dlg(end).string = 'Num Samples?';
dlg(end).value = '75';

dlg(end+1).fieldName = 'normResponses';
dlg(end).style = 'number';
dlg(end).string = 'Normalize responses to visual field coverage?';
dlg(end).value = '1';

dlg(end+1).fieldName = 'method';
dlg(end).style = 'number';
dlg(end).string = 'Use sum (0) or max (1) for pRF coverage map?';
dlg(end).value = 0;

dlg(end+1).fieldName = 'saveFlag';
dlg(end).style = 'checkbox';
dlg(end).string = 'Save Image?';
dlg(end).value = 0;

[resp ok] = generalDialog(dlg, mfilename);
if ~ok
	disp('User Aborted.')
    prf_size = 0;
    method = 0;
    normResponses = 0;
    fieldRange = 0;
    nSamples = 0;
    saveFlag = 0;
    return;
end

prf_size = resp.prf_size;
method = resp.method;
normResponses = resp.normResponses;
fieldRange = resp.fieldRange;
nSamples = resp.nSamples;
saveFlag = resp.saveFlag;

return


text:
function dx = repmat_der(~, varargin)
%REPMAT_DER
%   REPMAT_DER(X, N, DZDY)
%   REPMAT_DER(X, D1, D2, ..., DZDY)
%   REPMAT_DER(X, D, DZDY)
%   Derivative of REPMAT function, w.r.t. first input. Same syntax as
%   native REPMAT, plus derivative.

% Copyright (C) 2016 Joao F. Henriques.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

  dzdy = varargin{end} ;

  assert(numel(varargin) >= 2, 'Invalid syntax.') ;

  % collect repetitions for each dimension in a single vector
  if numel(varargin) == 2
    reps = varargin{1} ;
    if isscalar(reps)
      reps = [reps, reps] ;
    end
  else
    reps = [varargin{1:end-1}] ;
  end

  % iterate dimensions, summing derivative across all repetitions in that
  % dimension
  dx = dzdy ;
  for dim = 1:numel(reps)
    if reps(dim) > 1
      % split dimension in 2: one for the original, one for the repetitions
      sz = size(dx) ;
      dx = reshape(dx, [sz(1:dim-1), sz(dim) / reps(dim), reps(dim), sz(dim+1:end)]) ;

      % sum across the repetitions dimension
      dx = sum(dx, dim + 1) ;

      % merge the two dimensions again
      dx = reshape(dx, [sz(1:dim-1), sz(dim) / reps(dim), sz(dim+1:end)]) ;
    end
  end
end

repo: MFT


text:
function [uOutput] = scecimp(nFunction, sFilename)
% import SCEC (Caltech) Data Center format
% Filter function switchyard
%%%%     CHANGE THESE LINES %%%%%%%%%%%
if nFunction == FilterOp.getDescription
  uOutput = 'SCEC (Caltech) Data Center format';
elseif nFunction == FilterOp.getWebpage
  uOutput = 'SCEC.htm'; %%%%    DO NOT CHANGE %%%%%%%%%%%
elseif nFunction == FilterOp.importCatalog
    % Read formated data
    mData = textread(sFilename, '%s', 'delimiter', '\n', 'whitespace', '');
    % Create empty catalog
    uOutput = zeros(length(mData), 9);
    % Loop through all lines of catalog and convert them
    mData = char(mData);
    l = find( mData == ' ' );
    mData(l) = '0';
    errc = 1;

    try % this is the fast method ....
        i = 1:length(mData(:,1));
        [uOutput] = readvalues(mData,i,i);

    catch  % this is the line by line method ..
        for i = 1:length(mData(:,1))
            if rem(i,100) == 0 ; disp([ num2str(i) ' of ' num2str(length(mData)) ' events processed ']); end
            try
                uOutput(i,:) = readvalues(mData,i,1);
            catch
                errc = errc + 1;
                if errc == 100
                    if stoploop
                        return
                    end
                end
                msg.dbfprintf('Import: Problem in line %d of %s. Line ignored.\n',i, sFilename);
                uOutput(i,:)=nan;
            end
        end

    end
    l = isnan(uOutput(:,1));
    uOutput(l,:) = [];

end

%%%%%%%%%%%%%%%%%%%%%%
%%%%     CHANGE THESE LINES %%%%%%%%%%%

function [uOutput] = readvalues(mData,i,k)

  uOutput(k,1) = str2num(mData(i,34:37)) - str2num(mData(i,39:43))/60 ;  % Longitude
uOutput(k,2) = str2num(mData(i,26:27)) + str2num(mData(i,29:33))/60 ;  % Latitude
uOutput(k,3) = str2num(mData(i,1:4));    % Year
uOutput(k,4) = str2num(mData(i,6:7));    % Month
uOutput(k,5) = str2num(mData(i,9:10));   % Day
uOutput(k,6) = str2num(mData(i,47:49));  % Magnitude
uOutput(k,7) = str2num(mData(i,55:59));  % Depth
uOutput(k,8) = str2num(mData(i,13:14));  % Hour
uOutput(k,9) = str2num(mData(i,16:17));  % Minute
% convert to decimal years
uOutput(k,3) = decyear([uOutput(k,3) uOutput(k,4) uOutput(k,5) uOutput(k,8) uOutput(k,9)]);

% convert to decimal years
uOutput(k,3) = decyear([uOutput(k,3) uOutput(k,4) uOutput(k,5) uOutput(k,8) uOutput(k,9)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    DO NOT CHANGE %%%%%%%%%%%

function [mystop] = stoploop()

ButtonName=questdlg('More than 100 lines could not be read. Continue?', ...
    'Interrupt?', ...
    'Yes','No','Nope');

switch ButtonName
case 'Yes'
    disp('going on');
    mystop = 0;
case 'No'
    mystop = 1;
end % switch

author: CelsoReyes


text:
function D=dist_mat(P1, P2)

n1 = zeros(1,size(P1,1));
n2 = zeros(1,size(P2,1));
for i=1:size(P1,1)
    n1(i) = norm(P1(i,:));
end
for i=1:size(P2,1)
    n2(i) = norm(P2(i,:));
end
n1 = n1.^2;
n2 = n2.^2;

D = -2*P1*P2';
N = repmat(n1',[1,size(P2,1)]) + repmat(n2,[size(P1,1),1]);
D = D+N;

D = sqrt(D);

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
main_fig = figure(...
'Units','characters',...
'Position',[135 4 198 56],...
'Visible',get(0,'defaultfigureVisible'),...
'Color',get(0,'defaultfigureColor'),...
'IntegerHandle','off',...
'MenuBar','none',...
'Name','CNMF_E',...
'NumberTitle','off',...
'Resize','on',...
'PaperPosition',get(0,'defaultfigurePaperPosition'),...
'ScreenPixelsPerInchMode','manual',...
'ParentMode','manual',...
'HandleVisibility','callback',...
'Tag','figure1' );

author: zhoupc


text:
function ir = shorten_ir(ir,nsamples)
%SHORTEN_IR shortens an impulse response
%
%   Usage: ir = shorten_ir(ir,nsamples)
%
%   Input parameters:
%       ir          - impulse response with length x channels
%       nsamples    - length of the target impulse response
%
%   Output paramteres:
%       ir          - impulse response signal with nsamples x n
%
%   SHORTEN_IR(ir,nsamples) shortens a given impulse response to the given
%   number of samples nsamples and applying a 5% long hanning window.
%
%   See also: get_ir, reduce_ir

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 2;
nargmax = 2;
narginchk(nargmin,nargmax);
isargpositivescalar(nsamples);

%% ===== Computation ====================================================
% Window impulse response
win = hann_window(0,ceil(0.05*nsamples),nsamples);
ir = ir(1:nsamples,:) .* repmat(win,[1 size(ir,2)]);

author: sfstoolbox

repo: sfs-matlab


text:
classdef prtFeatSelSfs < prtFeatSel
% prtFeatSelSfs   Sequential forward feature selection object.
%
%    FEATSEL = prtFeatSelSfs creates a sequental forward feature selection
%    object.
%
%    FEATSEL = prtFeatSelSfs(PROPERTY1, VALUE1, ...) constructs a
%    prttFeatSelSfs object FEATSEL with properties as specified by
%    PROPERTY/VALUE pair
%
%    A prtFeatSelSfsobject has the following properties:
%
%    nFeatures             - The number of features to be selected
%    showProgressBar       - Flag indicating whether or not to show the
%                            progress bar during feature selection.
%    evaluationMetric      - The metric to be used to determine which
%                            features are selected. evaluationMetric must
%                            be a function handle. The function handle must
%                            be in the form:
%                            @(dataSet)prtEval(prtClass, dataSet, varargin)
%                            where prtEvak is a prtEval function, prtClass
%                            is a prt classifier object, and varargin
%                            represents optional input arguments to a
%                            prtEval function.
%
%    peformance            - The performance obtained by the using the
%                            features selected.
%    selectedFeatures      - The indices of the features selected that gave
%                            the best performance.
%
%   A prtFeatSelExhaustive object inherits the TRAIN and RUN methods
%   from prtClass.
%
%   Example:
%
%   dataSet = prtDataGenFeatureSelection;         % Generate a data set
%   featSel = prtFeatSelSfs;          % Create a feature selction object
%   featSel.nFeatures = 3;            % Select only one feature of the data
%   featSel = featSel.train(dataSet); % Train the feature selection object
%   outDataSet = featSel.run(dataSet);% Extract the data set with only the
%                                     % selected features
%
%   %   Change the scoring function to prtScorePdAtPf, and change the
%   %   classification method to prtClassMAP
%
%   featSel.evaluationMetric = @(DS)prtEvalPdAtPf(prtClassMap, DS, .9);
%
%   featSel = featSel.train(dataSet);
%   outDataSet = featSel.run(dataSet);
%
 % See Also:  prtFeatSelStatic, prtFeatSelExhaustive

    properties (SetAccess=private)
        name = 'Sequentual Feature Selection' % Sequentual Feature Selection
        nameAbbreviation = 'SFS' % SFS
    end

    properties
        % General Classifier Properties
        nFeatures = 3;                    % The number of features to be selected
        evaluationMetric = @(DS)prtEvalAuc(prtClassFld,DS);   % The metric used to evaluate performance
    end

    properties (SetAccess = protected)
        performance = [];        % The evalutationMetric for the selected features
        selectedFeatures = [];   % The integer values of the selected features
    end

    methods
        function Obj = prtFeatSelSfs(varargin)
            Obj.isCrossValidateValid = false;
            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end

        function Obj = set.nFeatures(Obj,val)
            if ~prtUtilIsPositiveScalarInteger(val);
                error('prt:prtFeatSelSfs','nFeatures must be a positive scalar integer.');
            end
            Obj.nFeatures = val;
        end

        function Obj = set.evaluationMetric(Obj,val)
            assert(isa(val, 'function_handle') && nargin(val)>=1,'prt:prtFeatSelExhaustive','evaluationMetric must be a function handle that accepts one input argument.');
            Obj.evaluationMetric = val;
        end

    end
    methods (Access=protected,Hidden=true)

        function Obj = trainAction(Obj,DS)

            nFeatsTotal = DS.nFeatures;
            nSelectFeatures = min(nFeatsTotal,Obj.nFeatures);

            Obj.performance = nan(1,nSelectFeatures);
            Obj.selectedFeatures = nan(1,nSelectFeatures);

            sfsSelectedFeatures = [];

            if Obj.showProgressBar
                h = prtUtilProgressBar(0,'Feature Selection - SFS','autoClose',true);
            end

            for j = 1:nSelectFeatures

                if Obj.showProgressBar && j == 1
                    h2 = prtUtilProgressBar(0,sprintf('Selecting Feature %d',j),'autoClose',false);
                elseif Obj.showProgressBar
                    h2.titleStr = sprintf('Selecting Feature %d',j);
                    h2.update(0);
                end

                availableFeatures = setdiff(1:nFeatsTotal,sfsSelectedFeatures);
                cPerformance = nan(size(availableFeatures));
                for i = 1:length(availableFeatures)
                    currentFeatureSet = cat(2,sfsSelectedFeatures,availableFeatures(i));
                    tempDataSet = DS.retainFeatures(currentFeatureSet);

                    cPerformance(i) = Obj.evaluationMetric(tempDataSet);

                    if Obj.showProgressBar
                        h2.update(i/length(availableFeatures));
                    end
                end

                if Obj.showProgressBar
                    % Make sure it's closed.
                    h2.update(1);
                end

                if all(~isfinite(cPerformance))
                    error('prt:prtFeatSelSfs','All evaluation matrics resulted in non-finite values. Check evalutionMetric');
                end

                % Randomly choose the next feature if more than one provide the same performance
                val = max(cPerformance);
                newFeatInd = find(cPerformance == val);
                newFeatInd = newFeatInd(max(1,ceil(rand*length(newFeatInd))));

                % In the (degenerate) case when rand==0, set the index to the first one
                Obj.performance(j) = val;
                sfsSelectedFeatures = cat(2,sfsSelectedFeatures,availableFeatures(newFeatInd));

                if Obj.showProgressBar
                    h.update(j/nSelectFeatures);
                end
            end
            Obj.selectedFeatures = sfsSelectedFeatures;
            if Obj.showProgressBar
                % Make sure it's closed.
                h.update(1);
            end
        end

        function DataSet = runAction(Obj,DataSet)
            if ~Obj.isTrained
                error('prt:prtFeatSelSfs','Attempt to run a prtFeatSel that is not trained');
            end
            DataSet = DataSet.retainFeatures(Obj.selectedFeatures);
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function c=ref_gdgt(f,g,a,M,c_t,c_f,c_w)
%REF_GDGT  Reference generalized DGT
%   Usage:  c=ref_dgtiv(f,g,a,M,c_t,c_f,c_w);
%
%   Linear algebra version of the algorithm. Create big matrix
%   containing all the basis functions and multiply with the transpose.

L=size(f,1);

N=L/a;

F=zeros(L,M*N);

l=(0:L-1).';

for n=0:N-1
  for m=0:M-1
    F(:,M*n+m+1)=exp(2*pi*i*(m+c_f)*(l+c_t)/M).*circshift(g,n*a+c_w);
  end;
end;

c=F'*f;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function x = isnr(ref, sig, obs)
%
% snr -- Compute Improvement Signal-to-Noise Ratio for images
%
% Usage:
%       x = isnr(ref, sig, obs)
%
% Input:
%       ref         Reference image
%       sig         Modified image
%       obs         Observed image
%
% Output:
%       x           SNR value
%
% Authors
%   Paul Rodriguez    prodrig@pucp.edu.pe
%   Brendt Wohlberg   brendt@tmail.lanl.gov
%

mse1 = mean((ref(:)-sig(:)).^2);
mse2 = mean((obs(:)-sig(:)).^2);
x = 10*log10(mse2/mse1);

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
% Create Noise Signal and embed sine wave in it at a random location
npts=5000; % # points in signal
srate=1000; % sample rate (Hz)
dur=npts/srate; % signal duration in sec
amp_n=3.5; % noise amplitude
amp_t=1.0; % sine wave amplitude
freq=100; % sine wave frequency
sinepts=400; % # points in sine wave

t=[0:npts-1]/srate; % signal time base
sine_t=[0:sinepts-1]/srate; % sine wave time base
noise=(rand(1,npts)-.5)*amp_n; %noise signal
sinewave=amp_t*(sin(2*pi*freq*sine_t)); % sine wave
% random location of sinewave
sinepos=floor(rand(1,1)*(npts-sinepts))+1;
signal = noise; % make signal equal to noise
endsine = sinepos + sinepts-1; %calc index end of sine wave
% add sinewave to signal starting at index sinepos
signal(sinepos:endsine) = signal(sinepos:endsine) + sinewave;

% Plot signal
subplot(5,1,1)
plot(t,signal);
hold on
%plot red dot on location of sine wave start
plot(sinepos/srate,0,'.r');
hold off
title('Original Signal');

% Step 1: Filter signal
hbandwidth=5; %half bandwidth
nyquist=srate/2;
ntaps=200;
hpf=(freq-hbandwidth)/nyquist;
lpf=(freq+hbandwidth)/nyquist;
[b,a]=fir1(ntaps, [hpf lpf]); %calc filter coefficients
signal_f=filter(b,a,signal); %filter signal
subplot(5,1,2)
plot(t,signal_f); %plot filtered signal
hold on
%plot red dot on location of sine wave
plot(sinepos/srate,0,'.r');
hold off
title('Step 1. Filter Signal');

signal_r=abs(signal_f); %abs value of filtered signal
subplot(5,1,3)
plot(t,signal_r); %plot filtered signal
hold on
%plot red dot on location of sine wave
plot(sinepos/srate,0,'.r');
hold off
title('Step 2. Rectify Signal');

% Step 3: Low-pass filter rectified signal
lpf=10/nyquist; % low-pass filter corner frequency
npoles=6;
[b,a]=butter(npoles, lpf); % Butterworth filter coeff
signal_env=filter(b,a,signal_r); %Filter signal
subplot(5,1,4)
plot(t,signal_env); %plot filtered signal
hold on
%plot red dot on location of sine wave
plot(sinepos/srate,0,'.r');
hold off
title('Step 3. Envelope Signal');

% Step 4: Threshold signal
threshold=amp_t/2;
gated=(signal_env>threshold);
subplot(5,1,5)
plot(t,signal_env); %plot filtered signal
hold on
%plot red dot on location of sine wave
plot(sinepos/srate,0,'.r');
plot(t, gated, 'r'); % plot threshold signal
hold off
title('Step 4. Threshold Signal');
xlabel('Time (s)');

d_gated=diff(gated);
plot(d_gated);

find(d_gated==1)

find(d_gated == -1)

abs(find(d_gated==1) - find(d_gated == -1))

author: yandld


text:
function display(this);

    disp(inputname(1));
    if length(this) == 1
        struct(this.EEG)
    else
        TMP = struct(this);

        TMP2 = TMP(1).EEG;
        fieldorder = fieldnames(TMP2);
        for index = 2:length(TMP)
            TMP2(index) = orderfields(TMP(index).EEG, fieldorder);
        end;
        TMP2
    end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function f = luflop (L, U)						    %#ok
%LUFLOP given L and U, computes # of flops required to compute them
%
% Example:
% f = luflop (L, U)
%
% Given an LU factorization, compute how many flops took to compute it.  This
% is the same as (assuming U has a zero-free diagonal):
%
%   Lnz = full (sum (spones (L))) - 1 ;
%   Unz = full (sum (spones (U')))' - 1 ;
%   f = 2*Lnz*Unz + sum (Lnz) ;
%
% except that no extra workspace is allocated for spones (L) and spones (U).
% L and U must be sparse.
%
% Note: the above expression has a subtle undercount when exact numerical
% cancelation occurs.  Try [L,U,P] = lu (sparse (ones (10))) and then
% luflop (L,U).
%
% See also LU

% Copyright 1995-2007 by Timothy A. Davis.

help luflop
error ('luflop mexFunction not found!  Use umfpack_make to compile luflop.') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Test file for @chebfun/logical.m.

function pass = test_logical(pref)

if (nargin < 1)
    pref = chebfunpref();
end

% Generate a few random points to use as test values.
seedRNG(6178);
x = 2 * rand(100, 1) - 1;

% Check a few basic cases.
f = chebfun(@(x) sin(2*x), [-1 1], pref);
g = logical(f);
pass(1) = isequal(g.pointValues, [1 0 1].') && all(feval(g, x) == 1);

f = chebfun(@(x) exp(x), [-1 -0.5 0.5 1], pref);
g = logical(f);
pass(2) = all(g.pointValues == 1) && all(feval(g, x) == 1);

f = chebfun(@(x) 0*x, [-1 -0.5 0.5 1], pref);
g = logical(f);
pass(3) = all(g.pointValues == 0) && all(feval(g, x) == 0);

f = chebfun(@(x) 0*x + 1, [-1 -0.5 0.5 1], pref);
g = logical(f);
pass(4) = all(g.pointValues == 1) && all(feval(g, x) == 1);

% Check complex chebfun.
f = chebfun(@(x) exp(1i*x).*sin(x), [-1 -0.5 0.5 1], pref);
g = logical(f);
pass(5) = isequal(g.pointValues, [1 1 0 1 1].') && all(feval(g, x) == 1);

% Check array-valued chebfun.
f = chebfun(@(x) [sin(x) exp(x)], [-1 -0.5 0.5 1], pref);
g = logical(f);
pass(6) = isequal(g.pointValues, [1 1 0 1 1 ; 1 1 1 1 1].') && ...
    all(all(feval(g, x) == 1));

%% Test for singular function:

% define the domain:
dom = [-2 7];

op = @(x) sin(30*x)./((x-dom(1)).*(x-dom(2)));
f = chebfun(op, dom, 'exps', [-1 -1], 'splitting', 'on');
h = logical(f);

% check values:

% Generate a few random points to use as test values:
x = diff(dom) * rand(100, 1) + dom(1);

fval = feval(h, x);
err = fval - 1;
pass(7) = ~any( err );

r = roots(f);
pass(8) = ~any( h(r) );

%% Test for function defined on unbounded domain:

% Functions on [-inf b]:

% Set the domain:
dom = [-Inf 0 3*pi ];
domCheck = [-1e6 3*pi];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

% Blow-up function:
op = @(x) 5+exp(x.^3);
opExact = @(x) logical(op(x));
f = chebfun({op 0}, dom);
g = logical(f);
gVals = feval(g, x);
gExact = opExact(x);
err = gVals - gExact;
pass(9) = ~any(err);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function linpack_s_test07 ( )

%*****************************************************************************80
%
%% TEST07 tests SGBFA and SGBSL.
%
%  Discussion:
%
%    SGBFA and SGBSL are for general banded matrices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 100;
  ml = 25;
  mu = 25;
  lda = 2*ml+mu+1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST07\n' );
  fprintf ( 1, '  For a general banded matrix,\n' );
  fprintf ( 1, '  SGBFA factors the matrix,\n' );
  fprintf ( 1, '  SGBSL solves a factored linear system.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Assign values to matrix A and right hand side B.
%
%  We want to try a problem with a significant bandwidth.
%
  m = ml + mu + 1;
  fprintf ( 1, '  The bandwidth of the matrix is %d\n', m );

  for j = 1 : n

    ilo = max ( 1, j - mu );
    ihi = min ( n, j + ml );

    temp = 0.0;
    for i = ilo : ihi
      a(i-j+m,j) = -1.0;
      temp = temp - 1.0;
    end

    temp = temp + 1.0;
    a(m,j) = 4.0 - temp;
    b(j) = 4.0;

  end
%
%  Force B to be a column vector.
%
  b = b';
%
%  Factor the matrix A.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Factor the matrix.\n' );

  [ a, ipivot, info ] = sgbfa ( a, lda, n, ml, mu );

  if ( info ~= 0 )
   fprintf ( 1, '  Error!  SGBFA returns INFO = %d\n', info );
    return
  end
%
%  Call SGBSL to solve the linear system.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the linear system.\n' );

  job = 0;
  b = sgbsl ( a, lda, n, ml, mu, ipivot, b, job );
%
%  Print the results.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The first and last 5 entries of the solution:\n' );
  fprintf ( 1, '  (All should be 1):\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    if ( i <= 5 | n-5 < i )
      fprintf ( 1, '  %6d  %14f\n', i, b(i) );
    end
    if ( i == 5 )
      fprintf ( 1, '  ......  ..............\n' );
    end
  end

  return
end


text:
function [cIX,gIX] = GrowClustersFromSeedsItr(thres_merge,thres_cap,thres_minsize,thres_reg,cIX,gIX,M_0)
disp('find ROIs')

% Set params

% may be determined by hist of distances between all cells (in a sample)
% in correlation distance, i.e. 1 - corr.coeff
% thres_merge = 0.4;
% thres_cap = 0.5;

% thres_minsize = 10; % cell number in final clusters

%%
gIX = SqueezeGroupIX(gIX);
M = M_0(cIX,:);
C = FindCentroid_Direct(gIX,M); % kmeans 20x20, 'supervoxels'
nFoxels = size(C,1);

%% Calculate correlation distance between all cluster-centroids
temp = pdist(C,'correlation');
Dist = squareform(temp);
for i = 1:nFoxels,
    Dist(i,i) = NaN;
end

%% ITERATION:
% initialize
ROI = [];
ROI(nFoxels).fxlist = []; % not capping
ROI(nFoxels).numcell = [];
ROIcount = 0;

% analogy ~ mask to store new ROI
BW = zeros(nFoxels,1);

%%
tic
while true, % loop through all qualified seeds
    % look for next seed
    [~,ix] = min(Dist(:));
    [I,J] = ind2sub(size(Dist),ix);
    % and add seed position to mask
    BW(I) = 1;
    BW(J) = 1;

    if Dist(I,J)>(1-thres_merge),
        break;
    end

    while true % grow ROI: expand mask and examine next neighbors
        IX_in = find(BW);
        IX_out = find(BW==0);

        % find next closest cluster
        D = 1-corr(C(IX_in,:)',C(IX_out,:)'); % distance matrix between new core and the original cluster means
        D2 = min(D,[],1);
        [a,ix] = min(D2);
        ix_pretend = IX_out(ix);

        BW_pretend = BW;
        BW_pretend(ix_pretend) = 1;

        % find correlation between potential cluster and new core
        list = find(BW_pretend);
        IX = [];
        for i_gIX = 1:length(list),
            IX = [IX;find(gIX==list(i_gIX))];
        end
        M_core = M(IX,:);
        [~,newCore] = kmeans(M_core,1,'distance','correlation');
        coredist = 1-corr(C(ix_pretend,:)',newCore');

        if a<(1-thres_merge) && coredist<(1-thres_cap),
            BW(ix_pretend) = 1;
        else
            break; % finished expanding this seed
        end
    end

    % save this ROI if bigger than thres
    list = find(BW);
    IX = [];
    for i_gIX = 1:length(list),
        IX = [IX;find(gIX==list(i_gIX))]; %#ok<AGROW>
    end
    numcell = length(IX);
    if numcell >= thres_minsize,
        ROIcount = ROIcount+1;
        ROI(ROIcount).fxlist = find(BW);
        ROI(ROIcount).numcell = numcell;
    end

    % update/reset
    Dist(list,:) = NaN;
    Dist(:,list) = NaN;
    C(list,:) = NaN;
    BW = zeros(nFoxels,1);
end

ROI(ROIcount+1:end) = [];
toc
%% Merge accordingly
for i = 1:ROIcount,
    fxlist = ROI(i).fxlist;
    for j = 2:length(fxlist),
        gIX(gIX==fxlist(j)) = fxlist(1);
    end
end
[gIX,numU] = SqueezeGroupIX(gIX);
disp(numU);

%% Regression with the centroid of each cluster, round 2
% disp('auto-reg');
% Reg = FindCentroid_Direct(gIX,M);
% [cIX,gIX] = AllCentroidRegression_direct(M_0,thres_reg,Reg);

%% size threshold
U = unique(gIX);
numU = length(U);
for i=1:numU,
    if length(find(gIX==U(i)))<thres_minsize,
        cIX(gIX==U(i)) = [];
        gIX(gIX==U(i)) = [];
    end
end
[gIX,numU] = SqueezeGroupIX(gIX);
disp(numU);
end

author: xiuyechen

repo: FishExplorer


text:
% LOOK FOR SINGULAR POINTS IN JOINT SPACE

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.

function singular_points = find_singular_points(param)

singular_points=[];

q=[0 0 0 0];

teta = eval(param.DH.theta);
d = eval(param.DH.d);
a = eval(param.DH.a);
alfa = eval(param.DH.alpha);

%J = eval(param.J)

q1=-pi:0.05:pi;
q2=-pi:0.05:pi;

for i=1:length(q1),
    fprintf('\nIteration %d of %d', i, length(q1))
   for j=1:length(q2),
        q=[q1(i) q2(j) 0];
        J = eval(param.J);

        val = det(J);
        if abs(val) < 0.0001
            singular_points = [singular_points q'];
        end
   end
end

repo: ARTE


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nickabattista@gmail.com) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Models the force interaction that is user-defined.
%
%           Note:
%                   1. User has complete control for artibtrary input file
%                      format and number of parameters
%                   2. Those parameters get passed to the matrix called
%                      "general_force"
%                   3. User gets options of including other things into
%                      force function outside of inputted parameters, e.g.,
%                      current_time, dt, previous location of Lag. Pts, etc.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [fx_genForce, fy_genForce] = give_Me_General_User_Defined_Force_Densities(ds,Nb,xLag,yLag,xLag_P,yLag_P,dt,current_time,general_force)

    %
    % INPUTS:
    %         ds: Lagragian spacing (defined by ds = 0.5*dx)
    %         Nb: # of Lagrangian Pts.
    %         xLag: current x-Lagrangian coordinate positions
    %         yLag: current y-Lagrangian coordinate positions
    %         xLag_P: previous x-Lagrangian coordinate positions
    %         yLag_P: previous y-Lagrangian coordinate positions
    %         dt: time-step value
    %         current_time: current time in simulation
    %         general_force: matrix containing all data from input file
    %

    %
    % NOTE: THIS EXAMPLE CREATES A USER-DEFINED FORCE THAT IS JUST AN
    % ORDINARY LINEAR SPRING
    %

    Nsprings = length(general_force(:,1));  % # of Springs
    sp_1 = general_force(:,1);              % Initialize storage for MASTER NODE Spring Connection
    sp_2 = general_force(:,2);              % Initialize storage for SLAVE NODE Spring Connection
    K_Vec = general_force(:,3);             % Stores spring stiffness associated with each spring
    RL_Vec = general_force(:,4);            % Stores spring resting length associated with each spring

    fx = zeros(Nb,1);                 % Initialize storage for x-forces
    fy = fx;                          % Initialize storage for y-forces

    %
    % Loops over all the master-nodes of the springs to compute forces
    %
    for i=1:Nsprings

        id_Master = sp_1(i);          % Master Node index
        id_Slave = sp_2(i);           % Slave Node index
        k_Spring = K_Vec(i);          % Spring stiffness of i-th spring
        L_r = RL_Vec(i);              % Resting length of i-th spring

        dx = xLag(id_Slave) - xLag(id_Master); % x-Distance btwn slave and master node
        dy = yLag(id_Slave) - yLag(id_Master); % y-Distance btwn slave and master node

        sF_x =  k_Spring * ( sqrt( dx^2 + dy^2 ) - L_r ) * ( dx / sqrt(dx^2+dy^2) ); % Compute x-Force
        sF_y =  k_Spring * ( sqrt( dx^2 + dy^2 ) - L_r ) * ( dy / sqrt(dx^2+dy^2) ); % Compute y-Force

        fx(id_Master,1) = fx(id_Master,1) + sF_x;  % Sum total forces for node, i in x-direction (this is MASTER node for this spring)
        fy(id_Master,1) = fy(id_Master,1) + sF_y;  % Sum total forces for node, i in y-direction (this is MASTER node for this spring)

        fx(id_Slave,1) = fx(id_Slave,1) - sF_x;    % Sum total forces for node, i in x-direction (this is SLAVE node for this spring)
        fy(id_Slave,1) = fy(id_Slave,1) - sF_y;    % Sum total forces for node, i in y-direction (this is SLAVE node for this spring)

    end

    % Store forces from artibrary force function onto desired Lagrangian Pts.
    fx_genForce = fx;
    fy_genForce = fy;

author: nickabattista