text:
function rrd = jplephem (et, ntarg, ncent)

% reads the jpl planetary ephemeris and gives
% the position and velocity of the point 'ntarg'
% with respect to point 'ncent'

% input

%   et    = julian ephemeris date at which interpolation is wanted

%   ntarg = integer number of 'target' point

%   ncent = integer number of center point

%   the numbering convention for 'ntarg' and 'ncent' is:

%        1 = mercury           8 = neptune
%        2 = venus             9 = pluto
%        3 = earth            10 = moon
%        4 = mars             11 = sun
%        5 = jupiter          12 = solar-system barycenter
%        6 = saturn           13 = earth-moon barycenter
%        7 = uranus           14 = nutations (longitude and obliq)
%                             15 = librations, if on ephemeris file

%        if nutations are wanted, set ntarg = 14.
%        for librations, set ntarg = 15. set ncent = 0.

% output

%   rrd = output 6-word array containing position and velocity
%         of point 'ntarg' relative to 'ncent'. the units are au and
%         au/day. for librations the units are radians and radians
%         per day. in the case of nutations the first four words of
%         rrd will be set to nutations and rates, having units of
%         radians and radians/day.

%         the option is available to have the units in km and km/sec.
%         for this, set km = 1 via global in the calling program.

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global cval ss au emrat ncon ipt np nv twot pc vc

global iephem ephname bary pvsun nrl fid lpt

rrd = zeros(6, 1);

list = zeros(12, 1);

et2 = zeros(2, 1);

% load time array

et2(1) = et;

et2(2) = 0;

% first entry?

if (iephem == 1)

    pvsun = zeros(6, 1);

    % read header file data

    fid = fopen(ephname, 'r');

    ttl = fread(fid, 252);

    cnam = fread(fid, 2400);

    ss = fread(fid, 3, 'double');

    ncon = fread(fid, 1, 'int');

    % astronomical unit

    au = fread(fid, 1, 'double');

    % earth-moon ratio

    emrat = fread(fid, 1, 'double');

    ipt = fread(fid, [3 12], 'int');

    numde = fread(fid, 1, 'int');

    lpt = fread(fid, 3, 'int');

    % move to next record

    status = fseek(fid, 8144, 'bof');

    % read "constant" values

    cval = fread(fid, 400, 'double');

    % initialization

    nrl = 0;

    bary = 0;

    pc(1) = 1;
    pc(2) = 0;

    vc(2) = 1;

    np = 2;
    nv = 3;

    twot = 0;

    iephem = 0;
end

if (ntarg == ncent)

    return;

end

%%%%%%%%%%%%%%%%%%%%%%
% nutations requested?
%%%%%%%%%%%%%%%%%%%%%%

if (ntarg == 14)

    if (ipt(2, 12) > 0)

        list(11) = 2;

        [pv, rrd] = state(et2, list);

        list(11) = 0;

        return;

    else

        fprintf('\n\njplephem - no nutations on this ephemeris file \n');

        return;

    end
end

%%%%%%%%%%%%%%%%%%%%%%%
% librations requested?
%%%%%%%%%%%%%%%%%%%%%%%

if (ntarg == 15)

    if (lpt(2) > 0)

        list(12) = 2;

        [pv, rrd] = state(et2, list);

        list(12) = 0;

        for i = 1:1:6
            rrd(i) = pv(i, 11);
        end

        return

    else

        fprintf('\n\n no librations on this ephemeris file \n');

        return;

    end

end

% force barycentric output by function 'state'

bsave = bary;

bary = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set up proper entries in 'list' array for state call
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:1:2

    k = ntarg;

    if (i == 2)
        k = ncent;
    end

    if (k <= 10)
        list(k) = 2;
    end

    if (k == 10)
        list(3) = 2;
    end

    if (k == 3)
        list(10) = 2;
    end

    if (k == 13)
        list(3) = 2;
    end

end

%%%%%%%%%%%%%%%%%%%%
% make call to state
%%%%%%%%%%%%%%%%%%%%

[pv, rrd] = state(et2, list);

if (ntarg == 11 || ncent == 11)
    for i = 1:1:6
        pv(i, 11) = pvsun(i);
    end
end

if (ntarg == 12 || ncent == 12)
    for i = 1:1:6
        pv(i, 12) = 0;
    end
end

if (ntarg == 13 || ncent == 13)
    for i = 1:1:6
        pv(i, 13) = pv(i, 3);
    end
end

if (ntarg * ncent == 30 && ntarg + ncent == 13)
    for i = 1:1:6
        pv(i, 3) = 0;
    end
else
    if (list(3) == 2)
        for i = 1:1:6
            pv(i, 3) = pv(i, 3) - pv(i, 10) / (1 + emrat);
        end
    end

    if (list(10) == 2)
        for i = 1:1:6
            pv(i, 10) = pv(i, 3) + pv(i, 10);
        end
    end
end

for i = 1:1:6
    rrd(i) = pv(i, ntarg) - pv(i, ncent);
end

bary = bsave;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [pv, nut] = state(et2, list)

% reads and interpolates the jpl planetary ephemeris file

% input

%   et2    2-word julian ephemeris epoch at which interpolation
%          is wanted.  any combination of et2(1) + et2(2) which falls
%          within the time span on the file is a permissible epoch.

%          a. for ease in programming, the user may put the
%             entire epoch in et2(1) and set et2(2) = 0.

%          b. for maximum interpolation accuracy, set et2(1) equal
%             to the most recent midnight at or before interpolation
%             epoch and set et2(2) equal to fractional part of a day
%             elapsed between et2(1) and epoch.

%          c. as an alternative, it may prove convenient to set
%             et2(1) = some fixed epoch, such as start of integration,
%             and et2(2) = elapsed interval between and epoch.

%   list   12-word integer array specifying what interpolation
%          is wanted for each of the bodies on the file.
%
%          list(i) = 0 => no interpolation for body i
%                  = 1 => position only
%                  = 2 => position and velocity

%          the designation of the astronomical bodies by i is:

%              i =  1 => mercury
%                =  2 => venus
%                =  3 => earth-moon barycenter
%                =  4 => mars
%                =  5 => jupiter
%                =  6 => saturn
%                =  7 => uranus
%                =  8 => neptune
%                =  9 => pluto
%                = 10 => geocentric moon
%                = 11 => nutations in longitude and obliquity
%                = 12 => lunar librations (if on file)

% output

%   pv   6 x 11 array that will contain requested interpolated
%        quantities. the body specified by list(i) will have its
%        state in the array starting at pv(1, i). (on any given
%        call, only those words in 'pv' which are affected by the
%        first 10 'list' entries (and by list(12) if librations are
%        on the file) are set. the rest of the 'pv' array
%        is untouched).  the order of components starting in
%        pv(1, i) is x, y, z, dx, dy, dz.

%        all output vectors are referenced to the earth mean
%        equator and equinox of j2000 if the de number is 200 or
%        greater; of b1950 if the de number is less than 200.

%        the moon state is always geocentric; the other nine states
%        are either heliocentric or solar-system barycentric,
%        depending on the setting of common flags (see below).

%        lunar librations, if on file, are put into pv(k, 11) if
%        list(12) is 1 or 2.
%
%   nut  4-word array that will contain nutations and rates,
%        depending on the setting of list(11). the order of
%        quantities in nut is:

%        d psi  (nutation in longitude)
%        d epsilon (nutation in obliquity)
%        d psi dot
%        d epsilon dot

% global

%   km    logical flag defining physical units of the output states
%         = 1 => kilometers and kilometers/second
%         = 0 => au and au/day

%         default value = 0 (km determines time unit
%         for nutations and librations. angle unit is always radians.)

%   bary  logical flag defining output center.
%         only the 9 planets are affected.
%         bary = 1 => center is solar-system barycenter
%              = 0 => center is sun
%         default value = 0

%   pvsun 6-word array containing the barycentric position and
%         velocity of the sun

global ss au ipt lpt

global nrl fid km bary pvsun coef

nut = zeros(4, 1);

pv = zeros(6, 11);

if (et2(1) == 0)
    return
end

s = et2(1) - 0.5;

tmp = split(s);

pjd(1) = tmp(1);

pjd(2) = tmp(2);

tmp = split(et2(2));

pjd(3) = tmp(1);

pjd(4) = tmp(2);

pjd(1) = pjd(1) + pjd(3) + 0.5;

pjd(2) = pjd(2) + pjd(4);

tmp = split(pjd(2));

pjd(3) = tmp(1);

pjd(4) = tmp(2);

pjd(1) = pjd(1) + pjd(3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% error return for epoch out of range
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (pjd(1) + pjd(4) < ss(1) || pjd(1) + pjd(4) > ss(2))
    fprintf('\n\n error in state - epoch out of range \n');

    return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate record number and relative time in interval
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nr = fix((pjd(1) - ss(1)) / ss(3)) + 2;

if (pjd(1) == ss(2))
    nr = nr - 1;
end

t(1) = ((pjd(1) - ((nr-2) * ss(3) + ss(1))) + pjd(4)) / ss(3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read correct record if not in core
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (nr ~= nrl)
    nrl = nr;

    status = fseek(fid, nr * 8144, 'bof');

    coef = fread(fid, 1018, 'double');
end

if (km == 1)
    t(2) = ss(3) * 86400;
    aufac = 1;
else
    t(2) = ss(3);
    aufac = 1 / au;
end

% interpolate barycentric state vector of sun

tmpv = zeros(3, 2);

ibuf = ipt(1, 11);
ncf = ipt(2, 11);
na = ipt(3, 11);

tmpv = interp(ibuf, t, ncf, 3, na, 2);

k = 0;

for j = 1:1:2
    for i = 1:1:3
        k = k + 1;
        pvsun(k) = tmpv(i, j) * aufac;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check and interpolate bodies requested
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tmpv = zeros(3, 2);

tempv = zeros(6, 1);

for i = 1:1:10
    if (list(i) ~= 0)

        ibuf = ipt(1, i);
        ncf = ipt(2, i);
        na = ipt(3, i);

        tmpv = interp(ibuf, t, ncf, 3, na, list(i));

        k = 0;

        for j = 1:1:2
            for m = 1:1:3
                k = k + 1;
                tempv(k) = tmpv(m, j);
            end
        end

        for j = 1:1:6
            if (i <= 9 && bary == 0)
                pv(j, i) = tempv(j) * aufac - pvsun(j);
            else
                pv(j, i) = tempv(j) * aufac;
            end
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do nutations if requested (and if on file)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (list(11) > 0 && ipt(2, 12) > 0)

    ibuf = ipt(1, 12);
    ncf = ipt(2, 12);
    na = ipt(3, 12);

    nut = interp(ibuf, t, ncf, 2, na, list(11));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get librations if requested (and if on file)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (list(12) > 0 && lpt(2) > 0)
    tmpv = interp(lpt(1), t, lpt(2), 3, lpt(3), list(12));

    for i = 1:1:3
        pv(i, 11) = tmpv(i, 1);

        pv(i + 3, 11) = tmpv(i, 2);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function pv = interp(ibuf, t, ncf, ncm, na, ifl)

% this function differentiates and interpolates a set
% of chebyshev coefficients to give position and velocity

% input

%   buf   1st location of array of chebyshev coefficients of position

%   t     t(1) is fractional time in interval covered by
%         coefficients at which interpolation is wanted
%         (0 <= t(1) <= 1). t(2) is length of whole
%         interval in input time units.

%   ncf   number of coefficients per component

%   ncm   number of components per set of coefficients

%   na    number of sets of coefficients in full array
%         (i.e., number of sub-intervals in full interval)

%   ifl   integer flag
%         = 1 for positions only
%         = 2 for pos and vel

% output

%   pv    interpolated quantities requested.  dimension
%         expected is pv(ncm, ifl)

global coef np nv twot pc vc

pv = zeros(6, 12);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get correct sub-interval number for this set of
% coefficients and get normalized chebyshev time
% within that subinterval.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dna = na;

dt1 = fix(t(1));

temp = dna * t(1);

ll = fix(temp - dt1) + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% tc is the normalized chebyshev time (-1 <= tc <= 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tc = 2 * (mod(temp, 1) + dt1) - 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check to see whether chebyshev time has changed,
% and compute new polynomial values if it has.
% (the element pc(2) is the value of t1(tc) and hence
% contains the value of tc on the previous call)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (tc ~= pc(2))
    np = 2;
    nv = 3;
    pc(2) = tc;
    twot = tc + tc;
end

% be sure that at least 'ncf' polynomials have been evaluated
% and are stored in the array 'pc'.

if (np < ncf)
    for i = np + 1:1:ncf
        pc(i) = twot * pc(i - 1) - pc(i - 2);
    end

    np = ncf;
end

bcoef = ncf * na * ncm;

cbody = zeros(bcoef, 1);

n = ibuf;

for m = 1:1:bcoef
    cbody(m) = coef(n);
    n = n + 1;
end

cbuf = zeros(ncf, ncm, na);

n = 0;

for l = 1:1:na
    for i = 1:1:ncm
        for j = 1:1:ncf
            n = n + 1;
            cbuf(j, i, l) = cbody(n);
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% interpolate to get position for each component
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:1:ncm
    pv(i, 1) = 0;

    for j = ncf:-1:1
        pv(i, 1) = pv(i, 1) + pc(j) * cbuf(j, i, ll);
    end
end

if (ifl <= 1)
    return
end

% if velocity interpolation is wanted, be sure enough
% derivative polynomials have been generated and stored.

vfac = (dna + dna) / t(2);

vc(3) = twot + twot;

if (nv < ncf)
    for i = nv + 1:1:ncf
        vc(i) = twot * vc(i - 1) + pc(i - 1) + pc(i - 1) - vc(i - 2);
    end

    nv = ncf;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% interpolate to get velocity for each component
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:1:ncm
    pv(i, 2) = 0;

    for j = ncf:-1:2
        pv(i, 2) = pv(i, 2) + vc(j) * cbuf(j, i, ll);
    end

    pv(i, 2) = pv(i, 2) * vfac;
end

%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

function fr = split(tt)

% this function breaks a number into a integer
% and a fractional part.

% input

%   tt = input number

% output

%   fr = 2-word output array
%        fr(1) contains integer part
%        fr(2) contains fractional part

%        for negative input numbers, fr(1) contains the next
%        more negative integer; fr(2) contains a positive fraction.

fr = zeros(2, 1);

fr(1) = fix(tt);

fr(2) = tt - fr(1);

if (tt >= 0 || fr(2) == 0)
    return
end

% make adjustments for negative input number

fr(1) = fr(1) - 1;

fr(2) = fr(2) + 1;

author: Sable


text:
function varargout = plot(varargin)
% plot method for GIfTI objects
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: plot.m 5888 2014-02-19 19:54:12Z guillaume $

% if ishandle(varargin{1})
%     h = figure(varargin{1});
% else
%     h = figure;
%     %axis equal;
%     %axis off;
%     %camlight;
%     %camlight(-80,-10);
%     %lighting phong;
% end
% cameramenu;

cdata = [];
ax = [];
if nargin == 1
    this = varargin{1};
    h = gcf;
else
    if ishandle(varargin{1})
        ax = varargin{1};
        h = figure(get(ax,'parent'));
        this = varargin{2};
    else
        this = varargin{1};
        h = gcf;
        cdata = subsref(varargin{2},struct('type','.','subs','cdata'));
    end
    if nargin > 2
        indc = varargin{3};
    else
        indc = 1;
    end
end

if isempty(ax), ax = axes('Parent',h); end
axis(ax,'equal');
axis(ax,'off');
hp = patch(struct(...
    'vertices',  subsref(this,struct('type','.','subs','vertices')),...
    'faces',     subsref(this,struct('type','.','subs','faces'))),...
    'FaceColor', 'b',...
    'EdgeColor', 'none',...
    'Parent',ax);

if ~isempty(cdata)
    set(hp,'FaceVertexCData',cdata(:,indc), 'FaceColor','interp')
end

axes(ax);
camlight;
camlight(-80,-10);
lighting phong;
axes(ax);
cameramenu;

if nargout
    varargout{1} = hp;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function A = create_3d_cubic_interp_matrix(T, varargin)
% A = create_3d_bilinear_interp_matrix(T, [szin])
% T should contain displacements in pixels
% size(T) = [sz1, sz2, sz3, 3]
%
% x_def = A * x(:);

sz = [size(T, 1), size(T, 2), size(T, 3)];
if numel(varargin) == 1
    szin = varargin{1};
else
    szin = sz;
end
npixin = szin(1) * szin(2) * szin(3);
npixout = sz(1) * sz(2) * sz(3);

[m1, m2, m3] = ndgrid(1:sz(1), 1:sz(2), 1:sz(3));
TT = T + cat(4, m1, m2, m3);
rowidx = floor(TT(:, :, :, 1));
colidx = floor(TT(:, :, :, 2));
dolidx = floor(TT(:, :, :, 3));
rowk = TT(:, :, :, 1) - rowidx;
colk = TT(:, :, :, 2) - colidx;
dolk = TT(:, :, :, 3) - dolidx;

Kr = koefs(rowk);
Kc = koefs(colk);
Kd = koefs(dolk);

for i = [-1, 0, 1, 2]
    for j = [-1, 0, 1, 2]
        for k = [-1, 0, 1, 2]
%             [i,j,k]
            r = rowidx(:) + i;
            c = colidx(:) + j;
            d = dolidx(:) + k;
            s = Kr{i+2} .* Kc{j+2} .* Kd{k+2};
            t = 1 : npixout;
            idx = (r <= szin(1)) & (r >= 1) & (c <= szin(2)) & (c >= 1) & (d >= 1) & (d <= szin(3));

            r = r(idx);
            c = c(idx);
            d = d(idx);
            s = s(idx);
            t = t(idx);

            sp_c_idx = r + szin(1) * (c - 1) + szin(1)*szin(2)* (d - 1);

            if i == -1 && j == -1 && k == -1
                A = sparse(t, sp_c_idx, s, npixout, npixin);
            else
                tmp = sparse(t, sp_c_idx, s, npixout, npixin);
                A = A + tmp;
            end
        end
    end
end
A = A/8;
end

function B = koefs(A)
    B{1} = A.*((2-A).*A-1);
    B{2} = A.^2.*(3*A-5)+2;
    B{3} = A.*((4-3*A).*A+1);
    B{4} = A.^2.*(A-1);
end

repo: pTVreg


text:
function [initXOptions, optionsAll] = hsvargplvmInitXOptions(Ytr, options, globalOpt)

stackedOpt = globalOpt.stackedInitOpt;

%--- Here we have the option of using Bayesian GPLVM or GPLVM for
% initialising the latent spaces. If this is the case, train the
% corresponding models
optionsAll = hsvargplvmCreateOptions(Ytr, options, globalOpt);
initXOptions = cell(1, options.H);
for h=1:options.H
    if strcmp(optionsAll.initX{h}, 'vargplvm') | strcmp(optionsAll.initX{h}, 'fgplvm')
        initXOptions{h}{1} = optionsAll;
        % DOn't allow the D >> N trick for layers > 1
        if h~=1
            if isfield(initXOptions{h}{1}, 'enableDgtN')
                initXOptions{h}{1}.enableDgtN = false;
            end
        end
        initXOptions{h}{1}.latentDim = optionsAll.Q{h};
        initXOptions{h}{1}.numActive = optionsAll.K{h}{1};
        initXOptions{h}{1}.kern = optionsAll.kern{h}{1};
        initXOptions{h}{1}.initX = 'ppca';
        initXOptions{h}{1}.initSNR = 90;
        initXOptions{h}{1}.numActive = 50;
        initXOptions{h}{2} = 160;
        initXOptions{h}{3} = 30;
        if ~isempty(stackedOpt)
            if isfield(stackedOpt, 'stackedInitVardistIters') && ~isempty(stackedOpt.stackedInitVardistIters)
                initXOptions{h}{2} = stackedOpt.stackedInitVardistIters;
            end
            if isfield(stackedOpt, 'stackedInitIters') && ~isempty(stackedOpt.stackedInitIters)
                initXOptions{h}{3} = stackedOpt.stackedInitIters;
            end
            if isfield(stackedOpt, 'stackedInitSNR') && ~isempty(stackedOpt.stackedInitSNR)
                initXOptions{h}{1}.initSNR = stackedOpt.stackedInitSNR;
            end
            if isfield(stackedOpt, 'stackedInitK') && ~isempty(stackedOpt.stackedInitK)
                initXOptions{h}{1}.numActive = stackedInitK;
            end
        end
    elseif ~isempty(stackedOpt) && (iscell(stackedOpt) && ~isempty(stackedOpt{h}))
        initXOptions{h} = stackedOpt{h};
    else
        initXOptions{h} = {};
    end
end

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
function surf = imJointSurfaceArea(img, L1, L2, varargin)
% Surface area of the interface between two labels.
%
%   S = imJointSurfaceArea(LBL, L1, L2)
%   Estimates the joint surface area between the two labels L1 and L2 in
%   the label image LBL.
%
%   S = imJointSurfaceArea(LBL, L1, L2, NDIRS)
%   Specifies the number of directions used for estimating surface area.
%   NDIRS can be either 3 or 13 (the default).
%
%   S = imJointSurfaceArea(..., RESOL)
%   Specifies image resolution. RESOL is a 1-by-3 row vector containing
%   resolution in the X, Y and Z direction (in that order).
%
%
%   Example
%     % generate a demo image
%     img = discreteBall(1:10, 1:100, 1:100, [50.12 50.23 50.34 40]);
%     % convert to image with two different labels
%     img2 = uint8(img + 1);
%     % compute joint surface area
%     imJointSurfaceArea(img2, 1, 2)
%     ans =
%         2.0102e+004
%
%   See also
%     imSurfaceArea

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-07-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% check image dimension and type
if ndims(img) ~= 3 || islogical(img)
    error('first argument should be a 3D image');
end

%% Process input arguments

% default number of directions
ndir = 13;

% default image resolution
delta = [1 1 1];

% Process user input arguments
while ~isempty(varargin)
    var = varargin{1};
    if ~isnumeric(var)
        error('option should be numeric');
    end

    % option is either connectivity or resolution
    if isscalar(var)
        ndir = var;
    else
        delta = var;
    end
    varargin(1) = [];
end

%% Initialisations

% distances between a pixel and its neighbours.
d1  = delta(1);
d2  = delta(2);
d3  = delta(3);

% volume of a voxel (used for computing line densities)
vol = d1 * d2 * d3;

%% Main processing for 3 directions

% number of transitions along the 3 main directions
n1a = sum(sum(sum(img(:,1:end-1,:)==L1 & img(:,2:end,:)==L2)));
n1b = sum(sum(sum(img(:,1:end-1,:)==L2 & img(:,2:end,:)==L1)));
n2a = sum(sum(sum(img(1:end-1,:,:)==L1 & img(2:end,:,:)==L2)));
n2b = sum(sum(sum(img(1:end-1,:,:)==L2 & img(2:end,:,:)==L1)));
n3a = sum(sum(sum(img(:,:,1:end-1)==L1 & img(:,:,2:end)==L2)));
n3b = sum(sum(sum(img(:,:,1:end-1)==L2 & img(:,:,2:end)==L1)));

if ndir == 3
    % compute surface area by averaging over the 3 main directions
    surf = 4/3 * ((n1a+n1b)/d1 + (n2a+n2b)/d2 + (n3a+n3b)/d3) / 2 * vol;
    return;
end

%% Additional processing for 13 directions

% Number of connected components along diagonals contained in the three
% main planes
n4a = sum(sum(sum(img(2:end,1:end-1,:)==L1   & img(1:end-1,2:end,:)==L2)));
n4b = sum(sum(sum(img(2:end,1:end-1,:)==L2   & img(1:end-1,2:end,:)==L1)));
n5a = sum(sum(sum(img(1:end-1,1:end-1,:)==L1 & img(2:end,2:end,:)==L2)));
n5b = sum(sum(sum(img(1:end-1,1:end-1,:)==L2 & img(2:end,2:end,:)==L1)));
n6a = sum(sum(sum(img(:,2:end,1:end-1)==L1   & img(:,1:end-1,2:end)==L2)));
n6b = sum(sum(sum(img(:,2:end,1:end-1)==L2   & img(:,1:end-1,2:end)==L1)));
n7a = sum(sum(sum(img(:,1:end-1,1:end-1)==L1 & img(:,2:end,2:end)==L2)));
n7b = sum(sum(sum(img(:,1:end-1,1:end-1)==L2 & img(:,2:end,2:end)==L1)));
n8a = sum(sum(sum(img(2:end,:,1:end-1)==L1   & img(1:end-1,:,2:end)==L2)));
n8b = sum(sum(sum(img(2:end,:,1:end-1)==L2   & img(1:end-1,:,2:end)==L1)));
n9a = sum(sum(sum(img(1:end-1,:,1:end-1)==L1 & img(2:end,:,2:end)==L2)));
n9b = sum(sum(sum(img(1:end-1,:,1:end-1)==L2 & img(2:end,:,2:end)==L1)));

%TODO: add the case of 9 directions ?

% Number of connected components along lines corresponding to diagonals of
% the unit cube
n10a = sum(sum(sum(img(1:end-1,1:end-1,1:end-1)==L1 & img(2:end,2:end,2:end)==L2)));
n10b = sum(sum(sum(img(1:end-1,1:end-1,1:end-1)==L2 & img(2:end,2:end,2:end)==L1)));
n11a = sum(sum(sum(img(2:end,1:end-1,1:end-1)==L1 & img(1:end-1,2:end,2:end)==L2)));
n11b = sum(sum(sum(img(2:end,1:end-1,1:end-1)==L2 & img(1:end-1,2:end,2:end)==L1)));
n12a = sum(sum(sum(img(1:end-1,2:end,1:end-1)==L1 & img(2:end,1:end-1,2:end)==L2)));
n12b = sum(sum(sum(img(1:end-1,2:end,1:end-1)==L2 & img(2:end,1:end-1,2:end)==L1)));
n13a = sum(sum(sum(img(2:end,2:end,1:end-1)==L1 & img(1:end-1,1:end-1,2:end)==L2)));
n13b = sum(sum(sum(img(2:end,2:end,1:end-1)==L2 & img(1:end-1,1:end-1,2:end)==L1)));

% space between 2 voxels in each direction
d12  = hypot(d1, d2);
d13  = hypot(d1, d3);
d23  = hypot(d2, d3);
d123 = sqrt(d1^2 + d2^2 + d3^2);

% Compute weights corresponding to surface fraction of spherical caps
% For isotropic case, weights correspond to:
% c1 = 0.04577789120476 * 2;  % Ox
% c2 = 0.04577789120476 * 2;  % Oy
% c3 = 0.04577789120476 * 2;  % Oz
% c4 = 0.03698062787608 * 2;  % Oxy
% c6 = 0.03698062787608 * 2;  % Oxz
% c8 = 0.03698062787608 * 2;  % Oyz
% c10 = 0.03519563978232 * 2;  % Oxyz
c = computeDirectionWeights3d13(delta);

% compute the weighted sum of each direction
% intersection count * direction weight / line density
surf = 4 * vol * (...
    (n1a+n1b)*c(1)/d1 + (n2a+n2b)*c(2)/d2 + (n3a+n3b)*c(3)/d3 + ...
    (n4a+n4b+n5a+n5b)*c(4)/d12 + (n6a+n6b+n7a+n7b)*c(6)/d13 + ...
    (n8a+n8b+n9a+n9b)*c(8)/d23 + ...
    (n10a + n10b + n11a + n11b + n12a + n12b + n13a + n13b) * c(10) / d123) / 2;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function covMa = cal_covMa4Block(Block,blockWin)
K = size(Block,1);
num = size(Block,4);
N = size(Block,2);
K_m = K/2+1;
X = zeros(size(Block));
for i = 1:num
    X(:,:,:,i) = fft(Block(:,:,i).*blockWin);
end
covMa = zeros(N,N,K_m,num);
for i = 1:num
    for j = 1:K_m
        temp = permute(X(j,:,:,i),[2 3 1 4]);
        covMa(:,:,j,i) = temp*temp'/size(X,3);
    end
end

repo: Voice-Separation-and-Enhancement


text:
%EXsparsity  Example script, deblurring without and with sparsity

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% clear workspace and window.
clear, clc

% Choose if you would like to see the results displayed in a single figure
% window ('subplots') or in multiple figure windows ('manyplots').
% dispres = 'subplots';
dispres = 'manyplots';

rng(0) % make sure this test is repeatable

% Test problem: Gaussian deblurring of a sparse 'dotk' image (starry night).
n = 256;                        % Image is n-by-n.
NoiseLevel = 0.10;              % Relative noise level in data.
PRoptions.trueImage = 'dotk';   % Sparse test image.
[A, b, x, ProbInfo] = PRblur(n,PRoptions);
bn = PRnoise(b,'gauss',NoiseLevel);

% Compute the CGLS reconstruction.
options = IRset('MaxIter', 80, 'x_true', x, 'NoStop', 'on');
[Xcgls, info_cgls] = IRcgls(A, bn, options);

% Compute ell1 sparse reconstructions with default (GCV) inner stopping rule.
[Xell1_GCV, info_ell1_GCV] = IRell1(A, bn, options);

% Compute ell1 sparse reconstruction with DP as inner stopping rule.
options = IRset(options, 'RegParam', 'discrep', 'NoiseLevel', NoiseLevel,...
                'eta', 1.1);
[Xell1_DP, info_ell1_DP] = IRell1(A, bn, options);

options.NoStopOut = 'on';
K = 80;
[Xirn_DP, info_irn_DP] = IRirn(A, bn, K, options);

% Display the reconstructions;
% uncomment as appropriate to avoid displaying titles and legends.
if strcmp(dispres, 'subplots')
    figure(1), clf
    subplot(2,3,1)
    imagesc(reshape(x,n,n)), axis image off
    title('Exact image','fontsize',16,'interpreter','latex')
    subplot(2,3,2)
    imagesc(reshape(b,n,n)), axis image off
    title('Blurred image','fontsize',16,'interpreter','latex')
    subplot(2,3,3)
    imagesc(reshape(max(0,info_cgls.BestReg.X),n,n)), axis image off
    title(['Best IRcgls after ',num2str(info_cgls.BestReg.It),' iterations'],'fontsize',16,'interpreter','latex')
    subplot(2,3,4)
    imagesc(reshape(max(0,info_ell1_GCV.BestReg.X),n,n)), axis image off
    title(['Best IRell1 w/ GCV after ',num2str(info_ell1_GCV.BestReg.It),...
          ' iterations'],'fontsize',16,'interpreter','latex')
    subplot(2,3,5)
    imagesc(reshape(max(0,info_ell1_DP.BestReg.X),n,n)), axis image off
    title(['Best IRell1 w/ DP ',num2str(info_ell1_DP.BestReg.It),...
          ' iterations'],'fontsize',16,'interpreter','latex')
    subplot(2,3,6)
    imagesc(reshape(max(0,info_irn_DP.BestReg.X),n,n)), axis image off
    title(['Best IRirn w/ DP after', num2str(info_irn_DP.BestReg.It),...
         ' iterations'],'fontsize',16,'interpreter','latex')
elseif strcmp(dispres, 'manyplots')
    figure(1), clf
    XX = reshape(x,n,n);
    imagesc(XX), axis image off
    title('Exact image','fontsize',16,'interpreter','latex')
    axes('units','normalized','position',[0.17 0.1 1.15 0.25]);
    imagesc(XX(227:256,227:256)), axis image, caxis([0 max(XX(:))])
    set(gca,'Xtick',[],'Ytick',[],'LineWidth',2,'Xcolor','red','Ycolor','red')
    %
    figure(2), clf
    BB = reshape((max(b,0)),n,n);
    imagesc(BB), axis image off
    title('Blurred image','fontsize',16,'interpreter','latex')
    axes('units','normalized','position',[0.17 0.1 1.15 0.25]);
    imagesc(BB(227:256,227:256)), axis image, caxis([0 max(BB(:))])
    set(gca,'Xtick',[],'Ytick',[],'LineWidth',2,'Xcolor','red','Ycolor','red')
    %
    figure(3), clf
    XCG = reshape((max(0,info_cgls.BestReg.X)),n,n);
    imagesc(XCG), axis image off
    title(['Best IRcgls after ',num2str(info_cgls.BestReg.It),' iterations'],'fontsize',16,'interpreter','latex')
    axes('units','normalized','position',[0.17 0.1 1.15 0.25]);
    imagesc(XCG(227:256,227:256)), axis image, caxis([0 max(XCG(:))])
    set(gca,'Xtick',[],'Ytick',[],'LineWidth',2,'Xcolor','red','Ycolor','red')
    %
    figure(4), clf
    I = reshape((max(0,info_ell1_DP.BestReg.X)),n,n);
    imagesc(I), axis image off
    title(['Best IRell1 w/ GCV after ',num2str(info_ell1_GCV.BestReg.It),...
          ' iterations'],'fontsize',16,'interpreter','latex')
    axes('units','normalized','position',[0.17 0.1 1.15 0.25]);
    imagesc(I(227:256,227:256)), axis image, caxis([0 max(I(:))])
    set(gca,'Xtick',[],'Ytick',[],'LineWidth',2,'Xcolor','red','Ycolor','red')
    %
    figure(5), clf
    II = reshape((max(0,info_ell1_DP.BestReg.X)),n,n);
    imagesc(II), axis image off
    title(['Best IRell1 w/ DP after ',num2str(info_ell1_DP.BestReg.It),...
          ' iterations'],'fontsize',16,'interpreter','latex')
    axes('units','normalized','position',[0.17 0.1 1.15 0.25]);
    imagesc(II(227:256,227:256)), axis image, caxis([0 max(II(:))])
    set(gca,'Xtick',[],'Ytick',[],'LineWidth',2,'Xcolor','red','Ycolor','red')
    %
    figure(6), clf
    III = reshape((max(0,info_irn_DP.BestReg.X)),n,n);
    imagesc(III), axis image off
    title(['Best IRirn w/ DP after ',num2str(info_irn_DP.BestReg.It),...
          ' iterations'],'fontsize',16,'interpreter','latex')
    axes('units','normalized','position',[0.17 0.1 1.15 0.25]);
    imagesc(III(227:256,227:256)), axis image, caxis([0 max(III(:))])
    set(gca,'Xtick',[],'Ytick',[],'LineWidth',2,'Xcolor','red','Ycolor','red')
end

return

% A number of instructions useful to save the displayed figures follow;
% the defualt is not to execute them. If you wish to save the displayed
% figures in the dedicated 'Results' folder, please comment the above
% return statement
oldcd = cd;
if strcmp(dispres, 'subplots')
    try
        cd('Results')
    catch
        mkdir('Results')
        cd('Results')
    end
    figure(1), print -dpng -r300 EXsparsity
elseif strcmp(dispres, 'manyplots')
    try
        cd('Results')
    catch
        mkdir('Results')
        cd('Results')
    end
    figure(1), print -depsc -r300 EXsparsity_a
    figure(2), print -depsc -r300 EXsparsity_b
    figure(3), print -depsc -r300 EXsparsity_c
    figure(4), print -depsc -r300 EXsparsity_d
    figure(5), print -depsc -r300 EXsparsity_e
    figure(6), print -depsc -r300 EXsparsity_f
end
cd(oldcd)

% Uncomment the following return statement if you wish to save the
% displayed figures as MATLAB figures

% return

oldcd = cd;
if strcmp(dispres, 'subplots')
    try
        cd('Results')
    catch
        mkdir('Results')
        cd('Results')
    end
    figure(1), saveas('EXsparsity.fig')
elseif strcmp(dispres, 'manyplots')
    try
        cd('Results')
    catch
        mkdir('Results')
        cd('Results')
    end
    saveas(figure(1), 'EXsparsity_a.fig')
    saveas(figure(2), 'EXsparsity_b.fig')
    saveas(figure(3), 'EXsparsity_c.fig')
    saveas(figure(4), 'EXsparsity_d.fig')
    saveas(figure(5), 'EXsparsity_e.fig')
    saveas(figure(6), 'EXsparsity_f.fig')
end
cd(oldcd)

repo: IRtools


text:
clear all; close all; clc

n=100;
L=20; x=linspace(-L,L,n); y=x;
[X,Y]=meshgrid(x,y);

Xd=[];
for j=1:100
u=tanh(sqrt(X.^2+Y.^2)).*cos(angle(X+i*Y)-(sqrt(X.^2+Y.^2))+j/10);
f=exp(-0.01*(X.^2+Y.^2));
uf=u.*f;
Xd(:,j)=reshape(uf,n^2,1);
pcolor(x,y,uf), shading interp, colormap(hot), caxis([-1 1]), drawnow
end

%%
[U,S,V]=svd(Xd,0);

figure(2)
subplot(4,1,3)
plot(100*diag(S)/sum(diag(S)),'ko','Linewidth',[2])
subplot(4,1,4)
semilogy(100*diag(S)/sum(diag(S)),'ko','Linewidth',[2])
subplot(2,1,1)
plot(V(:,1:4),'Linewidth',[2])
legend('mode1','mode2','mode3','mode4')
set(gca,'Fontsize',[15],'Xtick',[0 20 40 60 80 100])
subplot(4,1,3), set(gca,'Fontsize',[15],'Ylim',[0 60],'Ytick',[0 20 40 60],'Xlim',[0 40],'Xtick',[0 10 20 30 40])
subplot(4,1,4), set(gca,'Fontsize',[15],'Ylim',[10^(-20) 10^2],'Ytick',[10^(-20) 10^(-10) 10^2],'Xlim',[0 40],'Xtick',[0 10 20 30 40])

figure(3)
for j=1:4
  subplot(4,4,j)
  mode=reshape(U(:,j),n,n);
  pcolor(X,Y,mode), shading interp,caxis([-0.03 0.03]), colormap(gray)
  axis off
end

%%

figure(11)
u=tanh(sqrt(X.^2+Y.^2)).*cos(angle(X+i*Y)-(sqrt(X.^2+Y.^2)));
f=exp(-0.01*(X.^2+Y.^2));
uf=u.*f;
subplot(3,3,1),pcolor(x,y,uf), shading interp, caxis([-1 1]), axis off
subplot(3,3,2),pcolor(x,y,abs(uf)), shading interp, caxis([-1 1]), axis off
subplot(3,3,3),pcolor(x,y,uf.^5), shading interp, caxis([-1 1]), axis off
colormap(gray)

%% TRANSLATION
figure(5)
n=200; L=20; x=linspace(-L,L,n); y=x;  % space
m=41; T=10; t=linspace(0,T,m);         % time
c=3;   % wave speed

X=[];
for j=1:m
    X(:,j)=exp(-(x+15-c*t(j)).^2).';  % data snapshots
end
[U,S,V]=svd(X);  % SVD decomposition

%%
figure(6)
subplot(2,2,1)
waterfall(x,t,X.'),colormap([0 0 0])
view(20,75)
set(gca,'Xlim',[-20 20],'Xtick',[-20 -10 0 10 20],'Ylim',[0 10], ...
    'Ytick',[0 5 10],'Zlim',[0 1],'Ztick',[0  1],'Fontsize',[12])

[U2,S2,V2]=svd(X);

subplot(4,2,2)
plot(100*diag(S2)/sum(diag(S2)),'ko','Linewidth',[2])
set(gca,'Xlim',[0 40],'Xtick',0:10:40,'Ylim',[0 8],'Ytick',[0 4 8])
subplot(4,2,4)
semilogy(100*diag(S2)/sum(diag(S2)),'ko','Linewidth',[2])
grid on
set(gca,'Xlim',[0 40],'Xtick',0:10:40,'Ylim',[10^(-3) 2*10^1],'Ytick',[10^(-3) 10^(-2) 10^(-1) 10^0 10^1])

figure(8)
subplot(2,1,1)
plot(x,U2(:,1:4),'Linewidth',[2]);
legend('mode1','mode2','mode3','mode4','Location','SouthEast')
set(gca,'Fontsize',[15],'Ylim',[-0.15 0.15],'Ytick',[-0.15 0 0.15])
subplot(2,1,2)
plot(t,V2(:,1:4),'Linewidth',[2])
set(gca,'Fontsize',[15],'Ylim',[-.3 0.3],'Ytick',[-0.3 0 0.3])

author: dynamicslab


text:
function Gamma = initGamma_random(T,K,D,Pstructure,Pistructure,episodic,priorOFFvsON)
% Return a Gamma timecourse corresponding to random transitions with a given
% dirichletdiag (i.e. the transition probability prior)
% For ehmm , D must be a 2x2 matrix of prior transitions - or left
%   empty
%   first state is ON, second is OFF  (always start in OFF)
%
% Author: Romesh Abeysuriya, University of Oxford (2017)
%         Diego Vidaurre, Aarhus University (2020)

if nargin < 4, Pstructure = []; end
if nargin < 5, Pistructure = []; end
if nargin < 6, episodic = false; end
if nargin < 7, priorOFFvsON = 3; end

rng('default')
rng('shuffle') % make this "truly" random

% Form transition probability matrix
if episodic
    P = (D-1)*eye(2)+ones(2);
    P(:,2) = P(:,2) * priorOFFvsON;
    P = bsxfun(@rdivide,P,sum(P,2));
    Pi = [0 1];
    if ~isempty(Pstructure) || ~isempty(Pistructure)
       warning('initGamma_random: Pstructure and Postructure will be ignored')
    end
else
    P = (D-1)*eye(K)+ones(K);
    if nargin>=4, P(~Pstructure) = 0; end
    if nargin<5, Pistructure = true(1,K); end
    P = bsxfun(@rdivide,P,sum(P,2));
    Pi = zeros(1,K);
    Pi(Pistructure) = 1;
    Pi = Pi / sum(Pi);
end

% Preallocate
Gamma = zeros(sum(T),K);

if episodic
    for n = 1:length(T)
        for k = 1:K
            while true
                gamma = zeros(T(n),2);
                gamma(1,2) = 1;
                for t = 2:T(n)
                    gamma(t,:) = mnrnd(1,P(gamma(t-1,:)==1,:));
                end
                if all(sum(gamma)>0), break; end
            end
            gamma = gamma + 0.0001 * rand(T(n),2);
            gamma = gamma ./ repmat(sum(gamma,2),1,2);
            t0 = sum(T(1:n-1)); t1 = sum(T(1:n));
            Gamma(t0+1:t1,k) = gamma(:,1);
        end
    end
elseif isChain(P) % If it is a chain, it P must be ordered
    for n = 1:length(T)
        ch = cumsum(rand(1,K));
        ch = ch / ch(end);
        ch = round(T(n) * ch);
        ch = [0 ch(1:end-1)];
        gamma = zeros(T(n),K);
        for k = 1:K-1, gamma(ch(k)+1:ch(k+1),k) = 1; end
        gamma(ch(k+1)+1:end,K) = 1;
        t0 = sum(T(1:n-1)); t1 = sum(T(1:n));
        Gamma(t0+1:t1,:) = gamma;
    end
else
    for n = 1:length(T)
        gamma = zeros(T(n),K);
        if sum(Pi>0)==1, gamma(1,Pi>0) = 1;
        else, gamma(1,:) = mnrnd(1,Pi);
        end
        for t = 2:T(n)
            gamma(t,:) = mnrnd(1,P(gamma(t-1,:)==1,:));
        end
        gamma = gamma + 0.0001 * rand(T(n),K);
        t0 = sum(T(1:n-1)); t1 = sum(T(1:n));
        Gamma(t0+1:t1,:) = gamma ./ repmat(sum(gamma,2),1,K);
    end
end

end

function bool = isChain(P)
%bool = any(sum(P,2)==1);
P = P~=0;
bool = all(sum(P(1:end-1,:),2)==2) && sum(P(end,:))==1;
end

author: OHBA-analysis

repo: HMM-MAR


text:
function overlap=compute_overlap(results,ground_truth)
gt_boxes = [ground_truth(:,1:2), ground_truth(:,1:2) + ground_truth(:,3:4) - ones(size(ground_truth,1), 2)];
 results.gt = gt_boxes;
    %   compute the OP
    pd_boxes = results.res;
    pd_boxes = [pd_boxes(:,1:2), pd_boxes(:,1:2) + pd_boxes(:,3:4) - ones(size(pd_boxes,1), 2)  ];
    lenALL=size(ground_truth,1);
     OP = zeros(size(gt_boxes,1),1);
    for i=1:size(gt_boxes,1)
        b_gt = gt_boxes(i,:);
        b_pd = pd_boxes(i,:);
        OP(i) = computePascalScore(b_gt,b_pd);
    end
    thresholdSetOverlap = 0:0.05:1;
    successNumOverlap = zeros(length(thresholdSetOverlap));
     for tIdx=1:length(thresholdSetOverlap)
         successNumOverlap(idx,tIdx) = sum(OP >thresholdSetOverlap(tIdx));
     end
     aveSuccessRatePlot = successNumOverlap/(lenALL+eps);
     overlap=mean(aveSuccessRatePlot);
end

author: Daikenan

repo: ASRCF

save_path: github-repos/MATLAB/Daikenan-ASRCF


text:
function img_f = GT_filter(img)
    % plug-in your favorite filter here
    % e.g. L0 smoothing filter
    img_f = L0Smoothing(img, 0.02);
end

author: jimmy-ren


text:
%%*************************************************************************
%% parameters.m: set OPTIONS structure to specify default
%%               parameters for sqlp.m
%%
%% OPTIONS.vers        : version of direction to use.
%%                       1 for HKM direction
%%                       2 for NT direction
%%                       0 for the default (which uses the HKM direction if
%%                       semidefinite blocks exist; and NT direction of SOCP problms)
%% OPTIONS.gam         : step-length parameter,
%% OPTIONS.predcorr    : whether to use Mehrotra predictor-corrector.
%% OPTIONS.expon       : exponent in decrease of centering parameter sigma.
%% OPTIONS.gaptol      : tolerance for duality gap as a fraction of the
%%                       value of the objective functions.
%% OPTIONS.inftol      : tolerance for stopping due to suspicion of
%%                       infeasibility.
%% OPTIONS.steptol     : toloerance for stopping due to small steps.
%% OPTIONS.maxit       : maximum number of iteration allowed
%% OPTIONS.printlevel  : 3, if want to display result in each iteration,
%%                       2, if want to display only summary,
%%                       1, if want to display warning message,
%%                       0, no display at all.
%% OPTIONS.stoplevel   : 2, if want to automatically detect termination;
%%                       1, if want to automatically detect termination, but
%%                          restart automatically with a new iterate
%%                          when the algorithm stagnants because of tiny step-lengths.
%%                       0, if want the algorithm to continue forever except for
%%                          successful completion, maximum number of iterations, or
%%                          numerical failures. Note, do not change this field unless
%%                          you very sure.
%% OPTIONS.scale_data  : 1, if want to scale the data before solving the problem,
%%                          else = 0
%% OPTIONS.rmdepconstr : 1, if want to remove nearly dependent constraints,
%%                          else = 0.
%% OPTIONS.smallblkdim : block-size threshold determining what method to compute the
%%                       schur complement matrix corresponding to semidefintie block.
%%                       NOTE: this number should be small, say less than 20.
%% OPTIONS.parbarrier  : parameter values of the log-barrier terms in the SQLP problem.
%%                       Default = [],  meaning that the parameter values are all 0.
%% OPTIONS.schurfun    : [], if no user supplied routine to compute the Schur matrix,
%%                       else, it is a cell array where each cell is either [],
%%                       or contains a string that is the file name where the Schur matrix
%%                       of the associated block data is computed.
%%                       For example, if the SQLP data has the block structure
%%                       blk{1,1} = '1'; blk{1,2} = 10;
%%                       blk{2,1} = 's'; blk{2,2} = 50;
%%                       and
%%                       OPTIONS.schurfun{1} = [];
%%                       OPTIONS.schurfun{2} = 'myownschur', where
%%                             'myownschur' is a function with the calling sequence:
%%                              function   schur = myownschur(X2,Z2inv,schurfun_par(2,:));
%%                       This means that for the first block, the Schur
%%                       matrix is computed by the default method in SDPT3,
%%                       and for the second block, the user supplies the
%%                       routine to compute the associated Schur matrix.
%% OPTIONS.schurfun_par: [], if no user supplied routine to compute the Schur matrix,
%%                       else, it is a cell array where the p-th row is either [],
%%                       or is a cell array containing the parameters needed in
%%                       the user supplied Schur routine OPTIONS.schurfun{p}.
%%                       For example, for the block structure described
%%                       above, we may have:
%%                       OPTIONS.schurfun_par{1} = [];
%%                       OPTIONS.schurfun_par{2,1} = par1;
%%                       OPTIONS.schurfun_par{2,2} = par2;
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*************************************************************************

function OPTIONS = sqlparameters

OPTIONS.vers           = 0;
OPTIONS.gam            = 0;
OPTIONS.predcorr       = 1;
OPTIONS.expon          = 1;
OPTIONS.gaptol         = 1e-8;
OPTIONS.inftol         = 1e-8;
OPTIONS.steptol        = 1e-6;
OPTIONS.maxit          = 100;
OPTIONS.printlevel     = 3;
OPTIONS.stoplevel      = 1; %% do not change this field unless you very sure.
OPTIONS.scale_data     = 0;
OPTIONS.spdensity      = 0.4;
OPTIONS.rmdepconstr    = 0;
OPTIONS.smallblkdim    = 50;
OPTIONS.parbarrier     = [];
OPTIONS.schurfun       = [];
OPTIONS.schurfun_par   = [];
%%*************************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function plotCustom(v,pcmd,varargin)
%
% Syntax
%   plotcustom(v,@(x,y) drawCommand(x,y))  %
%
% Input
%  v  - @vector3d
%  s  - string
%
% Options
%
% Output
%
% See also

% initialize spherical plot
sP = newSphericalPlot(v,varargin{:});

for j = 1:numel(sP)

  % project data
  [x,y] = project(sP(j).proj,v,varargin{:});

  % plot custom
  for i = 1:length(x), pcmd{1}(sP(j).hgt,x(i),y(i)); end
end

author: mtex-toolbox

repo: mtex

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