text:
 function [Gauss,terms,logp,comptime] = laplacedegenerate_ep(y,X,K,varargin)

% specialized version of our fast variant of
% Expectation Propagation for Logistic Regression (2 classes) with a (correlated) Laplace prior
% in the degenerate case, i.e., when the number of features >> number of samples
% (also works when number of samples > number of features, but then less stable!!!)
%
% input:
% labels    = an N x 1 vector of class labels [1,2]
% examples  = an N x M matrix of input data
% K         = the prior precision matrix of size M x M
%
% Note: the bias term should be explicitly added to K and the examples!
%
% regression parameters refers to betas
% auxiliary variables refers to u and v whose precision matrix auxK = inv(Lambda) couples the features.
%
% we have a precision matrix of the form
%
% | K_beta         |
% |        K_u     |
% |            K_v |
%
% priorGauss: struct with fields
%
%     hatK       diagonal of precision matrix (number of samples x 1);
%                   (initially zero)
%     diagK      diagonal of precision matrix (number of features x 1)
%                   (initially zero)
%
%     precision matrix of regression parameters K_beta = A' hatK A + diagK
%
%     h          canonical mean (number of features x 1)
%                   (initially zero)
%
%     auxK       precision matrix of auxiliary variables (number of features x number of features; sparse)
%                   (contains the covariance structure of interest)
%
%     A          feature matrix (number of samples x number of features)
%
% terms: struct with fields
%
%     hatK       as priorGauss
%     diagK      as priorGauss
%     hath       canonical mean (number of samples x 1)
%     h          canonical mean (number of features x 1)
%                                canonical mean of regression parameters = h + A' hath
%     auxK       as priorGauss, but then only diagonal elements (number of features x 1)
%
% opt: struct with fields (all optional, defaults in brackets)
%
%     maxstepsize   maximum step size [1]
%     fraction      fraction or power for fractional/power EP [1]
%     niter         maximum number of iterations [100]
%     tol           convergence criterion [1e-5]
%     nweights      number of points for numerical integration [20]
%     temperature   temperature for simulated annealing [1]
%     verbose       give output during run [1]
%
% Gauss: struct with fields as in priorGauss plus
%
%     hatB       diagonal of projected covariance (number of samples x 1)
%                                                               hatB = A Covariance A'
%     hatn       projected mean (number of samples x 1)
%                                                               hatn = A m
%     diagC      diagonal of covariance matrix of regression parameters (number of features x 1)
%     auxC       diagonal of covariance matrix of auxiliary variables (number of features x 1)
%
%
% logp: estimated marginal loglikelihood (implementation doubtful...)
%
% comptime: computation time

%% initialization

% parse opt
opt = [];
for i=1:2:length(varargin)
  opt.(varargin{i}) = varargin{i+1};
end

if ~isfield(opt,'maxstepsize'), opt.maxstepsize = 1;    end
if ~isfield(opt,'fraction'),    opt.fraction    = 0.95; end
if ~isfield(opt,'niter'),       opt.niter       = 1000; end
if ~isfield(opt,'tol'),         opt.tol         = 1e-5; end
if ~isfield(opt,'nweights'),    opt.nweights    = 50;   end
if ~isfield(opt,'temperature'), opt.temperature = 1;    end
if ~isfield(opt,'lambda'),      opt.lambda      = 1;    end
if ~isfield(opt,'verbose'),     opt.verbose      = 1;    end

if opt.verbose
  fprintf('starting EP\n');
end

tic

[nsamples,nfeatures] = size(X);

%% create priorGauss and terms

A = X.*repmat(y,1,nfeatures);

% construct Gaussian representation
priorGauss.A     = A;
priorGauss.hatK  = zeros(nsamples,1);
priorGauss.h     = zeros(nfeatures,1);
priorGauss.diagK = zeros(nfeatures,1);
priorGauss.auxK  = K;

% compute additional terms for the EP free energy
[cholK,dummy,S]  = chol(K,'lower');
LogPriorRestAux2 = (2*sum(log(full(diag(cholK))))); % -nfeatures*log(2*pi)); % redundant MvG

% construct term representation
terms.hatK  = ones(nsamples,1)/10;
terms.hath  = zeros(nsamples,1);
terms.diagK = ones(nfeatures,1)./(10*opt.lambda);
terms.auxK  = zeros(nfeatures,1);
terms.h     = zeros(nfeatures,1);

%% precompute points and weights for numerical integration

[xhermite,whermite] = gausshermite(opt.nweights);
xhermite = xhermite(:);       % nhermite x 1
whermite = whermite(:);       % nhermite x 1

[xlaguerre,wlaguerre] = gausslaguerre(opt.nweights);
xlaguerre = xlaguerre(:);     % nlaguerre x 1
wlaguerre = wlaguerre(:);     % nlaguerre x 1

% divide all canonical parameters by the temperature
if opt.temperature ~= 1,
   [priorGauss,terms] = correct_for_temperature(priorGauss,terms,opt.temperature);
end

%% build initial Gauss

% add terms to prior
Gauss        = update_Gauss(priorGauss,terms);
Gauss        = canonical_to_moments(Gauss);
[myGauss,ok] = project_all(Gauss,terms,opt.fraction);
if ~ok,
   error('improper cavity distributions\n');
end

prior = Gauss; % save prior Gauss

%% enter the iterations

logp     = 0;
logpold  = 2*opt.tol;
change   = 0;
teller   = 0;
stepsize = opt.maxstepsize;

while abs(logp-logpold) > stepsize*opt.tol && teller < opt.niter,

   teller    = teller+1;
   logpold   = logp;
   oldchange = change;

   % compute the new term approximation by applying an adf update on all the cavity approximations
   [fullterms,logreglogz,crosslogz] = ...
      adfupdate_all(myGauss,terms,xhermite,whermite,xlaguerre,wlaguerre,opt.fraction,opt.temperature);

   ok  = 0;
   ok1 = 1;
   ok2 = 1;

   while ~ok,

      % try to replace the old term approximations by the new term approximations and check whether the new Gauss is still fine

     [newGauss,newterms,ok1] = try_update(Gauss,fullterms,terms,stepsize);
     %[newGauss,newterms,ok1] = try_update(Gauss,fullterms,terms,0.5);

      if ok1,
         % compute all cavity approximations needed for the next EP updates and check whether they are all fine
         [newGauss,myGauss,logZappx,ok2] = try_project(newGauss,newterms,opt.fraction);
      end

      ok = (ok1 & ok2);
      if ok,    % accept

        terms = newterms;
        Gauss = newGauss;
        stepsize = min(opt.maxstepsize,stepsize*1.9);

      else  % try with smaller stepsize

        stepsize = stepsize/2;
         if ~ok1,
            fprintf('improper full covariance: lowering stepsize to %g\n',stepsize');
         elseif ~ok2,
            fprintf('improper cavity covariance: lowering stepsize to %g\n',stepsize');
         end
         if stepsize < 1e-10,
            warning('Cannot find an update that leads to proper cavity approximations');
            teller = opt.niter;
            break;
         end

      end

   end

   % compute marginal moments

   if ok,

      % compute marginal loglikelihood

      %logp = sum(logz) + sum(crosslogz) + logdet/2 + contribQ;
%      CorrTermCross1 = (myGauss.m.^2)./myGauss.diagC - (newGauss.m.^2)./newGauss.diagC + log(myGauss.diagC/newGauss.diagC);
%	   CorrTermCross2 = log(myGauss.auxC/newGauss.auxC);
%      CorrTerm  = 0.5*CorrTermCross1 + 0.5*(CorrTermCross2 + CorrTermCross2);

      logp      = (sum(crosslogz) + sum(logreglogz))./opt.fraction + logZappx + LogPriorRestAux2;

      if opt.verbose
        fprintf('%d: %g (stepsize: %g)\n',teller,logp,stepsize);
      end

      % check whether marginal loglikelihood is going up and down, if so, lower stepsize

      change = logp-logpold;
      if change*oldchange < 0,   % possibly cycling
         stepsize = stepsize/2;
      end
      oldchange = change;
   end

end

comptime = toc;

if opt.verbose
  fprintf('EP finished in %s seconds\n',num2str(comptime));
end

Gauss.prior = prior;

%%% END MAIN

%%%%%%%%%
%
% compute the cavity approximations that result when subtracting a fraction of the term approximations

function [myGauss,ok] = project_all(Gauss,terms,fraction)

if nargin < 3,
   fraction = 1;
end

% take out and project in moment form

% (1) regression parameters

[myGauss.hatB,myGauss.hatn] = ...
   rank_one_update(Gauss.hatB,Gauss.hatn,-fraction*terms.hatK,-fraction*terms.hath);

% (2) cross terms between regression parameters and auxiliary parameters

[myGauss.diagC,myGauss.m] = ...
   rank_one_update(Gauss.diagC,Gauss.m,-fraction*terms.diagK,-fraction*terms.h);

 myGauss.auxC = rank_one_update(Gauss.auxC,[],-fraction*terms.auxK);

% check whether all precision matrices are strictly positive definite

if nargout > 1,
   ok = (all(myGauss.hatB > 0) & all(myGauss.diagC > 0) & all(myGauss.auxC > 0));
end

%%%%%%%%%
%
% compute the new term approximation by applying an adf update on all the cavity approximations

function [fullterms,logreglogz,crosslogz] = ...
   adfupdate_all(myGauss,terms,xhermite,whermite,xlaguerre,wlaguerre,fraction,temperature)

if nargin < 8,
   temperature = 1;
end
if nargin < 7,
   fraction = 1;
end

%% (1) regression parameters

oldm     = myGauss.hatn;
oldC     = myGauss.hatB;
sqrtC    = sqrt(oldC);
nsamples = length(oldm);
nhermite = length(whermite);

% translate and scale the sample points to get the correct mean and variance

x = repmat(oldm,1,nhermite) + sqrtC*xhermite';

% compute the terms at the sample points

g = logist(x);   % returns - log (1 + exp(-x)) with special attention for very small and very large x

% correct for fraction and temperature and incorporate the sample weights

g           = fraction*g/temperature + log(repmat(whermite',nsamples,1));
maxg        = max(g,[],2);
g           = g-repmat(maxg,1,nhermite);
expg        = exp(g);
denominator = sum(expg,2);
neww        = expg./repmat(denominator,1,nhermite);

% compute the moments

Ex   = sum(x.*neww,2);
Exx  = sum(x.^2.*neww,2);
newm = Ex;
newC = Exx-Ex.^2;

% derive the term approximation from the change in mean and variance

[fullterms.hatK,fullterms.hath,logzextra] = compute_termproxy(newC,newm,oldC,oldm,fraction);

% contributions to marginal loglikelihood
logreglogz = maxg + log(denominator) + logzextra;

%% (2) cross terms between regression parameters and auxiliary variables

oldm      = myGauss.m;
oldC      = myGauss.diagC;
oldlambda = myGauss.auxC;
nfeatures = length(oldm);
nlaguerre = length(wlaguerre);

% this part heavily relies on the accompanying note
% basic idea:
% - the cavity approximation on U is an exponential distribution
% - we have analytical formulas for the moments of x conditioned upon U
% - marginal moments can then be computed through numerical integration with Gauss-Laguerre

% translate and scale the sample points to get the correct mean

U  = 2*oldlambda*xlaguerre';     % nfeatures x nlaguerre

mm = repmat(oldm,1,nlaguerre);
CC = repmat(oldC,1,nlaguerre);

% compute the partition function (integral over x) given U and turn this into weights required for computing the marginal moments

g           = (1-fraction)*log(U)/2 - fraction*mm.^2./(U + fraction*CC)/2 - log(U+ fraction*CC)/2;

g           = bsxfun(@plus,g,log(wlaguerre'));
maxg        = max(g,[],2);
g           = bsxfun(@minus,g,maxg);
expg        = exp(g);

denominator = sum(expg,2);

neww        = bsxfun(@rdivide,expg,denominator);

% compute the marginal moments through numerical integration

ExgU = mm.*U./(U + fraction*CC);
Ex   = sum(ExgU.*neww,2);

ExxgU = ExgU.^2 + CC.*U./(U+fraction*CC);
Exx   = sum(ExxgU.*neww,2);
EU    = sum(U.*neww,2);

newm      = Ex;
newC      = Exx-Ex.^2;
newlambda = EU/2;

% derive the term approximation from the change in mean and variance

[fullterms.diagK,fullterms.h,logzextra1] = compute_termproxy(newC,newm,oldC,oldm,fraction);
[fullterms.auxK,dummy,logzextra2] = compute_termproxy(newlambda,zeros(nfeatures,1),oldlambda,zeros(nfeatures,1),fraction);

crosslogz = maxg + log(denominator) + logzextra1 + 2*logzextra2;
% multiplied the last term by two

%%%%%%%%%%
%
% compute the moments corresponding to the canonical parameters

function [Gauss,logp] = canonical_to_moments(Gauss)

[nsamples,nfeatures] = size(Gauss.A);

%% (1) regression parameters

if nsamples > nfeatures,   % in the non-degenerate case, this direct route is more stable and faster

   scaledA     = Gauss.A.*(repmat(Gauss.hatK,1,nfeatures));
   K           = Gauss.A'*scaledA + diag(Gauss.diagK);
   [C,logdet1] = invert_chol(K);
   Gauss.m     = C*Gauss.h;
   Gauss.hatB  = zeros(nsamples,1);          % only need diagonal
   for k=1:nsamples,
      Gauss.hatB(k) = Gauss.A(k,:)*C*Gauss.A(k,:)';
   end
   Gauss.diagC = diag(C);
   Gauss.hatn  = Gauss.A*Gauss.m;

else

   % this part heavily relies on the appendix of the accompanying note
   % basic idea:
   % - the precision matrix K is of the form A' hatK A + diagK, where both hatK and diagK are diagonal matrices
   % - apply Woodbury's formula to replace inverting an (nfeat x nfeat) matrix by an (nsample x nsample) alternative
   % - projections of the covariance matrix and the mean onto the feature matrix then follow immediately

   scaledA     = bsxfun(@rdivide,Gauss.A,Gauss.diagK');
   W           = Gauss.A*scaledA';
   W           = (W + W')/2;    % make symmetric
   [Q,logdet1] = invert_chol(diag(1./Gauss.hatK) + W);

   Gauss.hatB = zeros(nsamples,1);
   for k=1:nsamples,
      Gauss.hatB(k) = W(k,k) - W(k,:)*Q*W(:,k);
   end

   Gauss.m     = Gauss.h./Gauss.diagK - scaledA'*(Q*(scaledA*Gauss.h));
   Gauss.hatn  = Gauss.A*Gauss.m;

   Gauss.diagC = 1./Gauss.diagK;

%     for i=1:nfeatures,
%        Gauss.diagC(i) = Gauss.diagC(i) - scaledA(:,i)'*Q*scaledA(:,i);
%     end

   % adriana's recipe
   z = scaledA' * Q; for i=1:size(z), Gauss.diagC(i) = Gauss.diagC(i) - z(i,:) * scaledA(:,i); end

   logdet1 = logdet1 + sum(log(Gauss.diagK)) + sum(log(Gauss.hatK));

end

% compute quadratic term (BC)

qterm = sum(Gauss.m .* Gauss.diagK .* Gauss.m); % = m' * diagK * m
qterm = qterm + sum(Gauss.hatn .* Gauss.hatK .* Gauss.hatn);

logp1 = 0.5*(qterm - logdet1);

%% (2) auxiliary variables; i.e., wrt scale mixture representation of
% Laplace prior

% this is (by far) the most expensive step when nsamples << nfeatures
% and the precision matrix of the auxiliary variables is non-diagonal

[auxC,logdet2] = invert_chol(Gauss.auxK); % only need diagonal terms
Gauss.auxC     = full(diag(auxC));        % turn into full vector
logp2          = 0.5*( - logdet2);
logp           = logp1 + 2*logp2;

%%%%%%%%%%
%
% take out the old term proxies and add the new termproxies and check whether the resulting Gaussian is still normalizable

function [newGauss,newterms,ok] = try_update(Gauss,fullterms,terms,stepsize)

if nargin < 4,
   stepsize = 1;
end

% take out the old term proxies

newGauss = update_Gauss(Gauss,terms,-1);

% compute the new term proxies as a weighted combi of the old ones and the "full" (stepsize 1) term proxies

newterms = combine_terms(fullterms,terms,stepsize);

% add the new term proxies

newGauss = update_Gauss(newGauss,newterms,1);

[L,check,dummy] = chol(newGauss.auxK,'lower');   % check whether full covariance matrix is ok
                % note that this is bit inefficient, since we redo the Cholesky later when everything is fine

ok = (check == 0 & all(newGauss.hatK > 0) & all(newGauss.diagK > 0));  % perhaps a bit too strong???

%%%%%%%%%%%%
%
% compute the moment form of the current Gauss and all cavity approximations and check whether they are fine

function [Gauss,myGauss,logdet,ok] = try_project(Gauss,terms,fraction)

if nargin < 3,
   fraction = 1;
end

[Gauss,logdet] = canonical_to_moments(Gauss);
[myGauss,ok]   = project_all(Gauss,terms,fraction);

%%%%%%%%%%
%
% if we use a temperature < 1, to get closer to the MAP solution, we have to change the prior and initial term proxies accordingly

function [Gauss,terms] = correct_for_temperature(Gauss,terms,temperature)

% note: choose temperature small to implement MAP-like behavior

Gauss.hatK  = Gauss.hatK/temperature;
Gauss.h     = Gauss.h/temperature;
Gauss.auxK  = Gauss.auxK/temperature;
Gauss.diagK = Gauss.diagK/temperature;

terms.hatK  = terms.hatK/temperature;
terms.hath  = terms.hath/temperature;
terms.diagK = terms.diagK/temperature;
terms.auxK  = terms.auxK/temperature;
terms.h     = terms.h/temperature;

%%%%%%%%%%
%
% invert a positive definite matrix using Cholesky factorization

function [invA,logdet] = invert_chol(A)

if issparse(A)

  if 0 % matlab version; slower but useful in case of mex problems

    [L,dummy,S] = chol(sparse(A),'lower');   % now A = S*(L*L')*S' and (L*L') = S'*A*S

    n           = length(L);
    invdiagL2   = 1./spdiags(L,0).^2;

    invA = A;
    for i=n:-1:1,
      I      = i+find(L(i+1:n,i));
      invA(I,i) = -(invA(I,I)*L(I,i))/L(i,i);
      invA(i,I) = invA(I,i)';
      invA(i,i) = invdiagL2(i) - (invA(i,I)*L(I,i))/L(i,i);
    end

    invA = S*invA*S';

  else

    [L,dummy,S] = chol(sparse(A),'lower');   % now A = S*(L*L')*S' and (L*L') = S'*A*S

    if dummy
      error('matrix is not p.d.');
    end

    invA = fastinvc(L);
    invA = S*invA*S';
  end

else

  [L,dummy] = chol(A,'lower');

  if dummy
    error('matrix is not p.d.');
  end

  invA = inv(A);

end

if nargout > 1,
  logdet = 2*sum(log(full(diag(L))));
end

%%%%%%%%%%
%
% compute the term proxy when [oldC,oldm] changes to [newC,newm]

function [K,h,logz] = compute_termproxy(newC,newm,oldC,oldm,fraction)

if nargin < 5,
   fraction = 1;
end

K    = (1./newC - 1./oldC)/fraction;
h    = (newm./newC - oldm./oldC)/fraction;
logz = oldm.^2./oldC/2 - newm.^2./newC/2 + log(full(oldC./newC))/2 ;

%%%%%%%%%%
%
% Sherman-Morrison formula to compute the change from [oldC,oldm] to [newC,newm] when we add [K,h] to the corresponding canonical parameters

function [newC,newm] = rank_one_update(oldC,oldm,K,h)

dummy         = K.*oldC;
oneminusdelta = 1./(1+dummy);
newC          = oneminusdelta.*oldC;

if nargout > 1,
   newm = oneminusdelta.*(oldm + h.*oldC);
end

%%%%%%%%%%%
%
% general procedure for a weighted combi of the fields of two structures

function terms = combine_terms(terms1,terms2,stepsize)

names1 = fieldnames(terms1);
names2 = fieldnames(terms2);
names  = intersect(names1,names2);

terms = struct;
for i=1:length(names)
  terms.(names{i}) = stepsize*terms1.(names{i}) + (1-stepsize)*terms2.(names{i});
end

%%%%%%%%%%%
%
% updates the Gaussian representation with new term proxies

function Gauss = update_Gauss(Gauss,terms,const)

if nargin < 3,
   const = 1;
end

Gauss.h     = Gauss.h     + const*Gauss.A'*terms.hath + const*terms.h;
Gauss.hatK  = Gauss.hatK  + const*terms.hatK;
Gauss.diagK = Gauss.diagK + const*terms.diagK;

% get diagonal elements
diagidx             = 1:(size(Gauss.auxK,1)+1):numel(Gauss.auxK);
Gauss.auxK(diagidx) = Gauss.auxK(diagidx) + const*terms.auxK';

author: fieldtrip

repo: fieldtrip


text:
%% Copyright (C) 2014-2020, 2022-2023 Colin B. Macdonald
%% Copyright (C) 2017 NVS Abhilash
%% Copyright (C) 2017 Mike Miller
%% Copyright (C) 2022 Alex Vong
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deftypefn  Command {} sympref @var{cmd}
%% @deftypefnx Command {} sympref @var{cmd} @var{args}
%% @deftypefnx Function {@var{r} =} sympref ()
%% @deftypefnx Function {@var{r} =} sympref (@var{cmd})
%% @deftypefnx Function {@var{r} =} sympref (@var{cmd}, @var{args})
%% Preferences for the Symbolic package.
%%
%% @code{sympref} can set or get various preferences and
%% configurations.  The various choices for @var{cmd} and
%% @var{args} are documented below.
%%
%% Run @strong{diagnostics} on your system:
%% @example
%% @comment doctest: +SKIP
%% sympref diagnose
%%   @print{} ...
%% @end example
%%
%%
%% @strong{Display} of syms:
%%
%% @example
%% @group
%% sympref display
%%   @result{} ans = unicode
%%
%% @end group
%% @group
%% syms x
%% sympref display flat
%% sin(x/2)
%%   @result{} (sym) sin(x/2)
%%
%% sympref display ascii
%% sin(x/2)
%%   @result{} (sym)
%%          /x\
%%       sin|-|
%%          \2/
%%
%% sympref display unicode
%% sin(x/2)
%%   @result{} (sym)
%%          ⎛x⎞
%%       sin⎜─⎟
%%          ⎝2⎠
%%
%% sympref display default
%% @end group
%% @end example
%%
%% By default, a unicode pretty printer is used to display
%% symbolic expressions.  If that doesn't work (e.g., if you
%% see @code{?} characters) then try the @code{ascii} option.
%%
%%
%% @strong{Communication mechanism}:
%%
%% @example
%% @group
%% sympref ipc
%%   @result{} ans = default
%% @end group
%% @end example
%%
%% This default depends on your system.  If you have loaded the
%% @uref{https://gitlab.com/mtmiller/octave-pythonic, Pythonic package},
%% the default will be the @code{native} mechanism.
%% Otherwise, typically the @code{popen2} mechanism will be used,
%% which uses a pipe to communicate with Python.
%% If that doesn't work, try @code{sympref ipc system} which is
%% much slower, as a new Python process is started for each operation.
%%
%% Other options for @code{sympref ipc} include:
%% @itemize
%% @item @code{sympref ipc popen2}: force popen2 choice.
%% @item @code{sympref ipc native}: use the @code{py} interface to
%% interact directly with an embedded Python interpreter, e.g.,
%% provided by the Octave Pythonic package.
%% @item @code{sympref ipc system}: construct a long string of
%% the command and pass it directly to the python interpreter with
%% the @code{system()} command.  This typically assembles a multiline
%% string for the commands, except on Windows where a long one-line
%% string is used.
%% @item @code{sympref ipc systmpfile}: output the python commands
%% to a file in @code{tempdir} and then call that.
%% For debugging, will not be supported long-term.
%% @item @code{sympref ipc sysoneline}: put the python commands all
%% on one line and pass to @code{python -c} using a call to @code{system()}.
%% For debugging, will not be supported long-term.
%% @end itemize
%%
%% Except for @code{native}, all of these communication interfaces
%% depend on the current @strong{Python executable}, which can be queried:
%% @example
%% @comment doctest: +SKIP
%% sympref python
%%   @result{} ans = python
%% @end example
%% Changing this might help if you've installed
%% a local Python interpreter somewhere else on your system.
%% The value can be changed by setting the environment variable
%% @code{PYTHON}, which can be configured in the OS, or it can be
%% set within Octave using:
%% @example
%% @comment doctest: +SKIP
%% setenv PYTHON python3
%% setenv PYTHON $@{HOME@}/.local/bin/python
%% setenv PYTHON C:\Python\python.exe
%% sympref reset
%% @end example
%% If the environment variable is empty or not set, the package
%% uses a default setting (often @code{python}).
%%
%%
%% @strong{Reset}: reset the SymPy communication mechanism.  This can be
%% useful after an error occurs and the connection with Python
%% becomes confused.
%%
%% @example
%% @group
%% sympref reset                              % doctest: +SKIP
%% @end group
%% @end example
%%
%%
%% @strong{Default precision}: control the number of digits used by
%% variable-precision arithmetic (see also the @ref{digits} command).
%%
%% @example
%% @group
%% sympref digits          % get
%%   @result{} ans = 32
%% sympref digits 64       % set
%% sympref digits default
%% @end group
%% @end example
%%
%% Be @strong{quiet} by minimizing startup and diagnostics messages:
%% @example
%% @group
%% sympref quiet
%%   @result{} ans = 0
%% sympref quiet on
%% sympref quiet default
%% @end group
%% @end example
%%
%% Report the @strong{version} number:
%%
%% @example
%% @group
%% sympref version
%%   @result{} 3.1.1+
%% @end group
%% @end example
%%
%% @seealso{sym, syms}
%% @end deftypefn

function varargout = sympref(cmd, arg)

  persistent settings

  if (isempty(settings))
    settings = 42;
    sympref('defaults')
  end

  if (nargin == 0)
    varargout{1} = settings;
    return
  end

  if (isstruct (cmd))
    assert (isequal (sort (fieldnames (cmd)), ...
      sort ({'ipc'; 'display'; 'digits'; 'quiet'})), ...
      'sympref: structure has incorrect field names')
    settings = [];
    sympref ('quiet', cmd.quiet)
    sympref ('display', cmd.display)
    sympref ('digits', cmd.digits)
    sympref ('ipc', cmd.ipc)
    return
  end

  switch lower(cmd)
    case 'defaults'
      settings = [];
      settings.ipc = 'default';
      sympref ('display', 'default')
      sympref ('digits', 'default')
      sympref ('quiet', 'default')

    case 'version'
      assert (nargin == 1)
      varargout{1} = '3.1.1+';

    case 'display'
      if (nargin == 1)
        varargout{1} = settings.display;
      else
        arg = lower (arg);
        if (strcmp (arg, 'default'))
          arg = 'unicode';
          if (ispc () && (~isunix ()))
            % Unicode not working on Windows, Issue #83.
            arg = 'ascii';
          end
        end
        assert(strcmp(arg, 'flat') || strcmp(arg, 'ascii') || ...
               strcmp(arg, 'unicode'))
        settings.display = arg;
      end

    case 'digits'
      if (nargin == 1)
        varargout{1} = settings.digits;
      else
        if (ischar(arg))
          if (strcmpi(arg, 'default'))
            arg = 32;
          else
            arg = str2double(arg);
          end
        end
        arg = int32(arg);
        assert(arg > 0, 'precision must be positive')
        settings.digits = arg;
      end

    case 'snippet'
      warning ('OctSymPy:deprecated', ...
               'Debugging mode "snippet" has been removed');

    case 'quiet'
      if (nargin == 1)
        varargout{1} = settings.quiet;
      else
        if (strcmpi(arg, 'default'))
          settings.quiet = false;
        else
          settings.quiet = tf_from_input(arg);
        end
      end

    case 'python'
      if (nargin ~= 1)
        error('old syntax ''sympref python'' removed; use ''setenv PYTHON'' instead')
      end
      pyexec = getenv('PYTHON');
      if (isempty(pyexec))
        pyexec = defaultpython ();
      end
      varargout{1} = pyexec;

    case 'ipc'
      if (nargin == 1)
        varargout{1} = settings.ipc;
      else
        verbose = ~sympref('quiet');
        sympref('reset')
        settings.ipc = arg;
        switch arg
          case 'default'
            msg = 'Choosing the default [autodetect] communication mechanism';
          case 'native'
            msg = 'Forcing the native Python/C API communication mechanism';
          case 'system'
            msg = 'Forcing the system() communication mechanism';
          case 'popen2'
            msg = 'Forcing the popen2() communication mechanism';
          case 'systmpfile'
            msg = 'Forcing systmpfile ipc: warning: this is for debugging';
          case 'sysoneline'
            msg = 'Forcing sysoneline ipc: warning: this is for debugging';
          otherwise
            msg = '';
            if (~ ischar (arg))
              arg = num2str (arg);
            end
            warning('OctSymPy:sympref:invalidarg', ...
                    'Unsupported IPC mechanism ''%s'': hope you know what you''re doing', ...
                    arg)
        end
        if (verbose)
          disp(msg)
        end
      end

    case 'reset'
      verbose = ~sympref('quiet');
      python_ipc_driver('reset', []);

    %case 'path'
      %pkg_path = fileparts (mfilename ('fullpath'));
      % or
      %pkg_l = pkg ('list');
      %idx = strcmp ('octsympy', cellfun (@(x) x.name, pkg_l, 'UniformOutput', false));
      %if (~ any (idx))
      %  error ('the package %s is not installed', your_pkg);
      %end
      %pkg_path = pkg_l{idx}.dir

    case 'diagnose'
      if (strcmp (lower (sympref ('ipc')), 'default') && ...
          exist ('pyversion') && exist ('pyexec') && exist ('pyeval'))
        % TODO: see note in ipc_native
        assert_pythonic_and_sympy (true)
      else
        assert_have_python_and_sympy (sympref ('python'), true)
      end

    otherwise
      error ('sympref: invalid preference or command ''%s''', lower (cmd));
  end
end

function r = tf_from_input(s)

  if (~ischar(s))
    r = logical(s);
  elseif (strcmpi(s, 'true'))
    r = true;
  elseif (strcmpi(s, 'false'))
    r = false;
  elseif (strcmpi(s, 'on'))
    r = true;
  elseif (strcmpi(s, 'off'))
    r = false;
  elseif (strcmpi(s, '[]'))
    r = false;
  else
    r = str2double(s);
    assert(~isnan(r), 'invalid expression to convert to bool')
    r = logical(r);
  end
end

%!shared sympref_orig
%! sympref_orig = sympref ();

%!test
%! % test quiet, side effect of making following tests a bit less noisy!
%! sympref quiet on
%! a = sympref('quiet');
%! assert(a == 1)

%% Test for correct line numbers in Python exceptions.  These first tests
%% happen with the default ipc---usually popen2.  We don't explicitly test
%% popen2 because test suite should be portable.
%!error <line 1> pycall_sympy__ ('raise ValueError');
%!error <line 1> pycall_sympy__ ('raise ValueError', sym('x'));
%!error <line 1> pycall_sympy__ ('raise ValueError', sym([1 2 3; 4 5 6]));
%!error <line 1> pycall_sympy__ ('raise ValueError', {1; 1; 1});
%!error <line 1> pycall_sympy__ ('raise ValueError', struct('a', 1, 'b', 'word'));
%!error <line 2> pycall_sympy__ ( {'x = 1' 'raise ValueError'} );
%!error <line 3> pycall_sympy__ ( {'x = 1' 'pass' '1/0'} );
%!error <line 3> pycall_sympy__ ( {'a=1' 'b=1' 'raise ValueError' 'c=1' 'd=1'} );

%% Test for correct line error in Python exceptions.
%!error <raise ValueError> pycall_sympy__ ('raise ValueError');
%!error <raise ValueError> pycall_sympy__ ('raise ValueError', sym('x'));
%!error <raise ValueError> pycall_sympy__ ('raise ValueError', sym([1 2 3; 4 5 6]));
%!error <raise ValueError> pycall_sympy__ ('raise ValueError', {1; 1; 1});
%!error <raise ValueError> pycall_sympy__ ('raise ValueError', struct('a', 1, 'b', 'word'));
%!error <raise ValueError> pycall_sympy__ ( {'x = 1' 'raise ValueError'} );
%!error <1/0> pycall_sympy__ ( {'x = 1' 'pass' '1/0'} );
%!error <raise ValueError> pycall_sympy__ ( {'a=1' 'b=1' 'raise ValueError' 'c=1' 'd=1'} );

%!test
%! % system should work on all system, but just runs sysoneline on windows
%! sympref('ipc', 'system');
%! syms x

%!error <line 1> pycall_sympy__ ('raise ValueError')
%!error <line 1> pycall_sympy__ ('raise ValueError', sym('x'))
%!error <line 1> pycall_sympy__ ('raise ValueError', struct('a', 1, 'b', 'word'))

%!error <c=1; raise ValueError> pycall_sympy__ ({'a=1' 'b=1' 'c=1; raise ValueError' 'd=1'});

%!test
%! % sysoneline should work on all systems
%! sympref('ipc', 'sysoneline');
%! syms x

%!error <line 1> pycall_sympy__ ('raise ValueError')
%!error <line 1> pycall_sympy__ ('raise ValueError', struct('a', 1, 'b', 'word'))

%!test
%! sympref('ipc', 'systmpfile');
%! syms x

%!error <line 1> pycall_sympy__ ('raise ValueError')
%!error <line 1> pycall_sympy__ ('raise ValueError', struct('a', 1, 'b', 'word'))

%!test
%! s = warning ('off', 'OctSymPy:sympref:invalidarg');
%! sympref ('ipc', 'bogus');
%! assert (strcmp (sympref ('ipc'), 'bogus'))
%! warning (s)

%!error <invalid ipc mechanism> syms ('x')

%!test
%! sympref('defaults')
%! assert(strcmp(sympref('ipc'), 'default'))
%! sympref('quiet', 'on')

%!test
%! % restore sympref from structure
%! old = sympref ();
%! sympref ('display', 'ascii');
%! sympref ('digits', 64);
%! old = orderfields (old);  % re-ordering the fields should be ok
%! sympref (old);
%! new = sympref ();
%! assert (isequal (old, new))

%!error <incorrect field names>
%! s.a = 'hello';
%! s.b = 'world';
%! sympref (s)

%!test
%! syms x
%! sympref('reset');
%! % restore original sympref settings
%! sympref ('ipc',   sympref_orig.ipc);
%! syms x
%! sympref ('quiet', sympref_orig.quiet);

%!error <invalid preference or command> sympref ('nosuchsetting')
%!error <invalid preference or command> sympref ('nosuchsetting', true)

repo: octsympy


text:
function [soln,eqn,info] = Poisson(node,elem,bdFlag,pde,option)
%% POISSON Poisson equation: P1 linear element.
%
%   u = Poisson(node,elem,bdFlag,pde) produces the linear finite element
%   approximation of the Poisson equation
%
%       -div(d*grad(u))=f  in \Omega, with
%       Dirichlet boundary condition u=g_D on \Gamma_D,
%       Neumann boundary condition   d*grad(u)*n=g_N on \Gamma_N,
%       Robin boundary condition     g_R*u + d*grad(u)*n=g_N on \Gamma _R
%
%   The mesh is given by node and elem and the boundary edge is by
%   bdFlag. See meshdoc, bddoc for details. The data is given by the
%   structure pde which contains function handles f, g_D, g_N, g_R, or d.
%   For general elliptic equations with convection and reaction
%   coefficients, see ellipticpde.
%
%   soln = Poisson(node,elem,bdFlag,pde,option) specifies the options.
%
%   In the output, soln structure contains
%     - soln.u: solution u
%     - soln.Du: gradient of u
%
%   In the input, option structures contains
%    - option.dquadorder: quadrature order for diffusion coefficients
%    - option.fquadorder: quadrature order for computing right hand side f
%    - option.solver
%      'direct': the built in direct solver \ (mldivide)
%      'mg':     multigrid-type solvers mg is used.
%      'amg':    algebraic multigrid method is used.
%      'none':   only assemble the matrix equation but not solve.
%   The default setting is to use the direct solver for small size problems
%   and multigrid solvers for large size problems. For more options on the
%   multigrid solver mg, type help mg.
%
%   The function Poisson assembes the matrix equation AD*u = b and solves
%   it by the direct solver (small size <= 2e3) or the multigrid solver
%   (large size > 2e3). The Dirichlet boundary condition is built into the
%   matrix AD and the Neumann boundary condition is build into b.
%
%   The diffusion coefficient d is a scalar function or a column array with
%   the same length as the elem array.
%
%   When only one type of boundary condition is imposed, the input argument
%   bdFlag can be skipped. The boundary condition is implicitly given in
%   the pde structure by specifying g_D or g_N only. See examples below.
%
%   [soln,eqn] = Poisson(node,elem,bdFlag,pde) returns also the equation
%   structure eqn, which includes:
%     - eqn.AD:  modified stiffness matrix AD;
%     - eqn.b:   the right hand side.
%     - eqn.Lap: non-modified stiffness matrix
%
%   The solution u = AD\b. The matrix eqn.Lap can be used to evulate the
%   bilinear form a(u,v) = u'*eqn.Lap*v, especially the enery norm of a finite
%   element function u is given by by sqrt(u'*eqn.Lap*u).
%
%   [soln,eqn,info] = Poisson(node,elem,bdFlag,pde) returns also the
%   information on the assembeling and solver, which includes:
%     - info.assembleTime: time to assemble the matrix equation
%     - info.solverTime:   time to solve the matrix equation
%     - info.itStep:       number of iteration steps for the mg solver
%     - info.error:        l2 norm of the residual b - A*u
%     - info.flag:         flag for the mg solver.
%       flag = 0: converge within max iteration
%       flag = 1: iterated maxIt times but did not converge
%       flag = 2: direct solver
%       flag = 3: no solve
%
%   Example
%     squarePoisson; Poissonfemrate;
%
%   Example
%     clear all
%     node = [0,0; 1,0; 1,1; 0,1];
%     elem = [2,3,1; 4,1,3];
%     for k = 1:4
%       [node,elem] = uniformrefine(node,elem);
%     end
%     % Homogenous Dirichlet boundary condition
%     pde.f = inline('ones(size(p,1),1)','p');
%     pde.g_D = 0;
%     u = Poisson(node,elem,[],pde);
%     figure(1);
%     showresult(node,elem,u);
%
%   See also Poisson3, femPoisson, squarePoisson, Poissonfemrate
%
%   Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Preprocess
if ~exist('bdFlag','var'), bdFlag = []; end
if ~exist('option','var'), option = []; end
% important constants
N = size(node,1);
NT = size(elem,1);
Ndof = N;

%% Diffusion coefficient
time = cputime;  % record assembling time
if ~isfield(pde,'d'), pde.d = []; end
if ~isfield(option,'dquadorder'), option.dquadorder = 1; end
if ~isempty(pde.d) && isnumeric(pde.d)
   K = pde.d;                                 % d is an array
end
if ~isempty(pde.d) && ~isnumeric(pde.d)       % d is a function
    [lambda,weight] = quadpts(option.dquadorder);
    nQuad = size(lambda,1);
    K = zeros(NT,1);
    for p = 1:nQuad
		pxy = lambda(p,1)*node(elem(:,1),:) ...
			+ lambda(p,2)*node(elem(:,2),:) ...
			+ lambda(p,3)*node(elem(:,3),:);
        K = K + weight(p)*pde.d(pxy);
   end
end

%% Compute geometric quantities and gradient of local basis
[Dphi,area] = gradbasis(node,elem);

%% Assemble stiffness matrix
A = sparse(Ndof,Ndof);
for i = 1:3
    for j = i:3
        % $A_{ij}|_{\tau} = \int_{\tau}K\nabla \phi_i\cdot \nabla \phi_j dxdy$
        Aij = (Dphi(:,1,i).*Dphi(:,1,j) + Dphi(:,2,i).*Dphi(:,2,j)).*area;
        if ~isempty(pde.d)
            Aij = K.*Aij;
        end
        if (j==i)
            A = A + sparse(elem(:,i),elem(:,j),Aij,Ndof,Ndof);
        else
            A = A + sparse([elem(:,i);elem(:,j)],[elem(:,j);elem(:,i)],...
                           [Aij; Aij],Ndof,Ndof);
        end
    end
end
clear K Aij

%% Assemble the right hand side
b = zeros(Ndof,1);
if ~isfield(option,'fquadorder')
    option.fquadorder = 3;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if isreal(pde.f) % f is a real number or vector and not a function
   switch length(pde.f)
       case NT  % f is piecewise constant
         bt = pde.f.*area/3;
         b = accumarray(elem(:),[bt; bt; bt],[Ndof 1]);
       case N   % f is piecewise linear
         bt = zeros(NT,3);
         bt(:,1) = area.*(2*pde.f(elem(:,1)) + pde.f(elem(:,2)) + pde.f(elem(:,3)))/12;
         bt(:,2) = area.*(2*pde.f(elem(:,2)) + pde.f(elem(:,3)) + pde.f(elem(:,1)))/12;
         bt(:,3) = area.*(2*pde.f(elem(:,3)) + pde.f(elem(:,1)) + pde.f(elem(:,2)))/12;
         b = accumarray(elem(:),bt(:),[Ndof 1]);
       case 1   % f is a scalar e.g. f = 1
         bt = pde.f*area/3;
         b = accumarray(elem(:),[bt; bt; bt],[Ndof 1]);
   end
end
if ~isempty(pde.f) && ~isreal(pde.f)  % f is a function
    [lambda,weight] = quadpts(option.fquadorder);
    phi = lambda;                 % linear bases
	nQuad = size(lambda,1);
    bt = zeros(NT,3);
    for p = 1:nQuad
		% quadrature points in the x-y coordinate
		pxy = lambda(p,1)*node(elem(:,1),:) ...
			+ lambda(p,2)*node(elem(:,2),:) ...
			+ lambda(p,3)*node(elem(:,3),:);
		fp = pde.f(pxy);
        for i = 1:3
            bt(:,i) = bt(:,i) + weight(p)*phi(p,i)*fp;
        end
    end
    bt = bt.*repmat(area,1,3);
    b = accumarray(elem(:),bt(:),[Ndof 1]);
end
clear pxy bt

%% Set up boundary conditions
[AD,b,u,freeNode,isPureNeumann] = getbd(b);

%% Record assembling time
assembleTime = cputime - time;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
if isempty(freeNode), return; end
% Set up solver type
if isempty(option) || ~isfield(option,'solver')  || isfield(option,'mgoption')   % no option.solver
    if Ndof <= 2e3  % Direct solver for small size systems
        option.solver = 'direct';
    else            % MGCG  solver for large size systems
        option.solver = 'mg';
    end
end
if isPureNeumann
    option.solver = 'mg';
end
solver = option.solver;
% solve
switch solver
    case 'direct'
        t = cputime;
        u(freeNode) = AD(freeNode,freeNode)\b(freeNode);
        residual = norm(b - AD*u);
        info = struct('solverTime',cputime - t,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'none'
        info = struct('solverTime',[],'itStep',0,'err',[],'flag',3,'stopErr',[]);
    case 'mg'
        if ~isfield(option,'mgoption')   % no option.mgoption
            option.mgoption.x0 = u;
            option.mgoption.solver = 'CG';
        end
        [u,info] = mg(AD,b,elem,option.mgoption);
    case 'amg'
        if ~isfield(option,'amgoption')  % no option.amgoption
            option.amgoption.x0 = u;
            option.amgoption.solver = 'CG';
        end
        [u(freeNode),info] = amg(AD(freeNode,freeNode),b(freeNode),option.amgoption);
end
% post-process for pure Neumann problem
if isPureNeumann
    patchArea = accumarray(elem(:),[area;area;area]/3, [N 1]);
    uc = sum(u.*patchArea)/sum(area);
    u = u - uc;   % int u = 0
end

%% Compute Du
dudx =  u(elem(:,1)).*Dphi(:,1,1) + u(elem(:,2)).*Dphi(:,1,2) ...
      + u(elem(:,3)).*Dphi(:,1,3);
dudy =  u(elem(:,1)).*Dphi(:,2,1) + u(elem(:,2)).*Dphi(:,2,2) ...
      + u(elem(:,3)).*Dphi(:,2,3);
Du = [dudx, dudy];

%% Output
if nargout == 1
    soln = u;
else
    soln = struct('u',u,'Du',Du);
    eqn = struct('A',AD,'b',b,'freeNode',freeNode,'Lap',A);
    info.assembleTime = assembleTime;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [AD,b,u,freeNode,isPureNeumann] = getbd(b)
    %% Set up of boundary conditions.
    %
    % 1) Modify the matrix for Dirichlet boundary nodes, which are not degree
    % of freedom. Values at these nodes are evaluatation of pde.g_D. The
    % original stiffness matrix A is turn into the matrix AD by enforcing
    % AD(fixedNode,fixedNode)=I, AD(fixedNode,freeNode)=0, AD(freeNode,fixedNode)=0.
    %
    % 2) Modify the right hand side b. The Neumann boundary integral is added
    % to b. For Dirichlet boundary ndoes, b(fixedNode) is the evaluation of
    % pde.g_D.
    %
    % Special attentation should be given for the pure Neumann boundary
    % condition. To enforce the compatible condition, the vector b should have
    % mean value zero. To avoid a singular matrix, the 1st node is chosen as
    % fixedNode.
    %
    % The order of assigning Neumann and Dirichlet boundary condition is
    % important to get the right setting at the intersection nodes of Dirichlet
    % and Neumann boundary edges.
    %
    % Reference: Long Chen. Finite Element Methods and its Programming. Lecture
    % Notes.

    u = zeros(Ndof,1);
    %% Initial check
    if ~isfield(pde,'g_D'), pde.g_D = []; end
    if ~isfield(pde,'g_N'), pde.g_N = []; end
    if ~isfield(pde,'g_R'), pde.g_R = []; end

    %% Part 1: Modify the matrix for Dirichlet and Robin condition
    % Robin boundary condition
    Robin = [];
    isRobin = (bdFlag(:) == 3);
    if any(isRobin)
        allEdge = [elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])];
        Robin = allEdge(isRobin,:);
    end
    if ~isempty(Robin) && ~isempty(pde.g_R) && ~(isnumeric(pde.g_R) && (pde.g_R == 0))
        ve = node(Robin(:,1),:) - node(Robin(:,2),:);
        edgeLength = sqrt(sum(ve.^2,2));
        mid = (node(Robin(:,1),:) + node(Robin(:,2),:))/2;
        % use Simplson rule to compute int g_R phi_iphi_j ds
        ii = [Robin(:,1),Robin(:,1),Robin(:,2),Robin(:,2)];
        jj = [Robin(:,1),Robin(:,2),Robin(:,1),Robin(:,2)];
        temp = pde.g_R(mid).*edgeLength;
        ss = [1/3*temp, 1/6*temp, 1/6*temp, 1/3*temp];
        A = A + sparse(ii,jj,ss,Ndof,Ndof);
    end

    % Find Dirichlet boundary nodes: fixedNode
    fixedNode = []; freeNode = [];
    if ~isempty(bdFlag) % find boundary edges and boundary nodes
        [fixedNode,bdEdge,isBdNode] = findboundary(elem,bdFlag);
        freeNode = ~isBdNode;
    end
    if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
        % no bdFlag, only pde.g_D is given
        [fixedNode,bdEdge,isBdNode] = findboundary(elem);
        freeNode = ~isBdNode;
    end

    % Modify the matrix for different boundary conditions
    % Dirichlet boundary condition
    % Build Dirichlet boundary condition into the matrix AD by enforcing
    % AD(fixedNode,fixedNode)=I, AD(fixedNode,freeNode)=0, AD(freeNode,fixedNode)=0.
    if ~isempty(fixedNode)
        bdidx = zeros(Ndof,1);
        bdidx(fixedNode) = 1;
        Tbd = spdiags(bdidx,0,Ndof,Ndof);
        T = spdiags(1-bdidx,0,Ndof,Ndof);
        AD = T*A*T + Tbd;
    end
    % Neumann boundary condition
    isPureNeumann = false;
    if isempty(fixedNode) && isempty(Robin) % pure Neumann boundary condition
        isPureNeumann = true;
        AD = A;
        AD(1,1) = AD(1,1) + 1e-6;
%         fixedNode = 1;
%         freeNode = 2:Ndof;    % eliminate the kernel by enforcing u(1) = 0;
    end
    % Robin boundary condition
    if isempty(fixedNode) && ~isempty(Robin)
        AD = A;
    end

    %% Part 2: Find boundary edges and modify the right hand side b
    % Find boundary edges: Neumann
    Neumann = [];
    if ~isempty(bdFlag)  % bdFlag specifies different bd conditions
        Neumann = bdEdge;
    end
    if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
        % no bdFlag, only pde.g_N or pde.g_R is given in the input
        [tempvar,Neumann] = findboundary(elem); %#ok<ASGLU>
    end

    % Neumann boundary condition
    if  isnumeric(pde.g_N) && all(pde.g_N == 0)
        pde.g_N = [];
    end
    if ~isempty(Neumann) && ~isempty(pde.g_N)
        el = sqrt(sum((node(Neumann(:,1),:) - node(Neumann(:,2),:)).^2,2));
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 2;   % default order exact for linear gN
        end
        [lambdagN,weightgN] = quadpts1(option.gNquadorder);
        phigN = lambdagN;                 % linear bases
        nQuadgN = size(lambdagN,1);
        ge = zeros(size(Neumann,1),2);
        for pp = 1:nQuadgN
            % quadrature points in the x-y coordinate
            ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
                 + lambdagN(pp,2)*node(Neumann(:,2),:);
            gNp = pde.g_N(ppxy);
            for igN = 1:2
                ge(:,igN) = ge(:,igN) + weightgN(pp)*phigN(pp,igN)*gNp;
            end
        end
        ge = ge.*repmat(el,1,2);
        b = b + accumarray(Neumann(:), ge(:),[Ndof,1]);
    end
    % The case with non-empty Neumann edges but g_N=0 or g_N=[] corresponds to
    % the zero flux boundary condition on Neumann edges and no modification of
    % A,u,b is needed.

    % Dirichlet boundary condition
    if isnumeric(pde.g_D) && all(pde.g_D == 0)   % zero g_D
        pde.g_D = [];
    end
    if ~isPureNeumann && ~isempty(fixedNode) && ~isempty(pde.g_D)
        if isnumeric(pde.g_D)  % pde.g_D could be a numerical array
            u(fixedNode) = pde.g_D(fixedNode);
        else % pde.g_D is a function handle
            u(fixedNode) = pde.g_D(node(fixedNode,:));
        end
        b = b - A*u;
    end
    if ~isempty(fixedNode) % non-empty Dirichlet boundary condition
        b(fixedNode) = u(fixedNode);
    end
    % The case with non-empty Dirichlet nodes but g_D=0 or g_D=[] corresponds
    % to the zero Dirichlet boundary condition and no modification of u,b is
    % needed.

    % Pure Neumann boundary condition
    if isPureNeumann
        b = b - mean(b); % compatiable condition (f,1) + <gN,1> = 0
%         b(1) = 0;
    end
    end % end of getbd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end % end of Poisson

repo: ifem


text:
function L = evaluate_log_posterior(this, uv)
%EVALUATE_LOG_POSTERIOR computes the log-posterior (negative energy) of the
%   flow fields UV
%   Actually only proportional to the log posterior since the variance of neither the
%   spatial nor the data terms is considered
%
% This is a member function of the class 'classic_nl_optical_flow'.
%
% Authors: Deqing Sun, Department of Computer Science, Brown University
% Contact: dqsun@cs.brown.edu
% $Date: $
% $Revision: $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

% Spatial term
S = this.spatial_filters;
p = 0;

for i = 1:length(S)

    u_ = conv2(uv(:,:,1), S{i}, 'valid');
    v_ = conv2(uv(:,:,2), S{i}, 'valid');

    if isa(this.rho_spatial_u{i}, 'robust_function')

        p = p - sum(evaluate(this.rho_spatial_u{i}, u_(:)))...
            - sum(evaluate(this.rho_spatial_v{i}, v_(:)));

    elseif isa(this.rho_spatial_u{i}, 'gsm_density')

        p   = p + sum(evaluate_log(this.rho_spatial_u{i}, u_(:)'))...
                    + sum(evaluate_log(this.rho_spatial_v{i}, v_(:)'));

    else
        error('evaluate_log_posterior: unknown rho function!');
    end;
end;

% likelihood
It  = partial_deriv(this.images, uv, this.interpolation_method, this.deriv_filter);

if isa(this.rho_data, 'robust_function')

    l   = -sum(evaluate(this.rho_data, It(:)));

elseif isa(this.rho_data, 'gsm_density')

    l   = sum(evaluate_log(this.rho_data, It(:)'));

else
    error('evaluate_log_posterior: unknown rho function!');
end;

L = this.lambda*p + l;

if this.display
    fprintf('spatial\t%3.2e\tdata\t%3.2e\n', this.lambda*p, l);
end;

author: kristinbranson

repo: JAABA


text:
function [Gs, op, nodes] = mk_nbrs_of_dag(G0)
% MK_NBRS_OF_DAG Make all DAGs that differ from G0 by a single edge deletion, addition or reversal
% [Gs, op, nodes] = mk_nbrs_of_dag(G0)
%
% Gs{i} is the i'th neighbor.
% op{i} = 'add', 'del', or 'rev' is the operation used to create the i'th neighbor.
% nodes(i,1:2) are the head and tail of the operated-on arc.

Gs = {};
op = {};
nodes = [];

[I,J] = find(G0);
nnbrs = 1;
% all single edge deletions
for e=1:length(I)
  i = I(e); j = J(e);
  G = G0;
  G(i,j) = 0;
  Gs{nnbrs} = G;
  op{nnbrs} = 'del';
  nodes(nnbrs, :) = [i j];
  nnbrs = nnbrs + 1;
end

% all single edge reversals
for e=1:length(I)
  i = I(e); j = J(e);
  G = G0;
  G(i,j) = 0;
  G(j,i) = 1;
  if acyclic(G)
    Gs{nnbrs} = G;
    op{nnbrs} = 'rev';
    nodes(nnbrs, :) = [i j];
    nnbrs = nnbrs + 1;
  end
end

[I,J] = find(~G0);
% all single edge additions
for e=1:length(I)
  i = I(e); j = J(e);
  if i ~= j % don't add self arcs
    G = G0;
    G(i,j) = 1;
    if G(j,i)==0 % don't add i->j if j->i exists already
      if acyclic(G)
	Gs{nnbrs} = G;
	op{nnbrs} = 'add';
	nodes(nnbrs, :) = [i j];
	nnbrs = nnbrs + 1;
      end
    end
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [N,binEdges] = BF_SimpleBinner(xData,numBins)
% BF_SimpleBinner Generate a histogram from equally spaced bins.
% An alternative to the more opaque method of Matlab's internal histcounts function.
%
%---INPUTS:
% xData, a data vector.
% numBins, the number of bins.
%---OUTPUTS:
% N, the counts
% binEdges, the extremeties of the bins.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

minX = min(xData);
maxX = max(xData);

% Linearly spaced bins:
binEdges = linspace(minX,maxX,numBins+1);

N = zeros(numBins,1);
for i = 1:numBins
    if i < numBins
        N(i) = sum(xData >= binEdges(i) & xData < binEdges(i+1));
    else
        % the final bin
        N(i) = sum(xData >= binEdges(i) & xData <= binEdges(i+1));
    end
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function lambda = symm_random_eigenvalues ( n, d, key )

%*****************************************************************************80
%
%% SYMM_RANDOM_EIGENVALUES returns the eigenvalues of the SYMM_RANDOM matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real D(N), the desired eigenvalues for the matrix.
%
%    Input, integer KEY, a positive integer that selects the data.
%
%    Output, real LAMBDA(N), the eigenvalues.
%
  d = d(:);

  lambda = d;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% hohmann.m          July 9, 2013

% Hohmann two impulse orbit transfer between
% planar and non-coplanar circular orbits

% includes three-dimensional orbit graphics
% and graphical primer vector analysis

% Orbital Mechanics with MATLAB

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

global rtd dtr pvi pvdi

global mu req hn1 hn2 hn3 dinc

% astrodynamic and utility constants

om_constants;

% Brent root-finding tolerance

rtol = 1.0e-8;

% request inputs

clc; home;

fprintf('\nHohmann Orbit Transfer Analysis\n');

while (1)

    fprintf('\n\nplease input the initial altitude (kilometers)\n');

    alt1 = input('? ');

    if (alt1 > 0.0)
        break;
    end

end

while (1)

    fprintf('\n\nplease input the final altitude (kilometers)\n');

    alt2 = input('? ');

    if (alt2 > 0.0)
        break;
    end

end

while (1)

    fprintf('\n\nplease input the initial orbital inclination (degrees)');
    fprintf('\n(0 <= inclination <= 180)\n');

    inc1 = input('? ');

    if (inc1 >= 0.0 && inc1 <= 180.0)
        break;
    end

end

while (1)

    fprintf('\n\nplease input the final orbital inclination (degrees)');
    fprintf('\n(0 <= inclination <= 180)\n');

    inc2 = input('? ');

    if (inc2 >= 0.0 && inc2 <= 180.0)
        break;
    end

end

% convert orbit inclinations to radians

inc1 = inc1 * dtr;

inc2 = inc2 * dtr;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve the orbit transfer problem %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% calculate total inclination change (radians)

dinc = abs(inc2 - inc1);

% compute geocentric radii of initial and final orbits (kilometers)

r1 = req + alt1;

r2 = req + alt2;

% compute "normalized" radii

hn1 = sqrt(2.0 * r2 / (r2 + r1));

hn2 = sqrt(r1 / r2);

hn3 = sqrt(2.0 * r1 / (r2 + r1));

% compute "local circular velocity" of initial and final orbits (km/sec)

v1 = sqrt(mu / r1);

v2 = sqrt(mu / r2);

% compute transfer orbit semimajor axis (kilometers)

smat = 0.5 * (r1 + r2);

% compute transfer orbit eccentricity (non-dimensional)

ecct = (max(r1, r2) - min(r1, r2)) / (r1 + r2);

% compute transfer orbit perigee and apogee radii and velocities

rp = smat * (1.0 - ecct);

ra = smat * (1.0 + ecct);

vt1 = sqrt(2.0 * mu * ra / (rp * (rp + ra)));

vt2 = sqrt(2.0 * mu * rp / (ra * (rp + ra)));

% compute transfer orbit period (seconds)

taut = 2.0 * pi * sqrt(smat^3 / mu);

tof = 0.5 * taut;

if (abs(dinc) == 0)

    % coplanar orbit transfer

    if (r2 > r1)

        % higher-to-lower transfer

        dv1 = vt1 - v1;

        dv2 = v2 - vt2;

    else

        % lower-to-higher transfer

        dv1 = v1 - vt2;

        dv2 = vt1 - v2;

    end

    dinc1 = 0;

    dinc2 = 0;

    inct = inc1;

else

    % non-coplanar orbit transfer

    [xroot, froot] = brent('hohmfunc', 0, dinc, rtol);

    % calculate delta-v's

    dinc1 = xroot;

    dinc2 = dinc - dinc1;

    dv1 = v1 * sqrt(1.0 + hn1 * hn1 - 2.0 * hn1 * cos(dinc1));

    dv2 = v1 * sqrt(hn2 * hn2 * hn3 * hn3 + hn2 * hn2 ...
        - 2.0 * hn2 * hn2 * hn3 * cos(dinc2));

    if (inc2 > inc1)

        inct = inc1 + dinc1;

    else

        inct = inc1 - dinc1;

    end

end

% print results

clc; home;

fprintf('\nHohmann Orbit Transfer Analysis');
fprintf('\n-------------------------------\n\n');

fprintf('initial orbit altitude            %10.4f kilometers \n\n', alt1);

fprintf('initial orbit radius              %10.4f kilometers \n\n', alt1 + req);

fprintf('initial orbit inclination         %10.4f degrees \n\n', inc1 * rtd);

fprintf('initial orbit velocity            %10.4f meters/second \n\n\n', 1000.0 * v1);

fprintf('final orbit altitude              %10.4f kilometers \n\n', alt2);

fprintf('final orbit radius                %10.4f kilometers \n\n', alt2 + req);

fprintf('final orbit inclination           %10.4f degrees \n\n', inc2 * rtd);

fprintf('final orbit velocity              %10.4f meters/second \n', 1000.0 * v2);

fprintf('\n\nfirst inclination change          %10.4f degrees\n\n', dinc1 * rtd);

fprintf('second inclination change         %10.4f degrees\n\n', dinc2 * rtd);

fprintf('total inclination change          %10.4f degrees\n\n\n', rtd * (dinc1 + dinc2));

fprintf('first delta-v                     %10.4f meters/second \n\n', 1000.0 * dv1);

fprintf('second delta-v                    %10.4f meters/second \n\n', 1000.0 * dv2);

fprintf('total delta-v                     %10.4f meters/second \n\n\n', 1000.0 * (dv1 + dv2));

fprintf('transfer orbit semimajor axis     %10.4f kilometers \n\n', smat);

fprintf('transfer orbit eccentricity       %10.8f \n\n', ecct);

fprintf('transfer orbit inclination        %10.4f degrees \n\n', rtd * inct);

fprintf('transfer orbit perigee velocity   %10.4f meters/second \n\n', 1000.0 * vt1);

fprintf('transfer orbit apogee velocity    %10.4f meters/second \n\n', 1000.0 * vt2);

fprintf('transfer orbit coast time         %10.4f seconds \n', tof);

fprintf('                                  %10.4f minutes \n', tof / 60.0);

fprintf('                                  %10.4f hours \n\n', tof / 3600.0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create trajectory graphics %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% load orbital elements arrays, create state vectors and plot orbits

oevi(1) = r1;
oevi(2) = 0.0;
oevi(3) = inc1;
oevi(4) = 0.0;
oevi(5) = 0.0;

% determine correct true anomaly (radians)

if (alt2 > alt1)

    oevi(6) = 0.0;

else

    oevi(6) = 180.0 * dtr;

end

[ri, vi] = orb2eci(mu, oevi);

oevti(1) = smat;
oevti(2) = ecct;
oevti(3) = inct;
oevti(4) = 0.0;
oevti(5) = 0.0;

% determine correct true anomaly (radians)

if (alt2 > alt1)

    oevti(6) = 0.0;

else

    oevti(6) = 180.0 * dtr;

end

[rti, vti] = orb2eci(mu, oevti);

oevtf(1) = smat;
oevtf(2) = ecct;
oevtf(3) = inct;
oevtf(4) = 0.0;
oevtf(5) = 0.0;

% determine correct true anomaly (radians)

if (alt2 > alt1)

    oevtf(6) = 180.0 * dtr;

else

    oevtf(6) = 0.0;

end

[rtf, vtf] = orb2eci(mu, oevtf);

oevf(1) = r2;
oevf(2) = 0.0;
oevf(3) = inc2;
oevf(4) = 0.0;
oevf(5) = 0.0;

% determine correct true anomaly (radians)

if (alt2 > alt1)

    oevf(6) = 180.0 * dtr;

else

    oevf(6) = 0.0;

end

[rf, vf] = orb2eci(mu, oevf);

% compute orbital periods

period1 = 2.0 * pi * oevi(1) * sqrt(oevi(1) / mu);

period2 = 2.0 * pi * oevti(1) * sqrt(oevti(1) / mu);

period3 = 2.0 * pi * oevf(1) * sqrt(oevf(1) / mu);

deltat1 = period1 / 300;

simtime1 = -deltat1;

deltat2 = 0.5 * period2 / 300;

simtime2 = -deltat2;

deltat3 = period3 / 300;

simtime3 = -deltat3;

for i = 1:1:301

    simtime1 = simtime1 + deltat1;

    simtime2 = simtime2 + deltat2;

    simtime3 = simtime3 + deltat3;

    % compute initial orbit "normalized" position vector

    [rwrk, vwrk] = twobody2 (mu, simtime1, ri, vi);

    rp1_x(i) = rwrk(1) / req;

    rp1_y(i) = rwrk(2) / req;

    rp1_z(i) = rwrk(3) / req;

    % compute transfer orbit position vector

    [rwrk, vwrk] = twobody2 (mu, simtime2, rti, vti);

    rp2_x(i) = rwrk(1) / req;

    rp2_y(i) = rwrk(2) / req;

    rp2_z(i) = rwrk(3) / req;

    % compute final orbit position vector

    [rwrk, vwrk] = twobody2 (mu, simtime3, rf, vf);

    rp3_x(i) = rwrk(1) / req;

    rp3_y(i) = rwrk(2) / req;

    rp3_z(i) = rwrk(3) / req;

end

figure(1);

% create axes vectors

xaxisx = [1 1.5];
xaxisy = [0 0];
xaxisz = [0 0];

yaxisx = [0 0];
yaxisy = [1 1.5];
yaxisz = [0 0];

zaxisx = [0 0];
zaxisy = [0 0];
zaxisz = [1 1.5];

figure (1);

hold on;

grid on;

% plot earth

[x y z] = sphere(24);

h = surf(x, y, z);

colormap([127/255 1 222/255]);

set (h, 'edgecolor', [1 1 1]);

% plot coordinate system axes

plot3(xaxisx, xaxisy, xaxisz, '-g', 'LineWidth', 1);

plot3(yaxisx, yaxisy, yaxisz, '-r', 'LineWidth', 1);

plot3(zaxisx, zaxisy, zaxisz, '-b', 'LineWidth', 1);

% plot initial orbit

plot3(rp1_x, rp1_y, rp1_z, '-r', 'LineWidth', 1.5);

plot3(rp1_x(1), rp1_y(1), rp1_z(1), 'ob');

% plot transfer orbit

plot3(rp2_x, rp2_y, rp2_z, '-b', 'LineWidth', 1.5);

plot3(rp2_x(end), rp2_y(end), rp2_z(end), 'ob');

% plot final orbit

plot3(rp3_x, rp3_y, rp3_z, '-g', 'LineWidth', 1.5);

xlabel('X coordinate (ER)', 'FontSize', 12);

ylabel('Y coordinate (ER)', 'FontSize', 12);

zlabel('Z coordinate (ER)', 'FontSize', 12);

title('Hohmann Transfer: Initial, Transfer and Final Orbits', 'FontSize', 16);

axis equal;

view(50, 20);

rotate3d on;

print -depsc -tiff -r300 hohmann1.eps

%%%%%%%%%%%%%%%%%%%%%%%%
% create primer graphics
%%%%%%%%%%%%%%%%%%%%%%%%

dvi = (vti - vi)';

dvf = (vf - vtf)';

% perform primer vector initialization

pviniz(tof, rti, vti, dvi, dvf);

% number of graphic data points

npts = 300;

% plot behavior of primer vector magnitude

dt = tof / npts;

for i = 1:1:npts + 1

    t = (i - 1) * dt;

    if (t == 0)

        % initial value of primer magnitude and derivative

        pvm = norm(pvi);

        pvdm = dot(pvi, pvdi) / pvm;

    else

        % primer vector and derivative magnitudes at time t

        [pvm, pvdm] = pvector(rti, vti, t);

    end

    % load data array

    x1(i) = t / 60.0;

    y1(i) = pvm;

    y2(i) = pvdm;

end

figure(2);

hold on;

plot(x1, y1, '-r', 'LineWidth', 1.5);

plot(x1(1), y1(1), 'or');

plot(x1(end), y1(end), 'or');

title('Primer Vector Analysis of the Hohmann Transfer', 'FontSize', 16);

xlabel('simulation time (minutes)', 'FontSize', 12);

ylabel('primer vector magnitude', 'FontSize', 12);

grid;

% create eps graphics file with tiff preview

print -depsc -tiff -r300 primer.eps;

% plot behavior of magnitude of primer derivative

figure(3);

hold on;

plot(x1, y2, '-r', 'LineWidth', 1.5);

plot(x1(1), y2(1), 'or');

plot(x1(end), y2(end), 'or');

title('Primer Vector Analysis of the Hohmann Transfer', 'FontSize', 16);

xlabel('simulation time (minutes)', 'FontSize', 12);

ylabel('primer derivative magnitude', 'FontSize', 12);

grid;

% create eps graphics file with tiff preview

print -depsc -tiff -r300 primer_der.eps;

author: Sable

repo: mcbench-benchmarks


text:
function dcenter = ComputeDistFlyCenter(obj,n,fly1)

nflies = obj.nfliespermovie(n);
flyidx1 = obj.getFlyIdx(n,fly1);

% initialize
x_mm1 = obj.GetPerFrameData('x_mm',n,fly1);
y_mm1 = obj.GetPerFrameData('y_mm',n,fly1);
nframes1 = obj.nframes(flyidx1);
firstframe1 = obj.firstframes(flyidx1);
endframe1 = obj.endframes(flyidx1);
dcenter = nan(nflies,nframes1);

for fly2 = 1:nflies,
  if fly2 == fly1, continue; end

  flyidx2 = obj.getFlyIdx(n,fly2);
  firstframe2 = obj.firstframes(flyidx2);
  endframe2 = obj.endframes(flyidx2);

  % get start and end frames of overlap
  t0 = max(firstframe1,firstframe2);
  t1 = min(endframe1,endframe2);

  % no overlap
  if t1 < t0, continue; end

  % indices for these frames
  offi = firstframe1-1;
  offj = firstframe2-1;
  i0 = t0 - offi;
  i1 = t1 - offi;
  j0 = t0 - offj;
  j1 = t1 - offj;

  x_mm2 = obj.GetPerFrameData('x_mm',n,fly2);
  y_mm2 = obj.GetPerFrameData('y_mm',n,fly2);

  % centroid distance
  dx = x_mm2(j0:j1)-x_mm1(i0:i1);
  dy = y_mm2(j0:j1)-y_mm1(i0:i1);
  z = sqrt(dx.^2 + dy.^2);
  dcenter(fly2,i0:i1) = z;

end