text:
% Isoparametric Formulation Implementation
% clear memory
clear all
close all
clc
% E: modulus of elasticity
% A: area of cross section
% L: length of bar
E=8; L=4;
u_exact=@(x) (56-8*(x-2)-24*heaviside(x-5))/2/x;
hold on;
ezplot(u_exact,[2 6])
title('Exact solution v.s. FEM','interpreter','latex');
xlabel('x','interpreter','latex');
ylabel('Axial stress, $\it{\sigma}_{x}$','interpreter','latex');
for i=1:4  %For NEL = 1, 2, 4, 8
fprintf( '\nNumber of elements:%d\n\n',2^(i-1) );
% numberElements: number of elements
numberElements=2^(i-1);
% numberNodes: number of nodes
NNOD = 2;
numberNodes=numberElements+1;
% generation of coordinates and connectivities
elementNodes=[1:numberNodes-1;2:numberNodes]';
nodeCoordinates=linspace(2,L+2,numberNodes);
%Generate element length vector
Le=ones(1,numberElements)*L/numberElements;
A=(nodeCoordinates(1:end-1)+Le./2)*2;
%Evaluate area for each cross-section at the center of each element by A=2x;
% for structure:
   % displacements: displacement vector
   % force : force vector
   % stiffness: stiffness matrix
force=zeros(numberNodes,1);
stiffness=zeros(numberNodes,numberNodes);
% computation of the system stiffness matrix and force vector
for e=1:numberElements;
  % elementDof: element degrees of freedom (Dof)
  elementDof=elementNodes(e,:);
  detJacobian=Le(e)/2;
  invJacobian=1/detJacobian;
  ngp = 2;
  xc=0.5*(nodeCoordinates(elementDof(1))+nodeCoordinates(elementDof(end)));
  [w,xi]=gauss1d(ngp);
  if(nodeCoordinates(elementDof(1))==5)
     x=(5-xc)/detJacobian;
     [s,n]=shapeFunctionL2(x);
     force(elementDof)= force(elementDof)+...
         24*s';
  elseif(nodeCoordinates(elementDof(1))<5&&...
          nodeCoordinates(elementDof(end))>5)
     x=(5-xc)/detJacobian;
     [s,n]=shapeFunctionL2(x);
     force(elementDof)= force(elementDof)+...
         24*s';
  end
  for ip=1:ngp;
      [shape,naturalDerivatives]=shapeFunctionL2(xi(ip));
      B=naturalDerivatives*invJacobian;
      stiffness(elementDof,elementDof)=...
      stiffness(elementDof,elementDof)+ B'*B*w(ip)*detJacobian*E*A(e);
      force(elementDof)=force(elementDof)+...
          8*shape'*detJacobian*w(ip);
  end
end

% boundary conditions and solution
% prescribed dofs
prescribedDof=[1];
% solution
GDof=numberNodes;
displacements=solution(GDof,prescribedDof,stiffness,force);
% output displacements/reactions
outputDisplacementsReactionsPretty(displacements,stiffness, ...
    numberNodes,prescribedDof,force)
%stress/strain recover
fprintf('Axial stress\n')
fprintf('element\tgauss point\taxial stress\n')
ngp = 2;
elementNodeCoor=zeros(numberElements*NNOD,1);
elementNodeStr=zeros(numberElements*NNOD,1);
for e=1:numberElements;
  % elementDof: element degrees of freedom (Dof)
  elementDof=elementNodes(e,:) ;
  detJacobian=Le(e)/2;
  invJacobian=1/detJacobian;
  % Nodal coordiantes in natural coordinate
  xi=[-1 1];
  for ip=1:NNOD
      [shape,naturalDerivatives]=shapeFunctionL2(xi(ip));
      B=naturalDerivatives*invJacobian;
      elementNodeStr(ip+(e-1)*NNOD,1)=E*B*displacements(elementDof,1);
      fprintf('%2.0f\t%2.0fth ip\t%10.4e\n', e, ip,...
          elementNodeStr(ip+(e-1)*NNOD,1))
  end
  elementNodeCoor(1+(e-1)*ngp:2+(e-1)*ngp)=[nodeCoordinates(elementDof(1))...
      nodeCoordinates(elementDof(2))];
end

%post process
switch numberElements
    case 1
        str='r*--';
    case 2
        str='g*--';
    case 4
        str='k*--';
    case 8
        str='m*--';
end

plot(elementNodeCoor,elementNodeStr,str)
hold on;
end
legend('Exact solution','NEL=1','NEL=2','NEL=4','NEL=8','interpreter','latex');

repo: Aero-matlab


text:
%DEMO_MODELASSESMENT2  Demonstration for model assessment when the observation
%                      model is non-Gaussian
%
%  Description
%    We will consider the classification problem in demo_classific.
%    The analysis is conducted with full Gaussian process using
%    both probit and logit likelihood. The performance of these two
%    models are compared by evaluating the ten-fold
%    cross-validation, leave-one-out cross-validation, WAIC, DIC
%    and the effective number of parameters The inference will be
%    conducted using maximum a posterior (MAP) estimate for the
%    parameters using EP and Laplace approximation, via full Markov
%    chain Monte Carlo (MCMC) and with an integration approximation
%    (IA) for the parameters.
%
%    This demo is organised in two parts:
%     1) data analysis with with probit likelihood
%     2) data analysis with with logit likelihood
%
%  See also
%    DEMO_CLASSIFIC1, DEMO_MODELASSESMENT1

% Copyright (c) 2009-2010 Jarno Vanhatalo
% Copyright (c) 2010 Aki Vehtari

% 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.

% =====================================
% 1) data analysis with probit likelihood
% =====================================
disp('Data analysis with probit likelihood')
S = which('demo_classific');
L = strrep(S,'demo_classific.m','demodata/synth.tr');
x=load(L);
y=x(:,end);
y = 2.*y-1;
x(:,end)=[];
[n, nin] = size(x);

DIC=repmat(NaN,1,8);DIC2=repmat(NaN,1,8);DIC_latent=repmat(NaN,1,8);
p_eff=repmat(NaN,1,8);p_eff2=repmat(NaN,1,8);p_eff_latent=repmat(NaN,1,8);p_eff_latent2=repmat(NaN,1,8);

% Create covariance functions
gpcf = gpcf_sexp('lengthScale', [0.9 0.9], 'magnSigma2', 2);

% Set the prior for the parameters of covariance functions
pl = prior_logunif();
gpcf = gpcf_sexp(gpcf, 'lengthScale_prior', pl,'magnSigma2_prior', pl); %

% Create the GP structure
gp = gp_set('lik', lik_probit, 'cf', gpcf, 'jitterSigma2', 1e-4);

% ------- Laplace approximation --------
disp(['Probit with Laplace integration over the latent values '; ...
      'and MAP estimate for the parameters                    '])

% Set the approximate inference method
gp = gp_set(gp, 'latent_method', 'Laplace');

n=length(y);

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% Evaluate the effective number of parameters and DIC with focus on
% latent variables.
models{1} = 'pr_Laplace';
p_eff_latent(1) = gp_peff(gp, x, y);
[DIC_latent(1), p_eff_latent2(1)] = gp_dic(gp, x, y, 'focus', 'latent');
WAIC(1) = gp_waic(gp,x,y);

% Evaluate the 10-fold cross-validation results.
cvres = gp_kfcv(gp, x, y, 'display', 'fold');
mlpd_cv(1) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(gp, x, y);
mlpd_loo(1) = mean(lpy);

% ------- Expectation propagation --------
disp(['Probit with EP integration over the latent values and MAP '; ...
      'estimate for the parameters                               '])

% Set the approximate inference method
gp = gp_set(gp, 'latent_method', 'EP');

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the BFGS quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% Evaluate the effective number of parameters and DIC with focus on
% latent variables.
models{2} = 'pr_EP';
p_eff_latent(2) = gp_peff(gp, x, y) ;
[DIC_latent(2), p_eff_latent2(2)] = gp_dic(gp, x, y, 'focus', 'latent');
WAIC(2) = gp_waic(gp,x,y);

% Evaluate the 10-fold cross-validation results.
cvres = gp_kfcv(gp, x, y, 'display', 'fold');
mlpd_cv(2) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(gp, x, y);
mlpd_loo(2) = mean(lpy);

% ------- MCMC ---------------
disp(['Probit with MCMC integration over the latent values and '; ...
      'the parameters                                          '])

% Set the approximate inference method
gp = gp_set(gp, 'latent_method', 'MCMC');

% Sample
mcopt.nsamples=220;mcopt.display=20;
[rgp,gp]=gp_mc(gp, x, y, mcopt);
rgp=thin(rgp, 21, 2);

% Evaluate the effective number of parameters and DIC with focus on
% latent variables.
models{3} = 'pr_MCMC';
[DIC(3), p_eff(3)] =  gp_dic(rgp, x, y, 'focus', 'param');
[DIC2(3), p_eff2(3)] =  gp_dic(rgp, x, y, 'focus', 'all');
WAIC(3) = gp_waic(rgp,x,y);

% Evaluate the 10-fold cross-validation results.
mcopt.nsamples=50;mcopt.display=20;
cvres = gp_kfcv(gp, x, y, 'inf_method', 'MCMC', 'opt', mcopt, 'display', 'fold');
mlpd_cv(3) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(rgp, x, y);
mlpd_loo(3) = mean(lpy);

% --- Integration approximation approach ---
disp(['Probit with EP integration over the latent values and '; ...
      'grid integration over the parameters                  '])

% Use EP
gp = gp_set(gp, 'latent_method', 'EP');

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the BFGS quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% now perform the integration
clear opt
opt.int_method = 'grid';
opt.step_size = 2;
gp_array = gp_ia(gp, x, y, opt);

models{4} = 'pr_IA';
[DIC(4), p_eff(4)] =  gp_dic(gp_array, x, y, 'focus', 'param');
[DIC2(4), p_eff2(4)] =  gp_dic(gp_array, x, y, 'focus', 'all');
WAIC(4) = gp_waic(gp_array,x,y);

% Then the 10 fold cross-validation.
cvres = gp_kfcv(gp, x, y, 'inf_method', 'IA', 'opt', opt, 'display', 'fold');
mlpd_cv(4) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(gp_array, x, y);
mlpd_loo(4) = mean(lpy);

% =====================================
% 2) data analysis with logit likelihood
% =====================================
disp('Data analysis with logit likelihood')

S = which('demo_classific');
L = strrep(S,'demo_classific.m','demodata/synth.tr');
x=load(L);
y=x(:,end);
y = 2.*y-1;
x(:,end)=[];
[n, nin] = size(x);

% Create covariance functions
gpcf = gpcf_sexp('lengthScale', [0.9 0.9], 'magnSigma2', 2);

% Set the prior for the parameters of covariance functions
pl = prior_logunif();
gpcf = gpcf_sexp(gpcf, 'lengthScale_prior', pl,'magnSigma2_prior', pl); %

% Create the likelihood structure
lik = ('init');

% Create the GP structure
gp = gp_set('lik', lik_logit, 'cf', gpcf, 'jitterSigma2', 1e-4);

% ------- Laplace approximation --------
disp(['Logit with Laplace integration over the latent values and '; ...
      'MAP estimate for the parameters                           '])

% Set the approximate inference method
gp = gp_set(gp, 'latent_method', 'Laplace');

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the BFGS quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% Evaluate the effective number of parameters and DIC with focus on
% latent variables.
models{5} = 'lo_Laplace';
p_eff_latent(5) = gp_peff(gp, x, y);
[DIC_latent(5), p_eff_latent2(5)] = gp_dic(gp, x, y, 'focus', 'latent');
WAIC(5) = gp_waic(gp,x,y);

% Evaluate the 10-fold cross-validation results.
cvres = gp_kfcv(gp, x, y, 'display', 'fold');
mlpd_cv(5) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(gp, x, y);
mlpd_loo(5) = mean(lpy);

% ------- Expectation propagation --------
disp(['Logit with EP integration over the latent values and MAP'; ...
      'estimate for the parameters                             '])

% Set the approximate inference method
gp = gp_set(gp, 'latent_method', 'EP');

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the BFGS quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% Evaluate the effective number of parameters and DIC with focus on
% latent variables.
models{6} = 'lo_EP';
p_eff_latent(6) = gp_peff(gp, x, y) ;
[DIC_latent(6), p_eff_latent2(6)] = gp_dic(gp, x, y, 'focus', 'latent');
WAIC(6) = gp_waic(gp,x,y);

% Evaluate the 10-fold cross-validation results.
cvres = gp_kfcv(gp, x, y, 'display', 'fold');
mlpd_cv(6) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(gp, x, y);
mlpd_loo(6) = mean(lpy);

% ------- MCMC ---------------
disp(['Logit with MCMC integration over the latent values and '; ...
      'the parameters                                         '])

% Set the approximate inference method
gp = gp_set(gp, 'latent_method', 'MCMC');

% Sample
mcopt.nsamples=200;mcopt.display=20;
[rgp,gp] = gp_mc(gp, x, y, mcopt);
rgp=thin(rgp, 21, 2);

% Evaluate the effective number of parameters and DIC with focus on latent variables.
models{7} = 'lo_MCMC';
[DIC(7), p_eff(7)] =  gp_dic(rgp, x, y, 'focus', 'param');
[DIC2(7), p_eff2(7)] =  gp_dic(rgp, x, y, 'focus', 'all');
WAIC(7) = gp_waic(rgp,x,y);

% Evaluate the 10-fold cross-validation results.
mcopt.nsamples=50;mcopt.display=20;
cvres = gp_kfcv(gp, x, y, 'inf_method', 'MCMC', 'opt', mcopt, 'display', 'fold');
mlpd_cv(7) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(rgp, x, y);
mlpd_loo(7) = mean(lpy);

% --- Integration approximation approach ---
disp(['Logit with EP integration over the latent values and grid '; ...
      'integration over the parameters                           '])

% Use EP
gp = gp_set(gp, 'latent_method', 'EP');

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the BFGS quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% now perform the integration
clear opt
opt.int_method = 'grid';
opt.step_size = 2;
gp_array = gp_ia(gp, x, y, opt);

models{8} = 'lo_IA';
[DIC(8), p_eff(8)] =  gp_dic(gp_array, x, y, 'focus', 'param');
[DIC2(8), p_eff2(8)] =  gp_dic(gp_array, x, y, 'focus', 'all');
WAIC(8) = gp_waic(gp_array,x,y);

% Then the 10 fold cross-validation.
cvres = gp_kfcv(gp, x, y, 'inf_method', 'IA', 'opt', opt, 'display', 'fold');
mlpd_cv(8) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(gp_array, x, y);
mlpd_loo(8) = mean(lpy);

%========================================================
% PART 4 Print the results
%========================================================
disp('Summary of the results')

S = '       ';
for i = 1:length(models)
    S = [S '  ' models{i}];
end

S = sprintf([S '\n CV-mlpd  %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], mlpd_cv);
S = sprintf([S '\n LOO-mlpd %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], mlpd_loo);
S = sprintf([S '\n ']);
S = sprintf([S '\n WAIC     %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], WAIC);
S = sprintf([S '\n ']);
S = sprintf([S '\n DIC_h    %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], DIC);
S = sprintf([S '\n DIC_a    %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], DIC2);
S = sprintf([S '\n DIC_l    %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], DIC_latent);
S = sprintf([S '\n peff_h   %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], p_eff);
S = sprintf([S '\n peff_a   %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], p_eff2);
S = sprintf([S '\n peff_l   %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], p_eff_latent);
S = sprintf([S '\n peff_l2  %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], p_eff_latent2);
S = sprintf([S '\n ']);
S = sprintf([S '\n The notation is as follows:']);
S = sprintf([S '\n pr_*     = probit likelihood and inference method']);
S = sprintf([S '\n lo_*     = logit likelihood and inference method']);
S = sprintf([S '\n CV-mlpd  = mean log predictive density from the 10-fold CV. ']);
S = sprintf([S '\n LOO-mlpd = mean log predictive density from the 10-fold CV. ']);
S = sprintf([S '\n WAIC     = Widely applicable information criterion. ']);
S = sprintf([S '\n DIC_h    = DIC with focus on parameters. ']);
S = sprintf([S '\n DIC_a    = DIC with focus on parameters and laten variables (all). ']);
S = sprintf([S '\n DIC_l    = DIC with focus on latent variables. ']);
S = sprintf([S '\n peff_h   = effective number of parameters (latent variables marginalized). ']);
S = sprintf([S '\n peff_a   = effective number of parameters and latent variables. ']);
S = sprintf([S '\n peff_l   = effective number of latent variables evaluated with gp_peff. ']);
S = sprintf([S '\n peff_l2  = effective number of latent variables evaluated with gp_dic. ']);
S = sprintf([S '\n '])

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function ObjVal = step2(Chrom,mydata)
Dim=size(Chrom,2);
[Nind,Nvar] = size(Chrom);
Chrom=ceil(Chrom+0.5);
ObjVal = sum((Chrom .* Chrom)')';
return

author: Sable

repo: mcbench-benchmarks


text:
function [ a, b ] = p01_ab ( m )

%*****************************************************************************80
%
%% P01_AB evaluates the limits of the optimization region for problem 01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Output, real A(M), B(M), the lower and upper bounds.
%
  a(1:m) = - 5.0;
  b(1:m) = + 5.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function engine = bk_ff_hmm_inf_engine(bnet)
% BK_FF_HMM_INF_ENGINE Naive (HMM-based) implementation of fully factored form of Boyen-Koller
% engine = bk_ff_hmm_inf_engine(bnet)
%
% This is implemented on top of the forwards-backwards algo for HMMs,
% so it is *less* efficient than exact inference! However, it is good for educational purposes,
% because it illustrates the BK algorithm very clearly.

[persistent_nodes, transient_nodes] = partition_dbn_nodes(bnet.intra, bnet.inter);
assert(isequal(sort(bnet.observed), transient_nodes));
[engine.prior, engine.transmat] = dbn_to_hmm(bnet);

ss = length(bnet.intra);

engine.bel = [];
engine.bel_marginals = [];
engine.marginals = [];

engine = class(engine, 'bk_ff_hmm_inf_engine', inf_engine(bnet));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% Steady state transport solutions for one fixed value of Pe and various values of Da1
x = [0:0.01:1];
Pe = 1; Da1 = 16;
s = sqrt(0.25*Pe*Pe+Pe/Da1);
mu1 = 0.5*Pe+s; mu2 = 0.5*Pe-s;
s = mu2*exp(mu2)-mu1*exp(mu1);
c = (mu2*exp(mu2)*exp(mu1*x)-mu1*exp(mu1)*exp(mu2*x))./s;
for Da1 = [8 4 2 1 0.5 0.25 0.125];
    s = sqrt(0.25*Pe*Pe+Pe/Da1);
    mu1 = 0.5*Pe+s; mu2 = 0.5*Pe-s;
    s = mu2*exp(mu2)-mu1*exp(mu1);
    c = [c;(mu2*exp(mu2)*exp(mu1*x)-mu1*exp(mu1)*exp(mu2*x))./s];
end
plot (x,c);
legend('Da_1=16','Da_1=8','Da_1=4','Da_1=2','Da_1=1','Da_1=0.5','Da_1=0.25','Da_1=0.125');
xlabel ('x/L [-]'); ylabel ('c/c_{in} [-]');


text:
function [gradISE,gradJhr,gradJrr]=GaussMixISELTCovGrads(w1,mu1,P1,w2,mu2,P2,L2,PDFVals12,PDFVals22)
%%GAUSSMIXISELTCOVGRADS Compute the gradient of the non-normalized
%               integrated squared error (ISE) between two Gaussian mixture
%               PDFs with respect to the elements of square-roots of the
%               covariance matrices of the individual components of the
%               second Gaussian mixture. For a covariance matrix P, the
%               square root L is such that P=L*L'. However, L does not need
%               to be lower-triangular. This gradient can arise when
%               optimizing over the ISE with respect to the covariance
%               matrices. Optimizing over the square-root matrices rather
%               than the matrices themselves ensures that the covariance
%               matrices are always positive (semi)definite.
%
%INPUTS: w1 The N1X1 or 1XN1 vector of weights for the first Gaussian
%           mixture. All w1>0 and sum(w1)=1.
%       mu1 The xDimXN1 set of mean vectors for the first Gaussian mixture
%           distribution.
%        P1 The xDimXxDimXN1 set of positive definite covariance matrices
%           for the first Gaussian mixture distirbution.
% w2, mu2, P2, The length N2, xDimXN2, and xDimXxDimXN2 set of weights,
%           mean vectors and positive-definite covariance matrices for the
%           second Gaussian mixture distribution.
%        L2 The xDimXxDimXN2 set of square roots of the covariance matrices
%           in P2 such that P2(:,:,i)=L2(:,:,i)*L2(:,:,i)'.
% PDFVals12 A matrix such that the value in element (i,j) is
%           N(mu1(:,i);mu2(:,j),P1(:,:,i)+P2(:,:,j)), where N indicates the
%           multivariate Gaussian PDF evaluated at the first argument with
%           the second and third arguments being the mean and covarince
%           matrix. This parameter is returned by computeGaussMixISE.
% PDFVals22 A matrix such that the value in element (i,j) is
%           N(mu2(:,i);mu2(:,j),P2(:,:,i)+P2(:,:,j)). This parameter is
%           returned by computeGaussMixISE.
%
%OUTPUTS: gradISE The xDimXxDimXN2 set of derivatives of the ISE with
%                 respect to the elements of L2 (the square roots of P2).
% gradJhr, gradJrr In Chapter 3 of [1], the ISE is expressed in terms of
%                 Jhr and Jrr terms. These are the xDimXxDimXN2 gradients
%                 of those terms.
%
%Formule for gradJhr and gradJrr are Equation 3.45 in Section 3.3.3.3 of
%[1]. They relate to the ISE via Equation 3.20. See the function
%computeGaussMixISE to compute the ISE.
%
%EXAMPLE 1:
%In this example with a scalar PDF, we verify that the gradient obtained
%from this function is consistent with numerical differentiation.
% w1=[0.03,0.18,0.12,0.19,0.02,0.16,0.06,0.1,0.08,0.06];
% n1=length(w1);
% mu1=[1.45,2.20,0.67,0.48,1.49,0.91,1.01,1.42,2.77,0.89];
% P1=[0.0487,0.0305,0.1171,0.0174,0.0295,0.0102, 0.0323, 0.0380, 0.0115, 0.0679];
% P1=reshape(P1,[1,1,n1]);

% %The second PDF is the first with the five least-weight components deleted.
% w2=[0.18,0.12,0.19,0.16,0.1,0.08];
% w2=w2/sum(w2);
% n2=length(w2);
% mu2=[2.20,0.67,0.48,0.91,1.42,2.77];
% P2=[0.0305,0.1171,0.0174,0.0102,0.0380,0.0115];
% P2=reshape(P2,[1,1,n2]);
%
% [ISEVal,PDFVals12,PDFVals22]=computeGaussMixISE(w1,mu1,P1,w2,mu2,P2);
% L2=sqrt(P2);
% epsVal=1e-8;
% gradISENum=zeros(1,1,n2);
% for k=1:n2
%     L2Cur=L2;
%     L2Cur(1,1,k)=L2Cur(1,1,k)+epsVal;
%     ISEValCur=computeGaussMixISE(w1,mu1,P1,w2,mu2,L2Cur.^2);
%     gradISENum(k)=(ISEValCur-ISEVal)/epsVal;
% end
% gradISE=GaussMixISELTCovGrads(w1,mu1,P1,w2,mu2,P2,L2,PDFVals12,PDFVals22);
% RelErr=max(abs((gradISENum(:)-gradISE(:))./gradISENum(:)))
%The relative error will be about 7.162e-6, which indicates good numeric
%agreement.
%
%EXAMPLE 2:
%In this example with a bivariate PDF, we verify that the gradient obtained
%from this function is consistent with numerical differentiation.
% w1=[0.25;0.5;0.25];
% mu1=zeros(2,2);
% mu1(:,1)=[1;-1];
% mu1(:,2)=[-1;1];
% mu1(:,3)=[0;0];
% P1=zeros(2,2,2);
% P1(:,:,1)=[4/9,  14/45;
%            14/45,4/9];
% P1(:,:,2)=[4/9, 0;
%            0, 4/9];
% P1(:,:,3)=[2/9, -1/9;
%           -1/9, 3/9];
%
% %The second distribution just throws out the first component.
% w2=w1(2:3);
% w2=w2/sum(w2);
% mu2=mu1(:,2:3);
% P2=P1(:,:,2:3);
%
% n2=length(w2);
% [ISEVal,PDFVals12,PDFVals22]=computeGaussMixISE(w1,mu1,P1,w2,mu2,P2);
%
% numDim=size(mu2,1);
% L2=zeros(numDim,numDim,n2);
% for j=1:n2
%     L2(:,:,j)=chol(P2(:,:,j),'lower');
% end
%
% epsVal=1e-8;
% gradISENumDiff=zeros(numDim,numDim,n2);
% for j=1:n2
%     for curEl=1:(numDim^2)
%         [i1,i2]=ind2sub([numDim,numDim],curEl);
%
%         LCur=L2(:,:,j);
%         LCur(i1,i2)=LCur(i1,i2)+epsVal;
%
%         PCur=P2;
%         PCur(:,:,j)=LCur*LCur';
%         ISEValCur=computeGaussMixISE(w1,mu1,P1,w2,mu2,PCur);
%         gradISENumDiff(i1,i2,j)=(ISEValCur-ISEVal)/epsVal;
%     end
% end
% gradISE=GaussMixISELTCovGrads(w1,mu1,P1,w2,mu2,P2,L2,PDFVals12,PDFVals22);
% RelErr=max(max(abs((gradISENumDiff-gradISE)./gradISENumDiff)))
%The relative error will be about 6.99e-7, which indicates good numeric
%agreement.
%
%REFERENCES:
%[1] J. L. Williams, "Gaussian mixture reduction for tracking multiple
%    maneuvering targets in clutter," Master's thesis, Air Force Institute
%    of Technology, Mar. 2003. [Online].
%    Available: http://www.dtic.mil/srch/doc?collection=t3&id=ADA415317
%
%May 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

Nh=length(w1);
Nr=length(w2);
xDim=size(mu1,1);

%If only the covariance matrices are given.
if(nargin<7||isempty(L2))
    L2=zeros(xDim,xDim,Nr);
    for i=1:Nr
        L2(:,:,i)=chol(P2(:,:,i),'lower');
    end
end

%If only the lower-triangular Cholesky decompositions of the covariance
%matrices are given.
if(isempty(P2))
    P2=zeros(xDim,xDim,Nr);
    for i=1:Nr
        P2(:,:,i)=L2(:,:,i)*L2(:,:,i)';
    end
end

%The formulae for the gradient of Jhr and Jrr are given in Equation 3.45 of
%Section 3.3.3.3 of [1].
gradJhr=zeros(xDim,xDim,Nr);
gradJrr=zeros(xDim,xDim,Nr);

for j=1:Nr
    for i=1:Nh
        diff=mu1(:,i)-mu2(:,j);
        PSum=P1(:,:,i)+P2(:,:,j);
        PSumInv=inv(PSum);

        gradJhr(:,:,j)=gradJhr(:,:,j)+w1(i)*PDFVals12(i,j)*PSumInv*(diff*diff'-PSum)*PSumInv*L2(:,:,j);
    end
    gradJhr(:,:,j)=w2(j)*gradJhr(:,:,j);
end

for j=1:Nr
    for i=1:Nr
        diff=mu2(:,i)-mu2(:,j);
        PSum=P2(:,:,i)+P2(:,:,j);
        PSumInv=inv(PSum);

        gradJrr(:,:,j)=gradJrr(:,:,j)+w2(i)*PDFVals22(i,j)*PSumInv*(diff*diff'-PSum)*PSumInv*L2(:,:,j);
    end
    gradJrr(:,:,j)=2*w2(j)*gradJrr(:,:,j);
end

gradISE=gradJrr-2*gradJhr;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function value = givens_condition ( n )

%*****************************************************************************80
%
%% GIVENS_CONDITION returns the L1 condition of the GIVENS matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real VALUE, the L1 condition.
%
  a_norm = n * n;

  if ( n == 1 )
    b_norm = 1.0;
  else
    b_norm = 2.0;
  end

  value = a_norm * b_norm;

  return
end


text:
function F = coth(F, varargin)
%COTH   Hyperbolic cotangent of a CHEBFUN.
%   COTH(F) computes the hyperbolic cotangent of the CHEBFUN F.
%
%   COTH(F, PREF) does the same but uses the CHEBFUNPREF object PREF when
%   computing the composition.
%
% See also ACOTH.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Call the compose method:
F = compose(F, @coth, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function pi_values_test ( )

%*****************************************************************************80
%
%% PI_VALUES_TEST demonstrates the use of PI_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2004
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PI_VALUES_TEST:\n' );
  fprintf ( 1, '  PI_VALUES returns values of\n' );
  fprintf ( 1, '  the PI function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N         PI(N)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, fn ] = pi_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %8d  %10d\n', n, fn )

  end

  return
end


text:
function [x,rank,info] = gelsd(a,b_in,rcond)
%GELSD   Computing eigen information.
%   [X,INFO] = GELSD(A,B,RCOND)
%
%	No input parameter are optional. If you want an easy to use interface,
%	see GELS.

    x = [];
	rank = [] ;
	info = [];

	if (nargin~=3),
		disp('Wrong number of input parameters.');
		return;
	end;

	    %Do all the boring checking
	if (~isnumeric(a)),
		disp('The matrix is not composed of numeric values.');
		return;
	elseif (isinteger(a)),
		a=double(a);
	end;

	m=int32(size(a,1));
	n=int32(size(a,2));
	nrhs=int32(size(b_in,2));
	lda=int32(m);
	ldb=int32(max(m,n));

	if (isDouble(b_in)),
		b=zeros(double(ldb),double(nrhs));
		b(1:size(b_in,1) , 1:size(b_in,2)) = b_in;
	else
		b=single(zeros(double(ldb),double(nrhs)));
		b(1:size(b_in,1) , 1:size(b_in,2)) = b_in;
	end;
	rank=int32(0);

	if (isComplex(a) || isComplex(b)),
		if (isDouble(a) || isDouble(b))

			if (~isComplex(b)),
				b=complex(b);
			end;
			if (~isComplex(a)),
				a=complex(a);
			end;
			if (~isDouble(b)),
				b=double(b);
			end;
			if (~isDouble(a)),
				a=double(a);
			end;
			s=zeros(double(min(m,n)),1);
			work=complex(zeros(1));
			rwork=zeros(1);
			iwork=int32(zeros(1));
			lwork=int32(-1);
			info=int32(-1);
			[ m, n, nrhs, a, lda, x, ldb, s, rcond, rank, work, lwork, rwork, iwork, info]=lapack_zgelsd(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, rwork, iwork, info);
			lwork=int32(work(1));
			work=complex(zeros(double(lwork),1));
			lrwork=int32(rwork(1));
			smlsiz=25;
			nlvl = max( 0, int32( log2( double(min( m,n ))/double(smlsiz+1) ) ) + 1 );
			lrwork = 10*min(m,n) + 2*min(m,n)*smlsiz + 8*min(m,n)*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);

			liwork = max(1, 3*min(m,n)*nlvl + 11*min(m,n));
			rwork=zeros(double(lrwork),1);
			iwork=int32(zeros(double(liwork),1));
			[ m, n, nrhs, a, lda, x, ldb, s, rcond, rank, work, lwork, rwork, iwork, info]=lapack_zgelsd(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, rwork, iwork, info);

		else %complex single
			if (~isComplex(b)),
				b=complex(b);
			end;
			if (~isComplex(a)),
				a=complex(a);
			end;
			s=single(zeros(double(min(m,n)),1));

			work=single(complex(zeros(1,1)));
			work(1)=sqrt(-1)+1;
			lwork=int32(single(-1));
			rwork=single(zeros(1));
			iwork=int32(zeros(1));
			info=int32(single(-1));
			[ m, n, nrhs, a, lda, x, ldb, s, rcond, rank, work, lwork, rwork, iwork, info]=lapack_cgelsd(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, rwork, iwork, info);
			lwork=int32(work(1));
			work=single(complex(zeros(double(lwork),1)));
			lrwork=int32(rwork(1));
			smlsiz=25;
			nlvl = max( 0, int32( log2( double(min( m,n ))/double(smlsiz+1) ) ) + 1 );
			lrwork = 10*min(m,n) + 2*min(m,n)*smlsiz + 8*min(m,n)*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);

			liwork = max(1, 3*min(m,n)*nlvl + 11*min(m,n));
			rwork=single(zeros(double(lrwork),1));
			iwork=int32(zeros(double(liwork),1));
			[ m, n, nrhs, a, lda, x, ldb, s, rcond, rank, work, lwork, rwork, iwork, info]=lapack_cgelsd(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, rwork, iwork, info);

		end;
	else %real double
		if (isDouble(a) || isDouble(b)),
			smlsiz=int32(32);
			nlvl = max( 0, int32( log2( double(min( m,n ))/double((double(smlsiz)+1)) )  + 1) );
			lwork= int32(11*double(max(m,n)) + 2*double(max(m,n))*double(smlsiz) + 8*double(max(m,n))*double(nlvl) + double(max(m,n))*double(nrhs));
			liwork=int32(3 * double(min(m,n)) * double(nlvl) + 11 * double(min(m,n)));
			iwork=int32(zeros(double(liwork),1));
			if (~isDouble(b)),
				b=double(b);
			end;
			if (~isDouble(a)),
				a=double(a);
			end;
			s=zeros(double(min(m,n)),1);
			info=int32(-1);
			iwork=int32(zeros(double(liwork),1));
			work=zeros(double(lwork),1);
			[ m, n, nrhs, a, lda, x, ldb, s, rcond, rank, work, lwork, iwork, info]=lapack_dgelsd(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, iwork, info);
		else %real single
			smlsiz=int32(25);
			nlvl = max( 0, int32( log2( double(min( m,n ))/double((smlsiz+1)) ) ) + 1 );
			liwork=int32(3 * min(m,n) * nlvl + 11 * min(m,n));
			lwork= int32(11*max(m,n) + 2*max(m,n)*smlsiz + 8*max(m,n)*nlvl + max(m,n)*nrhs);
			iwork=int32(zeros(double(liwork),1));
			work=single(zeros(double(lwork),1));
			s=single(zeros(double(min(m,n)),1));

			info=int32(-1);
			[ m, n, nrhs, a, lda, x, ldb, s, rcond, rank, work, lwork, iwork, info]=lapack_sgelsd(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, iwork, info);
		end;
	end;

	x=x(1:double(n),1:double(nrhs));


text:
function [out1,out2,out3] = UserInputFunction(x,y)
if nargin==0; x = 0;y = 0; end
%[-------------------------------------]%
% EDIT THE FUNCTION TO YOUR OWN CHOICE%
  G = x + y;
%[-------------------------------------]%
if nargout==2,
   out1 = x; out2 = y;  out3 = G;
else
   out1 = G;
end


text:
function [alpha, beta, gamma, loglik, xi, gamma2] = fwdback_twoslice(engine, init_state_distrib, transmat, obslik, varargin)
% FWDBACK Compute the posterior probs. in an HMM using the forwards backwards algo.
%
% [alpha, beta, gamma, loglik, xi, gamma2] = fwdback(init_state_distrib, transmat, obslik, ...)
%
% Notation:
% Y(t) = observation, Q(t) = hidden state, M(t) = mixture variable (for MOG outputs)
% A(t) = discrete input (action) (for POMDP models)
%
% INPUT:
% init_state_distrib(i) = Pr(Q(1) = i)
% transmat(i,j) = Pr(Q(t) = j | Q(t-1)=i)
%  or transmat{a}(i,j) = Pr(Q(t) = j | Q(t-1)=i, A(t-1)=a) if there are discrete inputs
% obslik(i,t) = Pr(Y(t)| Q(t)=i)
%   (Compute obslik using eval_pdf_xxx on your data sequence first.)
%
% Optional parameters may be passed as 'param_name', param_value pairs.
% Parameter names are shown below; default values in [] - if none, argument is mandatory.
%
% For HMMs with MOG outputs: if you want to compute gamma2, you must specify
% 'obslik2' - obslik(i,j,t) = Pr(Y(t)| Q(t)=i,M(t)=j)  []
% 'mixmat' - mixmat(i,j) = Pr(M(t) = j | Q(t)=i)  []
%
% For HMMs with discrete inputs:
% 'act' - act(t) = action performed at step t
%
% Optional arguments:
% 'fwd_only' - if 1, only do a forwards pass and set beta=[], gamma2=[]  [0]
% 'scaled' - if 1,  normalize alphas and betas to prevent underflow [1]
% 'maximize' - if 1, use max-product instead of sum-product [0]
%
% OUTPUTS:
% alpha(i,t) = p(Q(t)=i | y(1:t)) (or p(Q(t)=i, y(1:t)) if scaled=0)
% beta(i,t) = p(y(t+1:T) | Q(t)=i)*p(y(t+1:T)|y(1:t)) (or p(y(t+1:T) | Q(t)=i) if scaled=0)
% gamma(i,t) = p(Q(t)=i | y(1:T))
% loglik = log p(y(1:T))
% xi(i,j,t-1)  = p(Q(t-1)=i, Q(t)=j | y(1:T))
% gamma2(j,k,t) = p(Q(t)=j, M(t)=k | y(1:T)) (only for MOG  outputs)
%
% If fwd_only = 1, these become
% alpha(i,t) = p(Q(t)=i | y(1:t))
% beta = []
% gamma(i,t) = p(Q(t)=i | y(1:t))
% xi(i,j,t-1)  = p(Q(t-1)=i, Q(t)=j | y(1:t))
% gamma2 = []
%
% Note: we only compute xi if it is requested as a return argument, since it can be very large.
% Similarly, we only compute gamma2 on request (and if using MOG outputs).
%
% Examples:
%
% [alpha, beta, gamma, loglik] = fwdback(pi, A, multinomial_prob(sequence, B));
%
% [B, B2] = mixgauss_prob(data, mu, Sigma, mixmat);
% [alpha, beta, gamma, loglik, xi, gamma2] = fwdback(pi, A, B, 'obslik2', B2, 'mixmat', mixmat);

if nargout >= 5, compute_xi = 1; else compute_xi = 0; end
if nargout >= 6, compute_gamma2 = 1; else compute_gamma2 = 0; end

[obslik2, mixmat, fwd_only, scaled, act, maximize, compute_xi, compute_gamma2] = process_options(varargin, 'obslik2', [], 'mixmat', [], 'fwd_only', 0, 'scaled', 1, 'act', [], 'maximize', 0, 'compute_xi', compute_xi, 'compute_gamma2', compute_gamma2);

[Q T] = size(obslik);

if isempty(obslik2)
  compute_gamma2 = 0;
end

if isempty(act)
  act = ones(1,T);
  transmat = { transmat } ;
end

scale = ones(1,T);

% scale(t) = Pr(O(t) | O(1:t-1)) = 1/c(t) as defined by Rabiner (1989).
% Hence prod_t scale(t) = Pr(O(1)) Pr(O(2)|O(1)) Pr(O(3) | O(1:2)) = Pr(O(1), ... ,O(T))
% or log P = sum_t log scale(t).
% Rabiner suggests multiplying beta(t) by scale(t), but we can instead
% normalise beta(t) - the constants will cancel when we compute gamma.

loglik = 0;

alpha = zeros(Q,T);
gamma = zeros(Q,T);
if compute_xi
  xi = zeros(Q,Q,T-1);
else
  xi = [];
end

%%%%%%%%% Forwards %%%%%%%%%%

t = 1;
alpha(:,1) = init_state_distrib(:) .* obslik(:,t);
if scaled
  %[alpha(:,t), scale(t)] = normaliseC(alpha(:,t));
  [alpha(:,t), scale(t)] = normalise(alpha(:,t));
end
if scaled, assert(approxeq(sum(alpha(:,t)),1)), end
for t=2:T
  %trans = transmat(:,:,act(t-1))';
  trans = transmat{act(t-1)};
  if maximize
    m = max_mult(trans', alpha(:,t-1));
    %A = repmat(alpha(:,t-1), [1 Q]);
    %m = max(trans .* A, [], 1);
  else
    m = trans' * alpha(:,t-1);
  end
  alpha(:,t) = m(:) .* obslik(:,t);
  if scaled
    %[alpha(:,t), scale(t)] = normaliseC(alpha(:,t));
    [alpha(:,t), scale(t)] = normalise(alpha(:,t));
  end
  if compute_xi & fwd_only  % useful for online EM
    %xi(:,:,t-1) = normaliseC((alpha(:,t-1) * obslik(:,t)') .* trans);
    xi(:,:,t-1) = normalise((alpha(:,t-1) * obslik(:,t)') .* trans);
  end
  if scaled, assert(approxeq(sum(alpha(:,t)),1)), end
end
if scaled
  if any(scale==0)
    loglik = -inf;
  else
    loglik = sum(log(scale));
  end
else
  loglik = log(sum(alpha(:,T)));
end

if fwd_only
  gamma = alpha;
  beta = [];
  gamma2 = [];
  return;
end

%%%%%%%%% Backwards %%%%%%%%%%

beta = zeros(Q,T);
if compute_gamma2
  M = size(mixmat, 2);
  gamma2 = zeros(Q,M,T);
else
  gamma2 = [];
end

beta(:,T) = ones(Q,1);
%gamma(:,T) = normaliseC(alpha(:,T) .* beta(:,T));
gamma(:,T) = normalise(alpha(:,T) .* beta(:,T));
t=T;
if compute_gamma2
  denom = obslik(:,t) + (obslik(:,t)==0); % replace 0s with 1s before dividing
  gamma2(:,:,t) = obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M]) ./ repmat(denom, [1 M]);
  %gamma2(:,:,t) = normaliseC(obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M])); % wrong!
end
for t=T-1:-1:1
  b = beta(:,t+1) .* obslik(:,t+1);
  %trans = transmat(:,:,act(t));
  trans = transmat{act(t)};
  if maximize
    B = repmat(b(:)', Q, 1);
    beta(:,t) = max(trans .* B, [], 2);
  else
    beta(:,t) = trans * b;
  end
  if scaled
    %beta(:,t) = normaliseC(beta(:,t));
    beta(:,t) = normalise(beta(:,t));
  end
  %gamma(:,t) = normaliseC(alpha(:,t) .* beta(:,t));
  gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
  if compute_xi
    %xi(:,:,t) = normaliseC((trans .* (alpha(:,t) * b')));
    xi(:,:,t) = normalise((trans .* (alpha(:,t) * b')));
    %xi(:,:,t) = (trans .* (alpha(:,t) * b'));
  end
  if compute_gamma2
    denom = obslik(:,t) + (obslik(:,t)==0); % replace 0s with 1s before dividing
    gamma2(:,:,t) = obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M]) ./ repmat(denom, [1 M]);
    %gamma2(:,:,t) = normaliseC(obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M]));
  end
end

% We now explain the equation for gamma2
% Let zt=y(1:t-1,t+1:T) be all observations except y(t)
% gamma2(Q,M,t) = P(Qt,Mt|yt,zt) = P(yt|Qt,Mt,zt) P(Qt,Mt|zt) / P(yt|zt)
%                = P(yt|Qt,Mt) P(Mt|Qt) P(Qt|zt) / P(yt|zt)
% Now gamma(Q,t) = P(Qt|yt,zt) = P(yt|Qt) P(Qt|zt) / P(yt|zt)
% hence
% P(Qt,Mt|yt,zt) = P(yt|Qt,Mt) P(Mt|Qt) [P(Qt|yt,zt) P(yt|zt) / P(yt|Qt)] / P(yt|zt)
%                = P(yt|Qt,Mt) P(Mt|Qt) P(Qt|yt,zt) / P(yt|Qt)
%


text:
% FT_POSTAMBLE_DEBUG is a helper script for debugging problems with FieldTrip
% functions. It is executed at normal termination of the calling function and
% disables the "onCleanup" function.
%
% Use as
%   ft_preamble debug
%   .... regular code goes here ...
%   ft_postamble debug
%
% See also FT_PREAMBLE, FT_POSTAMBLE, FT_PREAMBLE_DEBUG, DEBUGCLEANUP

% these variables are shared by the three debug handlers
global Ce9dei2ZOo_debug Ce9dei2ZOo_funname Ce9dei2ZOo_argin Ce9dei2ZOo_ws Ce9dei2ZOo_ns Ce9dei2ZOo_is Ce9dei2ZOo_ds

if isfield(cfg, 'verbose') && ischar(cfg.verbose)
  % restore the previous state of the notifications
  ft_warning(Ce9dei2ZOo_ws)
  ft_notice(Ce9dei2ZOo_ns);
  ft_info(Ce9dei2ZOo_is);
  ft_debug(Ce9dei2ZOo_ds);
end

if ~isfield(cfg, 'debug')
  % do not provide extra debugging facilities
  return
end

switch cfg.debug
  case {'displayonsuccess' 'saveonsuccess'}
    % do not clean up the global variables yet
    % these are still needed by the cleanup function

  otherwise
    % stack(1) is this script
    % stack(2) is the calling ft_postamble function
    % stack(3) is the main FieldTrip function that we are interested in
    Ce9dei2ZOx_funname = dbstack;
    Ce9dei2ZOx_funname = Ce9dei2ZOx_funname(3).name;

    if isequal(Ce9dei2ZOx_funname, Ce9dei2ZOo_funname)
      % this results in the cleanup function doing nothing
      Ce9dei2ZOo_debug   = 'no';
      Ce9dei2ZOo_funname = [];
      Ce9dei2ZOo_argin   = [];
    else
      % this happens for nested FieldTrip functions, e.g. when
      % ft_selectdata is called by another high-level function
    end
end


text:
function varargout = poly(varargin) %#ok<STOUT>
%POLY   POLY does not support SINGFUN objects.
%
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

error('CHEBFUN:SINGFUN:poly:noSupport', ...
      'POLY does not support SINGFUN objects.')

end


text:
%function cstr = strsplit(strarr, ch)
%
% FUNCTIONALITY
% Splits string <strarr> at positions where substring <ch>
% occurs. If <ch> is not given, then ascii character 10
% (newline) is the splitting position.
% This would split for example a read-in file into its lines.
%
% INPUT
%   strarr     string to split
%   ch         optional, character or substring at which to split,
%              default is character ascii 10 = 'newline'
%
% OUTPUT
%   cstr       cell array of strings, all leading and trailing
%              spaces are deleted. Substring ch is removed from string.
%
% EXAMPLES
%   a = strsplit('aaa bbb ccc ', ' ')
%   Result: a{1}='aaa', a{2}='bbb', a{3}='ccc'
%   a = strsplit('user@server', '@')
%   Result: a{1}='user', a{3}='server'
%   a = strsplit('aabbccddeeffgg', 'dd')
%   Result: a{1}='aabbcc', a{2} = 'eeffgg'
%
% RELATED
%   strfun, depreciated stringify.m

% Copyright (C) 2002 Marc Molinari, University of Southampton
% Author: Dr Marc Molinari <m.molinari@soton.ac.uk>
% $Revision: 1.1 $ $Date: 2005/04/15 17:12:14 $ $Tag$

author: Cloud-CV

repo: object-proposals


text:
function point = intersectLines(line1, line2, varargin)
%INTERSECTLINES Return all intersection points of N lines in 2D.
%
%   PT = intersectLines(L1, L2);
%   returns the intersection point of lines L1 and L2. L1 and L2 are 1-by-4
%   row arrays, containing parametric representation of each line (in the
%   form [x0 y0 dx dy], see 'createLine' for details).
%
%   In case of colinear lines, returns [Inf Inf].
%   In case of parallel but not colinear lines, returns [NaN NaN].
%
%   If each input is [N*4] array, the result is a [N*2] array containing
%   intersections of each couple of lines.
%   If one of the input has N rows and the other 1 row, the result is a
%   [N*2] array.
%
%   PT = intersectLines(L1, L2, EPS);
%   Specifies the tolerance for detecting parallel lines. Default is 1e-14.
%
%   Example
%   line1 = createLine([0 0], [10 10]);
%   line2 = createLine([0 10], [10 0]);
%   point = intersectLines(line1, line2)
%   point =
%       5   5
%
%   See also
%   lines2d, edges2d, intersectEdges, intersectLineEdge
%   intersectLineCircle

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2003-10-31
% Copyright 2003-2022 INRA - TPV URPOI - BIA IMASTE

%% Process input arguments

% extract tolerance
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

% check size of each input
N1 = size(line1, 1);
N2 = size(line2, 1);
N = max(N1, N2);
if N1 ~= N2 && N1*N2 ~= N
    error('matGeom:IntersectLines:IllegalArgument', ...
        'The two input arguments must have same number of lines');
end

%% Check parallel and colinear lines

% coordinate differences of origin points
dx = bsxfun(@minus, line2(:,1), line1(:,1));
dy = bsxfun(@minus, line2(:,2), line1(:,2));

% indices of parallel lines
denom = line1(:,3) .* line2(:,4) - line2(:,3) .* line1(:,4);
par = abs(denom) < tol;

% indices of colinear lines
col = abs(dx .* line1(:,4) - dy .* line1(:,3)) < tol & par ;

% initialize result array
x0 = zeros(N, 1);
y0 = zeros(N, 1);

% initialize result for parallel lines
x0(col) = Inf;
y0(col) = Inf;
x0(par & ~col) = NaN;
y0(par & ~col) = NaN;

% in case all line couples are parallel, return
if all(par)
    point = [x0 y0];
    return;
end

%% Extract coordinates of itnersecting lines

% indices of intersecting lines
inds = ~par;

% extract base coordinates of first lines
if N1 > 1
    line1 = line1(inds,:);
end
x1 =  line1(:,1);
y1 =  line1(:,2);
dx1 = line1(:,3);
dy1 = line1(:,4);

% extract base coordinates of second lines
if N2 > 1
    line2 = line2(inds,:);
end
x2 =  line2(:,1);
y2 =  line2(:,2);
dx2 = line2(:,3);
dy2 = line2(:,4);

% re-compute coordinate differences of origin points
dx = bsxfun(@minus, line2(:,1), line1(:,1));
dy = bsxfun(@minus, line2(:,2), line1(:,2));

%% Compute intersection points

denom = denom(inds);
x0(inds) = (x2 .* dy2 .* dx1 - dy .* dx1 .* dx2 - x1 .* dy1 .* dx2) ./ denom ;
y0(inds) = (dx .* dy1 .* dy2 + y1 .* dx1 .* dy2 - y2 .* dx2 .* dy1) ./ denom ;

% concatenate result
point = [x0 y0];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function Cw = getCameraCenter(camRtC2W)

Cw = camRtC2W(1:3,4);

author: jianxiongxiao

repo: ProfXkit


text:
function [X, R, S] = snapshot_gen_sto(design, doas, wavelength, t, ncov, scov)
%SNAPSHOT_GEN_STO Generates snapshots for the stochastic model.
%Syntax:
%   X = STO_SNAPSHOT_GEN(design, doas, wavelength[, t, ncov, scov]);
%   [X, R] = STO_SNAPSHOT_GEN(design, doas, wavelength[, t, ncov, scov]);
%   [X, R, S] = STO_SNAPSHOT_GEN(design, doas, wavelength[, t, ncov, scov]);
%Inputs:
%   design - Array design.
%   doas - DOA vector. For 2D DOAs, each column represents a DOA pair.
%   wavelength - Wavelength.
%   t - Number of snapshots.
%   ncov - Covariance matrix of the additive complex circular-symmetric
%          Gaussian noise. Can be a scalar, vector (for uncorrelated noise
%          with different powers), or a matrix.
%   scov - Covariance matrix of the source signals. Can be a scalar, vector
%          (for uncorrelated sources with different powers), or a matrix.
%Outputs:
%   X - Snapshots, where each columns is a single snapshot.
%   R - Sample covariance matrix (averaged by the number of snapshots).
%   S - A source_count x snapshot_count matrix consists of source signal
%       vectors.
if nargin <= 5
    scov = 1;
end
if nargin <= 4
    ncov = 1;
end
if nargin <= 3
    t = 1;
end
A = steering_matrix(design, wavelength, doas);
[m, k] = size(A);
S_internal = gen_ccsg(k, t, scov);
X = A * S_internal + gen_ccsg(m, t, ncov);
if nargout >= 2
    R = (X*X')/t;
    if nargout == 3
        S = S_internal;
    end
end
end

function X = gen_ccsg(m, n, cov)
X0 = randn(m, n) + 1j * randn(m, n);
if isscalar(cov)
    X = sqrt(cov/2) * X0;
elseif isvector(cov)
    X = bsxfun(@times, X0, cov(:));
else
    C = sqrtm(cov/2);
    X = C*X0;
end
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function [xUpdate,PUpdate,innov,Pzz,W]=QMCKalUpdateWithPred(z,R,zPred,PzPred,otherInfo)
%%QMCKALUPDATEWITHPRED Given the output of the measurement prediction step
%           from QMCKalMeasPred and a measurement, complete the measurement
%           update step of the quasi-Monte Carlo Kalman filter with
%           additive measurement noise. Separating the measurement
%           prediction step from the rest of the update step can make the
%           creation of multiple measurement association hypotheses from a
%           single target prediction more efficient. The full measurement
%           update function is QMCKalUpdate.
%
%INPUTS: z The zDimX1 vector measurement.
%        R The zDim X zDim measurement covariance matrix in the native
%          coordinate system of the measurement.
%    zPred The zDimXnumComp measurement predictions from the filter.
%   PzPred The zDimXzDimXnumComp covariance matrices associated with zPred.
% otherInfo The intermediate results returned in the otherInfo output of
%          the QMCKalMeasPred function.
%
%OUTPUTS: xUpdate The xDimXnumComp updated state vectors.
%         PUpdate The updated xDimXxDimXnumComp state covariance matrices
%                 associated with xUpdate.
%      innov, Pzz The zDimXnumComp innovations and the zDimXzDim innovation
%                 covariance matrices are returned in case one wishes to
%                 analyze the consistency of the estimator or use those
%                 values in gating or likelihood evaluation.
%               W The xDimXzDimXnumComp gains used in the update.
%
%See the comments to the function QMCKalMeasPred for an example of usage of
%this function. See the comments to QMCKalUpdate for more information on
%the algorithm.
%
%June 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

innovTrans=otherInfo.innovTrans;
stateTrans=otherInfo.stateTrans;
xPred=otherInfo.xPred;
PPred=otherInfo.PPred;
Pxz=otherInfo.Pxz;

xDim=size(xPred,1);
numComp=size(xPred,2);
zDim=size(z,1);

xUpdate=zeros(xDim,numComp);
PUpdate=zeros(xDim,xDim,numComp);
innov=zeros(zDim,numComp);
Pzz=zeros(zDim,zDim,numComp);
W=zeros(xDim,zDim,numComp);

for k=1:numComp
    Pzz(:,:,k)=PzPred(:,:,k)+R;

    %The innovation, transformed as necessary to keep values in a desired
    %range.
    innov(:,k)=innovTrans(z,zPred(:,k));

    %The filter gain
    W(:,:,k)=Pxz(:,:,k)/Pzz(:,:,k);

    %Updated state estimate
    xUpdate(:,k)=stateTrans(xPred(:,k)+W(:,:,k)*innov(:,k));

    %Updated state covariance matrix
    PUpdate(:,:,k)=PPred(:,:,k)-W(:,:,k)*Pzz(:,:,k)*W(:,:,k)';
end
end

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


text:
%% Converting sparse tensors to matrices and vice versa
% We show how to convert a sptensor to a matrix stored in _coordinate_
% format and with extra information so that it can be converted back to a
% sptensor.

%% Creating a sptenmat (sparse tensor as sparse matrix) object
% A sparse tensor can be converted to a sparse matrix. The matrix, however,
% is not stored as a MATLAB sparse matrix because that format is sometimes
% inefficient for converted sparse tensors. Instead, the row and column
% indices are stored explicitly.
%%
% First, we create a sparse tensor to be converted.
X = sptenrand([10 10 10 10],10) %<-- Generate some data.
%%
% All the same options for tenmat are available as for tenmat.
A = sptenmat(X,1) %<-- Mode-1 matricization.
%%
A = sptenmat(X,[2 3]) %<-- More than one mode is mapped to the columns.
%%
A = sptenmat(X,[2 3],'t') %<-- Specify column dimensions (transpose).
%%
A = sptenmat(X,1:4) %<-- All modes mapped to rows, i.e., vectorize.
%%
A = sptenmat(X,2) %<-- By default, columns are ordered as [1 3 4].
%%
A = sptenmat(X,2,[3 1 4]) %<-- Explicit column ordering.
%%
A = sptenmat(X,2,'fc') %<-- Foward cyclic.
%%
A = sptenmat(X,2,'bc') %<-- Backward cyclic.
%% Constituent parts of a sptenmat
A.subs %<-- Subscripts of the nonzeros.
%%
A.vals %<-- The corresponding nonzero values.
%%
A.tsize %<-- Size of the original tensor.
%%
A.rdims %<-- Dimensions that were mapped to the rows.
%%
A.cdims %<-- Dimensions that were mapped to the columns.
%% Creating a sptenmat from its constituent parts
B = sptenmat(A.subs,A.vals,A.rdims,A.cdims,A.tsize) %<-- Copies A
%%
B = sptenmat(double(A),A.rdims,A.cdims,A.tsize) %<-- More efficient to pass a matrix.
%% Creating a sptenmat with no nonzeros
A = sptenmat([],[],A.rdims,A.cdims,A.tsize) %<-- An empty sptenmat.
%% Creating an emtpy sptenmat
A = sptenmat %<-- A really empty sptenmat.
%% Use double to convert a sptenmat to a MATLAB sparse matrix
X = sptenrand([10 10 10 10],10); %<-- Create a tensor.
A = sptenmat(X,1) %<-- Convert it to a sptenmat
%%
B = double(A) %<-- Convert it to a MATLAB sparse matrix
%%
whos A B %<-- The storage for B (the sparse matrix) is larger than for A.
%%
C = B'; %<-- Transposing the result fixes the problem.
whos C
%% Use full to convert a sptenmat to a tenmat
B = sptenmat(sptenrand([3 3 3], 3), 1) %<-- Create a sptenmat
%%
C = full(B) %<-- Convert to a tenmat
%% Use sptensor to convert a sptenmat to a sptensor
Y = sptensor(A) %<-- Convert a sptenmat to a sptensor
%% Use size and tsize for the dimensions of a sptenmat
size(A) %<-- Matrix size
tsize(A) %<-- Corresponding tensor size
%% Subscripted reference for a sptenmat
% This is not supported beyond getting the constituent parts.
%% Subscripted assignment for a sptenmat
A(1:2,1:2) = ones(2) %<-- Replace part of the matrix.
%% Use end for the last index
% End is not supported.
%% Basic operations for sptenmat
norm(A) %<-- Norm of the matrix.
%%
+A %<-- Calls uplus.
%%
-A %<-- Calls uminus.
%% Use aatx to efficiently compute A * A' * x for a sptenmat
x = ones(10,1); %<-- Create vector
aatx(A,x) %<-- Compute A * A' * x
%%
double(A) * double(A)' * x %<-- Same as above but less efficient
%% Displaying a tenmat
% Shows the original tensor dimensions, the modes mapped to rows, the modes
% mapped to columns, and the matrix.
disp(A)

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function r = eq(p,q)
% this function compares the equality condition of x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if (isa(p, 'CellVector')&&isa(q, 'CellVector'))
    r=p;
    r.xvalue = p.xvalue==q.xvalue;
    r.yvalue = p.yvalue==q.yvalue;
    r.zvalue = p.zvalue==q.zvalue;
elseif isa(p, 'CellVector')
    r=p;
    r.xvalue = p.xvalue==q;
    r.yvalue = p.yvalue==q;
    r.zvalue = p.zvalue==q;
else
    r=q;
    r.xvalue = p==q.xvalue;
    r.yvalue = p==q.yvalue;
    r.zvalue = p==q.zvalue;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
%% Projection Based Shape Parameters
%
%%
% In this section we discuss shape parameters grains that depend on one
% dimensional projections, i.e.,
%
% || <grain2d.caliper.html |caliper|>  || caliper or Feret diameter in $\mu m$ || <grain2d.diameter.html |diameter|>  || diameter in $\mu m$ ||
%
% In order to demonstrate these parameters we first import a small sample
% data set.

% load sample EBSD data set
mtexdata forsterite silent

% reconstruct grains, discard boudnary grains and smooth them
[grains, ebsd.grainId] = calcGrains(ebsd('indexed'),'angle',5*degree);
ebsd(grains(grains.grainSize<5)) = [];
[grains, ebsd.grainId] = calcGrains(ebsd('indexed'),'angle',5*degree);
grains(grains.isBoundary) = [];

grains = smooth(grains('indexed'),10,'moveTriplePoints');

% plot all grains and highlight a specific one
plot(grains)

ind = 654;
hold on
plot(grains(ind).boundary,'lineWidth',5,'linecolor','blue')
hold off

%%
% The most well known projection based parameter is the
% <grain2d.diamter.html |diameter|> which refers to the longest distance
% between any two boundary points and is given in $\mu m$.

grains(ind).diameter

%%
% The diameter is a special case of the <grain2d.caliper.html |caliper|> or
% Feret diameter of a grain. By definition the caliper is the length of a
% grain when projected onto a line. We may trace the caliper with respect
% to projection direction

close all
omega = linspace(0,180);
plot(omega,grains(ind).caliper(omega*degree),'LineWidth',2)
ylabel('length in $\mu$m','Interpreter','latex')
xlabel('angle of the projection line in degree')
xlim([0,180])

%%
% We observe that that maximum caliper is about 7000 while the minimum
% caliper is about 2000. We may compute the exact direction and length of
% the maximum or minimum by passing the options |'longest'| or |'shortest'|
% to the function <grain2d.caliper.html |caliper|>. In this case the the
% output is of type @vector3d indicating the direction. The
% <vector3d.norm.html |norm|> of the vector coincides with the caliper for
% this projection direction. Hence, the |norm(grains.caliper('longest'))|
% coincides with the diameter.

plot(grains(ind),'micronbar','off')
legend('off')

norm(grains(ind).caliper('longest'))
norm(grains(ind).caliper('shortest'))

hold on
quiver(grains(ind),grains(ind).caliper('longest'),'noScaling')
quiver(grains(ind),grains(ind).caliper('shortest'),'noScaling')
hold off

%%
% The difference between the longest and the shortest caliper can be taken
% as a measure how round a grain is.

cMin = grains.caliper('shortest');
cMax = grains.caliper('longest');

plot(grains,(norm(cMax) - norm(cMin))./norm(cMax),'micronbar','off')
mtexColorbar('title','TODO')

%%
% This longest and shortest caliper are comparable to
% <grain2d.longAxis.html |grains.longAxis|> and <grain2d.shortAxis.html
% |grains.shortAxis|> computed from an ellipse fitted to the grain. In the
% case of rectangular particles, one might not primarily be interested in
% the longest caliper of a grain but rather in the direction normal to the
% shortest caliper. This is computed when specifying the option
% |'shortestPerp'|. If we imagine a very strong alignment of the long axes
% of orthorhombic particles, the maximum diameter may show a bimodal
% distribution (the two, roughly equally distributed diagonals of the
% particle).

% load some test grains
testgrains = mtexdata('testgrains');
testgrains = smooth(testgrains([6 8]),10);

% compute the longest caliper and the caliper perpendicular to the shortest
cMax = testgrains.caliper('longest');
cMinPerp = testgrains.caliper('shortestPerp');

% plot the grains and visualize the different long axes
plot(testgrains,'micronbar','off','lineWidth',2)
hold on
quiver(testgrains,cMax,'DisplayName','longest calliper','LineWidth',3)
quiver(testgrains,testgrains.longAxis,'DisplayName','long axis','LineWidth',3)
quiver(testgrains,cMinPerp,'DisplayName','perp to shortest','LineWidth',3)
hold off
legend('Location','east')

%% PAROR and SURFOR
%
% Another way of quantifying shape farbics is by making use of the
% cumulative projection function of the grains or the grain boundary
% segments. These methods are heavily inspired by
% <https://en.wikipedia.org/wiki/Flatland Edwin A. Abbotts 'Flatland - A
% romance of many dimensions' (1884)> and based on Panozzo, R., 1983,
% "Two-dimensional analysis of shape fabric using projections of digitized
% lines in a plane". Tectonophysics 95, 279-294. and Panozzo, R., 1984,
% "Two-dimensional strain from the orientation of lines in a plane." J.
% Struct. Geol. 6, 215-221. implemented in Mtex as <grain2d.paror.html
% |grains.paror|> and <grainBoundary.surfor.html |grainBoudnary.surfor|>
%
% As mentioned above the function <grain2d.caliper.html |caliper|> can be
% called with a list of angles and returns the projection length of all
% grains with respect to all angles.

% projection angle
omega = linspace(0,360*degree,361);
c = grains('Fo').caliper(omega);

subplot(1,2,1)
polarplot(omega,c,'LineWidth',2,'color',[0 0.25 0.5 0.25])
title('Forsterite')

% take the average
hold on
polarplot(omega,5*mean(c),'LineWidth',3,'color','k');
hold off

subplot(1,2,2)
c = grains('Enstatite').caliper(omega);

polarplot(omega,c,'LineWidth',2,'color',[0 0.25 0.5 0.25])
title('Enstatite')

% take the average
hold on
polarplot(omega,5*mean(c),'LineWidth',3,'color','k');
hold off

%%
% The above averaged caliper can be computed more directly by the function
% <grain2d.paror.html |grains.paror|> which returns the cumulative particle
% projection function normalized to 1. The projection angles can be
% regarded as the rotation angle of the particle (counterclockwise) while
% projecting from the y-axis onto the x-axis.

close all
cumplF = paror(grains('fo'),omega);
cumplE = paror(grains('en'),omega);

subplot(1,2,1)
polarplot(omega,cumplF,'LineWidth',3,'color','k')

subplot(1,2,2)
polarplot(omega,cumplE,'LineWidth',3,'color','k')

%%
% We can interpret the results in the following way. The minimum of the
% curve is a measure of the amplitude of the projection function and can be
% compared to an averaged axial ratio $b/a$ of the entire  fabric;
% isotropic fabrics would have a $b/a$ close to 1 while highly anisotropic
% fabrics can be identified by small $b/a$ values.

min(cumplF), min(cumplE)

%%
% The position of the maxima and minima of the projection function derived
% from <grain2d.paror.html |paror|> can be interpreted in the following
% way: the maximum position represents the preferred axis parallel to the
% longest projection and the normal the minimum position represents the
% preferred axis related to the normal to the shortest projection function.

% for the Forsterite
[~, id_max] = max(cumplF);
[~, id_min] = min(cumplF);

mod(omega(id_max)./degree,180)
mod(omega(id_min)./degree-90,180)

% for the Enstatite
[~, id_max] = max(cumplE);
[~, id_min] = min(cumplE);

mod(omega(id_max)./degree,180)
mod(omega(id_min)./degree-90,180)

%%
% The smaller the difference between these values, the closer the fabric is
% to an orthorhombic symmetry.
%
% Similarly to using the entire particle (the convex hull in case of the
% projection functions), we can use a distribution of lines
% <grainBoundary.surfor.html |grainBoudnary.surfor|>. This can be useful
% for the quantification of the grain boundary anisotropy or in general
% might be needed if we look at boundaries which do not form closed
% outlines, e.g. a list of subgrain or twin boundaries or the contact
% between certain phases.
%
% Let's compare the boundaries between the different unlike phases and
% between forsterite-forsterite in our sample:

close all
pairs = [1 1; nchoosek(1:3,2)];
phase = {'Fo' 'En' 'Di'};
for i=1:length(pairs)

  gB = grains.boundary(phase{pairs(i,:)});
  polarplot(omega, surfor(gB,omega), 'linewidth',2, ...
    'DisplayName',[phase{pairs(i,1)} '-' phase{pairs(i,2)}]);
  hold on

end
hold off
legend('Location','southoutside','Orientation','horizontal')

%%
% We can see that Forsterite-Forsterite boundaries form a fabric slightly
% more inclined with respect to the other phase boundariesand that the
% phase boundaries between the two pyroxenes (Enstatite and Diopside) show
% the lowest anisotropy.
%
%% Characteristic Shape
%
% The characteristic shape results from the cummulative sum of all grain
% boundary segements ordered by the angle of the segment direction. It can
% be regarded as to represent the average grain shape, however without the
% need to use closed areas such as it would be required when working with
% grains.
%
% Here we can compare the shape defined by Forterite-Forsterite,
% Enstatite-Enstatite and Forsterite-Enstatite boundaries

plotopts = {'normalize','linewidth',2, 'plain'};

shapeF = characteristicShape(grains.boundary('Fo','Fo'))
plot(shapeF,plotopts{:},'DisplayName','Fo-Fo')

hold on
shapeE = characteristicShape(grains.boundary('En','En'));
plot(shapeE,plotopts{:},'DisplayName','En-En')

hold on
shapeEF = characteristicShape(grains.boundary('En','Fo'));
plot(shapeEF,plotopts{:},'DisplayName','En-Fo')
hold off

legend('Location','southoutside','Orientation','horizontal')

%%
% The output of the command <grainBoundary.characteristicShape.html
% |characteristicShape|> is a <shape2d.shape2d.html |shape2d|> object which
% behaves very similar to a <grain2d.grain2d.html |grain2d|> object. Hence
% it is easy to derive things such as a long axis or e.g. the angle between
% the longest and the shortest caliper which can be regarded as a measure
% of asymmetry.

angle(shapeF.caliper('longest'),shapeF.caliper('shortest'))/degree
angle(shapeE.caliper('longest'),shapeE.caliper('shortest'))/degree
angle(shapeEF.caliper('longest'),shapeF.caliper('shortest'))/degree

author: mtex-toolbox

repo: mtex

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


text:
function im=seamcarving(im,k)
%% illustrative example of Seam carving for content aware image resizing
%
%
% usage: carvedimg=seamcarving(im,k)
%
% k is how many vertical seams to remove.
% im is the image.
%
% example:
%   img=imread('peppers.png')
%   carvedimg=seamcarving(img,50)
%   image([carvedimg img]);
%   axis equal;
%
% Author: Aslak Grinsted 2007...
% Based on ideas from Avidan & Shamir:
% http://video.google.com/videoplay?docid=-6221880321193117495
% Note i havent read their paper and they have probably lots of smart tricks
% for optimizations.
%

demo=nargout==0;
if nargin==0
    fex={'peppers.png' 'liftingbody.png' 'pears.png' 'trees.tif' 'football.jpg' 'onion.png'};
    fex=fex{ceil(rand*length(fex))};
    try
        [im,map]=imread(fex);
    catch
        [im,map]=imread('street1.jpg');
    end
    k=50;
end
im=im2double(im);

if demo
    close(findobj(0,'type','figure','tag','seam carving demo'));
    figure; set(gcf,'tag','seam carving demo','name','Seam Carving','NumberTitle','off')
    axes('position', [0 0 1 1]);
    if size(im,3)==1
        im=im/max(im(:));
        him=imagesc(im);
        colormap gray
    else
        him=image(im);
    end
    origim=im;
    axis equal
    axis off
end

for jj=1:k
    G=costfunction(im);
    %find shortest path in G
    Pot=G;
    for ii=2:size(Pot,1)
        pp=Pot(ii-1,:);
        ix=pp(1:end-1)<pp(2:end);
        pp([false ix])=pp(ix);
        ix=pp(2:end)<pp(1:end-1);
        pp(ix)=pp([false ix]);
        Pot(ii,:)=Pot(ii,:)+pp;
    end

    %Walk down hill
    pix=zeros(size(G,1),1);
    [mn,pix(end)]=min(Pot(end,:));
    pp=find(Pot(end,:)==mn);
    pix(end)=pp(ceil(rand*length(pp)));

    im(end,pix(end),:)=nan;
    for ii=size(G,1)-1:-1:1
        %[mn,gg]=min(Pot(ii,pix+(-1:1)));
        [mn,gg]=min(Pot(ii,max(pix(ii+1)-1,1):min(pix(ii+1)+1,end)));
        pix(ii)=gg+pix(ii+1)-1-(pix(ii+1)>1);
        im(ii,pix(ii),:)=bitand(ii,1);
%        G(ii,pix(ii))=1;
    end

    if demo
        set(him,'CDATA',im);
        %set(him,'CDATA',G,'CDataMapping','scaled');
        drawnow;
    end

    %remove seam from im & G:
    for ii=1:size(im,1)
%        G(ii,pix(ii):end-1)=G(ii,pix(ii)+1:end);
        im(ii,pix(ii):end-1,:)=im(ii,pix(ii)+1:end,:);
    end
    im(:,end,:)=[];
%    G(:,end)=[];

end

if demo
    set(him,'CDATA',[im origim])
    axis tight
end
if nargout==0
    clear im
end

    function G=costfunction(im) %%(xi,yi)
            G=zeros(size(im,1),size(im,2));
            for ii=1:size(im,3)
                %G=G+abs(filter2([1 0 -1],im(:,:,ii)))+abs(filter2([1;0;-1],im(:,:,ii)));
                G=G+(filter2([.5 1 .5; 1 -6 1; .5 1 .5],im(:,:,ii))).^2; %faster and reasonably good.
            end
    end
end


text:
function y=fixgaps(x)
% FIXGAPS Linearly interpolates gaps in a time series
% YOUT=FIXGAPS(YIN) linearly interpolates over NaN
% in the input time series (may be complex), but ignores
% trailing and leading NaN.
%

% R. Pawlowicz 6/Nov/99

y=x;

bd=isnan(x);
gd=find(~bd);

bd([1:(min(gd)-1) (max(gd)+1):end])=0;

y(bd)=interp1(gd,x(gd),find(bd));

end

author: matthieukomorowski


text:
function [lImFea] = extr_lIm_fea( lIm )

[nrow, ncol] = size(lIm);

lImFea = zeros([nrow, ncol, 4]);

% first order gradient filters
hf1 = [-1,0,1];
vf1 = [-1,0,1]';

lImFea(:, :, 1) = conv2(lIm, hf1, 'same');
lImFea(:, :, 2) = conv2(lIm, vf1, 'same');

% second order gradient filters
hf2 = [1,0,-2,0,1];
vf2 = [1,0,-2,0,1]';

lImFea(:, :, 3) = conv2(lIm,hf2,'same');
lImFea(:, :, 4) = conv2(lIm,vf2,'same');

author: thomas-koehler

repo: SupER


text:
function check = extreme_values_check ( a, b )

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

  check = 1;

  return
end


text:
function z = cvmgt(x, y, p)
%CVMGT  Conditional merge of two vectors
%
%   Z = CVMGT(X, Y, P) return a vector Z whose elements are X if P is true
%   and Y otherwise.  P, X, and Y should be the same shape except that X
%   and Y may be scalars.  CVMGT stands for conditional vector merge true
%   (an intrinsic function for the Cray fortran compiler).  It implements
%   the C++ statement
%
%     Z = P ? X : Y;

  z = zeros(size(p));
  if isscalar(x)
    z(p) = x;
  else
    z(p) = x(p);
  end
  if isscalar(y)
    z(~p) = y;
  else
    z(~p) = y(~p);
  end
end


text:
function [d] = EpipolarRelationshipCheck(p1,p2,F)

pp2 = p2;
ph2 = e2h(pp2);
l = F'*ph2;
l = l/sqrt(l(1)^2+l(2)^2);

pp1 = p1;
ph1 = e2h(pp1);

d = ph1'*l;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


function pars=extractpars(vars,default);
% function pars=extractpars(vars,default);
%
%

if(nargin<2)
  default=[];
end;

pars=default;
if(length(vars)==1)
 p=vars{1};
 s=fieldnames(p);
 for i=1:length(s)
   eval(['pars.' s{i} '=p.' s{i}]);
 end;

else

 for i=1:2:length(vars)
  if(isstr(vars{i}))
    if(i+1>length(vars)) error(sprintf('Parameter %s has no value\n',vars{i}));else val=vars{i+1};end;
    if(isstr(val))
     eval(['pars.' vars{i} '=''' val ''';']);
    else
     eval(['pars.' vars{i} '=' sprintf('%i',val) ';']);
    end;
  end;
 end;

end;


text:
function bdDistWeight = sc_get_bdDistWeight(mask)

[D,IDX] = bwdist(~mask);

end

repo: StructCompletion


function res = sin(x)
res = 9;
end % fucntion sin


text:
function path = vl_setup(varargin)
% VL_SETUP Add VLFeat Toolbox to the path
%   PATH = VL_SETUP() adds the VLFeat Toolbox to MATLAB path and
%   returns the path PATH to the VLFeat package.
%
%   VL_SETUP('NOPREFIX') adds aliases to each function that do not
%   contain the VL_ prefix. For example, with this option it is
%   possible to use SIFT() instead of VL_SIFT().
%
%   VL_SETUP('TEST') or VL_SETUP('XTEST') adds VLFeat unit test
%   function suite. See also VL_TEST().
%
%   VL_SETUP('QUIET') does not print the greeting message.
%
%   See also: VL_ROOT(), VL_HELP().

% Authors: Andrea Vedaldi and Brian Fulkerson

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

noprefix = false ;
quiet = true ;
xtest = false ;
demo = false ;

for ai=1:length(varargin)
  opt = varargin{ai} ;
  switch lower(opt)
    case {'noprefix', 'usingvl'}
      noprefix = true ;
    case {'test', 'xtest'}
      xtest = true ;
    case {'demo'}
      demo = true ;
    case {'quiet'}
      quiet = true ;
    case {'verbose'}
      quiet = false ;
    otherwise
      error('Unknown option ''%s''.', opt) ;
  end
end

% Do not use vl_root() to avoid conflicts with other VLFeat
% installations.

[a,b,c] = fileparts(mfilename('fullpath')) ;
[a,b,c] = fileparts(a) ;
root = a ;

addpath(fullfile(root,'toolbox'             )) ;
addpath(fullfile(root,'toolbox','aib'       )) ;
addpath(fullfile(root,'toolbox','geometry'  )) ;
addpath(fullfile(root,'toolbox','imop'      )) ;
addpath(fullfile(root,'toolbox','kmeans'    )) ;
addpath(fullfile(root,'toolbox','misc'      )) ;
addpath(fullfile(root,'toolbox','mser'      )) ;
addpath(fullfile(root,'toolbox','plotop'    )) ;
addpath(fullfile(root,'toolbox','quickshift')) ;
addpath(fullfile(root,'toolbox','sift'      )) ;
addpath(fullfile(root,'toolbox','special'   )) ;
addpath(fullfile(root,'toolbox','slic'      )) ;
addpath(fullfile(root,'toolbox','gmm'       )) ;
addpath(fullfile(root,'toolbox','vlad'      )) ;
addpath(fullfile(root,'toolbox','fisher'    )) ;

if vl_isoctave()
  addpath(genpath(fullfile(root,'toolbox','mex','octave'))) ;
  warning('off', 'Octave:possible-matlab-short-circuit-operator') ;
  pkg load image ;
else
  bindir = mexext ;
  if strcmp(bindir, 'dll'), bindir = 'mexw32' ; end
  addpath(fullfile(root,'toolbox','mex',bindir)) ;
end

if noprefix
  addpath(fullfile(root,'toolbox','noprefix')) ;
end

if xtest
  addpath(fullfile(root,'toolbox','xtest')) ;
end

if demo
  addpath(fullfile(root,'toolbox','demo')) ;
end

if ~quiet
  if exist('vl_version') == 3
    fprintf('VLFeat %s ready.\n', vl_version) ;
  else
    warning('VLFeat does not seem to be installed correctly. Make sure that the MEX files are compiled.') ;
  end
end

if nargout == 0
  clear path ;
end


text:
function geometry_test0368 ( )

%*****************************************************************************80
%
%% TEST0368 tests SEGMENT_POINT_DIST_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  test_num = 3;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0368\n' );
  fprintf ( 1, '  SEGMENT_POINT_NEAR_3D computes the nearest point\n' );
  fprintf ( 1, '    from a line segment to a point in 3D.\n' );

  for test = 1 : test_num

    [ p1, seed ] = r8vec_uniform_01 ( dim_num, seed );
    [ p2, seed ] = r8vec_uniform_01 ( dim_num, seed );
    [ p, seed ] = r8vec_uniform_01 ( dim_num, seed );

    [ pn, dist, t ] = segment_point_near_3d ( p1, p2, p );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  TEST = %d', test );
    fprintf ( 1, '  P1 =   %12f  %12f  %12f\n', p1(1:dim_num) );
    fprintf ( 1, '  P2 =   %12f  %12f  %12f\n', p2(1:dim_num) );
    fprintf ( 1, '  P =    %12f  %12f  %12f\n', p(1:dim_num) );
    fprintf ( 1, '  PN =   %12f  %12f  %12f\n', pn(1:dim_num) );
    fprintf ( 1, '  DIST = %12f\n', dist );
    fprintf ( 1, '  T =    %12f\n', t );

  end

  return
end


text:
function mv = mv_conditionsMenu(mv,hfig);
%
% mv = mv_conditionsMenu(mv,hfig);
%
% Add a menu for toggling on/off conditions
% in the multi voxel UI.
%
% ras, 04/05
% ras, 03/06: used to ignore baseline condition, now adds it.
if ieNotDefined('hfig')
    hfig = gcf;
end

if ieNotDefined('mv')
    mv = get(hfig,'UserData');
end

condNums = mv.trials.condNums;
condNames = mv.trials.condNames;
condColors = mv.trials.condColors;

% (callback will be as for the legend callback:
% umtoggle(gcbo);
% timeCourseUI;
mv.ui.condMenu = uimenu('ForegroundColor', 'k', 'Label', 'Conditions', 'Separator', 'on');
accelChars = '0123456789-=';
cbStr=['umtoggle(gcbo); multiVoxelUI;'];
for i = 1:length(condNames)
    if i < length(accelChars)
        accel = accelChars(i);
    else
        accel = '';
    end
    if(isempty(mv.trials.condNames{i}))
        mv.trials.condNames{i}=num2str(i);
    end

    hc(i) = uimenu(mv.ui.condMenu,'Label',condNames{i},...
                     'Separator','off',...
                     'Checked','on',...
                     'Accelerator',accel,...
                     'Tag',num2str(condNums(i)),...
                     'UserData',condColors{i},'Callback',cbStr);
end

set(hc(1), 'Checked', 'off'); % unselect baseline condition

mv.ui.condMenuHandles = hc;

set(hfig,'UserData',mv);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% load some data
load leq
%% plot it
plot(yy);, grid on
title('1-D signal')
%% embed it into a 3-D space with a delay of 10
N=10;
A=[yy(1:end-(N-1)),yy(N/2:end-N/2),yy(N:end)];
%% plot it again
figure
plot3(A(:,1),A(:,2),A(:,3))
grid on, axis tight
title('3-D signal')


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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the LEAF_LEAF-EXAMPLE geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Leaf_Leaf()

%-----------------------------------------------------
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%-----------------------------------------------------
Nx = 384         % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Ny = 64          % # of Eulerian Grid Pts. in y-Direction (MUST BE EVEN!!!)
Lx = 30.0        % Length of Eulerian Grid in x-Direction (cm)
Ly = 5.0         % Length of Eulerian Grid in y-Direction (cm)
ds= 0.5*Lx/Nx;   % Lagrangian spacing

%----------------------------------------------------------------
% AIR PROPERTIES FOR INPUT
%
%     Density:     [rho] = kg/m^3
%     Newton:      [N] = kg * m/s^2
%     Dyn. Visc.:  [N*s/m^2] = (kg*m/s^2)*(s/m^2) = kg / (m*s)
%
%----------------------------------------------------------------
% PROPERTIES IN MKS UNITS
%mu_Air = 18.03e-6; % N*s/m^2 dynamic viscosity of air at 18C (64.4F)
%rho_Air = 1.21;    % kg/m^3  density of air at 18C
mu_Air = 1.895e-5; % N*s/m^2 =  dynamic viscosity of air at 18C (64.4F)
rho_Air = 1.145;    % kg/m^3  density of air at 35C
%(from https://www.engineersedge.com/physics/viscosity_of_air_dynamic_and_kinematic_14483.htm)
%
% CONVERT TO CGS UNITS
mu_Air = mu_Air * (1000/1) * (1/100)     %(1000g/1kg) & (1m/100cm)
rho_Air = rho_Air * (1000/1) * (1/100)^3 %(1000g/1kg) & (1m/100cm)
%
% FINAL UNITS
% [mu_Air] = g/(cm*s)
% [rho_Air] = g/cm^3

%----------------------------------------------------
% LEAF GEOMETRY
%----------------------------------------------------
xLag = [];                     % Initialize Storage
yLag = [];                     % Initialize Storage
L_leaf = 1.0;                  % Length of tail off cylinder
struct_name = 'leaf'; % Name for .vertex, .spring, etc files
%
x0 = 3.0;                      % Starting x-Pt for LEAF
y0 = Ly/2;                     % Starting y_Pt for LEAF
%
xLag_L = x0:ds:x0+L_leaf;
yLag_L = y0*ones(1,length(xLag_L));

%----------------------------------------------------
% IF CHANNEL WALLS / CYLINDRICAL LEADING EDGE %
%----------------------------------------------------
r = 0.01;                      % Radii of Cylinder

%--------------------------------------------
% Call function to construct geometry
%       * No Channel Walls
%       * No Cylinder on Leading Edge
%--------------------------------------------
[xLag_Cy,yLag_Cy] = give_Me_Cylinder_Immersed_Boundary_Geometry(ds,r,x0,y0);
%
% % Plot Geometry to test
% plot(xLag_Ch(1:end/2),yLag_Ch(1:end/2),'r-'); hold on;
% plot(xLag_Ch(end/2+1:end),yLag_Ch(end/2+1:end),'r-'); hold on;
% plot(xLag_Cy,yLag_Cy,'r-'); hold on;
% plot(xLag_L,yLag_L,'m-'); hold on;
% axis([0 Lx 0 Ly]);
% %
% plot(xLag_Ch,yLag_Ch,'*'); hold on;
% plot(xLag_Cy,yLag_Cy,'g*'); hold on;
% xlabel('x'); ylabel('y');
% axis square;
%
% xLag = [xLag xLag_Ch xLag_Cy]; % Add xLagPts from Circle to xLag Pt. Vector (*no springs or beams*)
% yLag = [yLag yLag_Ch yLag_Cy]; % Add xLagPts from Circle to xLag Pt. Vector (*no springs or beams*)
%
% Nbefore = length(xLag);   % # of Lagrangian Pts Before Leaf to COUNT for...
%                           % ... Cylinder / Channel Points have TARGET PTS
%

xLag = [xLag xLag_Cy]; % Add xLagPts from Circle to xLag Pt. Vector (*no springs or beams*)
yLag = [yLag yLag_Cy]; % Add xLagPts from Circle to xLag Pt. Vector (*no springs or beams*)

Nbefore = length(xLag); % # of Lag Pts Before Leaf Starts

xLag = [xLag xLag_L];   % Update xLag Pts to Include Leaf
yLag = [yLag yLag_L];   % Update xLag Pts to Include Leaf

Ntot = length(xLag);    % Total # of Lag Pts

%------------------------
% PLOT GEOMETRY
%------------------------
plot(xLag,yLag,'.'); hold on;
xlabel('x'); ylabel('y');
axis([0 Lx 0 Ly]);

% Prints .vertex file!
print_Lagrangian_Vertices(xLag,yLag,struct_name);

% Prints .spring file!
k_Spring = 2.5e4; %1e7
print_Lagrangian_Springs(xLag,yLag,k_Spring,struct_name,Nbefore,Ntot);

% Prints .beam file!
k_Beam = 2e8; C = 0.0; %2e12 -> 4e12
print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name,Nbefore,Ntot);

% Prints .target file!
k_Target = 2.5e4;
print_Lagrangian_Target_Pts(xLag,k_Target,struct_name,Nbefore);

% Prints .mass file!
k_Mass = 1e0;         % 'spring' stiffness parameter for tethering
Mass = 0.5e-3;        % "MASS" value for 'ghost' lag movement
print_Lagrangian_Mass_Pts(xLag,k_Mass,Mass,struct_name,Nbefore,Ntot);

% Call function for initial concentration
[Concentration]= give_Me_Initial_Concentration(Nx,Ny);

% Prints .concentration file!
print_Concentration_Info(Nx,Ny,Concentration,struct_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called <struct_name>.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag);

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints TARGET points to a file called <struct_name>.target
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name,Nbefore)

    %N = length(xLag);

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', Nbefore+3 );

    %Loops over any Channel/Cylinder Pts AND First 3 Pts Along LEAF
    for s = 1:Nbefore+3
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAMs (Torsional Springs) to a file called <struct_name>.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name,Nbefore,Ntot)

    % k_Beam: beam stiffness
    % C: beam curvature

    if Nbefore ~= 0
        N = Ntot-Nbefore-1;
    else
        N = Ntot-Nbefore-2; % NOTE: Total number of beams
    end

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %------------------------------------------------
    % BEAM BETWEEN CYLINDER AND LEAF
    %       (including those TETHERED in place by
    %                                   TARGET PTS)
    %------------------------------------------------
    if Nbefore ~= 0
        s1 = Nbefore;
        s2 = Nbefore+1;
        s3 = Nbefore+2;
        fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s1, s2, s3, k_Beam, C);
    end

    %------------------------------------------------
    % BEAM BETWEEN VERTICES ALONG "LEAF"
    %       (including those TETHERED in place by
    %                                   TARGET PTS)
    %------------------------------------------------
    for s = Nbefore+2:Ntot-1
        s1 = s-1;
        s2 = s;
        s3 = s+1;
        fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s1, s2, s3, k_Beam, C);
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called <struct_name>.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,yLag,k_Spring,struct_name,Nbefore,Ntot)

    if Nbefore ~= 0
        N = Ntot-Nbefore;
    else
        N = Ntot-Nbefore-1;
    end

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %------------------------------------------------
    % SPRINGS BETWEEN CYLINDER AND LEAF
    %       (including those TETHERED in place by
    %                                   TARGET PTS)
    %------------------------------------------------
    if Nbefore ~= 0

        % springs between adjacent pts on cylinder
        for i=1:Nbefore
            if i<Nbefore
                x1 = xLag(i); y1 = yLag(i);
                x2 = xLag(i+1); y2 = yLag(i+1);
                ds = sqrt( (x1-x2)^2 + (y1-y2)^2 );
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', i, i+1, k_Spring, ds);
            else
                x1 = xLag(i); y1 = yLag(i);
                x2 = xLag(1); y2 = yLag(1);
                ds = sqrt( (x1-x2)^2 + (y1-y2)^2 );
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', i, 1, k_Spring, ds);

            end
        end

        % springs across opposite sides of cylinder
        for i=1:floor(Nbefore/2)
            x1 = xLag(i); y1 = yLag(i);
            x2 = xLag( i+floor(Nbefore/2) ); y2 = yLag( i+floor(Nbefore/2) );
            ds = sqrt( (x1-x2)^2 + (y1-y2)^2 );
            fprintf(spring_fid, '%d %d %1.16e %1.16e\n', i, i+floor(Nbefore/2), k_Spring, ds);
        end

        % spring between cylinder and leaf
        x1 = xLag(Nbefore); y1 = yLag(Nbefore);
        x2 = xLag(Nbefore+1); y2 = yLag(Nbefore+1);
        ds = sqrt( (x1-x2)^2 + (y1-y2)^2 );
        fprintf(spring_fid, '%d %d %1.16e %1.16e\n', Nbefore, Nbefore+1, k_Spring, ds);

    end

    %------------------------------------------------
    % SPRINGS BETWEEN VERTICES ALONG LEAF
    %       (including those TETHERED in place by
    %                                   TARGET PTS)
    %------------------------------------------------
    for s = Nbefore+1:Ntot-1

        x1 = xLag(s); y1 = yLag(s);
        x2 = xLag(s+1); y2 = yLag(s+1);

        ds = sqrt( (x1-x2)^2 + (y1-y2)^2 );

        fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds);

    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints TARGET points to a file called <struct_name>.mass
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Mass_Pts(xLag,kMass,Mass,struct_name,Nbefore,Ntot)

    %LOOP OVER LAG PTS FOR MASS PTS IN LAGRANGIAN INDEXING
    N = (Ntot-Nbefore);

    mass_fid = fopen([struct_name '.mass'], 'w');

    fprintf(mass_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = Nbefore+1:Ntot%Nbefore:Ntot
        fprintf(mass_fid, '%d %1.16e %1.16e\n', s, kMass, Mass );
    end

    fclose(mass_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry for channel
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Channel_Immersed_Boundary_Geometry(ds,L,w,Lx,Ly)

% The immsersed structure is a channel %
x = (Lx-L)/2:ds:(L+(Lx-L)/2);  %xPts
yBot = (Ly-w)/2;               %yVal for bottom of Channel
yTop = Ly - (Ly-w)/2;          %yVal for top of Channel

xLag = [x x];
yLag = [yBot*ones(1,length(x)) yTop*ones(1,length(x))];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry for cylinder
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLagN,yLagN] = give_Me_Cylinder_Immersed_Boundary_Geometry(ds,r,x0,y0)

% The immsersed structure is a cylinder %

dtheta = ds/ (4*r);
theta = 0; i=1;
while theta < 2*pi
   xLag(i) = (x0-r-ds) + r*cos(theta);
   yLag(i) = y0 + r*sin(theta);
   theta = theta + dtheta;
   i=i+1;
end

for i=1:length(xLag)
   xLagN(i) = xLag(end-(i-1));
   yLagN(i) = yLag(end-(i-1));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints initial concentration to a file called
% <struct_name>.concentration
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Concentration_Info(Nx,Ny,C,struct_name)

    con_fid = fopen([struct_name '.concentration'], 'w');

    for i=1:Ny
        for j=1:Nx
            fprintf(con_fid, '%1.16e ', C(i,j) );
        end
        fprintf(con_fid,'\n');
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: gives initial concentration gradient inside channel
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [C] = give_Me_Initial_Concentration(Nx,Ny)

C = zeros(Ny,Nx);

author: nickabattista


text:
function line_fekete_rule_test03 ( m )

%*****************************************************************************80
%
%% LINE_FEKETE_RULE_TEST03 seeks Fekete points in [-1,+1].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the dimension of the polynomial space.
%
  n = 5001;
  a = -1.0;
  b = +1.0;
  if ( nargin < 1 )
    m = 5;
  end
  x = linspace ( a, b, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_FEKETE_RULE_TEST03:\n' );
  fprintf ( 1, '  Seek Fekete points in [%g,%g]\n', a, b );
  fprintf ( 1, '  using %d equally spaced sample points\n', n );
  fprintf ( 1, '  for polynomials of degree M = %d\n', m );
  fprintf ( 1, '  with the Legendre basis and uniform weight.\n' );

  [ nf, xf, wf, vf ] = line_fekete_legendre ( m, a, b, n, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  NF = %d\n', nf );
  r8vec_print ( nf, xf, '  Estimated Fekete points XF:' );

  wf_sum = sum ( wf(1:nf) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sum(WF) = %g\n', wf_sum );

  return
end


text:
function varargout = transformPoint3d(pts, transfo, varargin)
%TRANSFORMPOINT3D Transform a point with a 3D affine transform.
%
%   PT2 = transformPoint3d(PT1, TRANS);
%   PT2 = transformPoint3d(X1, Y1, Z1, TRANS);
%   where PT1 has the form [xp yp zp], and TRANS is a 3-by-3, 3-by-4, or
%   4-by-4 matrix, returns the point transformed according to the affine
%   transform specified by TRANS.
%
%   The function accepts transforms given using the following formats:
%   [a b c]   ,   [a b c j] , or [a b c j]
%   [d e f]       [d e f k]      [d e f k]
%   [g h i]       [g h i l]      [g h i l]
%                                [0 0 0 1]
%
%   PT2 = transformPoint3d(PT1, TRANS)
%   also work when PT1 is a N-by-3-by-M-by-P-by-ETC array of double. In
%   this case, PT2 has the same size as PT1.
%
%   PT2 = transformPoint3d(X1, Y1, Z1, TRANS);
%   also work when X1, Y1 and Z1 are 3 arrays with the same size. In this
%   case, PT2 will be a 1-by-3 cell containing {X Y Z} outputs of size(X1).
%
%   [X2, Y2, Z2] = transformPoint3d(...);
%   returns the result in 3 different arrays the same size as the input.
%   This form can be useful when used with functions like meshgrid or warp.
%
%   MESH2 = transformPoint3d(MESH, TRANS)
%   transforms the field 'vertices' of the struct MESH and returns the same
%   struct with the transformed vertices.
%   (It is recommended to use the function 'transformMesh', within the
%   "meshes3d" module).
%
%   See also
%     points3d, transforms3d, transformMesh, createTranslation3d
%     createRotationOx, createRotationOy, createRotationOz, createScaling
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2005-02-10
% Copyright 2005-2022 INRA - TPV URPOI - BIA IMASTE

%% Parse input arguments

% Check special case: if first argument is a struct with a field named
% 'vertices', then the output will be the same struct, but with the
% transformed vertices.
if nargin == 2 && isstruct(pts) && isfield(pts, 'vertices')
    mesh = pts;
    mesh.vertices = transformPoint3d(mesh.vertices, transfo);
    varargout = {mesh};
    return;
end

% Parse x, y, and z coordinates of input points from input arguments
if nargin == 2
    % Point coordinates are given in a single N-by-3-by-M-by-etc argument.
    % Preallocate x, y, and z to size N-by-1-by-M-by-etc, then fill them in
    dim = size(pts);
    dim(2) = 1;
    [x, y, z] = deal(zeros(dim, class(pts)));
    x(:) = pts(:,1,:);
    y(:) = pts(:,2,:);
    z(:) = pts(:,3,:);

elseif nargin == 4
    % Point coordinates are given in 3 different arrays
    x = pts;
    y = transfo;
    z = varargin{1};
    transfo = varargin{2};
    dim = size(x);

else
    error('MatGeom:geom3d:WrongInputArgumentNumber', ...
        'Requires number of input arguments to be either 2 or 4');
end

%% Process transformation matrix

% extract the linear and the translation parts of the matrix
linear = transfo(1:3, 1:3)';
trans = [0 0 0];
if size(transfo, 2) > 3
    trans = transfo(1:3, 4)';
end

%% Main processing

% convert coordinates
try
    % vectorial processing, if there is enough memory.
    % same as:
    % res = (transfo * [x(:) y(:) z(:) ones(NP, 1)]')';
    res = bsxfun(@plus, [x(:) y(:) z(:)] * linear, trans);

    % Back-fill x,y,z with new result (saves calling costly reshape())
    x(:) = res(:,1);
    y(:) = res(:,2);
    z(:) = res(:,3);

catch ME
    disp(ME.message)
    % process each point one by one, writing in existing array
    NP = numel(x);
    for i = 1:NP
        res = [x(i) y(i) z(i)] * linear + trans;
        x(i) = res(1);
        y(i) = res(2);
        z(i) = res(3);
    end
end

% process output arguments
if nargout <= 1
    % results are stored in a unique array
    if length(dim) > 2 && dim(2) > 1
        warning('geom3d:shapeMismatch',...
            'Shape mismatch: Non-vector xyz input should have multiple x,y,z output arguments. Cell {x,y,z} returned instead.')
        varargout{1} = {x,y,z};
    else
        varargout{1} = [x y z];
    end

elseif nargout == 3
    % results are returned in three array with same size.
    varargout = {x, y, z};
end


text:
function r8poly_print ( n, a, title )

%*****************************************************************************80
%
%% R8POLY_PRINT prints out a polynomial.
%
%  Discussion:
%
%    The power sum form is:
%
%      p(x) = a(0) + a(1) * x + ... + a(n-1) * x^(n-1) + a(n) * x^(n)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of A.
%
%    Input, real A(1:N+1), the polynomial coefficients.
%    A(1) is the constant term and
%    A(N+1) is the coefficient of X^N.
%
%    Input, character TITLE(*), an optional title.
%
  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  fprintf ( 1, '\n' );

  n = r8poly_degree ( n, a );

  if ( a(n+1) < 0.0 )
    plus_minus = '-';
  else
    plus_minus = ' ';
  end

  mag = abs ( a(n+1) );

  if ( 2 <= n )
    fprintf ( 1, '  p(x) = %c%14f * x^%d\n', plus_minus, mag, n );
  elseif ( n == 1 )
    fprintf ( 1, '  p(x) = %c%14f * x\n', plus_minus, mag );
  elseif ( n == 0 )
    fprintf ( 1, '  p(x) = %c%14f\n', plus_minus, mag );
  end

  for i = n-1 : -1 : 0

    if ( a(i+1) < 0.0 )
      plus_minus = '-';
    else
      plus_minus = '+';
    end

    mag = abs ( a(i+1) );

    if ( mag ~= 0.0E+00 )

      if ( 2 <= i )
        fprintf ( 1, '         %c%14f * x^%d\n', plus_minus, mag, i );
      elseif ( i == 1 )
        fprintf ( 1, '         %c%14f * x\n', plus_minus, mag );
      elseif ( i == 0 )
        fprintf ( 1, '         %c%14f\n', plus_minus, mag );
      end
    end

  end

  return
end


text:
function x = block_levinson(y, L)
%BLOCK_LEVINSON Block Levinson recursion for efficiently solving
%symmetric block Toeplitz matrix equations.
%   BLOCK_LEVINSON(Y, L) solves the matrix equation T * x = y, where T
%   is a symmetric matrix with block Toeplitz structure, and returns the
%   solution vector x. The matrix T is never stored in full (because it
%   is large and mostly redundant), so the input parameter L is actually
%   the leftmost "block column" of T (the leftmost d columns where d is
%   the block dimension).

%   Author: Keenan Pepper
%   Last modified: 2007-12-23

%   References:
%     [1] Akaike, Hirotugu (1973). "Block Toeplitz Matrix Inversion".
%     SIAM J. Appl. Math. 24 (2): 234-241

s = size(L);
d = s(2);                 % Block dimension
N = s(1) / d;             % Number of blocks

B = reshape(L, [d,N,d]);  % This is just to get the bottom block row B
B = permute(B, [1,3,2]);  % from the left block column L
B = flipdim(B, 3);
B = reshape(B, [d,N*d]);

f = L(1:d,:)^-1;          % "Forward" block vector
b = f;                    % "Backward" block vector
x = f * y(1:d);           % Solution vector

for n = 2:N
    ef = B(:,(N-n)*d+1:N*d) * [f;zeros(d)];
    eb = L(1:n*d,:)' * [zeros(d);b];
    ex = B(:,(N-n)*d+1:N*d) * [x;zeros(d,1)];
    A = [eye(d),eb;ef,eye(d)]^-1;
    fn = [[f;zeros(d)],[zeros(d);b]] * A(:,1:d);
    bn = [[f;zeros(d)],[zeros(d);b]] * A(:,d+1:end);
    f = fn;
    b = bn;
    x = [x;zeros(d,1)] + b * (y((n-1)*d+1:n*d) - ex);
end

author: Sable