text:
function C=multiply4X5And5X5Matrices(A,B)
%%MULTIPLY4X4And5X5Matrices This just multiplies the 4X4 matrix A by the
%           5X5 matrix B to get a 4X5 matrix C. This demonstrates the
%           algorithm developed in [1], which minimizes the number of
%           scalar multiplication operations. This could be used as a
%           template for efficient implementation in other programming
%           languages. In Matlab, this is not faster than the built-in
%           matrix multiplication operation.
%
%INPUTS: A A 4X5 matrix.
%        B A 5X5 matrix.
%
%OUTPUTS: C A 4X5 matrix.
%
%EXAMPLE:
%This just shows that this function produces the same result as Matlab
%within finite precision limitations.
% A=randn(4,5);
% B=randn(5,5);
% C=multiply4X5And5X5Matrices(A,B);
% C1=A*B;%Matlab's way.
% RelErr=max(max(abs(abs((C-C1)./C1))))
%
%REFERENCES:
%[1] A. Fawzi, M. Balog, A. Huang, T. Hubert, B. Romera-Paredes,
%    M. Barekatain, A. Novikov, F. J. R. Ruiz, J. Schrittwieser, G.
%    Swirszcz, D. Silver, D. Hassabis, and P. Kohli, "Discovering faster
%    matrix multiplication algorithms with reinforcement learning," Nature,
%    vol. 610, pp. 47-53, 2022.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

h1=A(3,2)*(-B(2,1)-B(2,5)-B(3,1));
h2=(A(2,2)+A(2,5)-A(3,5))*(-B(2,5)-B(5,1));
h3=(-A(3,1)-A(4,1)+A(4,2))*(-B(1,1)+B(2,5));
h4=(A(1,2)+A(1,4)+A(3,4))*(-B(2,5)-B(4,1));
h5=(A(1,5)+A(2,2)+A(2,5))*(-B(2,4)+B(5,1));
h6=(-A(2,2)-A(2,5)-A(4,5))*(B(2,3)+B(5,1));
h7=(-A(1,1)+A(4,1)-A(4,2))*(B(1,1)+B(2,4));
h8=(A(3,2)-A(3,3)-A(4,3))*(-B(2,3)+B(3,1));
h9=(-A(1,2)-A(1,4)+A(4,4))*(B(2,3)+B(4,1));
h10=(A(2,2)+A(2,5))*B(5,1);
h11=(-A(2,1)-A(4,1)+A(4,2))*(-B(1,1)+B(2,2));
h12=(A(4,1)-A(4,2))*B(1,1);
h13=(A(1,2)+A(1,4)+A(2,4))*(B(2,2)+B(4,1));
h14=(A(1,3)-A(3,2)+A(3,3))*(B(2,4)+B(3,1));
h15=(-A(1,2)-A(1,4))*B(4,1);
h16=(-A(3,2)+A(3,3))*B(3,1);
h17=(A(1,2)+A(1,4)-A(2,1)+A(2,2)-A(2,3)+A(2,4)-A(3,2)+A(3,3)-A(4,1)+A(4,2))*B(2,2);
h18=A(2,1)*(B(1,1)+B(1,2)+B(5,2));
h19=-A(2,3)*(B(3,1)+B(3,2)+B(5,2));
h20=(-A(1,5)+A(2,1)+A(2,3)-A(2,5))*(-B(1,1)-B(1,2)+B(1,4)-B(5,2));
h21=(A(2,1)+A(2,3)-A(2,5))*B(5,2);
h22=(A(1,3)-A(1,4)-A(2,4))*(B(1,1)+B(1,2)-B(1,4)-B(3,1)-B(3,2)+B(3,4)+B(4,4));
h23=A(1,3)*(-B(3,1)+B(3,4)+B(4,4));
h24=A(1,5)*(-B(4,4)-B(5,1)+B(5,4));
h25=-A(1,1)*(B(1,1)-B(1,4));
h26=(-A(1,3)+A(1,4)+A(1,5))*B(4,4);
h27=(A(1,3)-A(3,1)+A(3,3))*(B(1,1)-B(1,4)+B(1,5)+B(3,5));
h28=-A(3,4)*(-B(3,5)-B(4,1)-B(4,5));
h29=A(3,1)*(B(1,1)+B(1,5)+B(3,5));
h30=(A(3,1)-A(3,3)+A(3,4))*B(3,5);
h31=(-A(1,4)-A(1,5)-A(3,4))*(-B(4,4)-B(5,1)+B(5,4)-B(5,5));
h32=(A(2,1)+A(4,1)+A(4,4))*(B(1,3)-B(4,1)-B(4,2)-B(4,3));
h33=A(4,3)*(-B(3,1)-B(3,3));
h34=A(4,4)*(-B(1,3)+B(4,1)+B(4,3));
h35=-A(4,5)*(B(1,3)+B(5,1)+B(5,3));
h36=(A(2,3)-A(2,5)-A(4,5))*(B(3,1)+B(3,2)+B(3,3)+B(5,2));
h37=(-A(4,1)-A(4,4)+A(4,5))*B(1,3);
h38=(-A(2,3)-A(3,1)+A(3,3)-A(3,4))*(B(3,5)+B(4,1)+B(4,2)+B(4,5));
h39=(-A(3,1)-A(4,1)-A(4,4)+A(4,5))*(B(1,3)+B(5,1)+B(5,3)+B(5,5));
h40=(-A(1,3)+A(1,4)+A(1,5)-A(4,4))*(-B(3,1)-B(3,3)+B(3,4)+B(4,4));
h41=(-A(1,1)+A(4,1)-A(4,5))*(B(1,3)+B(3,1)+B(3,3)-B(3,4)+B(5,1)+B(5,3)-B(5,4));
h42=(-A(2,1)+A(2,5)-A(3,5))*(-B(1,1)-B(1,2)-B(1,5)+B(4,1)+B(4,2)+B(4,5)-B(5,2));
h43=A(2,4)*(B(4,1)+B(4,2));
h44=(A(2,3)+A(3,2)-A(3,3))*(B(2,2)-B(3,1));
h45=(-A(3,3)+A(3,4)-A(4,3))*(B(3,5)+B(4,1)+B(4,3)+B(4,5)+B(5,1)+B(5,3)+B(5,5));
h46=-A(3,5)*(-B(5,1)-B(5,5));
h47=(A(2,1)-A(2,5)-A(3,1)+A(3,5))*(B(1,1)+B(1,2)+B(1,5)-B(4,1)-B(4,2)-B(4,5));
h48=(-A(2,3)+A(3,3))*(B(2,2)+B(3,2)+B(3,5)+B(4,1)+B(4,2)+B(4,5));
h49=(-A(1,1)-A(1,3)+A(1,4)+A(1,5)-A(2,1)-A(2,3)+A(2,4)+A(2,5))*(-B(1,1)-B(1,2)+B(1,4));
h50=(-A(1,4)-A(2,4))*(B(2,2)-B(3,1)-B(3,2)+B(3,4)-B(4,2)+B(4,4));
h51=A(2,2)*(B(2,1)+B(2,2)-B(5,1));
h52=A(4,2)*(B(1,1)+B(2,1)+B(2,3));
h53=-A(1,2)*(-B(2,1)+B(2,4)+B(4,1));
h54=(A(1,2)+A(1,4)-A(2,2)-A(2,5)-A(3,2)+A(3,3)-A(4,2)+A(4,3)-A(4,4)-A(4,5))*B(2,3);
h55=(A(1,4)-A(4,4))*(-B(2,3)+B(3,1)+B(3,3)-B(3,4)+B(4,3)-B(4,4));
h56=(A(1,1)-A(1,5)-A(4,1)+A(4,5))*(B(3,1)+B(3,3)-B(3,4)+B(5,1)+B(5,3)-B(5,4));
h57=(-A(3,1)-A(4,1))*(-B(1,3)-B(1,5)-B(2,5)-B(5,1)-B(5,3)-B(5,5));
h58=(-A(1,4)-A(1,5)-A(3,4)-A(3,5))*(-B(5,1)+B(5,4)-B(5,5));
h59=(-A(3,3)+A(3,4)-A(4,3)+A(4,4))*(B(4,1)+B(4,3)+B(4,5)+B(5,1)+B(5,3)+B(5,5));
h60=(A(2,5)+A(4,5))*(B(2,3)-B(3,1)-B(3,2)-B(3,3)-B(5,2)-B(5,3));
h61=(A(1,4)+A(3,4))*(B(1,1)-B(1,4)+B(1,5)-B(2,5)-B(4,4)+B(4,5)-B(5,1)+B(5,4)-B(5,5));
h62=(A(2,1)+A(4,1))*(B(1,2)+B(1,3)+B(2,2)-B(4,1)-B(4,2)-B(4,3));
h63=(-A(3,3)-A(4,3))*(-B(2,3)-B(3,3)-B(3,5)-B(4,1)-B(4,3)-B(4,5));
h64=(A(1,1)-A(1,3)-A(1,4)+A(3,1)-A(3,3)-A(3,4))*(B(1,1)-B(1,4)+B(1,5));
h65=(-A(1,1)+A(4,1))*(-B(1,3)+B(1,4)+B(2,4)-B(5,1)-B(5,3)+B(5,4));
h66=(A(1,1)-A(1,2)+A(1,3)-A(1,5)-A(2,2)-A(2,5)-A(3,2)+A(3,3)-A(4,1)+A(4,2))*B(2,4);
h67=(A(2,5)-A(3,5))*(B(1,1)+B(1,2)+B(1,5)-B(2,5)-B(4,1)-B(4,2)-B(4,5)+B(5,2)+B(5,5));
h68=(A(1,1)+A(1,3)-A(1,4)-A(1,5)-A(4,1)-A(4,3)+A(4,4)+A(4,5))*(-B(3,1)-B(3,3)+B(3,4));
h69=(-A(1,3)+A(1,4)-A(2,3)+A(2,4))*(-B(2,4)-B(3,1)-B(3,2)+B(3,4)-B(5,2)+B(5,4));
h70=(A(2,3)-A(2,5)+A(4,3)-A(4,5))*(-B(3,1)-B(3,2)-B(3,3));
h71=(-A(3,1)+A(3,3)-A(3,4)+A(3,5)-A(4,1)+A(4,3)-A(4,4)+A(4,5))*(-B(5,1)-B(5,3)-B(5,5));
h72=(-A(2,1)-A(2,4)-A(4,1)-A(4,4))*(B(4,1)+B(4,2)+B(4,3));
h73=(A(1,3)-A(1,4)-A(1,5)+A(2,3)-A(2,4)-A(2,5))*(B(1,1)+B(1,2)-B(1,4)+B(2,4)+B(5,2)-B(5,4));
h74=(A(2,1)-A(2,3)+A(2,4)-A(3,1)+A(3,3)-A(3,4))*(B(4,1)+B(4,2)+B(4,5));
h75=-(A(1,2)+A(1,4)-A(2,2)-A(2,5)-A(3,1)+A(3,2)+A(3,4)+A(3,5)-A(4,1)+A(4,2))*B(2,5);
h76=(A(1,3)+A(3,3))*(-B(1,1)+B(1,4)-B(1,5)+B(2,4)+B(3,4)-B(3,5));

C=zeros(4,5);
C(1,1)=-h10+h12+h14-h15-h16+h53+h5-h66-h7;
C(2,1)=h10+h11-h12+h13+h15+h16-h17-h44+h51;
C(3,1)=h10-h12+h15+h16-h1+h2+h3-h4+h75;
C(4,1)=-h10+h12-h15-h16+h52+h54-h6-h8+h9;
C(1,2)=h13+h15+h20+h21-h22+h23+h25-h43+h49+h50;
C(2,2)=-h11+h12-h13-h15-h16+h17+h18-h19-h21+h43+h44;
C(3,2)=-h16-h19-h21-h28-h29-h38+h42+h44-h47+h48;
C(4,2)=h11-h12-h18+h21-h32+h33-h34-h36+h62-h70;
C(1,3)=h15+h23+h24+h34-h37+h40-h41+h55-h56-h9;
C(2,3)=-h10+h19+h32+h35+h36+h37-h43-h60-h6-h72;
C(3,3)=-h16-h28+h33+h37-h39+h45-h46+h63-h71-h8;
C(4,3)=h10+h15+h16-h33+h34-h35-h37-h54+h6+h8-h9;
C(1,4)=-h10+h12+h14-h16+h23+h24+h25+h26+h5-h66-h7;
C(2,4)=h10+h18-h19+h20-h22-h24-h26-h5-h69+h73;
C(3,4)=-h14+h16-h23-h26+h27+h29+h31+h46-h58+h76;
C(4,4)=h12+h25+h26-h33-h35-h40+h41+h65-h68-h7;
C(1,5)=h15+h24+h25+h27-h28+h30+h31-h4+h61+h64;
C(2,5)=-h10-h18-h2-h30-h38+h42-h43+h46+h67+h74;
C(3,5)=-h10+h12-h15+h28+h29-h2-h30-h3+h46+h4-h75;
C(4,5)=-h12-h29+h30-h34+h35+h39+h3-h45+h57+h59;

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [samples, logp, diagn] = hmc_nuts(f, theta0, opt)
%HMC_NUTS No-U-Turn Sampler (NUTS)
%
%  Description
%    [SAMPLES, LOGP, DIAGN] = HMC_NUTS(f, theta0, opt)
%    Implements the No-U-Turn Sampler (NUTS), specifically,
%    algorithm 6 from the NUTS paper (Hoffman & Gelman, 2011). Runs
%    opt.Madapt steps of burn-in, during which it adapts the step
%    size parameter epsilon, then starts generating samples to
%    return.
%
%    f(theta) should be a function that returns the log probability its
%    gradient evaluated at theta. I.e., you should be able to call
%    [logp grad] = f(theta).
%
%    opt.epsilon is a step size parameter.
%    opt.M       is the number of samples to generate.
%    opt.Madapt  is the number of steps of burn-in/how long to run
%                the dual averaging algorithm to fit the step size
%                epsilon. Note that there is no need to provide
%                opt.epsilon if doing adaptation.
%    opt.theta0  is a 1-by-D vector with the desired initial setting
%                of the parameters.
%    opt.delta   should be between 0 and 1, and is a target HMC
%                acceptance probability. Defaults to 0.8 if
%                unspecified.
%
%
%    The returned variable "samples" is an (M+Madapt)-by-D matrix
%    of samples generated by NUTS, including burn-in samples.
%
%    Note that when used from gp_mc, opt.M and opt.Madapt are both 0 or
%    1 (hmc_nuts returns only one sample to gp_mc). Number of epsilon
%    adaptations should be set in hmc options structure hmc_opt.nadapt, in
%    gp_mc(... ,'hmc_opt', hmc_opt).
%
%    The returned structure diagn includes step-size vector
%    epsilon, number of rejected samples and dual averaging
%    parameters so its possible to continue adapting step-size
%    parameter.

%      Copyright (c) 2011, Matthew D. Hoffman
%      Copyright (c) 2012, Ville Tolvanen
%      All rights reserved.
%

% Redistribution and use in source and binary forms, with or
% without modification, are permitted provided that the following
% conditions are met:
%
% Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
%
% Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the
% distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
% CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
% INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
% SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
% LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
% USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
% AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
% LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
% ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

global nfevals;
nfevals = 0;

if ~isfield(opt, 'delta')
  delta = 0.8;
else
  delta = opt.delta;
end
if ~isfield(opt, 'M')
  M = 1;
else
  M = opt.M;
end
if ~isfield(opt, 'Madapt')
  Madapt = 0;
else
  Madapt = opt.Madapt;
end

diagn.rej = 0;
assert(size(theta0, 1) == 1);

D = length(theta0);
samples = zeros(M+Madapt, D);

[logp grad] = f(theta0);
samples(1, :) = theta0;

% Parameters to the dual averaging algorithm.
gamma = 0.05;
t0 = 10;
kappa = 0.75;
% Initialize dual averaging algorithm.
epsilonbar = 1;
Hbar = 0;

if isfield(opt, 'epsilon') && ~isempty(opt.epsilon)
  epsilon = opt.epsilon(end);

  % Hbar & epsilonbar are needed when doing adaptation of step-length
  if isfield(opt, 'Hbar') && ~isempty(opt.Hbar)
    Hbar = opt.Hbar;
  end
  if isfield(opt, 'epsilonbar') && ~isempty(opt.epsilonbar)
    epsilonbar=opt.epsilonbar;
  else
    epsilonbar=opt.epsilon;
  end
  mu = log(10*opt.epsilon(1));
else
  % Choose a reasonable first epsilon by a simple heuristic.
  epsilon = find_reasonable_epsilon(theta0, grad, logp, f);
  mu = log(10*epsilon);

  opt.epsilon = epsilon;
  opt.epsilonbar = epsilonbar;
  opt.Hbar = Hbar;
end

for m = 2:M+Madapt+1,
%     m
    % Resample momenta.
    r0 = randn(1, D);
    % Joint log-probability of theta and momenta r.
    joint = logp - 0.5 * (r0 * r0');
    % Resample u ~ uniform([0, exp(joint)]).
    % Equivalent to (log(u) - joint) ~ exponential(1).
    logu = joint - exprnd(1);
    % Initialize tree.
    thetaminus = samples(m-1, :);
    thetaplus = samples(m-1, :);
    rminus = r0;
    rplus = r0;
    gradminus = grad;
    gradplus = grad;
    % Initial height j = 0.
    j = 0;
    % If all else fails, the next sample is the previous sample.
    samples(m, :) = samples(m-1, :);
    % Initially the only valid point is the initial point.
    n = 1;
    rej = 0;
    % Main loop---keep going until the criterion s == 0.
    s = 1;
    while (s == 1)
        % Choose a direction. -1=backwards, 1=forwards.
        v = 2*(rand() < 0.5)-1;
        % Double the size of the tree.
        if (v == -1)
            [thetaminus, rminus, gradminus, tmp, tmp, tmp, thetaprime, gradprime, logpprime, nprime, sprime, alpha, nalpha] = ...
                build_tree(thetaminus, rminus, gradminus, logu, v, j, epsilon, f, joint);
        else
            [tmp, tmp, tmp, thetaplus, rplus, gradplus, thetaprime, gradprime, logpprime, nprime, sprime, alpha, nalpha] = ...
                build_tree(thetaplus, rplus, gradplus, logu, v, j, epsilon, f, joint);
        end
        % Use Metropolis-Hastings to decide whether or not to move to a
        % point from the half-tree we just generated.
        if ((sprime == 1) && (rand() < nprime/n))
          samples(m, :) = thetaprime;
          logp = logpprime;
          grad = gradprime;
        else
          rej = rej + 1;
        end
        % Update number of valid points we've seen.
        n = n + nprime;
        % Decide if it's time to stop.
        s = sprime && stop_criterion(thetaminus, thetaplus, rminus, rplus);
        % Increment depth.
        j = j + 1;
    end

    % Do adaptation of epsilon if we're still doing burn-in.
    eta = 1 / (length(opt.epsilon) + t0);
    Hbar = (1 - eta) * Hbar + eta * (delta - alpha / nalpha);
    if (m <= Madapt+1)
        epsilon = exp(mu - sqrt(m-1)/gamma * Hbar);
        eta = (length(opt.epsilon))^-kappa;
        epsilonbar = exp((1 - eta) * log(epsilonbar) + eta * log(epsilon));
    else
        epsilon = epsilonbar;
    end
    opt.epsilon(end+1) = epsilon;
    opt.epsilonbar = epsilonbar;
    opt.Hbar = Hbar;
    diagn.rej = diagn.rej + rej;
end

diagn.opt = opt;
end

function [thetaprime, rprime, gradprime, logpprime] = leapfrog(theta, r, grad, epsilon, f)
rprime = r + 0.5 * epsilon * grad;
thetaprime = theta + epsilon * rprime;
[logpprime, gradprime] = f(thetaprime);
rprime = rprime + 0.5 * epsilon * gradprime;
global nfevals;
nfevals = nfevals + 1;
end

function criterion = stop_criterion(thetaminus, thetaplus, rminus, rplus)
thetavec = thetaplus - thetaminus;
criterion = (thetavec * rminus' >= 0) && (thetavec * rplus' >= 0);
end

% The main recursion.
function [thetaminus, rminus, gradminus, thetaplus, rplus, gradplus, thetaprime, gradprime, logpprime, nprime, sprime, alphaprime, nalphaprime] = ...
                build_tree(theta, r, grad, logu, v, j, epsilon, f, joint0)
if (j == 0)
    % Base case: Take a single leapfrog step in the direction v.
    [thetaprime, rprime, gradprime, logpprime] = leapfrog(theta, r, grad, v*epsilon, f);
    joint = logpprime - 0.5 * (rprime * rprime');
    % Is the new point in the slice?
    nprime = logu < joint;
    % Is the simulation wildly inaccurate?
    sprime = logu - 1000 < joint;
    % Set the return values---minus=plus for all things here, since the
    % "tree" is of depth 0.
    thetaminus = thetaprime;
    thetaplus = thetaprime;
    rminus = rprime;
    rplus = rprime;
    gradminus = gradprime;
    gradplus = gradprime;
    % Compute the acceptance probability.
    alphaprime = min(1, exp(logpprime - 0.5 * (rprime * rprime') - joint0));
    nalphaprime = 1;
else
    % Recursion: Implicitly build the height j-1 left and right subtrees.
    [thetaminus, rminus, gradminus, thetaplus, rplus, gradplus, thetaprime, gradprime, logpprime, nprime, sprime, alphaprime, nalphaprime] = ...
                build_tree(theta, r, grad, logu, v, j-1, epsilon, f, joint0);
    % No need to keep going if the stopping criteria were met in the first
    % subtree.
    if (sprime == 1)
        if (v == -1)
            [thetaminus, rminus, gradminus, tmp, tmp, tmp, thetaprime2, gradprime2, logpprime2, nprime2, sprime2, alphaprime2, nalphaprime2] = ...
                build_tree(thetaminus, rminus, gradminus, logu, v, j-1, epsilon, f, joint0);
        else
            [tmp, tmp, tmp, thetaplus, rplus, gradplus, thetaprime2, gradprime2, logpprime2, nprime2, sprime2, alphaprime2, nalphaprime2] = ...
                build_tree(thetaplus, rplus, gradplus, logu, v, j-1, epsilon, f, joint0);
        end
        % Choose which subtree to propagate a sample up from.
        if (rand() < nprime2 / (nprime + nprime2))
            thetaprime = thetaprime2;
            gradprime = gradprime2;
            logpprime = logpprime2;
        end
        % Update the number of valid points.
        nprime = nprime + nprime2;
        % Update the stopping criterion.
        sprime = sprime && sprime2 && stop_criterion(thetaminus, thetaplus, rminus, rplus);
        % Update the acceptance probability statistics.
        alphaprime = alphaprime + alphaprime2;
        nalphaprime = nalphaprime + nalphaprime2;
    end
end
end

function epsilon = find_reasonable_epsilon(theta0, grad0, logp0, f)
epsilon = 0.1;
r0 = randn(1, length(theta0));
% Figure out what direction we should be moving epsilon.
[tmp, rprime, tmp, logpprime] = leapfrog(theta0, r0, grad0, epsilon, f);
acceptprob = exp(logpprime - logp0 - 0.5 * (rprime * rprime' - r0 * r0'));

% Here we presume that energy function returns NaN, if energy cannot be
% evaluated at the suggested hyperparameters so that we need smalled epsilon
if isnan(acceptprob)
  acceptprob=0;
end
a = 2 * (acceptprob > 0.5) - 1;
% Keep moving epsilon in that direction until acceptprob crosses 0.5.
while (acceptprob^a > 2^(-a))
    epsilon = epsilon * 2^a;
    [tmp, rprime, tmp, logpprime] = leapfrog(theta0, r0, grad0, epsilon, f);
    acceptprob = exp(logpprime - logp0 - 0.5 * (rprime * rprime' - r0 * r0'));
end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function ncV=grComp(E,n)
% Function ncV=grComp(E,n) find all components of the graph.
% Input parameter:
%   E(m,2) - the edges of graph;
%     1st and 2nd elements of each row is numbers of vertexes;
%     m - number of edges.
%   n - number of vertexes ( optional, by default n=max(max(E)) ).
%     This input parameter is needed, if last vertexes is isolated.
% Output parameter:
%   ncV(n,1) - the the vector-column with the number of component
%     for each vertex;
% Author: Sergiy Iglin
% e-mail: siglin@yandex.ru
% personal page: http://iglin.exponenta.ru

if nargin<1,
  error('There are no input data!')
end
[m,n1,E1] = grValidation(E); % data validation
E2=[E1(:,1:2);E1(:,[2 1])]; % all arrows and vice versa
[Dec,Ord]=grDecOrd(E2); % the components
ncV=sum(Dec*diag([1:size(Dec,2)]),2); % the numbers of components
if (nargin>1)&(n>n1), % last isolated vertexes
  ncV=[ncV;[1:n-n1]'+max(ncV)];
end
return

author: HuangCongQing


text:
function [vFrequency, vAmplitude] = fastfft(vData, SampleRate, Plot)

%FASTFFT   Create useful data from an FFT operation.
%   Usage: [vFrequency, vAmplitude] = fastfft(vData, SampleRate, [Plot])
%
%   (no plot will be shown if the last input == 0 or is not included)
%
%   This function inputs 'vData' as a vector (row or column),
%   'SampleRate' as a number (samples/sec), 'Plot' as anything,
%   and does the following:
%
%     1: Removes the DC offset of the data
%     2: Puts the data through a hanning window
%     3: Calculates the Fast Fourier Transform (FFT)
%     4: Calculates the amplitude from the FFT
%     5: Calculates the frequency scale
%     6: Optionally creates a Bode plot
%
%   Created 7/22/03, Rick Auch, mekaneck@campbellsville.com

%Make vData a row vector
if size(vData,2)==1
    vData = vData';
end

%Calculate number of data points in data
n = length(vData);

%Remove DC Offset
vData = vData - mean(vData);

%Put data through hanning window using hanning subfunction
vData = hanning(vData);

%Calculate FFT
vData = fft(vData);

%Calculate amplitude from FFT (multply by sqrt(8/3) because of effects of hanning window)
vAmplitude = abs(vData)*sqrt(8/3);

%Calculate frequency scale
vFrequency = linspace(0,n-1,n)*(SampleRate/n);

%Limit both output vectors due to Nyquist criterion
DataLimit = ceil(n/2);
vAmplitude = vAmplitude(1:DataLimit);
vFrequency = vFrequency(1:DataLimit);

if exist('Plot', 'var')==1 & Plot~=0
    plot(vFrequency, vAmplitude);
    title('Bode Plot');
    xlabel('Frequency (Hz)');
    ylabel('Amplitude');
end

%------------------------------------------------------------------------------------------
%Hanning Subfunction
function vOutput = hanning(vInput)
% This function takes a vector input and outputs the same vector,
% multiplied by the hanning window function

%Determine the number of input data points
n = length(vInput);

%Initialize the vector
vHanningFunc = linspace(0,n-1,n);

%Calculate the hanning funtion
vHanningFunc = .5*(1-cos(2*pi*vHanningFunc/(n-1)));

%Output the result
vOutput = vInput.*vHanningFunc;

author: Sable

repo: mcbench-benchmarks


text:
function outputImage = zoom(this, zoom_factor)
% Geometric shear of MrImage
%
%   Y = MrImage()
%   Y.zoom(inputs)
%
% This is a method of class MrImage.
%
% NOTE: This is a method of MrImage rather than MrImageGeometry, because
%       the latter is composed on the fly from affineTransformation and
%       dimInfo to integrate both information and sustain consistency.
%       Thus, the effect will only be visible using the world space plot
%       options, i.e. Y.plot('plotType', 'spmi').
% NOTE: The rotation is *added* to the rotation already defined in the
%       affine transformation matrix, equivalent to the rotation in SPM
%       Display.
% IN
%   zoom_factor            [1,3] zoom factor, i.e., [zx, zy, zz]
%
% OUT
%
% EXAMPLE
%   zoom
%
%   See also MrImage

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2019-12-09
% Copyright (C) 2019 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

outputImage = this.copyobj();

[outputImage.dimInfo, outputImage.affineTransformation] = ...
    this.geometry.perform_world_space_operation('zoom', zoom_factor, this.dimInfo);
end

author: translationalneuromodeling

repo: tapas


text:
function variance = f_variance ( m, n )

%*****************************************************************************80
%
%% F_VARIANCE returns the variance of the F central PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the parameters of the PDF.
%    1 <= M,
%    1 <= N.
%    Note, however, that the variance is not defined unless 5 <= N.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  if ( n < 5 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'F_VARIANCE - Fatal error!\n' );
    fprintf ( 1, '  The variance is not defined for N < 5.\n' );
    error ( 'F_VARIANCE - Fatal error!' );
  end

  variance = 2 * n * n * ( m + n - 2 ) / ( m * ( n - 2 )^2 * ( n - 4 ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%BCSFRONTIER BlueChipStock rolling efficient frontiers

addpath ./source

load BlueChipBacktest0
load BlueChipBacktest

% Plot 3D efficient frontiers

figure(1);
surf(X0,Y0,Z0,'FaceColor','interp','EdgeColor','none','FaceLighting','phong');
ylabel('\bfStd.Dev. of Returns');
zlabel('\bfMean of Returns');
title('\bfRolling Efficient Frontiers (Absolute Total Return)');
camlight right;
view(30,30);

i = input('Continue >');

figure(1);
surf(X,Y,Z,'FaceColor','interp','EdgeColor','none','FaceLighting','phong');
ylabel('\bfStd.Dev. of Returns');
zlabel('\bfMean of Returns');
title('\bfRolling Efficient Frontiers (Relative Total Return vs DJIA)');
camlight right;
view(30,30);


text:
function r4lib_test ( )

%*****************************************************************************80
%
%% R4LIB_TEST tests the R4LIB library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R4LIB_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the R4LIB library.\n' );

  r4lib_test001 ( );
  r4lib_test023 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R4LIB_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
%requires https://uk.mathworks.com/matlabcentral/fileexchange/29774-stoichiometry-tools

%             TODO
%             Error using parse_formula>parse_formula_ (line 191)
%             Could not parse formula:
%             O5C5H3A2
%                   ^^

formula = 'O5C5H3A2';
hillformulaExample = hillformula(formula);
assert(strcmp(hillformulaExample,'C5H3A2O5'))

author: opencobra

repo: cobratoolbox


text:
function [ Grid,varargout ] = CavityGridOperators( N )
%CAVITYGRIDOPERATORS generates the grid and operators
[Grid]=CollocationGrid_q(N);    % grid coordiuantes and indices

if nargout==2
    varargout{1}= CreateOperators_psi( Grid.D,eye(N+1));  %  derivative matrices
end

[~,Grid.wc] = clencurt(N);
 Grid.W = kron(Grid.wc,Grid.wc);   % integration coefficients for computation of Kinetic Energy

end

function [Grid]=CollocationGrid_q(N)

  [Grid.D,Grid.x] = cheb(N);
  [xx,yy]=meshgrid(Grid.x,Grid.x); Grid.xx=xx(:);Grid.yy=yy(:);

% Index positions for walls and interior
  Grid.Walls = find(abs(xx)==1 | abs(yy)==1);
  Grid.LeftWall = find(xx==-1); Grid.RightWall=find(xx==1);
  Grid.BottomWall=find(yy==-1); Grid.TopWall=find(yy==1);
  Grid.Interior = setdiff((1:length(Grid.xx)),Grid.Walls);

  % operators
  I = eye(N+1);
  Grid.D2 = (diag(1-Grid.x.^2)*Grid.D^2 - 4*diag(Grid.x)*Grid.D - 2*I);
  Grid.D4 = (diag(1-Grid.x.^2)*Grid.D^4 - 8*diag(Grid.x)*Grid.D^3 - 12*Grid.D^2); % Trefethen style

     % Propagation function
   Grid.S0 = (1-Grid.xx.^2).*(1-Grid.yy.^2);
   Grid.S1 = Grid.S0;
   Grid.S1(Grid.Walls)=1;

   % plot on cheb grid
   Grid.cplot = @(X) contourf(reshape(xx,N+1,N+1),reshape(yy,N+1,N+1),reshape(X,N+1,N+1),50,'LineStyle','None');

end

  function [x,w] = clencurt(N)
  theta = pi*(0:N)'/N; x = cos(theta);
  w = zeros(1,N+1); ii = 2:N; v = ones(N-1,1);
  if mod(N,2)==0
    w(1) = 1/(N^2-1); w(N+1) = w(1);
    for k=1:N/2-1, v = v - 2*cos(2*k*theta(ii))/(4*k^2-1); end
    v = v - cos(N*theta(ii))/(N^2-1);
  else
    w(1) = 1/N^2; w(N+1) = w(1);
    for k=1:(N-1)/2, v = v - 2*cos(2*k*theta(ii))/(4*k^2-1); end
  end
  w(ii) = 2*v/N;
  end

 function [ Operators_psi ] = CreateOperators_psi( D,I)
%CREATEOPERATORS_q creates the Operator acting on psi(x)
Operators_psi.Dx = kron(D,I);
Operators_psi.Dy = kron(I,D);
D2= D^2; D4=D2^2;
Operators_psi.del2 = kron(D2,I)+kron(I,D2);

Operators_psi.del4 = kron(D4,I)+kron(I,D4)+2*kron(D2,I)*kron(I,D2);

end

author: arbabiha


text:
function M = lbs_matrix(V,W)
  % LBS_MATRIX  construct a matrix that when multiplied against a column of
  % affine transformation entries computes new coordinates of the vertices
  %
  % M = lbs_matrix(V,W)
  %
  % Input:
  %   V  #V by dim list of vertex rest positions
  %   W  #W by #handles list of correspondence weights
  % Output:
  %   M  #V * dim by #handles * dim * (dim+1) matrix such that
  %     new_V(:) = LBS(V,W,A) = reshape(M * A,size(V)), where A is a column
  %     vectors formed by the entries in each handle's dim by dim+1
  %     transformation matrix. Specifcally, A =
  %       reshape(permute(Astack,[3 1 2]),n*dim*(dim+1),1)
  %     or A = [Lxx;Lyx;Lxy;Lyy;tx;ty], and likewise for other dim
  %     if Astack(:,:,i) is the dim by (dim+1) transformation at handle i
  %
  % Example:
  %  MLBS = lbs_matrix(V,W);
  %  n = size(V,1);
  %  dim = size(V,2);
  %  m = size(W,2);
  %  % stack of identity transformations
  %  Astack = repmat([eye(dim,dim) zeros(dim,1)],[1 1 m]);
  %  % collect transformations into column
  %  A = reshape(permute(Astack,[3 1 2]),m*dim*(dim+1),1);
  %  % apply transformations
  %  new_V = MLBS*A;
  %  new_V = reshape(new_V,[n dim]);
  %
  %  % Alternative:
  %  % Q is #W by 4 list of quats, T is #W by 3 list of translations
  %  % A is #W*4 by 3 stack of transposed affine transformations
  %  A = reshape(cat(2,permute(quat2mat(Q),[2 1 3]),permute(T,[2 3 1])),3,[])';
  %  % M is #V by #W*4 skinning matrix
  %  % M = (𝟙ᵀ⊗ [V 𝟙]) ⊙ (W ⊗ 𝟙ᵀ)
  %  %M = kron(ones(1,size(W,2)),[V ones(size(V,1),1)]).* ...
  %  %    kron(W,ones(1,size(V,2)+1));
  %  M = reshape([V ones(size(V,1),1)].*permute(W,[1 3 2]),size(V,1),[]);
  %  U = M*A;
  %

  % number of mesh (domain) vertices
  n = size(V,1);
  assert(n == size(W,1));
  % dimension of mesh
  dim = size(V,2);
  % number of handles
  m = size(W,2);

  % M = zeros(V*dim,m*dim*(dim+1));

  % repeat vertex positions so that VV(:,:,i) gives #V by #handles matrix where
  % each column is ith coordinates of vertices
  VV = permute(repmat(V,[1 1 m]),[1 3 2]);
  % multiply each column in VV by respective weights for that handle
  VVW = VV.*repmat(W,[1 1 dim]);
  % matrix of zeros
  Z = zeros(n,m);
  switch dim
  case 2
    M = [ ...
      VVW(:,:,1) Z          VVW(:,:,2) Z          W Z; ...
      Z          VVW(:,:,1) Z          VVW(:,:,2) Z W];
  case 3
    M = [ ...
      VVW(:,:,1) Z          Z          VVW(:,:,2) Z          Z          VVW(:,:,3) Z          Z          W Z Z; ...
      Z          VVW(:,:,1) Z          Z          VVW(:,:,2) Z          Z          VVW(:,:,3) Z          Z W Z; ...
      Z          Z          VVW(:,:,1) Z          Z          VVW(:,:,2) Z          Z          VVW(:,:,3) Z Z W];
  otherwise
    error('Only dim=2 or dim=3 supported');
  end

end

author: alecjacobson

repo: gptoolbox


text:
% You can create a decoupled network and convert it to a format that can be
% viewed in Samiam by running this script.

pedigree = struct('parents', [0,0;1,3;0,0;1,3;2,6;0,0;2,6;4,9;0,0]);
pedigree.names = {'Ira','James','Robin','Eva','Jason','Rene','Benjamin','Sandra','Aaron'};
phenotypeList = {'CysticFibrosis', 'NoCysticFibrosis'};
alleleFreqsThree = [0.1; 0.7; 0.2];
alleleListThree = {'F', 'f', 'n'};
alphaListThree = [0.8; 0.6; 0.1; 0.5; 0.05; 0.01];
positionsGeneCopy = [1040, 600, 1170, 600, 1105, 500; 1300, 400, 1430, 400, 1365, 300; 780, 600, 910, 600, 845, 500; 520, 400, 650, 400, 585, 300; 1560, 200, 1690, 200, 1625, 100; 2080, 400, 2210, 400, 2145, 300; 1820, 200, 1950, 200, 1885, 100; 260, 200, 390, 200, 325, 100; 0, 400, 130, 400, 65, 300];

% This will construct a decoupled Bayesian network and convert it into a
% file that can be viewed in SamIam.

factorListDecoupled = constructDecoupledGeneticNetwork(pedigree, alleleFreqsThree, alphaListThree);
sendToSamiamGeneCopy(pedigree, factorListDecoupled, alleleListThree, phenotypeList, positionsGeneCopy, 'cysticFibrosisBayesNetGeneCopy');

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [D, dist, aLongestString] = LCS(X,Y)
%%%Calculates the longest common substring between to strings.
%%%Code written by David Cumin
%%%email: d.cumin@auckland.ac.nz
%%%INPUT
%%%X, Y - both are strings e.g. 'test' or 'stingtocompare'
%%%OUTPUT
%%%D is the substring over the length of the shortest string
%%%dist is the length of the substring
%%%aLongestString is a sting of length dist (only one of potentially many)

%%%For example
%%% X = 'abcabc';
%%% Y = 'adcbac';
%%% [D dist str] = LCS(X,Y);
%%% results in:
%%% D = 0.6667
%%% dist = 4
%%% str = acbc
%%% this is seen for X: 'a-c-bc' and Y: 'a-cb-c'

%%%Make matrix
n =length(X);
m =length(Y);
L=zeros(n+1,m+1);
L(1,:)=0;
L(:,1)=0;
b = zeros(n+1,m+1);
b(:,1)=1;%%%Up
b(1,:)=2;%%%Left

for i = 2:n+1
    for j = 2:m+1
        if (X(i-1) == Y(j-1))
            L(i,j) = L(i-1,j-1) + 1;
            b(i,j) = 3;%%%Up and left
        else
            L(i,j) = L(i-1,j-1);
        end
        if(L(i-1,j) >= L(i,j))
            L(i,j) = L(i-1,j);
            b(i,j) = 1;%Up
        end
        if(L(i,j-1) >= L(i,j))
            L(i,j) = L(i,j-1);
            b(i,j) = 2;%Left
        end
    end
end
L(:,1) = [];
L(1,:) = [];
b(:,1) = [];
b(1,:) = [];
dist = L(n,m);

D = (dist / min(m,n));
if(dist == 0)
    aLongestString = '';
else
    %%%now backtrack to find the longest subsequence
    i = n;
    j = m;
    p = dist;
    aLongestString = {};
    while(i>0 && j>0)
        if(b(i,j) == 3)
            aLongestString{p} = X(i);
            p = p-1;
            i = i-1;
            j = j-1;
        elseif(b(i,j) == 1)
            i = i-1;
        elseif(b(i,j) == 2)
            j = j-1;
        end
    end

    if ischar(aLongestString{1})
        aLongestString = char(aLongestString)';
    else
        aLongestString = cell2mat(aLongestString);
    end
end

end


text:
function [ x, seed ] = r4vec_uniform_ab ( n, a, b, seed )

%*****************************************************************************80
%
%% R4VEC_UNIFORM_AB returns a scaled pseudorandom R4VEC.
%
%  Discussion:
%
%    Each dimension ranges from A to B.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real A, B, the range of the pseudorandom values.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X(N), the vector of pseudorandom values.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  x = zeros ( n, 1 );

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4VEC_UNIFORM_AB - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R4VEC_UNIFORM_AB - Fatal error!' );
  end

  for i = 1 : n

    k = floor ( seed / 127773 );

    seed = 16807 * ( seed - k * 127773 ) - k * 2836;

    if ( seed < 0 )
      seed = seed + 2147483647;
    end

    x(i) = a + ( b - a ) * seed * 4.656612875E-10;

  end

  return
end


text:
%% FACE_POOL runs the example interactively with MATLABPOOL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 June 2010
%
%  Author:
%
%    Gene Cliff, John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FACE_POOL\n' );
  fprintf ( 1, '  Evaluate FACE_FUN directly with MATLABPOOL.\n' );
%
%  Load the MAT file containing the grid coordinates as the variable
%  "Grid" and the element connectivity as the variable "e_conn".
%
  load GRID
%
%  We will be looking for tetrahedrons which have one face lying
%  on a plane defining one side of the box.  To define this box
%  we must specify:
%
%  II, the index of the "I" variable which is constant on the plane.
%  If the constant variable is X, then II is 1, and so on.
%
%  VAL_II, the value of the constant "I" variable.
%
%  XJ_LB, XJ_UB, the upper and lower bounds for the "J" nonconstant variable.
%
%  XK_LB, XK_UB, the upper and lower bounds for the "K" nonconstant variable.
%
%  The "J" and "K" variables are the first and second variables to follow
%  the "I" variable, assuming that we are counting with wrap around.
%  There are thus three possibilities for I, J and K:
%
%    I  J  K
%   -- -- --
%    1  2  3
%    2  3  1
%    3  1  2
%
%  Our choices below specify that we are looking at the plane Z = 3.0,
%  and more specifically, the box in that plane for which
%  -1.0 <= X <= +1.0, and -1.0 <= Y <= +1.0.
%
  ii = 3;
  val_ii = 3.0;
  xj_lb = -1.0;
  xj_ub = +1.0;
  xk_lb = -1.0;
  xk_ub = +1.0;
%
%  Call FACE_FUN to compute a list of the triangular faces of the
%  tetrahedral elements which lie on the boundary plane.
%
%  Get 4 "labs" to work on the job.
%
  matlabpool open local 4
  F_conn = face_fun ( Grid, e_conn, ii, val_ii, xj_lb, xj_ub, xk_lb, xk_ub );
  matlabpool close

  face_num = size ( F_conn, 2 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of triangles on boundary is %d\n', face_num );
%
%  Display the surface triangulation.
%
  i_pl = setdiff ( 1:3, ii );

  triplot ( F_conn', Grid(:, i_pl(1)), Grid(:, i_pl(2)) );

  if ( ii == 1 )
    title ( 'Triangular faces on X boundary plane' );
    xlabel ( '- Y axis -' );
    ylabel ( '- Z axis -' );
  elseif ( ii == 2 )
    xlabel ( '- X axis -' );
    title ( 'Triangular faces on Y boundary plane' );
    ylabel ( '- Z axis -' );
  elseif ( ii == 3 )
    xlabel ( '- X axis -' );
    ylabel ( '- Y axis -' );
    title ( 'Triangular faces on Z boundary plane' );
  end

  axis equal
  axis tight


text:
function [SDANN, SDNNI] = CalcSDANN(windows_all, tNN, NN, HRVparams)
% [SDANN, SDNNI] = CalcSDANN(windows_all, tNN, NN, HRVparams)
%
%   OVERVIEW:  Calculates SDANN and SDNNindex for all segments
%              of a specified constant length with no overlap.
%
%   INPUT:     windows_all :
%              NN          : a single column of NN (normal normal) interval
%                            data in seconds
%              tNN         : the time indices of the rr interval data (seconds)
%              HRVparams   : settings struct, including segment length - the
%                            length of the interval requested
%                            (it is defaulted at 5minutes == 300 seconds)
%              NOTE: seglength must be specified in the same units as <times>
%
%   OUTPUT:     SDANN : standard deviation of the averages of values
%               SDNNI : mean of the standard deviations of all values
%
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Main script written by Adriana N. Vest
%       Dependent scripts written by various authors
%       (see functions for details)
%	COPYRIGHT (C) 2016
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%%

if nargin<3
    error('not enough input arguments!');
end

segmentlength = HRVparams.windowlength;

for i = 1:length(windows_all)
    if ~isnan(windows_all(i))
        idx = find(tNN >= windows_all(i) & tNN < windows_all(i) + segmentlength);

        nn_win = NN(idx);

        sm(i) = mean(nn_win);     % mean of each segment
        sstd(i) = std(nn_win);    % stdev of each segment
    else
        sm(i) = NaN;
        sstd(i) = NaN;
    end
end

% Remove NaN values from calculation
seg_mean = sm;
seg_stdev = sstd;
idxrem = find(isnan(sm));
seg_mean(idxrem) = [];
seg_stdev(idxrem) = [];

SDANN = std(seg_mean);              % stdev of means of each segment
SDNNI = mean(seg_stdev);            % mean of stdevs from each segment

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function []=panel6disp(n,N,m,p,k,T,d1,d2,d3,d4,d5,d,Ymat,Xtilde,Units,endo,exo,const,Xi,theta_median,theta_std,theta_lbound,theta_ubound,sigma_median,D_estimates,gamma_estimates,alpha0,delta0,gama,a0,b0,rho,psi,acceptrate,startdate,enddate,forecast_record,forecast_estimates,Fcperiods,stringdates3,Fstartdate,Fcenddate,Feval,Fcomp,data_endo_c,IRF,IRFt,pref,names)

% recover a point estimate (the median) of the VAR coefficients
betatilde=Xi*theta_median(:,:,end);
Btilde=reshape(betatilde,k,N*n);

% check whether the model is stationary
[stationary,eigmodulus]=bear.checkstable(betatilde,N*n,p,k);

% estimate the in-sample evaluation criteria

% obtain a point estimate thetatilde of the structural factors, which is the median
Thetatilde=reshape(theta_median,T*d,1);
% compute fitted values
Ytilde=full(reshape(Xtilde*Thetatilde,N*n,T)');
% reshape for convenience
Ytilde=reshape(Ytilde,T,n,N);
Ymat=reshape(Ymat,T,n,N);

% loop over units
for ii=1:N
% estimate the residuals for this unit
EPS(:,:,ii)=Ymat(:,:,ii)-Ytilde(:,:,ii);

% Compute then the sum of squared residuals
% compute first the RSS matrix, defined in (1.9.5)
RSS(:,:,ii)=EPS(:,:,ii)'*EPS(:,:,ii);
% retain only the diagonal elements to get the vector of RSSi values
rss(:,:,ii)=diag(RSS(:,:,ii));

% Go on calculating R2
% generate Mbar
Mbar=eye(T)-ones(T,T)/T;
% then compute the TSS matrix, defined in (1.9.8)
TSS(:,:,ii)=Ymat(:,:,ii)'*Mbar*Ymat(:,:,ii);
% generate the R2 matrix in (1.9.9)
R2(:,:,ii)=eye(n)-RSS(:,:,ii)./TSS(:,:,ii);
% retain only the diagonal elements to get the vector of R2 values
r2(:,:,ii)=diag(R2(:,:,ii));

% then calculate the adjusted R2, using (1.9.11)
R2bar(:,:,ii)=eye(n)-((T-1)/(T-k))*(eye(n)-R2(:,:,ii));
% retain only the diagonal elements to get the vector of R2bar values
r2bar(:,:,ii)=diag(R2bar(:,:,ii));

end

% estimate the forecast evaluation criteria
% first note that forecast evaluation can only be conducted if it is activated, and if there is some observable data after the beginning of the forecast
if Feval==1 && Fcomp==1

% generate the elements required for the evaluation
forecast_record=reshape(forecast_record,N*n,1);
forecast_estimates=reshape(forecast_estimates,N*n,1);
data_endo_c=reshape(data_endo_c,Fcperiods,N*n);

% compute forecast evaluation
[RMSE,MAE,MAPE,Ustat,CRPS_estimates]=bear.panel6feval(N,n,forecast_record,forecast_estimates,Fcperiods,data_endo_c);

end

% start displaying and saving the general results

% preliminary task: create and open the txt file used to save the results

filelocation=fullfile(pref.results_path, [pref.results_sub '.txt']);
fid=fopen(filelocation,'wt');

% print toolbox header

fprintf('%s\n','');
fprintf(fid,'%s\n','');

% print the list of contributors
bear.printcontributors(fid);

% print then estimation results

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

toolboxinfo='BEAR toolbox estimates';
fprintf('%s\n',toolboxinfo);
fprintf(fid,'%s\n',toolboxinfo);

time=clock;
datestring=datestr(time);
dateinfo=['Date: ' datestring(1,1:11) '   Time: ' datestring(1,13:17)];
fprintf('%s\n',dateinfo);
fprintf(fid,'%s\n',dateinfo);

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

VARtypeinfo='Panel VAR: structural factor (dynamic)';
fprintf('%s\n',VARtypeinfo);
fprintf(fid,'%s\n',VARtypeinfo);

if IRFt==1
SVARinfo='structural decomposition: none';
fprintf('%s\n',SVARinfo);
fprintf(fid,'%s\n',SVARinfo);
elseif IRFt==2
SVARinfo='structural decomposition: choleski factorisation';
fprintf('%s\n',SVARinfo);
fprintf(fid,'%s\n',SVARinfo);
elseif IRFt==3
SVARinfo='structural decomposition: triangular factorisation';
fprintf('%s\n',SVARinfo);
fprintf(fid,'%s\n',SVARinfo);
end

temp='units: ';
for ii=1:N
temp=[temp ' ' Units{ii,1} ' '];
end
unitinfo=temp;
fprintf('%s\n',unitinfo);
fprintf(fid,'%s\n',unitinfo);

temp='endogenous variables: ';
for ii=1:n
temp=[temp ' ' endo{ii,1} ' '];
end
endoinfo=temp;
fprintf('%s\n',endoinfo);
fprintf(fid,'%s\n',endoinfo);

temp='exogenous variables: ';
if const==0 && m==0
temp=[temp ' none'];
elseif const==1 && m==1
temp=[temp ' constant '];
elseif const==0 && m>0
   for ii=1:m-1
   temp=[temp ' ' exo{ii,1} ' '];
   end
elseif const==1 && m>1
temp=[temp ' constant '];
   for ii=1:m-1
   temp=[temp ' ' exo{ii,1} ' '];
   end
end
exoinfo=temp;
fprintf('%s\n',exoinfo);
fprintf(fid,'%s\n',exoinfo);

sampledateinfo=['estimation sample: ' startdate '-' enddate];
fprintf('%s\n',sampledateinfo);
fprintf(fid,'%s\n',sampledateinfo);

samplelengthinfo=['sample size (omitting initial conditions): ' num2str(T)];
fprintf('%s\n',samplelengthinfo);
fprintf(fid,'%s\n',samplelengthinfo);

laginfo=['number of lags included in regression: ' num2str(p)];
fprintf('%s\n',laginfo);
fprintf(fid,'%s\n',laginfo);

hyperparam1='hyperparameters:';
fprintf('%s\n',hyperparam1);
fprintf(fid,'%s\n',hyperparam1);

hyperparam2=['IG shape on residual variance (alpha0):         ' num2str(alpha0)];
fprintf('%s\n',hyperparam2);
fprintf(fid,'%s\n',hyperparam2);

hyperparam3=['IG scale on residual variance (delta0):         ' num2str(delta0)];
fprintf('%s\n',hyperparam3);
fprintf(fid,'%s\n',hyperparam3);

hyperparam4=['AR coefficient on residual variance (gamma):    ' num2str(gama)];
fprintf('%s\n',hyperparam4);
fprintf(fid,'%s\n',hyperparam4);

hyperparam5=['IG shape on factor variance (a0):               ' num2str(a0)];
fprintf('%s\n',hyperparam5);
fprintf(fid,'%s\n',hyperparam5);

hyperparam6=['IG scale on factor variance (b0):               ' num2str(b0)];
fprintf('%s\n',hyperparam6);
fprintf(fid,'%s\n',hyperparam6);

hyperparam7=['AR coefficient on factors (rho):                ' num2str(rho)];
fprintf('%s\n',hyperparam7);
fprintf(fid,'%s\n',hyperparam7);

hyperparam8=['variance of Metropolis draw (psi):              ' num2str(psi)];
fprintf('%s\n',hyperparam8);
fprintf(fid,'%s\n',hyperparam8);

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% Metropolis Hastings summary

MHinfo1=['The acceptance rate of the Metropolis-Hastings step during the Gibbs sampler was ' num2str(acceptrate,'%.2f') '%.'];
fprintf('%s\n',MHinfo1);
fprintf(fid,'%s\n',MHinfo1);

MHinfo2='A good acceptance rate is typically comprised between 20% and 30%.';
fprintf('%s\n',MHinfo2);
fprintf(fid,'%s\n',MHinfo2);

if acceptrate<=20

MHinfo3='Your acceptance rate is therefore low, which implies that some values may be repeated';
fprintf('%s\n',MHinfo3);
fprintf(fid,'%s\n',MHinfo3);
MHinfo4='in the estimation of the posterior distribution for Z.';
fprintf('%s\n',MHinfo4);
fprintf(fid,'%s\n',MHinfo4);
MHinfo5='Your may consider decreasing the value of psi in order to remedy to this situation.';
fprintf('%s\n',MHinfo5);
fprintf(fid,'%s\n',MHinfo5);

elseif acceptrate>=30

MHinfo3='Your acceptance rate is therefore high, which implies that only a small portion of the';
fprintf('%s\n',MHinfo3);
fprintf(fid,'%s\n',MHinfo3);
MHinfo4='support of the posterior distribution of Z may have been covered by the algorithm.';
fprintf('%s\n',MHinfo4);
fprintf(fid,'%s\n',MHinfo4);
MHinfo5='Your may consider increasing the value of psi in order to remedy to this situation.';
fprintf('%s\n',MHinfo5);
fprintf(fid,'%s\n',MHinfo5);

elseif (acceptrate>20 && acceptrate<30)

MHinfo3='This is satisfied in this case. Hence, the value of psi that has been selected seems';
fprintf('%s\n',MHinfo3);
fprintf(fid,'%s\n',MHinfo3);
MHinfo4='suitable for an efficient estimation of the posterior distribution of Z.';
fprintf('%s\n',MHinfo4);
fprintf(fid,'%s\n',MHinfo4);

end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display factor estimates

factorinfo=['Structural factors (final sample period):'];
fprintf('%s\n',factorinfo);
fprintf(fid,'%s\n',factorinfo);

fprintf('%s\n','');
fprintf(fid,'%s\n','');

factorheader=fprintf('%35s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');
factorheader=fprintf(fid,'%35s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');

% common component

fprintf('%s\n','theta1 (common component)');
fprintf(fid,'%s\n','theta1 (common component)');
values=[theta_median(1,1,end) theta_std(1,1,end) theta_lbound(1,1,end) theta_ubound(1,1,end)];
fprintf('%-35s %15.3f %15.3f %15.3f %15.3f\n','common component',values);
fprintf(fid,'%-35s %15.3f %15.3f %15.3f %15.3f\n','common component',values);
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% unit component

fprintf('%s\n','theta2 (unit-specific component)');
fprintf(fid,'%s\n','theta2 (unit component)');
for ii=1:d2
values=[theta_median(d1+ii,1,end) theta_std(d1+ii,1,end) theta_lbound(d1+ii,1,end) theta_ubound(d1+ii,1,end)];
fprintf('%-35s %15.3f %15.3f %15.3f %15.3f\n',['unit ' int2str(ii) ' component'],values);
fprintf(fid,'%-35s %15.3f %15.3f %15.3f %15.3f\n',['unit ' int2str(ii) ' component'],values);
end
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% variable component

fprintf('%s\n','theta3 (variable-specific component)');
fprintf(fid,'%s\n','theta3 (endogenous variable component)');
for ii=1:d3
values=[theta_median(d1+d2+ii,1,end) theta_std(d1+d2+ii,1,end) theta_lbound(d1+d2+ii,1,end) theta_ubound(d1+d2+ii,1,end)];
fprintf('%-35s %15.3f %15.3f %15.3f %15.3f\n',['variable ' int2str(ii) ' component'],values);
fprintf(fid,'%-35s %15.3f %15.3f %15.3f %15.3f\n',['variable ' int2str(ii) ' component'],values);
end
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% lag component (if applicable)

if d4~=0
fprintf('%s\n','theta4 (lag-specific component)');
fprintf(fid,'%s\n','theta4 (lag component)');
for ii=1:d4
values=[theta_median(d1+d2+d3+ii,1,end) theta_std(d1+d2+d3+ii,1,end) theta_lbound(d1+d2+d3+ii,1,end) theta_ubound(d1+d2+d3+ii,1,end)];
fprintf('%-35s %15.3f %15.3f %15.3f %15.3f\n',['lag ' int2str(ii) ' component'],values);
fprintf(fid,'%-35s %15.3f %15.3f %15.3f %15.3f\n',['lag ' int2str(ii) ' component'],values);
end
fprintf('%s\n','');
fprintf(fid,'%s\n','');
end

% exogenous component (if applicable)

if d5~=0
fprintf('%s\n','theta5 (exogenous variable component)');
fprintf(fid,'%s\n','theta5 (exogenous component)');
% initiate equation count
eqcount=0;
% initiate exogenous count
exocount=0;
for ii=1:d5
   if exocount==m
   exocount=0;
   end
   if exocount==0
   eqcount=eqcount+1;
   end
exocount=exocount+1;
values=[theta_median(d1+d2+d3+d4+ii,1,end) theta_std(d1+d2+d3+d4+ii,1,end) theta_lbound(d1+d2+d3+d4+ii,1,end) theta_ubound(d1+d2+d3+d4+ii,1,end)];
fprintf('%-35s %15.3f %15.3f %15.3f %15.3f\n',['equation ' int2str(eqcount) ', exogenous ' int2str(exocount) ' component'],values);
fprintf(fid,'%-35s %15.3f %15.3f %15.3f %15.3f\n',['equation ' int2str(eqcount) ', exogenous ' int2str(exocount) ' component'],values);
end
fprintf('%s\n','');
fprintf(fid,'%s\n','');
end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display VAR stability results
eigmodulus=reshape(eigmodulus,p,N*n);
stabilityinfo1=['Roots (modulus) of the characteristic polynomial for the final sample period:'];
fprintf('%s\n',stabilityinfo1);
fprintf(fid,'%s\n',stabilityinfo1);
for jj=1:p
temp=num2str(eigmodulus(jj,1),'%.3f');
   for kk=2:N*n
   temp=[temp,'  ',num2str(eigmodulus(jj,kk),'%.3f')];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end
if stationary==1;
stabilityinfo2=['No root lies outside the unit circle.'];
stabilityinfo3=['The estimated VAR model satisfies the stability condition'];
fprintf('%s\n',stabilityinfo2);
fprintf(fid,'%s\n',stabilityinfo2);
fprintf('%s\n',stabilityinfo3);
fprintf(fid,'%s\n',stabilityinfo3);
else
stabilityinfo2=['Warning: at leat one root lies on or outside the unit circle.'];
stabilityinfo3=['The estimated VAR model will not be stable'];
end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display posterior for sigma
% reshape sigma
sigma_median=reshape(sigma_median(:,:,end),N*n,N*n);
% start displaying
sigmainfo=['sigma (residual covariance matrix): posterior estimates for the final sample period'];
fprintf('%s\n',sigmainfo);
fprintf(fid,'%s\n',sigmainfo);
% calculate the (integer) length of the largest number in sigma, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(sigma_median))))));
% add a separator, a potential minus sign, and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:N*n
temp=[];
   for jj=1:N*n
   % convert matrix entry into string
   number=num2str(sigma_median(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end

fprintf('%s\n','');
fprintf(fid,'%s\n','');

% then display the results for D and gamma, if a structural decomposition was selected

if IRF==1 && IRFt~=1
svarinfo1=['D (structural decomposition matrix): posterior estimates'];
fprintf('%s\n',svarinfo1);
fprintf(fid,'%s\n',svarinfo1);

% recover D
D=reshape(D_estimates(:,:,end),N*n,N*n);
% calculate the (integer) length of the largest number in D, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(D))))));
% add a separator, a potential minus sign and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:N*n
temp=[];
   for jj=1:N*n
   % convert matrix entry into string
   number=num2str(D(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end

fprintf('%s\n','');
fprintf(fid,'%s\n','');

svarinfo2=['gamma (structural disturbances covariance matrix): posterior estimates'];
fprintf('%s\n',svarinfo2);
fprintf(fid,'%s\n',svarinfo2);

% recover gamma
gamma=reshape(gamma_estimates(:,:,end),N*n,N*n);
% calculate the (integer) length of the largest number in D, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(D))))));
% add a separator, a potential minus sign and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:N*n
temp=[];
   for jj=1:N*n
   % convert matrix entry into string
   number=num2str(gamma(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end

end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% finally display in-sample and forecast evaluation measures
% if forecast evaluation is activated and possible, display the results
if Feval==1 && Fcomp==1
bear.panel6fprint(Units,N,n,endo,rss,r2,r2bar,RMSE,MAE,MAPE,Ustat,CRPS_estimates,stringdates3,Fstartdate,Fcenddate,Fcperiods,fid);

% if forecast evaluation is activated but not possible, return a message to signal it
elseif Feval==1 && Fcomp==0

finfo1=['Forecast evaluation cannot be conducted.'];
fprintf('%s\n',finfo1);
fprintf(fid,'%s\n',finfo1);
finfo2=['Forecasts start in ' Fstartdate ', while observable data is available only until ' names{end,1} '.'];
fprintf('%s\n',finfo2);
fprintf(fid,'%s\n',finfo2);
finfo3=['To obtain forecast evaluation, the forecast start date must be anterior to the end of the data set.'];
fprintf('%s\n',finfo3);
fprintf(fid,'%s\n',finfo3);

% if forecast evaluation is not activated altogether, do not do anything

end

author: european-central-bank

repo: BEAR-toolbox


text:
%% Importing Crystal Orientations
%
%%
% In order to import orientation data from a text file we first need to
% defined the corresponding <CrystalSymmetries.html crystal symmetry>, e.g.
% by a cif file

cs = crystalSymmetry.load('quartz.cif')

%%
% In the second step we may use the command <orientation.load.html
% |orientation.load|> to load the data. This function requires that one
% specifies the meaning of the column of the import file by the option
% |columnNames|.

fname = fullfile(mtexDataPath,'orientation','Tongue_Quartzite_Bunge_Euler');
ori = orientation.load(fname,'columnNames',{'phi1','Phi','phi2'},cs)

%%
% This creates a variable of type @orientation which can be used for
% further analysis, e.g. to plot pole figures

plotPDF(ori,Miller({0,0,0,1},{1,0,-1,0},cs))

author: mtex-toolbox

repo: mtex

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


text:
function [h mu ul ll] = circ_mtest(alpha, dir, xi, w, d)
%
% [pval, z] = circ_mtest(alpha, dir, w, d)
%   One-Sample test for the mean angle.
%   H0: the population has mean dir.
%   HA: the population has not mean dir.
%
%   Note: This is the equvivalent to a one-sample t-test with specified
%         mean direction.
%
%   Input:
%     alpha	sample of angles in radians
%     dir   assumed mean direction
%     [xi   alpha level of the test]
%     [w		number of incidences in case of binned angle data]
%     [d    spacing of bin centers for binned data, if supplied
%           correction factor is used to correct for bias in
%           estimation of r, in radians (!)]
%
%   Output:
%     h     0 if H0 can not be rejected, 1 otherwise
%     mu    mean
%     ul    upper (1-xi) confidence level
%     ll    lower (1-xi) confidence level
%
% PHB 7/6/2008
%
% References:
%   Biostatistical Analysis, J. H. Zar
%
% Circular Statistics Toolbox for Matlab

% By Philipp Berens, 2009
% berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html

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

if nargin<3
  xi = 0.05;
end

if nargin<4
  % if no specific weighting has been specified
  % assume no binning has taken place
	w = ones(size(alpha));
else
  if size(w,2) > size(w,1)
    w = w';
  end
  if length(alpha)~=length(w)
    error('Input dimensions do not match.')
  end
end

if nargin<5
  % per default do not apply correct for binned data
  d = 0;
end

% compute ingredients
mu = circ_mean(alpha,w);
t = circ_confmean(alpha,xi,w,d);
ul = mu + t;
ll = mu - t;

% compute test via confidence limits (example 27.3)
h = abs(circ_dist2(dir,mu)) > t;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [gen, ave1, ave2] = fusion_strategy(features_a, features_b, source_a, source_b, unit)

[m,n] = size(features_a);
[m1,n1] = size(source_a);
ave_temp1 = zeros(m1,n1);
ave_temp2 = zeros(m1,n1);
weight_ave_temp1 = zeros(m1,n1);
weight_ave_temp2 = zeros(m1,n1);

for i=2:m-1
    for j=2:n-1
        A1 =sum(sum(features_a(i-1:i+1,j-1:j+1)))/9;
        A2 =sum(sum(features_b(i-1:i+1,j-1:j+1)))/9;

        % weight average
        weight_ave_temp1(((i-2)*unit+1):((i-1)*unit),((j-2)*unit+1):((j-1)*unit)) = A1/(A1+A2);
        weight_ave_temp2(((i-2)*unit+1):((i-1)*unit),((j-2)*unit+1):((j-1)*unit)) = A2/(A1+A2);
        ave_temp1(((i-2)*unit+1):((i-1)*unit),((j-2)*unit+1):((j-1)*unit)) = A1;
        ave_temp2(((i-2)*unit+1):((i-1)*unit),((j-2)*unit+1):((j-1)*unit)) = A2;
%         % choose max
%         if A1>A2
%             gen_temp(((i-2)*unit+1):((i-1)*unit),((j-2)*unit+1):((j-1)*unit)) = ones(unit,unit);
%             temp_mask(i,j) = 1;
%         end
    end
end
% figure;imshow(temp_mask);
weight_ave_temp1 = weight_ave_temp1(1:m1,1:n1);
weight_ave_temp2 = weight_ave_temp2(1:m1,1:n1);
% figure;imshow(weight_ave_temp1);
% figure;imshow(weight_ave_temp2);

gen = source_a.*weight_ave_temp1 + source_b.*weight_ave_temp2;
% figure;imshow(gen);

ave1 = ave_temp1;
ave2 = ave_temp2;
end

author: thfylsty


text:
% SUMMARY:  This is a Discrete Hidden Markov code
%           This code is inspired by Murphy's PMTK3 toolbox. Using EM
%           algorithm. Details are from PRML
% AUTHOR:   QIUQIANG KONG, Queen Mary University of London
% Created:  17-09-2015
% Modified: 17-11-2015
% Ref      Chap 13. <Pattern Analysis and Machine Learning>
% -----------------------------------------------------------
% input
%   Data      cell of data
%   state_num state num
%   mix_num   multinominal num
% varargin input:
%   p_start0  p(z1), size: Q*1
%   A         p(zn|zn-1), transform matrix, size: Q*Q
%   phi0:     emission probability para
%       B       p(xn|zn), emission matrix, size: p*Q
%   iter_num  how many time the EM should run (default: 100)
%   converge  (default: 1+1e-4)
% output
%   p_start  p(z1), dim 1: Q
%   A        p(zn|zn-1), transform matrix, size: Q*Q
%   mu       p(xn|zn), emission matrix, size: p*Q
% ===========================================================
function [p_start, A, phi, loglik] = Dhmm(Data, state_num, mix_num, varargin)
for i1 = 1:2:length(varargin)
    switch varargin{i1}
        case 'p_start0'
            p_start = varargin{i1+1};
        case 'A0'
            A = varargin{i1+1};
        case 'phi0'
            phi = varargin{i1+1};
        case 'iter_num'
            iter_num = varargin{i1+1};
        case 'converge'
            converge = varargin{i1+1};
    end
end
Q = state_num;
M = mix_num;
if (~exist('p_start'))
    tmp = rand(1,Q);
    p_start = tmp / sum(tmp);
end
if (~exist('A'))
    tmp = rand(Q,Q);
    A = bsxfun(@rdivide, tmp, sum(tmp,2));
end
if (~exist('phi'))
    phi.B = ones(M,Q) / M;
end
if (~exist('iter_num'))
    iter_num = 100;
end
if (~exist('converge'))
    converge = 1 + 1e-4;
end

obj_num = length(Data);     % sequences num
[M,Q] = size(phi.B);        % multimominal num, state num
pre_ll = -inf;

for k = 1:iter_num
    % E STEP
    for r = 1:obj_num
        logp_xn_given_zn = Discrete_logp_xn_given_zn(Data{r}, phi);
        [LogGamma{r}, LogKsi{r}, Loglik{r}] = LogForwardBackward(logp_xn_given_zn, p_start, A);
    end

    % convert loggamma to gamma, logksi to ksi, substract the max
    [Gamma, Ksi] = UniformLogGammaKsi(LogGamma, LogKsi);

    % M STEP common
    [p_start, A] = M_step_common(Gamma, Ksi);

    % M STEP for Multinominal distribution
    B_numer = zeros(M,Q);
    B_denom = zeros(1,Q);
    for r = 1:obj_num
        xr = Data{r};
        Nr = length(xr);
        Xr = zeros(Nr,M);
        Xr(sub2ind([Nr,M], 1:Nr, xr')) = 1;
        B_numer = B_numer + Xr' * Gamma{r};
        B_denom = B_denom + sum(Gamma{r},1);
    end
    phi.B = bsxfun(@rdivide, B_numer, B_denom);

    % calculate loglik
    loglik = 0;
    for r = 1:obj_num
        loglik = loglik + Loglik{r};
    end
    if (loglik-pre_ll<log(converge)) break;
    else pre_ll = loglik; end
end

author: qiuqiangkong

repo: matlab-hmm


text:
function varargout = mul_const(x, h, der_y)

% x is [m1, m2, p, b]
% h is [m1, m2]
% der_y is same size as x

if nargin < 3
    der_y = [];
end

if isempty(der_y)
    y = bsxfun(@times, x, h);
    varargout = {y};
else
    der_x = bsxfun(@times, der_y, h);
    der_h = sum(sum(der_y .* x, 3), 4);
    varargout = {der_x, der_h};
end

end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
function [T]=txtfile2cell(fileName)
% function [T]=txtfile2cell(fileName)
% ------------------------------------------------------------------------
% This function read the text file specified by fileName (path and name)
% whereby each line is read into a seperate entry in the output cell array
% T.
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2015/102/10
%------------------------------------------------------------------------

fid=fopen(fileName);
T=textscan(fid,'%s','delimiter', '\n','Whitespace','');
T=T{1,1};
fclose(fid);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [net, info] = cnn_train(net, imdb, getBatch, varargin)
%CNN_TRAIN  An example implementation of SGD for training CNNs
%    CNN_TRAIN() is an example learner implementing stochastic
%    gradient descent with momentum to train a CNN. It can be used
%    with different datasets and tasks by providing a suitable
%    getBatch function.
%
%    The function automatically restarts after each training epoch by
%    checkpointing.
%
%    The function supports training on CPU or on one or more GPUs
%    (specify the list of GPU IDs in the `gpus` option). Multi-GPU
%    support is relatively primitive but sufficient to obtain a
%    noticable speedup.

% Copyright (C) 2014-15 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.numEpochs = 300 ;
opts.gpus = [] ; % which GPU devices to use (none, one, or more)
opts.learningRate = 0.001 ;
opts.continue = false ;
opts.expDir = fullfile('data','exp') ;
opts.conserveMemory = false ;
opts.backPropDepth = +inf ;
opts.sync = false ;
opts.prefetch = false ;
opts.cudnn = true ;
opts.weightDecay = 0.0005 ;
opts.momentum = 0.9 ;
opts.errorFunction = 'multiclass' ;
opts.errorLabels = {} ;
opts.plotDiagnostics = false ;
opts.memoryMapFile = fullfile(tempdir, 'matconvnet.bin') ;
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isnan(opts.train), opts.train = [] ; end

% -------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------

evaluateMode = isempty(opts.train) ;

if ~evaluateMode
  for i=1:numel(net.layers)
    if isfield(net.layers{i}, 'weights')
      J = numel(net.layers{i}.weights) ;
      for j=1:J
        net.layers{i}.momentum{j} = zeros(size(net.layers{i}.weights{j}), 'single') ;
      end
      if ~isfield(net.layers{i}, 'learningRate')
        net.layers{i}.learningRate = ones(1, J, 'single') ;
      end
      if ~isfield(net.layers{i}, 'weightDecay')
        net.layers{i}.weightDecay = ones(1, J, 'single') ;
      end
    end
    % Legacy code: will be removed
    if isfield(net.layers{i}, 'filters')
      net.layers{i}.momentum{1} = zeros(size(net.layers{i}.filters), 'single') ;
      net.layers{i}.momentum{2} = zeros(size(net.layers{i}.biases), 'single') ;
      if ~isfield(net.layers{i}, 'learningRate')
        net.layers{i}.learningRate = ones(1, 2, 'single') ;
      end
      if ~isfield(net.layers{i}, 'weightDecay')
        net.layers{i}.weightDecay = single([1 0]) ;
      end
    end
  end
end

% setup GPUs
numGpus = numel(opts.gpus) ;
if numGpus > 1
  if isempty(gcp('nocreate')),
    parpool('local',numGpus) ;
    spmd, gpuDevice(opts.gpus(labindex)), end
  end
elseif numGpus == 1
  gpuDevice(opts.gpus)
end
if exist(opts.memoryMapFile), delete(opts.memoryMapFile) ; end

% setup error calculation function
if isstr(opts.errorFunction)
  switch opts.errorFunction
    case 'none'
      opts.errorFunction = @error_none ;
    case 'multiclass'
      opts.errorFunction = @error_multiclass ;
      if isempty(opts.errorLabels), opts.errorLabels = {'top1e', 'top5e'} ; end
    case 'binary'
      opts.errorFunction = @error_binary ;
      if isempty(opts.errorLabels), opts.errorLabels = {'bine'} ; end
    otherwise
      error('Uknown error function ''%s''', opts.errorFunction) ;
  end
end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('resuming by loading epoch %d\n', start) ;
  load(modelPath(start), 'net', 'info') ;
end

for epoch=start+1:opts.numEpochs

  % train one epoch and validate
  learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  train = opts.train(randperm(numel(opts.train))) ; % shuffle
  val = opts.val ;
  if numGpus <= 1
    [net,stats.train] = process_epoch(opts, getBatch, epoch, train, learningRate, imdb, net) ;
    [~,stats.val] = process_epoch(opts, getBatch, epoch, val, 0, imdb, net) ;
  else
    spmd(numGpus)
      [net_, stats_train_] = process_epoch(opts, getBatch, epoch, train, learningRate, imdb, net) ;
      [~, stats_val_] = process_epoch(opts, getBatch, epoch, val, 0, imdb, net_) ;
    end
    net = net_{1} ;
    stats.train = sum([stats_train_{:}],2) ;
    stats.val = sum([stats_val_{:}],2) ;
  end

  % save
  if evaluateMode, sets = {'val'} ; else sets = {'train', 'val'} ; end
  for f = sets
    f = char(f) ;
    n = numel(eval(f)) ;
    info.(f).speed(epoch) = n / stats.(f)(1) * max(1, numGpus) ;
    info.(f).objective(epoch) = stats.(f)(2) / n ;
    info.(f).error(:,epoch) = stats.(f)(3:end) / n ;
  end
  if ~evaluateMode, save(modelPath(epoch), 'net', 'info') ; end

  figure(1) ; clf ;
  hasError = isa(opts.errorFunction, 'function_handle') ;
  subplot(1,1+hasError,1) ;
  if ~evaluateMode
    semilogy(1:epoch, info.train.objective, '.-', 'linewidth', 2) ;
    hold on ;
  end
  semilogy(1:epoch, info.val.objective, '.--') ;
  xlabel('training epoch') ; ylabel('energy') ;
  grid on ;
  h=legend(sets) ;
  set(h,'color','none');
  title('objective') ;
  if hasError
    subplot(1,2,2) ; leg = {} ;
    if ~evaluateMode
      plot(1:epoch, info.train.error', '.-', 'linewidth', 2) ;
      hold on ;
      leg = horzcat(leg, strcat('train ', opts.errorLabels)) ;
    end
    plot(1:epoch, info.val.error', '.--') ;
    leg = horzcat(leg, strcat('val ', opts.errorLabels)) ;
    set(legend(leg{:}),'color','none') ;
    grid on ;
    xlabel('training epoch') ; ylabel('error') ;
    title('error') ;
  end
  drawnow ;
  print(1, modelFigPath, '-dpdf') ;
end

% -------------------------------------------------------------------------
function err = error_multiclass(opts, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
[~,predictions] = sort(predictions, 3, 'descend') ;

% be resilient to badly formatted labels
if numel(labels) == size(predictions, 4)
  labels = reshape(labels,1,1,1,[]) ;
end

% skip null labels
mass = single(labels(:,:,1,:) > 0) ;
if size(labels,3) == 2
  % if there is a second channel in labels, used it as weights
  mass = mass .* labels(:,:,2,:) ;
  labels(:,:,2,:) = [] ;
end

error = ~bsxfun(@eq, predictions, labels) ;
err(1,1) = sum(sum(sum(mass .* error(:,:,1,:)))) ;
err(2,1) = sum(sum(sum(mass .* min(error(:,:,1:5,:),[],3)))) ;

% -------------------------------------------------------------------------
function err = error_binaryclass(opts, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
error = bsxfun(@times, predictions, labels) < 0 ;
err = sum(error(:)) ;

% -------------------------------------------------------------------------
function err = error_none(opts, labels, res)
% -------------------------------------------------------------------------
err = zeros(0,1) ;

% -------------------------------------------------------------------------
function  [net_cpu,stats,prof] = process_epoch(opts, getBatch, epoch, subset, learningRate, imdb, net_cpu)
% -------------------------------------------------------------------------

% move CNN to GPU as needed
numGpus = numel(opts.gpus) ;
if numGpus >= 1
  net = vl_simplenn_move(net_cpu, 'gpu') ;
else
  net = net_cpu ;
  net_cpu = [] ;
end

% validation mode if learning rate is zero
training = learningRate > 0 ;
if training, mode = 'training' ; else, mode = 'validation' ; end
if nargout > 2, mpiprofile on ; end

numGpus = numel(opts.gpus) ;
if numGpus >= 1
  one = gpuArray(single(1)) ;
else
  one = single(1) ;
end
res = [] ;
mmap = [] ;
stats = [] ;
start = tic ;

for t=1:opts.batchSize:numel(subset)
  fprintf('%s: epoch %02d: batch %3d/%3d: ', mode, epoch, ...
          fix(t/opts.batchSize)+1, ceil(numel(subset)/opts.batchSize)) ;
  batchSize = min(opts.batchSize, numel(subset) - t + 1) ;
  numDone = 0 ;
  error = [] ;
  for s=1:opts.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+opts.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : opts.numSubBatches * numlabs : batchEnd) ;
    [im, labels] = getBatch(imdb, batch) ;

    if opts.prefetch
      if s==opts.numSubBatches
        batchStart = t + (labindex-1) + opts.batchSize ;
        batchEnd = min(t+2*opts.batchSize-1, numel(subset)) ;
      else
        batchStart = batchStart + numlabs ;
      end
      nextBatch = subset(batchStart : opts.numSubBatches * numlabs : batchEnd) ;
      getBatch(imdb, nextBatch) ;
    end

    if numGpus >= 1
      im = gpuArray(im) ;
    end

    % evaluate CNN
    net.layers{end}.class = labels ;
    if training, dzdy = one; else, dzdy = [] ; end
    res = vl_simplenn(net, im, dzdy, res, ...
                      'accumulate', s ~= 1, ...
                      'disableDropout', ~training, ...
                      'conserveMemory', opts.conserveMemory, ...
                      'backPropDepth', opts.backPropDepth, ...
                      'sync', opts.sync, ...
                      'cudnn', opts.cudnn) ;

    % accumulate training errors
    error = sum([error, [...
      sum(double(gather(res(end).x))) ;
      reshape(opts.errorFunction(opts, labels, res),[],1) ; ]],2) ;
    numDone = numDone + numel(batch) ;
  end

  % gather and accumulate gradients across labs
  if training
    if numGpus <= 1
      [net,res] = accumulate_gradients(opts, learningRate, batchSize, net, res) ;
    else
      if isempty(mmap)
        mmap = map_gradients(opts.memoryMapFile, net, res, numGpus) ;
      end
      write_gradients(mmap, net, res) ;
      labBarrier() ;
      [net,res] = accumulate_gradients(opts, learningRate, batchSize, net, res, mmap) ;
    end
  end

  % print learning statistics

  time = toc(start) ;
  stats = sum([stats,[0 ; error]],2); % works even when stats=[]
  stats(1) = time ;
  n = (t + batchSize - 1) / max(1,numlabs) ;
  speed = n/time ;
  fprintf('%.1f Hz%s\n', speed) ;

  fprintf(' obj:%.3g', stats(2)/n) ;
  for i=1:numel(opts.errorLabels)
    fprintf(' %s:%.3g', opts.errorLabels{i}, stats(i+2)/n) ;
  end
  fprintf(' [%d/%d]', numDone, batchSize);
  fprintf('\n') ;

  % debug info
  if opts.plotDiagnostics && numGpus <= 1
    figure(2) ; vl_simplenn_diagnose(net,res) ; drawnow ;
  end
end

if nargout > 2
  prof = mpiprofile('info');
  mpiprofile off ;
end

if numGpus >= 1
  net_cpu = vl_simplenn_move(net, 'cpu') ;
else
  net_cpu = net ;
end

% -------------------------------------------------------------------------
function [net,res] = accumulate_gradients(opts, lr, batchSize, net, res, mmap)
% -------------------------------------------------------------------------
for l=numel(net.layers):-1:1
  for j=1:numel(res(l).dzdw)
    thisDecay = opts.weightDecay * net.layers{l}.weightDecay(j) ;
    thisLR = lr * net.layers{l}.learningRate(j) ;

    % accumualte from multiple labs (GPUs) if needed
    if nargin >= 6
      tag = sprintf('l%d_%d',l,j) ;
      tmp = zeros(size(mmap.Data(labindex).(tag)), 'single') ;
      for g = setdiff(1:numel(mmap.Data), labindex)
        tmp = tmp + mmap.Data(g).(tag) ;
      end
      res(l).dzdw{j} = res(l).dzdw{j} + tmp ;
    end

    if isfield(net.layers{l}, 'weights')
      net.layers{l}.momentum{j} = ...
        opts.momentum * net.layers{l}.momentum{j} ...
        - thisDecay * net.layers{l}.weights{j} ...
        - (1 / batchSize) * res(l).dzdw{j} ;
      net.layers{l}.weights{j} = net.layers{l}.weights{j} + thisLR * net.layers{l}.momentum{j} ;
    else
      % Legacy code: to be removed
      if j == 1
        net.layers{l}.momentum{j} = ...
          opts.momentum * net.layers{l}.momentum{j} ...
          - thisDecay * net.layers{l}.filters ...
          - (1 / batchSize) * res(l).dzdw{j} ;
        net.layers{l}.filters = net.layers{l}.filters + thisLR * net.layers{l}.momentum{j} ;
      else
        net.layers{l}.momentum{j} = ...
          opts.momentum * net.layers{l}.momentum{j} ...
          - thisDecay * net.layers{l}.biases ...
          - (1 / batchSize) * res(l).dzdw{j} ;
        net.layers{l}.biases = net.layers{l}.biases + thisLR * net.layers{l}.momentum{j} ;
      end
    end
  end
end

% -------------------------------------------------------------------------
function mmap = map_gradients(fname, net, res, numGpus)
% -------------------------------------------------------------------------
format = {} ;
for i=1:numel(net.layers)
  for j=1:numel(res(i).dzdw)
    format(end+1,1:3) = {'single', size(res(i).dzdw{j}), sprintf('l%d_%d',i,j)} ;
  end
end
format(end+1,1:3) = {'double', [3 1], 'errors'} ;
if ~exist(fname) && (labindex == 1)
  f = fopen(fname,'wb') ;
  for g=1:numGpus
    for i=1:size(format,1)
      fwrite(f,zeros(format{i,2},format{i,1}),format{i,1}) ;
    end
  end
  fclose(f) ;
end
labBarrier() ;
mmap = memmapfile(fname, 'Format', format, 'Repeat', numGpus, 'Writable', true) ;

% -------------------------------------------------------------------------
function write_gradients(mmap, net, res)
% -------------------------------------------------------------------------
for i=1:numel(net.layers)
  for j=1:numel(res(i).dzdw)
    mmap.Data(labindex).(sprintf('l%d_%d',i,j)) = gather(res(i).dzdw{j}) ;
  end
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

author: huangzehao

repo: caffe-vdsr


text:
function  [h, compUpV, compUp] =  lfmComputeH4VP(gamma1_p, gamma1_m, sigma2, t1, ...
    preFactor, preExp, mode)

% LFMCOMPUTEH4VP Helper function for computing part of the LFMVXLFM kernel.
% FORMAT
% DESC computes a portion of the LFMVXLFM kernel.
% ARG gamma1 : Gamma value for first system.
% ARG gamma2 : Gamma value for second system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG mode: indicates in which way the vectors t1 and t2 must be transposed
% RETURN h : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio A. Alvarez, 2010
%
% SEEALSO : lfmComputeH4Hat, lfmXlfmKernCompute

% KERN

% This could also be used with str2func changing between 'lfm' and 'lfmvp'

if mode==0
    if nargout > 1
        [compUpV{1}, compUp{1}] = lfmvpComputeUpsilonVector(gamma1_p,sigma2, t1, 0);
        [compUpV{2}, compUp{2}] = lfmvpComputeUpsilonVector(gamma1_m,sigma2, t1, 0);
        h =  compUpV{1}*( preExp(:,1)/preFactor(1) - preExp(:,2)/preFactor(2)).' ...
            + compUpV{2}*( preExp(:,2)/preFactor(3) - preExp(:,1)/preFactor(4)).';
    else
        h =  lfmvpComputeUpsilonVector(gamma1_p,sigma2, t1, 0)*( preExp(:,1)/preFactor(1) - preExp(:,2)/preFactor(2)).' ...
            + lfmvpComputeUpsilonVector(gamma1_m,sigma2, t1, 0)*( preExp(:,2)/preFactor(3) - preExp(:,1)/preFactor(4)).';
    end
else
    if nargout > 1
        compUp{1} = lfmComputeUpsilonVector(gamma1_p,sigma2, t1);
        compUp{2} = lfmComputeUpsilonVector(gamma1_m,sigma2, t1);
        h =  compUp{1}*(  preExp(:,2)/preFactor(2) - preExp(:,1)/preFactor(1)).' ...
            + compUp{2}*( preExp(:,1)/preFactor(4) - preExp(:,2)/preFactor(3)).';
        compUpV = compUp;
    else
        h =  lfmComputeUpsilonVector(gamma1_p,sigma2, t1)*(  preExp(:,2)/preFactor(2) - preExp(:,1)/preFactor(1)).' ...
            + lfmComputeUpsilonVector(gamma1_m,sigma2, t1)*( preExp(:,1)/preFactor(4) - preExp(:,2)/preFactor(3)).';
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [metabolite_structure,IDsAdded] = getSeed2Kegg(metabolite_structure)
% This function parses the file: ftp://ftp.kbase.us/assets/KBase_Reference_Data/Biochemistry/compounds.xls
% first column contains seed ID, 5th col contains Kegg ID. This file is
% provided in /data/ as 'compounds.xlsx'
%
% INPUT
% metabolite_structure  metabolite structure
%
% OUTPUT
% metabolite_structure  Updated metabolite structure
%
%
% Ines Thiele, 2020-2021

annotationSource = 'Seed mapping to Kegg';
annotationType = 'automatic';

Mets = fieldnames(metabolite_structure);
for i = 1 : length(Mets)
    SeedId{i,1} = Mets{i};
    SeedId{i,2} = metabolite_structure.(Mets{i}).seed;
end

a = 1;
IDsAdded = '';

[NUM,TXT,RAW]=xlsread('compounds.xlsx');
for i = 1 : size(SeedId,1)
    % no keggId exists
    if isempty(metabolite_structure.(SeedId{i,1}).keggId) || ~isempty(find(isnan(metabolite_structure.(SeedId{i,1}).keggId),1))
        match = strmatch(SeedId{i,2},RAW(:,1),'exact');
        if ~isempty(match)
            if contains(RAW(match,5),'|') % take only first kegg Id
                [RAW{match,5},rem]=strtok(RAW{match,5},'|');
            end
            if ~isempty(RAW{match,5})
                metabolite_structure.(SeedId{i,1}).keggId = RAW{match,5};
                metabolite_structure.(SeedId{i,1}).keggId_source =  [annotationSource,':',annotationType,':',datestr(now)];
                IDsAdded{a,1} = SeedId{i,1};
                IDsAdded{a,2} = 'keggId';
                IDsAdded{a,3} =  metabolite_structure.(SeedId{i,1}).keggId;
                a = a + 1;
            end
        end
    end
end


text:
function bestp = faceGA
%
% A silly face example of a genetic algorithm.
% By Tor Wager, June 2002
%
% calls functions:
% drawFace.m
% breedorgs.m   (from OptimizeDesign10)

gensize = 100;
numgens = 100;
numopt = 17;
change_every = 30;

% ------------------------------------------
% * set up figures
% ------------------------------------------

fgen = figure;set(gcf,'MenuBar','none');set(gcf,'Name','Faces in This Generation')
set(gcf,'Position',[30          70         975        1040])
set(gcf,'NumberTitle','off');set(gcf,'Color','w')

% for movie single figure
for i = 1:30, h(i) = subplot(6,5,i);, mypos(i,:) = get(h(i),'Position');,end
mypos(:,1) = mypos(:,1) - .1;
mypos(:,3) = mypos(:,3) * .7;
mypos(:,1) = mypos(:,1) * .7;
clf
for i = 1:30, nh(i) = axes('position',mypos(i,:));, end
set(gcf,'Position',[30          70        1514        1040])
hgen = axes('position',[.591 .5347 .36 .36]);
hideal = axes('position',[.591 .11 .36 .36]);

mov = avifile('facega.avi','Quality',75,'Compression','none','Fps',10);

%bgen = figure;set(gcf,'MenuBar','none');set(gcf,'Name','Best of This Generation')
%set(gcf,'Position',[1021         719         516         393])
%set(gcf,'NumberTitle','off');set(gcf,'Color','w')
%tgen = figure;set(gcf,'MenuBar','none');set(gcf,'Name','Target Face')
%set(gcf,'Position',[1024         195         516         393])
%set(gcf,'NumberTitle','off');set(gcf,'Color','w')

% ------------------------------------------
% * set up lists
% ------------------------------------------
p = rand(gensize,20);
idealface = [1.0000    1.0000    1.0000    1.0000         0    0.4486    0.5244    0.1715    0.5000    1.0000    0.2000    0.1414    0.4570    1.0000 0.8000    0.2248    0.9089    0.0073    0.5887    0.5421];
idealfit = repmat(idealface,gensize,1);
%figure(tgen); drawFace(idealface)
axes(hideal); cla; drawFace(idealface); title('Face To Match','FontSize',14)

for gen = 1:numgens

    if mod(gen,change_every) == 0
        idealface = rand(1,20);
        idealfit = repmat(idealface,gensize,1);
        axes(hideal); cla; drawFace(idealface)
    end

    for g = 1:gensize

        if g < 31

            axes(nh(g))
            %figure(fgen)
            %subplot(6,5,g)
            drawFace(p(g,:))
            H = gcf; mov = addframe(mov,H);
        end

    end

    % compute fitness and save best
    f = -1 * (sum(((idealfit(:,1:numopt) - p(:,1:numopt)) .^ 2),2));
    bestp = p(f == max(f),:); bestp = bestp(1,:);
    %figure(bgen); clf; drawFace(bestp)
    axes(hgen); cla; drawFace(bestp); title('Best of This Generation','FontSize',14)

    fm = f(f == max(f));
    text(.8,.8,1,num2str(round(fm(1)*1000)),'FontSize',18)

    f = (f - mean(f)) ./ std(f);

    % crossbreed
    p = breedorgs(f',p',2.1)';

    H = gcf; mov = addframe(mov,H);
end

mov = close(mov);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef LvdSensorSet < matlab.mixin.SetGet
    %LvdSensorSet Summary of this class goes here
    %   Detailed explanation goes here

    properties
        sensors(1,:) AbstractSensor

        lvdData LvdData
    end

    methods
        function obj = LvdSensorSet(lvdData)
            obj.lvdData = lvdData;
        end

        function addSensor(obj, sensor)
            obj.sensors(end+1) = sensor;
        end

        function removeSensor(obj, sensor)
            obj.sensors([obj.sensors] == sensor) = [];
        end

        function [listBoxStr, sensors] = getListboxStr(obj)
            listBoxStr = {};

            for(i=1:length(obj.sensors))
                listBoxStr{end+1} = obj.sensors(i).getListboxStr(); %#ok<AGROW>
            end

            sensors = obj.sensors;
        end

        function indPoints = getSensorsForInds(obj, inds)
            indPoints = obj.sensors(inds);
        end

        function inds = getIndsForSensors(obj, indSensors)
            inds = [];
            for(i=1:length(indSensors))
                for(j=1:length(obj.sensors))
                    if(indSensors(i) == obj.sensors(j))
                        inds(end+1) = j; %#ok<AGROW>
                        break;
                    end
                end
            end
        end

        function indSensor = getSensorAtInd(obj, ind)
            if(ind > 0 && ind <= length(obj.sensors))
                indSensor = obj.sensors(ind);
            else
                indSensor = AbstractSensor.empty(1,0);
            end
        end

        function tf = usesGroundObj(obj, groundObj)
            tf = false;
            for(i=1:length(obj.sensors))
                tf = tf || obj.sensors(i).usesGroundObj(groundObj);
            end
        end

        function tf = usesGeometricPoint(obj, point)
            tf = false;
            for(i=1:length(obj.sensors))
                tf = tf || obj.sensors(i).usesGeometricPoint(point);
            end
        end

        function tf = usesGeometricVector(obj, vector)
            tf = false;
            for(i=1:length(obj.sensors))
                tf = tf || obj.sensors(i).usesGeometricVector(vector);
            end
        end

        function tf = usesGeometricCoordSys(obj, coordSys)
            tf = false;
            for(i=1:length(obj.sensors))
                tf = tf || obj.sensors(i).usesGeometricCoordSys(coordSys);
            end
        end

        function tf = usesGeometricRefFrame(obj, refFrame)
            tf = false;
            for(i=1:length(obj.sensors))
                tf = tf || obj.sensors(i).usesGeometricRefFrame(refFrame);
            end
        end

        function tf = usesGeometricAngle(obj, angle)
            tf = false;
            for(i=1:length(obj.sensors))
                tf = tf || obj.sensors(i).usesGeometricAngle(angle);
            end
        end

        function tf = usesGeometricPlane(obj, plane)
            tf = false;
            for(i=1:length(obj.sensors))
                tf = tf || obj.sensors(i).usesGeometricPlane(plane);
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [y,labels] = generateGMM(n,lambda,mu,gamma,verbose)
%
% function y = generateGMM(n,lambda,mu)
% This function generates n observations following a MoG distribution
% given by the class frequencies lambda, their mean mu and their variance gamma.

try, verbose = ~~verbose; catch, verbose = 1; end

p = size(mu,1);
k = 1;
sumLambdak = lambda(1);

if verbose && (p == 2 || p == 3)
    col = repmat('bgrcmyk',1,ceil(1+length(gamma)./7));
    pl = figure;
    hold on
end

labels = zeros(n,length(gamma));

for i = 1:n

    ratio = i./n;

    if ratio>=sumLambdak && ~isequal(i,n)
        k = k+1;
        sumLambdak = sumLambdak+lambda(k);
    end

    labels(i,k) = 1;
    y(:,i) = mu(:,k) + sqrt(gamma(k)).*randn(p,1);

    if verbose && p == 2
        figure(pl),plot(y(1,i),y(2,i),[col(k),'.']);
    elseif verbose && p == 3
        figure(pl),plot3(y(1,i),y(2,i),y(3,i),[col(k),'.']);
    end
%     drawnow;
end

if verbose &&(p==2 || p==3)
    axis square
end

author: MBB-team

repo: VBA-toolbox


text:
function y = pinvNx2(x)

% PINVNX2 computes a pseudo-inverse of the slices of an Nx2xM real-valued matrix.
% Output has dimensionality 2xNxM. This implementation is generally faster
% than calling pinv in a for-loop, once M > 2

siz = [size(x) 1];
xtx = zeros([2,2,siz(3:end)]);
xtx(1,1,:,:) = sum(x(:,1,:,:).^2,1);
xtx(2,2,:,:) = sum(x(:,2,:,:).^2,1);
tmp          = sum(x(:,1,:,:).*x(:,2,:,:),1);
xtx(1,2,:,:) = tmp;
xtx(2,1,:,:) = tmp;
ixtx         = inv2x2(xtx);
y            = mtimes2xN(ixtx, permute(x, [2 1 3:ndims(x)]));

function [d] = inv2x2(x)

% INV2X2 computes inverse of matrix x, where x = 2x2xN, using explicit analytic definition

adjx  = [x(2,2,:,:) -x(1,2,:,:); -x(2,1,:,:) x(1,1,:,:)];
denom = x(1,1,:,:).*x(2,2,:,:) - x(1,2,:,:).*x(2,1,:,:);
d     = adjx./denom([1 1],[1 1],:,:);

function [z] = mtimes2xN(x, y)

% MTIMES2XN computes x*y where x = 2x2xM and y = 2xNxM
% and output dimensionatity is 2xNxM

siz   = size(y);
z     = zeros(siz);

for k = 1:siz(2)
  z(1,k,:,:) = x(1,1,:,:).*y(1,k,:,:) + x(1,2,:,:).*y(2,k,:,:);
  z(2,k,:,:) = x(2,1,:,:).*y(1,k,:,:) + x(2,2,:,:).*y(2,k,:,:);
end

author: spm


text:
function h = ksizeROT(npd,noIQR)
% "Rule of Thumb" estimate (Silverman)
%    Estimate is based on assumptions of Gaussian data and kernel
%    Actually the multivariate version in Scott ('92)
%  Use ksizeROT(X,1) to force use of stddev. instead of min(std,C*iqr)
%       (iqr = interquartile range, C*iqr = robust stddev estimate)
%

% Copyright (C) 2003 Alexander Ihler; distributable under GPL -- see README.txt

  X = getPoints(npd);
  N = size(X,2);  dim = size(X,1);
  if (nargin<2) noIQR=0; end;

  Rg = .282095; Mg=1;                     % See ksizeCalcUseful for derivation
  Re = .6;      Me = .199994;             %   this is the canonical kernel adjustment
  Rl = .25;     Ml = 1.994473;            %   for product kernels of these types
  switch(npd.type),
      case 0, prop = 1.0;                 % Approximate; 1D prop = 1.059224; % Gaussian
      case 1, prop = ((Re/Rg)^dim / (Me/Mg)^2 )^(1/(dim+4)); % 1D prop = 2.344944; % Epanetchnikov
      case 2, prop = ((Rl/Rg)^dim / (Ml/Mg)^2 )^(1/(dim+4)); % 1D prop = 0.784452; % Laplacian
  end;

  sig = std(X,0,2);            % estimate sigma (standard)
  if (noIQR)
    h = prop*sig*N^(-1/(4+dim));
  else
    iqrSig = .7413*iqr(X')';     % find interquartile range sigma est.
    if (max(iqrSig)==0) iqrSig=sig; end;
    h = prop * min(sig,iqrSig) * N^(-1/(4+dim));
  end;

%  if (min(h) == 0) warning('Near-zero covariance => Kernel size set to 0'); end;

author: ShapeNet

repo: RenderForCNN


text:
function coords = combineCoords(coords1, coords2, action)
%
% coords = combineCoords(coords1, coords2, action)
%
% Combines coords1 and coords2, removing duplicates, used for
% example to merge ROI coordinates into one big ROI.
%
% Possible values of action:
%	'intersect' / 'and': Take only coordinates that lie in coords1 AND
%						 coords2.
%	'union' / 'or': Take coordinates that lie in coords1 OR coords2.
%
%	'xor': Take coordinates that only lie in EITHER coords1 or coords2,
%			but not both.
%	'a not b' / 'setdiff': Take coordinates that lie in coords1 BUT NOT
%			coords2.
%
% coords, coords1, and coords2: 3xN arrays of (y,x,z) coordinates
% dims is size of volume
%
% rmk 10/30/98
% djh, 2/2001, updated to use intersect(coords1,coords2,'rows')
% ras, 04/10/05, deals w/ empty coords
% ras, 03/28/07, added some comments, more flexibility in specifying the
% action.

% Matlab functions work on rows, not cols
coords1 = coords1';
coords2 = coords2';

% made a little more complex, in case we
% run into empty sets of coords:
if isempty(coords1)
	switch lower(action)
      case {'intersection' 'intersect' 'and'}
         coords = [];
      case {'union' 'or'}
         coords = coords2;
      case 'xor'
         coords = coords2;
      case {'a not b' 'anotb' 'setdiff'}
         coords = [];
	end
elseif isempty(coords2)
	switch lower(action)
      case {'intersection' 'intersect' 'and'}
         coords = [];
      case {'union' 'or'}
         coords = coords1;
      case 'xor'
         coords = coords1;
      case {'a not b' 'anotb' 'setdiff'}
         coords = coords1;
	end
else
	switch lower(action)
      case {'intersection' 'intersect' 'and'}
         coords = intersect(coords1, coords2, 'rows');
      case {'union' 'or'}
         coords = union(coords1, coords2,'rows');
      case 'xor'
         coords = setxor(coords1, coords2, 'rows');
      case {'a not b' 'anotb' 'setdiff'}
         coords = setdiff(coords1, coords2, 'rows');
	end
end

% Transpose back to 3xN
coords = coords';

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft