text:
function [ci,bootstat]  = stat_bootci(nboot,bootfun,varargin)
%BOOTCI Bootstrap Confidence Interval
%   CI = BOOTCI(NBOOT,BOOTFUN,...) computes the 95 percent BCa bootstrap
%   confidence interval of the statistic defined by the function BOOTFUN.
%   NBOOT is a positive integer indicating the number of bootstrap data
%   samples used in the computation. BOOTFUN is a function handle specified
%   with @. The third and later input arguments to BOOTCI are data
%   (scalars, column vectors, or matrices) that are used to create inputs
%   to BOOTFUN. BOOTCI creates each bootstrap sample by sampling with
%   replacement from the rows of the non-scalar data arguments (these must
%   have the same number of rows). Scalar data are passed to BOOTFUN
%   unchanged.
%
%   If BOOTFUN returns a scalar, CI is a vector containing the lower and
%   upper bounds of the confidence interval. If BOOTFUN returns a vector of
%   length M, CI is an array of size 2-by-M, where CI(1,:) are lower bounds
%   and CI(2,:) are upper bounds. If BOOTFUN returns an array of size
%   M-by-N-by-P-by-..., CI is an array of size 2-by-M-by-N-by-P-by-...,
%   where CI(1,:,:,:,...) is an array of lower bounds and CI(2,:,:,:,...)
%   is an array of upper bounds.
%
%   CI = BOOTCI(NBOOT,{BOOTFUN,...},'alpha',ALPHA) computes the 100*(1-ALPHA)
%   percent BCa bootstrap confidence interval of the statistic defined by the
%   function BOOTFUN. ALPHA is a scalar between 0 and 1. The default value of
%   ALPHA is 0.05. A cell array groups BOOTFUN and the arguments used to create
%   inputs to it. ALPHA and any other arguments to BOOTCI appear outside the
%   cell array.
%
%   CI = BOOTCI(NBOOT,{BOOTFUN,...},...,'type',TYPE) computes the bootstrap
%   confidence interval of the statistic defined by the function BOOTFUN.
%   TYPE is the confidence interval type, specifying different methods of
%   computing the confidence interval. TYPE is a string chosen from
%       'norm' or 'normal':               normal approximated interval with
%                                         bootstrapped bias and standard
%                                         error;
%       'per' or 'percentile':            basic percentile method;
%       'cper' or 'corrected percentile': bias corrected percentile method;
%       'bca' :                           bias corrected and accelerated
%                                         percentile method;
%       'stud' or 'student':              studentized confidence interval.
%   The default value of TYPE is 'bca'.
%
%   CI = BOOTCI(NBOOT,{BOOTFUN,...},...,'type','stud','nbootstd',NBOOTSTD)
%   computes the studentized bootstrap confidence interval of the statistic
%   defined by the function BOOTFUN. The standard error of the bootstrap
%   statistics is estimated using bootstrap with NBOOTSTD bootstrap data
%   samples. NBOOTSTD is a positive integer value. The default value of
%   NBOOTSTD is 100.
%
%   CI = BOOTCI(NBOOT,{BOOTFUN,...},...,'type','stud','stderr',STDERR)
%   computes the studentized bootstrap confidence interval of statistics
%   defined by the function BOOTFUN. The standard error of the bootstrap
%   statistics is evaluated by the function STDERR. STDERR is a function
%   handle created using @. STDERR should take the same arguments as
%   BOOTFUN and return the standard error of the statistic computed by
%   BOOTFUN.
%
%   CI = BOOTCI(NBOOT,{BOOTFUN,...},...,'Weights',WEIGHTS) specifies
%   observation weights. WEIGHTS must be a vector of non-negative numbers
%   with at least one positive element. The number of elements in WEIGHTS
%   must be equal to the number of rows in non-scalar input arguments to
%   BOOTFUN. To obtain one bootstrap replicate, BOOTSTRP samples N out of N
%   with replacement using these weights as multinomial sampling
%   probabilities.
%
%   CI = BOOTCI(NBOOT,{BOOTFUN,...},...,'Options',OPTIONS) contains options
%   that specify whether to compute bootstrap iterations in parallel, and how
%   to use random numbers during the bootstrap sampling. OPTIONS is a struct
%   that can be created by a call to STATSET. BOOTCI uses the following fields:
%       'UseParallel'
%       'UseSubstreams'
%       'Streams'
%   For information on these fields see PARALLELSTATS.
%   NOTE: if 'UseParallel' is 'always' and 'UseSubstreams'
%   is 'never', then the length of Streams must equal the number
%   of processors used by BOOTCI. There are two possibilities.
%   If a MATLAB pool is open, then Streams is the same length as
%   the size of the MATLAB pool. If a MATLAB pool is not open,
%   then Streams must supply a single random number stream.
%
%   [CI,BOOTSTAT] = BOOTCI(...) also returns the bootstrapped statistic
%   computed for each of the NBOOT bootstrap replicate samples.  Each row of
%   BOOTSTAT contains the results of applying BOOTFUN to one bootstrap sample.
%   If BOOTFUN returns a matrix or array, then this output is converted to a
%   row vector for storage in BOOTSTAT.
%
%   Example:
%     Compute the confidence interval for the capability index in
%     statistical process control:
%          y = normrnd(1,1,30,1);                  % simulated process data
%          LSL = -3;  USL = 3;                     % process specifications
%          capable = @(x) (USL-LSL)./(6* std(x));  % process capability
%          bootci(2000,capable, y)                 % Bca confidence interval
%          bootci(2000,{capable, y},'type','per')  % basic percentile method
%
%   See also: BOOTSTRP, JACKKNIFE, STATSET, STATGET,
%   RANDSAMPLE, PARFOR, PARALLELSTATS.

% The BCa method is described in the following references:
%
%     T.J. DiCicio and B. Efron (1996), "Bootstrap confidence intervals,"
%     Statistical Science, v. 11, n. 3, pp. 189-228.
%
%     B. Efron and R.J. Tibshirani (1993), An Introduction to the
%     Bootstrap, Chapman & Hall, New York.
%
% Their formula involves a Z0 factor that is computed using the proportion
% of bootstrap values less than the original sample value.  In order to get
% reasonable results when the sample is lumpy, we include half of the
% bootstrap values that are tied with the original sample value when we
% compute Z0.

% Copyright 2005-2010 The MathWorks, Inc.
%   $Revision: 1.1.8.8 $  $Date: 2011/07/20 00:08:12 $

if nargin<2
    error(message('stats:bootci:TooFewInputs'));
end;
if nboot<=0 || nboot~=round(nboot)
    error(message('stats:bootci:BadNboot'))
end;

if ~iscell(bootfun) % default syntax
    type = 'bca';
    alpha = .05;
    fun = bootfun;
    data = varargin;
    weights = [];
    bootstrpOptions = statset('bootstrp');
else % syntax with optional type, alpha, nbootstd, and stderrfun name/value pairs
    fun = bootfun{1};
    data = bootfun(2:end);
    pnames = {'type', 'alpha', 'stderr', 'nbootstd', 'weights', 'options'};
    dflts =  {'bca', .05, [], 100, [], statset('bootstrp')};
    [type,alpha,stderrfun,nbootstd,weights,bootstrpOptions] = ...
                           internal.stats.parseArgs(pnames, dflts, varargin{:});
end

% error check for the bootfun
try
    obsstat = fun(data{:});
catch ME
    m = message('stats:bootci:BadBootfun',func2str(fun));
    throw(addCause(MException(m.Identifier,'%s',getString(m)),ME));
end
if any(~isfinite(obsstat))
        error(message('stats:bootci:NonfiniteBootfun'));
end;

% Save the original size of stat and reshape it
sz = [];
if ~isvector(obsstat)
    sz = size(obsstat);
end
obsstat = obsstat(:)'; % turn into a row-vector

% call subfunctions to compute the intervals
switch (lower(type))
    case {'norm','normal'}
        [ci,bootstat] = bootnorm(obsstat,nboot,fun,alpha,weights,bootstrpOptions,data{:});
    case {'per','percentile'}
        [ci,bootstat] = bootper(nboot,fun,alpha,weights,bootstrpOptions,data{:});
    case {'cper', 'corrected percentile'}
        [ci,bootstat] = bootcper(obsstat,nboot,fun,alpha,weights,bootstrpOptions,data{:});
    case 'bca'
        [ci,bootstat] = bootbca(obsstat,nboot,fun,alpha,weights,bootstrpOptions,data{:});
    case {'stud','student'}
        [ci,bootstat] = bootstud(obsstat,nboot,fun,alpha,nbootstd,stderrfun,...
            weights,bootstrpOptions,data{:});
    otherwise
        error(message('stats:bootci:BadType'))
end;

% Reshape
if ~isempty(sz)
    ci = reshape(ci,[2 sz]);
end

end   % bootci()

%-------------------------------------------------------------------------
function [ci,bstat] = bootnorm(stat,nboot,bootfun,alpha,weights,bootstrpOptions,varargin)
% normal approximation interval
% A.C. Davison and D.V. Hinkley (1996), p198-200

bstat = bootstrp(nboot,bootfun,varargin{:},'weights',weights,'Options',bootstrpOptions);

se = std(bstat,0,1);   % standard deviation estimate
bias = mean(bsxfun(@minus,bstat,stat),1);
za = norminv(alpha/2);   % normal confidence point
lower = stat - bias + se*za; % lower bound
upper = stat - bias - se*za;  % upper bound

% return
ci = [lower;upper];
end   % bootnorm()

%-------------------------------------------------------------------------
function [ci,bstat] = bootper(nboot,bootfun,alpha,weights,bootstrpOptions,varargin)
% percentile bootstrap CI

bstat = bootstrp(nboot,bootfun,varargin{:},'weights',weights,'Options',bootstrpOptions);

pct1 = 100*alpha/2;
pct2 = 100-pct1;
lower = prctile(bstat,pct1,1);
upper = prctile(bstat,pct2,1);

% return
ci =[lower;upper];
end % bootper()

%-------------------------------------------------------------------------
function [ci,bstat] = bootcper(stat,nboot,bootfun,alpha,weights,bootstrpOptions,varargin)
% corrected percentile bootstrap CI
% B. Efron (1982), "The jackknife, the bootstrap and other resampling
% plans", SIAM.

bstat = bootstrp(nboot,bootfun,varargin{:},'weights',weights,'Options',bootstrpOptions);

% stat is transformed to a normal random variable z0.
% z0 = invnormCDF[ECDF(stat)]
z_0 = fz0(bstat,stat);
z_alpha = norminv(alpha/2); % normal confidence point

% transform z0 back using the invECDF[normCDF(2z0-za)] and
% invECDF[normCDF(2z0+za)]
pct1 = 100*normcdf(2*z_0-z_alpha);
pct2 = 100*normcdf(2*z_0+z_alpha);

% inverse ECDF
m = numel(stat);
lower = zeros(1,m);
upper = zeros(1,m);
for i=1:m
    lower(i) = prctile(bstat(:,i),pct2(i),1);
    upper(i) = prctile(bstat(:,i),pct1(i),1);
end

% return
ci = [lower;upper];
end % bootcper()

%-------------------------------------------------------------------------
function [ci,bstat] = bootbca(stat,nboot,bootfun,alpha,weights,bootstrpOptions,varargin)
% corrected and accelerated percentile bootstrap CI
% T.J. DiCiccio and B. Efron (1996), "Bootstrap Confidence Intervals",
% statistical science, 11(3)

bstat = bootstrp(nboot,bootfun,varargin{:},'weights',weights,'Options',bootstrpOptions);

% same as bootcper, this is the bias correction
z_0 = fz0(bstat,stat);

% apply jackknife
try
    jstat = jackknife(bootfun,varargin{:},'Options',bootstrpOptions);
catch ME
    m = message('stats:bootci:JackknifeFailed',func2str(bootfun));
    throw(addCause(MException(m.Identifier,'%s',getString(m)),ME));
end
N = size(jstat,1);
if isempty(weights)
    weights = repmat(1/N,N,1);
else
    weights = weights(:);
    weights = weights/sum(weights);
end

% acceleration finding, see DiCiccio and Efron (1996)
mjstat = sum(bsxfun(@times,jstat,weights),1); % mean along 1st dim.
score = bsxfun(@minus,mjstat,jstat); % score function at stat; ignore (N-1) factor because it cancels out in the skew
iszer = all(score==0,1);
skew = sum(bsxfun(@times,score.^3,weights),1) ./ ...
    (sum(bsxfun(@times,score.^2,weights),1).^1.5) /sqrt(N); % skewness of the score function
skew(iszer) = 0;
acc = skew/6;  % acceleration

% transform back with bias corrected and acceleration
z_alpha1 = norminv(alpha/2);
z_alpha2 = -z_alpha1;
pct1 = 100*normcdf(z_0 +(z_0+z_alpha1)./(1-acc.*(z_0+z_alpha1)));
pct1(z_0==Inf) = 100;
pct1(z_0==-Inf) = 0;
pct2 = 100*normcdf(z_0 +(z_0+z_alpha2)./(1-acc.*(z_0+z_alpha2)));
pct2(z_0==Inf) = 100;
pct2(z_0==-Inf) = 0;

% inverse of ECDF
m = numel(stat);
lower = zeros(1,m);
upper = zeros(1,m);
for i=1:m
    lower(i) = prctile(bstat(:,i),pct2(i),1);
    upper(i) = prctile(bstat(:,i),pct1(i),1);
end

% return
ci = sort([lower;upper],1);
end % bootbca()

%-------------------------------------------------------------------------
function [ci,bstat] = bootstud(stat,nboot,bootfun,alpha,nbootstd,stderrfun,weights,bootstrpOptions,varargin)
% studentized bootstrap CI with bootstrp to estimate the se
% T.J. DiCiccio and B. Efron (1996), "Bootstrap Confidence Intervals",
% statistical science, 11(3)

% Should we use stderrfun to compute st. dev. or an inner bootstrap loop?
if isempty(stderrfun) % studentized with bootstrap error
    if nbootstd<=0 || nbootstd~=round(nbootstd)
        error(message('stats:bootci:BadNbootstd'))
    end
else  % studentized with stderrfun fun
    % error check for stderrfun
    try
        out=stderrfun(varargin{:});
    catch ME
        m = message('stats:bootci:BadStderr',func2str(stderrfun));
        throw(addCause(MException(m.Identifier,'%s',getString(m)),ME));
    end
    if any(~isfinite(out))
        error(message('stats:bootci:NonfiniteStderr'));
    end
    if isvector(out)
        out = out(:)';
    end
    if size(stat)~=size(out)
        error(message('stats:bootci:BadStderrSize'));
    end
end;

% bootstrap
[bstat,bootsam] = bootstrp(nboot,bootfun,varargin{:}, ...
    'weights',weights,'Options',bootstrpOptions); % bootstrap statistics

% find non-scalar data in varargin
la = length(varargin);
scalard = zeros(la,1);
for k = 1:la
   [row,col] = size(varargin{k});
   if max(row,col) == 1
      scalard(k) = 1;
   end
   if row == 1 && col ~= 1
      varargin{k} = varargin{k}(:);
   end
end

% bootstrap generated bootstrap replica to get student errors
N = size(bootsam,1);
sd_t = zeros(nboot,numel(stat));
for b=1:nboot
    db = cell(la,1);
    for k = 1:la
        % store the bootstrap data samples in a cell array
        if scalard(k) == 0
            db{k} = varargin{k}(bootsam(:,b),:);
        else
            db{k} = varargin{k};
        end
    end
    if isempty(weights)
        w = [];
    else
        w = weights(bootsam(:,b));
    end
    if isempty(stderrfun)
        bstatstd = bootstrp(nbootstd,bootfun,db{:},'weights',w,...
            'Options',bootstrpOptions);
        sd_t(b,:) = std(bstatstd,0,1);
    else
        bstatstd = stderrfun(db{:});
        sd_t(b,:) = reshape(bstatstd,1,numel(stat));
    end
end

% studentized statistics
above0 = sd_t>0;
mbstat = bsxfun(@minus,bstat,stat);
tstat = zeros(size(mbstat));
tstat(above0) = mbstat(above0)./sd_t(above0);

% percentiles for the studentized stats are computed.
lower = prctile(tstat,100*alpha/2,1);
upper = prctile(tstat,100*(1-alpha/2),1);

% back to the original stats from the studentized stats
lower = lower.*std(bstat,0,1) + stat;
upper = upper.*std(bstat,0,1) + stat;

% return
ci = [lower; upper];
end

% -------------------------
function z0=fz0(bstat,stat)
% Compute bias-correction constant z0
z0 = norminv(mean(bsxfun(@lt,bstat,stat),1) + mean(bsxfun(@eq,bstat,stat),1)/2);
end   % fz0()

repo: REST


text:
function surfacedata = quarticssurface

surfacedata = struct('phi',@phi,'gradient', @gradient, 'unitoutnormal',...
    @unitoutnormal, 'project', @project, 'initmesh',@initmesh,...
    'meancurvature', @meancurvature, 'u', @u, 'f', @f, 'gra_u',@gra_u);

function [node, elem] = initmesh
M = load('meshdata/quartics4528.mat');
node = M.node;
elem = M.elem;
end

function z = phi(p)
% level set function
x = p(:,1); y = p(:,2); z = p(:,3);
z = (x.^2 - 1).^2 + (y.^2 - 1).^2 + (z.^2 - 1).^2 - 1.05;
end

function n = unitoutnormal(p)
%
n = gradient(p);
l = sqrt(sum(n.^2,2));
n = n./repmat(l,1,3);
end

function n = gradient(p)
% gradient of phi
x = p(:,1); y = p(:,2); z = p(:,3);
n = [4*(x.^2 - 1).*x,4*(y.^2 - 1).*y, 4*(z.^2 - 1).*z];
end

function [node,dist] = project(p)
% projection function

s = sign(phi(p));

node = p;

normalAtNode = gradient(node);
valueAtNode = phi(node);
node = node - valueAtNode*ones(1,3).*normalAtNode./(dot(normalAtNode,normalAtNode,2)*ones(1,3));

vector = (-s)*ones(1,3).*(node - p);

d = s.*sqrt(dot(vector,vector,2));

normalAtNode = gradient(node);

node = p - d*ones(1,3).*normalAtNode./(sqrt(dot(normalAtNode,normalAtNode,2))*ones(1,3));

valueAtNode = phi(node);

normalAtNode = gradient(node);

vector = (-s)*ones(1,3).*(node - p);

d = s.*sqrt(dot(vector,vector,2));

e1 = normalAtNode./(sqrt(dot(normalAtNode, normalAtNode,2))*ones(1,3))-vector./(sqrt(dot(vector,vector,2))*ones(1,3));
error=sqrt(valueAtNode.^2./(dot(normalAtNode, normalAtNode,2))+dot(e1,e1,2));

k=1;
while max(abs(error)) > 1e-6 && k<200

    k=k+1;

    node = node - valueAtNode*ones(1,3).*normalAtNode./(dot(normalAtNode,normalAtNode,2)*ones(1,3));

    vector = -s*ones(1,3).*(node - p);
    d = s.*sqrt(dot(vector,vector,2));
    normalAtNode = gradient(node);

    node = p - d*ones(1,3).*normalAtNode./(sqrt(dot(normalAtNode,normalAtNode,2))*ones(1,3));

    valueAtNode = phi(node);

    normalAtNode = gradient(node);

    vector = (-s)*ones(1,3).*(node - p);

    d = s.*sqrt(dot(vector,vector,2));
    e1 = normalAtNode./(sqrt(dot(normalAtNode, normalAtNode,2))*ones(1,3))-vector./(sqrt(dot(vector,vector,2))*ones(1,3));
    error=sqrt(valueAtNode.^2./(dot(normalAtNode, normalAtNode,2))+dot(e1,e1,2));
end

if nargout == 2
    dist = d;
end
end

function mc = meancurvature(p)
x = p(:,1); y = p(:,2); z = p(:,3);
s = (x.^2-1).^2.*x.^2 + (y.^2-1).^2.*y.^2 + (z.^2-1).^2.*z.^2;
s1 =(3 .* z .^ 2 - 2 + 3 .* y .^ 2) .* x .^ 6 + (-6 .* y .^ 2 - 6 .* z .^ 2 + 4) .* x .^ 4 ...
    + (-2 + 6 .* y .^ 2 + 3 .* z .^ 6 - 6 .* z .^ 4 - 6 .* y .^ 4 + 6 .* z .^ 2 + 3 .* y .^ 6) .* x .^ 2 ...
+ (3 .* z .^ 2 - 2) .* y .^ 6 + (4 - 6 .* z .^ 2) .* y .^ 4 + ...
    (-2 + 6 .* z .^ 2 + 3 .* z .^ 6 - 6 .* z .^ 4) .* y .^ 2 - 2 .* z .^ 6 - 2 .* z .^ 2 + 4 .* z .^ 4;

mc = 0.5 * s1 ./ (s.^(3/2));
end

function z = f(p)
p = project(p);
x = p(:,1); y = p(:,2); z = p(:,3);
x2 = x.^2; y2 = y.^2; z2 = z.^2;
x4 = x2.^2; y4 = y2.^2; z4 = z2.^2;

x31 = x2.*x - x;
y31 = y2.*y - y;
z31 = z2.*z - z;

d1 = x31.^2 + y31.^2 + z31.^2;
d2 = x31 + y31 + z31;
t = exp(x+y+z);

d3 = (z2 + y2 - 2/3).*x4.*(x2 - 2);
d4 = (z4.*(z2-2) + 2*z2 + y4.*(y2 -2) + 2*y2 -2/3).*x2;
d5 = (z2 - 2/3).*y4.*(y2 - 2);
d6 = (z4.*(z2 -2) + 2*z2 - 2/3).*y2;

z = -2*t.*((-(y31+z31).*x31 - z31.*y31 + d1).*d1 -1.5*(d3 + d4 + d5 + d6 - 2/3*z31.^2).*d2)./d1.^2;
end

function z = u(p)
% exact solution
p = project(p);
z = exp(sum(p,2));
end

function z = gra_u(p)

p = project(p);
x = p(:,1); y = p(:,2); z = p(:,3);
x2 = x.^2; y2 = y.^2; z2 = z.^2;

x31 = x2.*x - x;
y31 = y2.*y - y;
z31 = z2.*z - z;

d1 = x31.^2 + y31.^2 + z31.^2;
t = exp(x+y+z);
z1 = (y31+z31).*x31 - y31.^2 - z31.^2;
z2 = (x31+z31).*y31 - x31.^2 - z31.^2;
z3 = (x31+y31).*z31 - x31.^2 - y31.^2;

z = -(t./d1)*ones(1,3).*[z1,z2,z3];
end

end

repo: ifem


text:
function disprog(i,N,steps);
%DISPROG Display progression of a loop.
%	DISPROG(i,N,steps) displays the progression of a loop.
%
%	I     : loop variable
%	N     : final value of i
%	STEPS : number of displayed steps.
%
%	Example:
%	 N=16; for i=1:N, disprog(i,N,5); end;

%	F. Auger, August, December 1995.
%       from an idea of R. Settineri.
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

global begin_time_disprog ;

if (i==1),
 begin_time_disprog=cputime;
end;

if (i==N),
 fprintf('100 %% complete in %g seconds.\n', cputime-begin_time_disprog);
 clear begin_time_disprog;
elseif (floor(i*steps/N)~=floor((i-1)*steps/N)),
 fprintf('%g ', floor(i*steps/N)*ceil(100/steps));
end;

author: HeLiangHIT


text:
classdef lds < dml.method
%LDS linear dynamical system.
%
%   DESCRIPTION
%   data X is represented as trials x features x timepoints
%
%   State can be partially observed/unobserved during training.
%   partial observability of multiple observations is also supported
%   observations are normally distributed conditional on the state.
%   multiple observation sequences are supported
%
%   bias term is *NOT* automatically added to the model
%
%   K = number of states
%   M = number of observations
%   T = number of timesteps
%
%   X and Y are swapped wrt the Kalman filter conventions
%
%   EXAMPLE
%   rand('seed',2); randn('seed',2);
%
%   nsamples = 1000; ncov = 2; ncycles = 10; ntrials = 2;
%   Y = sin(ncycles * 2 * pi * (1:nsamples) ./ nsamples);
%   Y = repmat(reshape(Y,[1 1 numel(Y)]),[ntrials 1 1]);
%   X = repmat(Y,[1 ncov 1]) + 0.5*randn(ntrials,ncov,nsamples);
%
%   k = dml.lds('inference','smooth','verbose',true,'indims',[ncov nsamples]);
%   k = k.train(X,Y);
%   Z = k.test(X);
%   figure
%   plot(squeeze(Y(1,:,:))','k');
%   hold on;
%   plot(squeeze(Z(1,:,:))','r');
%   legend('real','predicted');
%
%   k = dml.lds('verbose',true);
%   k = k.train(X,nan(size(Y))); % hidden state estimation
%   U = repmat(Y,[1 ncov 1]) + 0.1*randn(ntrials,ncov,nsamples);
%   Z = k.test(U);
%   figure
%   plot(squeeze(Y(1,:,:))','k--','LineWidth',2);
%   hold on;
%   plot(zscore(squeeze(Z(1,:,:))'),'r','LineWidth',2);
%   plot(squeeze(X(1,:,:))','bo');
%   legend('real','predicted','observed');
%
%   k = dml.lds('inference','smooth','verbose',true);
%   k = k.train(X,[Y nan(size(Y))]); % mixture of hidden + observed states
%   U = repmat(Y,[1 ncov 1]) + 0.1*randn(ntrials,ncov,nsamples);
%   Z = k.test(X);
%   figure
%   plot(zscore(squeeze(Y(1,:,:))'),'k');
%   hold on;
%   plot(zscore(squeeze(Z(1,:,:))'),'r');
%   plot(zscore(squeeze(X(1,:,:))'),'bo');
%   legend('real','predicted','observed');
%
%   REFERENCES
%   Pattern Recognition and Machine Learning, Bishop
%   A unifying review of linear dynamical systems, Gharamani
%
%   DEVELOPER
%   Marcel van Gerven (m.vangerven@donders.ru.nl)
%   Ali Bahramisharif (ali@cs.ru.nl)

  properties

    maxiter = 100; % maximum number of EM iterations
    thresh = 1e-4; % EM convergence threshold
    diagQ = 0; % regularize state noise to diagonal (0 <= diagQ <=1)
    diagR = 0; % regularize measurement noise to diagonal (0 <= diagR <=1)

    inference = 'smooth'; % filter / smooth

    epsilon = 0 % added to the diagonal of the covariance matrices for numerical stability (e.g. 1e-7);

    A % K x K transition matrix for the unobservable state
    C % M x K emission matrix for the observations
    mu0 % K x 1 the initial mean of the hidden state
    V0 % K x K the initial hidden noise covariance
    R % M x M measurement noise covariance
    Q % K x K state noise covariance

    loglik % log likelihood

    nhidden = 1 % number of hidden states; automatically determined if Y is given

  end

  methods

    function obj = lds(varargin)

      obj = obj@dml.method(varargin{:});
    end

    function obj = train(obj,X,Y)

      % resize data if indims is specified
      if ~isempty(obj.indims)
        X = reshape(X,[size(X,1) obj.indims]);
      end

      % cast to cell array (multiple observation sequences)
      % representation as variables * timepoints
      if ~iscell(X)
        if ndims(X) ~= 3, error('input should be of size trials x features x time'); end
        Xt = cell(1,size(X,1));
        for c=1:length(Xt)
          Xt{c} = reshape(X(c,:,:),[size(X,2) size(X,3)]);
        end
        X = Xt; clear Xt;
      else
        for c=1:length(X)
          X{c} = X{c};
        end
      end

      if nargin < 3

        % assume nhidden signals underlying the observations
        Y = cell(1,length(X));
        for c=1:length(X)
          Y{c} = nan(obj.nhidden,size(X{c},2));
        end

      else

        if ~iscell(Y)
          if ndims(Y) ~= 3, error('input should be of size trials x features x time'); end
          Yt = cell(1,size(Y,1));
          for c=1:length(Yt)
            Yt{c} = reshape(Y(c,:,:),[size(Y,2) size(Y,3)]);
          end
          Y = Yt; clear Yt;
        else
          for c=1:length(Y)
            Y{c} = Y{c}';
          end
        end

      end

      % remove empty elements
      eidx = cellfun(@(x)(isempty(x)),X);
      eidx = eidx & cellfun(@(x)(isempty(x)),Y);
      X = X(~eidx);
      Y = Y(~eidx);

      K = size(Y{1},1);
      M = size(X{1},1);
      N = length(X);

      if isempty(obj.A), obj.A = 1e-3*randn(K,K); end
      if isempty(obj.C), obj.C = 1e-3*randn(M,K); end
      if isempty(obj.mu0), obj.mu0 = 1e-3*randn(K,1); end
      if isempty(obj.V0), obj.V0 = 1e-3*eye(K); end
      if isempty(obj.R), obj.R = 1e-3*eye(M); end
      if isempty(obj.Q), obj.Q = 1e-3*eye(K); end

%       % DEBUG
%       [obj.A, obj.C, obj.Q, obj.R, obj.mu0, obj.V0, obj.loglik] = learn_kalman(X,obj.A,obj.C,obj.Q,obj.R,obj.mu0,obj.V0);
%       return;

      oldLL = inf;
      LL = 0;
      loglik = [];
      iter = 0;
      while abs(LL - oldLL) > obj.thresh && iter < obj.maxiter

        oldLL = LL;

        % E step
        LL=0;
        for c=1:N % iterate over sequences

          [G1t,G2t,G3t,G4t,G5t,G6t,mu0t,V0t,LLt] = obj.Estep(X{c},Y{c});

          if c==1

            G1=G1t;
            G2=G2t;
            G3=G3t;
            G4=G4t;
            G5=G5t;
            G6=G6t;
            mu0=mu0t;
            V0=V0t + mu0t*mu0t';
            LL=LLt;

          else

            G1 = G1 + G1t;
            G2 = G2 + G2t;
            G3 = G3 + G3t;
            G4 = G4 + G4t;
            G5 = G5 + G5t;
            G6 = G6 + G6t;

            mu0 = cat(2,mu0,mu0t);
            V0 = V0 + V0t + mu0t*mu0t';

            LL = LL + LLt;

          end

        end

        loglik = [loglik LL];

        if obj.verbose && ~isnan(LL)
          fprintf('EM step: %d; log likelihood: %g\n',iter,LL);
        end

        if iter && (LL < oldLL), fprintf('non-decreasing log likelihood!\n'); end

        % M step
        T = sum(cellfun(@(x)(size(x,2)),X));

        obj.C = G6 / G1;

        obj.R = (G5 - obj.C * G6') ./ T;
        obj.R = (obj.R + obj.R') ./ 2;

        if obj.diagR
          obj.R = (1-obj.diagR) * obj.R + obj.diagR * diag(diag(obj.R));
        end

        obj.A = G4 / G3;

        obj.Q = (G2 - obj.A * G4') ./ (T-N);
        obj.Q = (obj.Q + obj.Q') ./ 2;

        if obj.diagQ
          obj.Q = (1-obj.diagQ) * obj.Q + obj.diagQ * diag(diag(obj.Q));
        end
        obj.mu0 = mean(mu0,2);
        obj.V0 = V0/N - obj.mu0*obj.mu0'+(mu0-repmat(obj.mu0,[1 N]))*(mu0-repmat(obj.mu0,[1 N]))'/N;

        % symmetricize and make positive semidefinite
        obj.V0 = (obj.V0 + obj.V0') ./ 2;
        obj.V0 = obj.V0 + obj.epsilon*eye(size(obj.V0));

        iter = iter + 1;

        if ~any(isnan(Y{1}(:)))
          break;
        end

      end

      if obj.verbose
        fprintf('EM step: %d; log likelihood: %g\n',iter,LL);
      end

      obj.loglik = loglik(2:end);

    end

    function [mu,V,loglik] = test(obj,X)
     % LDS inference

     if ~iscell(X)
       if ndims(X) ~= 3, error('input should be of size trials x features x time'); end
       Xt = cell(1,size(X,1));
       for c=1:length(Xt)
         Xt{c} = squeeze(X(c,:,:));
       end
       X = Xt; clear Xt;
     end

     nsets = length(X);

     mu = cell(1,nsets);
     for c=1:nsets

       [m,V,loglik] = obj.filter(X{c});

       if strcmp(obj.inference,'smooth')
         [m,V] = obj.smooth(m,V);
       end

       mu{c} = m';

     end

     % reshape if all sequences are of the same size
     if numel(unique(cellfun(@(x)(size(x,1)),mu)))==1
      Z = zeros(length(mu),size(mu{1},2),size(mu{1},1));
      for z=1:length(mu)
        Z(z,:,:) = mu{z}';
      end
      mu = Z;
     end

    end

    function [mu,V,J,VP] = smooth(obj,mu1,V1,Y)
      % Kalman smoother; uses filtered means and variances

      T = numel(V1);

      nandim=[];
      dim=[];
      if exist('Y','var')
          for i=1:size(Y,1)
              if isnan(Y(i,1))
                  nandim=[nandim i];
              else
                  dim=[dim,i];
              end
          end

          A = obj.A(nandim,nandim);
          C = obj.C(:,nandim);
          Q = obj.Q(nandim,nandim);
      else

          A   = obj.A;
          Q   = obj.Q;
          C   = obj.C;
      end
      % P(t,t-1) at horizon
      P = A * V1{T-1} * A' + Q;
      PC = P * C';
      K = PC / (C * PC + obj.R);
      VP=cell(1,T);
      KC=K*C;
      VP{T}=(eye(size(KC))-KC)*A*V1{T-1};

      % RTS equations

      mu = mu1;
      V = V1;
      J = cell(1,T-1); % needed to compute transition matrix A in EM
      for n=T:-1:2

          P = A * V1{n-1} * A' + Q;
        if det(P)<0
            disp('bad predictive covariance')
            lambda=max(svd(P));
            P=P+(lambda+obj.epsilon)*eye(size(P));
        end
        if ~all(V{n}(:)==0)

          J{n-1} = (V1{n-1} * A') / P;

        if numel(nandim)>0
          mu(:,n-1) = mu1(:,n-1) + J{n-1} * (mu(:,n) - A * mu1(:,n-1)-obj.A(nandim,dim)*Y(dim,n-1));
        else
          mu(:,n-1) = mu1(:,n-1) + J{n-1} * (mu(:,n) - A * mu1(:,n-1));
        end

          V{n-1} = V1{n-1} + J{n-1} * (V{n} - P)*J{n-1}';

          if n<T
            VP{n}=V1{n}*J{n-1}'+J{n}*(VP{n+1}-A*V1{n})*J{n-1}';
          end

        else

          % fix when variance is zero (deterministic state)
          J{n} = zeros(size(V{n}));
          mu(:,n) = mu1(:,n);
          V{n} = zeros(size(V{n}));
          VP{n} = zeros(size(V{n}));

        end

        % symmetricize and make positive semidefinite
        V{n} = (V{n} + V{n}') ./ 2;
        V{n} = V{n} + obj.epsilon*eye(size(V{n}));

      end

    end

    function [mu,V,LL] = filter(obj,X,varargin)
      % Kalman filter algorithm; returns mean mu and covariance matrix V of
      % the states, the log likelihood and the observations Xobs with missing
      % values imputed (the latter remains to be done)

      nandim=[];
      dim=[];
      if ~isempty(varargin)
          Y=varargin{1};
      end
      if exist('Y','var')
          for i=1:size(Y,1)
              if isnan(Y(i,1))
                  nandim=[nandim i];
              else
                  dim=[dim,i];
              end
          end

          A = obj.A(nandim,nandim);
          C = obj.C(:,nandim);
          mu0 = obj.mu0(nandim);
          V0 = obj.V0(nandim,nandim);
          R = obj.R;
          Q = obj.Q(nandim,nandim);
          X1=X;
          X=X-obj.C(:,dim)*Y(dim,:);
      else
          A = obj.A;
          C = obj.C;
          mu0 = obj.mu0;
          V0 = obj.V0;
          R = obj.R;
          Q = obj.Q;
          Y=zeros(size(A,1),size(X,2));
          X1=X;
      end

      K = size(A,1);
      T = size(X,2);

      mu = nan(K,T);
      V = cell(1,T);
      I = eye(K);

      Xpred=mu;
      Vpred=V;

      % take measurements into account
      if isempty(X) || all(isnan(X(:,1)))

        mu(:,1) = mu0;
        V{1} = V0;

      else
        Xpred(:,1)=mu0;
        Vpred{1}=V0;
        obs = ~isnan(X(:,1));
        Cb = C(obs,:); % emission matrix for observed measurements

        % Kalman gain matrix
        K = (V0 * Cb') / (Cb * V0 * Cb' + R(obs,obs));

        mu(:,1) = mu0 + K * (X(obs,1) - Cb * mu0);
        V{1}  = (I - K * Cb) * V0;

      end

      % symmetricize and make positive semidefinite
      V{1} = (V{1} + V{1}') ./ 2;
      V{1} = V{1} + obj.epsilon*eye(size(V{1}));

      for t=2:T

        if numel(nandim)>0
          AM = A * mu(:,t-1)+obj.A(nandim,dim)*Y(dim,t-1);
        else
            AM = A * mu(:,t-1);
        end
        P = A * V{t-1} * A' + Q;

        if isempty(X) || all(isnan(X(:,t)))

          mu(:,t) = AM;
          V{t} = P;

        else
            Xpred(:,t)=AM;
            Vpred{t}=P;
          obs = ~isnan(X(:,t));
          Cb = C(obs,:); % emission matrix for observed measurements

          PC = P * Cb';

          K = PC / (Cb * PC + R(obs,obs));

          mu(:,t) = AM + K * (X(obs,t) - Cb * AM);
          V{t}  = (I - K * Cb) * P;

        end

        % symmetricize and make positive semidefinite
        V{t} = (V{t} + V{t}') ./ 2;
        V{t} = V{t} + obj.epsilon*eye(size(V{t}));

      end

      % compute log likelihood using error decomposition ignoring constant terms
      if nargout >= 3
          V1 = cell(1,T);
          K=size(obj.A,1);
          Y1=zeros(K,T);
          Y1(dim,:)=Y(dim,:);
          if numel(nandim)<1
              Y1=mu;
              V1=V;
          else
              Y1(nandim,:)=Xpred;
              i=1;
              V1{i}=zeros(K);
              V1{i}(dim,dim)=obj.Q(dim,dim);
              V1{i}(dim,nandim)=obj.Q(dim,nandim);
              V1{i}(nandim,dim)=obj.Q(nandim,dim);
              V1{i}(nandim,nandim)=Vpred{i};
              for i=2:T
                  V1{i}=zeros(K);
                  V1{i}(dim,dim)=obj.Q(dim,dim);
                  V1{i}(dim,nandim)=obj.Q(dim,nandim);
                  V1{i}(nandim,dim)=obj.Q(nandim,dim);
                  V1{i}(nandim,nandim)=Vpred{i};
              end
          end
          LL = obj.compute_loglik(X1,Y1,V1);
      end

    end

    function ll = likelihood(obj,X)
      % returns log likelihood of a sequence of observations

      if ~iscell(X)
        if ndims(X) ~= 3, error('input should be of size trials x features x time'); end
        Xt = cell(1,size(X,1));
        for c=1:length(Xt)
          Xt{c} = squeeze(X(c,:,:));
        end
        X = Xt; clear Xt;
      end

      ll = nan(length(X),1);

      for c=1:length(X)

        [mu,V,ll(c)] = obj.filter(X{c});

      end

    end

  end

  methods(Access=protected)

    function [G1,G2,G3,G4,G5,G6,mu0,V0,LL] = Estep(obj,X,Y)

      K = size(Y,1);
      M = size(X,1);
      T = size(X,2);

      nandim=[];
      dim=[];
      if exist('Y','var')
          for i=1:K
              if isnan(Y(i,1))
                  nandim=[nandim i];
              else
                  dim=[dim,i];
              end
          end
      end

      if any(isnan(Y(:))) % hidden state

          [mu1,V1,LL] = obj.filter(X,Y);
          [Y1,V1,J1,VP1] = obj.smooth(mu1,V1,Y);

          Y(nandim,:)=Y1;
          VP=cell(1,T);
          V = cell(1,T);
          %J = cell(1,T-1);
          for i=1:T
              if i>1
                  VP{i}=zeros(K);
                  VP{i}(nandim,nandim)=VP1{i};
              end
              V{i}=zeros(K);
              V{i}(dim,dim)=V{i}(dim,dim)+obj.epsilon*eye(length(dim));
              V{i}(nandim,nandim)=V1{i};
              %if i<T
              %    J{i}=zeros(K);
              %    J{i}(nandim,nandim)=J1{i};
              %end
          end

      else

        V = repmat({0},size(Y));
        %J = repmat({0},[size(Y,1) size(Y,2)-1]);
        VP = repmat({0},size(Y));
        LL = nan;

      end

      G1 = zeros(K,K);
      G4 = zeros(K,K);
      for t=1:T

        G1  = G1 + V{t} + Y(:,t)*Y(:,t)';

        if t>1

          G4 = G4 + Y(:,t)*Y(:,t-1)' + VP{t};
          %G4 = G4 + J{t-1}*V{t} + Y(:,t)*Y(:,t-1)'; % bishop

        end

      end

      G2 = G1 - Y(:,1) * Y(:,1)' - V{1};
      G3 = G1 - Y(:,T) * Y(:,T)' - V{T};

      mu0 = Y(:,1);
      V0 = V{1};

      % deal with partially observed observations
      if any(isnan(X(:)))

        C = obj.C;
        R = obj.R;

        G5 = zeros(M,M);
        for t=1:T

          XX = X(:,t) * X(:,t)';
          nidx = isnan(XX(:));
          if any(nidx)
            E = X(:,t) * Y(:,t)' * C' + R;
            XX(nidx) = E(nidx);
          end
          nidx = isnan(XX(:));
          if any(nidx)
            E = C * Y(:,t) * X(:,t)' + R;
            XX(nidx) = E(nidx);
          end
          nidx = isnan(XX(:));
          if any(nidx)
            E = C * (V{t} + Y(:,t)*Y(:,t)') * C' + R;
            XX(nidx) = E(nidx);
          end

          G5 =  G5 + XX;

        end

        G6 = zeros(M,K);
        for t=1:T

          XY = X(:,t) * Y(:,t)';
          EXY = C * (V{t} + Y(:,t)*Y(:,t)');
          nidx = isnan(XY(:));
          XY(nidx) = EXY(nidx);
          G6 = G6 + XY;

        end

      else
        G5 = X * X';
        G6 = X * Y';
      end

    end

    function LL = compute_loglik(obj,X,Y,V)

      T = size(X,2);

      C = obj.C;
      R = obj.R;

      LL = 0;
      for t=1:T

        obs = ~isnan(X(:,t));

        e = X(obs,t) - C(obs,:)*Y(:,t); % innovation (measurement error)
        S = C(obs,:)*V{t}*C(obs,:)' + R(obs,obs); % covariance of the innovation
        LL = LL - log(det(S)) ./ 2 - ((e' / S) * e) ./ 2;
        if ~isreal(LL)
            LL=real(LL);
            disp('bad likelihood')
        end
        LL = LL - (numel(obs)/2) * log(2*pi); %length of the feature vector

      end

    end

  end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function f = compose(f, op, g, data, pref)
%COMPOSE   Composition of UNBNDFUN objects.
%   H = COMPOSE(F, OP) returns a UNBNDFUN representing OP(F) where F is also a
%   UNBNDFUN object, and OP is a function handle.
%
%   H = COMPOSE(F, OP, G) returns OP(F, G) where F and G are UNBNDFUN objects,
%   and OP is a function handle. Note that F and G are assumed to have the same
%   domain. No warning will be given/error thrown if that is not the case, but
%   the output of the method will not be useful.
%
%   H = COMPOSE(F, G) returns a FUN representing G(F), where at least one of F
%   and G is an UNBNDFUN object. If the range of F is not included in the domain
%   of G, an error is thrown. Note that the domain of H will be the same as the
%   domain of F.
%
%   H = COMPOSE(F, OP, G, DATA, PREF) or H = COMPOSE(F, OP, [], DATA, PREF),
%   where F (and G) are UNBNDFUN objects, and OP is a function handle, uses the
%   constructor data in the structure DATA and the options passed by the
%   CHEBFUNPREF or preferences structure PREF to build the returned UNBNDFUN.
%   In particular, if the underlying tech class supports it, one can use PREF
%   to alter the constructor's behavior to take advantage of the fact that F
%   (and possibly OP or G) has additional structure beyond just being an
%   UNBNDFUN object.

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

% Parse inputs:
if ( nargin == 2 )

    if ( isa(op, 'unbndfun') )
        % Composition of two UNBNDFUN objects ( out = g(f) ):

        % OP is a UNBNDFUN! Rename it to g for clarity:
        g = op;

        % For composition to work with two FUN objects, the range of F must be
        % within the domain of G. However, we allow this method to assume that
        % it is the case.

        % When we compose F and G, we need to find the appropriate part of the
        % domain of G on which the range of F lies. Since we then want to work
        % with the ONEFUN of G, we apply the inverse map for the domain of G to
        % find the (sub)interval on [-1, 1] on which the range of F (and hence
        % its ONEFUN) lives.
        fMapped = g.mapping.Inv(f.onefun);

        % Try to do the composition, will only work if the range of F lies in
        % the domain of G.
        f.onefun = compose(fMapped, g.onefun);

    else
        % Standard composition ( out = OP(F) )

        % OP is an operator!
        f.onefun = compose(f.onefun, op);
    end

else
    % out = OP(F,G).

    if ( nargin < 5 )
        pref = chebfunpref();
    else
        pref = chebfunpref(pref);
    end

    if ( nargin < 4 )
        data = struct();
    end

    % Call ONEFUN/COMPOSE():
    if ( isempty(g) )
        f.onefun = compose(f.onefun, op, g, data, pref);
    else
        f.onefun = compose(f.onefun, op, g.onefun, data, pref);
    end

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function c = cardinal_cos ( j, m, t )

%*****************************************************************************80
%
%% CARDINAL_COS evaluates the J-th cardinal cosine basis function.
%
%  Discussion:
%
%    The base points are T(I) = pi * I / ( M + 1 ), 0 <= I <= M + 1.
%    Basis function J is 1 at T(J), and 0 at T(I) for I /= J
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    John Boyd,
%    Exponentially convergent Fourier-Chebyshev quadrature schemes on
%    bounded and infinite intervals,
%    Journal of Scientific Computing,
%    Volume 2, Number 2, 1987, pages 99-109.
%
%  Parameters:
%
%    Input, integer J, the index of the basis function.
%    0 <= J <= M + 1.
%
%    Input, integer M, indicates the size of the basis set.
%
%    Input, real T(*), one or more points in [0,pi] where the
%    basis function is to be evaluated.
%
%    Output, real C, the value of the function at T.
%
  if ( j == 0 || j == m + 1 )
    cj = 2.0;
  else
    cj = 1.0;
  end

  tj = pi * j / ( m + 1 );

  c = ( -1.0 ) ^ mod ( j + 1, 2 ) ...
    * sin ( t ) ...
    .* sin ( ( m + 1 ) * t ) ...
    / cj ...
    / ( m + 1 ) ...
    ./ ( cos ( t ) - cos ( tj ) );

  i = find ( abs ( t - tj ) < eps );
  c(i) = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef ( Hidden, Sealed ) FigureObserver < handle
    %uix.FigureObserver  Figure observer
    %
    %  A figure observer raises an event FigureChanged when the figure
    %  ancestor of a subject changes.

    %  Copyright 2014-2015 The MathWorks, Inc.
    %  $Revision: 1435 $ $Date: 2016-11-17 17:50:34 +0000 (Thu, 17 Nov 2016) $

    properties( SetAccess = private )
        Subject % subject
        Figure % figure ancestor
    end

    properties( Access = private )
        PreSetListeners % listeners to Parent PreGet
        PostSetListeners % listeners to Parent PreGet
        OldFigure = gobjects( 0 ) % previous figure ancestor
    end

    events( NotifyAccess = private )
        FigureChanged
    end

    methods

        function obj = FigureObserver( subject )
            %uix.FigureObserver  Create figure observer
            %
            %  o = uix.FigureObserver(s) creates a figure observer for the
            %  subject s.

            % Check
            validateattributes( subject, {'matlab.graphics.Graphics'}, ...
                {'scalar'}, '', 'subject' )

            % Store subject
            obj.Subject = subject;

            % Set up object
            obj.update()

        end % constructor

    end % structors

    methods( Access = private )

        function update( obj )
            %update  Update listeners and Figure property

            % Create fresh listeners
            obj.PreSetListeners = event.proplistener.empty( [1 0] ); % clear
            obj.PostSetListeners = event.proplistener.empty( [1 0] ); % clear
            o = obj.Subject;
            while ~isempty( o ) && ~isa( o, 'matlab.ui.Figure' )
                obj.PreSetListeners(end+1) = event.proplistener( o, ...
                    findprop( o, 'Parent' ), 'PreSet', @obj.onParentPreSet );
                obj.PostSetListeners(end+1) = event.proplistener( o, ...
                    findprop( o, 'Parent' ), 'PostSet', @obj.onParentPostSet );
                o = o.Parent;
            end

            % Store figure
            obj.Figure = o;

        end % update

        function onParentPreSet( obj, ~, ~ )
            %onParentPreSet  Event handler

            % Store old figure
            obj.OldFigure = obj.Figure;

        end % onParentPreSet

        function onParentPostSet( obj, ~, ~ )
            %onParentPostSet  Event handler

            % Update object
            obj.update()

            % Raise event
            oldFigure = obj.OldFigure;
            newFigure = obj.Figure;
            if ~isequal( oldFigure, newFigure )
                notify( obj, 'FigureChanged', ...
                    uix.FigureData( oldFigure, newFigure ) )
            end

            % Clear old figure
            obj.OldFigure = gobjects( 0 );

        end % onParentPostSet

    end % private methods

end % classdef

author: goGPS-Project


text:
function [info] = decode_fif(orig)

% DECODE_FIF is a helper function for real-time processing of Neuromag data. This
% function is used to decode the content of the optional neuromag_fif chunk(s).
%
% See also DECODE_RES4, DECODE_NIFTI1, SAP2MATLAB

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

% check that the required low-level toolbox is available
ft_hastoolbox('mne', 1);

global FIFF
if isempty(FIFF)
  FIFF = fiff_define_constants();
end

if isfield(orig, 'neuromag_header')
  % The binary blob was created on the little-endian Intel Linux acquisition
  % computer, whereas the default for fiff files is that they are stored in
  % big-endian byte order. MATLAB is able to swap the bytes on the fly by specifying
  % 'le" or "be" to fopen. The normal MNE fiff_open function assumes that it is big
  % endian, hence here we have to open it as little endian.

  filename = tempname;
  % write the binary blob to disk, byte-by-byte to avoid any swapping between little and big-endian content
  F = fopen_or_error(filename, 'w');
  fwrite(F, orig.neuromag_header, 'uint8');
  fclose(F);
  % read the content of the file using the standard reading functions
  [info, meas] = read_header(filename);
  % clean up the temporary file
  delete(filename);
end

% Typically, at the end of acquisition, the isotrak and hpiresult information
% is stored in the neuromag fiff container which can then (offline) be read by
% fiff_read_meas_info. However, for the purpose of head position monitoring
% (see Stolk et al., Neuroimage 2013) during acquisition, this crucial
% information requires to be accessible online. read_isotrak and read_hpiresult
% can extract information from the additionally chunked (neuromag2ft) files.

if isfield(orig, 'neuromag_isotrak')
  filename = tempname;
  % write the binary blob to disk, byte-by-byte to avoid any swapping between little and big-endian content
  F = fopen_or_error(filename, 'w');
  fwrite(F, orig.neuromag_isotrak, 'uint8');
  fclose(F);
  % read the content of the file using the standard reading functions
  [info.dig] = read_isotrak(filename);
  % clean up the temporary file
  delete(filename);
end

if isfield(orig, 'neuromag_hpiresult')
  filename = tempname;
  % write the binary blob to disk, byte-by-byte to avoid any swapping between little and big-endian content
  F = fopen_or_error(filename, 'w');
  fwrite(F, orig.neuromag_hpiresult, 'uint8');
  fclose(F);
  % read the content of the file using the standard reading functions
  [info.dev_head_t, info.ctf_head_t] = read_hpiresult(filename);
  % clean up the temporary file
  delete(filename);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [info, meas] = read_header(filename)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% open and read the file as little endian
[fid, tree] = fiff_open_le(filename); % open as little endian
[info, meas] = fiff_read_meas_info(fid, tree);
fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dig] = read_isotrak(filename)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global FIFF

% open the isotrak file (big endian)
% (typically stored in meas_info dir during acquisition, no fif extension required)
[fid, tree] = fiff_open(filename);

% locate the Polhemus data
isotrak = fiff_dir_tree_find(tree,FIFF.FIFFB_ISOTRAK);

dig=struct('kind',{},'ident',{},'r',{},'coord_frame',{});
coord_frame = FIFF.FIFFV_COORD_HEAD;
if length(isotrak) == 1
  p = 0;
  for k = 1:isotrak.nent
    kind = isotrak.dir(k).kind;
    pos  = isotrak.dir(k).pos;
    if kind == FIFF.FIFF_DIG_POINT
      p = p + 1;
      tag = fiff_read_tag(fid,pos);
      dig(p) = tag.data;
    else
      if kind == FIFF.FIFF_MNE_COORD_FRAME
        tag = fiff_read_tag(fid,pos);
        coord_frame = tag.data;
      elseif kind == FIFF.FIFF_COORD_TRANS
        tag = fiff_read_tag(fid,pos);
        dig_trans = tag.data;
      end
    end
  end
end
for k = 1:length(dig)
  dig(k).coord_frame = coord_frame;
end

if exist('dig_trans','var')
  if (dig_trans.from ~= coord_frame && dig_trans.to ~= coord_frame)
    clear('dig_trans');
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dev_head_t, ctf_head_t] = read_hpiresult(filename)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global FIFF

% open the hpiresult file (big endian)
% (typically stored in meas_info dir during acquisition, no fif extension required)
[fid, tree] = fiff_open(filename);

% locate the transformation matrix
dev_head_t=[];
ctf_head_t=[];
hpi_result = fiff_dir_tree_find(tree,FIFF.FIFFB_HPI_RESULT);
if length(hpi_result) == 1
  for k = 1:hpi_result.nent
    kind = hpi_result.dir(k).kind;
    pos  = hpi_result.dir(k).pos;
    if kind == FIFF.FIFF_COORD_TRANS
      tag = fiff_read_tag(fid,pos);
      cand = tag.data;
      if cand.from == FIFF.FIFFV_COORD_DEVICE && ...
          cand.to == FIFF.FIFFV_COORD_HEAD
        dev_head_t = cand;
      elseif cand.from == FIFF.FIFFV_MNE_COORD_CTF_HEAD && ...
          cand.to == FIFF.FIFFV_COORD_HEAD
        ctf_head_t = cand;
      end
    end
  end
end


text:
function pdf = inverse_gaussian_pdf ( x, a, b )

%*****************************************************************************80
%
%% INVERSE_GAUSSIAN_PDF evaluates the Inverse Gaussian PDF.
%
%  Discussion:
%
%    The Inverse Gaussian PDF is also known as the Wald PDF
%    and the Inverse Normal PDF.
%
%    PDF(X)(A,B)
%      = SQRT ( B / ( 2 * PI * X**3 ) )
%        * EXP ( - B * ( X - A )**2 / ( 2.0D+00 * A**2 * X ) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%    0.0 < X
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < A,
%    0.0 < B.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x <= 0.0 )
    pdf = 0.0;
  else
    pdf = sqrt ( b / ( 2.0 * pi * x^3 ) ) * ...
      exp ( - b * ( x - a )^2 / ( 2.0 * a^2 * x ) );
  end

  return
end


text:
function [ objects] = room3DFromHyp( points, floor_height )
%ROOM3DFROMHYP Popup room layout hypothesis to 3D
%   points: room corner in sequence in room layout
%   floor_height: assumed floor height, e.g. -160

% uv = xyz2uvN(points);
% high_height = floor_height * tan(uv(1:4,2))./tan(uv(5:8,2));
% ceil_height = mean(high_height);
%
% point3D = zeros(8,3);
% for i = 5:8
%     point3D(i,:) = LineFaceIntersection([0 0 floor_height], [0 0 1], [0 0 0], points(i,:));
% end
COUNTER_CLOCKWISE = true;
FITTINGRULES;

[pt I] = sortXYZ(points(1:4,:));
pd = points(5:8,:); pd = pd(I,:);
sort_points = [pd;pt];
ann_points = sort_points(POINTMAPPING{3},:);

objects = struct('out_points_w',[],'name',[],'type',[],'points',[],'x_w',[]);
[ out_point, x, fval, succ ] = initial_cuboid_3D( ann_points, 3, floor_height, 3, true );

objects.points = ann_points;
objects.type = 3;
objects.name = 'room';
objects.out_points_w = out_point;
objects.x_w = x;

end

author: zouchuhang

repo: LayoutNet


text:
%DEMGAUSS Demonstrate sampling from Gaussian distributions.
%
%	Description
%
%	DEMGAUSS provides a simple illustration of the generation of data
%	from Gaussian distributions. It first samples from a one-dimensional
%	distribution using RANDN, and then plots a normalized histogram
%	estimate of the distribution using HISTP together with the true
%	density calculated using GAUSS.
%
%	DEMGAUSS then demonstrates sampling from a Gaussian distribution in
%	two dimensions. It creates a mean vector and a covariance matrix, and
%	then plots contours of constant density using the function GAUSS. A
%	sample of points drawn from this distribution, obtained using the
%	function GSAMP, is then superimposed on the contours.
%
%	See also
%	GAUSS, GSAMP, HISTP
%

%	Copyright (c) Ian T Nabney (1996-2001)

clc
mean = 2; var = 5; nsamp = 3000;
xmin = -10; xmax = 10; nbins = 30;
disp('Demonstration of sampling from a uni-variate Gaussian with mean')
dstring = [num2str(mean), ' and variance ', num2str(var), '.  ', ...
    num2str(nsamp), ' samples are taken.'];
disp(dstring);
x = mean + sqrt(var)*randn(nsamp, 1);
fh1 = figure;
histp(x, xmin, xmax, nbins);
hold on;
axis([xmin xmax 0 0.2]);
plotvals = linspace(xmin, xmax, 200)';
probs = gauss(mean, var, plotvals);
plot(plotvals, probs, '-r');
xlabel('X')
ylabel('Density')

disp(' ')
disp('Press any key to continue')
pause;
mu = [3 2];
lam1 = 0.5;
lam2 = 5.0;
Sigma = lam1*[1,1]'*[1,1] + lam2*[1,-1]'*[1,-1];
disp(' ')
disp('Demonstration of sampling from a bi-variate Gaussian.  The mean is')
dstring = ['[', num2str(mu(1)), ', ', num2str(mu(2)), ...
      '] and the covariance matrix is'];
disp(dstring)
disp(Sigma);
ngrid = 40;
cmin = -5; cmax = 10;
cvals = linspace(cmin, cmax, ngrid);
[X1, X2] = meshgrid(cvals, cvals);
XX = [X1(:), X2(:)];
probs = gauss(mu, Sigma, XX);
probs = reshape(probs, ngrid, ngrid);

fh2 = figure;
contour(X1, X2, probs, 'b');
hold on

nsamp = 300;
dstring = [num2str(nsamp), ' samples are generated.'];
disp('The plot shows the sampled data points with a contour plot of their density.')
samples = gsamp(mu, Sigma, nsamp);
plot(samples(:,1), samples(:,2), 'or');
xlabel('X1')
ylabel('X2')
grid off;

disp(' ')
disp('Press any key to end')
pause;
close(fh1);
close(fh2);
clear all;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function F = erf(F, varargin)
%ERF   Error function of a CHEBFUN.
%   ERF(X) is the error function of the real-valued CHEBFUN X.
%
%   The error function is defined as:
%       erf(X)(s) = 2/sqrt(pi) * integral from 0 to X(s) of exp(-t^2) dt.
%
% See also ERFC, ERFCX, ERFINV, ERFCINV.

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

% Input must be real:
if ( ~isreal(F) )
    error('CHEBFUN:CHEBFUN:erf:notreal', 'Input must be real.');
end

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

end


text:
function [ a, seed ] = r8cb_random ( n, ml, mu, seed )

%*****************************************************************************80
%
%% R8CB_RANDOM randomizes a R8CB matrix.
%
%  Discussion:
%
%    The R8CB storage format is appropriate for a compact banded matrix.
%    It is assumed that the matrix has lower and upper bandwidths ML and MU,
%    respectively.  The matrix is stored in a way similar to that used
%    by LINPACK and LAPACK for a general banded matrix, except that in
%    this mode, no extra rows are set aside for possible fillin during pivoting.
%    Thus, this storage mode is suitable if you do not intend to factor
%    the matrix, or if you can guarantee that the matrix can be factored
%    without pivoting.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative, and no greater than N-1.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real A(ML+MU+1,N), the R8CB matrix.
%
%    Output, integer SEED, an updated seed for the random number generator.
%

%
%  Set the entries that correspond to matrix elements.
%
  for j = 1 : n

    ilo = max ( 1, j - mu );
    ihi = min ( n, j + ml );

    for i = j - mu : 0
      a(i-j+mu+1,j) = 0.0;
    end

    for i = ilo : ihi
      [ a(i-j+mu+1,j), seed ] = r8_uniform_01 ( seed );
    end

    for i = n + 1 : j + ml
      a(i-j+mu+1,j) = 0.0;
    end

  end

  return
end


text:
function [ y, m, d ] = day_borrow_julian ( y, m, d )

%*****************************************************************************80
%
%% DAY_BORROW_JULIAN borrows days from months in a Julian date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, integer M, integer D, a year, month, and day
%    representing a date.  On input, D might be negative.
%
%    Output, integer Y, integer M, integer D, a year, month, and day
%    representing a date.  On output,
%    M should have decreased by one month, and D gone up by the
%    number of days in the month we "cashed in".  Y may be affected
%    if the input value of M was 1.
%
  while ( d <= 0 )

    m = m - 1;

    [ y, m ] = month_borrow_julian ( y, m );

    days = month_length_julian ( y, m );

    d = d + days;

  end

  return
end


% params
totalNum = 6159264;
bValue = 1000;
lamda = 0;
qFile = '/home/yinhuan/Data/mapModel/yq/weightVector.txt';
visFilesDir = '/home/yinhuan/Data/mapModel/yq/visMatrix/';
maxQ = 179;
minQ = 1;

splitLength_0 = 50;
saveResultsDir_0 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/0/';
saveReIndexFile_0 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/index/0.txt';
saveNewQFile_0 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/weight/0.txt';

splitLength_1 = 100;
saveResultsDir_1 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/1/';
saveReIndexFile_1 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/index/1.txt';
saveNewQFile_1 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/weight/1.txt';

splitLength_2 = 200;
saveResultsDir_2 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/2/';
saveReIndexFile_2 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/index/2.txt';
saveNewQFile_2 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/weight/2.txt';

splitLength_3 = 400;
saveResultsDir_3 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/3/';
saveReIndexFile_3 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/index/3.txt';
saveNewQFile_3 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/weight/3.txt';

splitLength_4 = 800;
saveResultsDir_4 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/4/';
saveReIndexFile_4 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/index/4.txt';
saveNewQFile_4 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/weight/4.txt';

splitLength_5 = 1600;
saveResultsDir_5 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/5/';
saveReIndexFile_5 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/index/5.txt';
saveNewQFile_5 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/weight/5.txt';

splitLength_6 = 3200;
saveResultsDir_6 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/6/';
saveReIndexFile_6 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/index/6.txt';
saveNewQFile_6 = '/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/weight/6.txt';

%% iteration of optimization

visCells_origin = fromVisDirToCells(visFilesDir);

[ epsilon_soft_0, time_sum_0  ] = Uniform_loopCompress_Cells( lamda, qFile, visCells_origin, splitLength_0, totalNum, bValue, saveResultsDir_0 );
min_cost_0 = get_min_cost(minQ, maxQ, saveResultsDir_0, qFile, lamda, epsilon_soft_0);
[totalNum_0, visCells_0] = deleteZeroPoints_Cells( qFile, visCells_origin, saveResultsDir_0, saveReIndexFile_0, saveNewQFile_0 );

[ epsilon_soft_1, time_sum_1  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_0, visCells_0, splitLength_1, totalNum_0, bValue, saveResultsDir_1 );
min_cost_1 = get_min_cost(minQ, maxQ, saveResultsDir_1, saveNewQFile_0, lamda, epsilon_soft_1);
[totalNum_1, visCells_1] = deleteZeroPoints_Cells( saveNewQFile_0, visCells_0, saveResultsDir_1, saveReIndexFile_1, saveNewQFile_1 );

[ epsilon_soft_2, time_sum_2  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_1, visCells_1, splitLength_2, totalNum_1, bValue, saveResultsDir_2 );
min_cost_2 = get_min_cost(minQ, maxQ, saveResultsDir_2, saveNewQFile_1, lamda, epsilon_soft_2);
[totalNum_2, visCells_2] = deleteZeroPoints_Cells( saveNewQFile_1, visCells_1, saveResultsDir_2, saveReIndexFile_2, saveNewQFile_2 );

[ epsilon_soft_3, time_sum_3  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_2, visCells_2, splitLength_3, totalNum_2, bValue, saveResultsDir_3 );
min_cost_3 = get_min_cost(minQ, maxQ, saveResultsDir_3, saveNewQFile_2, lamda, epsilon_soft_3);
[totalNum_3, visCells_3] = deleteZeroPoints_Cells( saveNewQFile_2, visCells_2, saveResultsDir_3, saveReIndexFile_3, saveNewQFile_3 );

[ epsilon_soft_4, time_sum_4  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_3, visCells_3, splitLength_4, totalNum_3, bValue, saveResultsDir_4 );
min_cost_4 = get_min_cost(minQ, maxQ, saveResultsDir_4, saveNewQFile_3, lamda, epsilon_soft_4);
[totalNum_4, visCells_4] = deleteZeroPoints_Cells( saveNewQFile_3, visCells_3, saveResultsDir_4, saveReIndexFile_4, saveNewQFile_4 );

[ epsilon_soft_5, time_sum_5  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_4, visCells_4, splitLength_5, totalNum_4, bValue, saveResultsDir_5 );
min_cost_5 = get_min_cost(minQ, maxQ, saveResultsDir_5, saveNewQFile_4, lamda, epsilon_soft_5);
[totalNum_5, visCells_5] = deleteZeroPoints_Cells( saveNewQFile_4, visCells_4, saveResultsDir_5, saveReIndexFile_5, saveNewQFile_5 );

[ epsilon_soft_6, time_sum_6  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_5, visCells_5, splitLength_6, totalNum_5, bValue, saveResultsDir_6 );
min_cost_6 = get_min_cost(minQ, maxQ, saveResultsDir_6, saveNewQFile_5, lamda, epsilon_soft_6);
[totalNum_6, visCells_6] = deleteZeroPoints_Cells( saveNewQFile_5, visCells_5, saveResultsDir_6, saveReIndexFile_6, saveNewQFile_6 );

% save the results

compressIndex_5 = anal_reindex_last(saveResultsDir_6, saveReIndexFile_5);
compressIndex_4 = anal_reindex_middle(compressIndex_5, saveReIndexFile_4);
compressIndex_3 = anal_reindex_middle(compressIndex_4, saveReIndexFile_3);
compressIndex_2 = anal_reindex_middle(compressIndex_3, saveReIndexFile_2);
compressIndex_1 = anal_reindex_middle(compressIndex_2, saveReIndexFile_1);
compressIndex_0 = anal_reindex_middle(compressIndex_1, saveReIndexFile_0);

dlmwrite('/home/yinhuan/Data/mapModel/yq/iter_b_1000_0/compressIndex_final.txt', compressIndex_0, 'precision', '%d');


text:
function eclipVec=ICRS2Eliptic(vec,TT1,TT2,method)
%%ICRS2ECLIPTIC Convert a location vector from the International
%               Celestial Reference System (ICRS) to eliptic coordinates
%               either using the IAU 2006 precession model or the Vondrak
%               400 millennia precession model.
%
%INPUTS: x The NXnumVec collection of vectors in the ICRS to convert. N
%          can be 2, or 3. If the vectors are 2D, then they are assumed to
%          be azimuth and elevation in radians. 3D vectors are assumed to
%          be Cartesian position.
% Jul1, Jul2 Two parts of a Julian date given in terrestrial time (TT).
%          The units of the date are days. The full date is the sum of
%          both terms. The date is broken into two parts to provide more
%          bits of precision. It does not matter how the date is split.
%   method An optional parameter specifying which algorithm is to be used.
%          Possible values are
%          0 (The default if omitted or an empty matrix is passed) Use the
%            IAU 2006 precession model.
%          1 Use the long-term (Vondrak) precession model.
%
%OUTPUTS: xG The vectors rotated into the ecliptic coordinate system. If
%            the input was 2D azimuth and elevation, the output will be
%            the same. If the input was Cartesian, then the output will be
%            Cartesian.
%
%This function is a Matlab interface for the relevant functions in the
%International Astronomical Union's (IAU) Standard's of Fundamental
%Astronomy library.
%
%The ecliptic is defined in the IERS Conventions [1] to be the "the
%plane perpendicular to the mean heliocentric orbital angular momentum
%vector of the Earth-Moon barycentre in the BCRS".
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%eclipVec=ICRS2Eliptic(vec,TT1,TT2,method);
%
%REFERENCES:
%[1] G. Petit and B. Luzum, IERS Conventions (2010), International Earth
%    Rotation and Reference Systems Service Std. 36, 2010.
%
%July 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

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:
classdef Logger < handle
% A class that writes messages to log files (and/or standard out).
% To use the class, put the following two lines in your script:
%   global logger
%   logger = Logger();
% and then add and remove logfiles as required e.g.
%   logger.add_stdout(logger.All);
%   logger.add_filename('errorfile.txt',logger.Error);
%   logger.log(logger.All,'1 + 2 = %d\n',3);
%   logger.log(logger.Error,'An error occurred.\n');
%   logger.rm_filename('errorfile.txt');
%   fid = fopen('warningfile.txt','w');
%   logger.add_fid(fid);
%   logger.log(logger.Warning,'This is a warning.\n');
%   logger.clear()
%   fclose(fid);

properties (Access=private)
  fidlist = []
  fidlevels = [];
  filenamelist = {}
  filenamefids = []
  filenamelevels = []
end

properties (Access=public,Constant=true)
  None = 0
  Major_Error = 1
  Error = 2
  Warning = 4
  Minor_Warning = 8
  Info = 16
  Minor_Info = 32
  Debug_1 = 64
  Debug_2 = 128
  All = 255;
end

methods (Access=public)
  log(logger,level,message,varargin)
  add_fid(logger,fid,level)
  rm_fid(logger,fid)
  clear(logger)
  add_filename(logger,fname,level)
  rm_filename(logger,fname)
  add_stdout(logger,level)
  rm_stdout(logger)
  add_stderr(logger,level)
  rm_stderr(logger)
end

end

author: nesl


text:
function fg = fgCreate(varargin)
% Create a fiber group structure with default values
%
%   fg = fgCreate;
%
% Fields that can be set are
%   name = 'FG-1';
%   thickness = -0.5;        % Means ...
%   visibleFlag = 1;
%   colorRgb = [20 90 200];
%   seeds = [];
%   seedRadius = 0;
%   seedVoxelOffsets = [];
%   params = {};
%   fibers = fibers;
%   query_id = -1;
%   coordspace = struct('acpc', eye(4), 'img', [], 'mni', []);
%
% Example:
%   fg = fgCreate;
%   fg = fgCreate('name','myGroup');
%   fg = fgCreate('name','myGroup','colorRgb',[255 255 255]);
%   fg(2) = fgCreate;
%
% See also: fgSet, fgGet
%
% (c) Stanford VISTA Team

%% Set up the defaults
name = 'FG-1';
type = 'fibergroup';

thickness = -0.5;
visibleFlag = 1;
color = [20 90 200];
fibers = {};

fg.name = name;
fg.type = type;
fg.colorRgb = color;
fg.thickness = thickness;
fg.visible = visibleFlag;
fg.seeds = [];
fg.seedRadius = 0;
fg.seedVoxelOffsets = [];
fg.params = {};
fg.fibers = fibers;
fg.query_id = -1;

% Identifiers for the fibers in fg.fibers
% This is a list of possible names (e.g., arcuate, optic radiation) for
% each of the fibers in the cell array, fg.fibers{ii}.  There might be
% something like 20 different names when it is the Mori groups, for
% example.
fg.fiberNames = [];

% A numerical index, normally of the same length as fg.fibers, where each
% index points to the name associated with that fiber.  So, for the ith
% fiber cell array, fg.fibers{ii}, the index is fg.fiberNumbers(ii) and the
% name is fg.fiberNames(fg.fiberIndex(ii));
% The fiberIndex normally has the same length as fg.fibers.  There should
% be a smaller number of fiberNames.  To indicate that you don't know which
% group a fiber belongs to, use a non positive integer value, say -1.
fg.fiberIndex = [];

% Which space are the fiber coordinates represented in?
% Image space ('img') are basically indices into the data voxels.
% ACPC ('acpc') space is relative to the middle of the acpc.
% MNI ('mni') space is another option.
% This struct contains the transformations required to get from the
% fiber coordinates to the requested coordinate space (if this xform is
% known). If the fiber coordinates are given in a particular coordinate
% space. That matrix for that coordspace will be the identity matrix.
% Unknown transformations (that have not been set in any way) will be set
% to an empty matrix. The default space is 'acpc':
fg.coordspace = struct('acpc', eye(4), 'img', [], 'mni', []);

% If the user has set some inputs, overwrite the defaults here
fg = mrVarargin(fg,varargin);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function m = max(G)
% return maximum value

m = max([G.points]);

author: mtex-toolbox

repo: mtex

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


text:
function [ num_int, p ] = circles_imp_int_2d ( r1, center1, r2, center2 )

%*****************************************************************************80
%
%% CIRCLES_IMP_INT_2D: finds the intersection of two implicit circles in 2D.
%
%  Discussion:
%
%    Two circles can intersect in 0, 1, 2 or infinitely many points.
%
%    The 0 and 2 intersection cases are numerically robust; the 1 and
%    infinite intersection cases are numerically fragile.  The routine
%    uses a tolerance to try to detect the 1 and infinite cases.
%
%    An implicit circle in 2D satisfies the equation:
%
%      ( X - CENTER(1) )**2 + ( Y - CENTER(2) )**2 = R**2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R1, the radius of the first circle.
%
%    Input, real CENTER1(2,1), the center of the first circle.
%
%    Input, real R2, the radius of the second circle.
%
%    Input, real CENTER2(2,1), the center of the second circle.
%
%    Output, integer NUM_INT, the number of intersecting points found.
%    NUM_INT will be 0, 1, 2 or 3.  3 indicates that there are an infinite
%    number of intersection points.
%
%    Output, real P(2,2), if NUM_INT is 1 or 2,
%    the coordinates of the intersecting points.
%
  dim_num = 2;

  tol = eps;

  p(1:dim_num,1:2) = 0.0;;
%
%  Take care of the case in which the circles have the same center.
%
  t1 = ( abs ( center1(1,1) - center2(1,1) ) ...
       + abs ( center1(2,1) - center2(2,1) ) ) / 2.0;

  t2 = ( abs ( center1(1,1) ) + abs ( center2(1,1) ) ...
       + abs ( center1(2,1) ) + abs ( center2(2,1) ) + 1.0 ) / 5.0;

  if ( t1 <= tol * t2 )

    t1 = abs ( r1 - r2 );
    t2 = ( abs ( r1 ) + abs ( r2 ) + 1.0 ) / 3.0;

    if ( t1 <= tol * t2 )
      num_int = 3;
    else
      num_int = 0;
    end

    return

  end

  distsq = sum ( ( center1(1:2,1) - center2(1:2,1) ).^2 );

  root = 2.0 * ( r1 * r1 + r2 * r2 ) * distsq - distsq * distsq ...
    - ( r1 - r2 ) * ( r1 - r2 ) * ( r1 + r2 ) * ( r1 + r2 );

  if ( root < -tol )
    num_int = 0;
    return
  end

  sc1 = ( distsq - ( r2 * r2 - r1 * r1 ) ) / distsq;

  if ( root < tol )
    num_int = 1;
    p(1:dim_num,1) = center1(1:dim_num,1) ...
      + 0.5 * sc1 * ( center2(1:dim_num,1) - center1(1:dim_num,1) );
    return
  end

  sc2 = sqrt ( root ) / distsq;

  num_int = 2;

  p(1,1) = center1(1,1) + 0.5 * sc1 * ( center2(1,1) - center1(1,1) ) ...
                        - 0.5 * sc2 * ( center2(2,1) - center1(2,1) );
  p(2,1) = center1(2,1) + 0.5 * sc1 * ( center2(2,1) - center1(2,1) ) ...
                        + 0.5 * sc2 * ( center2(1,1) - center1(1,1) );

  p(1,2) = center1(1,1) + 0.5 * sc1 * ( center2(1,1) - center1(1,1) ) ...
                        + 0.5 * sc2 * ( center2(2,1) - center1(2,1) );
  p(2,2) = center1(2,1) + 0.5 * sc1 * ( center2(2,1) - center1(2,1) ) ...
                        - 0.5 * sc2 * ( center2(1,1) - center1(1,1) );

  return
end


text:
function openCalcParamWindowSLI(obj)
		%Here a window for the calculation parameters should be created
		%I suggest using the inputsdlg function

		Title = 'SLIDER Decluster ';
		Prompt={'Selection Method:', 'Sli_Type';...
			'Min. Mainshock Magnitude:','MainMag';...
			'Time before Eq (d):', 'BackTime';...
			'Time after Eq (d):', 'ForeTime';...
			'Min. search radius (km):', 'MinRad';...
			'Mag. radius scaler:', 'MagRadScale';...
			'Mag. radius constant:','MagRadConst'};

		LabelList2 =  {	'Magnitude';...
				'Time'};

		%Sli Method
		Formats(1,1).type='list';
		Formats(1,1).style='popupmenu';
		Formats(1,1).items=LabelList2;

		%min MainMag
		Formats(2,1).type='edit';
		Formats(2,1).format='float';
		Formats(2,1).limits = [-99 99];

		%Time before eq
		Formats(3,1).type='edit';
		Formats(3,1).format='float';
		Formats(3,1).limits = [0 9999];

		%Time after eq
		Formats(4,1).type='edit';
		Formats(4,1).format='float';
		Formats(4,1).limits = [0 9999];

		%Min Search Radius
		Formats(5,1).type='edit';
		Formats(5,1).format='float';
		Formats(5,1).limits = [0 9999];

		%Mag Search Radius Scaler
		Formats(6,1).type='edit';
		Formats(6,1).format='float';
		Formats(6,1).limits = [0 9999];

		%Mag Search Radius Constant
		Formats(7,1).type='edit';
		Formats(7,1).format='float';
		Formats(7,1).limits = [0 9999];

		%%%% SETTING DIALOG OPTIONS
		Options.WindowStyle = 'modal';
		Options.Resize = 'on';
		Options.Interpreter = 'tex';
		Options.ApplyButton = 'off';

		%open the dialog window
		[NewParameter,Canceled] = inputsdlg(Prompt,Title,Formats,obj.CalcParameter,Options);

		obj.CalcParameter=NewParameter;

		if Canceled==1
			obj.StartCalc=false;
		else
			obj.StartCalc=true;
		end

end

author: CelsoReyes


text:
%{
    Copyright (C) Cambridge Electronic Design Limited 2014
    Author: James Thompson
    Web: www.ced.co.uk email: james@ced.co.uk, softhelp@ced.co.uk

    This file is part of CEDS64ML, a MATLAB interface to the SON64 library.

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

    CEDS64ML 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 CEDS64ML.  If not, see <http://www.gnu.org/licenses/>.
%}

function [ dTBaseOut ] = CEDS64TimeBase( fhand, dTBaseIn )
%CEDS64TIMEBASE Gets and sets the file time base
%   [ dTBaseOut ] = CEDS64TimeBase( fhand, dTBaseIn )
%   Inputs
%   fhand - An integer handle to an open file
%   dTBaseIn - (Optional) The new time base as a double
%   Outputs
%   dTBaseOut - The old time base as a double, or a number <= 0 if an error.

if (nargin == 1 || nargin == 2) % always get the old time bas
    dTBaseOut = calllib('ceds64int', 'S64GetTimeBase', fhand);
else
    dTBaseOut = -22;
end

% has there been an error?
if dTBaseOut < 0
    return;
end

% if not set the new time bas if we're given one
if (nargin == 2)
    calllib('ceds64int', 'S64SetTimeBase', fhand, dTBaseIn);
end

end

author: brainstorm-tools


text:
%% Copyright (C) 2016-2019, 2022 Colin B. Macdonald
%%
%% 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
%% @defun  zeta (@var{x})
%% @defunx zeta (@var{n}, @var{x})
%% Numerical zeta function.
%%
%% Example:
%% @example
%% @group
%% zeta (1.1)
%%   @result{} ans = 10.584
%% @end group
%% @end example
%%
%% @strong{Note} this function may be slow for large numbers of inputs.
%% This is because it is not a native double-precision implementation
%% but rather the numerical evaluation of the Python @code{mpmath} function
%% @code{zeta}.
%%
%% TODO: The two-argument form is not yet implemented.
%%
%% @seealso{@@sym/zeta}
%% @end defun

function y = zeta (n, x)
  if (nargin ~= 1 && nargin ~= 2)
    print_usage ();
  end

  if (nargin == 1)
    x = n;
    cmd = { 'L = _ins[0]'
            'A = [complex(mpmath.zeta(x)) for x in L]'
            'return A,' };
    c = pycall_sympy__ (cmd, num2cell (x(:)));
    y = reshape (cell2mat (c), size (x));
    return
  end

  error ('zeta: two input arguments not implemented');
end

%!error zeta (1, 2, 3)
%!assert (isnan (zeta (nan)))

%!test
%! x = 1.1;
%! y = sym(11)/10;
%! A = zeta (x);
%! B = double (zeta (y));
%! assert (A, B, -4*eps);

%!test
%! y = [2 3 sym(pi); exp(sym(1)) 5 6];
%! x = double (y);
%! A = zeta (x);
%! B = double (zeta (y));
%! assert (A, B, -4*eps);

%!test
%! % maple:
%! % > A := [1+2*I, -2 + 5*I, 100, 10*I, -1e-4 + 1e-6*I, -20 + I];
%! % > for a in A do evalf(Zeta(a)) end do;
%! x = [1+2i; -2+5i; 100; 10i; -1e-4 + 1e-6*1i; -20-1i];
%! A = [  0.59816556976238173670 - 0.35185474521784529050*1i
%!        0.21425967567391921717 + 0.52503846985036050707*1i
%!        1.0
%!        1.7564685929749629608 - 0.10151198543617116894*1i
%!       -0.49990811617645824900 - 0.91873792757763831501e-6*1i
%!        175.09070083717643866 - 71.512541417467273425*1i ];
%! B = zeta (x);
%! assert (A, B, -eps)

%!assert (zeta (inf), 1.0)

repo: octsympy


text:
function [r,ElapsedTime] = RunCaisFile(FileName,AnsysPath)
%FileName: File of .cai type, only File Name, No Extension!!!!!
%%
%if nargin<2
%    AnsysPath='"E:\CloudCache\3E46DDDB3C282A64C62650819BF418D0A3B645F6\0\##PROGRAM_FILES##\ANSYS Inc\v110\ANSYS\bin\intel\ansys110"';
%end
%FileName='test';
%% .cai to .inp
r1=Cai2Ansys(FileName);
%% run .inp
[r2 ElapsedTime]=RunAnsys(AnsysPath,FileName);
%% Return Value
r=r1+r2;

author: Sable

repo: mcbench-benchmarks


text:
% KM_DEMO_AKCCA Kernel-Based Identification of Hammerstein Systems. This
% demo generates some Hammerstein system input and output data and performs
% identification of the system using KIHAM.
%
% Author: Steven Van Vaerenbergh (steven *at* gtas.dicom.unican.es), 2013.
%
% The algorithm in this file is based on the following publication:
% S. Van Vaerenbergh and L. A. Azpicueta-Ruiz, "Kernel-Based Identification
% of Hammerstein Systems for Nonlinear Acoustic Echo-Cancellation", 2014
% IEEE International Conference on Acoustics, Speech, and Signal Processing
% (ICASSP), Florence, Italy, May 2014.
%
% This file is part of the Kernel Methods Toolbox for MATLAB.
% https://github.com/steven2358/kmbox

%% STATE

close all; clear
% rs = sum(100*clock);
rs = 1; % seed for random generator
rng('default')
rng(rs)

%% PARAMETERS

% system
f = @(x) tanh(x);
L = 256; % linear channel length
N = 1024; % number of samples
SNR = 50;

% kiham
kerneltype = 'gauss';
kernelpar = .5;
M = 25; % dictionary size
ca = .001;
ch = 1;
it_max = 50; % maximum number of iterations
it_stop = 1E-6; % stop if change in cost is smaller than this

%% GENERATE DATA

H = [1 randn(1,L-1)]';

x = 5*(2*rand(N,1)-1);
% x = randn(N,1);
z = f(x);
z_mem = km_memfill(z,L);
y_ref = z_mem*H;

sig = y_ref;
sigvar = var(sig);
noisevar = 10^(-SNR/10)*sigvar;
noise = sqrt(noisevar)*randn(size(sig));
y = y_ref + noise;

%% PROGRAM
tic

xdict = linspace(min(x),max(x),M)';

[alpha,h] = km_kiham(x,y,xdict,L,kerneltype,kernelpar,ca,ch,it_max,it_stop);

fprintf('Training: elapsed time is %.2f seconds.\n',toc)
%% OUTPUT

x_test = linspace(min(x),max(x),100)';
z_test_true = f(x_test);
K_test = km_kernel(x_test,xdict,kerneltype,kernelpar);
z_test_est = K_test*alpha;
scale = norm(h)/norm(H)*sign(h(1))*sign(H(1)); % internal ambiguity
z_test_est = z_test_est*scale;

figure; hold all
plot(x_test,z_test_true,'.');
% grid on
% legend('true nonlinearity')
plot(x_test,z_test_est,'.','Color',[0 .5 0]);
legend({'true nonlinearity','estimated nonlinearity'},'Location','SE')
grid on

figure; hold all
stem(H);
stem(h/scale)
legend('true linear channel','estimated linear channel')

fprintf('Test MSE = %.4f dB\n',10*log10(mean((z_test_true-z_test_est).^2)))

author: epfl-lasa


text:
function ps_shapefile(filename,data,step)
%
%   12/2013 MA: Initial version
%
%

if nargin < 3
  step=':';
end

% get velocity
%if ~exist('ps_plot_v-d.mat')
% ps_plot('v-d',-1)
%end
%load ps_plot_v-d  ph_disp
%data= ph_disp;

disp('Loading lonlat matrix...')
load ps2 lonlat

step
if isinteger(step)
 step=1:step:size(lonlat,1);
end

lon=lonlat(step,:,1);
lat=lonlat(step,:,2);
data=double(data(step,:));

shape=struct('Geometry', 'Point','X',num2cell(lon),'Y',num2cell(lat),'V',num2cell(V));

shapewrite(shape,filename);

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
%  Figure 3.32      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
% script to generate Fig. 3.32
% Pole-zero effects
%
close all;
clear all;
Den1=[1 0.2 1.01];
Den2=[1 1.0];
Den11=conv(Den1,Den2);
alpha=0.1;
beta=1.0;
Num1=(1.01*1.0/(alpha^2+beta^2))*([1 2*alpha alpha^2+beta^2]);

t=0:.1:20;
sys=tf(Num1,Den11);
[y1]=step(sys,t)

Den1=[1 0.2 1.01];
Den2=[1 1.0];
Den22=conv(Den1,Den2);
% alpha=1;
alpha1=0.25;
beta=1.0;
Num2=(1.01*1.0/(alpha1^2+beta^2))*([1 2*alpha1 alpha1^2+beta^2]);

t=0:.1:20;
sys2=tf(Num2,Den22);
[y2]=step(sys2,t)

Den1=[1 0.2 1.01];
Den2=[1 1.0];
Den33=conv(Den1,Den2);
alpha2=0.5;
beta=1.0;
Num3=(1.01*1.0/(alpha2^2+beta^2))*([1 2*alpha2 alpha2^2+beta^2]);

t=0:.1:20;
sys2=tf(Num3,Den33);
[y3]=step(sys2,t)
plot(t,y1,t,y2,t,y3);
xlabel('Time (sec)');
ylabel('Unit step response');
text(2, 1.36, '\alpha = 0.5');
text(2,1.14,'\alpha = 0.25');
text(2,0.8,'\alpha = 0.1');
title('Fig. 3.32: Step responses');
%%%%%%%%%%%%%%%%%%%%%

%grid
nicegrid


text:
function Z = spinv(A, B)
% SPINV    Evaluate the sparse inverse matrix
%
% z = sinv(A)  returns the elements of inv(A)_ij, for which A_ij
%      is different from zero.
%
% z = sinv(LD, 1)  returns the elements of inv(A)_ij, for which A_ij
%     is different from zero, and where LD is the LDL cholesky
%     decomposition of A. LD has to be in the form returned by ldlchol
%     in SuiteSparse by Tim Davis.
%
%   Note! If z = sinv(LD, 1) is used LD must not be modified in Matlab
%   after ldlchol. Matlab destroys the symbolic sparsity structure in
%   the Cholesky decomposition, which is needed in the spinv
%   algorithm. If LD is modified the worst scenario is memory corruption.
%
%
%     For details, see:
%     Jarno Vanhatalo and Aki Vehtari (2008). Modelling local and
%     global phenomena with sparse Gaussian processes. Proceedings of
%     the 24th Conference on Uncertainty in Artificial Intelligence
%

% Copyright (c) 2008-2010      Jarno Vanhatalo

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

warning('Using M-file. Compile MEX-file!');

    n = size(A,1);

    if nargin == 1
        [LD, p, q] = ldlchol(A);
    else
        LD = A;
    end

    [I,J,ld] = find(LD);
    temp = [I(:) J(:) ; J(:) I(:)];
    temp = sortrows(unique(temp,'rows'),2);
    Iz = temp(:,1); Jz = temp(:,2);

    % Find the column starting points
    a1=zeros(n,1);
    a2 = cumsum(histc(J,1:n));
    a1(1) = 1; a1(2:end) = a2(1:end-1) + 1;
    az1=zeros(n,1);
    az2 = cumsum(histc(Jz,1:n));
    az1(1) = 1; az1(2:end) = az2(1:end-1) + 1;

    for j=1:n
        indaz{j} = az1(j):az2(j);
        indIz{j} = Iz(indaz{j})';
    end

    % Evaluate the sparse inverse
    z = zeros(size(Iz));
    z(end) = 1./ld(end);
    % Allocate memory
    cindit=zeros(n,1);
    for jj = n-1:-1:1
        fil = ld(a1(jj)+1:a1(jj+1)-1);
        fi = I(a1(jj)+1:a1(jj+1)-1);
        lfi = length(fi);
        Zt = zeros(lfi,lfi);
        indz = cumsum(histc(indIz{jj},[0 ; fi]));
        indz = az1(jj) + indz(1:end-1);

        i4=0;
        for i1 = 1:lfi
            cind1=indaz{fi(i1)};
            Icind1=indIz{fi(i1)};
            indfi = lfi;
            i2=length(Icind1);
            go = true;
            while go
                if Icind1(i2)==jj  % Find the indeces for the jj'th rows in fi columns
                    i4=i4+1;
                    cindit(i4)=cind1(i2);
                    go = false;
                end
                if indfi >= 1 && fi(indfi) == Icind1(i2) % Find the indeces for the fi'th rows in i2'nd columns
                    Zt(indfi,i1) = z(cind1(i2));
                    indfi = indfi-1;
                end
                i2 = i2-1;
            end
        end
        % remove extras
        cindi=cindit(1:i4);

        zij = -fil'*Zt;
        z(cindi) = zij;
        z(indz) = zij;
        zij = 1./ld(a1(jj)) - fil'*z(indz);
        z(az1(jj)-1+find(indIz{jj}==jj,1)) = zij;
    end

    Z = sparse(Iz,Jz,z);

    if nargin == 1
        r(q) = 1:n;
        Z = Z(r,r);
    end
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [p] = loadPlant_QuadRotor3d()
%
% Convenience function to load 3d quadcopter.
% Provided to simplify main function and modularize plant model definition.
%
% Depends:
%   - definePropulsionModel.m

% Enviromental params
p.g = -9.81 ; % World Coords is XYZ = [East, North, Up], i.e. gravity is a negative number
p.rho = 1.225 ; % air density during flight (kg/m^3)

% Inertial params
p.m = 5 ;
p.I = [0.625 0 0; 0 0.625 0; 0 0 1.25] ; % inertia tensor coords:
p.cg = [0 0 0] ; % (m) location of center of gravity

% control params
p.uMax = 1 ; % maximum throttle setting

% Propulsion system params - shared for all motors:
qRP.d_prop = 0.305*ones(4,1) ; % propeller diameter (m)
qRP.maxThrust = 25*ones(4,1) ; % thrust at 100% throttle (N)
qRP.maxRPM = 10000*ones(4,1) ; % RPM at 100% throttle (RPM)
qRP.maxTorque = ones(4,1) ;  % torque at 100% throttle (Nm)
qRP.thrustLocations = [0.5 0 0; 0 0.5 0; -0.5 0 0; 0 -0.5 0]; % motor locations (each row one motor in coords: [port, nose, top]
qRP.thrustAxes = repmat([0 0 1],4,1) ; % thrust axes of each motor in coords port, nose, top.
qRP.isSpinDirectionCCW = [1; 0; 1; 0] ; % bool to reverse motor spin direction around 'thrustAxes'.
plotflag = 0 ;
[p.propulsion] = definePropulsionModel(qRP,plotflag);

author: MatthewPeterKelly

repo: OptimTraj