text:
%% Setup
p = get_default_truck_trailer_params();
p.noiseLevel = 0.01;
p.trailerWheelbase = 11;
forwardTarget = [20 2 0.3];
xx = [0 -5 0 0];
opt = odeset('Events', @obstacleEvents);
dt = 0.2;

%% setup viz
figure(101)
ax(1) = subplot(121);
tp1 = truck_trailer_plot(p, ax(1));
title('No learning')
ax(2) = subplot(122);
tp2 = truck_trailer_plot(p, ax(2));
title('Learning via ADP')
%% No learning
numPullups = 5;

tsave = [0];
ysave = [xx];

for ct = 1:numPullups

%     if ct <=3
%         p.forwardTarget = [20 3 0.1];
%     else
%         p.forwardTarget = [20 0 0];
%     end
%forward
if ysave(end,2) >= -0.3
   p.forwardTarget = [20 0 0];
else
   p.forwardTarget =forwardTarget;
end

p.velocity = 1;
[t,y] = ode45(@(t, x)sys_truck_trailer_wrapper(t, x, p), ...
    tsave(end)+[0 40], [ysave(end,1:4) zeros(1,10)], opt);

tsave = [tsave; t];
ysave = [ysave; y(:,1:4)];

%backward
p.velocity = -1;
[t,y] = ode45(@(t, x)sys_truck_trailer_wrapper(t, x, p), ...
    tsave(end)+[0 40], [ysave(end,1:4) zeros(1,10)], opt);

tsave = [tsave; t];
ysave = [ysave; y(:,1:4)];

end

ts = (tsave(1):dt:tsave(end))';
[~, ix] = unique(tsave);
ys = interp1(tsave(ix), ysave(ix,:), ts, 'linear', 'extrap');
ts1 = ts;
ys1 = ys;
%%
% for jj = 1:numel(ts)
%     tp1.updateFig(ys(jj,:))
%     drawnow
% end

%% With ADP

tsave = [0];
ysave = [xx];
psave = zeros(3,3,0);
ksave = zeros(3,1,0);

for ct = 1:numPullups

%     if ct <=3
%         p.forwardTarget = [20 3 0.1];
%     else
%         p.forwardTarget = [20 0 0];
%     end

if ysave(end,2) >= -0.3
   p.forwardTarget = [20 0 0];
else
   p.forwardTarget = forwardTarget;
end
%forward
p.velocity = 1;
[t,y] = ode45(@(t, x)sys_truck_trailer_wrapper(t, x, p), ...
    tsave(end)+[0 40], [ysave(end,1:4) zeros(1,10)], opt);
tsave = [tsave; t];
ysave = [ysave; y(:,1:4)];

%backward
p.velocity = -1;
xbase = kron(xx(2:4), xx(2:4));
kbase = [0 0 0];
sbase = [0 0 0 0 0 0];
qbase = 0;

for ii = 1:30
    [t, y] = ode45(@(t, x)sys_truck_trailer_wrapper(t, x, p), ...
        tsave(end)+[0 0.5], [y(end,1:4) zeros(1,10)]);

    tsave = [tsave; t];
    ysave = [ysave; y(:,1:4)];

    xbase = [xbase;
        kron(y(end, 2:4), y(end, 2:4))];
    qbase = [qbase;
        y(end, 5)];
    kbase = [kbase;
        y(end, 6:8)];
    sbase = [sbase;
        y(end, 9:14)];
end

[t, y] = ode45(@(t, x)sys_truck_trailer_wrapper(t, x, p), ...
    tsave(end)+[0 20], [y(end,1:4) zeros(1,10)], opt);
tsave = [tsave; t];
ysave = [ysave; y(:,1:4)];

% Learning
xbase = xbase(2:end,:) - xbase(1:end-1,:);
kbase = kbase(2:end,:);
qbase = qbase(2:end,:);
sbase = sbase(2:end,:);
l = [xbase(:,[1 2 3 5 6 9]) kbase sbase] \ -qbase;
Ve = [l(1) l(2)/2 l(3)/2;
    l(2)/2 l(4) l(5)/2;
    l(3)/2 l(5)/2 l(6)]

Ke = (p.R \ l(7:9)) * 0.5

D1 = 3.0;     % trailer wheelbase
D2 = 11.0;   % tractor wheelbase
A = [  0    -1     0;
    0     0     0;
    0     0    1/p.trailerWheelbase];
B = [0;
    -1/p.truckWheelbase;
    1/p.truckWheelbase];

P = lyap((A-B*p.feedbackGain)',p.Q + p.feedbackGain'*p.R*p.feedbackGain)
Knext = p.R \ B'*P

psave(:,:,end+1) = P;
ksave(:,:,end+1) = p.feedbackGain;
p.feedbackGain = Knext;

end

%
% Resample t
ts = (tsave(1):dt:tsave(end))';
[~, ix] = unique(tsave);
ys = interp1(tsave(ix), ysave(ix,:), ts, 'linear', 'extrap');
ts2 =ts;
ys2 =ys;
% %%
% for jj = 1:numel(ts)
%     tp2.updateFig(ys(jj,:))
%     drawnow
% end

[Ps, ~, Ks] = care(A,B,p.Q, p.R);
%% Parallel plot
%v = VideoWriter('newfile.mp4','MPEG-4');
%open(v)
for jj = 1:numel(ts1)
    tp1.updateFig(ys1(jj,:))
    if jj <= numel(ts2)
        tp2.updateFig(ys2(jj,:))
    end
    drawnow
 %   frame = getframe(gcf);
 %   writeVideo(v,frame);
end
%close(v)

%%
figure(301)
plot(1:5, norm(squeeze(ksave - Ks')))

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [f, df] = mySparseFun(X0, x, y, xu, s, sn)
if nargin == 0
    clc; close all;
    x = [-20:0.1:60]';
    y = rand(length(x), 1);
    xu = [-20:0.1:60]';
    n = length(x);
    s = 1.0;
    sn = 0.1;
    X0 = [0.1 0.1 log([2.0 2.0])]';
end
a = X0(1);
b = X0(2);
l = exp(X0(3));
l2 = l*l;
sf = exp(X0(4));
sf2 = sf*sf;
sn2 = sn*sn;

m = length(xu);
n = length(x);
[Kuu, dKuu] = CalKFun(xu, xu, X0(3:4));
Kuu = Kuu + sf2*1e-2*eye(m);    % stablize.
save('./KInfo.mat', 'Kuu', 'X0', 'xu');
if rank(Kuu) < m
    Kuu = Kuu + 1e-6*eye(m);
end
L0 = chol(Kuu, 'lower');
iL0 = inv(L0);
iKuu = iL0'*iL0;
logDetKuu = 2*sum(log(diag(L0)));

[Kfu, dKfu] = CalKFun(x, xu, X0(3:4));
Kuf = Kfu';
dKuf = {dKfu{1}', dKfu{2}'};

Qff = Kfu*iKuu*Kuf;
DiagKff = sf2*eye(n);
Gama = sn2*eye(n) + s*(DiagKff - diag(diag(Qff)));
iGama = diag( 1.0 ./ diag(Gama) );
logDetGama = sum( log(diag(Gama)) );
ApxK = Gama + Qff;   % use Woodbury to calculate its inverse and determinants: Z = Gama, U = V = Kfu, W = iKuu.
S = Kuu + Kfu'*iGama*Kfu;
% S = S + 1e-10*eye(m);  % for numerical stabitity.
% if rank(S) < m
%     S = S + 1e-6*eye(m);
% end
L1 = chol(S, 'lower');
logDetS = 2*sum( log(diag(L1)) );
iL1 = inv(L1);
iS = iL1'*iL1;
iApxK = iGama - iGama*Kfu*iS*Kfu'*iGama;
tmp = iApxK*ApxK;
tt = diag(tmp);
% max(tt) - min(tt)

logDetK = logDetGama - logDetKuu + logDetS;   % use Woodbury equation.
mx = a*x + b;
f = -0.5*(y-mx)'*iApxK*(y-mx) - 0.5*logDetK-0.5*n*log(2*pi);
% ff = -0.5*(y-mx)'*pinv(ApxK)*(y-mx) - 0.5*log(det(ApxK)) - 0.5*n*log(2*pi);
% -0.5*(y-mx)'*pinv(ApxK)*(y-mx)+0.5*(y-mx)'*iApxK*(y-mx)
%%%%%%%%%% calculate gradient.
dQ = {};
dGama = {};
dDiagK = {zeros(n, n), 2*DiagKff};
dK = {};
Alpha = iApxK*(y-mx);
df2m = [x'; ones(1, n)] * Alpha;
A = Alpha * Alpha' - iApxK;
df2c = zeros(2, 1);
for i = 1 : 1 : 2
     dQ{i} = dKfu{i}*iKuu*Kuf + Kfu*iKuu*dKuf{i}+Kfu*(-iKuu*dKuu{i}*iKuu)*Kuf;
     dGama{i} = s*dDiagK{i} - s*diag(diag(dQ{i}));
     dK{i} = dQ{i} + dGama{i};
     nSum = 0;
     for id = 1 : 1 : n
         nSum = nSum + A(id, :) * dK{i}(:, id);
     end
     df2c(i) = 0.5*nSum;
end
% TestWoodFun(Gama, iKuu, Kfu, Kfu);
%%%%%%%%% Test dQ.
% Eps = 1e-2;
% for i = 1 : 1 : 2
%     v = zeros(2, 1);
%     v(i) = Eps;
%     XNew = X0(3:4) + v;
%     [KfuNew, ~] = CalKFun(x, xu, XNew);
%     [KuuNew, ~] = CalKFun(xu, xu, XNew);
%     QffNew = KfuNew*inv(KuuNew)*KfuNew';
%     dQ_Est = (QffNew - Qff) / Eps;
%     dKuuNew = (KuuNew - Kuu)/Eps;
%     tmpDiff = dQ{i} - dQ_Est;
%     norm(tmpDiff(:))
% end

df = [df2m; df2c];

f = -f;
df = -df;
bTest = 1;
end

function [K, dK] = CalKFun(x, y, p)
l = exp(p(1));
sf = exp(p(2));
l2 = l*l;
sf2 = sf*sf;
m = length(x);
n = length(y);
dK = {};
X = repmat(x, 1, n);
Y = repmat(y', m, 1);
R = (X - Y) .*(X - Y);
K = sf2*exp(-0.5/l2*R);
dK{1} = 1/l2 * R.*K;
dK{2} = 2*K;
% tmp = zeros(m,n);
% for i = 1 : 1 : m
%     for j = 1 : 1 : n
%         r = x(i) - y(j);
%         r2 = r*r;
%         tmp(i, j) = sf2*exp(-0.5*r2/l2);
%     end
% end
% bTest = 1;
end

function [K] = TestWoodFun(Z, W, U, V)
K = Z + U * W * V';
InvZ = pinv(Z);
InvW = pinv(W);
S = InvW+V'*InvZ*U;
InvK = InvZ - InvZ*U*pinv(S)*V'*InvZ;
detK = det(Z)*det(W)*det(InvW + V'*InvZ*U);

A = K*InvK;
figure;
tt = diag(A);
plot(tt, 'b.' );
max(tt) - min(tt)
bTest = 1;
end

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
function [C,W,G,V2W] = qslim(V,F,t,varargin)
  % QSLIM Very simple wrapper for qslim
  %
  % [C,W,G,V2W] = qslim(V,F,t,'ParameterName',ParameterValue)
  %
  % Inputs:
  %   V  #V by 3 input mesh vertex positions
  %   F  #F by 3 input mesh triangle indices (1-indexed)
  %   t  target number of faces {0 for all collapse}
  %   Optional:
  %     'QSlimFlags' followed by qslim flags
  % Outputs:
  %   C  #collapse list of edge collapse structs with fields
  %     .a index into V of vertex collapsed to
  %     .b index into V of vertex collapsed from
  %     .da  1 by 3 displacement of vertex a
  %     .rm list of indices in F of faces removed by this collapse
  %     .mod list of indices in F of faces modified by this collapse
  %   W  #W by 3 collapsed mesh vertex positions
  %   G  #G by 3 collapsed mesh triangle indices (1-indexed)
  %   V2W  #V list of indices into W
  %
  % See also: readLOG, perform_edge_collapse, path_to_qslim
  %

  qslim_flags = '';
  if nargin <3
    t = 0;
  end

  v = 1;
  while v<=numel(varargin)
    switch varargin{v}
    case 'QSlimFlags'
      assert((v+1)<=numel(varargin));
      v = v+1;
      qslim_flags = varargin{v};
    otherwise
      error(['Unsupported parametername: ' varargin{v}]);
    end
    v = v+1;
  end

  prefix = tempprefix();
  input = [prefix '.smf'];
  output = [prefix '.log'];

  % Write input to file
  writeSMF(input,V,F);

    % prepare command string
    command = sprintf('%s -t %d -M log -q %s %s >%s', ...
      path_to_qslim,t,qslim_flags,input,output);
  try
    [status,result] = system(command);
    if status ~= 0
      error(result);
    end

    [VV,FF,C] = readLOG(output);
    % should match input
    assert(size(VV,1)==size(V,1));
    assert(size(FF,1)==size(F,1));
    assert(isequal(F,FF));

    delete(input);
    delete(output);
  catch ME
    fprintf('%s\n',command);
    rethrow(ME);
  end

  % This is expensive so only perform if required by output
  if nargout>1
    [W,G,V2W] = perform_edge_collapse(V,F,C,t);
  end

end

author: alecjacobson

repo: gptoolbox


text:
function [res] = lmc(x)

%LMC calculates the left medcouple, a robust measure of
%left tail weight
%
% The left medcouple is described in:
%    Brys, G., Hubert, M. and Struyf, A. (2006),
%    "Robust Measures of Tail Weight",
%    Computational Statistics and Data Analysis,
%    50 (No 3), 733-759.
%
% For the up-to-date reference, please consult the website:
%    wis.kuleuven.be/stat/robust.html
%
% Required input arguments:
%    x : Data matrix (rows=observations, columns=variables)
%
% I/O:
%    result=lmc(x);
%
% Example:
%    result = lmc([chi2rnd(5,1000,1) trnd(3,1000,1)]);
%
% The output of LMC is a vector containing the left medcouple
%     for each column of the data matrix x
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Guy Brys
% Last Update: 17/03/2006

if (nargin<1)
    error('No input arguments')
end
if (size(x,1)==1)
    x = x';
end
for (i=1:size(x,2))
    res(i) = -mc(x(x(:,i)<=prctile(x(:,i),50),i));
end

author: marianux

repo: ecg-kit

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


text:
function mono_total_next_grlex_test ( )

%*****************************************************************************80
%
%% MONO_TOTAL_NEXT_GRLEX_TEST tests MONO_TOTAL_NEXT_GRLEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MONO_TOTAL_NEXT_GRLEX_TEST\n' );
  fprintf ( 1, '  MONO_TOTAL_NEXT_GRLEX can list the monomials\n' );
  fprintf ( 1, '  in M variables, of total degree N,\n' );
  fprintf ( 1, '  one at a time, in graded lexicographic order.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We start the process with (0,0,...,0,N).\n' );
  fprintf ( 1, '  The process ends with (N,0,...,0,0)\n' );

  n = 3;
  m = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Let M = %d\n', m );
  fprintf ( 1, '      N = %d\n', n );
  fprintf ( 1, '\n' );

  x = [ 0, 0, n ];
  i = 1;

  while ( 1 )

    fprintf ( 1, '  %2d:', i );
    for j = 1 : m
      fprintf ( 1, '  %1d', x(j) );
    end
    fprintf ( 1, '\n' );

    if ( x(1) == n )
      break
    end

    x = mono_total_next_grlex ( m, n, x );
    i = i + 1;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [c, ceq] = ma_XenonPropMassConstraint(stateLog, eventID, lbProp, ubProp, bodyIDApply, celBodyData, maData)
%ma_liqFuelOxPropMassConstraint Summary of this function goes here
%   Detailed explanation goes here
    normFact = 1;

    if(ischar(eventID) && strcmpi(eventID,'final'))
        eventNum = max(stateLog(:,13));
    else
%         hMAMainGUI = findall(0,'tag','ma_MainGUI');
%         maData = getappdata(hMAMainGUI,'ma_data');
        [~, eventNum] = getEventByID(eventID, maData.script);
    end

    eventLog = stateLog(stateLog(:,13)==eventNum,:);

    bodyEventLog = eventLog(eventLog(:,8)==bodyIDApply,:);
    if(isempty(bodyEventLog))
        finalEntry = eventLog(end,:);
    else
        finalEntry = bodyEventLog(end,:);
    end

    propMass = finalEntry(12);
    bodyID = finalEntry(8);

    if(bodyID == bodyIDApply || bodyIDApply==-1)
        if(lbProp == ubProp)
            c = [0 0];
            ceq(1) = propMass - ubProp;
        else
            c(1) = lbProp - propMass;
            c(2) = propMass - ubProp;
            ceq = [0];
        end
        c = c/normFact;
        ceq = ceq/normFact;
    else
        c = [0 0];
        ceq = [0];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef SPX_CGLeastSquare < handle
% Conjugate gradient algorithm for least squares  or
% normal equations implementation
% We are solving the problem minimize \| A x - b \|_2^2.
% A is of size M x N.
% M > N.
% A^T A is of size N x N.
% b is of size M x 1.
% x is of size N x 1.
% A is full rank. Rank(A) = N.
% We assume that A is well conditioned.
% We can simultaneously solve S problems.

    properties
        % Settable and gettable properties

        % Maximum number of iterations for which the algorithm can run
        MaxIterations
        % Threshold of norm (in terms of percentage)
        NormThreshold
    end

    properties(SetAccess=private)
        % Gettable properties

        % The full rank tall matrix operator M x N
        A
        % The vectors to be solved M x S
        B
        % Measurements dimension
        M
        % Problem dimension
        N
        % Number of equations to be solved
        S
        % Number of iterations taken for solving the problem
        Iterations
        % The solution vectors
        X
        % Residual vectors in the CG algorithm
        CGResiduals
        % Residual norms at the end of the corresponding CG problem.
        CGResidualNorms
        % Residual vectors of the LS problem
        LSResiduals
        % Norms of the LS residuals at the end
        LSResidualNorms
        % Indicates if the problems converged
        % This would be false only if the algorithm crossed max iterations
        Converged
    end

    methods
        % Public methods
        function self = SPX_CGLeastSquare(A, B)
            % Constructor
            if isa(A, 'spx.dict.Operator')
                self.A = A;
            elseif ismatrix(A)
                self.A = spx.dict.MatrixOperator(A);
            else
                error('Unsupported operator.');
            end
            [self.M, self.N] = size(A);
            if self.M < self.N
                error('Only over-determined systems are supported')
            end
            % Ideally we should check the rank too.
            % But that would be too time consuming.
            [bm, self.S] = size(B);
            if bm ~= self.M
                error('Dimensions mismatch.');
            end
            self.B = B;
            self.MaxIterations = self.N ;
            self.NormThreshold = 1e-6;
        end

        function result = solve(self)
            aa = self.A; % M x N
            % We compute A' B in advance.
            atbb = self.A.apply_ctranspose(self.B); % N x S
            % Initial estimate vectors are all zeros
            xx = zeros(self.N, self.S); % N x S
            result = xx;
            self.Iterations = zeros(self.S, 1);
            self.CGResidualNorms = zeros(self.S, 1);
            self.Converged = false(self.S, 1);
            % Initial residual vectors
            rr = atbb  - aa.apply_ctranspose(aa.apply(xx)); % N x S
            self.CGResiduals  = rr;
            % Norm squared of initial residual vectors
            deltas = spx.norm.inner_product_cw(rr, rr);
            % The factor with which the norm needs to be reduced
            epsilon = self.NormThreshold;
            % Target limits on norm squared of residuals
            limits = epsilon^2 * deltas;
            % Maximum number of iterations for which the algorithm
            % is allowed to run
            imax = self.MaxIterations;
            % Number of problems being solved
            ns = self.S;
            for s=1:ns
                % Initialize iteration counter
                i = 0;
                % The quantities for this problem
                limit = limits(s);
                delta = deltas(s);
                % First residual
                r = rr(:, s); % N x 1.
                % First estimate
                x = xx(:, s); % N x 1.
                % Target b
                atb = atbb(:, s); % N x 1.
                % First direction
                d = r;
                while i < imax && delta > limit
                    % Compute the intermediate variable
                    p = aa.apply(d); % M x 1
                    q = aa.apply_ctranspose(p); % N x 1
                    % the line search scale factor in current direction
                    % Note that d' * q = p' * p.
                    alpha = delta / (p' * p);
                    % Update estimate in current direction
                    x = x + alpha * d;
                    if mod(i , 50) == 0
                        % In order to avoid propagation of floating point
                        % errors, we will recompute the value of residual
                        r = atb - aa.apply_ctranspose(aa.apply(x));
                    else
                        % Otherwise we use a shortcut
                        r = r  - alpha * q;
                    end
                    % hold the current residual norm squared
                    delta_old = delta;
                    % Update residual norm squared
                    delta = r' * r;
                    % Compute the ratio
                    beta = delta / delta_old;
                    % choose the new direction
                    d = r + beta * d;
                    % Increase iteration counter
                    i = i + 1;
                end
                % The problem has been solved
                result(:, s) = x;
                % Number of iterations taken to solve this problem
                self.Iterations(s) = i;
                self.CGResidualNorms(s) = sqrt(delta);
                self.CGResiduals(:, s) = r;
                self.Converged(s) = delta <= limit;
            end
            % Maintain the result for reference
            self.X = result;
            % Compute the corresponding LS residuals
            self.LSResiduals = self.B - aa.apply(self.X);
            self.LSResidualNorms = spx.norm.norms_l2_cw(self.LSResiduals);
        end

        function result = hasConverged(self)
            % Returns if all the solutions have converged
            result = all(self.Converged);
        end

        function printResults(self)
            ns = self.S;
            nn = self.N;
            mm = self.N;
            for s = 1:ns
                fprintf('Problem: %d\n', s);
                fprintf('Iterations: %d\n', self.Iterations(s));
                fprintf('CG Residual norm: %.2f, LS Residual norm: %.2f, Converged: %d\n', ...
                    self.CGResidualNorms(s), self.LSResidualNorms(s), self.Converged(s));
                if nn < 10
                    % We will print the solutions too
                    fprintf('Solution vector: ');
                    fprintf('%.4f ', self.X(:, s));
                    fprintf('\n');
                    fprintf('CG Residual vector: ');
                    fprintf('%.4f ', self.CGResiduals(:, s));
                    fprintf('\n');
                end
                if mm < 10
                    fprintf('LS Residual vector: ');
                    fprintf('%.4f ', self.LSResiduals(:, s));
                    fprintf('\n');
                end
            end
        end
    end

    methods(Access=private)
        % Private methods

    end

    methods(Static)
        % Public static methods

    end

end

author: indigits

repo: sparse-plex


text:
%GETNEXTCLIQUES Find a pair of cliques ready for message passing
%   [i, j] = GETNEXTCLIQUES(P, messages) finds ready cliques in a given
%   clique tree, P, and a matrix of current messages. Returns indices i and j
%   such that clique i is ready to transmit a message to clique j.
%
%   We are doing clique tree message passing, so
%   do not return (i,j) if clique i has already passed a message to clique j.
%
%	 messages is a n x n matrix of passed messages, where messages(i,j)
% 	 represents the message going from clique i to clique j.
%   This matrix is initialized in CliqueTreeCalibrate as such:
%      MESSAGES = repmat(struct('var', [], 'card', [], 'val', []), N, N);
%
%   If more than one message is ready to be transmitted, return
%   the pair (i,j) that is numerically smallest. If you use an outer
%   for loop over i and an inner for loop over j, breaking when you find a
%   ready pair of cliques, you will get the right answer.
%
%   If no such cliques exist, returns i = j = 0.
%
%   See also CLIQUETREECALIBRATE
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function [i, j] = GetNextCliques(P, messages)

% initialization
% you should set them to the correct values in your code
i = 0;
j = 0;

N = length(messages);
mess = ones(size(messages));

for k = 1:N
	for l = 1:N
		mess(k,l) = length(messages(k,l).var)!=0;
	end
end
%mess=messages;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k = 1:N
	for l = 1:N
		if P.edges(k,l)==0 || mess(k,l)==1
			continue;
		end
		kmess = mess(:,k)';
	 	edges = P.edges(:,k);
		edges(l) = 0;
		kmess(l) = 0;
		%if sum(kmess) == sum(P.edges(k,:))-1
		if sum(kmess' == edges)==N
			i = k;
			j = l;
			return;
		end
		if i!=0
			break;
		end
	end
	if i!=0
		break;
	end
end
return;

author: anhncs

repo: Probabilistic-Graphical-Models


text:
% %function [xml, rxml] = LoadXml_old(FileBase)
% %loads the xml file using xmltools_olds (have to have it in the path)
% % rxml returns it's original layout - very messy structure but contains all
% % the xml file contents.
% % xml - is the ouput structure which is backwards compatible to LoadPar
% % output, so you can use it instead ..also loads some usefull stuff -
% % Anatomoical groups with Skips , Spike electrode groups
% % more can be added later (e.g. parameters of the process scripts)
% % this script is written for xml version 1.1 .. older version doesn't work.
% % additions are welcome
% %
% % By default, fbasename is the name of the directory where the data are, so
% % the exact xml file path should not be put in argument
% % By calling
% % xml = LoadXml(xmlfile,'raw')
% % the xml file specified in argument will be loaded
% % (added by A Peyrache)
%

function [xml, rxml] = LoadXml_old(fbasename,varargin)
% xml = struct;
warning(['this function is now deprecated from /buzcode and a '...
        ' backwards compatible version of LoadParameters.m will run now..\n'])
xml = LoadParameters(fbasename);
rxml =  xml;

%if isempty(varargin)
%    [fbasename mergedir rootdir] = extractfbasename(fbasename);
%elseif strcmpi(varargin,'raw')
%    fprintf('Loading directly xml file')
%else
%    error('Unrecognized option')
%end
% %if xml was in the filebase by chance
% xmli = strfind(fbasename,'.xml');
% if ~isempty(xmli)
%     fbasename = fbasename(1:xmli-1);
% end
% rxml = xmltools_old([fbasename '.xml']);
%
% rxml = rxml.child(2);
%
% % from this level all children are the different parameters fields
%
% xml.FileName = fbasename;
%
% for i=1:length(rxml.child)
%
%     switch rxml.child(i).tag
%
%         case 'generalInfo'
%             xml.Date = rxml.child(i).child(1).value; % date of xml file creation?
%
%         case 'acquisitionSystem'
%             xml.nBits = str2num(rxml.child(i).child(1).value); % number of bits of the file
%             xml.nChannels = str2num(rxml.child(i).child(2).value);
%             xml.SampleRate = str2num(rxml.child(i).child(3).value);
%             xml.SampleTime = 1e6/xml.SampleRate; %to make backwards compatible
%             xml.VoltageRange = str2num(rxml.child(i).child(4).value);
%             xml.Amplification = str2num(rxml.child(i).child(5).value);
%             xml.Offset =  str2num(rxml.child(i).child(6).value);
%
%         case 'fieldPotentials'
%             xml.lfpSampleRate = str2num(rxml.child(i).child.value);
%
%         case 'anatomicalDescription'
%             tmp = rxml.child(i).child.child;
%             for grpI =1:length(tmp)
%                 for chI=1:length(tmp(grpI).child)
%                     xml.AnatGrps(grpI).Channels(chI) = str2num(tmp(grpI).child(chI).value);
%                     xml.AnatGrps(grpI).Skip(chI) = str2num(tmp(grpI).child(chI).attribs.value);
%                 end
%             end
%
%         case 'spikeDetection'
%             if ~isempty(rxml.child(i).child)
%                 tmp =rxml.child(i).child.child;
%                 for grpI =1:length(tmp)
%                     for chI=1:length(tmp(grpI).child(1).child)
%                         xml.SpkGrps(grpI).Channels(chI) = str2num(tmp(grpI).child(1).child(chI).value);
%                     end
%                     if length(tmp(grpI).child)>1
%                         xml.SpkGrps(grpI).nSamples = str2num(tmp(grpI).child(2).value);
%                         xml.SpkGrps(grpI).PeakSample = str2num(tmp(grpI).child(3).value);
%                         xml.SpkGrps(grpI).nFeatures = str2num(tmp(grpI).child(4).value);
%                     end
%                     %backwards compatibility
%                     xml.nElecGps = length(tmp);
%                     xml.ElecGp{grpI} = xml.SpkGrps(grpI).Channels;
%                 end
%             else
%                 xml.nElecGps = 0;
%             end
%
%
%         case 'programs'
%             tmp = rxml.child(i).child;
%             for i=1:length(tmp)
%                 if strcmp(tmp(i).child(1).value,'process_mhipass')
%                     for j=1:length(tmp(i).child(2).child )
%                         if strcmp(tmp(i).child(2).child(j).child(1).value,'frequency')
%                             xml.HiPassFreq = str2num(tmp(i).child(2).child(j).child(2).value);
%                             break
%                         end
%                     end
%                 end
%             end
%     end
%
%
% end
%
%
% % general recursive parsing will have to wait.

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef MivContextDef < handle
    % MivContextDef. Defines contexts for MIV
    %
    %
    %     There are a number of contexts, which each represent particular
    %     regions of the lung. For example, the OriginalImage context is the
    %     entire image, whereas the LungROI is the parrallipiped region
    %     containing the lungs and airways. The LeftLung and RightLung comprise
    %     just the voumes of the left and right lung respectively.
    %
    %     Context sets describe a collection of related contexts. For example,
    %     the SingleLung context set contains LeftLung and RightLung.
    %
    %     Each plugin specifies its context set, which is the domain of the
    %     results produced by the plugin. Some plugins operate over the whole
    %     LungROI region, while some operate on individual lungs.
    %
    %     This class manages the situations where a result is requested for a
    %     particualar context (e.g. LungROI) but the plugin defines a different
    %     context set (e.g. SingleLung). The LungROI can be built from the two
    %     contexts in the SingleLung set. Therefore in this case, the plugin is run twice,
    %     once for the left and once for the right lung. Then the resutls are
    %     combined to produce the result for the LungROI context.
    %
    %     Thie class defines this heirarchy of contexts and context sets. In
    %     this way plugins can operate on whichever context is appropriate,
    %     while results can be requested for any context, and the conversions
    %     are handled automatically.
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = private)
        Contexts
        ContextSets
    end

    methods
        function obj = MivContextDef
            obj.CreateContexts;
        end

        function context = GetContexts(obj)
            context = obj.Contexts;
        end

        function context_sets = GetContextSets(obj)
            context_sets = obj.ContextSets;
        end

        function context = GetDefaultContext(~)
            context = PTKContext.OriginalImage;
        end

        function context = GetDefaultContextSet(~)
            context = PTKContextSet.OriginalImage;
        end

        function matches = ContextSetMatches(~, plugin_context_set, requested_context_set)
            % Returns true if the plugin_context_set can be used to
            % generate the context_set without conversion

            matches = (plugin_context_set == requested_context_set) || (plugin_context_set == PTKContextSet.Any);
        end

        function output_context = ChooseOutputContext(~, context)
            % If a specific context has been specified in the plugin, we use
            % this (note this is not normally the case, as plugins usually
            % specify a PTKContextSet rather than a PTKContext)
            if isa(context, 'PTKContext')
                output_context = plugin_info.Context;

                % If the plugin specifies a PTKContextSet of type
                % PTKContextSet.OriginalImage, then we choose to return a context
                % of PTKContext.OriginalImage
            elseif context == PTKContextSet.OriginalImage
                output_context = PTKContext.OriginalImage;

            % In all other cases we choose a default context of the lung ROI
            else
                output_context = PTKContext.OriginalImage;
            end
        end

    end

    methods (Access = private)
        function CreateContexts(obj)
            % Create the hierarchy of context types
            obj.ContextSets = containers.Map;
            full_set =  MimContextSetMapping(PTKContextSet.OriginalImage, []);
            any_set = MimContextSetMapping(PTKContextSet.Any, []);
            obj.ContextSets(char(PTKContextSet.OriginalImage)) = full_set;
            obj.ContextSets(char(PTKContextSet.Any)) = any_set;

            % Create the hierarchy of contexts
            obj.Contexts = containers.Map;
            full_context =  MimContextMapping(PTKContext.OriginalImage, full_set, 'PTKGetContextForOriginalImage', 'PTKOriginalImage', []);

            obj.Contexts(char(PTKContext.OriginalImage)) = full_context;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [grRules,rxnGeneMat,indexes2check] = standardizeGrRules(model,embedded)
% standardizeGrRules
%   Standardizes gene-rxn rules in a model according to the following
%       - No overall containing brackets
%       - Just enzyme complexes are enclosed into brackets
%       - ' and ' & ' or ' strings are strictly set to lowercases
%
%   A rxnGeneMat matrix consistent with the standardized grRules is created.
%
%   model        a model structure
%   embedded     true if this function is called inside of another
%                RAVEN function (opt, default false)
%
%   grRules      [nRxns x 1] cell array with the standardized grRules
%   rxnGeneMat   [nRxns x nGenes]Sparse matrix consitent with the
%                standardized grRules
%
%   If this function is going to be used in a model reconstruction or
%   modification pipeline it is recommended to run this function just
%   at the beginning of the process.
%
%   Usage: [grRules,rxnGeneMat,indexes2check]=standardizeGrRules(model,embedded)

%Preallocate fields
n          = length(model.rxns);
[g,~]      = size(model.genes);
rxnGeneMat = sparse(n,g);
grRules    = cell(n,1);
genes      = model.genes;
if nargin<2
    embedded = false;
end

if isfield(model,'grRules')
    originalGrRules=model.grRules;
    originalGrRules=grRulesPreparation(originalGrRules);
    %Search for potential logical errors in the grRules field
    indexes2check = findPotentialErrors(originalGrRules,embedded,model);

    for i=1:length(originalGrRules)
        originalSTR = originalGrRules{i};
        grRules{i,:}  = originalSTR;
        %Non-empty grRules are splitted in all their different isoenzymes
        genesSets   = getSimpleGeneSets(originalSTR);
        rxnGeneMat  = modifyRxnGeneMat(genesSets,genes,rxnGeneMat,i);
        %standardize the non-conflicting grRules
        if ~ismember(i,indexes2check)
            newSTR = [];
            if ~isempty(genesSets)
                %For each simple genes set in the rule
                for j=1:length(genesSets)
                    simpleSet  = genesSets{j};
                    %Enclose simpleSet in brackets
                    if length(genesSets)>1
                        if ~isempty(strfind(simpleSet,' and '))
                            simpleSet = horzcat('(',simpleSet,')');
                        end
                    end
                    %Separate genesSets in the substring (in case of
                    %isoenzymes)
                    if j<length(genesSets)
                        newSTR = [newSTR, simpleSet, ' or '];
                        %Add the last simpleSet
                    else
                        newSTR = [newSTR, simpleSet];
                    end
                end
                %Update grRule
                grRules{i} = char(newSTR);
            end
        end
    end
end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Function that gets a cell array with all the simple geneSets in a given
%grRule string
function genesSets = getSimpleGeneSets(originalSTR)
genesSets  = [];
%If gene rule is not empty split in all its different isoenzymes
if ~isempty(originalSTR)
    originalSTR = strtrim(originalSTR);
    %Remove all brackets
    originalSTR = strrep(originalSTR,'(','');
    originalSTR = strrep(originalSTR,')','');
    %Split all the different genesSets
    genesSets  = transpose(strsplit(originalSTR,' or '));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Function that gets a cell array of simple genes sets (single genes or
%enzyme complexes) associated with the i-th reaction and modifies the
%correspondent row in the rxnGeneMat accordingly.
function rxnGeneMat = modifyRxnGeneMat(genesSets,modelGenes,rxnGeneMat,i)

if ~isempty(genesSets)
    for j=1:length(genesSets)
        simpleSet  = genesSets{j};
        %        rxnGeneMat = modifyRxnGeneMat(simpleSet,model.genes,rxnGeneMat,i);
        %Get individual genes
        STR   = strrep(simpleSet,') and (',' and ');
        genes = strsplit(STR,' ');
        for k=1:length(genes)
            if ~strcmpi(genes(k),' and ')
                %Get gene index
                genePos = find(strcmpi(modelGenes,genes(k)));
                if ~isempty(genePos)
                    rxnGeneMat(i,genePos) = 1;
                    %else
                    %In this case the gene should be added to the
                    %genes field (and to all of its dependencies)
                end
            end
        end
    end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Function that gets the model field grRules and returns the indexes of the
%rules in which the pattern ") and (" is present.
function indexes2check = findPotentialErrors(grRules,embedded,model)
indxs_l       = find(~cellfun(@isempty,strfind(grRules,') and (')));
indxs_l_L     = find(~cellfun(@isempty,strfind(grRules,') and')));
indxs_l_R     = find(~cellfun(@isempty,strfind(grRules,'and (')));
indexes2check = vertcat(indxs_l,indxs_l_L,indxs_l_R);
indexes2check = unique(indexes2check);

if ~isempty(indexes2check)

    if embedded
        EM = 'Potentially problematic ") AND (" in the grRules for reaction(s): ';
        dispEM(EM,false,model.rxns(indexes2check),true)
    else
        STR = 'Potentially problematic ") AND (", ") AND" or "AND ("relat';
        STR = [STR,'ionships found in\n\n'];
        for i=1:length(indexes2check)
            index = indexes2check(i);
            STR = [STR '  - grRule #' model.rxns{index} ': ' grRules{index} '\n'];
        end
        STR = [STR,'\n This kind of relationships should only be present '];
        STR = [STR,'in  reactions catalysed by complexes of isoenzymes e'];
        STR = [STR,'.g.\n\n  - (G1 or G2) and (G3 or G4)\n\n For these c'];
        STR = [STR,'ases modify the grRules manually, writing all the po'];
        STR = [STR,'ssible combinations e.g.\n\n  - (G1 and G3) or (G1 a'];
        STR = [STR,'nd G4) or (G2 and G3) or (G2 and G4)\n\n For other c'];
        STR = [STR,'ases modify the correspondent grRules avoiding:\n\n '];
        STR = [STR,' 1) Overall container brackets, e.g.\n        "(G1 a'];
        STR = [STR,'nd G2)" should be "G1 and G2"\n\n  2) Single unit en'];
        STR = [STR,'zymes enclosed into brackets, e.g.\n        "(G1)" s'];
        STR = [STR,'hould be "G1"\n\n  3) The use of uppercases for logi'];
        STR = [STR,'cal operators, e.g.\n        "G1 OR G2" should be "G'];
        STR = [STR,'1 or G2"\n\n  4) Unbalanced brackets, e.g.\n        '];
        STR = [STR,'"((G1 and G2) or G3" should be "(G1 and G2) or G3"\n'];
        warning(sprintf(STR))
    end
end
end

function grRules = grRulesPreparation(grRules)
%Remove unnecessary blanks
grRules=strrep(grRules,'  ',' ');
grRules=strrep(grRules,'( ','(');
grRules=strrep(grRules,' )',')');
% Make sure that AND and OR strings are in lowercase
grRules=strrep(grRules,' AND ',' and ');
grRules=strrep(grRules,' OR ',' or ');
end

author: SysBioChalmers

repo: RAVEN


text:
function angles = polyangles(x, y)
%POLYANGLES Computes internal polygon angles.
%   ANGLES = POLYANGLES(X, Y) computes the interior angles (in
%   degrees) of an arbitrary polygon whose vertices are given in
%   [X, Y], ordered in a clockwise manner.  The program eliminates
%   duplicate adjacent rows in [X Y], except that the first row may
%   equal the last, so that the polygon is closed.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.6 $  $Date: 2003/11/21 14:44:06 $

% Preliminaries.
[x y] = dupgone(x, y); % Eliminate duplicate vertices.
xy = [x(:) y(:)];
if isempty(xy)
   % No vertices!
   angles = zeros(0, 1);
   return;
end
if size(xy, 1) == 1 | ~isequal(xy(1, :), xy(end, :))
   % Close the polygon
   xy(end + 1, :) = xy(1, :);
end

% Precompute some quantities.
d = diff(xy, 1);
v1 = -d(1:end, :);
v2 = [d(2:end, :); d(1, :)];
v1_dot_v2 = sum(v1 .* v2, 2);
mag_v1 = sqrt(sum(v1.^2, 2));
mag_v2 = sqrt(sum(v2.^2, 2));

% Protect against nearly duplicate vertices; output angle will be 90
% degrees for such cases. The "real" further protects against
% possible small imaginary angle components in those cases.
mag_v1(~mag_v1) = eps;
mag_v2(~mag_v2) = eps;
angles = real(acos(v1_dot_v2 ./ mag_v1 ./ mag_v2) * 180 / pi);

% The first angle computed was for the second vertex, and the
% last was for the first vertex. Scroll one position down to
% make the last vertex be the first.
angles = circshift(angles, [1, 0]);

% Now determine if any vertices are concave and adjust the angles
% accordingly.
sgn = convex_angle_test(xy);

% Any element of sgn that's -1 indicates that the angle is
% concave. The corresponding angles have to be subtracted
% from 360.
I = find(sgn == -1);
angles(I) = 360 - angles(I);

%-------------------------------------------------------------------%
function sgn = convex_angle_test(xy)
%   The rows of array xy are ordered vertices of a polygon. If the
%   kth angle is convex (>0 and <= 180 degress) then sgn(k) =
%   1. Otherwise sgn(k) = -1. This function assumes that the first
%   vertex in the list is convex, and that no other vertex has a
%   smaller value of x-coordinate. These two conditions are true in
%   the first vertex generated by the MPP algorithm. Also the
%   vertices are assumed to be ordered in a clockwise sequence, and
%   there can be no duplicate vertices.
%
%   The test is based on the fact that every convex vertex is on the
%   positive side of the line passing through the two vertices
%   immediately following each vertex being considered.  If a vertex
%   is concave then it lies on the negative side of the line joining
%   the next two vertices. This property is true also if positive and
%   negative are interchanged in the preceding two sentences.

% It is assumed that the polygon is closed.  If not, close it.
if size(xy, 1) == 1 | ~isequal(xy(1, :), xy(end, :))
   xy(end + 1, :) = xy(1, :);
end

% Sign convention: sgn = 1 for convex vertices (i.e, interior angle > 0
% and <= 180 degrees), sgn = -1 for concave vertices.

% Extreme points to be used in the following loop.  A 1 is appended
% to perform the inner (dot) product with w, which is 1-by-3 (see
% below).
L = 10^25;
top_left = [-L, -L, 1];
top_right = [-L, L, 1];
bottom_left = [L, -L, 1];
bottom_right = [L, L, 1];

sgn = 1; % The first vertex is known to be convex.

% Start following the vertices.
for k = 2:length(xy) - 1
   pfirst= xy(k - 1, :);
   psecond = xy(k, :); % This is the point tested for convexity.
   pthird = xy(k + 1, :);
   % Get the coefficients of the line (polygon edge) passing
   % through pfirst and psecond.
   w = polyedge(pfirst, psecond);

   % Establish the positive side of the line w1x + w2y + w3 = 0.
   % The positive side of the line should be in the right side of the
   % vector (psecond - pfirst).  deltax and deltay of this vector
   % give the direction of travel. This establishes which of the
   % extreme points (see above) should be on the + side. If that
   % point is on the negative side of the line, then w is replaced by -w.

   deltax = psecond(:, 1) - pfirst(:, 1);
   deltay = psecond(:, 2) - pfirst(:, 2);
   if deltax == 0 & deltay == 0
      error('Data into convexity test is 0 or duplicated.')
   end
   if deltax <= 0  & deltay >= 0 % Bottom_right should be on + side.
      vector_product = dot(w, bottom_right); % Inner product.
      w = sign(vector_product)*w;
   elseif deltax <= 0 & deltay <= 0 % Top_right should be on + side.
      vector_product = dot(w, top_right);
      w = sign(vector_product)*w;
   elseif deltax >= 0 & deltay <= 0  % Top_left should be on + side.
      vector_product = dot(w, top_left);
      w = sign(vector_product)*w;
   else % deltax >= 0 & deltay >= 0, so bottom_left should be on + side.
      vector_product = dot(w, bottom_left);
      w = sign(vector_product)*w;
   end
   % For the vertex at psecond to be convex, pthird has to be on the
   % positive side of the line.
   sgn(k) = 1;
   if (w(1)*pthird(:, 1) + w(2)*pthird(:, 2) + w(3)) < 0
      sgn(k) = -1;
   end
end

%-------------------------------------------------------------------%
function w = polyedge(p1, p2)
%   Outputs the coefficients of the line passing through p1 and
%   p2. The line is of the form w1x + w2y + w3 = 0.

x1 = p1(:, 1);  y1 = p1(:, 2);
x2 = p2(:, 1);  y2 = p2(:, 2);
if x1==x2
   w2 = 0;
   w1 = -1/x1;
   w3 = 1;
elseif y1==y2
   w1 = 0;
   w2 = -1/y1;
   w3 = 1;
elseif x1 == y1 & x2 == y2
   w1 = 1;
   w2 = 1;
   w3 = 0;
else
   w1 = (y1 - y2)/(x1*(y2 - y1) - y1*(x2 - x1) + eps);
   w2 = -w1*(x2 - x1)/(y2 - y1);
   w3 = 1;
end
w = [w1, w2, w3];

%-------------------------------------------------------------------%
function [xg, yg] = dupgone(x, y)
% Eliminates duplicate, adjacent rows in [x y], except that the
% first and last rows can be equal so that the polygon is closed.

xg = x;
yg = y;
if size(xg, 1) > 2
   I = find((x(1:end-1, :) == x(2:end, :)) & ...
            (y(1:end-1, :) == y(2:end, :)));
   xg(I) = [];
   yg(I) = [];
end

author: HuangCongQing


text:
function [x,resL2,qualMeasOut]= LSMR(proj,geo,angles,niter,varargin)

% LSMR solves the CBCT problem using LSMR.
%
%  LSMR(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem
%   using the projection data PROJ taken over ALPHA angles, corresponding
%   to the geometry descrived in GEO, using NITER iterations.
%
%  LSMR(PROJ,GEO,ANGLES,NITER,OPT,VAL,...) uses options and values for solving. The
%   possible options in OPT are:
%
%  'lambda'  Value of parameter lambda, default 0.
%  'Init'    Describes diferent initialization techniques.
%             * 'none'     : Initializes the image to zeros (default)
%             * 'FDK'      : intializes image to FDK reconstrucition
%             * 'multigrid': Initializes image by solving the problem in
%                            small scale and increasing it when relative
%                            convergence is reached.
%             * 'image'    : Initialization using a user specified
%                            image. Not recomended unless you really
%                            know what you are doing.
%  'InitImg'    an image for the 'image' initialization. Avoid.
%  'groundTruth'  an image as grounf truth, to be used if quality measures
%                 are requested, to plot their change w.r.t. this known
%                 data.
%  'restart'  true or false. By default the algorithm will restart when
%             loss of ortogonality is found.
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Malena Sabate Landman, Ander Biguri
%--------------------------------------------------------------------------
%%

[verbose,x,QualMeasOpts,gpuids,lambda,gt,restart]=parse_inputs(proj,geo,angles,varargin);

measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:)));
if ~any(isnan(gt(:)))
    QualMeasOpts{end+1}='error_norm';
    x0=gt;
    clear gt
end
if nargout<3 && measurequality
    warning("Image metrics requested but none catched as output. Call the algorithm with 3 outputs to store them")
    measurequality=false;
end
qualMeasOut=zeros(length(QualMeasOpts),niter);
resL2 = zeros(1,niter);

% David Chin-Lung Fong and Michael Saunders //doi.org/10.1137/10079687X

% Enumeration as given in the paper for 'Algorithm LSMR'
iter=0;
remember=0;
while iter<niter
    % (1) Initialize
    u=proj-Ax(x,geo,angles,'Siddon','gpuids',gpuids);
    normr = norm(u(:),2);
    beta = normr;
    u = u/beta;

    v=Atb(u,geo,angles,'matched','gpuids',gpuids);
    alpha = norm(v(:),2); % msl: do we want to check if it is 0?
    v = v/alpha;

    alphabar = alpha;
    zetabar = alpha * beta;
    rho = 1;
    rhobar = 1;
    cbar = 1;
    sbar = 0;
    h = v;
    hbar = 0;

    % Compute the residual norm ||r_k||
    betadd = beta;
    betad = 0;
    rhod = 1;
    tautilda = 0;
    thetatilda = 0;
    zeta = 0;
    d = 0;

    % (2) Start iterations
    for ii=iter:niter
        iter=iter+1;
        if measurequality && ~strcmp(QualMeasOpts,'error_norm')
            x0 = x; % only store if necesary
        end
        if (iter==1 && verbose);tic;end

        % (3) Continue the bidiagonalization
        u = Ax(v,geo,angles,'Siddon','gpuids',gpuids) - alpha*u;
        beta = norm(u(:),2);
        u = u / beta;

        v = Atb(u,geo,angles,'matched','gpuids',gpuids) - beta*v;

        alpha = norm(v(:),2);
        v = v / alpha;

        % (4) Construct and apply rotation \hat{P}_k
        alphahat = sqrt(alphabar^2 + lambda^2);
        chat = alphabar/alphahat;
        shat = lambda/alphahat;

        % (5) Construct and apply rotation P_k
        rhopre = rho;
        rho = sqrt(alphahat^2 + beta^2);
        c = alphahat / rho;
        s = beta / rho;
        theta = s * alpha;
        alphabar = c * alpha;

        % (6) Construct and apply rotation \bar{P}_k
        thetabar = sbar * rho;
        rhobarpre = rhobar;
        rhobar = sqrt((cbar*rho)^2 + theta^2);
        cbar = cbar * rho / rhobar;
        sbar = theta / rhobar;
        zetapre = zeta;
        zeta = cbar * zetabar;
        zetabar = -sbar * zetabar;

        % (7) Update \bar{h}, x, h
        hbar = h - (thetabar*rho/(rhopre*rhobarpre))*hbar;
        x = x + (zeta / (rho*rhobar)) * hbar;
        h = v - (theta / rho) * h;

        % (8) Apply rotation \hat{P}_k, P_k
        betaacute = chat* betadd;
        betacheck = - shat* betadd;

        % Computing ||r_k||

        betahat = c * betaacute;
        betadd = -s * betaacute;

        % Update estimated quantities of interest.
        %  (9) If k >= 2, construct and apply \tilde{P} to compute ||r_k||
        rhotilda = sqrt(rhod^2 + thetabar^2);
        ctilda = rhod / rhotilda;
        stilda = thetabar / rhotilda;
        thetatildapre = thetatilda;
        thetatilda = stilda * rhobar;
        rhod = ctilda * rhobar;
        % betatilda = ctilda * betad + stilda * betahat; % msl: in the orinal paper, but not used
        betad = -stilda * betad + ctilda * betahat;

        % (10) Update \tilde{t}_k by forward substitution
        tautilda = (zetapre - thetatildapre*tautilda) / rhotilda;
        taud = (zeta - thetatilda*tautilda) / rhod;

        % (11) Compute ||r_k||
        d = d + betacheck^2;
        gamma_var = d + (betad - taud)^2 + betadd^2;
        aux = sqrt(gamma_var); % this is the residual teh algorithm follows, but we lose ortogonality, so we compute it explicitly

        % ||A^T r_k || is just |zetabar|

        % (6) Test for convergence.
        % msl: I still need to implement this.
        % msl: There are suggestions on the original paper. Let's talk about it!

        if measurequality
            qualMeasOut(:,iter)=Measure_Quality(x0,x,QualMeasOpts);
        end
        % The following should never happen, but the reallity is that if we use
        % the residual from the algorithm, it starts diverging from this explicit residual value.
        % This is an interesting fact that I believe may be caused either by
        % the mismatch of the backprojection w.r.t the real adjoint, or
        % numerical issues related to doing several order of magnitude
        % difference operations on single precission numbers.
        aux=proj-Ax(x,geo,angles,'Siddon','gpuids',gpuids);
        resL2(iter)=im3Dnorm(aux,'L2');
        if iter>1 && resL2(iter)>resL2(iter-1)
            % we lost orthogonality, lets restart the algorithm unless the
            % user asked us not to.

            % undo bad step.
            x=x-(zeta / (rho*rhobar)) * hbar;
            % if the restart didn't work.
            if remember==iter || ~restart
                disp(['Algorithm stoped in iteration ', num2str(iter),' due to loss of ortogonality.'])
                return;
            end
            remember=iter;
            iter=iter-1;
            if verbose
                disp(['Orthogonality lost, restarting at iteration ', num2str(iter) ])
            end
            break
        end

        if (iter==1 && verbose)
            expected_time=toc*niter;
            disp('LSMR');
            disp(['Expected duration   :    ',secs2hms(expected_time)]);
            disp(['Expected finish time:    ',datestr(datetime('now')+seconds(expected_time))]);
            disp('');
        end
    end
end
end

%% parse inputs'
function [verbose,x,QualMeasOpts,gpuids, lambda,gt,restart]=parse_inputs(proj,geo,angles,argin)
opts=     {'init','initimg','verbose','qualmeas','gpuids','lambda','groundtruth','restart'};
defaults=ones(length(opts),1);

% Check inputs
nVarargs = length(argin);
if mod(nVarargs,2)
    error('TIGRE:LSMR:InvalidInput','Invalid number of inputs')
end

% check if option has been passed as input
for ii=1:2:nVarargs
    ind=find(ismember(opts,lower(argin{ii})));
    if ~isempty(ind)
        defaults(ind)=0;
    else
        error('TIGRE:LSMR:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]);
    end
end

for ii=1:length(opts)
    opt=opts{ii};
    default=defaults(ii);
    % if one option isnot default, then extranc value from input
    if default==0
        ind=double.empty(0,1);jj=1;
        while isempty(ind)
            ind=find(isequal(opt,lower(argin{jj})));
            jj=jj+1;
        end
        if isempty(ind)
            error('TIGRE:LSMR:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]);
        end
        val=argin{jj};
    end

    switch opt
        case 'init'
            x=[];
            if default || strcmp(val,'none')
                x=zeros(geo.nVoxel','single');
                continue;
            end
            if strcmp(val,'FDK')
                x=FDK(proj,geo,angles);
                continue;
            end
            if strcmp(val,'multigrid')
                x=init_multigrid(proj,geo,angles);
                continue;
            end
            if strcmp(val,'image')
                initwithimage=1;
                continue;
            end
            if isempty(x)
                error('TIGRE:LSMR:InvalidInput','Invalid Init option')
            end
            % % % % % % % ERROR
        case 'initimg'
            if default
                continue;
            end
            if exist('initwithimage','var')
                if isequal(size(val),geo.nVoxel')
                    x=single(val);
                else
                    error('TIGRE:LSMR:InvalidInput','Invalid image for initialization');
                end
            end
            %  =========================================================================
        case 'qualmeas'
            if default
                QualMeasOpts={};
            else
                if iscellstr(val)
                    QualMeasOpts=val;
                else
                    error('TIGRE:LSMR:InvalidInput','Invalid quality measurement parameters');
                end
            end
        case 'verbose'
            if default
                verbose=1;
            else
                verbose=val;
            end
            if ~is2014bOrNewer
                warning('TIGRE:Verbose mode not available for older versions than MATLAB R2014b');
                verbose=false;
            end
        case 'gpuids'
            if default
                gpuids = GpuIds();
            else
                gpuids = val;
            end
        case 'lambda'
            if default
                lambda = 0;
            else
                lambda = val;
            end
        case 'groundtruth'
            if default
                gt=nan;
            else
                gt=val;
            end
        case 'restart'
            if default
                restart=true;
            else
                restart=val;
            end
        otherwise
            error('TIGRE:LSMR:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in LSMR()']);
    end
end

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
%% FEM2D_SCRIPT carries out the finite element assembly and solution procedures.
%
%  Discussion:
%
%    The BATCH command runs scripts, not functions.  So we have to write
%    this short script if we want to work with BATCH!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 March 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_SCRIPT\n' );
  fprintf ( 1, '  Call FEM2D_FUN to set up and solve the system.\n' );
%
%  Get the problem parameters.
%
  param = p_data ( );
%
%  Set up the mesh, the linear system, and solve.
%
  [ x, z_ss ] = fem2d_fun ( param );


text:
% DEMOPTIMISEGP Shows that there is an optimum for the covariance function length scale.
% DESC shows that by varying the length scale an artificial data
% set has different likelihoods, yet there is an optimum for which
% the likelihood is maximised.

% COPYRIGHT : Neil D. Lawrence, 2006, 2008

% GP

randn('seed', 1e5);
rand('seed', 1e5);

fillColor = [0.7 0.7 0.7];
markerSize = 20;
markerWidth = 2;
markerType = 'k.';
lineWidth = 2;

x = linspace(-1, 1, 6)';
trueKern = kernCreate(x, {'rbf', 'white'});
kern.comp{2}.variance = 0.001;
K = kernCompute(trueKern, x);
y = gsamp(zeros(1, 6), K, 1)';

xtest = linspace(-1.5, 1.5, 200)';
kern = trueKern;

lengthScale = [0.05 0.1 0.25 0.5 1 2 4 8 16];
counter = 0;

figure(1)
p = plot(x, y, markerType);
set(p, 'markersize', markerSize, 'lineWidth', markerWidth);
set(gca, 'fontname', 'times')
set(gca, 'fontsize', 18)
set(gca, 'ylim', [-2 1])
set(gca, 'xlim', [-1.5 1.5])

zeroAxes(gca);
fileName = ['demOptimiseGp' num2str(counter)];
if exist('printDiagram') && printDiagram
  printPlot(fileName, '../tex/diagrams', '../html');
end
clf

void = semilogx(NaN, NaN, 'k.-');
set(gca, 'fontname', 'times')
set(gca, 'fontsize', 18)
set(gca, 'ylim', [-12 -4])
set(gca, 'xlim', [0.025 32])
grid on
ylabel('log-likelihood')
xlabel('length scale')
fileName = ['demOptimiseGp' num2str(counter) '0'];
if exist('printDiagram') && printDiagram
  printPlot(fileName, '../tex/diagrams', '../html');
end

clf

for i = 1:length(lengthScale)
  kern.comp{1}.inverseWidth = 1/(lengthScale(i)*lengthScale(i));
  K = kernCompute(kern, x);
  [invK, U] = pdinv(K);
  logDetK = logdet(K, U);
  ll(i) = -0.5*(logDetK + y'*invK*y + size(y, 1)*log(2*pi));
  llLogDet(i) = -.5*(logDetK+size(y, 1)*log(2*pi));
  llFit(i) = -.5*y'*invK*y;
  Kx = kernCompute(kern, x, xtest);
  ypredMean = Kx'*invK*y;
  ypredVar = kernDiagCompute(kern, xtest) - sum((Kx'*invK).*Kx', 2);

  counter = counter + 1;
  figure(counter)
  clf
  fill([xtest; xtest(end:-1:1)], ...
       [ypredMean; ypredMean(end:-1:1)] ...
         + 2*[ypredVar; -ypredVar], ...
       fillColor,'EdgeColor',fillColor)
  hold on;
  t = plot(xtest, ypredMean, 'k-');

  p = plot(x, y, markerType);
  set(p, 'markersize', markerSize, 'lineWidth', markerWidth);
  set(t, 'linewidth', lineWidth);
  set(gca, 'fontname', 'times')
  set(gca, 'fontsize', 18)
  set(gca, 'ylim', [-2 1])

  zeroAxes(gca);
  fileName = ['demOptimiseGp' num2str(counter)];
  if exist('printDiagram') && printDiagram
    printPlot(fileName, '../tex/diagrams', '../html');
  end

  counter = counter + 1;
  figure(counter)
  t = semilogx(lengthScale(1:i), ll(1:i), 'k.-');
  hold on
  t = [t; semilogx(lengthScale(1:i), llLogDet(1:i), 'k.:')];
  t = [t; semilogx(lengthScale(1:i), llFit(1:i), 'k.--')];
  set(t, 'markersize', markerSize, 'lineWidth', markerWidth);
  set(gca, 'fontname', 'times')
  set(gca, 'fontsize', 18)
  set(gca, 'ylim', [-15 5])
  set(gca, 'xlim', [0.025 32])
  grid on
  ylabel('log-likelihood')
  xlabel('length scale')
  fileName = ['demOptimiseGp' num2str(counter)];
  if exist('printDiagram') && printDiagram
    printPlot(fileName, '../tex/diagrams', '../html');
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = WorstOfCall(S1,S2)
    nsim = size(S1,1);
    nt = size(S1,2);
    A = [(S1(:,nt)-S1(:,1))./S1(:,1) (S2(:,nt)-S2(:,1))./S2(:,1)];
    val = [min(A,[],2) zeros(nsim,1)];
    y = mean(max(val,[],2));
end

author: Sable

repo: mcbench-benchmarks


text:
function [rmag, dec] = geodet4 (lat, alt)

% geodetic to geocentric coordinates

% input

%  lat  = geodetic latitude (radians)
%         (+north, -south; -pi/2 <= lat <= +pi/2)
%  alt  = geodetic altitude (kilometers)

% output

%  rmag = geocentric position magnitude (kilometers)
%  dec  = geocentric declination (radians)
%         (+north, -south; -pi/2 <= dec <= +pi/2)

% global constants

%  req  = equatorial radius (kilometers)
%  flat = flattening factor (non-dimensional)

% Orbital Mechanics with Matlab

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

global req flat

% "normalize" altitude

hhat = alt / req;

hp1 = hhat + 1;

% calculate trig terms

s2lat = sin(2 * lat);

c2lat = cos(2 * lat);

s4lat = sin(4 * lat);

c4lat = cos(4 * lat);

% geocentric distance (kilometers)

rho = hp1 + (0.5 * (c2lat - 1)) * flat ...
      + ((1/(4 * hp1) + (1/16)) * (1 - c4lat)) * flat * flat;

rmag = req * rho;

% geocentric declination (radians)

dec = lat + (-s2lat/hp1) * flat  ...
      + (-s2lat/(2 * (hp1 * hp1)) + (1/(4 * hp1 * hp1) + 1/(4 * hp1))...
      * s4lat) * flat * flat;


text:
function alpha = sphericalAngle(p1, p2, p3)
%SPHERICALANGLE Compute angle between points on the sphere
%
%   ALPHA = sphericalAngle(P1, P2, P3)
%   Computes angle (P1, P2, P3), i.e. the angle, measured at point P2,
%   between the direction (P2, P1) and the direction (P2, P3).
%   The result is given in radians, between 0 and 2*PI.
%
%   Points are given either as [x y z] (there will be normalized to lie on
%   the unit sphere), or as [phi theta], with phi being the longitude in [0
%   2*PI] and theta being the elevation on horizontal [-pi/2 pi/2].
%
%
%   NOTE:
%   this is an 'oriented' version of the angle computation, that is, the
%   result of sphericalAngle(P1, P2, P3) equals
%   2*pi-sphericalAngle(P3,P2,P1). To have the more classical relation
%   (with results given betwen 0 and PI), it suffices to take the minimum
%   of angle and 2*pi-angle.
%
%   Examples
%     % Use inputs as cartesian coordinates
%     p1 = [0 1 0];
%     p2 = [1 0 0];
%     p3 = [0 0 1];
%     alpha = sphericalAngle(p1, p2, p3)
%     alpha =
%         1.5708
%
%     % Use inputs as spherical coordinates
%     sph1 = [.1 0];
%     sph2 = [0 0];
%     sph3 = [0 .1];
%     alphas = sphericalAngle(sph1, sph2, sph3)
%     alphas =
%         1.5708
%
%
%   See also:
%   geom3d, angles3d, spheres, sph2cart

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/02/2005.
%

%   HISTORY
%   23-05-2006 fix bug for points with angle from center > pi/2
%   05-06-2013 fix bug for points given as spherical coordinates, better
%       support for multiple inputs

% test if points are given as matlab spherical coordinates
if size(p1, 2) == 2
    [x, y, z] = sph2cart(p1(:,1), p1(:,2), ones(size(p1,1), 1));
    p1 = [x y z];
    [x, y, z] = sph2cart(p2(:,1), p2(:,2), ones(size(p2,1), 1));
    p2 = [x y z];
    [x, y, z] = sph2cart(p3(:,1), p3(:,2), ones(size(p3,1), 1));
    p3 = [x y z];
end

% normalize points
p1  = normalizeVector3d(p1);
p2  = normalizeVector3d(p2);
p3  = normalizeVector3d(p3);

% create the plane tangent to the unit sphere and containing central point
plane = createPlane(p2, p2);

% project the two other points on the plane
pp1 = planePosition(projPointOnPlane(p1, plane), plane);
pp3 = planePosition(projPointOnPlane(p3, plane), plane);

% compute angle on the tangent plane
pp2 = zeros(max(size(pp1, 1), size(pp3,1)), 2);
alpha = angle3Points(pp1, pp2, pp3);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [ fea, out ] = ex_magnetostatics5( varargin )
%EX_MAGNETOSTATICS5 Magnetic field test model.
%
%   [ FEA, OUT ] = EX_MAGNETOSTATICS5( VARARGIN ) Magnetic field test model.
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       sfun        string {sflag2}        Shape function for pressure
%       hmax        scalar {0.1}           Grid size
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { 'sfun',     'sflag2';
            'hmax',     0.1;
            'iplot',    1;
            'tol',      1e-2;
            'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Geometry and grid generation.
fea.sdim = { 'x' 'y' };
fea.grid = rectgrid( round(1/opt.hmax) );
fea.grid.s( selcells(fea,'y<=(0.5+sqrt(eps))') ) = 2;

% Problem definition.
fea = addphys( fea, @magnetostatics );
fea.phys.ms.eqn.coef{3,end} = { 1 0 };
fea.phys.ms.eqn.coef{4,end} = { 0 1 };
fea.phys.ms.sfun = { opt.sfun };

% Parse and solve problem.
fea       = parsephys( fea );
fea       = parseprob( fea );
fea.sol.u = solvestat( fea, 'icub', 2, 'fid', opt.fid );   % Call to stationary solver.

% Postprocessing.
if( opt.iplot>0 )
  figure
  postplot( fea, 'surfexpr', fea.phys.ms.eqn.vars{2,2}, ...
            'isoexpr', fea.phys.ms.eqn.vars{2,2}, 'isolev', 25, ...
            'arrowexpr', fea.phys.ms.eqn.vars{8,2}, 'arrowcolor', 'w', 'arrowspacing', [45 30] )
  title( 'Magnetic field' )
end

% Error checking.
Az1  = intsubd( fea.phys.ms.eqn.vars{1,2}, fea, 1 );
Az2  = intsubd( fea.phys.ms.eqn.vars{1,2}, fea, 2 );
Mf1  = intsubd( fea.phys.ms.eqn.vars{2,2}, fea, 1 );
Mf2  = intsubd( fea.phys.ms.eqn.vars{2,2}, fea, 2 );
Scb1 = intbdr(  fea.phys.ms.bdr.vars{2,2}, fea, 1 );
Scb2 = intbdr(  fea.phys.ms.bdr.vars{2,2}, fea, 2 );
Scb3 = intbdr(  fea.phys.ms.bdr.vars{2,2}, fea, 3 );
Scb4 = intbdr(  fea.phys.ms.bdr.vars{2,2}, fea, 4 );
out.err = [ abs(Az1+3.178289e-8)/3.178289e-8 ;
            abs(Az2+3.178289e-8)/3.178289e-8 ;
            abs(Mf1-0.389419)/0.389419 ;
            abs(Mf2-0.520161)/0.520161 ;
            abs(Scb1-0.16233)/0.16233 ;
            abs(Scb2-0.836552)/0.836552 ;
            abs(Scb3+0.83767)/0.83767 ;
            abs(Scb4+0.16343)/0.16343 ];
out.pass = all( out.err < opt.tol );

if ( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function w = dualtree3D(x, J, Faf, af)

% 3D Dual-Tree Discrete Wavelet Transform
%
% USAGE:
%   w = dualtree3D(x, J, Faf, af)
% INPUT:
%   x - 3-D array
%   J - number of stages
%   Faf - first stage filters
%   af - filters for remaining stages
% OUPUT:
%   w{j}{i}{d} - wavelet coefficients
%        j = 1..J, i = 1..4, d = 1..7
%   w{J+1}{i} - lowpass coefficients
%        i = 1..4
% EXAMPLE:
%   x = rand(64,64,64);
%   J = 3;
%   [Faf, Fsf] = FSfarras;
%   [af, sf] = dualfilt1;
%   w = dualtree3D(x, J, Faf, af);
%   y = idualtree3D(w, J, Fsf, sf);
%   err = x - y;
%   max(max(max(abs(err))))
%
% WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
% http://taco.poly.edu/WaveletSoftware/

% normalization
x = x/2;

M = [
    1 1 1
    2 2 1
    2 1 2
    1 2 2
];

for i = 1:4
    f1 = M(i,1);
    f2 = M(i,2);
    f3 = M(i,3);
    [xi w{1}{i}] = afb3D(x, Faf{f1}, Faf{f2}, Faf{f3});
    for k = 2:J
        [xi w{k}{i}] = afb3D(xi, af{f1}, af{f2}, af{f3});
    end
    w{J+1}{i} = xi;
end

for k = 1:J
    for m = 1:7
        [w{k}{1}{m} w{k}{2}{m} w{k}{3}{m} w{k}{4}{m}] = ...
            pm4(w{k}{1}{m}, w{k}{2}{m}, w{k}{3}{m}, w{k}{4}{m});
    end
end

author: thomaskuestner


text:
function [x, infos] = fro_mu_partial_nmf(V, rank, in_options)
%
% This file is part of NMFLibrary
%
% Created by H.Kasai on Feb. 16, 2017
%
% Change log:
%
%       June 16, 2022 (Hiroyuki Kasai): Initial version.
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.alg       = 'mu';
    local_options.norm_h    = 0;
    local_options.norm_w    = 1;
    local_options.alpha     = 2;
    local_options.delta     = 0.1;
    local_options.myeps     = 1e-16;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    if ~strcmp(options.alg, 'mu')
        fprintf('Invalid algorithm: %s. Therfore, we use mu (i.e., multiplicative update).\n', options.alg);
        options.alg = 'mu';
    end

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    W = init_factors.W;
    H = init_factors.H;

    % initialize
    method_name = sprintf('MU-Partial (%s:%s)', options.alg, options.metric_type);
    epoch = 0;
    R_zero = zeros(m, n);
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, W, H, R_zero, options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('MU-Partial (%s:%s): Epoch = 0000, cost = %.16e, optgap = %.4e\n', options.alg, options.metric, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        if strcmp(options.alg, 'mu')

            % update H
            H = H .* (W' * V) ./ (W' * W * H);
            H = H + (H<options.myeps) .* options.myeps;

            % update W
            W(:, options.updateW) = W(:, options.updateW) .* (V * H(options.updateW, :)') ...
                ./ (W * (H * H(options.updateW, :)'));
            W = W + (W<options.myeps) .* options.myeps;

        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, W, H, R_zero, options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = W;
    x.H = H;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
% RES = histoMatch(MTX, N, X)
%
% Modify elements of MTX so that normalized histogram matches that
% specified by vectors X and N, where N contains the histogram counts
% and X the histogram bin positions (see histo).

% Eero Simoncelli, 7/96.

function res = histoMatch(mtx, N, X)

if ( exist('histo') == 3 )
  [oN, oX] = histo(mtx(:), size(X(:),1));
else
  [oN, oX] = hist(mtx(:), size(X(:),1));
end

oStep = oX(2) - oX(1);
oC = [0, cumsum(oN)]/sum(oN);
oX = [oX(1)-oStep/2, oX+oStep/2];

N = N(:)';
X = X(:)';
N = N + 1e-10;   %% HACK: no empty bins ensures nC strictly monotonic

nStep = X(2) - X(1);
nC = [0, cumsum(N)]/sum(N);
nX = [X(1)-nStep/2, X+nStep/2];

nnX = interp1(nC, nX, oC, 'linear');

if ( exist('pointOp') == 3 )
  res = pointOp(mtx, nnX, oX(1), oStep);
else
  res = reshape(interp1(oX, nnX, mtx(:)),size(mtx,1),size(mtx,2));
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function GammaExposureTMOv(hdrv, filenameOutput, tmo_gamma, tmo_fstop, tmo_quality, tmo_video_profile)
%
%
%       GammaExposureTMOv(hdrv, filenameOutput, tmo_gamma, tmo_fstop, tmo_quality, tmo_video_profile)
%
%       This function applies gamma + exposure correction for each frame of the HDR stream
%
%       Input:
%           -hdrv: a HDR video structure; use hdrvread to create a hdrv
%           structure
%           -filenameOutput: output filename (if it has an image extension,
%           single files will be generated)
%           -tmo_gamma: gamma for encoding the frame. If it is negative,
%           sRGB econding is applied
%           -tmo_fstop: f-stop value
%           -tmo_quality: the output quality in [1,100]. 100 is the best quality
%           1 is the lowest quality.%
%           -tmo_video_profile: the compression profile (encoder) for compressing the stream.
%           Please have a look to the profile of VideoWriter from the MATLAB
%           help. Depending on the version of MATLAB some profiles may be not
%           be present.
%
%       Output:
%           -frameOut: the tone mapped frame
%
%     Copyright (C) 2016  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%
%     This function applies a static TMO to an operator without taking into
%     account for temporal coherency
%

if(~exist('filenameOutput', 'var'))
    date_str = strrep(datestr(now()), ' ', '_');
    date_str = strrep(date_str, ':', '_');
    filenameOutput = ['static_tmo_output_', date_str, '.avi'];
end

if(~exist('tmo_gamma', 'var'))
    tmo_gamma = 2.2;
end

if(~exist('tmo_fstop', 'var'))
    tmo_fstop = 0.0;
end

if(~exist('tmo_quality', 'var'))
    tmo_quality = 95;
end

if(~exist('tmo_video_profile', 'var'))
    tmo_video_profile = 'MPEG-4';
end

if(tmo_gamma < 0)
    bsRGB = 1;
else
    bsRGB = 0;
end

name = RemoveExt(filenameOutput);
ext = fileExtension(filenameOutput);

bVideo = 0;
writerObj = 0;

if(strcmp(ext, 'avi') == 1 | strcmp(ext, 'mp4') == 1)
    bVideo = 1;
    writerObj = VideoWriter(filenameOutput, tmo_video_profile);
    writerObj.FrameRate = hdrv.FrameRate;
    writerObj.Quality = tmo_quality;
    open(writerObj);
end

hdrv = hdrvopen(hdrv);

exposure = 2^tmo_fstop;

disp('Tone Mapping...');
for i=1:hdrv.totalFrames
    disp(['Processing frame ',num2str(i)]);
    [frame, hdrv] = hdrvGetFrame(hdrv, i);

    %only physical values
    frame = RemoveSpecials(frame);
    frame(frame < 0) = 0;

    %gamma/sRGB encoding
    if(bsRGB)
        frameOut = ClampImg(ConvertRGBtosRGB(frame * exposure, 0), 0, 1);
    else
        frameOut = ClampImg(GammaTMO(frame, tmo_gamma, tmo_fstop, 0), 0, 1);
    end

    if(bVideo)
        writeVideo(writerObj, frameOut);
    else
        nameOut = [name, sprintf('%.10d',i), '.', ext];
        imwrite(frameOut, nameOut);
    end

end
disp('OK');

if(bVideo)
    close(writerObj);
end

hdrvclose(hdrv);

end

author: banterle


text:
%--------------------------------------------------------------------------
% C = lrsc(A,tau)
% Low Rank Subspace Clustering algorithm for clean data lying in a
% union of subspaces
%
% C = argmin |C|_* + tau/2 * |A - AC|_F^2 s.t. C = C'
%
% A: clean data matrix whose columns are points in a union of subspaces
% tau: scalar parameter
%--------------------------------------------------------------------------
% Adapted from LRSC by @ Rene Vidal, November 2012
%--------------------------------------------------------------------------

function C = clean_relaxed(A,tau)
% Make an estimate of tau if necessary
    if nargin < 2
        tau = 100/norm(A)^2;
    end
    threshold = 1/sqrt(tau)
    options = struct;
    options.lambda = threshold;
    options.tolerance =  16*eps;
    M = size(A, 1);
    options.p0 = ones(M, 1);
    [~, S, V, details] = spx.fast.lansvd(A, options);
    r = numel(S);
    C = V * (eye(r) - diag(1./(S.^2)/tau)) * V';
end


text:
function [k,sigma,mu] = myDistFit(Temperature)
%MYDISTFIT    Create plot of datasets and fits
%   MYDISTFIT(TEMPERATURE)
%   Creates a plot, similar to the plot in the main distribution fitting
%   window, using the data that you provide as input.  You can
%   apply this function to the same data you used with dfittool
%   or with different data.  You may want to edit the function to
%   customize the code and this help message.
%
%   Number of datasets:  1
%   Number of fits:  1

% This function was automatically generated on 04-Dec-2006 08:27:07

% Data from dataset "Temperature data":
%    Y = Temperature

% Force all inputs to be column vectors
Temperature = Temperature(:);

% Set up figure to receive datasets and fits
f_ = clf;
figure(f_);
set(f_,'Units','Pixels','Position',[654 334.5 680 468.45]);
legh_ = []; legt_ = {};   % handles and text for legend
ax_ = newplot;
set(ax_,'Box','on');
hold on;

% --- Plot data originally in dataset "Temperature data"
t_ = ~isnan(Temperature);
Data_ = Temperature(t_);
[F_,X_] = ecdf(Data_,'Function','cdf'...
              );  % compute empirical cdf
Bin_.rule = 1;
[C_,E_] = dfswitchyard('dfhistbins',Data_,[],[],Bin_,F_,X_);
[N_,C_] = ecdfhist(F_,X_,'edges',E_); % empirical pdf from cdf
h_ = bar(C_,N_,'hist');
set(h_,'FaceColor','none','EdgeColor',[0.333333 0 0.666667],...
       'LineStyle','-', 'LineWidth',1);
xlabel('Data');
ylabel('Density')
legh_(end+1) = h_;
legt_{end+1} = 'Temperature data';

% Nudge axis limits beyond data limits
xlim_ = get(ax_,'XLim');
if all(isfinite(xlim_))
   xlim_ = xlim_ + [-1 1] * 0.01 * diff(xlim_);
   set(ax_,'XLim',xlim_)
end

x_ = linspace(xlim_(1),xlim_(2),100);

% --- Create fit "Generalized Extreme Value Fit"

% Fit this distribution to get parameter values
t_ = ~isnan(Temperature);
Data_ = Temperature(t_);
% To use parameter estimates from the original fit:
%     p_ = [ -0.1140337456229, 30.99486601834, 23.95766494166];
p_ = gevfit(Data_, 0.05);
y_ = gevpdf(x_,p_(1), p_(2), p_(3));
h_ = plot(x_,y_,'Color',[1 0 0],...
          'LineStyle','-', 'LineWidth',2,...
          'Marker','none', 'MarkerSize',6);
legh_(end+1) = h_;
legt_{end+1} = 'Generalized Extreme Value Fit';

hold off;
leginfo_ = {'Orientation', 'vertical', 'Location', 'NorthEast'};
h_ = legend(ax_,legh_,legt_,leginfo_{:});  % create legend
set(h_,'Interpreter','none');

% ---- Added code after generation
k = p_(1);
sigma = p_(2);
mu = p_(3);


text:
function volume = torus_volume_3d ( r1, r2 )

%*****************************************************************************80
%
%% TORUS_VOLUME_3D computes the volume of a torus in 3D.
%
%  Discussion:
%
%    A torus with radii R1 and R2 is the set of points (X,Y,Z)
%    satisfying:
%
%    ( sqrt ( X**2 + Y**2 ) - R1 )**2 + Z**2 <= R2**2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R1, R2, the "inner" and "outer" radii of the
%    torus.
%
%    Output, real VOLUME, the volume of the torus.
%
  volume = 2.0 * pi * pi * r1 * r2 * r2;

  return
end

author: johannesgerer


text:
function ap = context_rescore(train_set, train_year)
% Train context rescoring SVMs and rescore the test predictions.
%   ap = context_rescore(train_set, train_year)
%
% Return value
%   ap            AP scores for all classes after context rescoring
%
% Arguments
%   train_set     Training dataset
%   train_year    Training dataset year

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
% Copyright (C) 2008, 2009, 2010 Pedro Felzenszwalb, Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

if nargin < 2
  conf = voc_config();
  train_year = conf.pascal.year;
  if nargin < 1
    train_set = conf.training.train_set_fg;
  end
end

context_train(train_set, train_year);
ap = context_test();

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Doyen Sahoo
% Contributors: Bin LI, Steven C.H. Hoi
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [] = resultManager(results, job)
%% The results analysis of the Algorithm Analyser can be seen here

    % r.returns_daily     = returns;
    % r.returns           = cumprod(returns_daily+1);
    % r.portfolio         = portfolio;
    % r.stats             = stats;
    % r.benchmarks        = benchmakrs;
    % r.benchmarks_daily  = benchmarks_daily;
    % r.chosenStrategy    = chosenStrategy;
    % r.dataFrequency     = dataFrequency;

    [results.r, results.c]  = size(results.portfolio);

    menuId = 6; % Result Manager for Algorithm Analyser
    displayMenu(menuId);
    prompt = 'Please enter your choice (1-6):';
    choice = input(prompt);

    switch(choice),
        case 1,
            % Print table of results
            displayTable(results);
            resultManager(results, job);

        case 2,
            % Plot the graphs of returns (both cumulative return and
            % log(cumulative return) of all algorithms
            cumReturns      = [results.benchmarks results.returns];
            dailyReturns    = [(results.benchmarks_daily-1) results.returns];

            figure;
            plot(cumReturns);
            xlim([0 results.r+1]);
            legend('Uniform BAH', 'Uniform CRP', 'Best Stock', 'Best CRP', results.chosenStrategy, 'Location', 'Best');
            grid on;
            title('Cumulative Returns of the Algorithms');
            xlabel('Time');
            ylabel('Portfolio Value');

            figure;
            semilogy(cumReturns);
            xlim([0 results.r+1]);
            legend('Uniform BAH', 'Uniform CRP', 'Best Stock', 'Best CRP', results.chosenStrategy, 'Location', 'Best');
            grid on;
            title('LOG Cumulative Returns of the Algorithms');
            xlabel('Time');
            ylabel('Portfolio Value');

            resultManager(results, job);

        case 3,
            % Produce all the risk plots. Tell user to chaneg the window
            % size in the config if needed
            load ../GUI/config/config.mat
            disp('Generating all the Risk Plots');
            msg = 'windowSize = ';
            msg = strcat(msg, num2str(windowRisk));
            disp(msg);
            disp('(You can modify the window size by changing windowRisk variable in the config)');

            % Produce all risk plots based on windowRisk
            window = windowRisk;
            if isnan(window) || window < 2
                errorMessage('Window size must be a number greater than or equal to 2');
            else
                window = round(window);

                % 1. Annualized Sharpe Ratio
                addpath('../GUI/lib');

                [r, c] = size(results.returns_daily);
                sharpeRatios = zeros(r,1);
                for i = 50:1:r
                    start = max(i-window+1, 1);
                    finish = i;
                    sharpeRatios(i) = sharpe(results.returns_daily(start:finish), results.returns(finish)/results.returns(start), results.dataFrequency);
                end
                figure;
                plot(sharpeRatios);
                xlim([0 r+1]);
                legend('Sharpe Ratio', 'Location', 'Best');
                title('Sharpe Ratio of the Algorithm');
                xlabel('Time');
                ylabel('Sharpe Ratio');

                % 2. Calmar Ratio
                calmarRatios = zeros(r,1);
                for i = 50:1:r
                    start = max(i-window+1, 1);
                    finish = i;
                    mdd = maxDD(results.returns(start:finish));
                    calmarRatios(i) = calmar(results.returns(start:finish), mdd, results.dataFrequency);
                end
                figure;
                plot(calmarRatios);
                xlim([0 r+1]);
                legend('Calmar Ratio', 'Location', 'Best');
                title('Calmar Ratio of the Algorithm');
                xlabel('Time');
                ylabel('Calmar Ratio');

                % 3. Sortino Ratio
                sortinos = zeros(r,1);
                for i = 50:1:r
                    start = max(i-window+1, 1);
                    finish = i;
                    sortinos(i) = sortino(results.returns_daily(start:finish),0);
                end

                figure;
                plot(sortinos);
                xlim([0 r+1]);
                legend('Sortino Ratio', 'Location', 'Best');
                title('Sortino Ratio of the Algorithm');
                xlabel('Time');
                ylabel('Sortino Ratio');

                % 4. Value at Risk
                vars = zeros(r,1);
                for i = 50:1:r
                    start = max(i-window+1, 1);
                    finish = i;
                    vars(i) = var5(results.returns_daily(start:finish));
                end

                figure;
                plot(vars);
                xlim([0 r+1]);
                legend('Value at Risk', 'Location', 'Best');
                title('Value at Risk of the Algorithm');
                xlabel('Time');
                ylabel('Value at Risk');

                % 5. Maximum Draw Down
                mdds = zeros(r,1);
                for i = 50:1:r
                    start = max(i-window+1, 1);
                    finish = i;
                    mdds(i) = maxDD(results.returns(start:finish));
                end

                figure;
                plot(mdds);
                xlim([0 r+1]);
                legend('Maximum Draw Down', 'Location', 'Best');
                title('Maximum Draw Down of the Algorithm');
                xlabel('Time');
                ylabel('Maximum Draw Down (%)');

                rmpath('../GUI/lib');
            end
            resultManager(results, job);

        case 4,
            % Portfolio allocation plots
            portfolio   = results.portfolio;
            expected    = mean(portfolio);
            deviation   = std(portfolio);

            errorbar(expected, deviation, 'xr');
            xlim([0 results.c+1]);

            title('Average Portfolio Allocation');
            xlabel('Assets');
            ylabel('Fraction of Portfolio');

            resultManager(results, job);

        case 5,
            % Save the results.
            prompt ='Enter Name of File to save results: ';
            filename = input(prompt);
            [ resultsTable ] = getTable(results);
            cd ../Log/Results/
            save(filename, 'results', 'job', 'resultsTable');
            disp('Results saved in /Log/Result ');
            cd ../../PGUI/
            resultManager(results, job);

        case 6,
            disp('Exiting Result Manager --> to Algorithm Analyser');
            [ job ] = jobInit();
            algorithmAnalyserMenu( job );

        otherwise,
            disp('ERROR: Please enter a valid input');
            resultManager(results, job);
    end

end

function [ algorithmJob ] = jobInit()
% Construct an instance of an algorithm analyser job

    % Read configuration
    load ../GUI/config/config.mat;
    algorithmJob.algorithmId    = 1;
    algorithmJob.datasetId      = 1;
    AL                          = char(algorithmList);
    algorithmJob.algorithm      = AL(algorithmJob.algorithmId,:);
    DL                          = char(dataList);
    algorithmJob.dataset        = DL(algorithmJob.datasetId,:);
    algorithmJob.parameters     = cell2mat(defaultParameters(algorithmJob.algorithmId,:));
end

function [ resultsTable ] = getTable(results)

    cumReturns      = [results.benchmarks results.returns];
    dailyReturns    = [results.benchmarks_daily-1 results.returns_daily];

    % Compute the statsand display the important numbers in table using the
    % library functions
    addpath('../GUI/lib');

    % Get final Values
    [r c] = size(cumReturns);
    results.finalValues = cumReturns(r,:);

    % Get the mean returnfor every day - This is a simple average
    results.meanReturns = mean(dailyReturns);

    % Get annualised returns
    denominator = 252/results.dataFrequency;
    Y = r/denominator;
    results.annualisedReturns = results.finalValues.^(1/Y)-1;

    % Get standard deviation- a measureof risk
    results.standardDeviation = std(dailyReturns);

    % Get annualised standard deviation
    results.annualisedStandardDeviation = results.standardDeviation * sqrt(denominator);

    % Get sharpe ratios
    results.sharpeRatios = sharpe(dailyReturns, results.finalValues,results.dataFrequency);

    % Get Sortino ratios
    results.sortinoRatios = sortino(dailyReturns, 0);

    % Get Value risks at level 5%
    results.valueAtRisks = var5(dailyReturns);

    % Get Maximum draw down
    results.mdds = maxDD_general(cumReturns);

    % Get Calmar ratios
    results.calmars = calmar(cumReturns, results.mdds, results.dataFrequency);

    % Fill up the tables

    tableData   = [results.finalValues; results.meanReturns; results.annualisedReturns; results.standardDeviation; results.annualisedStandardDeviation; results.sharpeRatios; results.calmars; results.sortinoRatios; results.valueAtRisks; results.mdds];

    Market      = tableData(:,1);
    Uniform     = tableData(:,2);
    BestStock   = tableData(:,3);
    BCRP        = tableData(:,4);
    Algorithm   = tableData(:,5);
    report      = {'Final Value','Mean Return','Annualised Return','Standard Deviation','Annualised Standard Deviation','Sharpe Ratio','Calmar Ratio','Sortino Ratio','Value at Risk','Maximum Draw Down'};
    % tableData   = table(Market,Uniform,BestStock,BCRP, Algorithm,
    % 'RowNames', report); % - works only in Matlab 2014a onwards

    finalDisplay = report';
    finalDisplay(2:end+1) = finalDisplay;
    finalDisplay{1} = '';

    finalDisplay{1,2} = 'Market';
    finalDisplay{1,3} = 'Uniform';
    finalDisplay{1,4} = 'BestStock';
    finalDisplay{1,5} = 'BCRP';
    finalDisplay{1,6} = 'Algorithm';

    for i = 1:1:10
        for j = 1:1:5
            finalDisplay{i+1,j+1} = tableData(i,j);
        end
    end

    resultsTable = finalDisplay;
end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function FR = FbApply2d( I, FB, shape, show )
% Applies each of the filters in the filterbank FB to the image I.
%
% To apply to a stack of images:
%  IFS = fevalArrays( images, @FbApply2d, FB, 'valid' );
%
% USAGE
%  FR = FbApply2d( I, FB, [shape], [show] )
%
% INPUTS
%  I       - 2D input array
%  FB      - filterbank - MxNxK set of K filters each of size MxN
%  shape   - ['full'] option for conv2 'full', 'same', 'valid'
%  show    - [0] first figure to use for optional display
%
% OUTPUTS
%  FR      - 3D set of filtered images
%
% EXAMPLE
%  load trees;  X=imresize(X,.5);  load FbDoG.mat;
%  FR = FbApply2d( X, FB, 'same', 1 );
%
% See also CONV2, FBMAKE
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<3 || isempty(shape)); shape = 'full'; end
if( nargin<4 || isempty(show)); show=0; end

nd=ndims(I);  ndf=ndims(FB);  nf=size(FB,3);
if( nd~=2  ); error('I must be an MxN array'); end
if( ndf~=2 && ndf~=3 ); error('FB must be an MxN or MxNxK array'); end
if( ~isa(I,'double')); I = double(I); end

% apply each filter to image
if( ndf==2 )
  FR = conv2( I, FB, shape );
else
  FR = repmat( conv2(I,FB(:,:,1),shape), [1 1 nf] );
  for i=2:nf; FR(:,:,i)=conv2(I,FB(:,:,i),shape); end
end

% optionally display
if( show )
  figure(show); im(I);
  figure(show+1); montage2(FB,struct('extraInfo',1));
  figure(show+2); montage2(FR,struct('extraInfo',1));
end

author: zouchuhang

repo: LayoutNet


text:
function obj_io_test04 ( output_file_name )

%*****************************************************************************80
%
%% TEST04 tests OBJ_WRITE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2008
%
%  Author:
%
%    John Burkardt
%
  face_num = 12;
  node_num = 8;
  normal_num = 0;
  order_max = 3;

  face_node = [ ...
    1, 3, 2; ...
    2, 3, 4; ...
    1, 6, 5; ...
    1, 2, 6; ...
    3, 7, 4; ...
    4, 7, 8; ...
    5, 6, 8; ...
    5, 8, 7; ...
    1, 5, 7; ...
    1, 7, 3; ...
    2, 4, 6; ...
    6, 4, 8 ]';
  face_order = [ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3 ];
  node_xyz = [ ...
    0.0, 0.0, 0.0; ...
    1.0, 0.0, 0.0; ...
    0.0, 1.0, 0.0; ...
    1.0, 1.0, 0.0; ...
    0.0, 0.0, 1.0; ...
    1.0, 0.0, 1.0; ...
    0.0, 1.0, 1.0; ...
    1.0, 1.0, 1.0 ]';
  normal_vector = [];
  vertex_normal = [];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04\n' );
  fprintf ( 1, '  OBJ_WRITE writes an ASCII OBJ file.\n' );
  fprintf ( 1, '  Here, we do NOT supply any normal vectors.\n' );

  obj_write ( output_file_name, node_num, face_num, normal_num, ...
    order_max, node_xyz, face_order, face_node, normal_vector, vertex_normal );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics data was written to the OBJ file "%s".\n', ...
    output_file_name );

  return
end

repo: jburkardt-m


text:
%TESTP Error estimation of Parzen classifier
%
% 	E = TESTP(A,H,T)
% 	E = TESTP(A,H)
%
% INPUT
%     A    input dataset
%     H    matrix smoothing parameters (optional, def: determined via
%            parzenc)
%     T    test dataset (optional)
%
% OUTPUT
%     E    estimated error rate
%
% DESCRIPTION
% Tests a dataset T on dataset A using a Parzen classification and returns
% the classification error E.  Returns the leave-one-out error estimate. If
% H is not given, it is determined by PARZENC.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, PARZEN_MAP, PARZENML, PARZENC.

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% may be to be merged with parzen_map, see also testk

% $Id: testp.m,v 1.4 2009/09/08 21:27:51 duin Exp $

function [e,d] = testp(a,h,t)

	isvaldfile(a,2,2);
	a = testdatasize(a);
	a = testdatasize(a,'objects');

	if nargin < 2, [W,h] = parzenc(a); end
	[m,k,c] = getsize(a);
	nlab = getnlab(a);
	lablist = getlablist(a);
	p = getprior(a);

	if length(h) == 1, h = h*ones(1,c); end
	if length(h) ~= c, error('Wrong number of smoothing parameters'); end

	% if no test dataset is specified
	if nargin <= 2
		% find for each sample cross-validated estimate
		% of aposteriori probability.
		d = classp(a,nlab,h,p);
		[dmax,J] = max(d',[],1);
		e = nlabcmp(lablist(J,:),lablist(nlab,:)) / m;
        % if the validation dataset is given
	elseif nargin == 3
		lablistt = getlablist(t);
		[n,kt] = size(t);
		nlabt = getnlab(t);
		if k ~= kt
			error('Data sizes do not match');
		end
		d = classp(a,nlab,h,p,t);
		[dmax,J] = max(d',[],1);
		e = nlabcmp(lablistt(J,:),lablistt(nlabt,:)) / n;
	end
return

%CLASSP estimate of Parzen density (if t is not specified,
% a leave-one-out error estimate for a is returned)
function F = classp(a,nlab,h,p,t)

	[m,k] = size(a);
	maxa = max(max(abs(a)));
	a = a/maxa;
	h = h/maxa;
	if nargin < 5
		mt = m;
	else
		[mt,kt] = size(t);
		t = t/maxa;
	end

	c = max(nlab);
	alf=sqrt(2*pi)^k; % density normalization factor
	[num,n] = prmem(mt,m); % use batches to avoid excessive memory usage
	F = ones(mt,c);
  s = sprintf('Compute distance matrix in %i batches: ',num);
  prwaitbar(num,s);
	for i = 0:num-1
    prwaitbar(num,i+1,[s int2str(i+1)]);
		if i == num-1
			nn = mt - num*n + n;
		else
			nn = n;
		end
		range = [i*n+1:i*n+nn];
		if nargin <= 4
			D = +distm(a,a(range,:));
			D(i*n+1:m+1:i*n+nn*m) = inf*ones(1,nn); % set distances to itself at inf
		else
			D = +distm(a,t(range,:));
		end
		for i=1:c
			I = find(nlab == i);
			if length(I) > 0
				F(range,i) = p(i).*sum(exp(-D(I,:)*0.5./(h(i).^2)),1)'./(length(I)*alf*h(i)^k);
			end
		end
  end
  prwaitbar(0);
	F = F + realmin;
	F = F ./ (sum(F')'*ones(1,c));

return

author: marianux

repo: ecg-kit