text:
function y = vl_nnrelu(x,varargin)
%VL_NNRELU CNN rectified linear unit.
%   Y = VL_NNRELU(X) applies the rectified linear unit to the data
%   X. X can have arbitrary size.
%
%   DZDX = VL_NNRELU(X, DZDY) computes the derivative of the block
%   projected onto DZDY. DZDX and DZDY have the same dimensions as
%   X and Y respectively.
%
%   VL_NNRELU(...,'OPT',VALUE,...) takes the following options:
%
%   `Leak`:: 0
%      Set the leak factor, a non-negative number. Y is equal to X if
%      X is not smaller than zero; otherwise, Y is equal to X
%      multipied by the leak factor. By default, the leak factor is
%      zero; for values greater than that one obtains the leaky ReLU
%      unit.
%
%   ADVANCED USAGE
%
%   As a further optimization, in the backward computation it is
%   possible to replace X with Y, namely, if Y = VL_NNRELU(X), then
%   VL_NNRELU(X,DZDY) gives the same result as VL_NNRELU(Y,DZDY).
%   This is useful because it means that the buffer X does not need to
%   be remembered in the backward pass.

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

if ~isempty(varargin) && ~ischar(varargin{1})  % passed in dzdy
  dzdy = varargin{1} ;
  varargin(1) = [] ;
else
  dzdy = [] ;
end

opts.leak = 0 ;
opts = vl_argparse(opts, varargin, 'nonrecursive') ;

if opts.leak == 0
  if nargin <= 1 || isempty(dzdy)
    y = max(x, 0) ;
  else
    y = dzdy .* (x > 0) ;
  end
else
  if nargin <= 1 || isempty(dzdy)
    y = x .* (opts.leak + (1 - opts.leak) * (x > 0)) ;
  else
    y = dzdy .* (opts.leak + (1 - opts.leak) * (x > 0)) ;
  end
end

repo: LapSRN


text:
function [Prefix, ViewName, Ext] = AfniPrefix (Prefix)
%
%   [Prefix, View, Ext] = AfniPrefix (Name)
%
%Purpose:
%   returns the afni prefix, view and extension from a brick name
%   No checking is done to make sure that View or Ext are acceptable
%   see PrefixStatus instead
%
%Input Parameters:
%   Name : 'elvis' or 'elvis+orig' or 'elvis+orig.HEAD'
%
%
%Output Parameters:
%   Prefix: 'elvis'
%   View: '+orig'
%   Ext: .HEAD
%
%Key Terms:
%
%More Info :
%   PrefixStatus
%
%
%
%     Author : Ziad Saad
%     Date : Tue Dec 4 15:24:09 EST 2001
%     SSCC/NIMH/ National Institutes of Health, Bethesda Maryland

%Define the function name for easy referencing
FuncName = 'AfniPrefix';

[Prefix, Ext] = RemoveExtension(Prefix, '.HEAD|.BRIK');
[Prefix, ViewName] = RemoveExtension(Prefix, '+orig|+acpc|+tlrc');

err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
clear
close all
clc;

%time
dt=0.01;

%coeffs
Cg=1.3;
Cd=1;
minSpeed=20;

%green
greenheight=double(imread('testgreen.tiff'));
[gradX,gradY]=gradient(greenheight(:,:,1));
forceFun=@(R,RL) (RL-R)/dt*Cd + [interp2(-gradX,R(1,:),R(2,:));interp2(-gradY,R(1,:),R(2,:))]*Cg;
%forceFun=@(R,RL) [interp2(-gradX,R(1,:),R(2,:));interp2(-gradY,R(1,:),R(2,:))];
l=size(greenheight,1);

%graphics
%sph=sphere(30)
alt=70;
az=0;

% balls
startLoc=[64;400];

angles=linspace(0,2*pi,360/4+1);
angles=angles(1:end-1);
speeds=linspace(minSpeed,1000,20);

numParticles=length(speeds)*length(angles);

speeds=reshape(speeds,1,1,length(speeds));
angleVectors=[cos(angles);sin(angles)];
startConditions=bsxfun(@times,speeds,angleVectors);
startConditions=reshape(startConditions,2,length(speeds)*length(angles));

r=repmat(startLoc,1,numParticles);
rl=r-startConditions*dt;

i=0;
while sum(sum(r~=rl))
    i=i+1;
    rn=2*r-rl+(forceFun(r,rl))*dt^2;
    rl=r;
    r=rn;

    %static friction
    dr=rl-r;
    s=sqrt(dr(1,:).^2+dr(2,:).^2)/dt;

    %walls
    haltBallsEdge=(r>l)|(r<1);
    haltBalls=(haltBallsEdge(1,:)|haltBallsEdge(2,:))|(s<minSpeed);

    r(:,haltBalls)=rl(:,haltBalls);

    if mod(i,10)==0
        figure(1)
        surf(greenheight,'edgecolor','none')
        hold on
%         scatter3(r(1,1),r(2,1),greenheight(floor(r(2,1)),floor(r(1,1))),'wo','filled')
%         hold off
%         az=az+1;
%         view([az alt])
        imshow(greenheight,[min(min(greenheight)),max(max(greenheight))],'colormap',colormap('parula'))
        hold on
        scatter(r(1,:),r(2,:),'w.')
        hold off
        title(num2str(s(1)))
    end
end

author: BrianHaidet

repo: AlphaPhoenix


text:
function [ x, seed ] = sample_sphere_uniform ( m, n, seed )

%*****************************************************************************80
%
%% SAMPLE_SPHERE_UNIFORM samples points inside the unit sphere.
%
%  Discussion:
%
%    The sphere has center 0 and radius 1.
%
%    We first generate a point ON the sphere, and then distribute it
%    IN the sphere.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Russell Cheng,
%    Random Variate Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley, 1998, pages 168.
%
%    Reuven Rubinstein,
%    Monte Carlo Optimization, Simulation, and Sensitivity
%    of Queueing Networks,
%    Wiley, 1986, page 232.
%
%  Parameters:
%
%    Input, integer M, the dimension of the space.
%
%    Input, integer N, the number of points.
%
%    Input/output, integer SEED, a seed for the random number generator.
%
%    Output, real X(M,N), the points.
%
  exponent = 1.0 / m;

  for j = 1 : n
%
%  Fill a vector with normally distributed values.
%
    [ y, seed ] = r8ec_normal_01 ( m, seed );
%
%  Compute the length of the vector.
%
    norm = sqrt ( sum ( y(1:m).^2 ) );
%
%  Normalize the vector.
%
    y(1:m) = y(1:m) / norm;
%
%  Now compute a value to map the point ON the sphere INTO the sphere.
%
    [ r, seed ] = r8_uniform_01 ( seed );

    x(1:m,j) = r^exponent * y(1:m)';

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef TestMagnitude
    %TestMagnitude

    methods (Static)
        function test_1
            x = randn(5,6);
            y = randn(5,6);
            mag = cv.magnitude(x, y);
            validateattributes(mag, {class(x)}, {'size',size(x)});

            x = single(x);
            y = single(y);
            mag = cv.magnitude(x, y);
            validateattributes(mag, {class(x)}, {'size',size(x)});
        end

        function test_error_argnum
            try
                cv.magnitude();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function res = isign(x)
%ISIGN        Interval sign, internal function for slope abs
%

%sign of interval array:
%  +1      for x>=0
%  -1      for x<=0
%  [-1,1]  for 0 in x

% written  12/06/98     S.M. Rump
% modified 12/27/02     S.M. Rump  resinf,ressup changed to type double: fix of Matlab 6.5 discrepancy to 6.0
%

  xinf = inf(x);
  xsup = sup(x);
  resinf = double( xinf>=0 );
  resinf(xsup<=0) = -1;
  ressup = resinf;
  index = ( xinf<0 ) & ( xsup>0 );
  if any(index(:))
    resinf(index) = -1;
    ressup(index) = 1;
  end
  res = infsup(resinf,ressup);

repo: OpenMAS


text:
function [x,y,ind] = spm_MDP_VB_ERP(MDP,FACTOR,T)
% auxiliary routine for hierarchical electrophysiological responses
% FORMAT [x,y] = spm_MDP_VB_ERP(MDP,FACTOR,T)
%
% MDP    - structure (see spm_MDP_VB)
% FACTOR - hidden factors (at high and low level) to plot
% T      - flag to return cell of expectations (at time T; usually 1)
%
% x      - simulated ERPs (high-level) (full lines)
% y      - simulated ERPs (low level)  (dotted lines)
% ind    - indices or bins at the end of each (synchronised) epoch
%
% This routine combines first and second level hidden expectations by
% synchronising them; such that first level updating is followed by an
% epoch of second level updating - during which updating is suspended
% (and expectations are held constant). The ensuing spike rates can be
% regarded as showing delay period activity. In this routine, simulated
% local field potentials are band pass filtered spike rates (between eight
% and 32 Hz).
%
% Graphics are provided for first and second levels, in terms of simulated
% spike rates (posterior expectations), which are then combined to show
% simulated local field potentials for both levels (superimposed).
%
% At the lower level, only expectations about hidden states in the first
% epoch are returned (because the number of epochs can differ from trial
% to trial).
%
% see also: spm_MDP_VB_LFP (for single level belief updating)
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MDP_VB_ERP.m 7656 2019-08-26 14:00:36Z karl $

% defaults: assume the first factor is of interest
%==========================================================================
try, f1 = FACTOR(1); catch, f1 = 1; end
try, f2 = FACTOR(2); catch, f2 = 1; end

% and T = 1
%--------------------------------------------------------------------------
if nargin < 3, T = 1; end

for m = 1:numel(MDP)

    % dimensions
    %----------------------------------------------------------------------
    xn  = MDP(m).xn{f1};      % neuronal responses
    Nb  = size(xn,1);         % number of time bins per epochs
    Nx  = size(xn,2);         % number of states
    Ne  = size(xn,3);         % number of epochs

    % expected hidden states
    %======================================================================
    x     = cell(Ne,Nx);
    y     = cell(Ne);
    for k = 1:Ne
        for j = 1:Nx
            x{k,j} = xn(:,j,T,k);
        end
        if isfield(MDP,'mdp')
            y{k}   = spm_MDP_VB_ERP(MDP(m).mdp(k),f2,1);
        else
            y{k}   = [];
        end
    end

    if nargin > 2, return, end

    % synchronise responses
    %----------------------------------------------------------------------
    u   = {};
    v   = {};
    uu  = spm_cat(x(1,:));
    for k = 1:Ne

        % low-level
        %------------------------------------------------------------------
        v{end + 1,1} = spm_cat(y{k});
        if k > 1
            u{end + 1,1} = ones(size(v{end,:},1),1)*u{end,1}(end,:);
        else
            u{end + 1,1} = ones(size(v{end,:},1),1)*uu(1,:);
        end

        % time bin indices
        %------------------------------------------------------------------
        ind(k) = size(u{end},1);

        % high-level
        %------------------------------------------------------------------
        u{end + 1,1} = spm_cat(x(k,:));
        v{end + 1,1} = ones(size(u{end,:},1),1)*v{end,1}(end,:);

        % time bin indices
        %------------------------------------------------------------------
        ind(k) = ind(k) + size(u{end},1);

    end

    % accumulate over trials
    %----------------------------------------------------------------------
    U{m,1} = u;
    V{m,1} = v;

end

% time bin (seconds)
%--------------------------------------------------------------------------
u  = spm_cat(U);
v  = spm_cat(V);
dt = 1/64;
t  = (1:size(u,1))*dt;

% bandpass filter between 8 and 32 Hz
%--------------------------------------------------------------------------
c  = 1/32;
x  = log(u + c);
y  = log(v + c);
x  = spm_conv(x,2,0) - spm_conv(x,16,0);
y  = spm_conv(y,2,0) - spm_conv(y,16,0);

if nargout > 2, return, end

% simulated firing rates and the local field potentials
%==========================================================================

% higher-level unit responses
%--------------------------------------------------------------------------
factor = MDP(1).label.factor{f1};
name   = MDP(1).label.name{f1};

subplot(4,1,1), image(t,1:(size(u,2)),64*(1 - u')), ylabel('Unit')
title(sprintf('Unit reponses : %s',factor),'FontSize',16)
if numel(name) < 16
    grid on, set(gca,'YTick',1:numel(name))
    set(gca,'YTickLabel',name)
end

% lower-level unit responses
%--------------------------------------------------------------------------
factor = MDP(1).MDP(1).label.factor{f2};
name   = MDP(1).MDP(1).label.name{f2};

subplot(4,1,2), image(t,1:(size(v,2)),64*(1 - v')), ylabel('Unit')
title(sprintf('Unit reponses : %s',factor),'FontSize',16)
if numel(factor) < 16
    grid on, set(gca,'YTick',1:numel(name))
    set(gca,'YTickLabel',name)
end

% event related responses at both levels
%--------------------------------------------------------------------------
subplot(4,1,3), plot(t,x',t,y','-.')
title('Local field potentials','FontSize',16)
ylabel('Depolarisation'),spm_axis tight
grid on, xlabel('time (seconds)')

author: spm


text:
run ~/src/vlfeat/toolbox/vl_setup
run matconvnet/matlab/vl_setupnn

%model = '32s' ;
%model = '16s' ;
model = '8s' ;

blobs = load(sprintf('../testcaffe/blobs%s.mat', model)) ;
opts.modelPath = sprintf('matconvnet/data/models/pascal-fcn%s-dag.mat',model) ;

% load and fix model
net = dagnn.DagNN.loadobj(load(opts.modelPath)) ;
net.mode = 'test' ;
net.conserveMemory = false ;

filename = '/home/vedaldi/src/deep-seg/data/voc11/JPEGImages/2008_000073.jpg' ;
im = single(imread(filename)) ;
im = bsxfun(@minus, single(im), ...
            reshape(net.meta.normalization.averageImage,1,1,3)) ;

net.eval({'data', im}) ;

ok = ismember({net.vars.name}, fieldnames(blobs)) ;
names = {net.vars(ok).name} ;

for i=1:numel(names) ;
  name = names{i} ;
  matches = cellfun(@(x) ~isempty(regexp(x,['^' name 'x*$'])), {net.vars.name}) ;
  j = max(find(matches)) ;
  if isempty(j), continue ; end

  name_ = net.vars(j).name ;
  a=net.vars(j).value ;
  b=blobs.(name);
  b=permute(b, [3 4 2 1]) ;
  if i==1, b = b(:,:,[3 2 1]) ; end % BGR RGB

  del=  max(abs((a(:) - b(:)))) ;
  str=sprintf('%d %s vs %s = %g',i,name,name_,del);
  disp(str);
end

author: vlfeat

repo: matconvnet-fcn


text:
function tests = test_drawSphericalPolygon
% Test suite for the file drawSphericalPolygon.
%
%   Test suite for the file drawSphericalPolygon
%
%   Example
%   test_drawSphericalPolygon
%
%   See also
%     drawSphericalPolygon

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2022-09-20,    using Matlab 9.9.0.1570001 (R2020b) Update 4
% Copyright 2022 INRAE - BIA-BIBS.

tests = functiontests(localfunctions);

function test_Simple(testCase) %#ok<*DEFNU>
% Test call of function without argument.

p1 = [0 -1 0];
p2 = [0 0 1];
p3 = [-1 0 1];
p4 = [-1 -1 1];
p5 = [-1 -1 0];
poly = [p1 ; p2 ; p3;p4;p5];
hFig = figure; hold on;

h = drawSphericalPolygon([0 0 0 1], poly, 'LineWidth', 2);
assertTrue(testCase, all(ishandle(h)));

close(hFig);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function y = amexpo1s(N,t0,T);
%AMEXPO1S Generate one-sided exponential amplitude modulation.
%	Y = AMEXPO1S(N,T0,T) generates a one-sided exponential
%	amplitude modulation centered on a time T0, and with a
%	spread proportional to T.
%	This modulation is scaled such that Y(T0)=1.
%
%	N  : number of points.
%	T0 : arrival time of the exponential	(default : N/2).
%	T  : time spreading			(default : 2*sqrt(N)).
%	Y  : signal.
%
%	Examples:
%	 z=amexpo1s(160);plot(z);
%	 z=amexpo1s(160,20,40);plot(z);
%
%	See also AMEXPO2S, AMGAUSS, AMRECT, AMTRIANG.

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

if (nargin == 0),
 error ( 'The number of parameters must be at least 1.' );
elseif (nargin == 1),
 t0=N/2; T=2*sqrt(N);
elseif (nargin ==2),
 T=2*sqrt(N);
end;

if (N<=0),
 error('N must be greater or equal to 1.');
else
 tmt0=(1:N)'-t0;
 y = exp(-sqrt(pi)*tmt0/T).*(tmt0>=0.0);
end;

author: HeLiangHIT


text:
function [ net,res,opts ] = rnn_ff( net,opts )
%NET_FF Summary of this function goes here
%   Detailed explanation goes here

    n_frames=opts.parameters.n_frames;
    n_hidden_nodes=opts.parameters.n_hidden_nodes;
    batch_size=opts.parameters.batch_size;

    if opts.use_gpu
        opts.input_data=gpuArray(single(opts.input_data));
        if isfield(opts,'inputs_labels')
            opts.input_labels=gpuArray(single(opts.input_labels));
        end
        if isfield(opts,'input_predicts')
            opts.input_predicts=gpuArray(single(opts.input_predicts));
        end
    end
    if ~isfield(opts.parameters,'Id_w')
       opts.parameters.Id_w=1;
    end

    for f=1:n_frames
        if isfield(opts,'input_predicts')
            res.Fit{f}(1).predicts=opts.input_predicts(:,:,f);
        end
        if isfield(opts,'input_labels')
            res.Fit{f}(1).class=opts.input_labels(:,f);
        end
    end

    opts.loss=zeros(1,n_frames,'like',opts.input_data);

    if isfield(opts,'input_labels')
        opts.err=zeros(2,n_frames,'like',opts.input_data);
    end

    %%%%
    res.Hidden{1}=zeros(n_hidden_nodes,batch_size,'like',opts.input_data);

    for f=1:n_frames
        %Process inputs
        res.Input{f}(1).x=[res.Hidden{f};opts.input_data(:,:,f)];%inputs

        %Input transform
        [ net{1},res.Input{f},opts ] = net_ff( net{1},res.Input{f},opts );

        %Update hidden nodes;
        res.Hidden{f+1}=opts.parameters.Id_w.*res.Hidden{f} + res.Input{f}(end).x;

        %Data fitting transform
        res.Fit{f}(1).x=res.Hidden{f+1};

        [ net{2},res.Fit{f},opts ] = net_ff( net{2},res.Fit{f},opts );

    end

    %stats
    for f=1:n_frames
        if isfield(opts,'input_labels')
            opts.err(:,f)=error_multiclass(res.Fit{f}(1).class,res.Fit{f});
        end
        opts.loss(:,f)=mean(res.Fit{f}(end).x(:));
    end

    if isfield(opts,'input_labels')
        opts.err=mean(opts.err,2)./opts.parameters.batch_size;
    end
    opts.loss=mean(opts.loss(:));
end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function partition=new_partitions(sigma,r,x)
%sigma=1;
%r=8.4;
T=r/2/sigma;
N=ceil(r);
%width=sigma*(2*N/r-1);
%x=-width:1/200:width;

partition=zeros(N,length(x));

for k=1:N
  if k==1
    for j=1:length(x)
      if x(j)<=sigma-N/T
        partition(k,j)=0;
      elseif sigma-N/T<x(j) & x(j)<-sigma
        partition(k,j)=rho((-x(j)-sigma)/(N/T-2*sigma));
      elseif -sigma<=x(j) & x(j)<=(1-N)/T+sigma
        partition(k,j)=1;
      elseif (1-N)/T+sigma<x(j) & x(j)<1/T-sigma
	partition(k,j)=rho((x(j)-((1-N)/T+sigma))/(N/T-2*sigma));
%      elseif (1-N)/T+sigma<x(j) & x(j)<(1-N/2)/T
%        partition(k,j)=1-rho((-x(j)+(1-N/2)/T)/(N/T/2-sigma))/2;
%      elseif x(j)==(1-N/2)/T
%        partition(k,j)=1/2;
%      elseif (1-N/2)/T<x(j) & x(j)<1/T-sigma
%        partition(k,j)=rho((x(j)-(1-N/2)/T)/(N/T/2-sigma))/2;
      else
        partition(k,j)=0;
      end
    end
  elseif k==N
    for j=1:length(x)
      if x(j)<=sigma-1/T
        partition(k,j)=0;
      elseif sigma-1/T<x(j) & x(j)<(N-1)/T-sigma
        partition(k,j)=1-rho((x(j)-(sigma-1/T))/(N/T-2*sigma));
%      elseif sigma-1/T<x(j) & x(j)<(N/2-1)/T
%        partition(k,j)=rho((-x(j)+(N/2-1)/T)/(N/T/2-sigma))/2;
%      elseif x(j)==(N/2-1)/T
%        partition(k,j)=1/2;
%      elseif (N/2-1)/T<x(j) & x(j)<(N-1)/T-sigma
%        partition(k,j)=1-rho((x(j)-(N/2-1)/T)/(N/T/2-sigma))/2;
      elseif (N-1)/T-sigma<=x(j) & x(j)<=sigma
        partition(k,j)=1;
      elseif sigma<x(j) & x(j)<N/T-sigma
        partition(k,j)=rho((x(j)-sigma)/(N/T-2*sigma));
      else
        partition(k,j)=0;
      end
    end
  else
    for j=1:length(x)
      if x(j)<=(-N+k-1)/T+sigma
        partition(k,j)=0;
      elseif (-N+k-1)/T+sigma<x(j) & x(j)<(k-1)/T-sigma
        partition(k,j)=1-rho((x(j)-((-N+k-1)/T+sigma))/(N/T-2*sigma));
%      elseif (-N+k-1)/T+sigma<x(j) & x(j)<(k-1-N/2)/T
%        partition(k,j)=rho((-x(j)+(k-1-N/2)/T)/(N/T/2-sigma))/2;
%      elseif x(j)==(k-1-N/2)/T
%        partition(k,j)=1/2;
%      elseif (k-1-N/2)/T<x(j) & x(j)<(k-1)/T-sigma
%        partition(k,j)=1-rho((x(j)-(k-1-N/2)/T)/(N/T/2-sigma))/2;
      elseif (k-1)/T-sigma<=x(j) & x(j)<=(-N+k)/T+sigma
        partition(k,j)=1;
      elseif (-N+k)/T+sigma<x(j) & x(j)<k/T-sigma
        partition(k,j)=rho((x(j)-((-N+k)/T+sigma))/(N/T-2*sigma));
%      elseif (-N+k)/T+sigma<x(j) & x(j)<(k-N/2)/T
%        partition(k,j)=1-rho((-x(j)+(k-N/2)/T)/(N/T/2-sigma))/2;
%      elseif x(j)==(k-N/2)/T
%        partition(k,j)=1/2;
%      elseif (k-N/2)/T<x(j) & x(j)<k/T-sigma
%        partition(k,j)=rho((x(j)-(k-N/2)/T)/(N/T/2-sigma))/2;
      else
        partition(k,j)=0;
      end
    end
  end
%  plot(x,partition(k,:))
%  axis([min(x) max(x) 0 1])
%  [k N]
%  pause
end
%plot(x,partition'); %'
%pause

%total=zeros(size(x));
%for j=1:N
%total=total+partition(j,:);
%end
%spy(1-total);

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function [c,E,exitVal]=minimaxPolyFit(f,xSpan,n,RelTol,maxIter)
%MINIMAXPOLYFIT Find an interpolating polynomial of order n to minimize
%               the maximum error approximating a scalar function f in the
%               range [xSpan(1);xSpan(2)]. The Remez algorithm is used.
%
%INPUTS: f A handle to the function that is to be approximated. The use
%          y=f(x) take a scalar x and returns a value y.
%    xSpan The range of inputs of x over which a minimax interpolating
%          polynomial is desired. The format is [xMin;xMax]. If this is
%          omitted or an empty matrix is passed, the values [-1;1] are
%          used.
%        n The degree of the polynomial approximation desired. If this is
%          omitted or an empty matrix is passed, the default value of n=5
%          is used. If too small a value is used and the function
%          oscillates a lot, the output polynomial will not be a true
%          minimax polynomial, because the algorithm assumes that there is
%          one maxima between regions chosen for selecting points. Also,
%          if f is an mth order polynomial, but is symmetric, then n
%          should be m+1 due to assumptions in the algorithm (the
%          resulting output will have the extra coefficient zero).
%   RelTol The relative tolerance of the maximum error used for
%          determining convergence. Iterations adjust the maximum error
%          observed with respect to points at which the function is
%          evaluated. When the maximum error stops changing across
%          iterations, convergence is assumed. The relative maximum error
%          is abs(E-EPrev)/abs(E); where E is the error and EPRev is the
%          error from the previous iteration. If this parameter is omitted
%          or an empty matrix is passed, a value of 1e-12 is used.
%  maxIter The maximum number of iterations for convergence. If omitted, a
%          value of 50 is used.
%
%OUTPUTS: c An interpolating polynomial such that polyval(c,xRel) gives the
%           interpolated value of f where xRel is the fractional distance
%           between xSpan(1) and xSpan(2). The value xRel ranges from 0 to
%           1. The development of the interpolation across 0 and 1 is often
%           numerically better than explicitly making the polynomial for
%           values from xSpan(1) to xSpan(2).
%         E The maximum error of the interpolating polynomial.
%   exitVal A value indicating how the algorithm terminated. Possible
%           values are:
%           0: The tolerance goal was met.
%           1: the maximum number of iterations was met.
%
%The algorithm used to fit a polynomial to an arbitrary nonlinear function
%is Remez algorithm. This is described ine Remez algorithm is well known
%and is described in Chapter 6 of [1].
%
%To initalize the Remez algorithm, one must choose a set of points in the
%region to be interpolated. One might be tempted to unifromely space the
%points, but as noted in [2] and Chaper 8.3 of [3], that can be a very bad
%idea. Rather, the use of Chebychev nodes is better and has much better
%asymptotic properties as the number of nodes increases.
%
%Given a set of points in xSpan, the algorithm find an interpolating
%polynomial that mades the magnitude of the error at all of the points the
%same and the sign of the error alternate. This is because there is a
%theorem that such an alternating equierror solution is a minimax solution.
%However, only a minimax solution for the initial set of points is given.
%Thus, the algorithm tries to move the points to maximize the error.
%
%Since the sign of the error between any two points flips, there is at
%least one zero between the points. This algorithm assumes that there is
%exactly one zero between the points, which should be reasonable if n is
%chosen to be large enough. In such an instance, the true error maxima of
%the function are between the zeros (and the outer bounds) of the function.
%Thus, the points can be moved there. Basically by maximizing the error at
%the chosen points, and getting an interpolating polynomial with equal
%error at the points, one finds the minimax approximation.
%
%As an example, consider interpolating the exponential function between 0
%and 2 with a sixth-order polynomial:
% f=@(x)exp(x);
% xSpan=[0;2];
% n=6;
% [c,E]=minimaxPolyFit(f,xSpan,n)
% %The maximum error E should be about 8.7e-6.
% %The resulting interpolating polynomial can be plotted as a fraction of
% %the distance between the points in xSpan. We will plot the error.
% numPoints=500;
% x=linspace(xSpan(1),xSpan(2),numPoints);
% xFrac=(x-xSpan(1))/(xSpan(2)-xSpan(1));
% error=polyval(c,xFrac)-exp(x);
% plot(x,error)
%
%REFERENCES:
%[1] V. K. Dzyadyk and I. A. Shevchuk, Theory of Uniform Approximation
%    of Functions by Polynomials, 1st ed. Berlin: Walter de Gruyter, 2008.
%[2] R. B. Platte, L. N. Trefethen, and A. B. Kuijlaars, "Impossibility of
%    fast stable approximation of analytic functions from equispaced
%    samples," SIAM Review, vol. 53, no. 2, pp. 308-318, 2011.
%[3] R. L. Burden and J. D. Faires, Numerical Analysis, 9th ed. Australia:
%    Brooks/ Cole, 2011.
%
%July 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(maxIter))
   maxIter=50;
end

if(nargin<4||isempty(RelTol))
   RelTol=1e-12;
end

if(nargin<3||isempty(n))
    n=5;
end

if(nargin<2||isempty(xSpan))
    xSpan=[-1;1];
end

%Step 1: Choose an initial set of nodes for interpolation.
%n+2 nodes are needed for an order-n polynomial approximation. The initial
%guess is set to the roots of the Chebyshev polynomial of the first kind.
%Such a choice is standard in most implementations of the Remez algorithm,
%because the Chebyshev polynomials distribute more nodes to the outside of
%the interval then to the center, reducing Runge's phenomenon.
k=(n+2):-1:1;
%The Chebyshev nodes for the (-1,1) span.
xNodes=cos((2*k(:)-1)/(2*(n+2))*pi);
%We are reducing the span of the nodes from -1/2 to 1/2 so that the
%interpolating polynomial on the output of the function takes an input from
%0 to 1.
xNodes=xNodes/2;

curIter=1;
EPrev=Inf;%This is used to determine convergence of the error estimate E.
while(1)
    %Step 2: Solve for the maximum error E and the coefficients of the
    %        interpolating polynomial P given the current set of nodes.
    %        This involves solving a linear system of equations.

    %Allocate space for the matrix.
    hMat=zeros(n+2,n+2);

    xDelta=(xNodes+1/2);
    hMat(:,1)=1;%The zeroth-order terms.
    for curCol=2:(n+1)
        hMat(:,curCol)=hMat(:,curCol-1).*xDelta;%This is xDelta.^(curCol-1)
    end

    %The last column corresponds to the E error term. This last column is
    %why the matrix is not a Vandermonde matrix.
    hMat(:,n+2)=(-1).^(0:(n+1));

    %Map the (-1,1) Chebyshev nodes into the range given by xSpan.
    xFNodes=(xSpan(2)-xSpan(1))*(xNodes+1/2)+xSpan(1);

    %Solve the system for the polynomial coefficients c and the error E
    coeffs=hMat\f(xFNodes);

    %The flipud is to make the ordering the same as that used by the
    %polyval function.
    c=flipud(coeffs(1:(n+1)));
    E=coeffs(n+2);

    relTolVal=abs(E-EPrev)/abs(E);

    %Check for convergence. Convergence is determined is the E value stops
    %changing significantly.
    %The check for the relative tolerance value being a nan deals with the
    %case where E=0 and E-EPrev=0.
    if(isnan(relTolVal)||relTolVal<RelTol)
        E=abs(E);
        exitVal=0;
        break;
    end

    %Return if the maximum number of iterations has passed.
    if(curIter==maxIter)
        E=abs(E);
        exitVal=1;
        break;
    end

    EPrev=E;

    %The coefficients c form the order-n minimax approximation of the
    %function at the n nodes. The error at those nodes has magnitude E.
    %However, these are most likely not the maximum error points for the
    %full function f. We want to move the points xNodes so that we are
    %minimizing the maximum error with respect to f in the interval.

    %Step 3: The Exchange Step. Given the current polynomial interpolation,
    %move the points to where it maximizes the magnitude of the error
    %function.

    %Step 3.1: Find the zeros of the error function between the control
    %          points.
    costFunZero=@(x)errorFun(x,xSpan,c,f);
    %There are n+1 points; the extra 2 points allocated are to make finding
    %the extrema in the next part of the step simpler.
    xIntervals=zeros(n+3,1);
    for curPoint=1:(n+1)
       %Find the zero of the error function between points curPoint and
       %curPoint+1.
       %Precision is arbitrary, but since X is from -1/2 to 1/2, it makes
       %sense to hard-code in the tolerance on X (and use a higher
       %tolernace than is the default)
       options=optimset('TolX',1e-14);
       xIntervals(curPoint+1)=fzero(costFunZero,[xNodes(curPoint),xNodes(curPoint+1)],options);
    end

    %Step 3.2: Find the extreme points and make them the new nodes.
    %The extrema are bracketed between each pair of roots. There are also
    %two more extrema located between the start of the interval xSpan(1)
    %and the first root and between the last root and the end of the
    %interval xSpan(2). We want to find all of the extrema. We will do line
    %searches to find those points. Those points are placed into xNodes as
    %they will form the new nodes for the next iteration.
    xIntervals(1)=-1/2;
    xIntervals(end)=1/2;

    costFunMax=@(x)negAbsErrorFun(x,xSpan,c,f);
    for curPoint=1:(n+2)
        %Find the maximum of the error function in the bounded interval
        %using a line search algorithm.
        options=optimset('TolX',1e-14);
        xNodes(curPoint)=fminbnd(costFunMax,xIntervals(curPoint),xIntervals(curPoint+1),options);
    end

    %Step 4: Iterate.
    curIter=curIter+1;
end
end

function val=errorFun(x,xSpan,c,f)
   %Map the (-1/2,1/2) Chebyshev node into the range given by xSpan.
   xf=(xSpan(2)-xSpan(1))*(x+1/2)+xSpan(1);

   val=polyval(c,x+1/2)-f(xf);
end

function val=negAbsErrorFun(x,xSpan,c,f)
   %Map the (-1/2,1/2) Chebyshev node into the range given by xSpan.
   xf=(xSpan(2)-xSpan(1))*(x+1/2)+xSpan(1);

   val=-abs(polyval(c,x+1/2)-f(xf));
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [rfAmp,rfPhase,rfFreq,rfCoil,rfTime]=rfSLR(p)
%create a SLR rf pulse starting from tStart
%  Inputs for dzrf are:
%    np -- number of points.         (required)
%    tb -- time-bandwidth product    (required)
%    ptype -- pulse type.  Options are:
%      st  -- small tip angle         (default)
%      ex  -- pi/2 excitation pulse
%      se  -- pi spin-echo pulse
%      sat -- pi/2 saturation pulse
%      inv -- inversion pulse
%    ftype -- filter design method.  Options are:
%      ms  -- Hamming windowed sinc (an msinc)
%      pm  -- Parks-McClellan equal ripple
%      ls  -- Least Squares           (default)
%      min -- Minimum phase (factored pm)
%      max -- Maximum phase (reversed min)
%    d1 -- Passband ripple        (default = 0.01)
%    d2 -- Stopband ripple        (default = 0.01)
%    pclsfrac -- pcls tolerance   (default = 1.5)

tStart=p.tStart;
tEnd=p.tEnd;
dt=p.dt;
FA=p.FA;
TBP=p.TBP;
rfFreq=p.rfFreq;
rfPhase=p.rfPhase;
PRipple=p.PRipple;
SRipple=p.SRipple;
PulseType=p.PulseType;
FilterType=p.FilterType;
rfCoil=p.CoilID;
Duplicates=max(1,p.Duplicates);
DupSpacing=max(0,p.DupSpacing);

rfTime=linspace(tStart,tEnd,ceil((tEnd-tStart)/dt)+1);
rf = dzrf(ceil((tEnd-tStart)/dt)+1,TBP,PulseType,FilterType,PRipple,SRipple,1.5);
rfAmp = real(rf);
rfAmp = DoB1Scaling(rfAmp,dt,FA)*rfAmp; %B1 Scaling

rfPhase=(rfPhase)*ones(size(rfTime));
rfFreq=(rfFreq)*ones(size(rfTime));
rfCoil=(rfCoil)*ones(size(rfTime));
rfAmp(1)=0;
rfAmp(end)=0;
rfFreq(1)=0;
rfFreq(end)=0;
rfPhase(1)=0;
rfPhase(end)=0;

% Create Duplicates
if Duplicates~=1 & DupSpacing ~=0
    rfAmp=repmat(rfAmp,[1 Duplicates]);
    rfFreq=repmat(rfFreq,[1 Duplicates]);
    rfPhase=repmat(rfPhase,[1 Duplicates]);
    rfCoil=repmat(rfCoil,[1 Duplicates]);
    TimeOffset = repmat(0:DupSpacing:(Duplicates-1)*DupSpacing,[length(rfTime) 1]);
    rfTime=repmat(rfTime,[1 Duplicates]) + (TimeOffset(:))';
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function coords = gsp_lle(G, dim, param)
%GSP_LLE Local Linear Embedding
%   Usage: coords = gsp_lle(G, dim);
%          coords = gsp_lle(G, dim, param);
%
%   Input parameters
%         G          : Graph
%         dim        : Dimensionality of the embedding
%         param      : Structure of optional parameters
%
%   Output parameters
%         coords     : Coordinates of the embedding
%
%   This function uses the weight matrix of a graph G, in order to compute
%   a *dim* -dimensional embedding (output coordinates). The algorithm used
%   is Locally Linear Embedding (LLE). Warning, this function might not
%   work if the graph is not connected.
%
%   *param* is a structure with optional parameters:
%
%   * *param.tol*    : Tolerance for the spectral gap (default 1e-6).
%   * *param.kernel* : The kernel used to create the graph weight matrix:
%     + 'exp'        : exponential kernel ($e^(frac{-d_{ij}}{sigma^2})$)
%     + '1/x'        : inverse of x kernel ($frac{1}{sigma+d_{ij}}$)
%     + '1/x^2'      : inverse of x^2 kernel ($frac{1}{(sigma+d_{ij})^2}$)
%     + 'resistance' : Resistance distance.
%   * *param.k*      : Max number of nearest neighbors. If not defined, the
%
%   number of neighbors varies from node to node since the algorithm
%   considers all the columns $j$ of the weight matrix that have non
%   zero values on the line $i$ as the neighbors of $i$.
%
%   References: saul2000introduction
%
%   See also: gsp_update_coordinates gsp_isomap gsp_laplacian_eigenmaps
%
%   Demo: gsp_demo_graph_embedding

% Authors : Dion O. E. Tzamarias
% Date    : 20/11/2015

%% TODO Fix bug param.k
%%

if nargin<3
    param = struct;
end

if ~isfield(param,'tol'), param.tol = 1e-6; end
if ~isfield(param,'kernel'), param.kernel = 'exp'; end

% Handling default parameters

W = zeros(G.N,G.N);

% if ~isfield(param,'k')
    [D] = gsp_weight2distance(G);
    [rows, cols ] = find(D);
    idx = cell(cols(end),1);
    for jj=1:cols(end)
        idx{jj} = rows(cols==jj).';
    end
% else
%     [D,idx] = gsp_weight2distance(G, param.kernel, param.k);
% end
D = D.^2; % D contains the squared distances
for ii=1:length(D)

    %   C = 0.5*(repmat(sum(D(idx{ii},idx{ii}),2),1,length(idx{ii}))/length(idx{ii})...
    %       + repmat(sum(D(idx{ii},idx{ii}),2).',length(idx{ii}),1)/length(idx{ii}) ...
    %       - D(idx{ii},idx{ii}) - repmat(sum(sum(D(idx{ii},idx{ii}),1)),...
    %       length(idx{ii}),length(idx{ii}))/length(idx{ii}).^2);
    C = 0.5*(bsxfun(@plus,sum(D(idx{ii},idx{ii}),2),sum(D(...
        idx{ii},idx{ii}),1))/length(idx{ii})-bsxfun(@plus,...
        D(idx{ii},idx{ii}),sum(sum(D(idx{ii},idx{ii}),1))...
        /length(idx{ii}).^2));

    if abs(det(C))< param.tol  %   C is singular
        C = C + (param.tol*trace(C)) / length(idx{ii})...
            * eye(size(C));
    end
%     if abs(det(C))< param.tol*1e-3
%         warning('C is singular')
%     end

    W(ii,idx{ii}) = C\ones(length(idx{ii}),1);
    W(ii,:) = W(ii,:)./sum(W(ii,:));
end

M = (eye(G.N,G.N) - W)' * (eye(G.N,G.N) - W);
M = sparse(M);

% Compute first dim eigenvectors of M
options.disp = 0;
options.isreal = 1;
options.issym = 1;
% only need bottom (no_dims + 1) eigenvectors
tol = param.tol;
[coords, e] = eigs(M, dim + 1, tol, options);
e = diag(e);

if nnz(e<param.tol) > 1
    disp('Multiple zero eigenvalues')
end

[~ , ind] = sort(e, 'ascend');
coords = coords(:,ind);
coords = coords(:,2:end);
end

%   The input parameter type refers to the type of the problem:
%   * 'weight_matrix' : Solve problem using the weight matrix of the graph G
%     Use the pairwise distance equation in order to reconstruct data from
%     neighbours.
%   * 'coords'       : Slove the problem using the coordinates. This is the
%     typical LLE aproach.

% switch lower(type)
%     case 'coords'
%         for ii=1:G.N
%
%             k_near_indx = knnsearch(G.coords,G.coords(ii,:),'k',n_neighbor+1);
%             k_nearest = G.coords(k_near_indx(2:end),:);
%
%             C = (repmat(G.coords(ii,:),n_neighbor,1) - k_nearest)* ...
%                 (repmat(G.coords(ii,:),n_neighbor,1) - k_nearest)' ;
%
%             if abs(det(C))<tol  % C is not singular
%                 C = C + ((tol*trace(C))^2) / n_neighbor * eye(size(C));
%             end
%             if abs(det(C))<tol
%                 warning('C is not singular')
%             end
%
%             W(ii,k_near_indx(2:end)) = C\ones(n_neighbor,1);
%             W(ii,:) = W(ii,:)./sum(W(ii,:));
%         end
%
%     case 'weight_matrix'

repo: gspbox


text:
%imu model parameter should be defined in cont. time
function [STM Qd]=mdl_ned_dcm(pos_n, vel_n, Cbn, acc, Aimu, Qimu, Cimu, Rimu, dt)

%continious model
[Anav N]=sys_ned_dcm_v000(pos_n, vel_n, Cbn, acc, 0, []);
nst_imu=size(Aimu,1);
nst=9+nst_imu;
Ac=[Anav N*Cimu;zeros(nst_imu,9) Aimu];
Qc=[N*Rimu*N' zeros(9,nst_imu);zeros(nst_imu,9) Qimu];

%Discretize the model
mx_a = dt*[-Ac,Qc;zeros(nst),Ac'];
mx_b = expm(mx_a);
STM = mx_b(nst+1:2*nst,nst+1:2*nst)';
Qd = STM*mx_b(1:nst,nst+1:2*nst);

repo: ekfmonoslam


text:
function [data, XYZvoxSphere, XYZmmSphere] = iimg_sphere_timeseries(images, XYZmmCenter, radius)
% :Usage:
% ::
%
%     function [data, XYZvoxSphere, XYZmmSphere] = iimg_sphere_timeseries(images, XYZmm, radius)
%
% :Inputs:
%
%   **images:**
%        list of image files
%
%   **XYZmm:**
%        [3 x n] array of mm coords
%
%   **radiu:**
%        radius in mm of sphere to generate
%
% :Outputs:
%
%   **data:**
%        voxel data
%
%   **XYZvoxSphere:**
%        voxel

    Vimages = spm_vol(images);
    [XYZvox(1,:) XYZvox(2,:) XYZvox(3,:)] = ind2sub(Vimages(1).dim(1:3), 1:prod(Vimages(1).dim(1:3)));
    XYZmm = Vimages(1).mat(1:3, :)*[XYZvox; ones(1, size(XYZvox, 2))];

    dist = [XYZmm(1,:) - XYZmmCenter(1); XYZmm(2,:) - XYZmmCenter(2); XYZmm(3,:) - XYZmmCenter(3)];
    whVoxelsInSphere = find(sum(dist.^2) <= radius^2);

    XYZmmSphere = unique(XYZmm(:,whVoxelsInSphere)', 'rows')';
    XYZvoxSphere = unique(XYZvox(:,whVoxelsInSphere)', 'rows')';

    data = spm_get_data(Vimages, XYZvoxSphere);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [axis] = extractHeatmapAxisData(obj,axis_data,axisName)
%extract information related to each axis
%   axis_data is the data extrated from the figure, axisName take the
%   values 'x' 'y' or 'z'

%-------------------------------------------------------------------------%

%-axis-side-%
% axis.side = eval(['axis_data.' axisName 'AxisLocation;']);

%-------------------------------------------------------------------------%

%-axis zeroline-%
axis.zeroline = false;

%-------------------------------------------------------------------------%

%-axis autorange-%
axis.autorange = false;

%-------------------------------------------------------------------------%

%-axis exponent format-%
axis.exponentformat = obj.PlotlyDefaults.ExponentFormat;

%-------------------------------------------------------------------------%

%-axis tick font size-%
axis.tickfont.size = axis_data.FontSize;

%-------------------------------------------------------------------------%

%-axis tick font family-%
axis.tickfont.family = matlab2plotlyfont(axis_data.FontName);

%-------------------------------------------------------------------------%

tl = eval(['axis_data.' axisName 'Data;']);
tl = length(tl);

w = axis_data.Position(4);
h = axis_data.Position(3);

ticklength = min(obj.PlotlyDefaults.MaxTickLength,...
    max(tl*w*obj.layout.width,tl*h*obj.layout.height));

%-axis ticklen-%
axis.ticklen = 0.1; %ticklength;

%-------------------------------------------------------------------------%

% col = eval(['255*axis_data.' axisName 'Color;']);
axiscol = 'rgb(150, 150, 150)';

%-axis linecolor-%
axis.linecolor = axiscol;
%-axis tickcolor-%
axis.tickcolor = axiscol;
%-axis tickfont-%
axis.tickfont.color = 'black';
%-axis grid color-%
axis.gridcolor = 'rgb(0, 0, 0)';

%-------------------------------------------------------------------------%

axis.showgrid = true;

%-------------------------------------------------------------------------%

lw = 0.5;
linewidth = max(1,lw*obj.PlotlyDefaults.AxisLineIncreaseFactor);

%-axis line width-%
axis.linewidth = linewidth;
%-axis tick width-%
axis.tickwidth = linewidth;
%-axis grid width-%
axis.gridwidth = linewidth*1.2;

%-------------------------------------------------------------------------%

%-setting ticks-%
axis.ticks = 'inside';
axis.mirror = true;

labels = eval(['axis_data.' axisName 'DisplayLabels;']);
vals = eval(['axis_data.' axisName 'DisplayData;']);

axis.showticklabels = true;
axis.type = 'category';
axis.autorange = true;
axis.ticktext = labels;
axis.tickvals = vals;
axis.autotick = false;
axis.tickson = 'boundaries';

%-------------------------------LABELS------------------------------------%

label = eval(['axis_data.' axisName 'Label;']);

%-------------------------------------------------------------------------%

%-title-%
axis.title = label;

%-------------------------------------------------------------------------%

%-axis title font color-%
axis.titlefont.color = 'black';

%-------------------------------------------------------------------------%

%-axis title font size-%
axis.titlefont.size = axis_data.FontSize*1.3;
axis.tickfont.size = axis_data.FontSize*1.15;

%-------------------------------------------------------------------------%

%-axis title font family-%
axis.titlefont.family = matlab2plotlyfont(axis_data.FontName);
axis.tickfont.family = matlab2plotlyfont(axis_data.FontName);

%-------------------------------------------------------------------------%

if strcmp(axis_data.Visible,'on')
    %-axis showline-%
    axis.showline = true;
else
    %-axis showline-%
    axis.showline = false;
    %-axis showticklabels-%
    axis.showticklabels = false;
    %-axis ticks-%
    axis.ticks = '';
end

%-------------------------------------------------------------------------%
end

author: plotly


text:
function area = circle_segment_area_from_angle ( r, theta )

%*****************************************************************************80
%
%% CIRCLE_SEGMENT_AREA_FROM_ANGLE computes the area of a circle segment.
%
%  Discussion:
%
%    Begin with a circle of radius R.  Choose two points P1 and P2 on the
%    circle, and draw the chord P1:P2.  This chord divides the circle
%    into two pieces, each of which is called a circle segment.
%    Consider one of the pieces.  The "angle" of this segment is the angle
%    P1:C:P2, where C is the center of the circle.  Let Q be the point on
%    the chord P1:P2 which is closest to C.  The "height" of the segment
%    is the distance from Q to the perimeter of the circle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%    0 < R.
%
%    Input, real THETA, the angle of the circle segment.
%
%    Output, real AREA, the area of the circle segment.
%
  area = r * r * ( theta - sin ( theta ) ) / 2.0;

  return
end


text:
function printMetricsExt(metrics, metricsInfo, dispHeader,dispMetrics,padChar)
% print metrics
%
% ...
%

% default names
if nargin==1
    metricsInfo.names.long = {'IDF1', 'IDP', 'IDR',    'Recall','Precision','False Alarm Rate', ...
        'GT Tracks','Mostly Tracked','Partially Tracked','Mostly Lost', ...
        'False Positives', 'False Negatives', 'ID Switches', 'Fragmentations', ...
        'MOTA','MOTP', 'MOTA Log'};

    metricsInfo.names.short = {'IDF1', 'IDP', 'IDR',    'Rcll','Prcn','FAR', ...
        'GT','MT','PT','ML', ...
        'FP', 'FN', 'IDs', 'FM', ...
        'MOTA','MOTP', 'MOTAL'};

    metricsInfo.widths.long = [5 4 4   6 9 16 9 14 17 11 15 15 11 14 5 5 8];
    metricsInfo.widths.short = [5 4 4  5 5 5 4 4 4 4 6 6 5 5 5 5 5];

    metricsInfo.format.long = {'.1f','.1f','.1f', ...
        '.1f','.1f','.2f', ...
        'i','i','i','i', ...
        'i','i','i','i', ...
        '.1f','.1f','.1f'};

    metricsInfo.format.short=metricsInfo.format.long;
end

namesToDisplay=metricsInfo.names.long;
widthsToDisplay=metricsInfo.widths.long;
formatToDisplay=metricsInfo.format.long;

namesToDisplay=metricsInfo.names.short;
widthsToDisplay=metricsInfo.widths.short;
formatToDisplay=metricsInfo.format.short;

if nargin<3, dispHeader=1; end
if nargin<4
    dispMetrics=1:length(metrics);
end
if nargin<5
    padChar={' ',' ','|',' ',' ','| ','',' ',' ','|','',' ',' ','| ',' ',' ',' ',' '};
end

if dispHeader
    for m=dispMetrics
        printString=sprintf('fprintf(''%%%is%s'',char(namesToDisplay(m)))',widthsToDisplay(m),char(padChar(m)));
        eval(printString)
    end
    fprintf('\n');
end

for m=dispMetrics
    printString=sprintf('fprintf(''%%%i%s%s'',metrics(m))',widthsToDisplay(m),char(formatToDisplay(m)),char(padChar(m)));
    eval(printString)
end

% if standard, new line
if nargin<4
    fprintf('\n');
end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function d = degree(varargin)

F = varargin{1};
F = flatten(F);
d = subdegree(F.clauses);

function d = subdegree(clauses)
d = -inf;
for i = 1:length(clauses)
    d = max(d,degree(sdpvar(clauses{i}.data)));
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [axes,mult] = elements(cs,multiplicity)
% extract symmetry elements by multiplicity
%
% Syntax
%
%   axes = elements(cs,multiplicity) % rotational axes with fixed multiplicity
%   [axes,multiplicity] = elements(cs) % all rotational axes with multiplicity
%
% Input
%  cs - @crystalSymmetry
%  multiplicity - double
%
% Output
%  axes - rotational axes vector3d
%  multiplicity - double
%

rot = cs.rot(cs.rot.angle>1*degree);
axes =  rot.axis;
mult = round(2*pi ./ rot.angle);
[axes, ~, id] = unique(axes,'tolerance',1e-3,'antipodal');
mult = accumarray(id,mult,[],@max);

axes = [axes;-axes];
mult = [mult;mult];

if nargin == 2
  axes = axes(mult == multiplicity);
end

author: mtex-toolbox

repo: mtex

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


text:
% Copyright (C) 2018  Symeon Symeonidis, Stefanos Tsantilas, Stelios Mitilineos
% simos421@gmail.com, steftsantilas@gmail.com, smitil@gmail.com
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.

function CstPickMidPoint(mws,Name,id)

Pick = invoke(mws,'Pick');
invoke(Pick,'PickMidpointFromId',strcat('component1:',Name),num2str(id));
end

repo: CST-MATLAB-API


text:
%DEMO_MODELASSESMENT1 Demonstration for model assessment with WAIC,
% DIC, number of effective parameters and ten-fold cross validation
%
%  Description
%    We will consider the regression problem in demo_regression1.
%    The analysis is conducted with full Gaussian process, and FIC
%    and PIC sparse approximations. The performance of these models
%    are compared by evaluating ten-fold cross validation,
%    leave-one-out cross-validation, WAIC, DIC and the effective
%    number of parameters. The inference will be conducted using
%    maximum a posterior (MAP) estimate for the parameters, via
%    full Markov chain Monte Carlo (MCMC) and with an integration
%    approximation (IA) for the parameters.
%
%    This demo is organised in three parts:
%     1) data analysis with full GP model
%     2) data analysis with FIC approximation
%     3) data analysis with PIC approximation
%
%  See also DEMO_REGRESSION1, DEMO_REGRESSION_SPARSE1

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

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

%========================================================
% PART 1 data analysis with full GP model
%========================================================
disp('Full GP model with Gaussian noise model')

% Load the data
S = which('demo_regression1');
L = strrep(S,'demo_regression1.m','demodata/dat.1');
data=load(L);
x = [data(:,1) data(:,2)];
y = data(:,3);
[n, nin] = size(x);

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

% ---------------------------
% --- Construct the model ---
gpcf = gpcf_sexp('lengthScale', [1 1], 'magnSigma2', 0.2^2);
lik = lik_gaussian('sigma2', 0.2^2);
gp = gp_set('lik', lik, 'cf', gpcf, 'jitterSigma2', 1e-9);

% -----------------------------
% --- Conduct the inference ---
%
% We will make the inference first by finding a maximum a posterior estimate
% for the parameters via gradient based optimization. After this we will
% perform an extensive Markov chain Monte Carlo sampling for the parameters.
%
disp('MAP estimate for the parameters')

% --- MAP estimate ---
%     (see gp_optim for more details)

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

% Evaluate the effective number of parameters, DIC and WAIC with focus on
% latent variables.
models{1} = 'full_MAP';
p_eff_latent = gp_peff(gp, x, y);
% For easier comparison to other methods, compute mean log
% predictive density (mlpd) instead of deviance (-2n*mlpd)
[DIC_latent, p_eff_latent2] = gp_dic(gp, x, y, 'focus', 'latent', 'output', 'mlpd');
WAICV(1) = gp_waic(gp,x,y);
WAICG(1) = gp_waic(gp,x,y, 'method', 'G');

% Evaluate the 10-fold cross-validation results.
disp('MAP estimate for the parameters - k-fold-CV')
cvres =  gp_kfcv(gp, x, y, 'display', 'fold');
mlpd_cv(1) = cvres.mlpd_cv;
rmse_cv(1) = cvres.mrmse_cv;

disp('MAP estimate for the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(gp, x, y);
mlpd_loo(1) = mean(lpy);
rmse_loo(1) = sqrt(mean((y-Ef).^2));

% --- MCMC approach ---
disp('MCMC integration over the parameters')

% Do the sampling (this takes about 1 minute)
[rfull,g,opt] = gp_mc(gp, x, y, 'nsamples', 220, 'display', 20);

% After sampling delete the burn-in and thin the sample chain
rfull = thin(rfull, 21, 2);

% Evaluate the effective number of parameters, DIC and WAIC.
models{2} = 'full_MCMC';
% For easier comparison to other methods, compute mean log
% predictive density (mlpd) instead of deviance (-2n*mlpd)
[DIC(2), p_eff(2)] =  gp_dic(rfull, x, y, 'focus', 'param', 'output', 'mlpd');
[DIC2(2), p_eff2(2)] =  gp_dic(rfull, x, y, 'focus', 'all', 'output', 'mlpd');
WAICV(2) = gp_waic(rfull,x,y);
WAICG(2) = gp_waic(rfull,x,y, 'method', 'G');

% Evaluate the 10-fold cross validation results.
%
% We reduce the number of samples so that the sampling takes less time.
% 50 is too small sample size, though, and for reliable results the 10-CV
% should be run with larger sample size.
disp('MCMC integration over the parameters - k-fold-CV')
opt.nsamples= 50;
cvres =  gp_kfcv(gp, x, y, 'inf_method', 'MCMC', 'opt', opt, 'rstream', 1, 'display', 'fold');
mlpd_cv(2) = cvres.mlpd_cv;
rmse_cv(2) = cvres.mrmse_cv;

disp('MCMC integration over the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(rfull, x, y);
mlpd_loo(2) = mean(lpy);
rmse_loo(2) = sqrt(mean((y-Ef).^2));

% --- Integration approximation approach ---
disp('Grid integration over the parameters')

gp_array = gp_ia(gp, x, y, 'int_method', 'grid');

models{3} = 'full_IA';
% For easier comparison to other methods, compute mean log
% predictive density (mlpd) instead of deviance (-2n*mlpd)
[DIC(3), p_eff(3)] =  gp_dic(gp_array, x, y, 'focus', 'param', 'output', 'mlpd');
[DIC2(3), p_eff2(3)] =  gp_dic(gp_array, x, y, 'focus', 'all', 'output', 'mlpd');
WAICV(3) = gp_waic(gp_array,x,y);
WAICG(3) = gp_waic(gp_array,x,y, 'method', 'G');

% Then the 10 fold cross-validation.
disp('Grid integration over the parameters - k-fold-CV')
clear opt
opt.int_method = 'grid';
cvres = gp_kfcv(gp, x, y, 'inf_method', 'IA', 'opt', opt, 'display', 'fold');
mlpd_cv(3) = cvres.mlpd_cv;
rmse_cv(3) = cvres.mrmse_cv;

disp('Grid integration over the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(gp_array, x, y);
mlpd_loo(3) = mean(lpy);
rmse_loo(3) = sqrt(mean((y-Ef).^2));

%========================================================
% PART 2 data analysis with FIC GP
%========================================================
disp('GP with FIC sparse approximation')

% ---------------------------
% --- Construct the model ---

% Here we conduct the same analysis as in part 1, but this time we
% use FIC approximation

% Initialize the inducing inputs in a regular grid over the input space
[u1,u2]=meshgrid(linspace(-1.8,1.8,6),linspace(-1.8,1.8,6));
X_u = [u1(:) u2(:)];

% Create the FIC GP structure
gp_fic = gp_set('type', 'FIC', 'lik', lik, 'cf', gpcf, 'jitterSigma2', 1e-6, 'X_u', X_u)

% -----------------------------
% --- Conduct the inference ---

% --- MAP estimate using scaled conjugate gradient algorithm ---
disp('MAP estimate for the parameters')

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

% Evaluate the effective number of parameters and DIC with focus on
% latent variables.
models{4} = 'FIC_MAP';
p_eff_latent(4) = gp_peff(gp_fic, x, y);
[DIC_latent(4), p_eff_latent2(4)] = gp_dic(gp_fic, x, y, 'focus', 'latent', 'output', 'mlpd');
WAICV(4) = gp_waic(gp_fic,x,y);
WAICG(4) = gp_waic(gp_fic,x,y, 'method', 'G');

% Evaluate the 10-fold cross validation results.
disp('MAP estimate for the parameters - k-fold-CV')
cvres = gp_kfcv(gp_fic, x, y, 'display', 'fold');
mlpd_cv(4) = cvres.mlpd_cv;
rmse_cv(4) = cvres.mrmse_cv;

disp('MAP estimate for the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(gp_fic, x, y);
mlpd_loo(4) = mean(lpy);
rmse_loo(4) = sqrt(mean((y-Ef).^2));

% --- MCMC approach ---
% (the inducing inputs are fixed)
disp('MCMC integration over the parameters')

% Do the sampling (this takes about 1 minute)
rfic = gp_mc(gp_fic, x, y, 'nsamples', 220, 'display', 20);

% After sampling we delete the burn-in and thin the sample chain
rfic = thin(rfic, 21, 2);

% Evaluate the effective number of parameters, DIC and WAIC. Note that
% the effective number of parameters as a second output, but here
% we use explicitly the gp_peff function
models{5} = 'FIC_MCMC';
[DIC(5), p_eff(5)] =  gp_dic(rfic, x, y, 'focus', 'param', 'output', 'mlpd');
[DIC2(5), p_eff2(5)] =  gp_dic(rfic, x, y, 'focus', 'all', 'output', 'mlpd');
WAICV(5) = gp_waic(rfic,x,y);
WAICG(5) = gp_waic(rfic,x,y, 'method', 'G');

% We reduce the number of samples so that the sampling takes less time.
% 50 is too small sample size, though, and for reliable results the 10-CV
% should be run with larger sample size. We also set the save option to 0.
clear opt
opt.nsamples= 50; opt.display=20;
disp('MCMC integration over the parameters - k-fold-CV')
cvres = gp_kfcv(gp_fic, x, y, 'inf_method', 'MCMC', 'opt', opt, 'display', 'fold');
mlpd_cv(5) = cvres.mlpd_cv;
rmse_cv(5) = cvres.mrmse_cv;

disp('MCMC integration over the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(rfic, x, y);
mlpd_loo(5) = mean(lpy);
rmse_loo(5) = sqrt(mean((y-Ef).^2));

% --- Integration approximation approach ---
disp('Grid integration over the parameters')
gpfic_array = gp_ia(gp_fic, x, y, 'int_method', 'grid');

models{6} = 'FIC_IA';
[DIC(6), p_eff(6)] =  gp_dic(gpfic_array, x, y, 'focus', 'param', 'output', 'mlpd');
[DIC2(6), p_eff2(6)] =  gp_dic(gpfic_array, x, y, 'focus', 'all', 'output', 'mlpd');
WAICV(6) = gp_waic(gpfic_array,x,y);
WAICG(6) = gp_waic(gpfic_array,x,y, 'method', 'G');

% Then the 10 fold cross-validation.
disp('Grid integration over the parameters - k-fold-CV')
clear opt
opt.int_method = 'grid';
cvres = gp_kfcv(gp_fic, x, y, 'inf_method', 'IA', 'opt', opt, 'display', 'fold');
mlpd_cv(6) = cvres.mlpd_cv;
rmse_cv(6) = cvres.mrmse_cv;

disp('Grid integration over the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(gpfic_array, x, y);
mlpd_loo(6) = mean(lpy);
rmse_loo(6) = sqrt(mean((y-Ef).^2));

%========================================================
% PART 3 data analysis with PIC approximation
%========================================================
disp('GP with PIC sparse approximation')

[u1,u2]=meshgrid(linspace(-1.8,1.8,6),linspace(-1.8,1.8,6));
X_u = [u1(:) u2(:)];

% Initialize test points
[p1,p2]=meshgrid(-1.8:0.1:1.8,-1.8:0.1:1.8);
p=[p1(:) p2(:)];

% set the data points into clusters. Here we construct two cell arrays.
%  trindex  contains the block index vectors for training data. That is
%           x(trindex{i},:) and y(trindex{i},:) belong to the i'th block.
b1 = [-1.7 -0.8 0.1 1 1.9];
mask = zeros(size(x,1),size(x,1));
trindex={};
for i1=1:4
  for i2=1:4
    ind = 1:size(x,1);
    ind = ind(: , b1(i1)<=x(ind',1) & x(ind',1) < b1(i1+1));
    ind = ind(: , b1(i2)<=x(ind',2) & x(ind',2) < b1(i2+1));
    trindex{4*(i1-1)+i2} = ind';
  end
end

% Create the PIC GP structure and set the inducing inputs and block indexes
gpcf = gpcf_sexp('lengthScale', [1 1], 'magnSigma2', 0.2^2);
lik = lik_gaussian('sigma2', 0.2^2);

gp_pic = gp_set('type', 'PIC', 'lik', lik, 'cf', gpcf, 'jitterSigma2', 1e-6, 'X_u', X_u);
gp_pic = gp_set(gp_pic, 'tr_index', trindex);

% -----------------------------
% --- Conduct the inference ---

% --- MAP estimate using scaled conjugate gradient algorithm ---
disp('MAP estimate for the parameters')

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

models{7} = 'PIC_MAP';
p_eff_latent(7) = gp_peff(gp_pic, x, y);
[DIC_latent(7), p_eff_latent2(7)] = gp_dic(gp_pic, x, y, 'focus', 'latent', 'output', 'mlpd');
WAICV(7) = gp_waic(gp_pic, x, y);
WAICG(7) = gp_waic(gp_pic, x, y, 'method', 'G');

% Evaluate the 10-fold cross validation results.
disp('MAP estimate for the parameters - k-fold-CV')
cvres = gp_kfcv(gp_pic, x, y, 'display', 'fold');
mlpd_cv(7) = cvres.mlpd_cv;
rmse_cv(7) = cvres.mrmse_cv;

disp('MAP estimate for the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(gp_pic, x, y);
mlpd_loo(7) = mean(lpy);
rmse_loo(7) = sqrt(mean((y-Ef).^2));

% --- MCMC approach ---
disp('MCMC integration over the parameters')

% Do the sampling (this takes about 1 minute)
rpic = gp_mc(gp_pic, x, y, 'nsamples', 220, 'display', 20);

% After sampling we delete the burn-in and thin the sample chain
rpic = rmfield(rpic, 'tr_index');
rpic = thin(rpic, 21, 2);
rpic.tr_index = trindex;

% Evaluate the effective number of parameters and DIC. Note that
% the effective number of parameters as a second output, but here
% we use explicitly the gp_peff function
models{8} = 'PIC_MCMC';
[DIC(8), p_eff(8)] =  gp_dic(rpic, x, y, 'focus', 'param', 'output', 'mlpd');
[DIC2(8), p_eff2(8)] =  gp_dic(rpic, x, y, 'focus', 'all', 'output', 'mlpd');
WAICV(8) = gp_waic(rpic, x, y);
WAICG(8) = gp_waic(rpic, x, y, 'method', 'G');

% We reduce the number of samples so that the sampling takes less time.
% 50 is too small sample size, though, and for reliable results the 10-CV
% should be run with larger sample size. We also set the save option to 0.
clear opt
opt.nsamples= 50; opt.display=20;
disp('MCMC integration over the parameters - k-fold-CV')
cvres = gp_kfcv(gp_pic, x, y, 'inf_method', 'MCMC', 'opt', opt, 'display', 'fold');
mlpd_cv(8) = cvres.mlpd_cv;
rmse_cv(8) = cvres.mrmse_cv;

disp('MCMC integration over the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(rpic, x, y);
mlpd_loo(8) = mean(lpy);
rmse_loo(8) = sqrt(mean((y-Ef).^2));

% --- Integration approximation approach ---
disp('Grid integration over the parameters')

gppic_array = gp_ia(gp_pic, x, y, 'int_method', 'grid');

models{9} = 'PIC_IA';
[DIC(9), p_eff(9)] =  gp_dic(gppic_array, x, y, 'focus', 'param', 'output', 'mlpd');
[DIC2(9), p_eff2(9)] =  gp_dic(gppic_array, x, y, 'focus', 'all', 'output', 'mlpd');
WAICV(9) = gp_waic(gppic_array, x, y);
WAICG(9) = gp_waic(gppic_array, x, y, 'method', 'G');

% Then the 10 fold cross-validation.
disp('Grid integration over the parameters - k-fold-CV')
clear opt
opt.int_method = 'grid';
cvres = gp_kfcv(gp_pic, x, y, 'inf_method', 'IA', 'opt', opt, 'display', 'fold');
mlpd_cv(9) = cvres.mlpd_cv;
rmse_cv(9) = cvres.mrmse_cv;

disp('Grid integration over the parameters - LOO-CV')
[Ef,Varf,lpy] =  gp_loopred(gppic_array, x, y);
mlpd_loo(9) = mean(lpy);
rmse_loo(9) = sqrt(mean((y-Ef).^2));

%========================================================
% PART 4 Print the results
%========================================================
disp('Summary of the results')
S = '      ';
for i = 1:length(models)
    S = [S '  ' models{i}];
end

S = sprintf([S '\n CV-mlpd  %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], mlpd_cv);
S = sprintf([S '\n CV-rmse  %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], rmse_cv);
S = sprintf([S '\n LOO-mlpd %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], mlpd_loo);
S = sprintf([S '\n LOO-rmse %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], rmse_loo);
S = sprintf([S '\n ']);
S = sprintf([S '\n WAIC_V   %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], WAICV);
S = sprintf([S '\n WAIC_G   %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], WAICG);
S = sprintf([S '\n ']);
S = sprintf([S '\n DIC_h    %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], DIC);
S = sprintf([S '\n DIC_a    %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], DIC2);
S = sprintf([S '\n DIC_l    %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], DIC_latent);
S = sprintf([S '\n peff_h   %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], p_eff);
S = sprintf([S '\n peff_a   %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], p_eff2);
S = sprintf([S '\n peff_l   %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], p_eff_latent);
S = sprintf([S '\n peff_l2  %5.2f     %5.2f     %5.2f     %5.2f    %5.2f     %5.2f    %5.2f    %5.2f    %5.2f'], p_eff_latent2);
S = sprintf([S '\n ']);
S = sprintf([S '\n ']);
S = sprintf([S '\n The notation is as follows:']);
S = sprintf([S '\n CV-mlpd  = mean log predictive density from the 10-fold CV. ']);
S = sprintf([S '\n CV-rmse  = root mean squared error from the 10-fold LOO-CV. ']);
S = sprintf([S '\n LOO-mlpd = mean log predictive density from the LOO-CV. ']);
S = sprintf([S '\n LOO-rmse = root mean squared error from the 10-fold CV. ']);
S = sprintf([S '\n WAIC_V   = WAIC via variance method ']);
S = sprintf([S '\n WAIC_G   = WAIC via Gibbs training utility method ']);
S = sprintf([S '\n DIC_h    = DIC with focus on parameters. ']);
S = sprintf([S '\n DIC_a    = DIC with focus on parameters and latent variables (all). ']);
S = sprintf([S '\n DIC_l    = DIC with focus on latent variables. ']);
S = sprintf([S '\n peff_h   = effective number of parameters (latent variables marginalized). ']);
S = sprintf([S '\n peff_a   = effective number of parameters and latent variables. ']);
S = sprintf([S '\n peff_l   = effective number of latent variables evaluated with gp_peff. ']);
S = sprintf([S '\n peff_l2  = effective number of latent variables evaluated with gp_dic. ']);
S = sprintf([S '\n '])

author: fieldtrip

repo: fieldtrip


text:
function i4_bit_hi1_test ( )

%*****************************************************************************80
%
%% I4_BIT_HI1_TEST tests I4_BIT_HI1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_BIT_HI1_TEST\n' );
  fprintf ( 1, '  I4_BIT_HI1 returns the location of the high 1 bit.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     I     I4_BIT_HI1(I)\n' );
  fprintf ( 1, '\n' );

  seed = 123456789;

  for test = 1 : 10

    [ i, seed ] = i4_uniform_ab ( 0, 100, seed );

    j = i4_bit_hi1 ( i );

    fprintf ( 1, '%6d %6d\n', i, j );

  end

  return
end

repo: jburkardt-m


text:
function new_slide = pptTitle(op, ttlText, fontSize, fontName);
%
% new_slide = pptTitle(op, ttlText, [fontSize=32], [fontName]);
%
% Add a title slide to a PowerPoint presentation, with the provided
% text. The optional fontSize and fontName fields allow you to specify
% the type and size of the text font. The default fontSize is 32, the
% default fontName is whatever is saved in the power point file's template.
%
% SEE ALSO: pptOpen, pptPaste, pptClose.
%
% ras, 07/2007.
if ~exist('fontSize', 'var'),	fontSize = 32;		end
if notDefined('fontName'),	fontName = '';		end

slide_count = get(op.Slides, 'Count');

% Add a the title slide:
slide_count = int32(double(slide_count)+1);
new_slide = invoke(op.Slides, 'Add', slide_count, 1);

% Set the title text, font
set(new_slide.Shapes.Title.TextFrame.TextRange, 'Text', ttlText);

if ~isempty(fontSize)
	set(new_slide.Shapes.Title.TextFrame.TextRange.Font, 'Size', fontSize);
end

if ~isempty(fontName)
	set(new_slide.Shapes.Title.TextFrame.TextRange.Font, 'Name', fontName);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [g, nu] = cmpndNoiseNuG(noise, mu, varSigma, y)

% CMPNDNOISENUG  Update nu and g parameters associated with compound noise model.

% NOISE

% NOISE

g = zeros(size(mu));
nu = zeros(size(g));
for i = 1:length(noise.comp)
  [g(:, i), nu(:, i)] = ...
      noiseUpdateNuG(noise.comp{i}, ...
                       mu(:, i), varSigma(:, i), ...
                       y(:, i));
  if any(nu(:, i)< 0)
    if noise.comp{i}.logconcave
      warning('nu less than zero in log concave model.')
    else
      fprintf('nu less than zero\n')
    end
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


x = [ 1 2 3 4 5 6];
y = [ 2 6 8 7 8 5];
barh(x,y);
title('\bfExample of a Horizontal Bar Plot');
xlabel('\bf\ity');
ylabel('\bf\itx');
axis([0 10 0 7]);


text:
function [mesh,elt2fce] = mshFace(mesh)
%+========================================================================+
%|                                                                        |
%|                 OPENMSH - LIBRARY FOR MESH MANAGEMENT                  |
%|           openMsh is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : mshFace.m                                     |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Submesh to face mesh                          |
%|  `---'  |                                                              |
%+========================================================================+

% Particles mesh
if (size(mesh.elt,2) == 1)
    error('mshFace.m : unavailable case')

% Edge mesh
elseif (size(mesh.elt,2) == 2)
    error('mshFace.m : unavailable case')

% Triangular mesh
elseif (size(mesh.elt,2) == 3)
    fce2vtx = mesh.elt;
    elt2fce = (1:size(mesh.elt,1))';
    col     = mesh.col;

% Tetrahedral mesh
elseif (size(mesh.elt,2) == 4)
    % All faces
    fce2vtx = [ mesh.elt(:,[2,3,4]) ; mesh.elt(:,[3,4,1]) ; ...
        mesh.elt(:,[4,1,2]) ; mesh.elt(:,[1,2,3]) ];

    % Faces unicity
    tmp           = sort(fce2vtx,2);
    [~,I,elt2fce] = unique(tmp,'rows');

    % Final faces
    fce2vtx = fce2vtx(I,:);
    elt2fce = reshape(elt2fce,size(mesh.elt,1),4);

    % Colours
    col             = zeros(size(fce2vtx,1),1);
    tmp             = mesh.col * ones(1,4);
    col(elt2fce(:)) = tmp(:);

    % Identify boundary
    bnd      = ( accumarray(elt2fce(:),1).*col ~= accumarray(elt2fce(:),tmp(:)) );
    col(bnd) = -1;

% Unknown type
else
    error('mshFace.m : unavailable case')
end

% Mesh format
mesh = msh(mesh.vtx,fce2vtx,col);
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function laguerre_polynomial_test ( )

%*****************************************************************************80
%
%% LAGUERRE_POLYNOMIAL_TEST tests the LAGUERRE_POLYNOMIAL library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGUERRE_POLYNOMIAL_TEST:\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the LAGUERRE_POLYNOMIAL library.\n' );

  laguerre_polynomial_test01 ( );
  laguerre_polynomial_test02 ( );
  laguerre_polynomial_test03 ( );
  laguerre_polynomial_test04 ( );
  laguerre_polynomial_test05 ( );
  laguerre_polynomial_test06 ( );

  p = 5;
  b = 0.0;
  laguerre_polynomial_test07 ( p, b );

  p = 5;
  b = 1.0;
  laguerre_polynomial_test07 ( p, b );

  p = 5;
  e = 0;
  laguerre_polynomial_test08 ( p, e );

  p = 5;
  e = 1;
  laguerre_polynomial_test08 ( p, e );
%
%  Make some plots.
%
  a = 0.0;
  b = 5.0;
  index = [ 0, 1, 2, 3, 4, 5, 10 ];
  filename = 'l_polynomial.png';
  l_polynomial_plot ( a, b, index, filename );

  a = 0.0;
  b = 5.0;
  index = [ 0, 1, 2, 3, 4, 5, 10 ];
  index2 = [ 1, 1, 1, 1, 1, 1, 1 ];
  filename = 'lm1_polynomial.png';
  lm_polynomial_plot ( a, b, index, index2, filename );

  a = 0.0;
  b = 5.0;
  index = [ 0, 1, 2, 3, 4, 5, 10 ];
  index2 = [ 2, 2, 2, 2, 2, 2, 2 ];
  filename = 'lm2_polynomial.png';
  lm_polynomial_plot ( a, b, index, index2, filename );

  a = 0.0;
  b = 5.0;
  index =  [ 0,   1,   2,   3,   4,   5,   10 ];
  index2 = [ 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5 ];
  filename = 'lf05_function.png';
  lf_function_plot ( a, b, index, index2, filename );

  close
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGUERRE_POLYNOMIAL_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer


text:
%% Total Variation Denoising
% Test for Rudin-Osher-Fatemi denoising (ROF) using FB-like method.

addpath('../');
addpath('../toolbox/');

%%
% Load image.

n = 256;
y = load_image('lena',n*2);
y = rescale(crop(y,n));
y = y + randn(n)*.06;

%%
% Display it.

clf;
imageplot(clamp(y));

%%
% We aim at minimising:

%%
% |min_x 1/2*norm(y-x,'fro')^2 + lambda*norm(K(x),1)|

%%
% Regularization parameter.

lambda = .2;

%%
% where |K| is a vectorial gradient and |norm(u,1)| is a vectorial L1
% norme.

K = @(x)grad(x);
KS = @(x)-div(x);

%%
% It can be put as the minimization of |F(K*x) + G(x)|

Amplitude = @(u)sqrt(sum(u.^2,3));
F = @(u)lambda*sum(sum(Amplitude(u)));
G = @(x)1/2*norm(y-x,'fro')^2;

%%
% The proximity operator of |F| is the vectorial soft thresholding.

Normalize = @(u)u./repmat( max(Amplitude(u),1e-10), [1 1 2] );
ProxF = @(u,tau)repmat( perform_soft_thresholding(Amplitude(u),lambda*tau), [1 1 2]).*Normalize(u);
ProxFS = compute_dual_prox(ProxF);

%%
% The proximity operator of G.

ProxG = @(x,tau)(x+tau*y)/(1+tau);

%%
% Function to record progression of the functional.

options.report = @(x)G(x) + F(K(x));

%%
% Run the ADMM algorihtm.

options.niter = 300;
[xAdmm,EAdmm] = perform_admm(y, K,  KS, ProxFS, ProxG, options);

%%
% Display image.

clf;
imageplot(xAdmm);

%%
% Since the functional to mimize is stricly convex, we can use a FB scheme
% on the dual problem.

GradGS = @(x)x+y;
L = 8;
options.method = 'fista';
[xFista,EFista] = perform_fb_strongly(y, K, KS, GradGS, ProxFS, L, options);

%%
% Compare the energy decays.

clf;
plot([EAdmm(:) EFista(:)]);
axis tight;
legend('ADMM', 'FISTA');
axis([1 length(EAdmm) EFista(end)*.9 2000]);

author: gpeyre

repo: matlab-toolboxes


text:
function searchlight_saveresults(modeldir)
% This function combines and saves searchlight analysis results.
%
% :Usage:
% ::
%
%     searchlight_saveresults(modeldir)
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2014  Wani Woo
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..
%
% :Input:
%
%   **modeldir:**
%        the directory where the searchlight result mat files (e.g.,
%        searchlight_results_*_01Aug2014.mat)
%
% :Outputs:
%
%   This function saves the brain maps where each voxel contains a summary
%   stat value (e.g., accuracy or outcome correlation) in the modeldir. The
%   naming convention of the result maps are as follow:
%
%   results_searchlight_(a)_(b)_dat(c).nii
%      1. will be a number - the number of test results
%      2. acc: accuracy, r: outcome correlation, p: p-values, se: standard error,
%         and thr: threshold for the single-interval test
%      3. If the test was on one dataset, (c) will be empty, but if the
%         test was cross-prediction, (c) will be 11, 12, 21, or 22.
%           - 11: trained on the first dataset, and tested on the first dataset
%           - 12: trained on the first dataset, and tested on the second dataset
%           - 21: trained on the second dataset, and tested on the first dataset
%           - 22: trained on the second dataset, and tested on the second dataset
%
% :Output Examples:
% ::
%
%    results_searchlight_1_r_dat11.nii (outcome correlation)
%    results_searchlight_1_p_dat11.nii (p value for the correlation values)
%
%    results_searchlight_2_acc_dat11.nii (accuracy)
%    results_searchlight_2_p_dat11.nii (p value for the accuracy)
%    results_searchlight_2_se_dat11.nii (standard error for the accuracy)
%    results_searchlight_2_thr_dat11.nii (threshold for the accuracy test)

load(fullfile(modeldir, 'searchlight_dream_variables.mat'));
data = load(fullfile(modeldir, 'searchlight_dream_variables.mat'), 'varargin');

do_cross = false;

% parsing varargin to know it's cross_classify or not
for i = 1:numel(data.varargin)
    if ischar(data.varargin{i})
        switch data.varargin{i}
            case 'dat2'
                do_cross = true;
        end
    end
end

res_files = filenames(fullfile(modeldir, 'searchlight_results_*mat'));

% load one test_results to set up results variables
load(res_files{1}, 'test_results');

% image frame
if do_cross
    dat.dat = zeros(size(dat.dat,1),4);
else
    dat.dat = zeros(size(dat.dat,1),1);
end

% results fields
for m = 1:numel(test_results)
    res_fields{m} = fields(test_results{m});
    for j = 1:numel(res_fields{m}), eval(['dat_' num2str(m) '_'  res_fields{m}{j} ' = dat;']); end
end

% combining results
for ii = 1:numel(res_files)
    load(res_files{ii}, 'test_results');

    for m = 1:numel(test_results)
        for j = 1:numel(res_fields{m})
            eval(['test_results{' num2str(m) '}.' res_fields{m}{j} '(isnan(test_results{' num2str(m) '}.' res_fields{m}{j} ')) = 0;']);
            eval(['dat_' num2str(m) '_' res_fields{m}{j} '.dat = dat_' num2str(m) '_' res_fields{m}{j} '.dat + test_results{' num2str(m) '}.' res_fields{m}{j} ';']);
        end
    end
end

% image names for each column
if do_cross
    col_names = {'dat11', 'dat12', 'dat21', 'dat22'};
else
    col_names = {'dat'};
end

% write images
for m = 1:numel(test_results)
    for jj = 1:numel(res_fields{m})
        eval(['for ii = 1:size(dat_' num2str(m) '_' res_fields{m}{jj} '.dat,2), dat.dat = dat_' ...
            num2str(m) '_' res_fields{m}{jj} '.dat(:,ii); dat.fullpath = fullfile(modeldir, [''results_searchlight_' ...
            num2str(m) '_' res_fields{m}{jj} '_'' col_names{ii} ''.nii'']); write(dat); end']);
    end
end

end


text:
function drawRawLines(SenFig,Raw)

% DRAWRAWLINES  Draw raw lines.
%   DRAWRAWLINES(SENFIG,RAW)  redraws the 2D segments in RAW on figure
%   SenFig.

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

if ~isempty(Raw.data.segments.coord)

    m = size(Raw.data.segments.coord,2); % number of segments
    pn = {'xdata','ydata'};
    pv = mat2cell(Raw.data.segments.coord([1 3 2 4],:)',ones(1,m),[2 2]);

    set(SenFig.raw.segments(1:m),...
        pn,     pv,...
        'vis',  'on')
    set(SenFig.raw.segments(m+1:end),...
        'vis',  'off')

end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function dist = imDistance(varargin)
%IMDISTANCE Distance map computed from a set of points
%
%   DIST = imDistance(DIM, POINTS);
%   Create a distance map of the point given by POINTS, in the image with
%   size given by DIM = [NY NX].
%   POINTS is a N-by-2 array of double containing point coordinates
%
%   DIST = imDistance(LX, LY, POINTS)
%   Specify the position of vertices along each axis. See meshgrid for
%   details.
%
%   DIST = imDistance(..., N)
%   Use N point randomly and independently created inside the image bounds.
%
%   DIST = imDistance(..., EDGE)
%   Also specify edge condition. EDGE can be one of:
%   - 'free': regions touching image borders are removed
%   - 'periodic': copy the set of points to simulate periodic set
%   - 'remove': image edges will be set to distance 0, and neighbor pixels
%       set accordingly.
%
%   TODO: manage different types of distance (now is only euclidian).
%
%   Example
%   % basic example with integer coordinates
%     img = imDistance([100 100], [10 10;90 90;90 10;10 90;50 50]);
%     imshow(img, [0 max(img(:))]);
%
%   % specify floating coordinates, with periodic edge conditions
%     lx = 1:100;
%     ly = 1:100;
%     pts = [10.2 10.3;90.4 90.5; 10.8 90.8;50.1 50.0]
%     img = imDistance(lx, ly, pts, 'periodic');
%     figure; imshow(img, [0 max(img(:))]);
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 16/06/2004.
%

%   HISTORY
%   28/06/2004 use faster algorithm
%   22/05/2006 allocate d1 and d2 only once (should be faster)
%   29/05/2009 add possibility to specify grid with meshgrid-like syntax

%% Extract input arguments

% default parameters values
points      = [];       % points array
N           = 20;       % default number of germs
edgecond    = 'free';   % edge condition
lx = 1:100;
ly = 1:100;

% extraction of image dimensions
if ~isempty(varargin)
    var = varargin{1};
    if size(var, 1) > 1 && size(var, 2) > 2
        % case of a 2x3 matrix with starting position, increment, end
        % position for each coordinate
        lx = var(1,1):var(1,2):var(1,3);
        ly = var(2,1):var(2,2):var(2,3);
        varargin(1) = [];

    elseif size(var, 1) == 1 && size(var, 2) == 2
        % first argument contains the size of the output image
        lx = 1:var(2);
        ly = 1:var(1);
        varargin(1) = [];

    elseif length(varargin) > 1
        % first and second arguments contain vector for each coordinate
        % respectively
        lx = varargin{1};
        ly = varargin{2};
        varargin(1:2) = [];

    else
        error('wrong input arguments in imDistance');
    end
end

% extraction of points, or number of points
if ~isempty(varargin)
    var = varargin{1};

    if length(var)>1    % use an array of germs
        points = var;
    else                % get number of germs
        N = var;
    end

    varargin(1) = [];
end

% extraction of edge condition
if ~isempty(varargin)
    edgecond = varargin{1};
end

%% Initialisations

% create array of points, if does not exist
if isempty(points)
    % generation of germs array
    points = zeros(N, 2);
    for i = 1:N
        points(i,1) = rand * lx(end) + lx(1);
        points(i,2) = rand * ly(end) + ly(1);
    end
end

% number of points
N = size(points, 1);

% size of output image
Nx = length(lx);
Ny = length(ly);

% create periodic conditions if needed
if strcmp(edgecond, 'periodic')
    % compute size of the box
    wx = lx(end)  + lx(2) - 2*lx(1);
    wy = ly(end)  + ly(2) - 2*ly(1);

    points = repmat(points, [9 1]);
    for i = [1 2 3]
        points((i-1)*N+1:i*N, 1) = points((i-1)*N+1:i*N, 1) - wx;
    end

    for i = [7 8 9]
        points((i-1)*N+1:i*N, 1) = points((i-1)*N+1:i*N, 1) + wx;
    end

    for i = [1 4 7]
        points((i-1)*N+1:i*N, 2) = points((i-1)*N+1:i*N, 2) - wy;
    end

    for i = [3 6 9]
        points((i-1)*N+1:i*N, 2) = points((i-1)*N+1:i*N, 2) + wy;
    end

    % re-compute number of points
    N = size(points, 1);
end

%% Generation of distance function

% allocate memory
dist = Nx * Ny * ones([Ny Nx]);

% pixels coordinates
[x, y] = meshgrid(lx, ly);

% update distance for each point
for p = 1:N
    dist = min(dist, hypot(x - points(p,1), y - points(p,2)));
end

%old algorithm, slower
%dist = zeros(dim);
%for i=1:dim(1)
    %disp (int2str(x));
    %for j=1:dim(2)
        %mat = points - ones(N,1)*[i j];
        %dist(i, j) = min(sqrt(diag(mat*mat')));
        % end
    %end

%% Optional processing for borders

% distance from borders of image
if strcmp(edgecond, 'remove')
    d1 = min(repmat((1:Ny)'-1, [1 Nx]), repmat((Ny:-1:1)'-1, [1 Nx]));
    dist = min(dist, d1);
    d1 = min(repmat((1:Nx)-1,  [Ny 1]), repmat((Nx:-1:1)-1,  [Ny 1]));
    dist = min(dist, d1);
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [ AllCuboid ] = regionBasedHypothesisFromSS( SS1, imgH, imgW )
%REGIONBASEDHYPOTHESISFROMSS Convert selective search result to cuboid
%   SS1: selective search result: hierarchy grouping of over-segment
%   AllCuboid:
%       xyzbox: write coordinates of 4 corners to a vector
%       views: visible normal direction of surfaces
%       score: intersection/union with segmentation

% [imgH, imgW, ~] = size(rotImg);
vanishingPoint = [-1 0 0; 1 0 0; 0 -1 0; 0 1 0;0 0 -1; 0 0 1];
viewID = [1 2 3 4 5 6];
colorTypes = {'Hsv', 'Lab', 'RGI', 'H', 'Intensity'};

%% extract ss cuboid candidate
AllCuboid = struct('views',zeros(2000,3),'xyzBox',zeros(2000,36),'score',zeros(2000,1),'count',0);
minAngleDist = pi/18;

for colorid = 1:length(colorTypes)
    hBlobs1 = RecreateBlobHierarchyIndIm( SS1(colorid).blobIndIm, SS1(colorid).blobBoxes, SS1(colorid).hierarchy{1});
    hBlobs2 = RecreateBlobHierarchyIndIm( SS1(colorid).blobIndIm, SS1(colorid).blobBoxes, SS1(colorid).hierarchy{2});
    hBlobs = [hBlobs1; hBlobs2];
    [ hBlobs_valid, hBlobs_centroid ] = regionShapeValidation( hBlobs, 1024, 2048, 2500, 3, 0.65, true );

    hBlobs_ViewPoint = coords2uv( hBlobs_centroid, imgW, imgH );

    yValid = hBlobs_ViewPoint(:,2)/(minAngleDist);
    yValid(abs(yValid)>=1) = 0;
    hBlobs_ViewPoint(yValid~=0,2) = minAngleDist*sign(yValid(yValid~=0));

    for j = 1:size(hBlobs_ViewPoint,1)
        xClose = hBlobs_ViewPoint(j,1) - [-pi -pi/2 0 pi/2 pi];
        [B,I] = min(abs(xClose));
        if B<minAngleDist
            hBlobs_ViewPoint(j,1) = -pi + (I-1)*pi/2 + sign(xClose(I))*minAngleDist;
        end
    end

    fov = 2/3*pi;
    fovSize = 640;
    big_map = zeros(imgH, imgW);

    for bid = find(hBlobs_valid')
        rect = hBlobs{bid}.rect;
        map = big_map;
        map(rect(1):rect(3),rect(2):rect(4)) = hBlobs{bid}.mask;
        warped_image = imgLookAt(map, hBlobs_ViewPoint(bid,1), hBlobs_ViewPoint(bid,2), fovSize, fov );

        fprintf('%d/%d, ', bid, length(hBlobs));
        vp = uv2xyzN(hBlobs_ViewPoint(bid,:));
	fprintf('1');
        [out2D, valid, division] = projectPoint2SeparateView( vanishingPoint, vp, fov, fovSize );
	fprintf('2');
        vanishing = out2D(valid,:);
	fprintf('3');
        subViewID = viewID(valid);
	fprintf('4');
        [ cuboid ] = getSegHypothesisID( warped_image, 200, vanishing, 1);

	fprintf('Get cuboid\n');
        for k = 1:length(cuboid)
            ID = AllCuboid.count + k;
            AllCuboid.score(ID) = cuboid(k).score;
            if cuboid(k).type==1
                points = cuboid(k).points;
                startID = cuboid(k).startID;
                surf1 = points(rem(([startID(1):startID(1)+3]-1),6)+1,:);
                surf2 = points(rem(([startID(2):startID(2)+3]-1),6)+1,:);
                [ out3DNorm1, ~ ] = projectPointFromSeparateView(  surf1, vp, fov, fovSize );
                [ out3DNorm2, ~ ] = projectPointFromSeparateView(  surf2, vp, fov, fovSize );
                AllCuboid.xyzBox(ID,1:24) = reshape([out3DNorm1;out3DNorm2]', 1, []);
                AllCuboid.views(ID,1:2) = subViewID(cuboid(k).direct(1:2));
            elseif cuboid(k).type==3
                points = cuboid(k).points;
                startID = cuboid(k).startID;
                surf1 = [points(rem(([startID(1):startID(1)+2]-1),6)+1,:);cuboid(k).addPoint];
                surf2 = [points(rem(([startID(2):startID(2)+2]-1),6)+1,:);cuboid(k).addPoint];
                surf3 = [points(rem(([startID(3):startID(3)+2]-1),6)+1,:);cuboid(k).addPoint];
                [ out3DNorm1, ~ ] = projectPointFromSeparateView(  surf1, vp, fov, fovSize );
                [ out3DNorm2, ~ ] = projectPointFromSeparateView(  surf2, vp, fov, fovSize );
                [ out3DNorm3, ~ ] = projectPointFromSeparateView(  surf3, vp, fov, fovSize );
                AllCuboid.xyzBox(ID,1:36) = reshape([out3DNorm1;out3DNorm2;out3DNorm3]', 1, []);
                AllCuboid.views(ID,:) = subViewID(cuboid(k).direct);
            end

        end
        AllCuboid.count = AllCuboid.count + length(cuboid);

    end

end

AllCuboid.views = AllCuboid.views(1:AllCuboid.count,:);
AllCuboid.xyzBox = AllCuboid.xyzBox(1:AllCuboid.count,:);
AllCuboid.score = AllCuboid.score(1:AllCuboid.count);

end

author: zouchuhang

repo: LayoutNet


text:
function [minPoint,minValue,exitCode]=goldenSectionSearch(g,xSpan,XTol,maxIter)
%%GOLDENSECTIONSEARCH Perform a golden section search to find the minimum
%                     of a continuous, real, function that is unimodal in a
%                     given interval. The golden section search linearly
%                     decreases the search region. The algorithm is
%                     considered more efficient than a ternary search.
%                     Unlike Matlab's built-in fminbnd function, quadratic
%                     interpolation is not used to try to speed up the
%                     convergence rate.
%
%INPUTS: g The handle to a real function to minimize. The function takes a
%          scalar input and returns a real scalar output.
%    xSpan The 2X1 or 1X2 span of values in which the minimum is to be
%          found xSpan(1)<xSpan(2).
%     XTol The absolute tolerance on the argument of g. When the region of
%          uncertainty of the minimum goes below xTol, then convergence is
%          declared. If omitted or an empty matrix is passed, XTol=1e-9 is
%          used.
%  maxIter The maximum number of iterations to use. If this parameter is
%          omitted or an empty matrix is passed, then maxIter=100; is used.
%          One to 2 function evaluations are performed each iteration.
%
%OUTPUTS: minPoint The value such that g(minVal) is minimized after the
%                  search.
%         minValue The value g(minVal).
%         exitCode A value indicating how the function terminated. Possible
%                  values are:
%                  0 The tolerance XTol in the argument of g was achieved.
%                  1 The maximum number of iterations was reached.
%
%The golden section search algorithm is implemented as described in
%Appendix C of [1]. The convergence rate is linear. This function evaluates
%g at the bounds of the interval.
%
%As an example, consdier finding the minimum of the function
%g(x)=(x-3)*x^3*(x-6)^4.
%This function has multiple minima and maxima. However, if one minimum can
%be bounded, then the golden section search can be used to find it. If
%multiple minima are bounded, then the golden section search can be used,
%but it might not convergence to the best minimum. We will consider the
%the case of three ranges.
% %Example 1
% alphaSpan=[0;3];
% g=@(x)(x-3)*x^3*(x-6)^4;
% [minPoint,minValue]=goldenSectionSearch(g,alphaSpan);
% %Here we find a minPoint of about 1.7354 and a min value of about -2186.1.
% %Next, we consider bounding a second minimum.
% %Example 2
% alphaSpan=[3;7];
% [minPoint1,minValue1]=goldenSectionSearch(g,alphaSpan);
% %Here, we get a minimum point of 6 with a value of 0. Now, if we take a
% %span that covers both minimum:
% %Example 3:
% alphaSpan=[1;7];
% [minPoint2,minValue2]=goldenSectionSearch(g,alphaSpan);
% %We also get a minimum point of 6 even though that is not the absolute
% %minimum over the entire region. This is because the entire region is not
% %unimodal. If we used alphaSpan=[0;7]; then the first minimum (the lower
% one) in that span would have been found.
%
%REFERENCES:
%[1] D. P. Bertsekas, Nonlinear Programming, 2nd ed. Belmont, MA: Athena
%    Science, 1999.
%
%January 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(maxIter))
    maxIter=100;
end

if(nargin<3||isempty(XTol))
	XTol=1e-9;
end

tau=(3-sqrt(5))/2;

alpha=xSpan(1);
alphaBar=xSpan(2);

if(alphaBar<alpha)
   error('alphaSpan(2) is less than alphaSpan(1)');
end

b=alpha+tau*(alphaBar-alpha);
bBar=alphaBar-tau*(alphaBar-alpha);

gAlpha=g(alpha);
gAlphaBar=g(alphaBar);

gb=g(b);
gbBar=g(bBar);

exitCode=1;
for curIter=1:maxIter
    alphaPrev=alpha;
    alphaBarPrev=alphaBar;
    if(gb<gbBar)
        if(gAlpha<=gb)
            alphaBar=b;
            gAlphaBar=gb;

            bBar=alphaBar-tau*(alphaBar-alpha);
            gbBar=g(bBar);
        else
            alphaBar=bBar;
            gAlphaBar=gbBar;

            %Because of the use of the Golden ratio for tau, only one
            %function evaluation is needed in this instance.
            bBar=b;
            gbBar=gb;
        end

        b=alpha+tau*(alphaBar-alpha);
        gb=g(b);
    elseif(gb>gbBar)
        if(gbBar>=gAlphaBar)
            alpha=bBar;
            gAlpha=gbBar;

            b=alpha+tau*(alphaBar-alpha);
            gb=g(b);
        else
            alpha=b;
            gAlpha=gb;

            %Because of the use of the golden ratio for tau, only one
            %function evaluation is needed in this instance.
            b=bBar;
            gb=gbBar;
        end

        bBar=alphaBar-tau*(alphaBar-alpha);
        gbBar=g(bBar);
    else
        alpha=b;
        alphaBar=bBar;

        gAlpha=gb;
        gAlphaBar=gbBar;

        b=alpha+tau*(alphaBar-alpha);
        gb=g(b);
        bBar=alphaBar-tau*(alphaBar-alpha);
        gbBar=g(bBar);
    end

    %If the accuracy bound has been achieved. The second condition deals
    %with finite precision limitiations if XTol is very close to 0.
    if(abs(alphaBar-alpha)<=XTol||(alphaPrev==alpha&&alphaBarPrev==alphaBar))
        exitCode=0;
        break;
    end
end

%After the iterations, we have the function evaluated at alpha, b, bBar,
%and alphaBar. We will just return the point/ value that is the smallest.

points=[alpha;b;bBar;alphaBar];
values=[gAlpha;gb;gbBar;gAlphaBar];

[values,idx]=sort(values,'ascend');
minPoint=points(idx(1));
minValue=values(1);
end

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

repo: TrackerComponentLibrary


text:
function plot_GridPointAnalyse(params,nMode,nPlot)
% Example : plot_GridPointAnalyse(params,0,2)
% function to analyse grid points and plots distribution of rate change
% values
%
% Input
% params    : Input matrix with result
% nMode     : calculate gridpoint from synthetic PSQ (0) or ask for grid
%               point (1)
% nPlot     : plots Z (1), sigma(Z) (2), beta (3), sigma(beta) (4)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[vIn,vBnd] = inpoly(params.mPolygon,params.mPSQ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  controlling inpoly
figure;
plot(params.mPolygon(:,1),params.mPolygon(:,2),'.');
hold on;plot(params.mPSQ(:,1),params.mPSQ(:,2),'r')
hold on;plot(params.mPolygon(vIn,1),params.mPolygon(vIn,2),'ro')
hold on;plot(params.mPolygon(vBnd,1),params.mPolygon(vBnd,2),'ko')
%%

switch nPlot
    case 1
        mPlot=squeeze(params.mResult1(:,1,:));
    case 2
        mPlot=calc_ProbColorbar2Value(squeeze(params.mResult2(:,1,:)));
    case 3
        mPlot=squeeze(params.mResult3(:,1,:));
    case 4
        mPlot=calc_ProbColorbar2Value(squeeze(params.mResult4(:,1,:)));
end

mPlot=mPlot(vIn,:);
mPlot=reshape(mPlot,size(mPlot,1)*size(mPlot,2),1);
figure;
[n,xout] = hist(mPlot(~isinf(mPlot)));
plot(xout,n,'r');

author: CelsoReyes