text:
function display(c)

%tstoolbox/@core/display
%   Syntax:
%     * display(c)
%
%   Input Arguments:
%     * c - core object
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

disp(['  Dlens : ' num2str(c.dlens)])

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% fixedSegmentation = removeBridgesMex(segmentation, anatomicalData,
%			 averageWhiteIntensity, thresholdIntensity, averageGrayIntensity)
%
% Mex function implementing the bridge removal algorithm described
% in Kriegeskorte & Goebel (2001) An Efficient Algorithm for
% Topologically Correct Segmentation of the Cortical Sheet in Anatomical
% MR Volumes. Neuroimage, 14, 329-346.
%
% To compile, try:
% mex removeBridgesMex.cpp CRemoveBridges.cpp
%
% On linux, this might make faster code:
% mex -O COPTIMFLAGS='-O3 -march=i686 -DNDEBUG' removeBridgesMex.cpp CRemoveBridges.cpp
%
% HISTORY (from the source code):
% 2002.08.?? Niko & Rainer kindly gave us their code from BrainVoyager.
% 2002.09.?? Ian Spiro modified the BV code to build a DLL that takes
%            data rather than filenames. He also wrote a simple mex wrapper
%            to call this DLL.
% 2002.11.13 Bob Dougherty: overhauled the code to make it cleaner (no more
%            globals!) and platform-independent. To do this, I had to remove
%            all the GUI stuff and the threading. The code now runs serially
%            and uses printf for progress reporting (although this can be
%            easily changed by modifying the updateProgress method). I also
%            wrote a mex wrapper for this class.

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


%% CAVE-eitan_office
% %{
imfolder = 'images\CAVE-eitan_office';
im_n = 6;
imfile = cell(im_n,1);
for ii = 1:im_n
    imfile{ii} = sprintf('%s\\e-%02d.jpg', imfolder, ii);
end

im = cell(im_n,1);
for ii = 1:im_n
    im{ii} = imread(imfile{ii});
end

imsize = zeros(im_n,3);

for ii = 1:im_n
    imsize(ii,:) = size(im{ii});
    if imsize(ii,1) > 720
        scale = 720/size(im{ii}, 1);
        im{ii} = imresize(im{ii}, scale);

        imsize(ii,:) = size(im{ii});
    end
end

mosaic = REW_mosaic( im, [], 0, 'persp', 0.01, imfolder );
%}


%% BVAR tutorial: Inference with Minnesota Prior
% Author:   Filippo Ferroni
% Date:     27/02/2020

clear all
close all
clc
addpath ../../cmintools/
addpath ../../bvartools/
pkg load statistics
pkg load io

%% %=========================================================================
%%% INFERENCE %%%
%%=========================================================================

% load the data
load('../BVAR tutorial/Data.mat')
y= [IPI HICP CORE Euribor1Y M3 EXRATE];% collect the variables used in the VAR

%% Ex1/ Minnesota Prior

lags = 6;
options.max_minn_hyper  = 1;
options.minn_prior_tau  = 10;      % set tau
options.index_est       = [3 4];   % define the hyper-parameters over which to maximize
options.lb              = [0 0];   % sets the lower bounds
options.ub              = [20 20]; % sets the upper bounds
options.max_compute     = 3;       % optimization  by Matlab Simplex
BVAR                    = bvar(y,lags,options);

%% Ex2/ Minnesota Prior

clear options
lags = 6;
% setting the default values for the hyperparameters
hyperpara(1)    = 3;		  % tau
hyperpara(2)    = 0.5;		  % decay
hyperpara(3)    = 5;		  % lambda
hyperpara(4)    = 2;		  % mu
hyperpara(5)    = 2;		  % omega
% setting the options
options.index_est	   = 1:1;      % hyper-parameter over which maximize
options.max_compute    = 3;      % maximize  using Matlab fmincon function
options.lb             = [0.05]; % Lower bound
[postmode,logmlike,HH] = bvar_max_hyper(hyperpara,y,lags,options);

%% Ex3

hyperpara(1)            = postmode(1); % use as starting value previous mode
options.index_est       = 1:3; % set hyper-parameters over which maximize
options.lb              = [0.05 0.05 0.05]; % Lower bounds
options.ub              = [50 50 50];       % Upper bounds
[postmode,log_dnsty,HH] = bvar_max_hyper(hyperpara,y,lags,options);


text:
classdef BatchNorm < dagnn.ElementWise
  properties
    numChannels
    epsilon = 1e-5
    opts = {'NoCuDNN'} % ours seems slightly faster
  end

  properties (Transient)
    moments
  end

  methods
    function outputs = forward(obj, inputs, params)
      if strcmp(obj.net.mode, 'test')
        outputs{1} = vl_nnbnorm(inputs{1}, params{1}, params{2}, ...
                                'moments', params{3}, ...
                                'epsilon', obj.epsilon, ...
                                obj.opts{:}) ;
      else
        [outputs{1},obj.moments] = ...
            vl_nnbnorm(inputs{1}, params{1}, params{2}, ...
                       'epsilon', obj.epsilon, ...
                       obj.opts{:}) ;
      end
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      [derInputs{1}, derParams{1}, derParams{2}, derParams{3}] = ...
        vl_nnbnorm(inputs{1}, params{1}, params{2}, derOutputs{1}, ...
                   'epsilon', obj.epsilon, ...
                   'moments', obj.moments, ...
                   obj.opts{:}) ;
      obj.moments = [] ;
      % multiply the moments update by the number of images in the batch
      % this is required to make the update additive for subbatches
      % and will eventually be normalized away
      derParams{3} = derParams{3} * size(inputs{1},4) ;
    end

    % ---------------------------------------------------------------------
    function obj = BatchNorm(varargin)
      obj.load(varargin{:}) ;
    end

    function params = initParams(obj)
      params{1} = ones(obj.numChannels,1,'single') ;
      params{2} = zeros(obj.numChannels,1,'single') ;
      params{3} = zeros(obj.numChannels,2,'single') ;
    end

    function attach(obj, net, index)
      attach@dagnn.ElementWise(obj, net, index) ;
      p = net.getParamIndex(net.layers(index).params{3}) ;
      net.params(p).trainMethod = 'average' ;
      net.params(p).learningRate = 0.1 ;
    end
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [ov ov_n1 ov_n2] = calc_overlap(dres1, f1, dres2, f2)
%%f2 can be an array and f1 should be a scalar.
%%% this will find the overlap between dres1(f1) (only one) and all detection windows in dres2(f2(:))

f2=f2(:)';
n = length(f2);

cx1 = dres1.x(f1);
cx2 = dres1.x(f1)+dres1.w(f1)-1;
cy1 = dres1.y(f1);
cy2 = dres1.y(f1)+dres1.h(f1)-1;

gx1 = dres2.x(f2);
gx2 = dres2.x(f2)+dres2.w(f2)-1;
gy1 = dres2.y(f2);
gy2 = dres2.y(f2)+dres2.h(f2)-1;

ca = dres1.h(f1).*dres1.w(f1);  %% area
ga = dres2.h(f2).*dres2.w(f2);

%%% find the overlapping area
xx1 = max(cx1, gx1);
yy1 = max(cy1, gy1);
xx2 = min(cx2, gx2);
yy2 = min(cy2, gy2);
w = xx2-xx1+1;
h = yy2-yy1+1;

inds  = find((w > 0) .* (h > 0));  %% real overlap
ov    = zeros(1, n);
ov_n1 = zeros(1, n);
ov_n2 = zeros(1, n);
inter = w(inds).*h(inds);         %% area of overlap
u     = ca + ga(inds) - w(inds).*h(inds);     %% area of union
ov(inds)    = inter ./ u;                       %% intersection / union
ov_n1(inds) = inter / ca;                       %% intersection / area in dres1
ov_n2(inds) = inter ./ ga(inds);                  %% intersection / area in dres2

author: yuxng


text:
function obj = subsasgn(obj,subs,dat)
% Overloaded subsasgn function for file_array objects
%__________________________________________________________________________
% Copyright (C) 2005-2013 Wellcome Trust Centre for Neuroimaging

%
% $Id: subsasgn.m 6157 2014-09-05 18:17:54Z guillaume $

if isempty(subs), return; end

%-Subscript types '.' or '{}'
%==========================================================================
if ~strcmp(subs(1).type,'()')
    if strcmp(subs(1).type,'.')
        %error('Attempt to reference field of non-structure array.');
        if numel(struct(obj))~=1
            error('Can only change the fields of simple file_array objects.');
        end
        switch(subs(1).subs)
        case 'fname',      obj = asgn(obj,@fname,     subs(2:end),dat); %fname(obj,dat);
        case 'dtype',      obj = asgn(obj,@dtype,     subs(2:end),dat); %dtype(obj,dat);
        case 'offset',     obj = asgn(obj,@offset,    subs(2:end),dat); %offset(obj,dat);
        case 'dim',        obj = asgn(obj,@dim,       subs(2:end),dat); %obj = dim(obj,dat);
        case 'scl_slope',  obj = asgn(obj,@scl_slope, subs(2:end),dat); %scl_slope(obj,dat);
        case 'scl_inter',  obj = asgn(obj,@scl_inter, subs(2:end),dat); %scl_inter(obj,dat);
        case 'permission', obj = asgn(obj,@permission,subs(2:end),dat); %permission(obj,dat);
        otherwise, error(['Reference to non-existent field "' subs.subs '".']);
        end
        return;
    end
    if strcmp(subs(1).type,'{}')
        error('Cell contents reference from a non-cell array object.');
    end
end

%-Subscript type '()'
%==========================================================================
if numel(subs)~=1, error('Expression too complicated.'); end

dm   = size(obj);
sobj = struct(obj);

if length(subs.subs) < length(dm)
    l   = length(subs.subs);
    dm  = [dm(1:(l-1)) prod(dm(l:end))];
    if numel(sobj) ~= 1
        error('Can only reshape simple file_array objects.');
    end
    if numel(sobj.scl_slope)>1 || numel(sobj.scl_inter)>1
        error('Can not reshape file_array objects with multiple slopes and intercepts.');
    end
end

%-Index vectors
%--------------------------------------------------------------------------
dm   = [dm ones(1,16)];
di   = ones(1,16);
args = cell(1,length(subs.subs));
for i=1:length(subs.subs)
    if ischar(subs.subs{i})
        if ~strcmp(subs.subs{i},':')
            error('Invalid subscript value.');
        end
        args{i} = int32(1:dm(i));
    else
        args{i} = int32(subs.subs{i});
    end
    di(i) = length(args{i});
end

%-Check read/write permissions
%--------------------------------------------------------------------------
for j=1:length(sobj)
    if strcmp(sobj(j).permission,'ro')
        error('Array is read-only.');
    end
end

%-Perform the assignment
%--------------------------------------------------------------------------
if length(sobj)==1
    sobj.dim = dm;
    if numel(dat)~=1
        subfun(sobj,double(dat),args{:});
    else
        dat = repmat(double(dat),di);
        subfun(sobj,dat,args{:});
    end
else
    for j=1:length(sobj)
        ps  = [sobj(j).pos ones(1,length(args))];
        dm  = [sobj(j).dim ones(1,length(args))];
        siz = ones(1,16);
        for i=1:length(args)
            msk      = args{i}>=ps(i) & args{i}<(ps(i)+dm(i));
            args2{i} = find(msk);
            args3{i} = int32(double(args{i}(msk))-ps(i)+1);
            siz(i)   = numel(args2{i});
        end
        if numel(dat)~=1
            dat1 = double(subsref(dat,struct('type','()','subs',{args2})));
        else
            dat1 = double(dat) + zeros(siz);
        end
        subfun(sobj(j),dat1,args3{:});
    end
end

%==========================================================================
% function sobj = subfun(sobj,dat,varargin)
%==========================================================================
function sobj = subfun(sobj,dat,varargin)
va = varargin;

%-Get datatype
%--------------------------------------------------------------------------
dt  = datatypes;
ind = find(cat(1,dt.code)==sobj.dtype);
if isempty(ind), error('Unknown datatype'); end
if dt(ind).isint, dat(~isfinite(dat)) = 0; end

%-Apply DC offset
%--------------------------------------------------------------------------
if ~isempty(sobj.scl_inter)
    inter = sobj.scl_inter;
    if numel(inter) > 1
        inter = resize_scales(inter,sobj.dim,varargin);
    end
    dat = double(dat) - inter;
end

%-Apply scalefactor
%--------------------------------------------------------------------------
if ~isempty(sobj.scl_slope)
    slope = sobj.scl_slope;
    if numel(slope) > 1
        slope = resize_scales(slope,sobj.dim,varargin);
        dat   = double(dat) ./ slope;
    else
        dat   = double(dat) / slope;
    end
end

%-Convert data into output datatype
%--------------------------------------------------------------------------
if dt(ind).isint, dat = round(dat); end

%ws = warning('off'); % Avoid warning messages in R14 SP3
dat = feval(dt(ind).conv,dat);
%warning(ws);

%-Write data to file
%--------------------------------------------------------------------------
nelem = dt(ind).nelem;
if nelem==1
    mat2file(sobj,dat,va{:});

    %if sobj.be, swap = @(x) swapbytes(x); else swap = @(x) x; end
    %m = memmapfile(sobj.fname, 'Format', {dt(ind).prec, sobj.dim,'dat'}, ...
    %       'Offset', sobj.offset, 'Writable', true);
    %m.Data.dat = subsasgn(m.Data.dat, substruct('()',va), swap(dat));
elseif nelem==2
    sobj1       = sobj;
    sobj1.dim   = [2 sobj.dim];
    sobj1.dtype = dt(find(strcmp(dt(ind).prec,{dt.prec}) & (cat(2,dt.nelem)==1))).code;
    dat         = reshape(dat,[1 size(dat)]);
    dat         = [real(dat) ; imag(dat)];
    mat2file(sobj1,dat,int32([1 2]),va{:});
else
    error('Inappropriate number of elements per voxel.');
end

%==========================================================================
% function obj = asgn(obj,fun,subs,dat)
%==========================================================================
function obj = asgn(obj,fun,subs,dat)
if ~isempty(subs)
    tmp = feval(fun,obj);
    tmp = subsasgn(tmp,subs,dat);
    obj = feval(fun,obj,tmp);
else
    obj = feval(fun,obj,dat);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Test file for trigtech/imag.m

function pass = test_imag(pref)

% Get preferences.
if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();

%%
% Test a scalar-valued function.
f = testclass.make(@(x) cos(pi*x) + 1i*sin(pi*x), [], pref);
g = testclass.make(@(x) sin(pi*x), [], pref);
h = imag(f);
g = prolong(g, length(h));
pass(1) = norm(h.coeffs - g.coeffs, inf) < 10*vscale(h)*eps;

%%
% Test an array-valued function.
f = testclass.make(@(x) [cos(sin(pi*x)) + 1i*sin(cos(pi*x)), -exp(1i*pi*x)], [], pref);
g = testclass.make(@(x) [sin(cos(pi*x)), -imag(exp(1i*pi*x))], [], pref);
h = imag(f);
g = prolong(g, length(h));
pass(2) = norm(h.coeffs - g.coeffs, inf) < 100*max(vscale(h)*eps);

%%
% Test a real function.
f = testclass.make(@(x) cos(pi*x), [], pref);
g = imag(f);
pass(3) = numel(g.coeffs) == 1;

%%
% Test an array-valued real function.
f = testclass.make(@(x) [cos(pi*x), sin(pi*x)], [], pref);
g = imag(f);
pass(4) = all(size(g.coeffs) == [1, 2]) && all(g.coeffs == 0);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [Y, optinf] = celnet_gpu(D, S, lambda, mu, opt)

% celnet_gpu -- Convolutional Elastic Net (GPU version)
%
%         argmin_{x_m} (1/2)||\sum_m d_m * x_m - s||_2^2 +
%                      lambda \sum_m ||x_m||_1 + (mu/2) \sum_m ||x_m||_2^2
%
%         The solution is computed using an ADMM approach (see
%         boyd-2010-distributed) with efficient solution of the main
%         linear systems (see wohlberg-2016-efficient).
%
% Usage:
%       [Y, optinf] = celnet_gpu(D, S, lambda, mu, opt)
%
% Input:
%       D           Dictionary filter set (3D array)
%       s           Input image
%       lambda      Regularization parameter (l1)
%       mu          Regularization parameter (l2)
%       opt         Algorithm parameters structure
%
% Output:
%       Y           Dictionary coefficient map set (3D array)
%       optinf      Details of optimisation
%
%
% Options structure fields:
%   Verbose          Flag determining whether iteration status is displayed.
%                    Fields are iteration number, functional value,
%                    data fidelity term, l1 regularisation term, l2
%                    regularisation term, and primal and dual residuals
%                    (see Sec. 3.3 of boyd-2010-distributed). The value of
%                    rho is also displayed if options request that it is
%                    automatically adjusted.
%   MaxMainIter      Maximum main iterations
%   AbsStopTol       Absolute convergence tolerance (see Sec. 3.3.1 of
%                    boyd-2010-distributed)
%   RelStopTol       Relative convergence tolerance (see Sec. 3.3.1 of
%                    boyd-2010-distributed)
%   L1Weight         Weighting array for coefficients in l1 norm of X.
%                    Array should have the same dimensions as X, but the
%                    first two dimensions may be of unit size, corresponding
%                    to a weighting that varies with filter index but is
%                    spatially constant.
%   L2Weight         Weighting array for l2 norm of X. Array should have
%                    dimensions corresponding to the non-spatial dimensions
%                    of X since spatial weighting is no possible (i.e.
%                    weighting varies only with filter and sample index).
%   Y0               Initial value for Y
%   U0               Initial value for U
%   rho              ADMM penalty parameter
%   AutoRho          Flag determining whether rho is automatically updated
%                    (see Sec. 3.4.1 of boyd-2010-distributed)
%   AutoRhoPeriod    Iteration period on which rho is updated
%   RhoRsdlRatio     Primal/dual residual ratio in rho update test
%   RhoScaling       Multiplier applied to rho when updated
%   AutoRhoScaling   Flag determining whether RhoScaling value is
%                    adaptively determined (see wohlberg-2015-adaptive). If
%                    enabled, RhoScaling specifies a maximum allowed
%                    multiplier instead of a fixed multiplier.
%   RhoRsdlTarget    Residual ratio targeted by auto rho update policy.
%   StdResiduals     Flag determining whether standard residual definitions
%                    (see Sec 3.3 of boyd-2010-distributed) are used instead
%                    of normalised residuals (see wohlberg-2015-adaptive)
%   RelaxParam       Relaxation parameter (see Sec. 3.4.3 of
%                    boyd-2010-distributed)
%   NoBndryCross     Flag indicating whether all solution coefficients
%                    corresponding to filters crossing the image boundary
%                    should be forced to zero.
%   AuxVarObj        Flag determining whether objective function is computed
%                    using the auxiliary (split) variable
%   HighMemSolve     Use more memory for a slightly faster solution
%
%
% Authors: Brendt Wohlberg <brendt@lanl.gov>
%          Ping-Keng Jao <jpk7656@gmail.com>
% Modified: 2015-12-18
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

gS = gpuArray(S);
gD = gpuArray(D);
glambda = gpuArray(lambda);

if nargin < 5,
  opt = [];
end
if nargin < 4,
  gmu = gpuArray(0);
else
  gmu = gpuArray(mu);
end
checkopt(opt, defaultopts([]));
opt = defaultopts(opt);

% Set up status display for verbose operation
hstr = 'Itn   Fnc       DFid      l1        l2        r         s      ';
sfms = '%4d %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e';
nsep = 64;
if opt.AutoRho,
  hstr = [hstr '   rho  '];
  sfms = [sfms ' %9.2e'];
  nsep = nsep + 10;
end
if opt.Verbose && opt.MaxMainIter > 0,
  disp(hstr);
  disp(char('-' * ones(1,nsep)));
end

% Start timer
tstart = tic;

% Collapsing of trailing singleton dimensions greatly complicates
% handling of both SMV and MMV cases. The simplest approach would be
% if S could always be reshaped to 4d, with dimensions consisting of
% image rows, image cols, a single dimensional placeholder for number
% of filters, and number of measurements, but in the single
% measurement case the third dimension is collapsed so that the array
% is only 3d.
if size(S,3) > 1,
  xsz = [size(S,1) size(S,2) size(D,3) size(S,3)];
  % Insert singleton 3rd dimension (for number of filters) so that
  % 4th dimension is number of images in input s volume
  S = reshape(S, [size(S,1) size(S,2) 1 size(S,3)]);
else
  xsz = [size(S,1) size(S,2) size(D,3) 1];
end
% Compute filters in DFT domain
gDf = fft2(gD, size(S,1), size(S,2));
% Convolve-sum and its Hermitian transpose
gDop = @(x) sum(bsxfun(@times, gDf, x), 3);
gDHop = @(x) bsxfun(@times, conj(gDf), x);
% Compute signal in DFT domain
gSf = fft2(gS);
% S convolved with all filters in DFT domain
gDSf = gDHop(gSf);

% Set up l2 weight array
if isscalar(opt.L2Weight),
  gwl2 = gpuArray(opt.L2Weight);
else
  gwl2 = gpuArray(reshape(opt.L2Weight, [1 1 size(opt.L2Weight,1) ...
                      size(opt.L2Weight,2)]));
end

% Default lambda is 1/10 times the lambda value beyond which the
% solution is a zero vector
if nargin < 3 | isempty(glambda),
  gb = ifft2(gDHop(gSf), 'symmetric');
  glambda = 0.1*max(vec(abs(gb)));
end
% Set up algorithm parameters and initialise variables
grho = gpuArray(opt.rho);
if isempty(grho), grho = 50*glambda+1; end;
gmwr = gmu*gwl2 + grho;
if isempty(opt.RhoRsdlTarget),
  if opt.StdResiduals,
    opt.RhoRsdlTarget = 1;
  else
    opt.RhoRsdlTarget = 1 + (18.3).^(log10(glambda) + 1);
  end
end
if opt.HighMemSolve,
  gcn = bsxfun(@rdivide, gDf, gmwr);
  gcd = sum(gDf.*bsxfun(@rdivide, conj(gDf), gmu*gwl2 + grho), 3) + 1.0;
  gC = bsxfun(@rdivide, gcn, gcd);
  clear cn cd;
else
  C = [];
end
gNx = prod(gpuArray(xsz));
optinf = struct('itstat', [], 'opt', opt);
gr = gpuArray(Inf);
gs = gpuArray(Inf);
gepri = gpuArray(0);
gedua = gpuArray(0);

% Initialise main working variables
% X = [];
if isempty(opt.Y0),
  gY = gpuArray.zeros(xsz);
else
  gY = gpuArray(opt.Y0);
end
gYprv = gY;
if isempty(opt.U0),
  if isempty(opt.Y0),
    gU = gpuArray.zeros(xsz);
  else
    gU = (glambda/grho)*sign(gY);
  end
else
  gU = gpuArray(opt.U0);
end

% Main loop
k = 1;
while k <= opt.MaxMainIter & (gr > gepri | gs > gedua),

  % Solve X subproblem
  gXf = solvedbd_sm(gDf, gmwr, gDSf + grho*fft2(gY - gU), gC);
  gX = ifft2(gXf, 'symmetric');

  % See pg. 21 of boyd-2010-distributed
  if opt.RelaxParam == 1,
    gXr = gX;
  else
    gXr = opt.RelaxParam*gX + (1-opt.RelaxParam)*gY;
  end

  % Solve Y subproblem
  gY = shrink(gXr + gU, (glambda/grho)*opt.L1Weight);
  if opt.NoBndryCross,
    gY((end-size(gD,1)+2):end,:,:,:) = 0;
    gY(:,(end-size(gD,1)+2):end,:,:) = 0;
  end

  % Update dual variable
  gU = gU + gXr - gY;

  % Compute data fidelity term in Fourier domain (note normalisation)
  if opt.AuxVarObj,
    gYf = fft2(gY); % This represents unnecessary computational cost
    gJdf = sum(vec(abs(sum(bsxfun(@times,gDf,gYf),3)-gSf).^2)) / ...
           (2*xsz(1)*xsz(2));
    gJl1 = sum(abs(vec(bsxfun(@times, opt.L1Weight, gY))));
    gJl2 = sum(vec(gwl2.*sum(sum(gY.^2, 1),2)))/2;
  else
    gJdf = sum(vec(abs(sum(bsxfun(@times,gDf,gXf),3)-gSf).^2)) / ...
           (2*xsz(1)*xsz(2));
    gJl1 = sum(abs(vec(bsxfun(@times, opt.L1Weight, gX))));
    gJl2 = sum(vec(gwl2.*sum(sum(gX.^2, 1),2)))/2;
  end
  gJfn = gJdf + glambda*gJl1 + gmu*gJl2;

  gnX = norm(gX(:)); gnY = norm(gY(:)); gnU = norm(gU(:));
  if opt.StdResiduals,
    % See pp. 19-20 of boyd-2010-distributed
    gr = norm(vec(gX - gY));
    gs = norm(vec(grho*(gYprv - gY)));
    gepri = sqrt(gNx)*opt.AbsStopTol+max(gnX,gnY)*opt.RelStopTol;
    gedua = sqrt(gNx)*opt.AbsStopTol+grho*gnU*opt.RelStopTol;
  else
    % See wohlberg-2015-adaptive
    gr = norm(vec(gX - gY))/max(gnX,gnY);
    gs = norm(vec(gYprv - gY))/gnU;
    gepri = sqrt(gNx)*opt.AbsStopTol/max(gnX,gnY)+opt.RelStopTol;
    gedua = sqrt(gNx)*opt.AbsStopTol/(grho*gnU)+opt.RelStopTol;
  end

  % Record and display iteration details
  tk = toc(tstart);
  optinf.itstat = [optinf.itstat; [k gather(gJfn) gather(gJdf) gather(gJl1) ...
                   gather(gJl2) gather(gr) gather(gs) gather(gepri) ...
                      gather(gedua) gather(grho) tk]];
  if opt.Verbose,
    if opt.AutoRho,
      disp(sprintf(sfms, k, gather(gJfn), gather(gJdf), gather(gJl1), ...
                   gather(gJl2), gather(gr), gather(gs), gather(grho)));
    else
      disp(sprintf(sfms, k, gather(gJfn), gather(gJdf), gather(gJl1), ...
                   gather(gJl2), gather(gr), gather(gs)));
    end
  end

  % See wohlberg-2015-adaptive and pp. 20-21 of boyd-2010-distributed
  if opt.AutoRho,
    if k ~= 1 && mod(k, opt.AutoRhoPeriod) == 0,
      if opt.AutoRhoScaling,
        grhomlt = sqrt(gr/(gs*opt.RhoRsdlTarget));
        if grhomlt < 1, grhomlt = 1/grhomlt; end
        if grhomlt > opt.RhoScaling, grhomlt = opt.RhoScaling; end
      else
        rhomlt = opt.RhoScaling;
      end
      grsf = 1;
      if gr > opt.RhoRsdlTarget*opt.RhoRsdlRatio*gs, grsf = grhomlt; end
      if gs > (opt.RhoRsdlRatio/opt.RhoRsdlTarget)*gr, grsf = 1/grhomlt; end
      grho = grsf*grho;
      gU = gU/grsf;
      if opt.HighMemSolve && grsf ~= 1,
        gmwr = gmu*gwl2 + grho;
        gcn = bsxfun(@rdivide, gDf, gmwr);
        gcd = sum(gDf.*bsxfun(@rdivide, conj(gDf), gmwr), 3) + 1.0;
        gC = bsxfun(@rdivide, gcn, gcd);
        clear gcn gcd;
      end
    end
  end

  gYprv = gY;
  k = k + 1;

end

% Record run time and working variables
optinf.runtime = toc(tstart);
optinf.X = gather(gX);
optinf.Xf = gather(gXf);
optinf.Y = gather(gY);
optinf.U = gather(gU);
optinf.lambda = gather(glambda);
optinf.mu = gather(mu);
optinf.rho = gather(grho);
Y = gather(gY);
if opt.Verbose && opt.MaxMainIter > 0,
  disp(char('-' * ones(1,nsep)));
end

return

function u = vec(v)

  u = v(:);

return

function u = shrink(v, lambda)

  if isscalar(lambda),
    u = sign(v).*max(0, abs(v) - lambda);
  else
    u = sign(v).*max(0, bsxfun(@minus, abs(v), lambda));
  end

return

function opt = defaultopts(opt)

  if ~isfield(opt,'Verbose'),
    opt.Verbose = 0;
  end
  if ~isfield(opt,'MaxMainIter'),
    opt.MaxMainIter = 1000;
  end
  if ~isfield(opt,'AbsStopTol'),
    opt.AbsStopTol = 0;
  end
  if ~isfield(opt,'RelStopTol'),
    opt.RelStopTol = 1e-4;
  end
  if ~isfield(opt,'L1Weight'),
    opt.L1Weight = 1;
  end
  if ~isfield(opt,'L2Weight'),
    opt.L2Weight = 1;
  end
  if ~isfield(opt,'Y0'),
    opt.Y0 = [];
  end
  if ~isfield(opt,'U0'),
    opt.U0 = [];
  end
  if ~isfield(opt,'rho'),
    opt.rho = [];
  end
  if ~isfield(opt,'AutoRho'),
    opt.AutoRho = 1;
  end
  if ~isfield(opt,'AutoRhoPeriod'),
    opt.AutoRhoPeriod = 1;
  end
  if ~isfield(opt,'RhoRsdlRatio'),
    opt.RhoRsdlRatio = 1.2;
  end
  if ~isfield(opt,'RhoScaling'),
    opt.RhoScaling = 100;
  end
  if ~isfield(opt,'AutoRhoScaling'),
    opt.AutoRhoScaling = 1;
  end
  if ~isfield(opt,'RhoRsdlTarget'),
    opt.RhoRsdlTarget = [];
  end
  if ~isfield(opt,'StdResiduals'),
    opt.StdResiduals = 0;
  end
  if ~isfield(opt,'RelaxParam'),
    opt.RelaxParam = 1.8;
  end
  if ~isfield(opt,'NoBndryCross'),
    opt.NoBndryCross = 0;
  end
  if ~isfield(opt,'AuxVarObj'),
    opt.AuxVarObj = 0;
  end
  if ~isfield(opt,'HighMemSolve'),
    opt.HighMemSolve = 0;
  end

return

author: thfylsty


text:
function D = sldiff_pw(X1, X2, type)
%SLDIFF_PW Measures the pair-wise difference
%
% $ Synatx $
%   - D = sldiff_pw(X1, X2, type)
%
% $ Arguments $
%   - X1, X2:             the two sample matrices
%   - type:               the type of difference measurement
%                         default = 'abssum'
%
% $ Description $
%   - D = sldiff_pw(X1, X2, type) computes the measurment of differences
%     between the samples in X1 and those in X2 in a pairwise manner.
%     All samples should be stored in columns. And the samples in X1 and
%     X2 should be of the same dimension. If X1 and X2 are of sizes
%     dxn1 and dxn2 respectively. Then D is a matrix of size n1 x n2.
%
%   - The measurment types supported are listed below
%     \*
%     \t  Table 1. The difference measurement types                 \\
%     \h    name      &         description                         \\
%          'abssum'   &   sum of absolute values of differences     \\
%          'maxdiff'  &   maximum of absolute values of differences \\
%          'mindiff'  &   minimum of absolute values of differences \\
%     \*
%
% $ History $
%   - Created by Dahua Lin on Dec 06th, 2005
%   - Modified by Dahua Lin on Sep 10th, 2006
%       - Re-implement the core in C++: pwdiff_core
%       - The efficiency is increased by 10 times.
%

%% parse and verify input arguments

if nargin < 3 || isempty(type)
    type = 'abssum';
end

switch type
    case 'abssum'
        pdmcode = 1;
    case 'maxdiff'
        pdmcode = 2;
    case 'mindiff'
        pdmcode = 3;
    otherwise
        error('sltoolbox:invalidarg', ...
            'Invalid type of pwdiff computation: %s', type);
end

%% Compute

D = pwdiff_core(X1, X2, pdmcode);

author: lmthang


text:
function mv = mv_plotsMenu(mv,hfig);
%
% mv = mv_plotsMenu(mv,hfig);
%
% Attach menus providing plotting options for
% the MultiVoxel UI.
%
% ras, 04/05
if ieNotDefined('hfig')
    hfig = gcf;
end

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

mv.ui.plotMenu = uimenu('ForegroundColor', 'k', 'Label', 'Plot', ...
                        'Separator', 'on');

h(1) = uimenu(mv.ui.plotMenu,'Label','Whole Time Course For Each Voxel',...
        'Separator','off',...
        'Accelerator','J',...
        'CallBack','mv_selectPlotType(1); mv_plotWholeTc;');

h(2) = uimenu(mv.ui.plotMenu,'Label','Mean Amps: Voxel by Condition',...
        'Separator','off',...
        'Accelerator','M',...
        'CallBack','mv_selectPlotType(2); mv_plotAmps([],[2 3]);');

h(3) = uimenu(mv.ui.plotMenu,'Label','Mean Amps Subplots',...
        'Separator','off',...
        'Accelerator','B',...
        'CallBack','mv_selectPlotType(3); mv_ampsSubplots;');

h(4) = uimenu(mv.ui.plotMenu, 'Label', 'Mean Amps Subplots (remove means)', ...
        'Separator', 'off', ...
		'Accelerator', 'V', ...
        'CallBack','mv_selectPlotType(4); mv_ampsSubplots([],[],[],1);');

h(5) = uimenu(mv.ui.plotMenu,'Label','Mean Amps: Trial by Voxel by Condition',...
        'Separator','off',...
        'CallBack','mv_selectPlotType(5); mv_plotAmps([],[1 2 3]);');

h(6) = uimenu(mv.ui.plotMenu,'Label','Histograms',...
        'Separator','on',...
        'Accelerator','H',...
        'CallBack','mv_selectPlotType(6); mv_histograms; ');

h(7) = uimenu(mv.ui.plotMenu,'Label','Cross-Correlation Matrix',...
        'Separator','off',...
        'Accelerator','R',...
        'CallBack','mv_selectPlotType(7); mv_corrcoef; ');

h(8) = uimenu(mv.ui.plotMenu,'Label','Cross-Correlation Matrix (remove means)',...
        'Separator','off',...
        'Accelerator','T',...
        'CallBack','mv_selectPlotType(8); mv_corrcoef([], 1); ');

h(9) = uimenu(mv.ui.plotMenu,'Label','Visualize GLMs',...
        'Separator','on',...
        'Accelerator','G',...
        'CallBack','mv_selectPlotType(9); mv_visualizeGlm; ');

mv.ui.plotHandles = h;

return


text:
function dq = line2dq(v,r)

% LINE2DQ    Transforms a line position expressed in vector notation into
%            its dual quaternion representation.
%
%       DQ = LINE2DQ(V,R) transforms the line position, specified by the
%       line orientation V and the coordinates of any point of the line, R,
%       into its dual quaternion representation. V and R must have the same
%       size. V  (resp. R) is either a vector of size 3 or an array of size
%       3*N (column i represents the orientation (resp. coordinates of a
%       point) of Line i) where N is the number of lines. DQ is either a
%       vector of size 8, either an array of size (8*N) depending on V
%       format. Each column of DQ represents the dual quaternion
%       representation of the corresponding line position.
%
% See also POS2DQ, VEL2DQ, LINEVEL2DQ, DQ2LINE

sv = size(v);
sr = size(r);
if sv == [1 3]
    v = v.';
    sv = size(v);
end
if sr == [1 3]
    r = r.';
    sr = size(r);
end

% check that v and r have the same size
if sv ~= sr
    error('DualQuaternion:line2dquat:sizesDoNotMatch',...
        'Arrays v and r should be the same size. Size of v matrix is [%d %d] while size of r matrix is [%d %d]',sv(1),sv(2),sr(1),sr(2));
end
% if the format is wrong
if sv(1) ~= 3
    error('DualQuaternion:line2dquat:wrongsize',...
        '%d rows in the V and R array. It should be 3. ',sv(1));
end

% normalization of the axis vector (if necessary)
n = length(v(1,:));
n2 = sum(v.^2).^0.5;
n2 = repmat(n2,3,1);
v =v./n2;

% construction of the line dual quaternion
dq = sym(zeros(8,n));
dq(2:4,:) = v;
dq(6:8,:) = cross(r,v);

author: Sable

repo: mcbench-benchmarks


text:
function moebius_values_test ( )

%*****************************************************************************80
%
%% MOEBGIUS_VALUES_TEST demonstrates the use of MOEBIUS_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MOEBIUS_VALUES_TEST:\n' );
  fprintf ( 1, '  MOEBIUS_VALUES returns values of \n' );
  fprintf ( 1, '  the Moebius function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N         MU(N)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

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

    if ( n_data == 0 )
      break
    end

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

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [C, S] = QLiftDec2MinMin(X, N)
%-----------------------------------------------------------------------------
% QLiftDec2MinMin
% Multilevel 2-D decomposition by the lifting scheme and using quincunx grids
%
% The MinMin scheme has been proposed by Heijmans and Goutsias, see e.g.
%    H.J.A.M. Heijmans, J. Goutsias,
%    Multiresolution signal decomposition schemes.
%    Part 2: morphological wavelets.
%    CWI Report PNA-R9905, Amsterdam, 1999.
%    http://repository.cwi.nl:8888/cwi_repository/docs/IV/04/04625D.pdf
%
% Calls for: QLmaxlev,
%            storeQ1001, storeR,
%            getcolor01, getcolor10, getcolor00, getcolor11,
%            putcolor01, putcolor10, putcolor00, putcolor11.
% See also: QLiftRec2MaxMin
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: February 3, 2003.
% (c) 1999-2003 Stichting CWI, Amsterdam
%-----------------------------------------------------------------------------
%Firstly, check input data
%
if  isempty(X)
  error(' QLiftDec2MinMin - empty matrix ');
else
  if mod(N, 2) == 1
    error(' QLiftDec2MinMin - only an even number of levels is accepted ');
  end
  if QLmaxlev(size(X), 'minmin') < N
    error(' QLiftDec2MinMin - too many levels requested ');
  end
  if N < 2
    disp([' QLiftDec2MinMin - WARNING too few levels requested ' ...
          '-> empty decomposition ']);
  end
end
%
%Secondly, start decomposition
%
O = X; % For the sake of efficient use of memory this could be improved upon.
% We descend to coarser grids, integer lev indicates number of scale.
C = []; S = [];

for lev=1:2:N
%
   [nO, mO] = size(O);
   if ( nO < 3 ) || ( mO < 3)
     error(' QLiftDec2MinMin - too many levels ');
   end
   minO = min(min(O));
   maxO = max(max(O));
%  cmin = minO-(maxO-minO);
   cmax = maxO+(maxO-minO);
%
%  The Lifting Scheme proceeds from a rectangular grid
%  towards a quincunx grid.
%
%  Stage: predict
   A00 =getcolor00(O);
   A11 =getcolor11(O);
%  Quincunx grid Q0011 is the union of the values at .00 and .11: "even slots"
%  Quincunx grid Q1001 is the union of the values at .10 and .01: "odd slots"
   Q1001D01 = getcolor01(O) - synA01min(A11, A00, cmax);                % Y1
   Q1001D10 = getcolor10(O) - synA10min(A11, A00, cmax);                % Y1
%  At this point the union (quincunx) of Q1001D01 & Q1001D10
%  contains the DETAILS of O.
%
%  For the inverse transform Q1001D01 and Q1001D10 have to be stored:
   [C, S] = storeQ1001( Q1001D10, Q1001D01, lev, 'd', C, S);
%
%  Stage: update
   Q0011A00 = A00 + ...
        min(zeros(size(A00)), synA00min(Q1001D10, Q1001D01, cmax));     % X1
   clear A00;
   Q0011A11 = A11 + ...
        min(zeros(size(A11)), synA11min(Q1001D10, Q1001D01, cmax));     % X1
   clear A11 Q1001D10 Q1001D01;
%  At this point the union (quincunx) of Q0011A00 & Q0011A11
%  contains the updated APPROXIMATION of O, the DETAILS of O
%  were in the union (quincunx) of Q1001D01 & Q1001D10 (see above).
%
%  The Lifting Scheme proceeds by a subsequent step from quincunx
%  to rectangular grid.
%
%  Q0011 is split into the 11 colour with the "odd slots" and
%  the 00 colour with the "even slots".
%
%  Stage: predict
   DETAIL11 = Q0011A11 - synA11Qmin(Q0011A00, size(Q0011A11), cmax);        % Y2
   clear Q0011A11;
%  Stage: update
   APPROX00 = Q0011A00 + ...
     min(zeros(size(Q0011A00)), synA00Qmin(DETAIL11, size(Q0011A00), cmax));% X2
%
%  DETAIL11 presents the detail gridfunction w.r.t. Q0011
%  APPROX00 now represents the updated version of the approximation of Q0011
%
%  For the inverse transform DETAIL11 has to be stored:
   [C, S] = storeR( DETAIL11, lev+1, 'd', C, S);
   clear Q0011A00 DETAIL11;
%
%  At this point gridfunction DETAIL11 containing the DETAILS has been stored,
%  gridfunction APPROX00 contains the updated APPROXIMATION, on the (down-
%  sampled) rectangular grid and has to be stored as well if at the highest
%  scale.
%  Note that APPROX00 is downsampled onto a rectangular grid with dimensions of
%  half size of the original O.
   if lev+1 >= N
     [C, S] = storeR(APPROX00, lev+1, 'a', C, S);
%    It is obligatory that at least at one scale the Approximation has to be
%    stored or else the scheme cannot be inverted.
     clear APPROX00;
%    In the Lifting Scheme all scales have now been processed!
   else
%    We proceed to the next scale.
     O = APPROX00; clear APPROX00;
   end
end
%-----------------------------------------------------------------------------


text:
% planarcams ... alignment under assumption of planar camera arrangements
%
% [align,cam] = planarmove(in,cam,config,idx)
% in, cam, config ... see the main GOCAL script
% idx ... indexes of cameras which are suppose to be in one plane
%
% align ... structures aligned wit the specified world frame
%
% $Id: planarcams.m,v 1.1 2003/07/03 15:36:55 svoboda Exp $

function [align,cam] = planarcams(in,cam,config,idx)

% fit a plane to the reconstructed points and estimate normal

plane.n = planefit(in.Ce(:,idx)');

new.n = [0,0,1]'; % align the xy plane horizontally

rotaxis = cross(new.n,plane.n);
rotangle = acos( (plane.n'*new.n)/norm(plane.n)*norm(new.n) );

R = nfi2r(rotaxis,rotangle);
s = 3;
t = [0,0,1]' - s*R*mean(in.Xe(1:3,:)')';

align.simT.s = s;
align.simT.R = R;
align.simT.t = t;

[align.P, align.X]	= align3d(in.Pe,in.Xe,align.simT);
% save aligned data
if 1 % SAVE_STEPHI | SAVE_PGUHA
	[align.Cst,align.Rot] = savecalpar(align.P,config);
end
drawscene(align.X,align.Cst',align.Rot,61,'cloud','Graphical Output Validation: Aligned data, TopView',config.cal.cams2use);

set(gca,'CameraTarget',[0,0,0]);
set(gca,'CameraPosition',[0,0,1]);

figure(61),
% print -depsc graphevalaligned.eps
eval(['print -depsc ', config.paths.data, 'topview.eps'])

drawscene(align.X,align.Cst',align.Rot,62,'cloud','Graphical Output Validation: Aligned data, SideView',config.cal.cams2use);

set(gca,'CameraTarget',[0,0,0.9]);
set(gca,'CameraPosition',[2,0,0.9]);

figure(62),
% print -depsc graphevalaligned.eps
eval(['print -depsc ', config.paths.data, 'sideview.eps'])

return

author: strawlab

repo: MultiCamSelfCal


text:
function [f,g] = call_log_reg(v, y, X, lambda)
% Primal objective function (Omega = inv(Sigma))
% f = 0.5*(logdet(V) - logdet(Sigma) - tr(V*SigmaInv) - (m-mu)'*SigmaInv*(m-mu)
%     + L) - sum_d fb(mbar_d, vbar_d)
% where mbar = X*m, vbar = diag(X*V*X')
%
% Written by Emtiyaz,
% Modified by Wu Lin
  global m_cache;
  global V_cache;
  global iter_counter;
  global is_cache;
  global times_cache;
  iter_counter = iter_counter + 1;

  [D L] = size(X);

  Omega = diag(lambda);
  %Extract mean, Cholesky and bias
  m = v(1:L);
  idx = L + [1:L*(L+1)/2];
  U = triu(unpackcovariance(v(idx),L));

  % compute V
  V = U'*U;

  if is_cache==1
      m_cache(:,iter_counter) = m;
      V_cache(:,:,iter_counter) = V;
  end

  % compute kl and its gradient
  kl = 0.5*(2*sum(log(diag(U))) + sum(log(lambda(2:end))) - trace(V*Omega) -m'*(Omega*m) + L);

  % contribution from the bound
  mbar = X*m;
  vbar = sum(X.*(V*X')',2); % diag(X*V*X') efficient
  [fb, gmb, gvb] = E_log_p('bernoulli_logit', y, mbar, vbar, []);
  fb = -fb;
  gm_lvb = X'*(-gmb);
  gV_lvb = X'*bsxfun(@times, -gvb, X); %efficient X'*diag(gllp/2)*X;

  % final
  f = kl - sum(fb);
  gm = - Omega*m - gm_lvb;
  gU = diag(1./diag(U)) - triu(U*Omega) - triu(2*U*gV_lvb);

  g=[gm(:); gU(triu(ones(L))==1)];

  % return
  g=-g;
  f=-f;

  if iter_counter+1<=length(times_cache)
      times_cache(iter_counter+1) = toc;
  else
      times_cache(end) = toc;
  end

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
%                           solvePhaseMax.m
%
%  Implementation of the PhaseMax algorithm proposed in the paper using
%  FASTA. Note: For this code to run, the solver "fasta.m" must be in your
%  path.
%
%% I/O
%  Inputs:
%     A:    m x n matrix or a function handle to a method that
%           returns A*x.
%     At:   The adjoint (transpose) of 'A'. If 'A' is a function handle, 'At'
%           must be provided.
%     b0:   m x 1 real,non-negative vector consists of all the measurements.
%     x0:   n x 1 vector. It is the initial guess of the unknown signal x.
%     opts: A struct consists of the options for the algorithm. For details,
%           see header in solvePhaseRetrieval.m or the User Guide.
%
%     Note: When a function handle is used, the
%     value of 'At' (a function handle for the adjoint of 'A') must be
%     supplied.
%
%  Outputs:
%     sol:  n x 1 vector. It is the estimated signal.
%     outs: A struct consists of the convergence info. For details,
%           see header in solvePhaseRetrieval.m or the User Guide.
%
%
%  See the script 'testPhaseMaxGaussian.m' and 'testPhaseMaxFourier.m' for
%  two examples of proper usage of this function.
%
%% Notations
%
%
%% Algorithm Description
%  Solve the PhaseMax signal reconstruction problem
%         maximize <x0,x>
%         subject to |Ax|<=b0
%
%  The problem is solved by approximately enforcing the constraints using a
%  quadratic barrier function.  A continuation method is used to increase
%  the strength of the barrier function until a high level of numerical
%  accuracy is reached.
%  The objective plus the quadratic barrier has the form
%     <-x0,x> + 0.5*max{|Ax|-b,0}^2.
%
%  For a detailed explanation, see the PhaseMax paper referenced below. For
%  more details about FASTA, see the FASTA user guide, or the paper "A
%  field guide to forward-backward splitting with a FASTA implementation."
%
%% References
%  Paper Title:   PhaseMax: Convex Phase Retrieval via Basis Pursuit
%  Authors:       Tom Goldstein, Christoph Studer
%  arXiv Address: https://arxiv.org/abs/1610.07531
%
%  Copyright Goldstein & Studer, 2016.  For more details, visit
%  https://www.cs.umd.edu/~tomg/projects/phasemax/

function [sol, outs] = solvePhaseMax(A,At,b0,x0,opts)
    % Initialization
    m = length(b0);              % number of measurements
    n = length(x0);              % length of the unknown signal
    remainIters = opts.maxIters; % The remaining fasta iterations we have.
                                 % It's initialized to opts.maxIters.

    %  Normalize the initial guess relative to the number of measurements
    x0 = (x0/norm(x0(:)))*mean(b0(:))*(m/n)*100;

    %  re-scale the initial guess so that it approximately satisfies |Ax|=b
    sol = x0.* min(b0./abs(A(x0)));

    % Initialize values potentially computed at each round.
    ending = 0;       % Indicate whether any of the ending condition (except maxIters) has been met in FASTA.
    iter = 0;
    currentTime = [];
    currentResid = [];
    currentReconError = [];
    currentMeasurementError = [];

    % Initialize vectors for recording convergence information
    [solveTimes,measurementErrors,reconErrors,residuals] = initializeContainers(opts);

    %%  Define objective function components for the gradient descent method FASTA
    f = @(z) 0.5*norm(max(abs(z)-b0,0))^2;     % f(z) = 0.5*max{|x|-b,0}^2 : This is the quadratic penalty function
    gradf = @(z)  (sign(z).*max(abs(z)-b0,0)); % The gradient of the quadratic penalty

    % Options to hand to fasta
    fastaOpts.maxIters = opts.maxIters;
    fastaOpts.stopNow = @(x, iter, resid, normResid, maxResid, opts) ...
        processIteration(x, resid);  % Use customized stopNow in order to get
                                         % solveTime, residual and error at each iteration.
    fastaOpts.verbose=0;
    startTime = tic;                    % Start timer
    constraintError = norm(abs(A(sol))-b0);  % Keep track of the current error in the solution.
    while remainIters > 0 & ~ending     % Iterate over continuation steps
        g = @(x) -real(x0'*x);          % The linear part of the objective
        proxg = @(x,t) x+t*x0;          % The proximal operator of the linear objective
        fastaOpts.tol = norm(x0)/100;   % use a tighter tolerance when the solution is more exact
        % Call FASTA to solve the inner minimization problem
        [sol, fastaOuts] = fasta(A, At, f, gradf, g, proxg, sol, fastaOpts);  % Call a solver to minimize the quadratic barrier problem

        fastaOpts.tau = fastaOuts.stepsizes(end);     % Record the most recent stepsize for recycling.
        x0 = x0/10;                                   % do continuation - this makes the quadratic penalty stronger

        % Update the max number of iterations for fasta
        remainIters = remainIters - fastaOuts.iterationCount;
        fastaOpts.maxIters = min(opts.maxIters, remainIters);

        % Monitor convergence and check stopping conditions
        newConstraintError = norm(max(abs(A(sol))-b0,0));
        relativeChange = abs(constraintError-newConstraintError)/norm(b0);
        if relativeChange <opts.tol   % terminate when error reduction stalls
            break;
        end
        constraintError = newConstraintError;
    end

    % Create output according to the options chosen by user
    outs = generateOutputs(opts, iter, solveTimes, measurementErrors, reconErrors, residuals);

    % Display verbose output if specified
    if opts.verbose == 1
        displayVerboseOutput(iter, currentTime, currentResid, currentReconError, currentMeasurementError);
    end

    % Runs code upon each FASTA iteration. Returns whether FASTA should
    % terminate.
    function stop = processIteration(x, residual)
        iter = iter + 1;
        % Record convergence information and check stopping condition
        % If xt is provided, reconstruction error will be computed and used for stopping
        % condition. Otherwise, residual will be computed and used for stopping
        % condition.
        if ~isempty(opts.xt)
            xt = opts.xt;
            %  Compute optimal rotation
            alpha = (x(:)'*xt(:))/(x(:)'*x(:));
            x = alpha*x;
            currentReconError = norm(x-xt)/norm(xt);
            if opts.recordReconErrors
                reconErrors(iter) = currentReconError;
            end
        end

        if isempty(opts.xt)
            currentResid = residual;
        end

        if opts.recordResiduals
            residuals(iter) = residual;
        end

        currentTime = toc(startTime);                % Record elapsed time so far
        if opts.recordTimes
            solveTimes(iter) = currentTime;
        end
        if opts.recordMeasurementErrors
            currentMeasurementError = norm(abs(A(sol)) - b0) / norm(b0);
            measurementErrors(iter) = currentMeasurementError;
        end

        % Display verbose output if specified
        if opts.verbose == 2
            displayVerboseOutput(iter, currentTime, currentResid, currentReconError, currentMeasurementError);
        end

        % Test stopping criteria.
        stop = false;
        if currentTime >= opts.maxTime % Stop if we're run over the max runtime
            stop = true;
        end
        if ~isempty(opts.xt) % If true solution is specified, terminate when close to true solution
            assert(~isempty(currentReconError),'If xt is provided, currentReconError must be provided.');
            stop = stop || currentReconError < opts.tol;
            ending = stop;  % When true, this flag will terminate outer loop
        end
        stop = stop || residual < fastaOpts.tol;    % Stop FASTA is the tolerance is reached
    end

end

author: tomgoldstein

repo: phasepack-matlab


text:
function Asup = mg_ns_cd_setup(x,y,viscosity,flowsol, domain,lncoarse,lnfine,outbnd)
%mg_ns_cd_setup GMG convection-diffusion problem (Navier-Stokes)
%   Asup = mg_ns_cd_setup(x,y,viscosity,flowsol, domain,lncoarse,lnfine,outbnd)
%   input
%          x            x grid coordinates
%          y            y grid coordinates
%          viscosity    viscosity in NS equations
%          flowsol      current discrete velocity field
%          domain       domain parameter, 1 for cavity, 3 for step
%          lncoarse     coarse grid index, log2(nc) for nc x nc square grid
%          lnfine       finest grid index, log2(nf) for nf x nf square grid
%          outbnd       location of outflow boundary
%   output
%          Asup         discrete convection-diffusion operator for grid
%                       of level lncoarse
%
%   IFISS function: HCE; 18 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

% adaptation of AR mg_cd_setup allowing convection coefficients to come
% from flowsol rather than function call

if domain==1,
   n=length(x)-1; np=n/2; nq=n/4;
   nvtx=(n+1)*(n+1);
   [X,Y]=meshgrid(x,y);
   xx=reshape(X',nvtx,1);
   yy=reshape(Y',nvtx,1);
   xy=[xx(:),yy(:)];
% assembly process
   kx = 1;
   ky = 1;
   mel=0;
   for j=1:np
      for i=1:np
         mref=(n+1)*(ky-1)+kx;
         pref=(np+1)*(j-1)+i;
         mel=mel+1;
         nvv(1) = mref;
         nvv(2) = mref+2;
         nvv(3) = mref+2*n+4;
         nvv(4) = mref+2*n+2;
         nvv(5) = mref+1;
         nvv(6) = mref+n+3;
         nvv(7) = mref+2*n+3;
         nvv(8)=  mref+n+1;
         nvv(9)=  mref+n+2;
         npp(1) = pref;
         npp(2) = pref+1;
         npp(3) = pref+np+2;
         npp(4) = pref+np+1;
         mv(mel,1:9)=nvv(1:9);
         mp(mel,1:4)=npp(1:4);
         kx = kx + 2;
      end
      ky = ky + 2;
      kx = 1;
   end
%% compute boundary vertices and edges (four boundary edges)
   k1=find( xy(:,2)==-1 );
   e1=[]; for k=1:mel, if any(mv(k,5)==k1), e1=[e1,k]; end, end
   ef1=ones(size(e1));
%
   k2=find( xy(:,1)==1  & xy(:,2)<1   & xy(:,2) >-1);
   e2=[]; for k=1:mel, if any(mv(k,6)==k2), e2=[e2,k]; end, end
   ef2=2*ones(size(e2));
%
   k3=find( xy(:,2)==1 );
   e3=[]; for k=1:mel, if any(mv(k,7)==k3), e3=[e3,k]; end, end
   ef3=3*ones(size(e3));
%
   k4=find( xy(:,1)==-1 & xy(:,2)<1   & xy(:,2) >-1 );
   e4=[]; for k=1:mel, if any(mv(k,8)==k4), e4=[e4,k]; end, end
   ef4=4*ones(size(e4));
%
   bound=sort([k1;k2;k3;k4]);
   mbound=[e1',ef1';e2',ef2';e3',ef3';e4',ef4'];

elseif domain==3,
   n=length(y)-1; np=n/2; nq=n/4;
   nvtx= np*(np*1) + (outbnd*np+1)*(n+1);
% negative x-values
   xneg=x(1:n/2);yneg=y(1:n/2);
   xpos=x(n/2+1:end);ypos=y(n/2+1:end);
   [Xneg,Ypos]=meshgrid(xneg,ypos);
   xx=reshape(Xneg',np*(np+1),1);
   yy=reshape(Ypos',np*(np+1),1);
   xyleft=[xx(:),yy(:)];

% assembly process
   kx = 1;
   ky = 1;
   mel=0;
% loop over 2x2 macroelements
   for j=1:nq
      for i=1:nq
         mref=np*(ky-1)+kx;
         pref=nq*(j-1)+i;
         mel=mel+1;
         nvv(1) = mref;
         nvv(2) = mref+2;
         nvv(3) = mref+2*np+2;
         nvv(4) = mref+2*np;
         nvv(5) = mref+1;
         nvv(6) = mref+np+2;
         nvv(7) = mref+2*np+1;
         nvv(8)=  mref+np;
         nvv(9)=  mref+np+1;
         npp(1) = pref;
         npp(2) = pref+1;
         npp(3) = pref+nq+1;
         npp(4) = pref+nq;
         mv(mel,1:9)=nvv(1:9);
         mp(mel,1:4)=npp(1:4);
         kx = kx + 2;
      end
      ky = ky + 2;
      kx = 1;
   end

% correction along the internal boundary
   mref=2*np*(((outbnd+1)/2)*np+1)+1;
   pref=2*nq*(((outbnd+1)/2)*nq+1)+1;
   for mel=nq:nq:nq*nq;
      nvv=mv(mel,:);
      npp=mp(mel,:);
      nvv(2) = mref;
      nvv(3) = mref+(2*outbnd)*np+2;
      nvv(6) = mref+outbnd*np+1;
      npp(2) = pref;
      npp(3) = pref+outbnd*nq+1;
      mv(mel,1:9)=nvv(1:9);
      mp(mel,1:4)=npp(1:4);
      mref=mref+(2*outbnd)*np+2;
      pref=pref+outbnd*nq+1;
   end
% positive x_values
   [Xpos,Y]=meshgrid(xpos,y);
   xx=reshape(Xpos',(outbnd*np+1)*(n+1),1);
   yy=reshape(Y',(outbnd*np+1)*(n+1),1);
   xyright=[xx(:),yy(:)];
   xy=[xyleft;xyright];
%
   kx = 1;
   ky = 1;
   mel=nq*nq;
   for j=1:np
      for i=1:outbnd*nq
         mref = (outbnd*np+1)*(ky-1)+kx + np*(np+1);
         pref = (outbnd*nq+1)*(j-1)+i +  nq*(nq+1);
         mel=mel+1;
         nvv(1) = mref;
         nvv(2) = mref+2;
         nvv(3) = mref+(2*outbnd)*np+4;
         nvv(4) = mref+(2*outbnd)*np+2;
         nvv(5) = mref+1;
         nvv(6) = mref+outbnd*np+3;
         nvv(7) = mref+(2*outbnd)*np+3;
         nvv(8)=  mref+outbnd*np+1;
         nvv(9)=  mref+outbnd*np+2;
         npp(1) = pref;
         npp(2) = pref+1;
         npp(3) = pref+outbnd*nq+2;
         npp(4) = pref+outbnd*nq+1;
         mv(mel,1:9)=nvv(1:9);
         mp(mel,1:4)=npp(1:4);
         kx = kx + 2;
      end
      ky = ky + 2;
      kx = 1;
   end

%% compute boundary vertices and edges (five boundary edges)
   k1=find( xy(:,1) <0  & xy(:,2)==0 );
   e1=[]; for k=1:mel, if any(mv(k,5)==k1), e1=[e1,k]; end, end
   ef1=ones(size(e1));
%
   k2=find( xy(:,1)==0  & xy(:,2)<=0 );
   e2=[]; for k=1:mel, if any(mv(k,8)==k2), e2=[e2,k]; end, end
   ef2=4*ones(size(e2));
%
   k3=find( xy(:,1) >0  & xy(:,2)==-1);
   e3=[]; for k=1:mel, if any(mv(k,5)==k3), e3=[e3,k]; end, end
   ef3=ones(size(e3));
%
   k5=find( xy(:,2)==1 );
   e5=[]; for k=1:mel, if any(mv(k,7)==k5), e5=[e5,k]; end, end
   ef5=3*ones(size(e5));
%
   k6=find( xy(:,1)==-1 & xy(:,2)<1   & xy(:,2) >0 );
   e6=[]; for k=1:mel, if any(mv(k,8)==k6), e6=[e6,k]; end, end
   ef6=4*ones(size(e6));
%
   bound=sort([k1;k2;k3;k5;k6]);
   mbound=[e1',ef1';e2',ef2';e3',ef3';;e5',ef5';e6',ef6'];
end

%
%% set up matrices for Q1 approximation
[ev,ebound]=mg_q1grid(x,y,xy,mv,bound,mbound);
[A,N,M,epe,eph,epw,usolc,vsolc] = mg_ns_q1cd(xy,ev,flowsol,domain,lncoarse,lnfine,outbnd);
% SUPG
supg=0;
A = viscosity*A + N;
%% compute element peclet numbers
epe = epe/viscosity;
%% include streamline diffusion matrix (if necessary)
esupg=find(epe<=1); expe=epe;
if any(expe),
   expe=0.5*(1-1./expe);
   expe(esupg)=inf;
   epp=expe; epp(esupg)=0; epp=epp.*eph./epw;
   S=mg_ns_q1cd_supg(xy,ev,expe,eph,epw,usolc,vsolc);A=A+S;
end
%
%% impose zero boundary conditions
Asup = mg_zerobc(A,xy,bound);


text:
function dataset = vl_svmdataset(x, varargin)
% VL_SVMDATASET  Construct advanced SVM dataset structure
%   DATASET = VL_SVMDATASET(X) wraps the matrix X into a dataset
%   structure. All the formats supported by VL_SVMTRAIN() are
%   suppported here as well.
%
%   DATASET = VL_SVMDATASET(X, 'homkermap', HOM) where HOM is a
%   structure with zero or more of the following fields: HOM.kernel,
%   HOM.order, HOM.window, HOM.gamma, HOM.period, See VL_HOMKERMAP()
%   for a discussion of these parameters. The effect is to use the
%   homogeneous kernel map to expand the data on the fly. This is
%   equivalent (albeit somewhat slower) to running VL_HOMKERMAP()
%   prior to learning. The main advantage is the significant memory
%   saving as the expanded data needs not to be stored in memory.
%
%   See:: VL_SVMTRAIN(), VL_HOMKERMAP(), VL_HELP().

% Author: Daniele Perrone and Andrea Vedaldi

opts.homkermap = [] ;
opts = vl_argparse(opts, varargin) ;

if issparse(x)
  error('X is a sparse matrix (sparse data will support be added in a future version.)') ;
end

if ~isa(x,'single') && ~isa(x,'double')
  error('X is neither SINGLE nor DOUBLE.') ;
end

sz = size(x) ;
if numel(sz) > 2
  error('X has more than two dimensions.') ;
end

dataset.data = x ;

if isstruct(opts.homkermap)
  if numel(opts.homkermap) ~= 1
    error('HOMKERMAP is not a singleton.') ;
  end
  ok = all(ismember(fieldnames(opts.homkermap), ...
                    {'kernel', 'order', 'window', 'gamma', 'period'})) ;
  if ~ok
    error('HOMKERMAP contains an unsupported field name.') ;
  end
  dataset.homkermap = opts.homkermap ;
elseif ~isempty(opts.homkermap)
    error('HOMKERMAP is not a structure nor []') ;
end

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [ra, dec] = apstar (tjd, n, rai, deci, pmra, pmdec, parlax, radvel)

% apparent place of a star

%  tjd    = tt julian date for apparent place (in)

%  n      = body identification number for the earth (in) (no longer used)

%  rai    = icrs right ascension in hours (in)

%  deci   = icrs declination in degrees (in)

%  pmra   = icrs proper motion in ra in milliarcseconds/year (in)

%  pmdec  = icrs proper motion in dec in milliarcseconds/year (in)

%  parlax = parallax in milliarcseconds (in)

%  radvel = radial velocity in kilometers/second (in)

%  ra     = apparent right ascension in hours (out)

%  dec    = apparent declination in degrees (out)

% note: coordinate system for output ra and dec is equator and equinox of date

% ported from NOVAS 3.1

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

locatn = 0;

icoord = 1;

ttjd = tjd;

object = 'star';

star(1) = rai;

star(2) = deci;

star(3) = pmra;

star(4) = pmdec;

star(5) = parlax;

star(6) = radvel;

observ = zeros(6, 1);

skypos = place (ttjd, object, locatn, icoord, star, observ);

ra = skypos(4);

dec = skypos(5);

author: goGPS-Project


text:
function transforms2d(varargin)
%TRANSFORMS2D Description of functions operating on transforms
%
%   By 'transform' we mean an affine transform. A planar affine transform
%   can be represented by a 3x3 matrix.
%
%   Example
%   % create a translation by the vector [10 20]:
%   T = createTranslation([10 20])
%   T =
%        1     0    10
%        0     1    20
%        0     0     1
%
%
%   See also:
%   createTranslation, createRotation, createRotation90, createScaling
%   createHomothecy, createLineReflection, createBasisTransform
%   transformPoint, transformVector, transformLine, transformEdge
%   rotateVector, fitAffineTransform2d
%   polynomialTransform2d, fitPolynomialTransform2d
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2008-10-13,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

help('transforms2d');

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function ar=lpczz2ar(zz)
%LPCZZ2AR Convert z-place poles to ar coefficients AR=(ZZ)
% The complex poles must occur in complex conjugate pairs
% but the order is unimportant.

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: lpczz2ar.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p]=size(zz);
ar=zeros(nf,p+1);
for k=1:nf
  ar(k,:)=real(poly(zz(k,:)));
end

author: jtkim-kaist

repo: VAD

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


text:
function b = r83_vxm ( n, a, x )

%*****************************************************************************80
%
%% R83_VXM multiplies a vector by an R83 matrix.
%
%  Discussion:
%
%    The r83 storage format is used for a tridiagonal matrix.
%    The superdiagonal is stored in entries (1,2:N), the diagonal in
%    entries (2,1:N), and the subdiagonal in (3,1:N-1).  Thus, the
%    original matrix is "collapsed" vertically into the array.
%
%  Example:
%
%    Here is how a R83 matrix of order 5 would be stored:
%
%       *  A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54  *
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the linear system.
%
%    Input, real A(3,N), the R83 matrix.
%
%    Input, real X(N), the vector to be multiplied by A'.
%
%    Output, real B(N), the product A' * x.
%
  b(1:n)   =            a(2,1:n)   .* x(1:n);
  b(1:n-1) = b(1:n-1) + a(3,1:n-1) .* x(2:n);
  b(2:n)   = b(2:n)   + a(1,2:n)   .* x(1:n-1);

  return
end


text:
function disp(x)
%DISP         Display function for pop-up windows in debugger
%

% written  04/26/13     S.M. Rump
%

  display(x,longprecision);

repo: OpenMAS


text:
function ker = set_grid(ker, nthreads)

global data;

block_size = min(nthreads, 512);
grid_size = min([ceil(nthreads / block_size), 1], data.gpu.max_grid_size);
if any(ker.block_size ~= block_size) || any(ker.grid_size ~= grid_size)
    %gpu_sync();
    ker.block_size = block_size;
    ker.grid_size = grid_size;
    ker.ker.ThreadBlockSize = block_size;
    ker.ker.GridSize = grid_size;
    data.gpu.kernels.(ker.name) = ker;
end

author: jimmy-ren


text:
function [gamma, loglik, marginals, marginalsT] = bk_ff_fb(prior, transmat, obslik, filter_only, hnodes, ns)
% BK_FF_FB Fully factored Boyen-Koller version of forwards-backwards
% [gamma, loglik, marginals, marginalsT] = bk_ff_hmm(prior, transmat, obslik, filter_only, hnodes, ns)

ss  = length(ns);
S = length(prior);
T = size(obslik, 2);
marginals = cell(ss,T);
marginalsT = cell(ss,T);
scale = zeros(1,T);
alpha = zeros(S, T);

transmat2 = transmat';
for t=1:T
  if t==1
    [alpha(:,t), scale(t)] = normalise(prior(:) .* obslik(:,t));
  else
    [alpha(:,t), scale(t)] = normalise((transmat2 * alpha(:,t-1)) .* obslik(:,t));
  end
  [marginals(:,t), marginalsT(:,t)] = project_joint_onto_marginals(alpha(:,t), hnodes, ns);
  alpha(:,t) = combine_marginals_into_joint(marginalsT(:,t), hnodes, ns);
  %fprintf('alpha t=%d\n', t);
  %celldisp(marginals(1:8,t))
end
loglik = sum(log(scale));

if filter_only
  gamma = alpha;
  return;
end

beta = zeros(S,T);
gamma = zeros(S,T);
t = T;
beta(:,t) = ones(S,1);
gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
[marginals(:,t), marginalsT(:,t)] = project_joint_onto_marginals(gamma(:,t), hnodes, ns);

for t=T-1:-1:1
  b = beta(:,t+1) .* obslik(:,t+1);
  beta(:,t) = normalise((transmat * b));
  [junk, tempT] = project_joint_onto_marginals(beta(:,t), hnodes, ns);
  beta(:,t) = combine_marginals_into_joint(tempT, hnodes, ns);
  %gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
  %[marginals(:,t), marginalsT(:,t)] = project_joint_onto_marginals(gamma(:,t), hnodes, ns);
end

gamma2 = zeros(S,T);
for t=T-1:-1:1
  b = beta(:,t+1) .* obslik(:,t+1);
  xi(:,:,t) = normalise((transmat .* (alpha(:,t) * b')));
  if t==T-1
    gamma2(:,T) = sum(xi(:,:,T-1), 1)';
  end
  gamma2(:,t) = sum(xi(:,:,t), 2);
  [marginals(:,t), marginalsT(:,t)] = project_joint_onto_marginals(gamma2(:,t), hnodes, ns);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [handles,msh] = dtiMrMeshAddImages(handles,msh,imOrigin,xIm,yIm,zIm,textures)
%   Add the images to the mrMesh window
%
%   [handles,msh] = ...
%      dtiMrMeshAddImages(handles,msh,imOrigin,xIm,yIm,zIm,textures)
%
% HISTORY:
% Authors: Dougherty, Wandell
% 2005.01.17 RFD: fixed minor inefficiency in how mrMesh params were
% passed.
%
% Stanford VISTA Team

if ieNotDefined('handles'), error('dtiFiberUI handles required.'); end
if ieNotDefined('msh'), msh = dtiGet(handles,'mrMesh'); end
if ieNotDefined('imOrigin'), error('imOrigin required.'); end
if ~exist('xIm','var'), error('xIm required.'); end
if ~exist('yIm','var'), error('yIm required.'); end
if ~exist('zIm','var'), error('zIm required.'); end
if ieNotDefined('textures'), textures = []; end

if ~isempty(xIm) && ~iscell(xIm), tmp{1} = xIm; clear xIm; xIm = tmp; end
if ~isempty(yIm) && ~iscell(yIm), tmp{1} = yIm; clear yIm; yIm = tmp;  end
if ~isempty(zIm) && ~iscell(zIm), tmp{1} = zIm; clear zIm; zIm = tmp; end
if ~iscell(imOrigin), tmp = imOrigin; clear imOrigin; imOrigin{1} = tmp; end

xSize = length(xIm);
ySize = length(yIm);
zSize = length(zIm);

if xSize + ySize + zSize == 0, return, end
if xSize == 0, xSize = max([xSize,ySize,zSize]); end
if ySize == 0, ySize = max([xSize,ySize,zSize]); end
if zSize == 0, zSize = max([xSize,ySize,zSize]); end

if (xSize ~= ySize) || (xSize ~= zSize) || (ySize ~= zSize)
    error('You should have to same number of images to display');
end

nImages = xSize;

if (length(imOrigin) ~= nImages)
    error('You should have an origin specified for each image');
end

% sbOffset = 10;

% imSize must be 2^n (openGL requirement)
%
imSize = 0;
for i = 1:nImages
    if ~isempty(xIm), imSize = max([imSize,2.^ceil(log2(size(xIm{i}')))]); end;
    if ~isempty(yIm), imSize = max([imSize,2.^ceil(log2(size(yIm{i})))]); end;
    if ~isempty(zIm), imSize = max([imSize,2.^ceil(log2(size(zIm{i}')))]); end;
end

imSize = [imSize imSize];
imTemplate.class = 'image';
imTemplate.width = imSize(1);
imTemplate.height = imSize(2);
imTemplate.tex_width = imSize(2);
imTemplate.tex_height = imSize(1);

% 1 is opaque. 0 is completely transparent.
% We need to set this parameter in the GUI at some point.  Leave it like
% this for now so Bob and Michal see the effect and comment.
imTemplate.transparency = msh.transparency;

if ~isempty(xIm)
    for i = 1:nImages
        im = imTemplate;
        [id,s,r] = mrMesh(msh.host, msh.id, 'add_actor', im);
        msh.imgActors(3*i) = r.actor;
        im.actor = r.actor;
        sz = size(xIm{i}');
        pos = floor((imSize-sz)./2)+1;
        imData = zeros(imSize);
        imData(pos(1):pos(1)+sz(1)-1, pos(2):pos(2)+sz(2)-1) = xIm{i}';
        im.texture = repmat(imData(:)', 4, 1);

        if isempty(textures)
            im.texture(4,:) = im.transparency*ones(size(im.texture(1,:)));
        else
            textures.xIm{i} = textures.xIm{i}';
            if(any(size(textures.xIm{i})~=imSize))
                imData = zeros(imSize);
                imData(pos(1):pos(1)+sz(1)-1, pos(2):pos(2)+sz(2)-1) = textures.xIm{i};
                im.texture(4,:) = imData(:);
            else
                im.texture(4,:) = textures.xIm{i}(:);
            end
        end

        im.rotation = [0 0 1; 0 1 0; 1 0 0];
        im.origin = imOrigin{i}.x;
        [id,s,r] = mrMesh(msh.host, msh.id, 'set', im);
    end
end

if~isempty(yIm)
    for i = 1:nImages
        im = imTemplate;
        [id,s,r] = mrMesh(msh.host, msh.id, 'add_actor', im);
        msh.imgActors(2 + 3*(i - 1)) = r.actor;
        im.actor = r.actor;
        sz = size(yIm{i});
        pos = floor((imSize-sz)./2)+1;
        imData = zeros(imSize);
        imData(pos(1):pos(1)+sz(1)-1, pos(2):pos(2)+sz(2)-1) = yIm{i};
        im.texture = repmat(imData(:)', 4, 1);

        if isempty(textures)
            im.texture(4,:) = im.transparency*ones(size(im.texture(1,:)));
        else
            %textures.yIm{i} = textures.yIm{i}';
            if(any(size(textures.yIm{i})~=imSize))
                imData = zeros(imSize);
                imData(pos(1):pos(1)+sz(1)-1, pos(2):pos(2)+sz(2)-1) = textures.yIm{i};
                im.texture(4,:) = imData(:);
            else
                im.texture(4,:) = textures.yIm{i}(:);
            end
        end

        im.rotation = [1 0 0; 0 0 1; 0 1 0];
        im.origin = imOrigin{i}.y;
        [id,s,r] = mrMesh(msh.host, msh.id, 'set', im);
    end
end

if ~isempty(zIm)
    for i = 1:nImages
        im = imTemplate;
        [id,s,r] = mrMesh(msh.host, msh.id, 'add_actor', im);
        msh.imgActors(1 + 3*(i - 1)) = r.actor;
        im.actor = r.actor;
        sz = size(zIm{i}');        pos = floor((imSize-sz)./2)+1;
        imData = zeros(imSize);
        imData(pos(1):pos(1)+sz(1)-1, pos(2):pos(2)+sz(2)-1) = zIm{i}';
        im.texture = repmat(imData(:)', 4, 1);

        if isempty(textures)
            im.texture(4,:) = im.transparency*ones(size(im.texture(1,:)));
        else
            textures.zIm{i} = textures.zIm{i}';
            if(any(size(textures.zIm{i})~=imSize))
                imData = zeros(imSize);
                imData(pos(1):pos(1)+sz(1)-1, pos(2):pos(2)+sz(2)-1) = textures.zIm{i};
                im.texture(4,:) = imData(:);
            else
                im.texture(4,:) = textures.zIm{i}(:);
            end
        end

        im.rotation = [1 0 0; 0 1 0; 0 0 1];
        im.origin = imOrigin{i}.z;
        [id,s,r] = mrMesh(msh.host, msh.id, 'set', im);
        %msh.img(1) = im;
    end
end

return;

author: vistalab


text:
function [ cluster_center, seed ] = cluster_initialize_5 ( dim_num, ...
  point_num, cluster_num, point, seed )

%*****************************************************************************80
%
%% CLUSTER_INITIALIZE_5 initializes the cluster centers to random values.
%
%  Discussion:
%
%    In this case, each cluster center is a random convex combination
%    of the data points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the number of spatial dimensions.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer CLUSTER_NUM, the number of clusters.
%
%    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates
%    of the points.
%
%    Input, integer SEED, a seed for the random
%    number generator.
%
%    Output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM),
%    the coordinates of the cluster centers.
%
%    Output, integer SEED, a seed for the random
%    number generator.
%

%
%  Get a PxC block of random factors.
%
  [ factor, seed ] = r8mat_uniform_01 ( point_num, cluster_num, seed );
%
%  Make each column of factors have unit sum.
%
  for j = 1 : cluster_num
    column_sum = sum ( factor(1:point_num,j) );
    factor(1:point_num,j) = factor(1:point_num,j) / column_sum;
  end
%
%  Set centers = points * factors.
%
  cluster_center(1:dim_num,1:cluster_num) = ...
    point(1:dim_num,1:point_num) * factor(1:point_num,1:cluster_num);

  return
end


text:
function [F] = mci_nmm_r2p2_dfdp (x,u,P,M)
% Parameter Jacobian for two region, two parameter NMM
% FORMAT [F] = mci_nmm_r2p2_dfdp (x,u,P,M)
%
% x         State
% u         Inputs
% P         Parameters
% M         Model structure
%
% F         F(i,j) = df(x)_i/dp_j
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny and Biswa Sengupta
% $Id: mci_nmm_r2p2_dfdp.m 6548 2015-09-11 12:39:47Z will $

% 18 state variables, 2 parameters
F=zeros(18,2);

curr_P=M.can_P; % Canonical parameter set

w21=P(1);
w12=P(2);

% 2 free parameters
curr_P.A{1}(2,1)=P(1); % Forward connection, w_21
curr_P.A{2}(1,2)=P(2); % Backward connection, w_12

P=curr_P;

% default parameters
E = [1 1/2 1/8]*32;         % extrinsic rates (forward, backward, lateral)
D = [2 16];                 % propogation delays (intrinsic, extrinsic)
H = [4 32];                 % receptor densities (excitatory, inhibitory)
T = [8 16];                 % synaptic constants (excitatory, inhibitory)
R = [2 1]/3;                % parameters of static nonlinearity

% neuronal states into matrix form; x(r,:) for region r
x = spm_unvec(x,M.x);
% extrinsic delays
De = D(2).*exp(P.D)/1000;

% delayed pyramidal cell activity
pyr_ext = presynaptic (x(:,9)-De*(x(:,5)-x(:,6)),P,R);

Te    = T(1)/1000*exp(P.T(:,1));         % excitatory time constants
He    = H(1)*exp(P.G(:,1));              % excitatory receptor density

HeTe=He/Te;

% Effect of forward connection on region 2 stellate cells
F(8,1)=HeTe*pyr_ext(1)*E(1)*w21*exp(w12);
% Effect of backward connection on region 1 pyramidal cells
F(9,2)=HeTe*pyr_ext(2)*E(2)*w12*exp(w12);
% Effect of backward connection on region 1 inhibitory cells
F(15,2)=HeTe*pyr_ext(2)*E(2)*w12*exp(w12);
end

function [S] = presynaptic (x,P,R)

% pre-synaptic inputs: s(V)
%--------------------------------------------------------------------------
R     = R.*exp(P.S);
S     = 1./(1 + exp(-R(1)*(x - R(2)))) - 1./(1 + exp(R(1)*R(2)));

end

author: spm


text:
%LINES5D  Generates three 5-dimensional lines
%
%	A = LINES5D(N);
%
% Generates a data set of N points, on 3 non-crossing, non-parallel lines
% in 5 dimensions.
%
% If N is a vector of sizes, exactly N(I) objects are generated
% for class I, I = 1,2.Default: N = [50 50 50].
%
% See also DATASETS, PRDATASETS

% Copyright: E. Pekalska, R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: lines5d.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function data = lines5d(N)
			if nargin< 1, N = [50 50 50]; end

	N = genclass(N,ones(1,3)/3);
	n1 = N(1);
	n2 = N(2);
	n3 = N(3);

  s1 = [0 0 0 1 0];
  s2 = [1 1 1 0 0];
  s3 = [0 1 0 1 0];
  s4 = [1 1 1 1 1];
  s5 = [0 1 1 0 1];
  s6 = [1 0 1 1 1];
  c1 = [0:1/(n1-1):1]';
  c2 = [0:1/(n2-1):1]';
  c3 = [0:1/(n3-1):1]';
  a  = c1*s1 + (1-c1)*s2;
  a  = [a; c2*s3 + (1-c2)*s4];
  a  = [a; c3*s5 + (1-c3)*s6];

	data = prdataset(a,genlab(N));
	data = setname(data,'5D Lines');

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% Test file for chebtech/roots.m

function pass = test_roots(pref)

if ( nargin < 1 )
    pref = chebtech.techPref();
end

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    %% Test roots of a bessel function:
    map = @(x) (x+1)*50;
    f = testclass.make(@(x) besselj(0, map(x)), [], pref);
    r = map(roots(f));
    exact = [   2.40482555769577276862163; 5.52007811028631064959660
                8.65372791291101221695437; 11.7915344390142816137431
                14.9309177084877859477626; 18.0710639679109225431479
                21.2116366298792589590784; 24.3524715307493027370579
                27.4934791320402547958773; 30.6346064684319751175496
                33.7758202135735686842385; 36.9170983536640439797695
                40.0584257646282392947993; 43.1997917131767303575241
                46.3411883716618140186858; 49.4826098973978171736028
                52.6240518411149960292513; 55.7655107550199793116835
                58.9069839260809421328344; 62.0484691902271698828525
                65.1899648002068604406360; 68.3314693298567982709923
                71.4729816035937328250631; 74.6145006437018378838205
                77.7560256303880550377394; 80.8975558711376278637723
                84.0390907769381901578795; 87.1806298436411536512617
                90.3221726372104800557177; 93.4637187819447741711905
                96.6052679509962687781216; 99.7468198586805964702799 ];
    pass(n, 1) = norm(r-exact,Inf) < 1e1*length(f)*eps;

    %% Test roots of an oscillatory function:
    k = 500;
    f = testclass.make(@(x) sin(pi*k*x), [], pref);
    r = roots(f);
    pass(n, 2) = norm(r-(-k:k)'/k, inf) < length(f)*eps;

    %% Test a perturbed polynomial:
    f = testclass.make( @(x) (x-.1).*(x+.9).*x.*(x-.9) + 1e-14*x.^5, ...
        [], pref);
    r = roots(f);
    pass(n, 3) = length(r) == 4 && norm(feval(f, r), inf) < ...
        1e2*length(f)*eps;

    %% Test a some simple polynomials:
    f = testclass.make([-1 ; 1], [], pref);
    r = roots(f);
    pass(n, 4) = all( r == 0 );

    f = testclass.make([1 ; 0 ; 1]);
    r = roots(f);
    pass(n, 5) = numel(r) == 2 && (norm(r, inf) < eps);

    %% Test some complex roots:
    f = testclass.make(@(x) 1 + 25*x.^2, [], pref);
    r = roots(f, 'complex', 1);

    pass(n, 6) = norm( r - [1i ; -1i]/5, inf) < 10*eps;

    f = testclass.make(@(x) (1 + 25*x.^2).*exp(x), [], pref);
    r = roots(f, 'complex', 1, 'prune', 1);

    pass(n, 7) = norm( r - [1i ; -1i]/5, inf) < 10*length(f)*eps;

    f = testclass.make(@(x) sin(100*pi*x));
    r1 = roots(f, 'complex', 1, 'recurse', 0);
    r2 = roots(f, 'complex', 1);

    pass(n, 8) = numel(r1) == 201 && numel(r2) >= 213;

    %% Test an array-valued function:
    f = testclass.make(@(x) [sin(pi*x), cos(pi*x)], [], pref);
    r = roots(f);
    r2 = [-1 0 1 -.5 .5 NaN].';
    pass(n, 9) = all( r(:) - r2 < 10*length(f)*eps | isnan(r2) );

    % Adding test for 'qz' flag:
    f = testclass.make(@(x) 1e-10*x.^3 + x.^2 - 1e-12, [], pref);
    r = roots(f, 'qz', 1);
    pass(n, 10) = ~isempty( r );
    pass(n, 11) = norm(feval(f, r), inf) < 10*eps;

    % Add a rootfinding test for low degree non-even functions:
    f = testclass.make(@(x) (x-.5).*(x-1/3), [], pref);
    r = roots(f, 'qz', 1);
    pass(n, 12) = norm(feval(f, r), inf) < eps;
end

end

author: chebfun


text:
% this script compares hedging based on Black-Scholes deltas with Factors on Demand hedging
% see Meucci, A. (2010) "Factors on Demand", Risk, 23, 7, p. 84-89
% available at http://ssrn.com/abstract=1565134

clc; clear; close all;
%%%%%%%%%%%%%%%%%%%
% inputs
tau_tilde=5;    % estimation step (days)
tau=40;         % time to horizon (days)
Time2Mats=[100 150 200 250 300];    % current time to maturity of call options in days
Strikes = [850 880 910 940 970];    % strikes of call options, same dimension as Time2Mat

r_free=0.04;    % risk-free rate
J=10000;        % number of simulations

%%%%%%%%%%%%%%%%%%%
% load underlying and volatility surface
load('DB_ImplVol');
numCalls = length(Time2Mats);
timeLength = length(spot);
numSurfPoints = length(days2Maturity)*length(moneyness);

%%%%%%%%%%%%%%%%%%%
% estimate invariant distribution assuming normality
% variables in X are changes in log(spot) and changes in log(imp.vol)
% evaluated at the 'numSurfPoints' points on the vol surface (vectorized).
X = zeros(timeLength-1,numSurfPoints+1);
% log-changes of underlying spot
X(:,1) = diff(log(spot));

% log-changes of implied vol for different maturities
impVolSeries = reshape(impVol,timeLength,numSurfPoints);
for i = 1:numSurfPoints,
    X(:,i+1) = diff(log(impVolSeries(:,i)));
end
muX = mean(X);
SigmaX = cov(X,1);

%%%%%%%%%%%%%%%%%%%
% project distribution to investment horizon
muX = muX*tau/tau_tilde;
SigmaX = SigmaX*tau/tau_tilde;

%%%%%%%%%%%%%%%%%%%
% linearly interpolate the vol surface at the current time to obtain
% implied vol for the given calls today, and price the calls
spot_T = spot(end);
volSurf_T = squeeze(impVol(end,:,:));
time2Mat_T = Time2Mats;
moneyness_T = Strikes/spot_T;
impVol_T = interpne(volSurf_T,[time2Mat_T',moneyness_T'],{days2Maturity,moneyness})'; % function by John D'Errico
callPrice_T = BlackScholesCall(spot_T,Strikes,r_free,impVol_T,Time2Mats/252);

%%%%%%%%%%%%%%%%%%%
% generate simulations at horizon
X_ = mvnrnd(muX,SigmaX,J);

% interpolate vol surface at horizon for the given calls
spot_ = spot_T*exp(X_(:,1));
impVol_ = zeros(J,numCalls);
for j = 1:J,
    volSurf = volSurf_T.*exp(reshape(X_(j,2:end),length(days2Maturity),length(moneyness)));
    time2Mat_ = Time2Mats-tau;
    moneyness_ = Strikes/spot_(j);
    impVol_(j,:) = interpne(volSurf,[time2Mat_',moneyness_'],{days2Maturity,moneyness})';  % function by John D'Errico
end

% price the calls
callPrice_ = zeros(J,numCalls);
for i = 1:numCalls,
    callPrice_(:,i) = BlackScholesCall(spot_,Strikes(i),r_free,impVol_(:,i),time2Mat_(i)/252);
end

% linear returns of the calls
Rc = callPrice_./repmat(callPrice_T,J,1) - 1;
% linear returns of the underlying
Rsp = spot_./spot_T - 1;

%%%%%%%%%%%%%%%%%%%
% compute the OLS linear (affine) model: Rc = a + b*Rsp + U
Z = [ones(J,1), Rsp];
olsLoadings = (Rc'*Z)/(Z'*Z);
a = olsLoadings(:,1);
b = olsLoadings(:,2);

%%%%%%%%%%%%%%%%%%%
% compute Black-Scholes delta and cash held in replicating portfolio
[callPrice_T,delta_T,cash_T] = BlackScholesCall(spot_T,Strikes,r_free,impVol_T,Time2Mats/252);
a_bs = cash_T./callPrice_T*r_free*tau/252;
b_bs = (delta_T./callPrice_T.*spot_T)';
fprintf('OLS: a = [%s\t]\n',sprintf('\t%7.4f',a'))
fprintf('B-S: a = [%s\t]\n',sprintf('\t%7.4f',a_bs'))
fprintf('OLS: b = [%s\t]\n',sprintf('\t%7.4f',b'))
fprintf('B-S: b = [%s\t]\n',sprintf('\t%7.4f',b_bs'))

for i = 1:numCalls
    figure
    plot(Rsp,Rc(:,i),'.')
    xlabel('return underlying')
    ylabel('return call option')
end


text:
% delete mesh, algorithm, or everything at once
% Danil Kirsanov, 09/2007

function object = geodesic_delete(object)

global geodesic_library;
if ~libisloaded(geodesic_library)       %everything is already cleared
    object = [];
    return;
end

if nargin == 0          % the simplest way to delete everything is to unload the library
   unloadlibrary(geodesic_library);
else
    if libisloaded(geodesic_library)
        if strcmp(object.object_type, 'mesh')
            calllib(geodesic_library, 'delete_mesh', object.id);      % delete mesh and corresponding algorithms
        else
            calllib(geodesic_library, 'delete_algorithm', object.id); % delete a single algorithm
        end
    end
end
object = [];

author: alecjacobson

repo: gptoolbox


text:
function [mu, s2, MU, SIG2] = simulGPmc(hyp, inf, meanfunc, cov, lik, input, target, test, lag, Nsamples)
% simulGPmc - Simulation of the dynamic GP model, where the output variance is
% propagated using the Monte Carlo method
%
%% Syntax
%  [mu, s2, MU, SIG2] = simulGPmc(hyp, inf, mean, cov, lik, input, target, test,
%  lag, Nsamples)
%
%% Description
% Idea: at every time step the output of GP model is approximated with
% Nsamples samples, which are used as the future inputs of the GP model.
% Samples are re-used if necessary (ie. y(k-1) for y(k-2) if lag=2 etc.)
% Uses routines gpx and gmx_sample.
%
% Input:
% * hyp      ... the structure of optimized hyperparameters
% * inf      ... the function specifying the inference method
% * meanfunc ... the prior mean function
% * cov      ... the specified covariance function, see help covFun for more info
% * lik      ... the likelihood function
% * input    ... the input part of the training data,  NxD matrix
% * target   ... the output part of the training data (ie. target), Nx1 vector
% * test     ... the input matrix for simulation, kxD vector, see
%                construct.m for more info
% * lag      ... the order of the model (number of used lagged outputs)
% * Nsamples ... the number of samples used in algorithm (ie. runs of simulation)
%
% Output:
% * mu    ... the mean predicted output
% * s2    ... the associated variances (with noise variances)
% * MU    ... the matrix of all predicted means, kxNsamples
% * SIG2  ... the associated predicted variances
%
% See also:
% gpx, gmx_sample, simulGPnaive
%
% Examples:
% demo_example_gp_simulation
%
%%
% * Written by J. Prikryl, November 2010
% * Based on the work of C.E. Rasmussen, A. Girard, K. Azman.
%

% meanfunc ... 'mean' is used as a matlab core function in this file

Ndx = 800;
DSig = 3;
num_iters = length(test);
sum_time  = 0;

[N, D] = size(input);
PDF = zeros(Nsamples,lag);

% Preallocate mu and s2
mu = zeros ( num_iters, 1 );
s2 = zeros ( num_iters, 1 );

% 1st step - input is a point
test_ = test(1,:);
[mu(1), s2(1), post] = gpx(hyp, inf, meanfunc, cov, lik, input, target, test_);

MU(1,:) = mu(1)*ones(1,Nsamples);
SIG2(1,:) = s2(1)*ones(1,Nsamples);

pdf = gmx_sample(mu(1),s2(1),Nsamples);
PDF(:,lag) = pdf;
% instead of a cycle, create one monstrous test_ matrix
test_ = repmat( test(2,:), Nsamples, 1 );
test_(:,lag) = pdf;
[MU(2,:), SIG2(2,:), post] = gpx(hyp, inf, meanfunc, cov, lik, input, target, test_, post);
mu(2,1) = mean(MU(2,:));
s2(2,1) = mean(SIG2(2,:)) + mean((MU(2,:)-mu(2)).^2);

% steps from 3 on ...
for k=3:num_iters

    if(mod(k,50)==0 || k==3)
        disp(['simulGPmc, step: ',int2str(k),'/',int2str(length(test))]);
    end

    if(k>lag)
        % samples of previous distributions
        for jj=1:lag-1
            PDF(:,jj) = PDF(:,jj+1);
        end

    else  % k <= lag
        % part after 'else' not yet tested for lag>=3
        col0 = lag-(k-1);
        PDF(:,1:col0) = repmat(test(k,1:col0),Nsamples,1);
        for jj=col0+1:lag-1
            PDF(:,jj) = PDF(:,jj+1);
        end
    end

    pdf = gmx_sample(MU(k-1,:),SIG2(k-1,:),Nsamples);
    PDF(:,lag) = pdf;

    % simulate for all cases, again using the matrix version
    t0 = tic;
    test_ = repmat ( test(k,:), Nsamples, 1 );
    test_(:,1:lag) = PDF;
    [MU(k,:), SIG2(k,:), post] = gpx(hyp, inf, meanfunc, cov, lik, input, target, test_, post);
    calltime = toc(t0);
%     fprintf ( 'gpr_simul() ....... %f sec\n\n', calltime ); % uncomment
%     if you wish
    sum_time = sum_time + calltime;

    % approximate output distribution with gauss - calculate m and v
    mu(k,1) = mean(MU(k,:));
    s2(k,1) = mean(SIG2(k,:)) + mean((MU(k,:)-mu(k)).^2);

end

%     fprintf ( 'average time needed for gpr_simul() .... %f sec\n', sum_time/(num_iters-2) );
%     uncomment if you wish
return;

author: Dynamic-Systems-and-GP

repo: GPdyn


text:
function blend_test13 ( )

%*****************************************************************************80
%
%% BLEND_TEST13 tests BLEND_123.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2008
%
%  Author:
%
%    John Burkardt
%
  n1 = 3;
  n2 = 5;
  n3 = 4;

  d(1:n1,1:n2,1:n3) = 0.0;

  for i = 1 : n1
    for j = 1 : n2
      for k = 1 : n3
        if ( i == 1 | i == n1 | j == 1 | j == n2 | k == 1 | k == n3 )
          d(i,j,k) = ( i + j + k );
        end
      end
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BLEND_TEST13\n' );
  fprintf ( 1, '  BLEND_123 blends face values into a table.\n' );
  fprintf ( 1, '\n' );

  r8block_print ( n1, n2, n3, d, '  Initial data array' );

  for i = 1 : n1

    r = ( i - 1 ) / ( n1 - 1 );

    for j = 1 : n2

      s = ( j - 1 ) / ( n2 - 1 );

      for k = 1 : n3

        t = ( k - 1 ) / ( n3 - 1 );

        if ( i == 1 | i == n1 | j == 1 | j == n2 | k == 1 | k == n3 )
          continue
        end

        d(i,j,k) = blend_123 ( r, s, t, ...
          d(1,1,1), d(1,1,n3), d(1,n2,1), d(1,n2,n3), ...
          d(n1,1,1), d(n1,1,n3), d(n1,n2,1), d(n1,n2,n3), ...
          d(i,1,1), d(i,1,n3), d(i,n2,1), d(i,n2,n3), ...
          d(1,j,1), d(1,j,n3), d(n1,j,1), d(n1,j,n3), ...
          d(1,1,k), d(1,n2,k), d(n1,1,k), d(n1,n2,k), ...
          d(1,j,k), d(n1,j,k), d(i,1,k), d(i,n2,k), d(i,j,1), d(i,j,n3) )

      end

    end

  end

  r8block_print ( n1, n2, n3, d, '  Interpolated data array' );

  return
end


text:
function bestFits = ellipseDetection(img, params)
% ellipseDetection: Ellipse detection
%
% Overview:
% --------
% Fits an ellipse by examining all possible major axes (all pairs of points) and
% getting the minor axis using Hough transform. The algorithm complexity depends on
% the number of valid non-zero points, therefore it is beneficial to provide as many
% restrictions in the "params" input arguments as possible if there is any prior
% knowledge about the problem.
%
% The code is reasonably fast due to (optional) randomization and full code vectorization.
% However, as the algorithm needs to compute pairwise point distances, it can be quite memory
% intensive. If you get out of memory errors, either downsample the input image or somehow
% decrease the number of non-zero points in it.
% It can deal with big amount of noise but can have severe problem with occlusions (major axis
% end points need to be visible)
%
% Input arguments:
% --------
% img
%   - One-channel input image (greyscale or binary).
% params
%   - Parameters of the algorithm:
%       minMajorAxis: Minimal length of major axis accepted.
%       maxMajorAxis: Maximal length of major axis accepted.
%        rotation, rotationSpan: Specification of restriction on the angle of the major axis in degrees.
%                                If rotationSpan is in (0,90), only angles within [rotation-rotationSpan,
%                               rotation+rotationSpan] are accepted.
%       minAspectRatio: Minimal aspect ratio of an ellipse (in (0,1))
%       randomize: Subsampling of all possible point pairs. Instead of examining all N*N pairs, runs
%                  only on N*randomize pairs. If 0, randomization is turned off.
%       numBest: Top numBest to return
%       uniformWeights: Used to prefer some points over others. If false, accumulator points are weighted
%                       by their grey intensity in the image. If true, the input image is regarded as binary.
%       smoothStddev: In order to provide more stability of the solution, the accumulator is convolved with
%                     a gaussian kernel. This parameter specifies its standard deviation in pixels.
%
% Return value:
% --------
% Returns a matrix of best fits. Each row (there are params.numBest of them) contains six elements:
% [x0 y0 a b alpha score] being the center of the ellipse, its major and minor axis, its angle in degrees and score.
%
% Based on:
% --------
% - "A New Efficient Ellipse Detection Method" (Yonghong Xie Qiang , Qiang Ji / 2002)
% - random subsampling inspired by "Randomized Hough Transform for Ellipse Detection with Result Clustering"
%   (CA Basca, M Talos, R Brad / 2005)
%
% Update log:
% --------
% 1.1: More memory efficient code, better documentation, more parameters, more solutions possible, example code.
% 1.0: Initial version
%
%
% Author: Martin Simonovsky
% e-mail: <mys007@seznam.cz>
% Release: 1.1
% Release date: 25.7.2013
%
%
% --------
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

    % default values
    if nargin==1; params=[]; end
    % - parameters to contrain the search
    if ~isfield(params,'minMajorAxis');     params.minMajorAxis = 10; end
    if ~isfield(params,'maxMajorAxis');     params.maxMajorAxis = 200; end
    if ~isfield(params,'rotation');            params.rotation = 0; end
    if ~isfield(params,'rotationSpan');        params.rotationSpan = 0; end
    if ~isfield(params,'minAspectRatio');    params.minAspectRatio = 0.1; end
    if ~isfield(params,'randomize');        params.randomize = 2; end
    % - others
    if ~isfield(params,'numBest');            params.numBest = 3; end
    if ~isfield(params,'uniformWeights');   params.uniformWeights = true; end
    if ~isfield(params,'smoothStddev');        params.smoothStddev = 1; end

    eps = 0.0001;
    bestFits = zeros(params.numBest,6);
    params.rotationSpan = min(params.rotationSpan, 90);
    H = fspecial('gaussian', [params.smoothStddev*6 1], params.smoothStddev);

    [Y,X]=find(img);
    Y = single(Y); X = single(X);
    N = length(Y);

    fprintf('Possible major axes: %d * %d = %d\n', N, N, N*N);

    % compute pairwise distances between points (memory intensive!) and filter
    % TODO: do this block-wise, just appending the filtered results (I,J)
    distsSq = bsxfun(@minus,X,X').^2 + bsxfun(@minus,Y,Y').^2;
    [I,J] = find(distsSq>=params.minMajorAxis^2 & distsSq<=params.maxMajorAxis^2);
    idx = I<J;
    I = uint32(I(idx)); J = uint32(J(idx));

    fprintf('..after distance constraint: %d\n', length(I));

    % compute pairwise angles and filter
    if params.rotationSpan>0
        tangents = (Y(I)-Y(J)) ./ (X(I)-X(J));
        tanLo = tand(params.rotation-params.rotationSpan);
        tanHi = tand(params.rotation+params.rotationSpan);
        if tanLo<tanHi
            idx = tangents > tanLo & tangents < tanHi;
        else
            idx = tangents > tanLo | tangents < tanHi;
        end
        I = I(idx); J = J(idx);
        fprintf('..after angular constraint: %d\n', length(I));
    else
        fprintf('..angular constraint not used\n');
    end

    npairs = length(I);

    % compute random choice and filter
    if params.randomize>0
        perm = randperm(npairs);
        pairSubset = perm(1:min(npairs,N*params.randomize));
        clear perm;
        fprintf('..after randomization: %d\n', length(pairSubset));
    else
        pairSubset = 1:npairs;
    end

    % check out all hypotheses
    for p=pairSubset
        x1=X(I(p)); y1=Y(I(p));
        x2=X(J(p)); y2=Y(J(p));

        %compute center & major axis
        x0=(x1+x2)/2; y0=(y1+y2)/2;
        aSq = distsSq(I(p),J(p))/4;
        thirdPtDistsSq = (X-x0).^2 + (Y-y0).^2;
        K = thirdPtDistsSq <= aSq; % (otherwise the formulae in paper do not work)

        %get minor ax propositions for all other points
        fSq = (X(K)-x2).^2 + (Y(K)-y2).^2;
        cosTau = (aSq + thirdPtDistsSq(K) - fSq) ./ (2*sqrt(aSq*thirdPtDistsSq(K)));
        cosTau = min(1,max(-1,cosTau)); %inexact float arithmetic?!
        sinTauSq = 1 - cosTau.^2;
        b = sqrt( (aSq * thirdPtDistsSq(K) .* sinTauSq) ./ (aSq - thirdPtDistsSq(K) .* cosTau.^2 + eps) );

        %proper bins for b
        idxs = ceil(b+eps);

        if params.uniformWeights
            weights = 1;
        else
            weights = img(sub2ind(size(img),Y(K),X(K)));
        end
        accumulator = accumarray(idxs, weights, [params.maxMajorAxis 1]);

        %a bit of smoothing and finding the most busy bin
        accumulator = conv(accumulator,H,'same');
        accumulator(1:ceil(sqrt(aSq)*params.minAspectRatio)) = 0;
        [score, idx] = max(accumulator);

        %keeping only the params.numBest best hypothesis (no non-maxima suppresion)
        if (bestFits(end,end) < score)
            bestFits(end,:) = [x0 y0 sqrt(aSq) idx atand((y1-y2)/(x1-x2)) score];
            if params.numBest>1
                [~,si]=sort(bestFits(:,end),'descend');
                bestFits = bestFits(si,:);
            end
        end
    end
end

author: Sable


text:
function out = Csum(n,k, UBorLB, accuracy)
%sum of binomial coefficients from n choose 0 to n choose k.
%UBorLB = 'ub' or 'lb', upper or lower bound
%accuracy: 1 - more accurate, 0 - faster

%
%   Created in 2013 by Victoria Kostina (vkostina@caltech.edu)
%

if accuracy
    out = 0;
    for j = 0:k
        out = out + C(n, j, UBorLB);
    end
else
    mult = 1;
    if strcmpi(UBorLB, 'ub')
        if n <= 2*k
            mult = Inf;
        else
            mult = (n - k)/(n - 2*k);
        end
    end
    out = mult*C(n, k, UBorLB);
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre