text:
% ------------------------GOLDEN SECTION METHOD----------------------------
% -------------------------------------------------------------------------
% Copyright (c) 2009, Katarzyna Zarnowiec, all rights reserved
% mailto: katarzyna.zarnowiec@gmail.com
% -------------------------------------------------------------------------

figure; hold on;

a=0;                            % start of interval
b=2;                            % end of interval
epsilon=0.000001;               % accuracy value
iter= 50;                       % maximum number of iterations
tau=double((sqrt(5)-1)/2);      % golden proportion coefficient, around 0.618
k=0;                            % number of iterations

x1=a+(1-tau)*(b-a);             % computing x values
x2=a+tau*(b-a);

f_x1=f(x1);                     % computing values in x points
f_x2=f(x2);

plot(x1,f_x1,'rx')              % plotting x
plot(x2,f_x2,'rx')

while ((abs(b-a)>epsilon) && (k<iter))
    k=k+1;
    if(f_x1<f_x2)
        b=x2;
        x2=x1;
        x1=a+(1-tau)*(b-a);

        f_x1=f(x1);
        f_x2=f(x2);

        plot(x1,f_x1,'rx');
    else
        a=x1;
        x1=x2;
        x2=a+tau*(b-a);

        f_x1=f(x1);
        f_x2=f(x2);

        plot(x2,f_x2,'rx')
    end

    k=k+1;
end

% chooses minimum point
if(f_x1<f_x2)
    sprintf('x_min=%f', x1)
    sprintf('f(x_min)=%f ', f_x1)
    plot(x1,f_x1,'ro')
else
    sprintf('x_min=%f', x2)
    sprintf('f(x_min)=%f ', f_x2)
    plot(x2,f_x2,'ro')
end


text:
function [model,entriesNotAdded] = addFoodRxns2HM(model,foodRxns)
% Takes a set of food reactions to a whole-body metabolic model
%
% USAGE
%
%   Format:    [model,entriesNotAdded] = addFoodRxns2HM(model,foodRxns)
%   Example 1: [model,entriesNotAdded] = addFoodRxns2HM(model,'All Food')
%   Example 2: [model,entriesNotAdded] = addFoodRxns2HM(model,{'Food_EX_AvocadoSalad[d]','Food_EX_BakerPotatoes[d]'})

% INPUTS
%
%   model:     a whole-body metabolic model
%   foodRxns:  a set of instructions to tell the function which food
%               reactions should be added to the model. Valid foodRxns
%               entries are as follows:
%                   - 'All' : adds all food reactions in the vmh database
%                   including food category reactions
%                   - 'All Food': adds all food reactions in the vmh
%                   database but does not include food category reactions
%                   - 'All Cat': adds all food category reactions to the
%                   model
%                   - A cellular array containing a list of food
%                   reactions to include in the model.
%
% OUTPUTS
%
%   model:          returns an updated model with the desired food reactions
%   entriesNotAdded:    a list of reactions that were not in the vmh database
%                       and therefore were not added to the model
%
% AUTHORS
%   Bronson R. Weston 2022

% Set up foodItems array
load('fdTable.mat')
load('fdCategoriesTable.mat')
if ischar(foodRxns)
    if strcmp(foodRxns,'All')
        foodRxns=[fdTable.Properties.VariableNames(2:end),fdCategoriesTable.Properties.VariableNames(2:end)];
        foodRxns=strcat('Food_EX_',foodRxns);
        foodRxns=strcat(foodRxns,'[d]');
    elseif strcmp(foodRxns,'All Food')
        foodRxns=fdTable.Properties.VariableNames(2:end);
        foodRxns=strcat('Food_EX_',foodRxns);
        foodRxns=strcat(foodRxns,'[d]');
    elseif strcmp(foodRxns,'All Cat')
        foodRxns=fdCategoriesTable.Properties.VariableNames(2:end);
        foodRxns=strcat('Food_EX_',foodRxns);
        foodRxns=strcat(foodRxns,'[d]');
    else
        foodRxns={foodRxns};
    end
else
    foodRxns=foodRxns;
end
%Combine fdTable and fdCategoriesTable and make necessary adjustements
fdCategoriesTable.Var1=[];
foodTable=[fdTable,fdCategoriesTable];
pro_Dindex=find(contains(foodTable.Var1,'pro_D')); %for now, remove pro_D from table.
foodTable(pro_Dindex,:)=[];
EnergyIndex=find(contains(foodTable.Var1,'Energy in Kcal'));
foodTable.Var1(EnergyIndex)={'Energy_in_Kcal'};

% Remove food items that already have reactions present in model
[~,~,ib] = intersect(model.rxns,foodRxns);
foodRxns(ib)= [];

%Add 'Energy_in_Kcal[d]' metabolite if not present
if ~any(strcmp(model.mets,'Energy_in_Kcal[d]'))
    model=addMetabolite(model, 'Energy_in_Kcal[d]');
    model = addMultipleReactions(model, {'EX_DietEnergy'}, {'Energy_in_Kcal[d]'}, [-1], 'lb', 0, 'ub', 100000);
end

% Identify food reactions to add
augFoodRxns=regexprep(foodRxns,'Food_EX_','');
augFoodRxns=regexprep(augFoodRxns,'\[d\]','');

[~,ia,ib] = intersect(foodTable.Properties.VariableNames,augFoodRxns);
entriesNotAdded=foodRxns(setdiff(1:length(foodRxns),ib));
foodRxns=foodRxns(ib);

%Add all food category reactions
foodMetabolites= foodTable.Var1.';
foodMetabolites= strcat(foodMetabolites,'[d]');
sMatrix=-1*table2array(foodTable(:,ia));
model = addMultipleReactions(model, foodRxns, foodMetabolites, sMatrix, 'lb', zeros(1,length(foodRxns)), 'ub', zeros(1,length(foodRxns)));
end

author: opencobra

repo: cobratoolbox


text:
function lik = lik_liks(varargin)
%  LIKS creates a likelihood structure which is composed by many
%       different likelihoods.
%
%  Description:
%    LIK_LIK = LIK_LIKS(LIKELIHOODS, VALUE1, ...) creates a likelihood
%    structure that allows alternative likelihoods for latent variables.
%    The full likelihood is the product of independent likelihoods
%    (independent observations given the latent process),
%
%               __k     __
%      p(y|f) = ||      ||   L_j(y_i|f_i, th_j)
%                 j=1   i \in I_j
%
%      f = (f1' f2' ... fk')'
%      y = (y1' y2' ... yk')'
%
%    Here j is the index for different likelihoods and I_j is the index
%    vector telling for which observations likelihood j is used.
%
%    When using the Liks likelihood you must give the vector/matrix z as
%    an extra parameter to each function that requires also y. The matrix z
%    must include the covariates for independent likelihoods (see e.g.
%    LIK_POISSON) and a column telling the index of likelihood attached to
%    respective observation. For example:
%      lik = lik_liks('likelihoods', {lik1 lik2}, 'classVariables', 2);
%      gp = gp_set('lik', lik, 'cf', k, 'latent_method', 'EP');
%      z = [1 1 1 2 2 2]';
%      gp_optim(gp, x, y, 'z', z)
%    Here z is a vector telling that the first 3 observations are related
%    to likelihood 1 and the last 3 to likelihood 2.
%
%    Parameters for Liks likelihood are [default]
%      likelihoods      - array of likelyhood structures [ {} ]
%      classVariables   - a scalar telling which column of matrix z defines
%                         the likelihood indices
%
%    For the demonstration of the use of lik_liks see DEMO_MULTIVARIATEGP.
%
%  See also
%    GP_SET, LIK_*, PRIOR_*
%
% Copyright (c) 2011 Jaakko Riihimäki
% Copyright (c) 2011 Aki Vehtari
% Copyright (c) 2012 Ville Tolvanen
% Copyright (c) 2015-2017 Jarno Vanhatalo
% ------------- 2015-2017 Marcelo Hartmann
% Copyright (c) 2017 Ville Tolvanen

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip = inputParser;
  ip.FunctionName = 'LIK_LIKS';
  ip.addOptional('lik', [], @isstruct);
  ip.addParamValue('likelihoods', {}, @(x) ~isempty(x) && iscell(x));
  ip.addParamValue('classVariables', [], @(x) ~isempty(x) && isscalar(x));
  ip.parse(varargin{:});
  lik = ip.Results.lik;

  if isempty(lik)
      init = true;
      lik.type = 'Liks';
  else
      if ~isfield(lik, 'type') || ~isequal(lik.type, 'Liks')
          error('First argument does not seem to be a valid likelihood function structure')
      end
      init = false;
  end

  % Initialize likelihoods
  if init || ~ismember('likelihoods', ip.UsingDefaults)
      lik.liks = ip.Results.likelihoods;
  end

  % number of likelihood functions
  lik.nliks = length(lik.liks);
  if lik.nliks < 2
      error('use one likelihood structure')

  else
      % column of z where the classVariables can be found
      lik.classVariables = ip.Results.classVariables;
  end

  % Initialize prior structure
  % Even if there is no prior at all
  % see line 2065 in the file gpla_e.m, and you may see the light !!!
  if init
      lik.p = [];
  end

  if init
      % Set the function handles to the subfunctions
      lik.fh.pak = @lik_liks_pak;
      lik.fh.unpak = @lik_liks_unpak;
      lik.fh.lp = @lik_liks_lp;
      lik.fh.lpg = @lik_liks_lpg;
      lik.fh.ll = @lik_liks_ll;
      lik.fh.llg = @lik_liks_llg;
      lik.fh.llg2 = @lik_liks_llg2;
      lik.fh.llg3 = @lik_liks_llg3;
      lik.fh.tiltedMoments = @lik_liks_tiltedMoments;
      lik.fh.siteDeriv = @lik_liks_siteDeriv;
      lik.fh.predy = @lik_liks_predy;
      lik.fh.invlink = @lik_liks_invlink;
      lik.fh.recappend = @lik_liks_recappend;
  end

end

function [w, s, h] = lik_liks_pak(lik)
% LIK_LIKS_PAK  Combines each likelihood parameters into one vector.
%
%  Description:
%    W = LIK_LIKS_PAK(LIK) takes a likelihood structure LIK and
%    combines the parameters into a single row vector W. This is a
%    mandatory subfunction used for example in energy and gradient
%    computations.
%
%   See also
%   LIK_LIKS_UNPAK, GP_PAK

 w = []; s = {}; h = [];

 for j = 1:lik.nliks
     lik_j = lik.liks{j};
     [wj, sj, hj] = lik_j.fh.pak(lik_j);
     w = [w wj];
     s = [s; sj];
     h = [h hj];

 end

end

function [lik, w] = lik_liks_unpak(lik, w)
% LIK_LIKS_UNPAK  Extracts each likelihood parameters from the vector.
%
%  Description:
%    [LIK, W] = LIK_LIKS_UNPAK(W, LIK) takes each likelihood
%    structure inside LIK and extracts the parameters from the vector W
%    to the whole LIK structure. This is a mandatory subfunction used
%    for example in energy and gradient computations.
%
%   See also:
%   LIK_LIKS_PAK, GP_UNPAK

% Assignment is inverse parameter transformation

 for j = 1:lik.nliks
     lik_j = lik.liks{j};
     [lik_j, w] = lik_j.fh.unpak(lik_j, w);
     lik.liks{j} = lik_j;

 end

end

function lp = lik_liks_lp(lik, varargin)
% LIK_LIKS_LP  log(prior) of each the likelihood parameters
%
%  Description:
%    LP = LIK_LIKS_LP(LIK) takes the a likelihood structure LIK and
%    returns log(p(th)), where th collects the parameters. This
%    subfunction is needed when there are likelihood parameters.
%
%  See also:
%    LIK_LIKS_LLG, LIK_LIKS_LLG2, LIK_LIKS_LLG3, GPLA_E

 % If there is prior for the likelihood parameters
 lp = 0;

 for j = 1:lik.nliks
     if ~isempty(lik.liks{j}.fh.pak(lik.liks{j}))
         lp = lp + lik.liks{j}.fh.lp(lik.liks{j});
     end
 end

end

function lpg = lik_liks_lpg(lik)
% LIK_LIKS_LPG  dlog(prior)/dth of each the likelihood parameters th
%
%  Description:
%    E = LIK_LIKS_LPG(LIK) takes a likelihood structure LIK and
%    returns d log(p(th))/dth for each likelihood function,
%    where th collects the parameters.
%    This subfunction is needed when there are likelihood parameters.
%
%  See also:
%    LIK_LIKS_LLG, LIK_LIKS_LG3, LIK_LIKS_LLG2, GPLA_G

 lpg = [];

 for j = 1:lik.nliks
     if ~isempty(lik.liks{j}.fh.pak(lik.liks{j}))
         lpg = [lpg lik.liks{j}.fh.lpg(lik.liks{j})];
     end
 end

end

function ll = lik_liks_ll(lik, y, ff, z)
% LIK_LIKS_LL log-likelihood
%
%  Description:
%    LL = LIK_LIKS_LL(LIK, Y, F) takes a likelihood structure LIK.
%    Returns the log-likelihood, sum_{i=1}^{k} (log p_i(y|f)).
%    This subfunction is needed when using Laplace approximation
%    or MCMC for inference with non-Gaussian likelihoods. This
%    subfunction is also used in information criteria (DIC, WAIC)
%    computations.
%
%  See also:
%    LIK_LLG, LIK_LLG3, LIK_LLG2, GPLA_E

 n  = size(y, 1);
 indClass = 1:size(z,2)==lik.classVariables;
 zi = z(:, indClass);
 z  = z(:, ~indClass);

 indj  = unique(zi);
 nind = numel(indj);

 if n ~= numel(zi)
     error('row-length of y and z are different')
 end

 f = ff(:);
 ll = 0;
 for j = 1:nind
     ind = zi==indj(j);
     likj = lik.liks{indj(j)};

     yj = y(ind);
     fj = f(ind);
     if isempty(z)
         zj = z;
     else
         zj = z(ind);
     end
     ll = ll + likj.fh.ll(likj, yj, fj, zj);

 end

end

function llg = lik_liks_llg(lik, y, ff, param, z)
% LIK_LIKS_LLG  Gradient of the log-likelihood
%
%  Description:
%    LLG = LIK_LIKS_LLG(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, Returns the gradient of the log likelihood
%    with respect to PARAM for each likelihood function. At the moment PARAM
%    can be 'param' or 'latent'. This subfunction is needed when using
%    Laplace approximation or MCMC for inference with non-Gaussian
%    likelihoods.
%
%  See also:
%    LIK_LIKS_LL, LIK_LIKS_LLG2, LIK_LIKS_LLG3, GPLA_E

 n  = size(y, 1);
 indClass = 1:size(z,2)==lik.classVariables;
 zi = z(:, indClass);
 z  = z(:, ~indClass);

 indj  = unique(zi);
 nind = numel(indj);

 if n ~= numel(zi)
     error('row-length of y and z are different')
 end

 f = ff(:);

 switch param
     case 'param'
         llg = [];
     case 'latent'
         llg=zeros(size(f));
 end
 for j = 1:nind
     ind = zi==indj(j);
     likj = lik.liks{indj(j)};

     yj = y(ind);
     fj = f(ind);
     if isempty(z)
         zj = z;
     else
         zj = z(ind);
     end

     switch param
         case 'param'
             if ~isempty(lik.liks{j}.fh.pak(likj))
                 llg = [llg likj.fh.llg(likj, yj, fj, param, zj)];
             end
         case 'latent'
             llg(ind) = likj.fh.llg(likj, yj, fj, param, zj);
     end
 end

end

function llg2 = lik_liks_llg2(lik, y, ff, param, z)
% LIK_LIKS_LLG2  Second gradients of the log-likelihood
%
%  Description
%    LLG2 = LIK_LIKS_LLG2(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, Returns the hessian of the log likelihood
%    with respect to PARAM. At the moment PARAM can be only
%    'latent'. LLG2 is a vector with diagonal elements of the
%    Hessian matrix (off diagonals are zero). This subfunction
%    is needed when using Laplace approximation or EP for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_LIKS_LL, LIK_LIKS_LLG, LIK_LIKS_LLG3, GPLA_E

 n  = size(y, 1);
 indClass = 1:size(z,2)==lik.classVariables;
 zi = z(:, indClass);
 z  = z(:, ~indClass);

 indj = unique(zi);
 nind = numel(indj);

 if n ~= numel(zi)
     error('row-length of y and z are different')
 end

 f = ff(:);

 nlikpar = length(lik.fh.pak(lik));
 z0 = zeros(n, nlikpar);
 aux(1) = 0;

 switch param
     case 'latent'
         llg2=zeros(size(f));
     case 'latent+param'
         llg2 = [];
 end
 for j = 1:nind
     ind = zi==indj(j);
     likj = lik.liks{indj(j)};

     yj = y(ind);
     fj = f(ind);
     if isempty(z)
         zj = z;
     else
         zj = z(ind);
     end

     switch param
         case 'param'

         case 'latent'
             llg2(ind) = likj.fh.llg2(likj, yj, fj, param, zj);

         case 'latent+param'
             if ~isempty(likj.fh.pak(likj))
                 % take the column vectors
                 llg2_tmp = likj.fh.llg2(likj, yj, fj, param, zj);

                 % auxiliar indexes
                 aux(end + 1) = aux(end) + size(likj.fh.pak(likj), 2);

                 % auxiliar matrices for derivatives w.r.t. parameters in
                 % the specific likelihood j
                 z0(ind, (aux(end - 1) + 1) : aux(end)) = llg2_tmp;
                 llg2 = [llg2 z0(:, aux(end - 1) + 1 : aux(end))];
                 z0(ind, (aux(end - 1) + 1) : aux(end)) = 0;

             end
     end

 end

end

function llg3 = lik_liks_llg3(lik, y, ff, param, z)
% LIK_LIKS_LLG3  Third gradients of the log likelihood
%
%  Description:
%    LLG3 = LIK_LIKS_LLG3(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, returns the third gradients of the log
%    likelihood with respect to PARAM. At the moment PARAM can be
%    only 'latent'. LLG3 is a vector with third gradients. This
%    subfunction is needed when using Laplace approximation for
%    inference with non-Gaussian likelihoods.
%
%  See also:
%    LIK_LIKS_LL, LIK_LIKS_LLG, LIK_LIKS_LLG2, GPLA_E, GPLA_G

 n  = size(y, 1);
 indClass = 1:size(z,2)==lik.classVariables;
 zi = z(:, indClass);
 z  = z(:, ~indClass);

 indj  = unique(zi);
 nind = numel(indj);

 if n ~= numel(zi)
     error('row-length of y and z are different')
 end

 f = ff(:);

 switch param
     case 'latent'
         llg3=zeros(size(f));
     case 'latent2+param'
         llg3 = [];
 end
 % auxiliar matrix for derivatives of parameters w.r.t. many likelihoods
 nlikpar = length(lik.fh.pak(lik));
 z0 = zeros(n, nlikpar);
 aux(1) = 0;

 for j = 1:nind
     switch param
         case 'param'

         case 'latent'
             ind = zi==indj(j);
             likj = lik.liks{indj(j)};

             yj = y(ind);
             fj = f(ind);
             if isempty(z)
                 zj = z;
             else
                 zj = z(ind);
             end

             llg3(ind) = likj.fh.llg3(likj, yj, fj, param, zj);

         case 'latent2+param'
             if ~isempty(lik.liks{indj(j)}.fh.pak(lik.liks{indj(j)}))
                 % take indexes and respective observations for specific likelihood
                 ind = zi==indj(j);
                 likj = lik.liks{indj(j)};

                 yj = y(ind);
                 fj = f(ind);
                 if isempty(z)
                     zj = z;
                 else
                     zj = z(ind);
                 end

                 % take the column vectors
                 llg3_tmp = likj.fh.llg3(likj, yj, fj, param, zj);

                 % auxiliar indexes
                 aux(end + 1) = aux(end) + size(likj.fh.pak(likj), 2);

                 % auxiliar matrices for derivatives w.r.t parameters in
                 % the specific likelihood j
                 z0(ind, (aux(end - 1) + 1) : aux(end)) = llg3_tmp;
                 llg3 = [llg3 z0(:, aux(end - 1) + 1 : aux(end))];
                 z0(ind, (aux(end - 1) + 1) : aux(end)) = 0;

             end

     end

 end

end

function [logM_0, m_1, sigm2hati1] = lik_liks_tiltedMoments(lik, y, i1, sigm2_i, myy_i, z)
% LIK_LIKS_TILTEDMOMENTS  Returns the marginal moments for EP algorithm
%
%  Description
%    [M_0, M_1, M2] = LIK_LIKS_TILTEDMOMENTS(LIK, Y, I, S2,
%    MYY) takes a likelihood structure LIKS, the observation Y, index
%    I, cavity variance S2 and mean MYY. Returns the zeroth
%    moment M_0, mean M_1 and variance M_2 of the posterior
%    marginal (see Rasmussen and Williams (2006): Gaussian
%    processes for Machine Learning, page 55). This subfunction
%    is needed when using EP for inference with non-Gaussian
%    likelihoods.
%
%  See also
%    GPEP_E

 n  = size(y, 1);
 indClass = 1:size(z,2)==lik.classVariables;
 zi = z(:, indClass);
 z  = z(:, ~indClass);

 indj = unique(zi);
 nind = numel(indj);

 if n ~= numel(zi)
     error('row-length of y and z are different')
 end

 logM_0 = zeros(n, 1);
 m_1 = zeros(n, 1);
 sigm2hati1 = zeros(n, 1);

 for j = 1:nind
     likj = lik.liks{indj(j)};
     ind = zi==indj(j);
     if numel(sigm2_i)>1
         yj = y(ind);
         if isempty(z)
             zj = z;
         else
             zj = z(ind);
         end
         sigm2_ij = sigm2_i(ind);
         myy_ij = myy_i(ind);
         [logM_0(ind), m_1(ind), sigm2hati1(ind)] = ...
             likj.fh.tiltedMoments(likj, yj, 1:length(yj), sigm2_ij, myy_ij, zj);
     else
         if any(find(ind)==i1)
             [logM_0, m_1, sigm2hati1] = ...
                 likj.fh.tiltedMoments(likj, y, i1, sigm2_i, myy_i, z);
         end
     end

 end

end

function [g_i] = lik_liks_siteDeriv(lik, y, i1, sigm2_i, myy_i, z)
% LIK_LIKS_SITEDERIV  Evaluate the expectation of the gradient
%                 of the log likelihood term with respect
%                 to the multiple-likelihood parameters for EP
%
%  Description [M_0, M_1, M2] =
%    LIK_LIKS_SITEDERIV(LIK, Y, I, S2, MYY, Z) takes a
%    likelihood structure LIKS, incedence counts Y, expected
%    counts Z, index I and cavity variance S2 and mean MYY.
%    Returns E_f [d log p(y_i|f_i) /d a], where a is the
%    likelihood parameter and the expectation is over the
%    marginal posterior. This is done for all likelihoods.
%    This term is needed when evaluating the
%    gradients of the marginal likelihood estimate Z_EP with
%    respect to the likelihood parameters (see Seeger (2008):
%    Expectation propagation for exponential families). This
%    subfunction is needed when using EP for inference with
%    non-Gaussian likelihoods and there are likelihood parameters.
%
%  See also
%    GPEP_G

 n  = size(y, 1);
 indClass = 1:size(z,2)==lik.classVariables;
 zi = z(:, indClass);
 z  = z(:, ~indClass);

 indj = unique(zi);
 nind = numel(indj);

 if n ~= numel(zi)
     error('row-length of y and z are different')
 end

 nlikpar = length(lik.fh.pak(lik));
 g_i = zeros(1, nlikpar);
 aux = 0;

 for j = 1:nind
     if ~isempty(lik.liks{indj(j)}.fh.pak(lik.liks{indj(j)}))
         % auxiliar indexes for paramters
         aux = aux + 1;

         % indexes for that specific likelihood
         ind = find(zi==indj(j));
         likj = lik.liks{indj(j)};

         if any(ind == i1)
             % !!!! if some specific likelihood has more than one parameter this will not work
             g_i(1, aux) = likj.fh.siteDeriv(likj, y, i1, sigm2_i, myy_i, z);
         end
     end

 end

end

function [lpy, Ey, Vary] = lik_liks_predy(lik, Ef, Varf, yt, zt)
% LIK_LIKS_PREDY  Returns the predictive mean, variance and density of y
%
%  Description:
%    LPY = LIK_LIKS_PREDY(LIK, EF, VARF YT, ZT)
%    Returns logarithm of the predictive density PY of YT, that is
%    p(yt | zt) = \int p(yt | f, zt) p(f|y) df.
%    This requires the observations YT.
%    This subfunction is needed when computing posterior predictive
%    distributions for future observations.
%
%    [LPY, EY, VARY] = LIK_LIKS_PREDY(LIK, EF, VARF) takes a
%    likelihood structure LIK, posterior mean EF and posterior
%    Variance VARF of the latent variable and returns the
%    posterior predictive mean EY and variance VARY of the
%    observations related to the latent variables. This subfunction
%    is needed when computing posterior predictive distributions for
%    future observations.
%
%
%  See also:
%    GPLA_PRED, GPEP_PRED, GPMC_PRED1

% number of values to predict;
n = size(Ef, 1);
indClass = 1:size(zt,2)==lik.classVariables;
zi = zt(:, indClass);
zt  = zt(:, ~indClass);

% check some conditions
if ~issorted(zi)
    error('you need to give the class variable increasing downwards');
end

if max(zi) > lik.nliks
    error('more classes than the number of classes you are trying to model');
end

% getting the information of the classes in the data
indj = unique(zi);
nind = numel(indj);

% log-density
lpy = zeros(n, 1);
if nargout > 1
    Ey = zeros(n, 1);
    Vary = zeros(n, 1);

    for j = 1:nind
        ind = zi==indj(j);
        likj = lik.liks{indj(j)};

        if numel(yt) ~= 0
           [lpy(ind), Ey(ind), Vary(ind)] = ...
           likj.fh.predy(likj, Ef(ind), Varf(ind), yt(ind), zt(ind));
        else
           [~, Ey(ind), Vary(ind)] = likj.fh.predy(likj, Ef(ind), Varf(ind), yt, zt(ind));
        end
    end
else
    for j = 1:nind
        ind = zi==indj(j);
        likj = lik.liks{indj(j)};

        if isempty(zt)
            lpy(ind) = likj.fh.predy(likj, Ef(ind), Varf(ind), yt(ind), []);
        else
            lpy(ind) = likj.fh.predy(likj, Ef(ind), Varf(ind), yt(ind), zt(ind));
        end
    end
end

end

function mu = lik_liks_invlink(lik, f, z)
%LIK_LIKS_INVLINK  Returns values of inverse link function
%
%  Description
%    P = LIK_LIKS_INVLINK(LIK, F) takes a likelihood structure LIK and
%    latent values F and returns the values MU of inverse link function.
%    This subfunction is needed when using gp_predprctmu.
%
%     See also
%     LIK_POISSON_LL, LIK_POISSON_PREDY

 n  = size(f, 1);
 indClass = 1:size(z,2)==lik.classVariables;
 zi = z(:, indClass);
 z  = z(:, ~indClass);

 indj = unique(zi);
 nind = numel(indj);

 if n ~= size(zi,1)
     error('row-length of f and z are different')
 end

 for j = 1:nind
     likj = lik.liks{indj(j)};
     ind = zi==indj(j);
     fj = f(ind,:);
     if isempty(z)
         zj = z;
     else
         zj = z(ind);
     end

     mu(ind,:) = likj.fh.invlink(likj, fj, zj);
 end

end

function reclik = lik_liks_recappend(reclik, ri, lik)
% RECAPPEND Append the parameters to the record
%
%  Description:
%    RECLIK = LIK_LIKS_RECAPPEND(RECLIK, RI, LIK) takes a
%    likelihood record structure RECLIK, record index RI and
%    likelihood structure LIK with the current MCMC samples of
%    the parameters. Returns RECLIK which contains all the old
%    samples and the current samples from LIK. This subfunction
%    is needed when using MCMC sampling (gp_mc).
%
%  See also:
%    GP_MC

 if nargin == 2
     % Initialize the record
     reclik.type = 'Liks';

     % Initialize parameters
     nliks = length(ri.liks);
     for i = 1:nliks
         lik_i = ri.liks{i};
         reclik.liks{i} = lik_i.fh.recappend([], ri.liks{i});
     end

     % Set the function handles
     reclik.fh.pak = @lik_liks_pak;
     reclik.fh.unpak = @lik_liks_unpak;
     reclik.fh.lp = @lik_liks_lp;
     reclik.fh.lpg = @lik_liks_lpg;
     reclik.fh.ll = @lik_liks_ll;
     reclik.fh.llg = @lik_liks_llg;
     reclik.fh.llg2 = @lik_liks_llg2;
     reclik.fh.llg3 = @lik_liks_llg3;
     reclik.fh.tiltedMoments = @lik_liks_tiltedMoments;
     reclik.fh.siteDeriv = @lik_liks_siteDeriv;
     reclik.fh.predy = @lik_liks_predy;
     reclik.fh.recappend = @lik_liks_recappend;

     if isfield(ri, 'classVariables')
         reclik.classVariables = ri.classVariables;
         reclik.nliks = ri.nliks;
     end

 else
     % Append to the record
     % Loop over all of the likelihood functions
     nliks = length(lik.liks);
     for i = 1:nliks;
         lik_i = lik.liks{i};
         reclik.liks{i} = lik_i.fh.recappend(reclik.liks{i}, ri, lik_i);
     end

 end
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function [pf,f]=lpccc2pf(cc,np,nc,c0)
%LPCCC2PF Convert complex cepstrum to power spectrum PF=(CC,NP,NC)
%
%  Inputs: cc(nf,n)     Complex ceptral coefficients excluding c(0), one frame per row
%          np           Size of output spectrum is np+1 [n]
%                       Alternatively, np can be a vector of output frequencies in the range 0 to 0.5
%          nc           Highest cepstral coefficient to use [np or, if np is a vector, n]
%                       Set nc=-1 to use n coefficients
%          c0(nf,1)     Cepstral coefficient cc(0) [0]
%
% Outputs: pf(nf,np+1)  Power spectrum from DC to Nyquist
%          f(1,np+1)    Normalized frequencies (0 to 0.5)
%
% The "complex cepstral coefficients", cc(n), are the inverse discrete-time Fourier transform
% of the log of the complex-valued spectrum. The cc(n) are real-valued and, for n<0, cc(n)=0.
% The "real cepstral coeffcients", rc(n), are the inverse discrete-time Fourier transform
% of the log of the magnitude spectrum; rc(0)=cc(0) and rc(n)=0.5*cc(n) for n~=0.
% For highest speed, choose np to be a power of 2.

%      Copyright (C) Mike Brookes 1998-2014
%      Version: $Id: lpccc2pf.m 5026 2014-08-22 17:47:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nf,mc]=size(cc);
if nargin<2 || ~numel(np)
    if nargout
        np=mc;
    else
        np=128;
    end
end
if nargin>=3 && numel(nc)==1 && nc==-1 nc=mc; end
if nargin<4 || ~numel(c0) c0=zeros(nf,1); end
if numel(np)>1 || np(1)<1
    if nargin<3 || ~numel(nc) nc=mc; end
    f=np(:)';
    if nc==mc
        pf=exp(2*[c0 cc]*cos(2*pi*(0:mc)'*f));
    else
        pf=exp(2*[c0 lpccc2cc(cc,nc)]*cos(2*pi*(0:nc)'*f));
    end
else
    if nargin<3 || ~numel(nc) nc=np; end
    if nc==mc
        pf=exp(2*real(rfft([c0 cc].',2*np).'));
    else
        pf=exp(2*real(rfft([c0 lpccc2cc(cc,nc)].',2*np).'));
    end
    f=linspace(0,0.5,np+1);
end
if ~nargout
    plot(f,db(pf.')/2);
    xlabel('Normalized frequency f/f_s');
    ylabel('Gain (dB)');
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
%% Image Warping
% A demo program shows how perspective transformation applied on an image.
% Based on a sample code from
% <http://study.marearts.com/2015/03/image-warping-using-opencv.html MareArts blog>.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.3.1/samples/cpp/warpPerspective_demo.cpp>
%

function varargout = warp_perspective_demo_gui(im)
    % load source image and set initial ROI corners (4 points clockwise)
    roi = [];
    if nargin < 1
        if true
            im = fullfile(mexopencv.root(), 'test', 'books_right.jpg');
            roi = [360 109; 532 138; 460 417; 317 338];
        else
            im = fullfile(mexopencv.root(), 'test', 'box_in_scene.png');
            roi = [243 188; 328 151; 388 294; 319 351];
        end
        img = cv.imread(im, 'Color',true);
    elseif ischar(im)
        img = cv.imread(im, 'Color',true);
    else
        img = im;
    end
    if isempty(roi)
        roi = bsxfun(@rdivide, [size(img,2) size(img,1)], ...
            [1.7 4.2; 1.15 3.32; 1.33 1.1; 1.93 1.36]);
    end

    % create the UI and hook event handlers
    h = buildGUI(img, roi);
    if nargout > 0, varargout{1} = h; end
    opts = {'Interruptible','off', 'BusyAction','cancel'};
    set(h.fig, 'CloseRequestFcn',@(~,~) delete(h.fig), ...
        'WindowKeyPressFcn',@onType, opts{:});
    set(h.fig(1), 'WindowButtonDownFcn',@onMouseDown, opts{:});

    %% Callback Functions

    function redraw()
        %REDRAW  Apply transformation using current ROI and display results

        % reapply transformation
        out1 = drawROI(img, roi);
        [out2, sz] = warpROI(img, roi);

        % show results and ajust second plot to fit image
        set(h.img(1), 'CData',out1);
        set(h.img(2), 'CData',out2, 'XData',[1 sz(1)], 'YData',[1 sz(2)]);
        set(h.ax(2), 'XLim',[0 sz(1)]+0.5, 'YLim',[0 sz(2)]+0.5);
        pos = get(h.fig(2), 'Position');
        set(h.fig(2), 'Position',[pos(1:2) sz]);
        drawnow;
    end

    function onType(~,e)
        %ONTYPE  Event handler for key press on figure

        switch e.Key
            case 'h'
                helpdlg({
                    'Use your mouse to select a point and move it'
                    'to see transformation changes.'
                    ''
                    'Hot keys:'
                    'h - this help dialog'
                    'q - quit the program'
                    'r - change order of points to rotate transformation'
                    'i - change order of points to invert transformation '
                });
            case {'q', 'escape'}
                close(h.fig);
            case 'r'
                roi = circshift(roi, -1);
                redraw();
            case 'i'
                roi = roi([2 1 4 3],:);
                redraw();
        end
    end

    function onMouseDown(~,~)
        %ONMOUSEDOWN  Event handler for mouse down on figure

        % hit-test for closest ROI corner
        pt = getCurrentPoint(h.ax(1));
        d = sum(abs(bsxfun(@minus, roi, pt)), 2);
        [mn, pt_idx] = min(d);
        if mn < 20
            % attach event handlers, and change mouse pointer
            set(h.fig(1), 'Pointer','cross', ...
                'WindowButtonMotionFcn',{@onMouseMove, pt_idx}, ...
                'WindowButtonUpFcn',@onMouseUp);
        end
    end

    function onMouseMove(~,~,idx)
        %ONMOUSEMOVE  Event handler for mouse move on figure

        % move specified ROI corner
        pt = getCurrentPoint(h.ax(1));
        roi(idx,:) = pt;
        redraw();
    end

    function onMouseUp(~,~)
        %ONMOUSEUP  Event handler for mouse up on figure

        % detach event handlers, and restore mouse pointer
        set(h.fig(1), 'Pointer','arrow', ...
            'WindowButtonMotionFcn','', ...
            'WindowButtonUpFcn','');
    end
end

%% Helper Functions

function out = drawROI(img, pts)
    %DRAWROI  Show ROI corners with labels

    labels = {'TL', 'TR', 'BR', 'BL'};
    out = cv.polylines(img, pts, 'Closed',true, 'Color',[0 0 255], 'Thickness',2);
    out = cv.circle(out, pts, 5, 'Color',[0 255 0], 'Thickness',3);
    out = cv.putText(out, labels, pts, ...
        'FontScale',0.8, 'Color',[255 0 0], 'Thickness',2);
end

function [out, dsz] = warpROI(img, srcPts)
    %WARPROI  Compute warped image from ROI

    % map ROI corners to corresponding destination rectangle
    len = diff(srcPts([1:end 1],:), 1, 1);  % [TR-TL; BR-TR; BL-BR; TL-BL]
    len = cellfun(@norm, num2cell(len, 2)); % lengths of sides
    w = max(len(1), len(3));
    h = max(len(2), len(4));
    dstPts = [0 0; w 0; w h; 0 h];

    % compute homography between points and apply persepective transformation
    dsz = round([w h]);
    H = cv.findHomography(srcPts, dstPts);
    out = cv.warpPerspective(img, H, 'DSize',dsz);
end

function p = getCurrentPoint(ax)
    %GETCURRENTPOINT  Retrieve current mouse location

    p = get(ax, 'CurrentPoint');
    p = p(1,1:2) - 1;
end

function h = buildGUI(img, roi)
    %BUILDGUI  Creates the UI

    % apply initial perspective transformation
    out1 = drawROI(img, roi);
    out2 = warpROI(img, roi);
    sz1 = size(out1);
    sz2 = size(out2);

    % properties
    fprops = {'Menubar','none', 'Resize','on'};
    aprops = {'Units','normalized', 'Position',[0 0 1 1]};
    h = struct();

    % show input image + ROI corners
    h.fig(1) = figure('Name','Image', ...
        'Position',[100 200 sz1(2) sz1(1)], fprops{:});
    h.ax(1) = axes('Parent',h.fig(1), aprops{:});
    if mexopencv.isOctave()
        h.img(1) = imshow(out1);
    else
        h.img(1) = imshow(out1, 'Parent',h.ax(1));
    end

    % show warped image
    h.fig(2) = figure('Name','Warped', ...
        'Position',[200+sz1(2) 200 sz2(2) sz2(1)], fprops{:});
    h.ax(2) = axes('Parent',h.fig(2), aprops{:});
    if mexopencv.isOctave()
        h.img(2) = imshow(out2);
    else
        h.img(2) = imshow(out2, 'Parent',h.ax(2));
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function value = sum(SO3F, varargin)
% Calculates the integral of a SO3FunHarmonic based on
%
% $$ v = \int f(R) dR $$
%
% with $vol(SO(3)) = \int_{SO(3)} 1 dR = 8\pi^2$.
%
% If there is a second argument it sums up along a specified dimension of a
% vector-valued SO3FunHarmonic.
%
% Syntax
%   value = sum(SO3F)
%   SO3F = sum(SO3F, d)
%
% Input
%  SO3F - @SO3FunHarmonic
%  d    - dimension to take the sum value over
%
% Output
%  SO3F  - @SO3FunHarmonic
%  value - double
%
% Description
%
% SO3F is a 3x3 SO3Fun
% sum(SO3F) returns a 3x3 matrix with the integrals of each function
% sum(SO3F, 1) returns a 1x3 SO3Fun which contains the pointwise sums along the first dimension
%

if nargin == 1
  value = 8*pi^2*mean(SO3F,varargin{:});
else
  SO3F.fhat = sum(SO3F.fhat, varargin{1}+1);
  value = SO3F;
end

end

author: mtex-toolbox

repo: mtex

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


text:
function this = subsasgn(this, selector, value)
%SUBSASGN   Modify parameters of CLASSIC+NL optical flow algorithm
%   THIS = SUBSASGN(THIS, SEL, VAL) sets parameter of a HS optical
%   flow algorithm THIS using given selector SEL to value VAL.  The
%   modified object is returned.
%   Available selectors:
%     - ".optical_flow":   Flow base object
%     - ".pyramid_levels": Levels of the image pyramid
%     - ".warping_mode":   'forward' for forward warping and 'backward'
%                          for backward warping
%     - ".solver":         Selector for the linear equation solver
%                          ('backslash' for the builtin solver and 'sor'
%                          for a MEX-based SOR solver)
%     - ".max_iters":      Maximum number of iterations at each pyramid
%                          level
%     - ".limit_update":   Flag whether to restrict the incremental flow
%                          at each level to [-1, 1]
%
% This is a member function of the class 'classic_nl_optical_flow'.
%
% Authors: Deqing Sun, Department of Computer Science, Brown University
% Contact: dqsun@cs.brown.edu
% $Date: $
% $Revision: $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  % Allow full access to all fields of the structure
  this = builtin('subsasgn', this, selector, value);

author: kristinbranson

repo: JAABA


text:
function [paramtransformsettings, names] = whiteKernExtractParamTransformSettings(kern)

% WHITEKERNEXTRACTPARAM Extract parameter transform settings from the WHITE kernel structure.
% FORMAT
% DESC Extract parameters from the white noise kernel structure into a
% vector of parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% FORMAT
% DESC Extract parameters and parameter names from the white noise
% kernel structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
% RETURN names : cell array of strings giving paramter names.
%
% SEEALSO whiteKernParamInit, whiteKernExpandParam, kernExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006
%
% COPYRIGHT : Jaakko Peltonen, 2011
%
% KERN

% The "white" kernel uses just one transformation, for the variance
% parameter of the kernel.

paramtransformsettings=cell(1);
paramtransformsettings{1}=kern.transforms(1).transformsettings;
if nargout > 1
  names = {'variance'};
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function bessel_j1_values_test ( )

%*****************************************************************************80
%
%% BESSEL_J1_VALUES_TEST demonstrates the use of BESSEL_J1_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BESSEL_J1_VALUES_TEST:\n' );
  fprintf ( 1, '  BESSEL_J1_VALUES stores values of \n' );
  fprintf ( 1, '  the Bessel J1 function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X            J1(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = bessel_j1_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [z, logQ] = sampleStateSeq_RGS( logSoftEv, eta_ii, F_ii, ii, TargetPsi )
% Draw new assignments for the discrete hidden state sequence
%  for each time series object, under *restricted settings*.
% Uses a fast message backward, sample forwards algorithm.

ks = find( F_ii > 0 );
T  = size(logSoftEv,2);
K  = length(ks);

if K == 1
    z = ks(1)*ones( 1, Tii );
    continue;
end

pi_init = 1/K*ones(1,K);
pi      = eta_ii;

normC = max( logSoftEv, [], 1);
logSoftEv = bsxfun( @minus, logSoftEv, normC );
Lik = exp( logSoftEv );

SEED = randomseed();
randomseed( SEED+1 );
[z, logQ] = SampleHMMStateSeqWithQsC( pi, Lik, pi_init, SEED(1) );
z = ks(z);

end % main function

author: michaelchughes

repo: NPBayesHMM


text:
function p = myismember(a,A)
% MYISMEMBER Is 'a' an element of a set of positive integers? (much faster than built-in ismember)
% p = myismember(a,A)

%if isempty(A) | a < min(A) | a > max(A)  % slow

if length(A)==0
  p = 0;
elseif a < min(A)
  p = 0;
elseif a > max(A)
  p = 0;
else
  bits = zeros(1, max(A));
  bits(A) = 1;
  p = bits(a);
end
p = logical(p);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% Laplace Transform properties

% Differentiation in the time-domain

x=cos(t);

L=diff(x,t) ;
laplace(L,s)

X=laplace(x,s);
x0=1;
R=s*X-x0;
simplify(R)

author: Sable

repo: mcbench-benchmarks


text:
function displayKaplanMeier_HN_tempFig5(pathRF,nameOutcomes,nameSets,thresholds)
% thresholds: (optional). Defines the thresholds of the risk stratification groups.
%             Defaut is one threshold of 0.5 for two risk groups.

if nargin < 4
    thresholds = 0.5;
else
    if numel(thresholds) > 2
        error('Max 3 curves (risk groups) are allowed for the moment')
    end
    thresholds = thresholds./100;
end
if numel(thresholds) == 1
    symbols = {'-b','--r'};
elseif numel(thresholds) == 2
    symbols = {'-g','--b','.-r'};
end

startpath = pwd;
nOutcomes = numel(nameOutcomes);

cd(pathRF), load('testing')
thresholds = [0,thresholds,1];
nCurves = numel(thresholds) - 1;
figure
for o = 1:nOutcomes
    legendCell = cell(1,nCurves);
    nameOutcome = nameOutcomes{o}; fSet = nameSets{o};
    if strcmp(nameOutcome,'Locoregional')
        name = 'Recurrence-free probability';
    elseif strcmp(nameOutcome,'Distant')
        name = 'Metastasis-free probability';
    elseif strcmp(nameOutcome,'Death')
        name = 'Survival probability';
    elseif strcmp(nameOutcome,'DeathSign')
        name = 'Survival probability';
    end
    load(['testResultsRF_',fSet,'_',nameOutcome]) % results gets out of there
    prob = results.probResponse;
    if strcmp(nameOutcome,'DeathSign')
        outcome = testing.outcomes.Death;
        timeData =  testing.timeToEvents.Death;
    else
        outcome = testing.outcomes.(nameOutcome);
        timeData =  testing.timeToEvents.(nameOutcome);
    end
    censData = 1 - outcome;
    hFig = subplot(1,nOutcomes,o); hold on, maxTT = 0; handles = []; tables1 = cell(1,nCurves); tables2 = cell(1,nCurves); X = cell(1,nCurves);
    for c = 1:nCurves
        time = timeData(prob >= thresholds(c) & prob < thresholds(c+1));
        cens = censData(prob >= thresholds(c) & prob < thresholds(c+1));
        [tableT,tableTT,t,T,xcg,ycg,~] = kmplot([time,cens],0.05,1,0); tables1{c} = tableT; tables2{c} = tableTT; X{c} = [time,cens];
        S = stairs(t,T,symbols{c},'LineWidth',2); handles = [handles,S];
        h = plot(xcg,ycg,'k+','MarkerSize',6,'LineWidth',1);
        maxT = max(t);
        if maxT > maxTT
            maxTT = maxT;
        end
        legendCell{c} = [num2str(thresholds(c)),' <= prob_{RF} < ',num2str(thresholds(c+1))];
    end
%     if nCurves == 2
%         pVal = getPval(tables1{1},tables1{2},tables2{1},tables2{2},X{1},X{2});
%         annotation(hFig,'textbox',...
%         [0.2635 0.300000001206285 0.22142856536167 0.0642857130794299],...
%         'String',{['p - value = ',num2str(pVal)]},...
%         'LineStyle','none','FitBoxToText','off');
%     elseif nCurves == 3
%         % Between low-risk and medium-risk curves (see bottom-left value)
%         pVal = getPval(tables1{1},tables1{2},tables2{1},tables2{2},X{1},X{2});
%         annotation(hFig,'textbox',...
%         [0.2635 0.300000001206285 0.22142856536167 0.0642857130794299],...
%         'String',{['p - value = ',num2str(pVal)]},...
%         'LineStyle','none','FitBoxToText','off');
%
%         % Between medium-risk and high-risk curves (see bottom-right value)
%         pVal = getPval(tables1{2},tables1{3},tables2{2},tables2{3},X{2},X{3});
%         annotation(hFig,'textbox',...
%         [0.6135 0.300000001206285 0.22142856536167 0.0642857130794299],...
%         'String',{['p - value = ',num2str(pVal)]},...
%         'LineStyle','none','FitBoxToText','off');
%     end
    axis([0 maxTT 0 1.05]), axis square
    ind = strfind(nameOutcome,'Death');
    if ~isempty(ind)
        nameOutcome(ind:ind+4) = [];
        nameOutcome = ['Survival',nameOutcome];
    end
    set(gca,'FontSize',12)
    title([nameOutcome,'\newline(',nameSets{o},')'],'FontSize',14)
    xlabel('Time (days)','FontSize',14)
    ylabel(name,'FontSize',14)
    if nCurves == 3
        legendCell = {'Low-risk group','Medium-risk group','High-risk group'};
    else
    end
    legend(handles,legendCell,'Location','SouthWest')
    grid ON
    hold off
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [l,L_k,L_hm,L_beta] = invPinHoleAPlucker(Sk,hm,beta)

% INVPINHOLEAPLUCKER Retro-projects anchored plucker line
%   INVPINHOLEAPLUCKER(K,L,BETA) retro-projects the anchored Plucker line
%   from the homogeneous 2D line HM and a pin hole camera K=[u0;v0;au;av]
%   at the origin. BETA specifies the unobservable direction of the line.
%   BETA is a 2-vector expressed in the plane base given by
%   PLANEVEC2PLANEBASE.
%
%   [l,L_k,L_hm,L_beta] = ... returns Jacobians wrt K, HM and BETA.
%
%   See also INVPINHOLEPLUCKER.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if nargout == 1
    % Plucker line
    pl = invPinHolePlucker(Sk,hm,beta) ;

    % anchored Plucker line L
    l = anchorPlucker(pl,[0;0;0]);

else
    % Plk lin in Sensor Frame
    [pl, PL_k, PL_hm, PL_beta] = invPinHolePlucker(Sk,hm,beta) ;

    % L in sensor frame
    [l, L_pl] = anchorPlucker(pl,[0;0;0]);

    % chain rule
    L_k    = L_pl * PL_k;
    L_hm   = L_pl * PL_hm;
    L_beta = L_pl * PL_beta;

end

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

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
clear all; close all; clc

%% folder setup
saveFigFlag = 1;

outputDir = GetOutputDataDir;
saveDir1 = fullfile(outputDir,'multimotor_0228');
% saveDir2 = fullfile(outputDir,'motor_map_seed');
if ~exist(saveDir1, 'dir'), mkdir(saveDir1), end;
% if ~exist(saveDir2, 'dir'), mkdir(saveDir2), end;

%% init

hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);
% i_fish = 8;
% setappdata(hfig,'isMotorseed',0);

%% run fish
range_fish = GetFishRange;%[1:3,5:18];
% M_thres_reg = zeros(3,18);
% M_numTopCorr = zeros(1,18);
% M_motorseedRegs = cell(1,18);
% M_compareMotorCellNumber = zeros(2,18);

%%
stimflag = 'P';
ClusterIDs = [6,2];

% stimflag = []; % for default set
% ClusterIDs = [6,1];
M_stimrange = GetStimRange(stimflag);

range_fish_valid = [];
for i_fish = range_fish
    if ~isempty(M_stimrange{i_fish})
       range_fish_valid = [range_fish_valid,i_fish];
    end
end

tscriptstart = tic;
for i_fish = range_fish_valid
    %% load data for chosen stim range
    stimrange = M_stimrange{i_fish};
    [cIX,gIX_load,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs,stimrange);

    %% Method 1: stim+motor regression
    fishset = getappdata(hfig,'fishset');
    [~,~,reg_sens] = GetStimRegressor(stim,fishset,i_fish);
    [~,~,reg_motor] = GetMotorRegressor(behavior,i_fish);

    % cell based
    gIX = (1:length(cIX))';
    tic
    [stimcorr,motorcorr] = MotorSourceCorrelation(M,reg_sens,reg_motor);
    toc
    [fig1,fig2] = MultiMotorVisuals(hfig,stimcorr,motorcorr,cIX,gIX);

    saveDir = fullfile(saveDir1,'motor rank - cell based - 2D');
    figName = ['Fish' num2str(i_fish)];
    SaveFigureHelper(saveFigFlag, saveDir, figName,fig1);

    saveDir = fullfile(saveDir1,'motor rank - cell based - anat');
    figName = ['Fish' num2str(i_fish)];
    SaveFigureHelper(saveFigFlag, saveDir, figName,fig2);

    %% cluster based
    % plot 1
%     setappdata(hfig,'isMotorseed',0);
%     %     setappdata(hfig,'stimrange',1:2);
%
%     UpdateTimeIndex(hfig);
%     behavior = getappdata(hfig,'behavior');
%     [~,~,reg_motor] = GetMotorRegressor(behavior,i_fish);

    gIX = gIX_load;
    C = FindClustermeans(gIX,M);
    tic
    [stimcorr,motorcorr] = MotorSourceCorrelation(C,reg_sens,reg_motor);
    toc
    [fig1,fig2] = MultiMotorVisuals(hfig,stimcorr,motorcorr,cIX,gIX);

%     % plot 2
%     setappdata(hfig,'isMotorseed',1);
%     %     setappdata(hfig,'stimrange',1:2);
%
%     UpdateTimeIndex(hfig);
%     behavior = getappdata(hfig,'behavior');
%     [~,~,reg_motor] = GetMotorRegressor(behavior,i_fish);
%
%     gIX = gIX_load;
%     C = FindClustermeans(gIX,M);
%     tic
%     [stimcorr,motorcorr] = MotorSourceCorrelation(C,reg_sens,reg_motor);
%     toc
%     [fig3,fig4] = MultiMotorVisuals(hfig,stimcorr,motorcorr,cIX,gIX);
    %%
    saveDir = fullfile(saveDir1,'motor rank - cluster based - 2D');
    figName = ['Fish' num2str(i_fish)];
    SaveFigureHelper(saveFigFlag, saveDir, figName,fig1);

    saveDir = fullfile(saveDir1,'motor rank - cluster based - anat');
    figName = ['Fish' num2str(i_fish)];
    SaveFigureHelper(saveFigFlag, saveDir, figName,fig2);
    %%

% %  saveas(gcf, fn, 'png');
% %
% % %%
% %     [rawmotorcorr,IX] = max(corr(reg_motor',M'),[],1);
% %     thres_reg = 0.5;
% %     ix = find(rawmotorcorr>thres_reg);
    %% Method 2: stimAvr+motor regression
    % cluster based;
    gIX = gIX_load;
    C = FindClustermeans(gIX,M);
    reg_sens = GetStimAvrClusmean(hfig,C);
    [~,~,reg_motor] = GetMotorRegressor(behavior,i_fish);

    nUnit = size(C,1);
    stimcorr = zeros(nUnit,1);
    motorcorr = zeros(nUnit,1);
    tic
    for i = 1:size(C,1)
    [stimcorr(i),motorcorr(i)] = MotorSourceCorrelation(C(i,:),reg_sens(i,:),reg_motor);
    end
    toc

    [fig1,fig2] = MultiMotorVisuals(hfig,stimcorr,motorcorr,cIX,gIX);

    saveDir = fullfile(saveDir1,'stimAvr rank - cluster based - 2D');
    figName = ['Fish' num2str(i_fish)];
    SaveFigureHelper(saveFigFlag, saveDir, figName,fig1);

    saveDir = fullfile(saveDir1,'stimAvr rank - cluster based - anat');
    figName = ['Fish' num2str(i_fish)];
    SaveFigureHelper(saveFigFlag, saveDir, figName,fig2);

end
toc(tscriptstart)

author: xiuyechen

repo: FishExplorer


%==============================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% function v0 = getVelocityStartingGuess(omega,m,nt)
%
% returns starting guess for velocity, that is, zero vector of appropiate
% size, which depends on regularizer, spatial and temporal discretization
% size.
%
% Input:
%
%  omega   - spatial domain
%      m   - number of cells for velocity
%     nt   - number of time points
%
% Output:
%
%     v0   - vector of zeros. size depends on grid for regularizer and nt
%
% =========================================================================
function v0 = getVelocityStartingGuess(omega,m,nt)
if not(exist('nt','var')) || isempty(nt)
    nt = regularizer('get','nt');% look up if regularizer is configured
    if isempty(nt),  nt = 0; end;
end
dim = numel(omega)/2;
switch regularizer('get','grid')
    case 'nodal'
        v0 = zeros(dim*prod(m+1),1);
    case 'staggered'
        v0 = zeros(sum(prod(ones(dim,1)*m+eye(dim),2)),1);
    otherwise % assume cell-centered
        v0 = zeros(dim*prod(m),1);
end
v0 = repmat(v0,nt+1,1);


text:
function plotclustfir(infir,oXc)

colors = {'r' 'g' 'b' 'y' 'c' 'm' 'k'};
if max(oXc)>length(colors)
    error('not enough colors!');
end

for c=1:size(infir,1);
    for s=1:size(infir,2);
        fir(c,:,:,s)=infir{c,s};
    end
end

mfir=squeeze(mean(fir,4));
figure;
for c=1:max(oXc);
leg{c}=['clust',num2str(c)];
end

for n=1:size(fir,1);        % for each condition
    subplot(1,size(mfir,1),n);
    for c=1:max(oXc);
    hold on;
    plotfir=squeeze(mean(mfir(n,:,oXc==c),3));
    plotfir=squeeze(mfir(n,:,oXc==c));
    plot(plotfir,'color',colors{c},'LineWidth',3)
    end
    legend(leg)
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [disp_row, disp_col] = resp_newton(response, responsef, iterations, ky, kx, use_sz)

[max_resp_row, max_row] = max(response, [], 1);
[init_max_response, max_col] = max(max_resp_row, [], 2);
max_row_perm = permute(max_row, [2 3 1]);
col = max_col(:)';
row = max_row_perm(sub2ind(size(max_row_perm), col, 1:size(response,3)));

trans_row = mod(row - 1 + floor((use_sz(1)-1)/2), use_sz(1)) - floor((use_sz(1)-1)/2);
trans_col = mod(col - 1 + floor((use_sz(2)-1)/2), use_sz(2)) - floor((use_sz(2)-1)/2);
init_pos_y = permute(2*pi * trans_row / use_sz(1), [1 3 2]);
init_pos_x = permute(2*pi * trans_col / use_sz(2), [1 3 2]);
max_pos_y = init_pos_y;
max_pos_x = init_pos_x;

% pre-compute complex exponential
exp_iky = exp(bsxfun(@times, 1i * ky, max_pos_y));
exp_ikx = exp(bsxfun(@times, 1i * kx, max_pos_x));

% gradient_step_size = gradient_step_size / prod(use_sz);

ky2 = ky.*ky;
kx2 = kx.*kx;

iter = 1;
while iter <= iterations
    % Compute gradient
    ky_exp_ky = bsxfun(@times, ky, exp_iky);
    kx_exp_kx = bsxfun(@times, kx, exp_ikx);
    y_resp = mtimesx(exp_iky, responsef, 'speed');
    resp_x = mtimesx(responsef, exp_ikx, 'speed');
    grad_y = -imag(mtimesx(ky_exp_ky, resp_x, 'speed'));
    grad_x = -imag(mtimesx(y_resp, kx_exp_kx, 'speed'));
    ival = 1i * mtimesx(exp_iky, resp_x, 'speed');
    H_yy = real(-mtimesx(bsxfun(@times, ky2, exp_iky), resp_x, 'speed') + ival);
    H_xx = real(-mtimesx(y_resp, bsxfun(@times, kx2, exp_ikx), 'speed') + ival);
    H_xy = real(-mtimesx(ky_exp_ky, mtimesx(responsef, kx_exp_kx, 'speed'), 'speed'));
    det_H = H_yy .* H_xx - H_xy .* H_xy;

    % Compute new position using newtons method
    max_pos_y = max_pos_y - (H_xx .* grad_y - H_xy .* grad_x) ./ det_H;
    max_pos_x = max_pos_x - (H_yy .* grad_x - H_xy .* grad_y) ./ det_H;

    % Evaluate maximum
    exp_iky = exp(bsxfun(@times, 1i * ky, max_pos_y));
    exp_ikx = exp(bsxfun(@times, 1i * kx, max_pos_x));

    iter = iter + 1;
end
max_response = 1 / prod(use_sz) * real(mtimesx(mtimesx(exp_iky, responsef, 'speed'), exp_ikx, 'speed'));

% check for scales that have not increased in score
ind = max_response < init_max_response;
max_pos_y(ind) = init_pos_y(ind);
max_pos_x(ind) = init_pos_x(ind);
disp_row = (mod(max_pos_y(1,1,1) + pi, 2*pi) - pi) / (2*pi) * use_sz(1);
disp_col = (mod(max_pos_x(1,1,1) + pi, 2*pi) - pi) / (2*pi) * use_sz(2);
end

repo: AutoTrack


text:
 function z = fatrix2_embed(mask, dim, x)
%function z = fatrix2_embed(mask, dim, x)
%|
%| version of embed that handles an empty mask
%|
%| in
%|	x	[N L]
%|	dim	[]		size(mask) if mask is nonempty
%|	mask	[(dim)]		possibly empty
%|
%| out
%|	z	[(dim) L]
%|
%| Copyright 2010-12-21, Jeff Fessler, University of Michigan

if nargin < 3, ir_usage, end

if ~isempty(mask)
	z = embed(x, mask);
else
	z = reshape(x, [dim size(x,2)]);
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function fig_state= fig_set(varargin)
%FIG_SET - Setup a new figure
%
%Synposis:
% FIG_STATE= fig_set(FN)
% FIG_STATE= fig_set(FN, <OPT>)
% FIG_STATE= fig_set(<OPT>)
%
%Input:
% FN:  figure number
% OPT: property/value list or struct of options with fields/properties:
%  .Fn          - figure number, alternative to using the first argument
%  .Name        - name of figure window, default: ''
%  .Silent      - if true, doesn't bring figure to the foreground,
%                 default: false
%  .Hide        - if true, sets the figure 'Visible' property to false,
%                 default: false
%  .Toolsoff    - if true, turns off the toolbar of the figure window,
%                 default: true
%  .Clf         - if true, clears the figure, default: false
%  .Resize      - resize the automatically defined dimensions of the figure
%                 (width, height) by the factors [wf hf], default: [1 1]
%  .GridSize    - figures are successively arranged in a [M x N] grid that
%                 covers the whole screen, according to their figure
%                 number, default: [2 2]
%
%Output:
% FIG_STATE: structure with a 'WasVisible' field, to use with function
%            fig_publish

props= {'Fn'              1            '!INT'
        'Name'            ''           'CHAR'
        'Silent'          0            '!BOOL'
        'Hide'            0            '!BOOL'
        'ToolsOff'        1            '!BOOL'
        'Clf'             0            '!BOOL'
        'Resize'          [1 1]        '!DOUBLE[2]'
        'Square'          0            '!BOOL'
        'GridSize'        [2 2]        '!DOUBLE[2]'
        'Props'           {}           'PROPLIST'
        'ShiftUpwards'    1            '!BOOL'
        'DesktopBorder'   [0 30 0 0]   '!DOUBLE[4]'
        'WindowBorder'    [5 20]       '!DOUBLE[2]'
       };

if nargin>0 && isnumeric(varargin{1}),
  opt= opt_proplistToStruct(varargin{2:end});
  opt.Fn= varargin{1};
else
  opt= opt_proplistToStruct(varargin{:});
end
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

pos_screen= get(0, 'ScreenSize');
actualsize(1)= pos_screen(3) - opt.GridSize(2)*2*opt.WindowBorder(1);
actualsize(2)= pos_screen(4) - opt.GridSize(1)*sum(opt.WindowBorder);
actualsize= actualsize - sum(opt.DesktopBorder([1 2; 3 4]));

fig_size= floor(actualsize./fliplr(opt.GridSize));
iv= mod(opt.GridSize(1) - 1 - floor((opt.Fn-1)/opt.GridSize(2)), opt.GridSize(1));
ih= mod(opt.Fn-1, opt.GridSize(2));

incr= fig_size + [2*opt.WindowBorder(1) sum(opt.WindowBorder)];
fig_pos= opt.DesktopBorder([1 2]) + opt.WindowBorder([1 1]) + [ih iv].*incr;

fig_state= struct('WasVisible',[]);
if ishandle(opt.Fn),
  fig_state.WasVisible= get(opt.Fn, 'Visible');
end
if opt.Silent && ishandle(opt.Fn),
  set(0, 'CurrentFigure',opt.Fn);
else
  figure(opt.Fn);
end
if opt.Hide,
  set(opt.Fn, 'Visible','off');
end
if opt.ToolsOff,
  set(gcf, 'ToolBar','none', 'MenuBar','none');
end
fig_size_orig= fig_size;
fig_size= round(fig_size .* opt.Resize);
if opt.ShiftUpwards && fig_size(2)~=fig_size_orig(2),
  fig_pos(2)= fig_pos(2) + fig_size_orig(2) - fig_size(2);
end
drawnow;
set(opt.Fn, 'Position', [fig_pos fig_size]);
if opt.Square,
  oldUnits= get(gcf, 'Units');
  set(gcf, 'Units','Points');
  pos= get(gcf, 'Position');
  pos([3 4])= min(pos([3 4]));
  set(gcf, 'Position',pos);
  set(gcf, 'Units',oldUnits);
end

if ~isdefault.Name,
  set(opt.Fn, 'Name',opt.Name);
end
if ~isempty(opt.Props),
  set(opt.Fn, opt.Props{:});
end
if opt.Clf,
  clf;
end

if nargout==0,
  clear fig_state;
end

author: bbci


text:
nj = 10000;
ms = 210;
mt = 210;

xj = sort((rand(nj,1)*2-1)*pi);
yj = sort((rand(nj,1)*2-1)*pi);
cj = randn(nj,1)+1i*randn(nj,1);

eps=1e-12;

iflag = +1;

tic
fk = dirft2d1(nj,xj,yj,cj,iflag,ms,mt);
toc

tic
fk1 = nufft2d1(nj,xj,yj,cj,iflag,eps,ms,mt);
toc

abs_error=norm(fk-fk1,2)
rel_error=norm(fk-fk1,2)/norm(fk,2)

iflag = -1;

tic
fk = dirft2d1(nj,xj,yj,cj,iflag,ms,mt);
toc

tic
fk1 = nufft2d1(nj,xj,yj,cj,iflag,eps,ms,mt);
toc

abs_error=norm(fk-fk1,2)
rel_error=norm(fk-fk1,2)/norm(fk,2)

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


%==============================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% function [Sc,dS,d2S] = curvatureST(uc,omega,m,varargin)
%
% Matrix-free spatio-temporal curvature regularization energy for vc
% where vc is cell-centered
%
% Sv) = 0.5 * \int_{\omega}\int_0^1
%               alpha(1)*v(x,t)'*A*v(x,t)+ alpha(2)*v(x,t)'*B*v(x,t) dx dt,
%
% where A is the curvature operator and B the first-order time derivative
% operator.
%
% Input:
%
%   vc          instationary velocity field (cell-centered)
%   omega       spatial domain
%   m           number of discretization points in space
%   varargin    optional parameters (see below)
%
% Optional Input:
%   tspan       time span (default: [0 1])
%   nt          number of time points for velocity
%
%
% Output:
%
%   Sc          current value  (0.5 * hd * uc'* A *uc)
%   dS          derivative     (hd * uc'*A )
%   d2S         Hessian        A
%  if ~matrixFree,  d2S is sparse matrix; else, d2S is struct; endif
%
% see also curvature.m
% ==================================================================================

function [Sc,dS,d2S] = mbCurvatureST(vc,omega,m,varargin)
if nargin == 0
    help(mfilename);
    return;
end

if strcmp(vc,'para')
    Sc = 'cell-centered';       % grid
    dS = 0;                     % matrixFree
    d2S = 'pcg';  % solver
    return;
end

persistent A omegaOld mOld alphaOld ntOld tspanOld

if ~exist('mOld','var'),     mOld = [];     end;
if ~exist('omegaOld','var'), omegaOld = []; end;
if ~exist('alphaOld','var'), alphaOld = []; end;

alpha       = [1 1e-3];
tspan       = [0 1];
nt          = [];
for k=1:2:length(varargin), % overwrites default parameter
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

dim = numel(omega)/2;

if isempty(nt), % roughly estimate nt
    nt = round(numel(vc)/(prod(m)*dim))-1;
end

    build =  isempty(mOld) || isempty(omegaOld) ...
        || length(mOld) ~= length(m) || length(omegaOld) ~= length(omega) ...
        || any(mOld ~= m) || any(omegaOld~=omega)|| any(alphaOld~=alpha) || any(tspanOld~=tspan) || any(ntOld~=nt);
    if build,
        fprintf('%s - rebuilding op\n',mfilename);
        A = getCurvatureMatrixST(omega,tspan,m,nt,alpha);
        A = A'*A;
        mOld = m; omegaOld = omega; alphaOld = alpha; ntOld = nt; tspanOld = tspan;
    end;
    dS  = vc'*A;
    Sc  = 0.5*dS*vc;
    d2S = A;

function C = getCurvatureMatrixST(omega,tspan,m,nt,alpha)

h   = (omega(2:2:end)-omega(1:2:end))./m;
hd  = prod(h);
% compute time-stepsize
dt = abs(tspan(2)-tspan(1))/nt;

% build gradient matrix for one transformation
C =  getCurvatureMatrix(omega,m);
a = sqrt(alpha(1).*hd.*dt*[1/2;ones(nt-1,1);1/2]);
% apply spatial regularization to all transformations and sum up in
% time
C =  kron(sdiag(a(:)),C);

% get time regularization matrix
if alpha(2)>0
    b = sqrt(alpha(2)*hd.*dt);
    Bt = b * getTimeDiffusionMatrix(nt,dt,prod(m),length(omega)/2);
    C  = [C;Bt];
end

function D = sdiag(v)
D = diag(sparse(v(:)));

function A = getTimeDiffusionMatrix(nt,dt,n,dim)
Dt=spdiags(ones(nt+1,1)*[-1/dt 1/dt],0:1,nt,nt+1);
A = kron(Dt, speye(n*dim)); % 2nd deriveative over time for two components


text:
function [mclab, mcprc, allprc, trainw] = segmentation2labels(imsegs, smap)

nseg = max(smap);

mclab = zeros(nseg, 1);
mcprc = zeros(nseg, 1);
allprc = zeros(nseg, 7);
trainw = zeros(nseg, 1);

imsegs.npixels = imsegs.npixels(:);
imsegs.labels = imsegs.labels(:);

for s = 1:nseg
    ind = find(smap==s);
    for k = 1:7
        allprc(s, k) = allprc(s, k) + sum((imsegs.labels(ind)==k).*imsegs.npixels(ind));
    end
    allprc(s, :) = allprc(s, :) / max(sum(allprc(s, :)),1E-6);
    [mcprc(s), mclab(s)] = max(allprc(s, :));
    trainw(s) = sum(imsegs.npixels(ind))/sum(imsegs.npixels);
end

author: zouchuhang

repo: LayoutNet


text:
function movie_of_slice_timeseries(imgs, slicenumber, moviename, orientation)
% :Inputs:
%
%   **imgs:**
%        Image names
%
%   **slicenumber:**
%        which slice to view (in volume)
%
%   **moviename:**
%        something like: 'slice10_timeseries.avi'
%
%   **orientation:**
%        'axial' or 'sagittal'
%
% :Examples:
% ::
%
%    slicenumber = 10;
%    moviename = 'slice10_timeseries.avi';

    if(~exist('orientation', 'var') || isempty(orientation))
        orientation = 'axial';
    end

    disp('Mapping volumes.')
    V = spm_vol(imgs);

    disp('Loading data.')

    sl = timeseries_extract_slice(V, slicenumber, orientation);

    disp('Setting limits.')
    mymad = mad(sl(:));
    mymean = mean(sl(:));
    limits = [0 mymean + 3 * mymad];

    slice_fig = figure();

    nimgs = size(sl, 3);

    disp('Starting movie.')
    movlength = 60; % in seconds

    [mov, nframes, axh] = setup_movie([], movlength, moviename);

    for i = 1:nimgs
        H = gca;
        imagesc(sl(:, :, i));
        set(H, 'CLim', limits)
        colorbar
        axis image

        [dd, ff, ee] = fileparts(V(i).fname);
        title(sprintf('Image %3.0f: %s', i, ff), 'FontSize', 24)

        mov = add_a_frame(mov, H);
    end

    tmp = close(mov); %#ok;
    close(slice_fig);

    disp('Success!')
end

function [mov, nframes, axh] = setup_movie(mov, movlength, moviename)
    axh = gca;

    fps = 5;
    nframes = movlength .* fps;

    if isempty(mov)
        mov = avifile(moviename, 'Quality', 75, 'Compression', 'None', 'Fps', fps);
    end

    % add to existing
    %O = struct('add2movie', [], 'zoom', 1, 'azOffset', [], 'elOffset', [], 'timecourse', [], 'fps', 5, 'length', 6);

    axis image
end

function mov = add_a_frame(mov, H)
    drawnow

    lightRestoreSingle(H);
    try
        mov = addframe(mov, H);
    catch
        disp('Cannot write frame.  Failed to set stream format??')
        mov = close(mov);
    end
end


text:
function [bler, ber] = Simulation(max_iter, max_err, max_runs, resolution, ebno_vec, N, K)
R = K/N;
n = log2(N);
num_block_err_bp = zeros(length(ebno_vec), 1);
num_bit_err_bp = zeros(length(ebno_vec), 1);
num_runs = zeros(length(ebno_vec), 1);

%Indices for enc/decoding
[M_up, M_down] = index_Matrix(N);
lambda_offset = 2.^(0 : n);
llr_layer_vec = get_llr_layer(N);

%code construction
design_snr = 2.5;
sigma_cc = 1/sqrt(2 * R) * 10^(-design_snr/20);
[channels, ~] = GA(sigma_cc, N);
[~, channel_ordered] = sort(channels, 'descend');
info_bits = sort(channel_ordered(1 : K), 'ascend');
frozen_bits = ones(N , 1);
frozen_bits(info_bits) = 0;

tic
for i_run = 1 : max_runs
    if mod(i_run, ceil(max_runs/resolution)) == 1
        disp(['Sim iteration running = ', num2str(i_run)]);
        disp(['N = ' num2str(N) ' K = ' num2str(K) ' GA construction SNR = ' num2str(design_snr) 'dB' ' Max Iter Number = ' num2str(max_iter)])
        disp('BP BLER')
        disp(num2str([ebno_vec' num_block_err_bp./num_runs]));
        disp(' ')
    end
    u = zeros(N, 1);
    info = rand(K, 1) < 0.5 ;
    u(info_bits) = info;
    x = polar_encoder(u, lambda_offset, llr_layer_vec);
    bpsk = 1 - 2 * x;
    noise = randn(N, 1);
    for i_ebno = 1 : length(ebno_vec)
        if num_block_err_bp(i_ebno) > max_err
            continue;
        end
        num_runs(i_ebno) = num_runs(i_ebno) + 1;
        sigma = 1 / sqrt(2 * R) * 10^(-ebno_vec(i_ebno)/20);
        y = bpsk + sigma * noise;
        llr = 2/sigma^2 * y;
        [info_esti_bp, ~, ~, ~] = BP_Decoder_LLR(info_bits, frozen_bits, llr, max_iter, M_up, M_down);
        if any(info_esti_bp ~= info)
            num_block_err_bp(i_ebno) =  num_block_err_bp(i_ebno) + 1;
            num_bit_err_bp(i_ebno) = num_bit_err_bp(i_ebno) + sum(info ~= info_esti_bp);
        end

    end
end
toc
bler = num_block_err_bp./num_runs;
ber = num_bit_err_bp./num_runs/K;
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
function [ objectIndex ] = GetScenario_waypoint_90Degrees(varargin)
% This scenario is designed to present a waypoint following exercise.

fprintf('[SCENARIO]\tGetting the waypoint following exercise.\n');

% DEFAULT CONFIGURATION
defaultConfig = struct('file','scenario.mat',...
                       'agents',[],...
                       'agentVelocity',0,...
                       'noiseFactor',0,...
                       'plot',0);

% Instanciate the scenario builder
SBinstance = scenarioBuilder();
% PARSE THE USER OVERRIDES USING THE SCENARIO BUILDER
[inputConfig] = SBinstance.configurationParser(defaultConfig,varargin);
% AGENT CONDITIONS
agentIndex = inputConfig.agents;
agentNumber = numel(agentIndex);                    % Declare the number of agents

%% DEFINE THE AGENT CONFIGURATION
% MOVE THROUGH THE AGENTS AND INITIALISE WITH GLOBAL PROPERTIES
fprintf('[SCENARIO]\tAssigning agent definition...\n');
for i=1:agentNumber
    agentIndex{i}.SetGLOBAL('globalPosition',[0;0;0] + inputConfig.noiseFactor*randn(3,1));
    agentIndex{i}.SetGLOBAL('globalVelocity',[inputConfig.agentVelocity;0;0] + inputConfig.noiseFactor*randn(3,1));
    agentIndex{i}.SetGLOBAL('quaternion',[1;0;0;0]);
end

%% DEFINE THE WAYPOINT CONFIGURATION
fprintf('[SCENARIO]\tBuilding the new scenario...\n');
waypointRadius = 0.1;
nameString = sprintf('WP-%s',agentIndex{1}.name);
waypointPositions = [ 20 20;
                       0 20;
                       0 20];

% ASSIGN AGENT GLOBAL PROPERTIES, ONE SIDE OF THE RINGS TO THE OTHER
waypointIndex{1} = waypoint('radius',waypointRadius,'name',nameString);
waypointIndex{1}.SetGLOBAL('globalPosition',waypointPositions(:,1));
waypointIndex{1}.SetGLOBAL('globalVelocity',zeros(3,1));
waypointIndex{1}.SetGLOBAL('quaternion',[1;0;0;0]);                        % Append properties from the sphereical scenario

waypointIndex{2} = waypoint('radius',waypointRadius,'name',nameString);
waypointIndex{2}.SetGLOBAL('globalPosition',waypointPositions(:,2));
waypointIndex{2}.SetGLOBAL('globalVelocity',zeros(3,1));
waypointIndex{2}.SetGLOBAL('quaternion',[1;0;0;0]);                        % Append properties from the sphereical scenario

for index = 1:numel(waypointIndex)
   waypointIndex{index} = waypointIndex{index}.CreateAgentAssociation(agentIndex{1},1/index);  % Create waypoint with association to agent
end

% BUILD THE COLLECTIVE OBJECT INDEX
objectIndex = horzcat(agentIndex,waypointIndex);
% PLOT THE SCENE
if inputConfig.plot
    SBinstance.plotObjectIndex(objectIndex);
end
% CLEAR THE REMAINING VARIABLES
clearvars -except objectIndex
end

repo: OpenMAS


text:
function filename = dccnpath(filename)

% DCCNPATH manages the filename and path for test files. It helps to locate and read
% test file from Linux, Windows or macOS computers both inside and outside the DCCN.
%
% Use as
%  filename = dccnpath(filename)
%
% The default location for FieldTrip and its test data is '/home/common/matlab/fieldtrip'.
% This function will search-and-replace this string by the location that applies to
% your computer. It will replace '/home' by 'H:' and will replace forward by backward slashes.
%
% In case you have a local copy, you can override the default location by
%   global ft_default
%   ft_default.dccnpath = '/your/copy';
%
% Note that most test scripts expect data located at /home/common/matlab/fieldtrip/data/ftp
% or /home/common/matlab/fieldtrip/data/test, hence you should organize your local
% copy of the data under /your/copy/data/ftp and /your/copy/data/test.
%
% Copyright (C) 2012-2022, Donders Centre for Cognitive Neuroimaging, Nijmegen, NL

% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

global ft_default

% switch between linux and windows paths
if ispc
  filename = strrep(filename,'/home','H:');
  filename = strrep(filename,'/','\');
else
  filename = strrep(filename,'H:','/home');
  filename = strrep(filename,'\','/');
end

% alternative1 alllows to test with local files in the present working directory
% this is often convenient when initially setting up a new test script while the data is not yet uploaded
[p, f, x] = fileparts(filename);
alternative1 = [f x];

if strcmp(alternative1, 'test')
  % this should not be used when the filename is only "test", since that is too generic
  alternative1 = '';
elseif endsWith(pwd, alternative1)
  % exist(alternative1, 'file') also returns 7 in case the present working directory name matches alternative1
  % this is not the match we want
  alternative1 = '';
end

% alternative2 allows for the user to specify the path in the global ft_default variable
% see https://github.com/fieldtrip/fieldtrip/issues/1998
if isfield(ft_default, 'dccnpath')
  alternative2 = strrep(filename, '/home/common/matlab/fieldtrip', ft_default.dccnpath);
else
  alternative2 = '';
end

if ~exist(filename, 'file') && exist(alternative1, 'file')
  warning('using local copy %s instead of %s', alternative1, filename);
  filename = alternative1;
elseif ~exist(filename, 'file') && exist(alternative2, 'file')
  warning('using local copy %s instead of %s', alternative2, filename);
  filename = alternative2;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function cS = forsterite
% Olivine

cs_Fo = crystalSymmetry('mmm',[0.466 1 0.586]);
N = Miller({0,1,0},{0,0,1},{0,2,1},{1,1,0},{1,0,1},{1,2,0},cs_Fo);
dist = [0.4, 1.3, 1.4, 1.05, 1.85, 1.35];
cS = crystalShape(N./dist);

author: mtex-toolbox

repo: mtex


text:
function sF = approximation(nodes, y, varargin)
% computes a least square problem to get an approximation
% Syntax
%   sF = S2FunHarmonic.approximation(S2Grid, f)
%   sF = S2FunHarmonic.approximation(S2Grid, f, 'bandwidth', bandwidth, 'tol', TOL, 'maxit', MAXIT, 'weights', W)
%
% Input
%  S2Grid - grid on the sphere
%  f      - function values on the grid (may be multidimensional)
%
% Options
%  bandwidth  - maximum degree of the spherical harmonics used to approximate the function
%  to         - tolerance for lsqm
%  maxIt      - maximum number of iterations for lsqm
%  W          - weight w_n for the node nodes (default: voronoi weights)
%

% make points unique
s = size(y);
y = reshape(y, length(nodes), []);
[nodes,~,ind] = unique(nodes(:));

% take the mean over duplicated nodes
for k = 1:size(y,2)
  yy(:,k) = accumarray(ind,y(:,k),[],@nanmean); %#ok<AGROW>
end

% consider the case of symmetry
if isa(nodes,'Miller')

  nodes = nodes.symmetrise('noAntipodal').';
  nodes.antipodal = nodes.CS.isLaue;
  yy = repmat(yy,size(nodes,2),1);
  nodes = nodes(:);

end

y = reshape(yy, [length(nodes) s(2:end)]);

tol = get_option(varargin, 'tol', 1e-6);
maxit = get_option(varargin, 'maxit', 40);

if check_option(varargin, 'antipodal') || nodes.antipodal
  nodes.antipodal = true;
  bw = get_option(varargin, 'bandwidth', ceil(sqrt(length(nodes))));
  bw = floor(bw/2)*2; % make bandwidth even
  mask = sparse((bw+1)^2); % only use even polynomial degree
  for m = 0:2:bw
    mask((m^2+1):(m^2+2*m+1), (m^2+1):(m^2+2*m+1)) = speye(2*m+1);
  end
else
  bw = get_option(varargin, 'bandwidth', ceil(sqrt(length(nodes)/2)));
  mask = speye((bw+1)^2);
end

W = get_option(varargin, 'weights');
if isempty(W)
  W = nodes.calcVoronoiArea;
else
  W = W(ind);
end
W = sqrt(W(:));

% initialize nfsft
nfsftmex('precompute', bw, 1000, 1, 0);
plan = nfsftmex('init_advanced', bw, length(nodes), 1);
nfsftmex('set_x', plan, [nodes.rho'; nodes.theta']); % set vertices
nfsftmex('precompute_x', plan);

b = W.*y;

s = size(b);
b = reshape(b, s(1), []);
num = size(b, 2);

fhat = zeros((bw+1)^2, num);
for index = 1:num
  [fhat(:, index), flag] = lsqr(...
    @(x, transp_flag) afun(transp_flag, x, plan, W, mask), ...
    b(:, index), tol, maxit);
  fhat(:, index) = mask*fhat(:, index);
end
fhat = reshape(fhat, [(bw+1)^2 s(2:end)]);

% finalize nfsft
nfsftmex('finalize', plan);

sF = S2FunHarmonic(fhat);

% ensure symmetry if required
if isa(nodes,'Miller'), sF = sF.symmetrise(nodes.CS); end

end

function y = afun(transp_flag, x, plan, W, mask)

if strcmp(transp_flag, 'transp')

  x = x.*W;

  nfsftmex('set_f', plan, x);
  nfsftmex('adjoint', plan);
  y = nfsftmex('get_f_hat_linear', plan);

  y = mask*y;

elseif strcmp(transp_flag, 'notransp')

  x = mask*x;

  nfsftmex('set_f_hat_linear', plan, x);
  nfsftmex('trafo', plan);
  y = nfsftmex('get_f', plan);

  y = y.*W;

end

end


function bool = isOctave(varargin)
%ISOCTAVE Determines whether (a certain) version of Octave is being used
%
% See also: isEnvironment, isMATLAB
bool = isEnvironment('Octave', varargin{:});


text:
function push_cIX_gIX(hfig,cIX,gIX,numK)
% push cIX, gIX to hfig
if ~exist('numK','var')
    numK=max(gIX);
end

var_script={hfig,cIX,gIX,numK};
GUI_FishExplorer('_place_holder_',0,0,0,1,'push_cIX_gIX',var_script);

end


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcLoadSeg(Temp,Event,handles)
handles = guidata(handles.MU_matrix_display);

[filename,pathname]=uigetfile({'*.mat','MAT-files (*.mat)'},'MultiSelect','off');
if filename==0
    return;
end
load([pathname filename]);
try
    if ~isequal(size(handles.Mask),size(Mask))
        error('Segmentation mask has different matrix size.');
    end
    handles.V.Segs=Segs;
    handles.Mask=Mask;
    handles=MU_update_image(handles.Matrix_display_axes,{handles.TMatrix,handles.Mask},handles,0);
catch me
    error_msg{1,1}='ERROR!!! Loading segmentation fails.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);
    return;
end
guidata(handles.MU_matrix_display, handles);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function s = sct(V,varargin)
  % SCT
  %
  % s = sct(V,varargin)
  %
  % Shorthand for:
  %  s = scatter(V(:,1),V(:,2), ...

  switch size(V,2)
  case 3
    ss = scatter3(V(:,1),V(:,2),V(:,3),varargin{:});
  case 2
    ss = scatter(V(:,1),V(:,2),varargin{:});
  end

  if nargout>0, s = ss; end
end

author: alecjacobson

repo: gptoolbox


text:
function [centroids, idx] = runkMeans(X, initial_centroids, ...
                                      max_iters, plot_progress)
%RUNKMEANS runs the K-Means algorithm on data matrix X, where each row of X
%is a single example
%   [centroids, idx] = RUNKMEANS(X, initial_centroids, max_iters, ...
%   plot_progress) runs the K-Means algorithm on data matrix X, where each
%   row of X is a single example. It uses initial_centroids used as the
%   initial centroids. max_iters specifies the total number of interactions
%   of K-Means to execute. plot_progress is a true/false flag that
%   indicates if the function should also plot its progress as the
%   learning happens. This is set to false by default. runkMeans returns
%   centroids, a Kxn matrix of the computed centroids and idx, a m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set default value for plot progress
if ~exist('plot_progress', 'var') || isempty(plot_progress)
    plot_progress = false;
end

% Plot the data if we are plotting progress
if plot_progress
    figure;
    hold on;
end

% Initialize values
[m n] = size(X);
K = size(initial_centroids, 1);
centroids = initial_centroids;
previous_centroids = centroids;
idx = zeros(m, 1);

% Run K-Means
for i=1:max_iters

    % Output progress
    fprintf('K-Means iteration %d/%d...\n', i, max_iters);
    if exist('OCTAVE_VERSION')
        fflush(stdout);
    end

    % For each example in X, assign it to the closest centroid
    idx = findClosestCentroids(X, centroids);

    % Optionally, plot progress here
    if plot_progress
        plotProgresskMeans(X, centroids, previous_centroids, idx, K, i);
        previous_centroids = centroids;
        fprintf('Press enter to continue.\n');
        pause;
    end

    % Given the memberships, compute new centroids
    centroids = computeCentroids(X, idx, K);
end

% Hold off if we are plotting progress
if plot_progress
    hold off;
end

end

author: Ayatans

repo: Machine-Learning-homework


text:
function [mopt, vopt, avals, mvals, vvals] = ...
    scimat_optimal_intersecting_plane(scimat, m0, v0, params)
% SCIMAT_OPTIMAL_INTERSECTING_PLANE  Optimise intersection plane for SCIMAT
% segmentation mask.
%
% [MOPT, VOPT, AVALS, MVALS, VVALS] = ...
%      scimat_optimal_intersecting_plane(SCIMAT, M0, V0, PARAMS)
%
%   This function computes the plane that intersects a SCIMAT segmentation
%   mask in a way that minimizes the segmentation area intersected by the
%   plane. That is, in some sense in finds the plane more orthogonal to the
%   segmented volume.
%
%   (Note that the area is computed on the convex hull of the plane
%   intersection with the volume.)
%
%   SCIMAT is the struct with the volume that we want to intersect (see
%   "help scimat" for details).
%
%   M0 is the rotation centroid. This centroid will not change.
%
%   V0 is a 3-vector that describes the normal vector to the initial
%   intersecting plane. By default, the initial plane is horizontal.
%
%   PARAMS is a struct with optimisation parameters:
%
%   * PARAMS.TYPE is a string. 'local' (default) means that the
%   intersected area will be minimised using multidimensional unconstrained
%   nonlinear minimization (Nelder-Mead, fminsearch). 'global' means that
%   area values will be systematically computed over a range of plane
%   inclination values.
%
%   * PARAMS.RAD can be a scalar or 2-vector:
%
%     scalar:   RAD is the radius of a 2D smoothing disk. The 2D image
%               obtained from intersecting the 3D volume by the plane at
%               each optimization iteration will be smoothed out using the
%               disk before computing the area.
%     2-vector: The smoothing element is a 3D ball. Instead of smoothing
%               a 2D image at each iteration, the whole 3D image volume is
%               smoothed with the ball at the beginning. RAD(1) is the ball
%               radius in the XY-plane. RAD(2) is the ball height.
%
%   * PARAMS.RANGE (global optimisation only): Angular range. The azimuth
%     of V0 will be changed RANGE(1) rad in any direction to compute
%     area values. The elevation of V0 will be changed RANGE(2) rad.
%     (Default RANGE(1) = RANGE(2) = 0.5236 rad = 30º).
%
%   * PARAMS.N (global optimisation only): Number of azimuth (N(1)) or
%   elevation (N(2)) samples. (Default N(1) = N(2) = 61).
%
%   MOPT is the centroid of the optimal intersection.
%
%   VOPT is the normal vector to the optimal plane.
%
%   Note that because M0 and MOPT are both contained in the optimal plane,
%   the duples (M0, VOPT) and (MOPT, VOPT) define the same optimal plane.
%
%   AVALS is a vector with the record of area values from the optimization.
%
%   MVALS is a volume with the coordinates of M at the different tested
%   values. If optimisation is 'local', then MVALS is a matrix where each
%   column is the centroid at each iteration step. If optimisation is
%   'global', then MVALS is a volume where MVALS(T,P,:) is the centroid for
%   the normal vector VVALS(T,P,:).
%
%   VVALS is a volume like MVALS, only for the normal vectors.

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

%% Checks and initialization

% check arguments
narginchk(2, 4) ;
nargoutchk(0, 5) ;

% defaults
if (nargin < 3 || isempty(v0))
    v0 = [0 0 1];
end
if (nargin < 4 || isempty(params))
    params.type = 'local';
    params.rad = [];
    params.range = [30 30] / 180 * pi;
    params.n = [61 61];
    se = [];
end
if (~isfield(params, 'type'))
    params.type = 'local';
end
if (~isfield(params, 'rad'))
    params.rad = [];
    se = [];
end
if (~isfield(params, 'range'))
    params.range = [30 30] / 180 * pi;
end
if (~isfield(params, 'n'))
    params.n = [61 61];
end

% prevent user entering rotation matrix instead of initial vector by
% mistake
if (size(v0, 2) ~= 1 || size(v0, 1) ~= 3)
    error('V0 must be a column 3-vector')
end

% remove the dummy dimension and convert image data to double
scimat = scimat_squeeze(scimat, true);

% convert radius from real world size into number of pixels
params.rad = round(...
    params.rad ./ [scimat.axis(1:length(params.rad)).spacing]);

% create disk for dilation/erosion if we are going to smooth in 2D
if (length(params.rad) == 1)
    se = strel('disk', params.rad);
end

% 3D smoothing of the segmentation edges
if (length(params.rad) == 2)
    se = strel('ball', params.rad(1), params.rad(2));
    scimat.data = imdilate(scimat.data, se);
    scimat.data = imerode(scimat.data, se);
end

% generate 3D grid of coordinates
[x, y, z] = scimat_ndgrid(scimat);

% % DEBUG: compute intersection of SCIMAT volume with the initial plane
% % (if you want to visualize the image as in Seg3D, you need to do 'axis
% % xy')
% im = scimat_intersect_plane(scimat, m0, v0, x, y, z);

%% Optimisation of the intersection area

% convert Cartesian coordinates into spherical coordinates (length has to
% be one); note: we use phi for azimuth, and theta for elevation, contrary
% to Matlab's naming convention
[phi0, theta0] = cart2sph(v0(1), v0(2), v0(3));

% group spherical coordinates into vector
alpha0 = [phi0, theta0];

% init variables to keep track of the evolution of area values in the
% optimisation
avals = [];
mvals = [];
vvals = [];

% local or global optimisation
if strcmp(params.type, 'local')

    % run optimisation to find minimum area; note that v0 is the only
    % variable optimised, but the rest (scimat, x, y, z, m, rad, se) are
    % available to segmented_area_of_intersection() because the latter is a
    % subfunction
    alpha = fminsearch(@segmented_area_of_intersection, alpha0);

    % convert result from spherical to Carterian coordinates
    [aux1 aux2 aux3] = sph2cart(alpha(1), alpha(2), 1.0);
    vopt = [aux1 aux2 aux3];

    % final centroid of the intersecting plane
    mopt = mvals(:, end);

elseif strcmp(params.type, 'global')

    % interval of azimuth angle values, phi \in [-180º, 180º] or \in
    % [0, 360º]
    phimin = phi0 - abs(params.range(1));
    phimax = phi0 + abs(params.range(1));

    % interval of elevation angle values, theta \in [-90º, 90º]
    thmin = max(theta0 - abs(params.range(2)), -pi/2);
    thmax = min(theta0 + abs(params.range(2)), pi/2);

    % sample angle intervals
    phivals = linspace(phimin, phimax, params.n(1));
    thetavals = linspace(thmin, thmax, params.n(2));

    % create matrices to save outputs; note that for each area value we
    % need to save a 3-vector with the rotation point, and a 3-vector with
    % the normal plane
    avals = zeros(length(thetavals), length(phivals));
    mvals = zeros(length(thetavals), length(phivals), 3);
    vvals = zeros(length(thetavals), length(phivals), 3);

    % compute area for each combination of elevation and azimuth angles
    for T = 1:length(thetavals) % elevation
        for P = 1:length(phivals) % azimuth
            [a, mnew, v] = segmented_area_of_intersection(...
                [phivals(P) thetavals(T)]);

            % put values in output matrices
            avals(T, P) = a;
            mvals(T, P, :) = mnew;
            vvals(T, P, :) = v;

        end
    end

    % find minimum area
    [foo, idx] = min(avals(:));

    % convert linear index to multiple subscripts
    [T, P] = ind2sub(size(avals), idx);

    % output optimal plane
    mopt = squeeze(mvals(T, P, :));
    vopt = squeeze(vvals(T, P, :));

else
    error(['Optimisation type not implemented: ' params.type])
end

    %% Objective function (the function we are trying to minimise)

    % rotate plane, intersect with image, and compute segmented area
    function [a, mnew, v] = segmented_area_of_intersection(alpha)

        % convert spherical to Carterian coordinates
        [aux1 aux2 aux3] = sph2cart(alpha(1), alpha(2), 1.0);
        v = [aux1 aux2 aux3]';

        % vector cannot be zero
        if (norm(v) == 0)
            error('Normal vector to plane cannot be (0,0,0)')
        end

        % this function cannot deal with vertical planes, because of a
        % singularity
        if (v(3) == 0)
            error('Intersecting plane cannot be vertical')
        end

        % compute intersection of plane with volume
        [im, zp, xp, yp] = scimat_intersect_plane(scimat, m0, v, x, y, z);

        % 2D smoothing of the segmentation edges
        if (length(params.rad) == 1)
            im = imdilate(im, se);
            im = imerode(im, se);
        end

%         % DEBUG: plot rotated plane
%         hold off
%         plot3(xp(:), yp(:), zp(:), '.r')

        % find segmented voxels in the 2D cut
        idx = find(im);

        % get coordinates of segmented voxels
        xps = xp(idx);
        yps = yp(idx);
        zps = zp(idx);

%         % DEBUG: visualize intersection projected onto horizontal plane
%         hold off
%         imagesc(xp(:), yp(:), im > 0)
%         hold on
%         % DEBUG: compute and plot convex hull
%         idx2 = convhull(xps, yps);
%         vxs = xps(idx2);
%         vys = yps(idx2);
%         plot(vxs, vys, 'w')
%         xlabel('x (m)')
%         ylabel('y (m)')
%         pause

        % compute a rotation matrix from the Cartesian system to the
        % rotated plane
        rotmat = vec2rotmat(v);

        % we are now seeing the rotated plane projected onto the horizontal
        % plane, i.e. we see the segmentation mask in perspective.
        % In order to see the true area of the segmentation mask, we need
        % to change the system of coordinates so that the rotated plane
        % becames the XY plane

        % first, move segmented points so that centroid is at (0,0,0)...
        xps = xps - m0(1);
        yps = yps - m0(2);
        zps = zps - m0(3);

        % ...second, make the rotated plane horizontal, by inverting the
        % rotation...
        xyzps = [xps(:) yps(:) zps(:)] * rotmat;
        xps = xyzps(:, 1);
        yps = xyzps(:, 2);
        zps = xyzps(:, 3);

        % if everything has gone alright, then the z-coordinate of xyzp
        % should be zero (+numerical errors), because the rotated plane is
        % now the XY plane
        assert(abs(min(zps)) < 1e-10)
        assert(abs(max(zps)) < 1e-10)

%         % DEBUG: visualize segmentation mask in real world coordinates
%         hold off
%         plot(xps + m0(1), yps + m0(2), 'r*')
%         axis ij

        % compute convex hull (reuse idx2): note convex hull coordinates
        % are on projected space
        idx2 = convhull(xps, yps);
        vxs = xps(idx2);
        vys = yps(idx2);

        % compute x-,y-coordinates centroid and area of polygon
        [mnew, a] = polycenter(vxs, vys);
        mnew(3) = 0;

        % the centroid is now on projected coordinates, but we need to put
        % it back on the real world coordinates
        mnew = rotmat * mnew';
        mnew = (mnew' + m0)';

        % for the global algorithm, values are recorded in a different way
        if strcmp(params.type, 'local')
            % kept track of optimisation evolution
            avals = [avals a];
            vvals = [vvals v];
            mvals = [mvals mnew];
        end

    end
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function quadmom_test08 ( )

%*****************************************************************************80
%
%% QUADMOM_TEST08 integrates sin(x) against a lower truncated normal weight.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gene Golub, John Welsch,
%    Calculation of Gaussian Quadrature Rules,
%    Mathematics of Computation,
%    Volume 23, Number 106, April 1969, pages 221-230.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUADMOM_TEST08:\n' );
  fprintf ( 1, '  Integrate sin(x) against a lower truncated normal weight.\n' );

  mu = 0.0;
  sigma = 1.0;
  a = -3.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  MU = %g\n', mu );
  fprintf ( 1, '  SIGMA = %g\n', sigma );
  fprintf ( 1, '  A = %g\n', a );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   N    Estimate\n' );
  fprintf ( 1, '\n' );
%
%  N is the order of the quadrature rule.
%
  for n = 1 : 9
%
%  Compute the M = 2 * N + 1 moments.
%
    m = 2 * n + 1;

    moment = moments_truncated_normal_a ( m, mu, sigma, a );
%
%  Compute the N points and weights by the method of moments.
%
    [ x, w ] = moment_method ( n, moment );
%
%  Evaluate the quadrature rule when f(x) = sin(x).
%
    q = w' * sin ( x );

    fprintf ( 1, '  %2d  %14.6g\n', n, q );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = lp_value ( n, o, x )

%*****************************************************************************80
%
%% LP_VALUE evaluates a Legendre polynomial at several points X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, integer O, the degree of the polynomial.
%    0 <= O.
%
%    Input, real X(N,1), the evaluation points.
%
%    Output, real VALUE(N,1), the value of the Legendre polynomial
%    of degree O at the points X.
%
  x = x ( : );

  v = zeros ( n, o + 1 );

  v(1:n,1) = 1.0;

  if ( 1 <= o )

    v(1:n,2) = x(1:n,1);

    for j = 2 : o

      v(1:n,j+1) = ( ( 2 * j - 1 ) * x(1:n,1) .* v(1:n,j)   ...
                  -  (     j - 1 ) *             v(1:n,j-1) ) ...
                  /  (     j     );

    end

  end

  value(1:n,1) = v(1:n,o+1);

  return
end