text:
% nbins = 8;
% psize = 10;
% npatches = 8;
% optflowwinsig = 3;
% optflowsig = 2;
% optreliability = 1e-4;
% patchsz = psize*npatches;

expdir = '/groups/branson/home/kabram/bransonlab/data/katie/20170308_133311/20170308_133311_08';
moviefilestr = 'movie.avi';
frames = 2831:2900;
vidoutfile = sprintf('%s_From%d_To%d_HOF_%s.avi',expname,fly,t0,t1,datestr(now,'yyyymmdd'));
makeVideo = true;

%%
params = getParamsKatie;
npatchesx = params.npatchesx;
npatchesy = params.npatchesy;
psize = params.psize;
nbins = params.nbins;
scale = params.scale;

optflowwinsig = params.optflowwinsig ;
optflowsig = params.optflowsig ;
optreliability = params.optreliability ;

maxflow = 5;
fly_thres = 90;
flow_thres = params.flow_thres;

% Antennal grooming for andy

[~,expname] = fileparts(expdir);
moviefile = fullfile(expdir,moviefilestr);

[readframe,nframes] = get_readframe_fcn(moviefile);
for i = 1:numel(frames)
  curf = readframe(frames(i));
  nextf = readframe(frames(i)+1);

  curpatch = curf;
  curpatch2 = nextf;

  im1(:,:,i) = curpatch;
  im2(:,:,i) = curpatch2;
end

[nr,nc,~] = size(im1);
%
% figure out the histogram bins
tmptheta = (0:179)*pi/180;
res = nan(nbins,numel(tmptheta));
m = single(ones(10,10));
o = single(zeros(10,10));
for i = 1:numel(tmptheta),
%   fprintf('i = %d, theta = %f\n',i,tmptheta(i));
  o(:) = single(tmptheta(i));
  rescurr = gradientHist(m,o,1,nbins,1);
  res(:,i) = rescurr(1,1,:);
end

bincenters = nan(1,nbins);
for i = 1:nbins,
  bincenters(i) = tmptheta(argmax(res(i,:)));
end

% this seems to be what the centers correspond to
bincenters = linspace(0,pi,nbins+1);
bincenters = bincenters(1:nbins);

% mayank divides by 2
bincenters2 = bincenters*2;
dt = mean(diff(bincenters2));
binedges2 = [bincenters2(1)-dt/2,(bincenters2(1:end-1)+bincenters2(2:end))/2,bincenters2(end)+dt/2];

%% make a video

tr = 191;
hfig = 100;
t0 = frames(1);
t1 = frames(end);
colorpos = [1,0,0];
colorneg = [0,.3,1];

maxv2 = 0;
colors = hsv(nbins);
sz = size(im1);
bwimg = zeros(sz(1),sz(2));
ctr = [ceil( (sz(1)+1)/2),ceil( (sz(2)+1)/2)];
bwimg(ctr(1),ctr(2))=1;
dimg = bwdist(bwimg,'euclidean');
[xx,yy]= meshgrid(1:sz(2),1:sz(1));
aimg = atan2(-(yy-ctr(1)),-(xx-ctr(2)));

for t = t0:20:t1,

  im1curr = gaussSmooth(im1(:,:,t-t0+1),1,'same');
  im2curr = gaussSmooth(im2(:,:,t-t0+1),1,'same');
  uv = estimate_flow_interface(im1curr,im2curr,'hs-brightness',{'max_warping_iters',2});
  Vx = uv(:,:,1); Vy = uv(:,:,2);

  M = sqrt(Vx.^2 + Vy.^2);
  O = mod(atan2(Vy,Vx)/2,pi);
  O = min(O,pi-1e-6);
  H = gradientHist(single(M),single(O),psize,nbins,1);
  maxv2 = max(maxv2,max(H(:)));
end

%%

maxv2 = 4;

if makeVideo
  vid = VideoWriter(vidoutfile);
  open(vid);
end

for t = t0:t1

  im1curr = gaussSmooth(im1(:,:,t-t0+1),1,'same');
  im2curr = gaussSmooth(im2(:,:,t-t0+1),1,'same');

  imsz = size(im1curr);
  pairimg = zeros(2*imsz(1),imsz(2),3);
  ttimg = [];
  tt(:,:,1) = im2curr;
  tt(:,:,2) = im1curr;
  tt(:,:,3) = im2curr;
  pairimg(1:imsz(1),:,:) = tt;

  uv = estimate_flow_interface(im1curr,im2curr,'hs-brightness',...
    {'max_warping_iters',2 });
  uvorig = uv;
  pairimg(imsz(1)+(1:imsz(1)),1:imsz(2),:) = flowToColor(uv,maxflow);

  pairimg = uint8(pairimg);
  pairimg = imresize(pairimg,scale);
  Vx = uvorig(:,:,1); Vy = uvorig(:,:,2);
  M = sqrt(Vx.^2 + Vy.^2);
  O = mod(atan2(Vy,Vx)/2,pi);
  O = min(O,pi-1e-6);
  Horig = gradientHist(single(M),single(O),psize,nbins,1);

  if t == t0,
    figure(hfig);
    clf;
    hax = axes('Position',[0,0,1,1]);
    set(hfig,'Units','pixels','Position',get(0,'ScreenSize'));

    him = imshow(pairimg);
    axis image;
    truesize;
    colormap gray;
    hold on;
    axis off;
    htext = text(1,nr-1,sprintf('%.2f s',(t-1)/200),'Color','k','HorizontalAlignment','left','VerticalAlignment','bottom','FontSize',18);

  else
    delete(horig(ishandle(horig)));
    set(him,'CData',pairimg);
    set(htext,'String',sprintf('%.2f s',(t-1)/200));
  end

  horig = plotHofFeatures(Horig,hax,maxv2,colors,binedges2,psize*scale,[0 0]);

  drawnow;
  if makeVideo
    fr = getframe(hax);
    writeVideo(vid,fr);
  elseif t0<t1
    pause;
  end
end

if makeVideo,  close(vid); end

author: kristinbranson


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = MC_merton_sgs(S,r,sigma,T, a,b,lambda,NSim)
%Implements the SGS method for the Merton model
nu = r - lambda*(exp(a+0.5*b^2)-1)-0.5*sigma^2; % martingale correction

X = log(S)*ones(NSim,1);          % X is the log price path

for k=1:NSim
    t = 0;
    tau = [];
    %simulate the jump times first
    while t < T
        dt = -log(rand)/lambda;% jump time
        t = t + dt;            % add the jump times
        tau = [tau; dt];             % not good matlab but works!
    end
    tau(end) = T-(t-dt);
    N = length(tau);                    % N number of jumps + 1
    W1 = randn(1,N);                    % uniforms for diffusion
    if N > 1
        W2 = randn(1,N-1);              % uniforms for jumps
        for i = 1:N-1
            Z = nu*tau(i) + sigma*sqrt(tau(i)) * W1(i);
            lnY = a+b*W2(i);
            X(k) = X(k) + Z + lnY;
        end
    end

    dt_end = tau(end);
    Z = nu * dt_end + sigma * sqrt(dt_end) * W1(end);
    X(k) = X(k) + Z;
end
y = exp(X);
end

author: Sable

repo: mcbench-benchmarks


text:
function pass = test_definePoint(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

% Test a scalar-valued CHEBFUN object:
f = chebfun(@(x) x, [-1, 0, 1], pref);
% Use DEFINEPOINT:
f = definePoint(f, .5, 1);
pass(1) = all(f.domain == [-1, 0, .5, 1]) && feval(f, .5) == 1;
% Use SUBSASGN:
f(-.5) = 2;
pass(2) = all(f.domain == [-1, -.5, 0, .5, 1]) && feval(f, -.5) == 2;

% Test an array-valued CHEBFUN object:
f = chebfun(@(x) [x, x], [-1, 0, 1], pref);
% Use DEFINEPOINT:
f = definePoint(f, .5, [1, 2]);
pass(3) = all(f.domain == [-1, 0, .5, 1]) && all(feval(f, .5) == [1, 2]);
% Use SUBSASGN:
f(-.5) = [2, 3];
pass(4) = all(f.domain == [-1, -.5, 0, .5, 1]) && all(feval(f, -.5) == [2, 3]);

% Test an array-valued CHEBFUN object at a number of points:
f = chebfun(@(x) [x, x], [-1, 0, 1], pref);
% Use DEFINEPOINT:
f = definePoint(f, [-.25, .5], [1, 2 ; 3 4]);
pass(5) = all(size(f.pointValues) == [5, 2]) && ...
    all(all(f.pointValues == [-1 -1 ; 1 2 ; 0 0 ; 3 4 ; 1 1]));
% Use SUBSASGN:
f([-.25, .5]) = [1, 2 ; 3 4];
pass(6) = all(size(f.pointValues) == [5, 2]) && ...
    all(all(f.pointValues == [-1 -1 ; 1 2 ; 0 0 ; 3 4 ; 1 1]));

% Test an array-valued CHEBFUN object at a number of points (vector expansion):
f = chebfun(@(x) [x, x], [-1, 0, 1], pref);
% Use DEFINEPOINT:
f = definePoint(f, [-.25, .5], [1, 2]);
pass(7) = all(size(f.pointValues) == [5, 2]) && ...
    all(all(f.pointValues == [-1 -1 ; 1 2 ; 0 0 ; 1 2 ; 1 1]));
% Use SUBSASGN:
f([-.25, .5]) = [1, 2];
pass(8) = all(size(f.pointValues) == [5, 2]) && ...
    all(all(f.pointValues == [-1 -1 ; 1 2 ; 0 0 ; 1 2 ; 1 1]));

% Test an array-valued CHEBFUN object at a number of points (scalar expansion):
f = chebfun(@(x) [x, x], [-1, 0, 1], pref);
% Use DEFINEPOINT:
f = definePoint(f, [-.25, .5], 1);
pass(9) = all(size(f.pointValues) == [5, 2]) && ...
    all(all(f.pointValues == [-1 -1 ; 1 1 ; 0 0 ; 1 1 ; 1 1]));
% Use SUBSASGN:
f([-.25, .5]) = 1;
pass(10) = all(size(f.pointValues) == [5, 2]) && ...
    all(all(f.pointValues == [-1 -1 ; 1 1 ; 0 0 ; 1 1 ; 1 1]));

%% Test on singular function: piecewise smooth chebfun

% define the domain:
dom = [-2 -1 0 1];

op1 = @(x) sin(x);
op2 = @(x) 1./((1+x).^0.5);
op3 = @(x) x+1;
op = {op1, op2, op3};

f = chebfun(op, dom, 'exps', [0 0 -0.5 0 0 0]);
brkpts = zeros(1,2);
brkpts(1) = -0.7;
brkpts(2) = -0.5;
g = definePoint(f, brkpts(1), 1);
g(brkpts(2)) = 2;

% check values:
check = zeros(1,4);
check(1) = all(g.domain == unique([dom, brkpts]));
check(2) = feval(g, brkpts(1)) == 1;
check(3) = feval(g, brkpts(2)) == 2;
check(4) = all(g.pointValues == [f.pointValues(1:2); 1; 2; f.pointValues(3:4)]);

pass(11) = all( check );

%% Test for function defined on unbounded domain:

% Functions on [-inf inf]:

% Set the domain:
dom = [-Inf Inf];

% Blow-up function:
op = @(x) (1-exp(-x.^2));
f = chebfun(op, dom);
brkpts = zeros(1,2);
brkpts(1) = -7;
brkpts(2) = 2;
g = definePoint(f, brkpts(1), 3);
g(brkpts(2)) = 4;

% check values:
check = zeros(1,4);
check(1) = all(g.domain == unique([dom, brkpts]));
check(2) = feval(g, brkpts(1)) == 3;
check(3) = feval(g, brkpts(2)) == 4;
check(4) = all(g.pointValues == [f.pointValues(1); 3; 4; f.pointValues(end)]);

pass(12) = all( check );

%% Test assigning to an endpoint (#896)
f = chebfun(0, [-1, 1]);
f(1) = 1;
pass(13) = abs(feval(f, 1) - 1) < eps;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function plot_3d_odomap_without_groundtruth(odoMapFile, ...
    xOdoMapScaling, yOdoMapScaling, zOdoMapScaling, ...
    xOdoMapTrans, yOdoMapTrans, zOdoMapTrans, ...
    xGtScaling, yGtScaling, zGtScaling)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     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/>.

    % load ground truth data
%     [frameId, gt_x, gt_y, gt_z, gt_rx, gt_ry, gt_rz] = load_ground_truth_data(groundTruthFile);

    % load experience map data
    [id, odomap_x, odomap_y, odomap_z, odomap_yaw, odomap_pitch] = load_exp_map_data(odoMapFile);

    figure('color',[1 1 1]);
    hold on

    plot3(odomap_y * yOdoMapScaling + yOdoMapTrans, ...
        odomap_x * xOdoMapScaling + xOdoMapTrans, ...
        odomap_z * zOdoMapScaling + zOdoMapTrans, '.b');
%     plot3((gt_x - gt_x(1)) * xGtScaling,(gt_y - gt_y(1)) * yGtScaling, (gt_z - gt_z(1)) * zGtScaling, '.r');
    hold off;

%     view(3)
view(0, 90);
    grid on
%     xlabel('x');
%     ylabel('y');
    xl = xlabel('x');
    yl = ylabel('y');
    zlabel('z');
    set(xl,'Rotation',15);
    set(yl,'Rotation',-30);
%     title('3D Experience Map');
    %                     legend('Result','Truth' ,'1');
    % axis([-10 20 -10 20 -10 20]);
    %                     axis equal
%     legend('Experience map','Ground truth', 'location', '0');
    hh = legend('Odomap','Ground truth');

    set(hh,'position',[.76 .82 .05 .05]);
    axis on
    rotate3d on
end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
%COVM Compute covariance matrix for large datasets
%
% 	C = COVM(A)
%
% Similar to C = COV(A) this routine computes the covariance matrix
% for the datavectors stored in the rows of A. No large intermediate
% matrices are created. If class(A) is 'prdataset' then class(C) is
% 'double'.

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

% $Id: covm.m,v 1.3 2010/02/08 15:34:14 duin Exp $

function c = covm(a,n)
		[m,k] = size(a);
if nargin < 2, n = 0; end
if n ~= 1 & n ~= 0
	error('Second parameter should be either 0 or 1')
end
[loops,n0,n1] = prmem(m,k);
if loops == 1
	c = prcov(+a,n);
	c = (c+c')/2;
	return
end
c = zeros(k,k);
u = ones(n0,1)*mean(a);
for j = 1:loops
	if j == loops, n = n1; else n = n0; end
	nn = (j-1)*n0;
	b = +a(nn+1:nn+n,:) - u(1:n,:);
	c = c + b'*b;
end
c = (c + c')/2;
if n
	c = c/m;
else
	c = c/(m-1);
end

author: marianux

repo: ecg-kit

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


text:
function [r,c,V] = findnearest(srchvalue,srcharray,bias)

% Usage:
% Find the nearest numerical value in an array to a search value
% All occurances are returned as array subscripts
%
% Output:
%
% For 2D matrix subscripts (r,c) use:
%
%       [r,c] = findnearest(srchvalue,srcharray,gt_or_lt)
%
%
% To also output the found value (V) use:
%
%       [r,c,V] = findnearest(srchvalue,srcharray,gt_or_lt)
%
%
% For single subscript (i) use:
%
%         i   = findnearest(srchvalue,srcharray,gt_or_lt)
%
%
% Inputs:
%
%    srchvalue = a numerical search value
%    srcharray = the array to be searched
%    bias      = 0 (default) for no bias
%                -1 to bias the output to lower values
%                 1 to bias the search to higher values
%                (in the latter cases if no values are found
%                 an empty array is ouput)
%
%
% By Tom Benson (2002)
% University College London
% t.benson@ucl.ac.uk

if nargin<2
    error('Need two inputs: Search value and search array')
elseif nargin<3
    bias = 0;
end

% find the differences
srcharray = srcharray-srchvalue;

if bias == -1   % only choose values <= to the search value

    srcharray(srcharray>0) =inf;

elseif bias == 1  % only choose values >= to the search value

    srcharray(srcharray<0) =inf;

end

% give the correct output
if nargout==1 | nargout==0

    if all(isinf(srcharray(:)))
        r = [];
    else
        r = find(abs(srcharray)==min(abs(srcharray(:))));
    end

elseif nargout>1
    if all(isinf(srcharray(:)))
        r = [];c=[];
    else
        [r,c] = find(abs(srcharray)==min(abs(srcharray(:))));
    end

    if nargout==3
        V = srcharray(r,c)+srchvalue;
    end
end


text:
function f = sym2func(fun,args)
%SYM2FUNC  Convert a symbolic variable equation to function handle
%
%   f = sym2func(fun) converts the symbolic equation into a string, then
%   replaces any x1 to x(1), x2 to x(2), etc, and returns it as a function
%   handle suitable for standard evaluation.

%   Copyright (C) 2013 Jonathan Currie (IPL)

if(nargin < 2), args = 'x'; end

%Extract symbolic variables in the expression
sv = symvar(fun);
%For each symbolic var, convert x1 to x(1) as cell of strings
cv = cell(size(sv));
for i = 1:length(sv)
    cv{i} = symVar2MVar(sv(i));
end
%Substitute our indexed variables into the symbolic expression
sfun = subs(fun,sv,cv);

%Convert to string
str = char(sfun);
%Remove 'matrix' if it is found
k = strfind(str,'matrix');
if(k)
    str = str(k+7:end-1);
end
%Ensure we have right size
[r,c] = size(fun);
if(c==1 && r > 1)
    str = regexprep(str,',',';'); %force column
elseif(c>1 && r>1)
    str = regexprep(str,'],','];'); %keep matrix
end

if(iscell(args))
    a = args{1};
    for i = 2:length(args)
        a = sprintf('%s,%s',a,args{i});
    end
    args = a;
end

%Build function handle
f = eval(['@(' args ') ' str]);

function str = symVar2MVar(svar)
s = char(svar);
for i = 2:length(s)
    if(~isalpha(s(i)))
        str = [s(1:i-1) '(' s(i:end) ')'];
        return;
    end
end
%No index?, oh well... could be for e.g. just t
str = s;

function is = isalpha(str)
if((str>='a' && str<='z') || (str>='A' && str<='Z'))
    is = true;
else
    is = false;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function f = minus(f, g)
%-   Subtraction of two CHEBFUN3T objects.
%   F - G subtracts G from F, where F and G are CHEBFUN3T objects or
%   scalars.
%
% See also CHEBFUN3T/PLUS and CHEBFUN3T/UMINUS.

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

% f - g = f + (-g)
f = plus(f, uminus(g));

end


text:
function test_suite = test_index_unique
% tests for cosmo_index_unique
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_index_unique_()
    n=100;
    a=round(rand(n,1)*2)+1;
    b=round(rand(n,1)*2)+1;
    c=round(rand(n,1)*2)+1;

    abc=[a b c];
    [unq1,unused,idx1]=unique(abc,'rows');
    [idx2,vals2]=cosmo_index_unique({a,b,c});

    assertEqual(unq1,[vals2{:}]);

    nunq=size(unq1,1);
    assert(nunq==numel(idx2));

    % test some random elements
    rp=randperm(nunq);
    rp=rp(1:10);
    for j=rp
        assertEqual(find(idx1==j),idx2{j})
    end

    % test with some string labels
    labels={'a','bb','ccc'}';

    [idx3,vals3]=cosmo_index_unique({labels(a),b,labels(c)});
    assertEqual(idx3,idx2);
    v2to3={labels(vals2{1}),vals2{2},labels(vals2{3})};
    assertEqual(vals3,v2to3)

    % test with matrix input
    [idx4,vals4]=cosmo_index_unique(abc);
    assertEqual(idx4,idx2);
    assertEqual(vals4,[vals2{:}]);

    % test vector input
    rp=randperm(10);
    [idx5,vals5]=cosmo_index_unique(rp+20);
    assertEqual(idx5,{1});
    assertEqual(vals5,rp+20);

    [idx6,vals6]=cosmo_index_unique(rp'+20);
    assertEqual(cell2mat(idx6(rp))',1:numel(rp));
    assertEqual(vals6,(1:numel(rp))'+20);

    % test empty input

    % test wrong inputs
    n=100;
    a=round(rand(n,1)*2)+1;
    b=round(rand(n,1)*2)+1;
    c=round(rand(n,1)*2)+1;

    abc=[a b c];
    [unq1,unused,idx1]=unique(abc,'rows');
    [idx2,vals2]=cosmo_index_unique({a,b,c});

    assertEqual(unq1,[vals2{:}]);

    nunq=size(unq1,1);
    assert(nunq==numel(idx2));

    % test some random elements
    rp=randperm(nunq);
    rp=rp(1:10);
    for j=rp
        assertEqual(find(idx1==j),idx2{j})
    end

    % test with some string labels
    labels={'a','bb','ccc'}';

    [idx3,vals3]=cosmo_index_unique({labels(a),b,labels(c)});
    assertEqual(idx3,idx2);
    v2to3={labels(vals2{1}),vals2{2},labels(vals2{3})};
    assertEqual(vals3,v2to3)

    % test with matrix input
    [idx4,vals4]=cosmo_index_unique(abc);
    assertEqual(idx4,idx2);
    assertEqual(vals4,[vals2{:}]);

    % test vector input
    rp=randperm(10);
    [idx5,vals5]=cosmo_index_unique(rp+20);
    assertEqual(idx5,{1});
    assertEqual(vals5,rp+20);

    [idx6,vals6]=cosmo_index_unique(rp'+20);
    assertEqual(cell2mat(idx6(rp))',1:numel(rp));
    assertEqual(vals6,(1:numel(rp))'+20);

    % test empty input

    % test wrong inputs
    aet=@(x)assertExceptionThrown(@()cosmo_index_unique(x),'');
    aetp=@(x,m)assertExceptionThrown(@()cosmo_index_unique(x),m);
    aet({[1,2],[3]});
    aet({[1,2],'a'});

    if cosmo_wtf('is_matlab')
        v=cosmo_wtf('version');
        is_prior_to_2012b=str2num(v(1))<=7;

        if is_prior_to_2012b
            id_different_classes='MATLAB:CELL:UNIQUE:InputClass';
        else
            id_different_classes='MATLAB:UNIQUE:InputClass';
        end
    else
        id_different_classes='';
    end
    aetp({[1,2],{1,2}},id_different_classes);
    aet(ones([2 2 2]));

function test_unique_with_nans
    n_rows_half=100+ceil(rand()*20);
    %n_rows_half=2;
    n_rows=n_rows_half*2;
    x=ceil(rand(n_rows,2)*sqrt(n_rows_half));

    rp=get_non_identity_randperm(n_rows);

    rp1=rp(1:n_rows_half);
    rp2=rp(n_rows_half+(1:n_rows_half));

    % first column, half of the rows become NaN
    x(rp1,1)=NaN;
    x(rp1,2)=1:n_rows_half;
    x(rp2,1)=1:n_rows_half;
    x(rp2,2)=1:n_rows_half;

    rp_y=get_non_identity_randperm(n_rows);

    y=x(rp_y,:);

    % test with numeric input
    [x_idx,x_unq]=cosmo_index_unique(x);
    [y_idx,y_unq]=cosmo_index_unique(y);

    assertFalse(isequal(x_unq,y_unq));
    % due to NaNs, the indices must be different
    assertFalse(isequal([x_idx{:}],rp_y([y_idx{:}])));

    % test with cell input
    [xx_idx,xx_unq]=cosmo_index_unique({x(:,1),x(:,2)});
    [yy_idx,yy_unq]=cosmo_index_unique({y(:,1),y(:,2)});

    assertFalse(isequal(xx_unq,yy_unq));
    % due to NaNs, the indices must be different
    assertFalse(isequal([xx_idx{:}],rp_y([yy_idx{:}])));

function test_index_unique_empty()
    empty_args={{},{[],{}},{[]},{{},{}}};
    for k=1:numel(empty_args)
        empty_arg=empty_args{k};

        [idxs,unq]=cosmo_index_unique(empty_arg);
        assertTrue(iscell(idxs));
        assertTrue(numel(idxs)==1);
        assertEqual(numel(unq),numel(empty_arg));
        assertTrue(all(cellfun(@isempty,unq)));
    end

function rp=get_non_identity_randperm(n)
    while true
        rp=randperm(n);
        if ~isequal(rp,1:n)
            break;
        end
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function title = p01_title ( )

%*****************************************************************************80
%
%% P01_TITLE returns the title of problem P01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Problem 1, Enright and Pryce #A1, exponential decay.';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [selection,value] = listdlg(varargin)
%LISTDLG  List selection dialog box.
%   [SELECTION,OK] = LISTDLG('ListString',S) creates a modal dialog box
%   which allows you to select a string or multiple strings from a list.
%   SELECTION is a vector of indices of the selected strings (length 1 in
%   the single selection mode).  This will be [] when OK is 0.  OK is 1 if
%   you push the OK button, or 0 if you push the Cancel button or close the
%   figure.
%
%   Double-clicking on an item or pressing <CR> when multiple items are
%   selected has the same effect as clicking the OK button.  Pressing <CR>
%   is the same as clicking the OK button. Pressing <ESC> is the same as
%   clicking the Cancel button.
%
%   Inputs are in parameter,value pairs:
%
%   Parameter       Description
%   'ListString'    cell array of strings for the list box.
%   'SelectionMode' string; can be 'single' or 'multiple'; defaults to
%                   'multiple'.
%   'ListSize'      minimum [width height] of listbox in pixels; defaults
%                   to [160 300]. The maximum [width height] is fixed to
%                   [800 600].
%   'InitialValue'  vector of indices of which items of the list box
%                   are initially selected; defaults to the first item.
%   'Name'          String for the figure's title; defaults to ''.
%   'PromptString'  string matrix or cell array of strings which appears
%                   as text above the list box; defaults to {}.
%   'OKString'      string for the OK button; defaults to 'OK'.
%   'CancelString'  string for the Cancel button; defaults to 'Cancel'.
%
%   A 'Select all' button is provided in the multiple selection case.
%
%   Example:
%     d = dir;
%     str = {d.name};
%     [s,v] = listdlg('PromptString','Select a file:',...
%                     'SelectionMode','single',...
%                     'ListString',str)
 %
%  See also DIALOG, ERRORDLG, HELPDLG, INPUTDLG,
%    MSGBOX, QUESTDLG, WARNDLG.

%   Copyright 1984-2005 The MathWorks, Inc.
%   $Revision: 6641 $  $Date: 2005/10/28 15:54:55 $

%   'uh'            uicontrol button height, in pixels; default = 22.
%   'fus'           frame/uicontrol spacing, in pixels; default = 8.
%   'ffs'           frame/figure spacing, in pixels; default = 8.

% simple test:
%
% d = dir; [s,v] = listdlg('PromptString','Select a file:','ListString',{d.name});
%

figname = '';
smode = 2;   % (multiple)
promptstring = {};
liststring = [];
listsize = [160 300];
initialvalue = [];
okstring = 'OK';
cancelstring = 'Cancel';
fus = 8;
ffs = 8;
uh = 22;

if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    cfg_message('MATLAB:listdlg:InvalidArgument', 'Arguments to LISTDLG must come param/value in pairs.')
end
for i=1:2:length(varargin)
    switch lower(varargin{i})
     case 'name'
      figname = varargin{i+1};
     case 'promptstring'
      promptstring = varargin{i+1};
     case 'selectionmode'
      switch lower(varargin{i+1})
       case 'single'
        smode = 1;
       case 'multiple'
        smode = 2;
      end
     case 'listsize'
      listsize = varargin{i+1};
     case 'liststring'
      liststring = varargin{i+1};
     case 'initialvalue'
      initialvalue = varargin{i+1};
     case 'uh'
      uh = varargin{i+1};
     case 'fus'
      fus = varargin{i+1};
     case 'ffs'
      ffs = varargin{i+1};
     case 'okstring'
      okstring = varargin{i+1};
     case 'cancelstring'
      cancelstring = varargin{i+1};
     otherwise
      cfg_message('MATLAB:listdlg:UnknownParameter', ['Unknown parameter name passed to LISTDLG.  Name was ' varargin{i}])
    end
end

if ischar(promptstring)
    promptstring = cellstr(promptstring);
end

if isempty(initialvalue)
    initialvalue = 1;
end

if isempty(liststring)
    cfg_message('MATLAB:listdlg:NeedParameter', 'ListString parameter is required.')
end

liststring=cellstr(liststring);
lfont = cfg_get_defaults('cfg_ui.lfont');
bfont = cfg_get_defaults('cfg_ui.bfont');
tmpObj = uicontrol('style','listbox',...
    'max',100,...
    'Visible','off',...
    lfont{:});
lext = cfg_maxextent(tmpObj, liststring);
ex = get(tmpObj,'Extent');
ex = ex(4);
delete(tmpObj);
listsize = min([800 600],max(listsize, [max(lext)+16, ex*numel(liststring)]));

fp = get(0,'defaultfigureposition');
w = 2*(fus+ffs)+listsize(1);
h = 2*ffs+6*fus+ex*length(promptstring)+listsize(2)+uh+(smode==2)*(fus+uh);
fp = [fp(1) fp(2)+fp(4)-h w h];  % keep upper left corner fixed

fig_props = { ...
    'name'                   figname ...
    'color'                  get(0,'defaultUicontrolBackgroundColor') ...
    'resize'                 'off' ...
    'numbertitle'            'off' ...
    'menubar'                'none' ...
    'windowstyle'            'modal' ...
    'visible'                'off' ...
    'createfcn'              ''    ...
    'position'               fp   ...
    'closerequestfcn'        'delete(gcbf)' ...
            };

fig = figure(fig_props{:});

if ~isempty(promptstring)
    prompt_text = uicontrol('style','text','string',promptstring,...
        'horizontalalignment','left',...
        'position',[ffs+fus fp(4)-(ffs+fus+ex*length(promptstring)) ...
        listsize(1) ex*length(promptstring)]); %#ok
end

btn_wid = (fp(3)-2*(ffs+fus)-fus)/2;

listbox = uicontrol('style','listbox',...
                    'position',[ffs+fus ffs+uh+4*fus+(smode==2)*(fus+uh) listsize],...
                    'string',liststring,...
                    'backgroundcolor','w',...
                    'max',smode,...
                    'tag','listbox',...
                    'value',initialvalue, ...
                    'callback', {@doListboxClick}, ...
                    lfont{:});

ok_btn = uicontrol('style','pushbutton',...
                   'string',okstring,...
                   'position',[ffs+fus ffs+fus btn_wid uh],...
                   'callback',{@doOK,listbox},...
                   bfont{:});

cancel_btn = uicontrol('style','pushbutton',...
                       'string',cancelstring,...
                       'position',[ffs+2*fus+btn_wid ffs+fus btn_wid uh],...
                       'callback',{@doCancel,listbox},...
                       bfont{:});

if smode == 2
    selectall_btn = uicontrol('style','pushbutton',...
                              'string','Select all',...
                              'position',[ffs+fus 4*fus+ffs+uh listsize(1) uh],...
                              'tag','selectall_btn',...
                              'callback',{@doSelectAll, listbox});

    if length(initialvalue) == length(liststring)
        set(selectall_btn,'enable','off')
    end
    set(listbox,'callback',{@doListboxClick, selectall_btn})
end

set([fig, ok_btn, cancel_btn, listbox], 'keypressfcn', {@doKeypress, listbox});

set(fig,'position',getnicedialoglocation(fp, get(fig,'Units')));
% Make ok_btn the default button.
setdefaultbutton(fig, ok_btn);

% make sure we are on screen
movegui(fig)
set(fig, 'visible','on'); drawnow;

try
    % Give default focus to the listbox *after* the figure is made visible
    uicontrol(listbox);
    uiwait(fig);
catch
    if ishandle(fig)
        delete(fig)
    end
end

if isappdata(0,'ListDialogAppData__')
    ad = getappdata(0,'ListDialogAppData__');
    selection = ad.selection;
    value = ad.value;
    rmappdata(0,'ListDialogAppData__')
else
    % figure was deleted
    selection = [];
    value = 0;
end

%% figure, OK and Cancel KeyPressFcn
function doKeypress(src, evd, listbox) %#ok
switch evd.Key
 case 'escape'
  doCancel([],[],listbox);
end

%% OK callback
function doOK(ok_btn, evd, listbox) %#ok
if (~isappdata(0, 'ListDialogAppData__'))
    ad.value = 1;
    ad.selection = get(listbox,'value');
    setappdata(0,'ListDialogAppData__',ad);
    delete(gcbf);
end

%% Cancel callback
function doCancel(cancel_btn, evd, listbox) %#ok
ad.value = 0;
ad.selection = [];
setappdata(0,'ListDialogAppData__',ad)
delete(gcbf);

%% SelectAll callback
function doSelectAll(selectall_btn, evd, listbox) %#ok
set(selectall_btn,'enable','off')
set(listbox,'value',1:length(get(listbox,'string')));

%% Listbox callback
function doListboxClick(listbox, evd, selectall_btn) %#ok
% if this is a doubleclick, doOK
if strcmp(get(gcbf,'SelectionType'),'open')
    doOK([],[],listbox);
elseif nargin == 3
    if length(get(listbox,'string'))==length(get(listbox,'value'))
        set(selectall_btn,'enable','off')
    else
        set(selectall_btn,'enable','on')
    end
end

author: spm


text:
% ---------------------------------------------------------------
% |Jelle C.J. Veraa, Stock Bar chart plotter function, 3/28/2009 |
% |inspired by and some code from:                               |
% |Nagi Hatoum, Candlestick chart ploter (sic) function, 5/29/03 |
% ---------------------------------------------------------------
%
% Important: date plotting subroutine needs 'tlabel' by Carlos Adrian
%   Vargas Aguilera
%   Available at the MATLAB Central File Exchange site, File ID: #19314
%   http://www.mathworks.com/matlabcentral/fileexchange/19314
%
% function barChartPlot(O,H,L,C)
%          barChartPlot(O,H,L,C,date)
%          barChartPlot(O,H,L,C,date,colorLine)
%          barChartPlot([O H L C])
%          barChartPlot([O H L C],date)
%          barChartPlot([O H L C],date,colorLine)
%
% To use colors without a date, put '0' as date
%     i.e. barChartPlot(OHLC,0,'g');
%
% Required inputs: column vectors of O(pen), H(igh), L(ow)
% and C(lose) prices of commodity.
% Alternative: [rowsx4] OHLC matrix
%
% optional inputs [default]:
% - date: serial date number (make with 'datenum') [no dates, index#]
% - colorLine: Color for lines                     ['k']
%
% Note: identical inputs as required for cndlV2 except for colors

function barChartPlot(varargin)
% See if we have [OHLC] or seperate vectors and retrieve our
% required variables (Feel free to make this code more pretty ;-)
isMat = size(varargin{1},2);
indexShift = 0;
useDate = 0;

if isMat == 4,
    O = varargin{1}(:,1);
    H = varargin{1}(:,2);
    L = varargin{1}(:,3);
    C = varargin{1}(:,4);
else
    O = varargin{1};
    H = varargin{2};
    L = varargin{3};
    C = varargin{4};
    indexShift = 3;
end

if nargin+isMat < 7,
    colorLine = 'k';
else
    colorLine = varargin{3+indexShift};
end

if nargin+isMat < 6
    date = (1:length(O))';
else
    if varargin{2+indexShift} ~= 0
        date = varargin{2+indexShift};
        useDate = 1;
    else
        date = (1:length(O))';
    end
end

% Get the length of the data
l=length(O);
% Start plotting
hold on
% Draw line from Low to High%%%%%%%%%%%%%%%%%
for i=1:l
   line([date(i) date(i)], [L(i) H(i)],'Color',colorLine)
end

% w = Width of body, change multiplier to draw body thicker or thinner
% the 'min' ensures no errors on weekends ('time gap Fri. Mon.' > wanted
% spacing)
w=.4*min([(date(2)-date(1)) (date(3)-date(2))]);
% Draw open/close marks
 for i=1:l
    xOpen=[date(i)-w date(i)];
    xClose=[date(i)+w date(i)];
    yOpen=[O(i) O(i)];
    yClose=[C(i) C(i)];
    line(xOpen,yOpen,'Color',colorLine);
    line(xClose,yClose,'Color',colorLine);
 end

%%%%%%%%%% if dates supplied, put them on the x axis%%%%%%%%%%
if (nargin+isMat > 5) && useDate,
    tlabel('x');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hold off


text:
function c = integrate(c)

   %INTEGRATE   integrates each trace.
   %
   % C = INTEGRATE(C) integrates each trace. Prior to integrating, the trend
   % and mean of each trace is removed. Depending on the application, it is
   % recommended that users also consider applying a gentle high pass filter
   % to remove unwanted low frequencies that can dominate the integrated
   % waveforms.

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   c.traces = detrend(c.traces);
   c.traces = demean(c.traces);
   c.traces = integrate(c.traces);
end

author: geoscience-community-codes

repo: GISMO


text:
function [modelstring, typeflag, order] = modeltype(ar,ma)

%MODELTYPE provides a description for an ARMA-model
%   [modelstring, typeflag, order] = modeltype(AR,MA)
%
%   modelstring: AR, MA, ARMA
%   typeflag:    1   2   3
%   order: AR(p): p, MA(q): q, ARMA(r,r-1): r, NaN otherwise.
%   (ARMA setting for 'order' only for use with ARMAsel)
%
%   Example:
%   modeltype([1 .2 .5],[1 .3])
%   modeltype: ARMA(2,1)
%
%   See also: ARMAsel.

%S. de Waele, MARCH 2001

k = length(ar)-1; l = length(ma)-1;
if ~l,
   modelstring = ['AR(' int2str(k) ')']; typeflag = 1; order = k;
elseif ~k
   modelstring = ['  MA(' int2str(l) ')']; typeflag = 2; order = l;
else
   modelstring = ['ARMA(' int2str(k) ',' int2str(l) ')']; typeflag = 3;
   if l == k-1,
       order = k;
   else
       order = NaN;
   end
end

if ~nargout,
   disp(['modeltype: ' modelstring])
   clear modelstring typeflag
end


text:
classdef CoreTest < handle
    % CoreTest Base class for MTest unit tests
    %
    % Unit tests should inherit from this class
    %
    %     Licence
    %     -------
    %     Part of CoreMat. https://github.com/tomdoel/coremat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
    end

    methods
        function Assert(obj, condition, message)
            if ~condition
                CoreErrorUtilities.ThrowException('CoreTest:AssertonFailure', ...
                    ['Assertion failure: ' message]);

            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function extractSTPerframeFtrs(outdir,ftrs,stationary,flowname,params)

if ~exist(outdir,'dir')
  mkdir( outdir);
end

ff = fields(ftrs);

% Initialize the struct for features of all the frames
for fnum = 1:numel(ff)
  curf = ff{fnum};
  if strcmp(curf,'hogftrs')
    pfname = 'hf';
  elseif strcmp(curf,'flowftrs') && stationary
    pfname = [flowname 's'];
  elseif strcmp(curf,'flowftrs') && ~stationary
    pfname = flowname;
  else
    error('Unknown feature type');
  end
  ny = size(ftrs.(curf){1},1);
  nx = size(ftrs.(curf){1},2);
  nb = size(ftrs.(curf){1},3);
  psz = params.psize;
  for yy = 1:ny
    for xx = 1:nx
      for oo = 1:nb
        perframe = struct;
        perframe.units.num = {};
        perframe.units.den = {'s'};
        for fly = 1:numel(ftrs.(ff{fnum}))
          tt = ftrs.(ff{fnum}){fly}(yy,xx,oo,:);
          perframe.data{fly} = tt(:);
        end
        perframe_name = sprintf('st_%s_%02d_%02d_%d_ny%d_nx%d_ns%d',pfname,yy,xx,oo,ny,nx,psz);
        perframe.params = params;
        outfilename = fullfile(outdir,perframe_name);
        save(outfilename,'-struct','perframe','-v7.3');
      end
    end
  end
end

repo: JAABA


text:
function [f, lineSegs, theta] = fov(A, pref)
%FOV   Field of values (numerical range) of matrix A.
%   F = FOV(A), where A is a square matrix, returns a CHEBFUN F with domain [0
%   2*pi]. The image F([0 pi]) is a curve describing the extreme points of the
%   boundary of the field of values A, a convex region in the complex plane.
%
%   For a generic matrix, the boundary of the field of values is smooth and all
%   boundary points are extreme points. If A is normal, the field of values is
%   the convex hull of the eigenvalues of A, so the extreme points consist only
%   of the eigenvalues and hence F has one constant piece for each eigenvalue,
%   so it is not continuous on [0, 2*pi].
%
%   The numerical abscissa of A is equal to max(real(F)), though this is much
%   better computed as max(real(eig(A + A')))/2.
%
%   The algorithm used is that of C. R. Johnson, Numerical determination of the
%   field of values of a general complex matrix, SIAM J. Numer. Anal. 15 (1978),
%   595-602.
%
%   F = FOV(A, PREF) allows the preferences in the CHEBFUNPREF structure PREF to
%   be used in constructing F. Note that PREF.splitting will always be set to
%   TRUE by FOV and the domain will always be [0, 2*pi].
%
%   [F, LINESEGS, THETA] = FOV(A) also returns a quasimatrix LINESEGS whose
%   columns are CHEBFUN objects defining the line segments connecting up the
%   extreme points and a vector THETA specifying the values in [0, 2*pi] where
%   the discontinuities in F occur.
%
% Example 1 (smooth boundary)
%   A = randn(5);
%   F = fov(A);
%   hold off, plot(F, '-b')
%   e = eig(A);
%   hold on, plot(e, '*k'), hold off, axis equal
%
% Example 2 (boundary has a corner)
%   A = [0 1 0 ; 0 0 0 ; 0 0 1];
%   [F, lineSegs] = fov(A);
%   plot(F, '-b'), hold on
%   plot(lineSegs, '-r')
%   e = eig(A);
%   plot(complex(e), '*k'), hold off, axis equal

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

if ( nargin < 2 )
    % Obtain preferences:
    pref = chebfunpref();
end

% Set splitting to 'on' and the domain to [0, 2*pi].
pref.splitting = true;
pref.domain = [0, 2*pi];

% Construct a CHEBFUN of the FOV curve, and try to merge any unnecessary
% breakpoints out of the result:
f = chebfun(@(theta) fovCurve(theta, A), pref);
f = merge(f);

if ( nargout == 1 )
    return
end

%% LINESEGS - implemented by Michael Overton
% Look for possible discontinuities at any break points. The boundary of the
% field of values may have corners, but it is continuous: hence the need
% sometimes for linear interpolation connecting up the sets of extreme points.

ends = f.domain;      % Break points are between 0 and 2*pi.
delta = 1e-14;        % Must be bigger than machine precision but not much
                      % bigger: if it is too small spurious line segments will
                      % have tiny length, which hardly matters.
if ( any(diff(ends) < delta) )
    delta = min(diff(ends))/3;
end
m = length(ends) - 1; % Number of pieces

% Get the values to the left and right of break points. The first break point is
% zero, which wraps around to 2*pi and needs special treatment. Exclude the
% right end point 2*pi.
lVal = feval(f, ends(1:end-1), 'left');  % Values to left of break points.
lVal(1) = feval(f, ends(end),  'left');  % Value to left of 2*pi.
rVal = feval(f, ends(1:end-1), 'right'); % Values to right of break points.

% Determine which points correspond to discontinuities:
tol = 10*delta*max(abs([lVal, rVal]), [], 2);
discont = abs(lVal - rVal) > tol;

% Define additional CHEBFUNs for the line segments joining discontinuities.
% (First idea was to do this by inserting tiny intervals with steep slopes into
% the CHEBFUN f, but this leads to loss of accuracy.) Each column of lineSegs
% represents a line in the complex plane connecting the points lVal(j) and
% rVal(j) to each other, along with the corresponding theta value in [0,2*pi].
theta = ends(discont);
lineSegs = chebfun([lVal(discont) ; rVal(discont)], [-1 1], 'tech', @chebtech2);

end

function z = fovCurve(theta, A)
z = NaN(size(theta));
for j = 1:length(theta)
    r = exp(1i*theta(j));
    B = r*A;
    H = (B + B')/2;
    [X, D] = eig(H);
    [ignored, k] = max(diag(D));
    v = X(:,k);
    z(j) = v'*A*v/(v'*v);
end
end


text:
% SB2_CONTROLSETTINGS  Set parameters to control the SPARSEBAYES algorithm
%
% CONTROLS = SB2_CONTROLSETTINGS
%
% OUTPUT ARGUMENTS:
%
%	CONTROLS	A structure whose fields control various aspects of the
%				running of the SPARSEBAYES algorithm.
%
%	.ZeroFactor			Small number equivalent to zero for Q^2-S
%	.MinDeltaLogAlpha	Termination criterion for changes in log-alpha
%	.MinDeltaLogBeta	Termination criterion for changes in log-beta
%
%	.PriorityAddition	Prefer "addition" operations
%	.PriorityDeletion	Prefer "deletion" operations
%
%	.BetaUpdateStart	How many "fast start" beta updates
%	.BetaUpdateFrequency
%						How regularly to update beta after the above
%	.BetaMaxFactor		Minimum value control for noise estimate
%
%	.PosteriorModeFrequency
%						How regularly to re-find the posterior mode
%
%	.BasisAlignmentTest	Test for redundant basis vectors?
%	.AlignmentMax		Basis redundancy criterion
%
% NOTES:
%
% The various definitions in the file are effectively "fixed" and not
% modified elsewhere.
%
% The interested user may wish to experiment with the operation of the
% SPARSEBAYES algorithm by modifying the values the file directly. See the
% inline comments for hints on the various control settings.
%

%
% Copyright 2009, Vector Anomaly Ltd
%
% This file is part of the SPARSEBAYES library for Matlab (V2.0).
%
% SPARSEBAYES 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.
%
% SPARSEBAYES 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 SPARSEBAYES in the accompanying file "licence.txt"; if not, write to
% the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
% MA 02110-1301 USA
%
% Contact the author: m a i l [at] m i k e t i p p i n g . c o m
%
function CONTROLS = SB2_ControlSettings

%% Define parameters which influence the underlying operation of the
%% SparseBayes inference algorithm

% TOLERANCES
%
% Any Q^2-S "relevance factor" less than this is considered to be zero
%
CONTROLS.ZeroFactor			= 1e-12;
%
% If the change in log-alpha for the best re-estimation is less than this,
% we consider termination
%
CONTROLS.MinDeltaLogAlpha	= 1e-3;
%
% In the Gaussian case, we also require a beta update to change the value
% of log-beta (inverse noise variance) less than this to terminate
%
CONTROLS.MinDeltaLogBeta	= 1e-6;

% ADD/DELETE
%
% - preferring addition where possible will probably make the algorithm a
% little slower and perhaps less "greedy"
%
% - preferring deletion may make the model a little more sparse and the
% algorithm may run slightly quicker
%
% Note: both these can be set to 'true' at the same time, in which case
% both take equal priority over re-estimation.
%
CONTROLS.PriorityAddition	= false;
CONTROLS.PriorityDeletion	= true;

% (GAUSSIAN) NOISE
%
% When to update the noise estimate
%
% The number of iterations from the start for which we update it every
% iteration (to get in the right ball-park to begin with)
%
CONTROLS.BetaUpdateStart		= 10;
%
% After the above, we only regularly update it after
% a given number of iterations
%
CONTROLS.BetaUpdateFrequency	= 5;
%
% Prevent zero-noise estimate (perfect fit) problem
% -	effectively says the noise variance estimate is clamped to be no
%	lower than variance-of-targets / BetaMaxFactor.
%
CONTROLS.BetaMaxFactor			= 1e6;

% POSTERIORMODE
%
% How many alpha updates to do in between each full posterior mode
% computation in the non-Gaussian case
%
% In principle, this should be set to one (to update the posterior every
% iteration) but it may be more efficient to do several alpha updates before
% re-finding the posterior mode.
%
CONTROLS.PosteriorModeFrequency	= 1;

% REDUNDANT BASIS
%
% Check for basis vector alignment/correlation redundancy
%
CONTROLS.BasisAlignmentTest		= true;
%
ALIGNMENT_ZERO					= 1e-3;
%
% If BasisAlignmentTest is true, any basis vector with inner product more
% than MAX_ALIGNMENT with any existing model vector will not be added
%
CONTROLS.AlignmentMax			= 1 - ALIGNMENT_ZERO;

repo: REST


text:
function [Y,h]=glyph(X,r,c)
%Syntax: [Y,h]=glyph(X,r,c)
%__________________________
%
% 3D glyph visualization. The glyph is a convex deltahedron where each
% node is connected to its 6 nearest nodes.
%
% Y is the N-by-3 matrix whith the cartesian coordinates of the points on
%   the glyph.
% h returns a vector of tetrahedra handles. Each element of h is a handle
%   to the set of patches forming one tetrahedron. Type "help tetramesh"
%   for more info.
% X is the N-by-3 matrix whith the cartesian coordinates of the points on
%   the sphere.
% r is the range parameter.
% c is the color parameter.
%
%
% References:
%
% Sangole A., Knopf G. K. (2002): Representing high-dimensional data sets
% as close surfaces. Journal of Information Visualization 1: 111-119
%
% Sangole A., Knopf G. K. (2003): Geometric representations for
% high-dimensional data using a spherical SOFM. International Journal of
% Smart Engineering System Design 5: 11-20
%
% Sangole A., Knopf G. K. (2003): Visualization of random ordered numeric
% data sets using self-organized feature maps. Computers and Graphics 27:
% 963-976
%
% Sangole A. P. (2003): Data-driven Modeling using Spherical
% Self-organizing Feature Maps. Doctor of Philosophy (Ph.D.) Thesis.
% Department of Mechanical and Materials Engineering. Faculty of
% Engineering. The University of Western Ontario, London, Ontario, Canada.
%
%
% Remark:
%
% If no output is desired, the function plots the glyph.
%
%
% Archana P. Sangole, PhD., P.E. (TX chapter)
% School of Physical & Occupational Therapy
% McGill University
% 3654 Promenade Sir-William-Osler
% Montreal, PQ, H3G 1Y5
% e-mail: archana.sangole@mail.mcgill.ca
%
% CRIR, Rehabilitation Institute of Montreal
% 6300 Ave Darlington
% Montreal, PQ, H3S 2J5
% Tel: 514.340.2111 x2188
% Fax: 514.340.2154
%
%
% Alexandros Leontitsis, PhD
% Department of Education
% University of Ioannina
% 45110- Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 23-Mar-2006

% Add a center to the sphere, which is [0 0 0]
X=[X;zeros(1,3)];

% 3-dimensional Delaunay tessellation
tri=delaunay3(X(:,1),X(:,2),X(:,3));

% Remove the center
X(end,:)=[];

if nargin<2 | isempty(r)==1
    r=ones(length(X),1);
else
    % r must be a vector
    if min(size(r))>1
        error('r must be a vector.');
    end
    % The length of r should be equal to the length of X.
    if length(r)~=length(X)
        error('The length of r should be equal to the length of X.');
    end
    r=r(:);
end

if nargin<3 | isempty(c)==1
    c=r;
else
    % c must be a vector
    if min(size(c))>1
        error('c must be a vector.');
    end
    % The length of c should be equal to the length of X.
    if length(c)~=length(X)
        error('The length of c should be equal to the length of X.');
    end
    c=c(:);
end

% Go to spherical coordinates, ...
[theta,phi]=cart2sph(X(:,1),X(:,2),X(:,3));
% ... and compute the cartesian coordinates with the given r
[Y(:,1),Y(:,2),Y(:,3)]=sph2cart(theta,phi,r);

% Add a center to the glyph, which is [0 0 0]
Y=[Y;zeros(1,3)];
% Sort the vertices of each tereahedron ...
tri=sort(tri')';
% ... in order to calculate the color
i=1:length(tri);
cnew(i)=mean(c(tri(i,1:end-1))')';

% If no output is desired, plot the glyph
if nargout==0
    % Plot the glyph
    tetramesh(tri,Y,cnew,'FaceAlpha',1,'LineStyle','-');
    % Define the x axis
    xlim([-1.2 1.2])
    % Define the y axis
    ylim([-1.2 1.2])
    % Define the z axis
    zlim([-1.2 1.2])
    % Define the color axis
    caxis([1 1.2]);
    % Remove the axes
    axis off
    % Freeze aspect ratio properties to faciliate 3D rotation
    axis vis3d
% Else if the handle is desired
elseif nargout==2
    % Retrieve the handle
    h=tetramesh(tri,Y,cnew,'FaceAlpha',1,'LineStyle','none');
end


text:
function varargout = drawPolygon3d(varargin)
%DRAWPOLYGON3D Draw a 3D polygon specified by a list of vertex coords
%
%   drawPolygon3d(POLY);
%   packs coordinates in a single N-by-3 array.
%
%   drawPolygon3d(PX, PY, PZ);
%   specifies coordinates in separate numeric vectors (either row or
%   columns)
%
%   drawPolygon3d(..., PARAM, VALUE);
%   Specifies style options to draw the polyline, see plot for details.
%
%   H = drawPolygon3d(...);
%   also returns a handle to the list of created line objects.
%
%   Example
%     t = linspace(0, 2*pi, 100)';
%     xt = 10 * cos(t);
%     yt = 5 * sin(t);
%     zt = zeros(1,100);
%     figure; drawPolygon3d(xt, yt, zt, 'b');
%
%   See Also:
%   polygons3d, fillPolygon3d, drawPolyline3d
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2011-08-17 from drawPolyline3d, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% HISTORY

% check case we want to draw several curves, stored in a cell array
var = varargin{1};
if iscell(var)
    hold on;
    h = [];
    for i = 1:length(var(:))
        h = [h; drawPolygon3d(var{i}, varargin{2:end})]; %#ok<AGROW>
    end
    if nargout > 0
        varargout{1} = h;
    end
    return;
end

% extract curve coordinates
if min(size(var)) == 1
    % if first argument is a vector (either row or column), then assumes
    % first argument contains x coords, second argument contains y coords
    % and third one the z coords
    px = var;
    if length(varargin) < 3
        error('geom3d:drawPolygon3d:Wrong number of arguments in fillPolygon3d');
    end
    py = varargin{2};
    pz = varargin{3};
    varargin = varargin(4:end);
else
    % first argument contains both coordinate
    px = var(:, 1);
    py = var(:, 2);
    pz = var(:, 3);
    varargin = varargin(2:end);
end

%% draw the polygon

% check that the polygon is closed
if px(1) ~= px(end) || py(1) ~= py(end) || pz(1) ~= pz(end)
    px = [px(:); px(1)];
    py = [py(:); py(1)];
    pz = [pz(:); pz(1)];
end

% draw the closed curve
h = plot3(px, py, pz, varargin{:});

if nargout > 0
    varargout = {h};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [dataX, dataY] = pix2data( ax, pixelX, pixelY )
%PIX2DATA Translate a Figure Pixel position into an Axes data point
%  If no Axes given, it defaults to GCA.
%  [X, Y] = PIX2DATA( AXES, x, y )
%  [X, Y] = PIX2DATA( AXES, [x y] )
%  [X, Y] = PIX2DATA( x, y )
%  [X, y] = PIX2DATA( [x y] )
%
%  For scalar or vector X and Y, output will have X and Y in columns
%  POINT  = PIX2DATA( AXES, x, y )
%  POINT  = PIX2DATA( AXES, [x y] )
%  POINT  = PIX2DATA( x, y )
%  POINT  = PIX2DATA( [x y] )
%
%  See Also: DATA2PIX

% Chuck Packard, The Mathworks, Inc.
% 25 March 95

if nargin == 0
   error( 'Need AXES and X,Y position.' )
end

%See if there is a GCA, don't create one
if isempty( get(0,'children') )
   error( 'No Figure windows.' );
end
if isempty( get( gcf, 'children' ) )
   error( 'No Axes in current Figure.' );
end

if nargin == 1
   if size( ax, 2 ) ~= 2
      error( 'When defaulting to GCA, first argument should be [X Y].' )
   end
   pixelY = ax(:,2);
   pixelX = ax(:,1);
   ax = gca;

elseif nargin == 2
   if length(ax) == 1 & size(pixelX,2) == 2
      pixelY = pixelX(:,2);
      pixelX(:,2) = [];
   else
      pixelY = pixelX;
      pixelX = ax;
      ax = gca;
   end
end

%Check data, make X and Y column vectors
xSize = size(pixelX);
ySize = size(pixelY);
if any( ySize ~= xSize )
   error( 'Size of X and Y data not consistent.' );
end

%Get position of Axes in pixels
origAxesUnits = get(ax,'units');
set( ax, 'units', 'pixels' )
axPos = fix(get(ax,'position'));
set( ax, 'units', origAxesUnits )

%Find extent of Axes in data
xlimits = get(ax,'xlim');
ylimits = get(ax,'ylim');
dataWidth = xlimits(2) - xlimits(1);
dataHeight = ylimits(2) - ylimits(1);

%Find normalized position of points and determine pixel position
% as offset from lower left of Axes.
X = xlimits(1) + (((pixelX - axPos(1)) ./ axPos(3)) .* dataWidth);
Y = ylimits(1) + (((pixelY - axPos(2)) ./ axPos(4)) .* dataHeight);

if nargout < 2
   if min(ySize) ~= 1
      %What would it mean to concatenate two matrices for this function?
      dataX = X;
   elseif ySize(1) == 1 & ySize(2) > 1
      %Need to transpose data so it looks nice in column format
      dataX = [X' Y'];
   else
      dataX = [X Y];
   end
else
   dataX = X;
   dataY = Y;
end


text:
function cg_test023 ( )

%*****************************************************************************80
%
%% CG_TEST023 tests R83S_CG.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 July 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CG_TEST023\n' );
  fprintf ( 1, '  Test R83S_CG, applying CG to an R83S matrix.\n' );

  n = 10;
%
%  Let A be the -1 2 -1 matrix.
%
  a = r83s_dif2 ( n, n );
%
%  Choose a random solution.
%
  seed = 123456789;
  [ x1, seed ] = r8vec_uniform_01 ( n, seed );
%
%  Compute the corresponding right hand side.
%
  b = r83s_mv ( n, n, a, x1 );
%
%  Call the CG routine.
%
  x2 = ones ( n, 1 );
  x2 = r83s_cg ( n, a, b, x2 );
%
%  Compute the residual.
%
  r = r83s_resid ( n, n, a, x2, b );
  r_norm = r8vec_norm ( n, r );
%
%  Compute the error.
%
  e_norm = r8vec_diff_norm ( n, x1, x2 );
%
%  Report.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of variables N = %d\n', n );
  fprintf ( 1, '  Norm of residual ||Ax-b|| = %g\n', r_norm );
  fprintf ( 1, '  Norm of error ||x1-x2|| = %g\n', e_norm );

  return
end


text:
function halton_test ( )

%*****************************************************************************80
%
%% HALTON_TEST tests the HALTON library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2008
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HALTON_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the HALTON library.\n' );

  halton_test01 ( );
  halton_test0125 ( );
  halton_test0126 ();
  halton_test02 ( );
  halton_test025 ( );
  halton_test03 ( );
  halton_test04 ( );
  halton_test045 ( );
  halton_test05( );
  halton_test06 ( );
  halton_test07 ( );
  halton_test08 ( );
  halton_test09 ( );
  halton_test10 ( );

  halton_test11 ( );
  halton_test12 ( );
  halton_test13 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HALTON_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function test_values_test100 ( )

%*****************************************************************************80
%
%% TEST100 demonstrates the use of LAPLACE_CDF_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST100:\n' );
  fprintf ( 1, '  LAPLACE_CDF_VALUES stores values of\n' );
  fprintf ( 1, '  the Laplace Cumulative Density Function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      Mu        Beta          X                  CDF(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, mu, beta, x, fx ] = laplace_cdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %12f  %12f  %24.16f\n', mu, beta, x, fx );

  end

  return
end


text:
function T = findFutileCycle(model, cut,closedModel)
% This function attempts to find reactions involved in the atp-driving
% futile cycle
%
% INPUT
% model         model structure
% cut           cutoff value for reactions to be displayed
% closedModel   if 1 the model will be closed otw the applied medium
%               constraints count (DEFAULT:1);
% OUTPUT
% T         Table of reactions potentially involved
%
% Ines Thiele 03/2022

if ~exist('cut','var')
    cut = 250;
end

if ~exist('closedModel','var')
    closedModel = 1;
end

modelClosed = model;
if closedModel
    modelexchanges1 = strmatch('Ex_',modelClosed.rxns);
    modelexchanges4 = strmatch('EX_',modelClosed.rxns);
    modelexchanges2 = strmatch('DM_',modelClosed.rxns);
    modelexchanges3 = strmatch('sink_',modelClosed.rxns);
    selExc = (find( full((sum(abs(modelClosed.S)==1,1) ==1) & (sum(modelClosed.S~=0) == 1))))';

    modelexchanges = unique([modelexchanges1;modelexchanges2;modelexchanges3;modelexchanges4;selExc]);
    modelClosed.lb(ismember(modelClosed.rxns,modelClosed.rxns(modelexchanges)))=0;
end
modelClosed = changeObjective(modelClosed,'DM_atp_c_');

fba = optimizeCbModel(modelClosed,'max','zero')
fba.v(find(abs(fba.v)<=1e-6))=0;
tab = [modelClosed.lb modelClosed.ub fba.v];
a = printRxnFormula(model);

R = modelClosed.rxns(find(abs(fba.v)>cut));
L = modelClosed.lb(find(abs(fba.v)>cut));
U = modelClosed.ub(find(abs(fba.v)>cut));
F = fba.v(find(abs(fba.v)>cut));
A = a(find(abs(fba.v)>cut));
if ~isempty(A)
T = table(R, num2str(L),num2str(U), num2str(F),A);
else
    T = {};
end

author: opencobra

repo: cobratoolbox


text:
% read in a sample image -- also see letters.png, bagel.png
img = imread('mushroom.png');
imshow(img);

% the standard skeletonization:
imshow(bwmorph(img,'skel',inf));

% the new method:
imshow(bwmorph(skeleton(img)>35,'skel',Inf));

% in more detail:
[skr,rad] = skeleton(img);

% the intensity at each point is proportional to the degree of evidence
% that this should be a point on the skeleton:
imagesc(skr);
colormap jet
axis image off

% skeleton can also return a map of the radius of the largest circle that
% fits within the foreground at each point:
imagesc(rad)
colormap jet
axis image off

% thresholding the skeleton can return skeletons of thickness 2,
% so the call to bwmorph completes the thinning to single-pixel width.
skel = bwmorph(skr > 35,'skel',inf);
imshow(skel)
% try different thresholds besides 35 to see the effects

% anaskel returns the locations of endpoints and junction points
[dmap,exy,jxy] = anaskel(skel);
hold on
plot(exy(1,:),exy(2,:),'go')
plot(jxy(1,:),jxy(2,:),'ro')


text:
function y = A_fhp_rect(x, OMEGA, m, n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Shiqian Ma
% Date : 09/05/2007
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% yc = 1/n*fft2(reshape(x,m,n));

yc = 1/sqrt(m*n)*fft2(reshape(x,m,n));
%%%% y = [yc(1,1); sqrt(2)*real(yc(OMEGA)); sqrt(2)*imag(yc(OMEGA))];
y = [yc(1,1); real(yc(OMEGA)); imag(yc(OMEGA))];

author: thomaskuestner


text:
% Re-select te corners after calibration

if ~exist('n_ima')|~exist('fc'),
   fprintf(1,'No calibration data available.\n');
   return;
end;

check_active_images;

flag = 0;
for kk = ind_active,
   if ~exist(['y_' num2str(kk)]),
      flag = 1;
   else
      eval(['ykk = y_' num2str(kk) ';']);
      if isnan(ykk(1,1)),
	 flag = 1;
      end;
   end;
end;

if flag,
   fprintf(1,'Need to calibrate once before before recomputing image corners. Maybe need to load Calib_Results.mat file.\n');
   return;
end;

if n_ima == 0,
    fprintf(1,'No image data available\n');
    return;
end

%if ~exist(['I_' num2str(ind_active(1))]),
%   n_ima_save = n_ima;
%   active_images_save = active_images;
%   ima_read_calib;
%   n_ima = n_ima_save;
%   active_images = active_images_save;
%   check_active_images;
%   if no_image_file,
%      disp('Cannot extract corners without images');
%      return;
%   end;
%end;

fprintf(1,'\nRe-extraction of the grid corners on the images (after first calibration)\n');

if ~dont_ask,
    disp('Window size for corner finder (wintx and winty):');
    wintx = input('wintx ([] = 20) = ');
    if isempty(wintx), wintx = 20; end;
    wintx = round(wintx);
    winty = input('winty ([] = 20) = ');
    if isempty(winty), winty = 20; end;
    winty = round(winty);
else
    wintx = 20;
    winty = 20;
end;

fprintf(1,'Window size = %dx%d\n',2*wintx+1,2*winty+1);

if ~dont_ask,
    ima_numbers = input('Number(s) of image(s) to process ([] = all images) = ');
else
    ima_numbers = [];
end;

if isempty(ima_numbers),
   ima_proc = 1:n_ima;
else
   ima_proc = ima_numbers;
end;

if ~dont_ask,
    q_auto = input('Use the projection of 3D grid or manual click ([]=auto, other=manual): ','s');
else
    q_auto = [];
end;

fprintf(1,'Processing image ');

for kk = ima_proc;
    if active_images(kk),
        fprintf(1,'%d...',kk);

        if ~type_numbering,
            number_ext =  num2str(image_numbers(kk));
        else
            number_ext = sprintf(['%.' num2str(N_slots) 'd'],image_numbers(kk));
        end;
        ima_name = [calib_name  number_ext '.' format_image];

        if exist(ima_name),
            if format_image(1) == 'p';
                if format_image(2) == 'p',
                    I = double(loadppm(ima_name));
                else
                    I = double(loadpgm(ima_name));
                end;
            else
                if format_image(1) == 'r',
                    I = readras(ima_name);
                else
                    I = double(imread(ima_name));
                end;
            end;
            if size(I,3)>1,
                I = 0.299 * I(:,:,1) + 0.5870 * I(:,:,2) + 0.114 * I(:,:,3);
            end;
            [ny,nx,junk] = size(I);
            if isempty(q_auto),
                eval(['y = y_' num2str(kk) ';']);
                xc = cornerfinder(y+1,I,winty,wintx); % the four corners
                eval(['wintx_' num2str(kk) ' = wintx;']);
                eval(['winty_' num2str(kk) ' = winty;']);
                eval(['x_' num2str(kk) '= xc - 1;']);
            else
                fprintf(1,'\n');
                fprintf(1,'\nProcessing image %d...\n',kk);
                click_calib_fisheye_no_read;
            end;
        else
            fprintf(1,'Image %s not found!!!...',ima_name);
        end;
    else
        if ~exist(['omc_' num2str(kk)]),
            eval(['dX_' num2str(kk) ' = NaN;']);
            eval(['dY_' num2str(kk) ' = NaN;']);
            eval(['wintx_' num2str(kk) ' = NaN;']);
            eval(['winty_' num2str(kk) ' = NaN;']);
            eval(['x_' num2str(kk) ' = NaN*ones(2,1);']);
            eval(['X_' num2str(kk) ' = NaN*ones(3,1);']);
            eval(['n_sq_x_' num2str(kk) ' = NaN;']);
            eval(['n_sq_y_' num2str(kk) ' = NaN;']);
        end;
    end;
end;

% Recompute the error:
comp_error_calib_fisheye;
fprintf(1,'\ndone\n');

repo: ekfmonoslam


text:
function [CODE, Declarations, NewVars] = ParseTmpFile(FILENAME,Result_Name)

%This function is designed to read the output of ccode.m back into matlab
%and parse it into a more useful format. It assumes a very specific format
%for the file, and assumes that the last line is a result, and
%automatically assigns it Result_Num.

%Figure out what to call things
persistent Result_Var_Num
persistent Local_Var_Num

%Initialize the persistent variables
if isempty(Result_Var_Num), Result_Var_Num=0; end
if isempty(Local_Var_Num), Local_Var_Num=2; end

%What should local variables be called?
BaseName = 'tmp';

%Check what the result name should be
if nargin == 1
    Result_Name = ['Result_' num2str(Result_Var_Num)];
    Result_Var_Num = Result_Var_Num + 1;  %Increment the number
end

% First, open the file and count the number of lines
fid = fopen(FILENAME);
N_Lines = 0;
while 1
    tline = fgetl(fid);
    if ~ischar(tline), break, end
    N_Lines = N_Lines + 1;
%     disp(tline)
end
% disp(N_Lines)
fclose(fid);

%Next, initialize the outputs;
CODE = cell(N_Lines,1);
Variables = cell(N_Lines,1);
Declarations = cell(2,1);
Declarations{1} = 'double  ';  %initialize the declarations for the tmp variables
Declarations{2} = 'double  ';  %initialize the declaration for the results

%Now, do the useful things:
fid = fopen(FILENAME);
for i=1:N_Lines
    tline = fgetl(fid);
    %Figure out what to name things:
    if i==N_Lines
        %This is the last line of the file -> it stores the result
        VarName = Result_Name;
        index = strfind(tline,'=');  %Search for the equal sign
        CODE{i} = [VarName ' = ' tline((index+1):end)];
        Declarations{2} = [Declarations{2} VarName];   %declare the result
    else
        %This is a regular line, so get the name from the file
        index = strfind(tline,'=');  %Search for the equal sign
        VarName = tline(1:(index-1)); %Get the first part of the string
        VarName = VarName(~isspace(VarName));  %Remove the spaces
        CODE{i} = tline;
        if i==1
            %Then first line, don't need the comma
            Declarations{1} = [Declarations{1}  VarName ];
        else
        Declarations{1} = [Declarations{1} ', ' VarName];   %Add it to the tmp pile
        end
    end

    Variables{i} = VarName;
end

Declarations{1} = [Declarations{1} ';'];
Declarations{2} = [Declarations{2} ';'];

fclose(fid);  %Close the file

%% Rename local variables
% This function will be called multiple times, and the output stored in a
% single file for all function calls. Thus, if the same variable names are
% used (such as t2,t3,t4...) then they will cause a compile error. To fix
% this, the function uses a persistent variable to continue increasing the
% local variable ID numbers.

NewVars = Variables;   %Store the new variable names here
for i=1:(length(Variables)-1)
    NewVars{i} = [BaseName num2str(Local_Var_Num)];  %Create the new name
    Local_Var_Num = Local_Var_Num + 1;   %Increment the counter
end

%Now do a find and replace operation on all of the functions:
for i=1:length(CODE)
    for j = 1:length(NewVars)
        CODE{i} = strrep(CODE{i},Variables{j},NewVars{j});
    end
end

for i = 1:length(Declarations)
    for j = 1:length(NewVars)
        Declarations{i} = strrep(Declarations{i},Variables{j},NewVars{j});
    end
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function g = noiseGradientParam(noise, mu, varsigma, y)

% NOISEGRADIENTPARAM Gradient wrt the noise model's parameters.
% FORMAT
% DESC computes the gradient of the log Z of the given noise model
% with respect to the of functions with respect to the given
% noise model's parameters.
% ARG noise : the noise structure for which the gradients are being
% computed.
% ARG mu : the input means for which the gradients are being computed.
% ARG varSigma : the input variances for which the gradients are being computed.
% ARG y : the target values for the noise model.
% RETURN g : gradients of the log probability with respect to
% the noise parameters. The ordering of the vector should match
% that provided by the function noiseExtractParam.
%
%
% SEEALSO noiseCreate, noiseParamInit, noiseGradVals, noiseGradientParam
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

fhandle = str2func([noise.type 'NoiseGradientParam']);
g = fhandle(noise, mu, varsigma, y);

% check if there is a prior over parameters
if isfield(noise, 'priors')
  for i = 1:length(noise.priors)
    index = noise.priors(i).index;
    g(index) = g(index) + priorGradient(noise.priors(i), params(index));
  end
end

% Check if parameters are being optimised in a transformed space.
if isfield(noise, 'transforms')
  fhandle = str2func([noise.type 'NoiseExtractParam']);
  params = fhandle(noise);
  for i = 1:length(noise.transforms)
    index = noise.transforms(i).index;
    fhandle = str2func([noise.transforms(i).type 'Transform']);
    g(index) = g(index).*fhandle(params(index), 'gradfact');
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef Optimizer_PG < Optimizer_Unconstrained

    properties  (GetAccess = public, SetAccess = private)
        optimality_tol
    end

    properties (GetAccess = public, SetAccess = protected)
        type = 'PROJECTED GRADIENT'
    end

    properties (Access = private)
       upperBound
       lowerBound
    end

    methods (Access = public)

        function obj = Optimizer_PG(cParams)
            obj@Optimizer_Unconstrained(cParams);
            obj.upperBound = cParams.ub;
            obj.lowerBound = cParams.lb;
        end

        function compute(obj)
            x_n      = obj.designVariable.value;
            gradient = obj.objectiveFunction.gradient;
            x_new = x_n-obj.lineSearch.value*gradient;
            ub = obj.upperBound*ones(length(x_n(:,1)),1);
            x2 = reshape(x_n,[],2); norm(x2(:,1)-x2(:,2))
            lb = obj.lowerBound*ones(length(x_n(:,1)),1);
            x_new = max(min(x_new,ub),lb);
            obj.designVariable.update(x_new);
            l2Norm = obj.designVariable.computeL2normIncrement();
            obj.optimalityCond = sqrt(l2Norm);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function noOffset = removeOffset(readings,refFreq)
%this function removes the offset from frequency readings by calculating
%the mean, subtracting the ref freq by the mean, and then adding the diff
%to the array of readings

offset = refFreq - mean(readings);
noOffset = readings + offset;

repo: mcbench-benchmarks


text:
function det = r8ge_np_det ( n, a_lu )

%*****************************************************************************80
%
%% R8GE_NP_DET computes the determinant of a matrix factored by R8GE_NP_FA.
%
%  Discussion:
%
%    The R8GE storage format is used for a general M by N matrix.  A storage
%    space is made for each logical entry.  The two dimensional logical
%    array is mapped to a vector, in which storage is by columns.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, real A_LU(N,N), the LU factors from R8GE_NP_FA.
%
%    Output, real DET, the determinant of the matrix.
%
  det = 1.0E+00;

  for i = 1 : n
    det = det * a_lu(i,i);
  end

  return
end


text:
function ME = myfunc(x)
% Break Test

for i = 1:100
    try
        x = x - 1;
        if x < 0
            error('x has become negative')
        end
    catch ME
        fprintf('x become negative')
        break
    end
end

repo: Aero-matlab


text:
function res = le(a,b)
%LE           Implements  a <= b  for Taylor, compares only a.x and b.x
%

% written  05/21/09     S.M. Rump
%

  if ~isa(a,'taylor')
    res = ( a<=reshape(b.t(1,:),b.size) );
  elseif ~isa(b,'taylor')
    res = ( reshape(a.t(1,:),a.size)<=b );
  else
    res = ( reshape(a.t(1,:),a.size)<=reshape(b.t(1,:),b.size) );
  end

repo: OpenMAS


text:
% --------------------------------------------------------
% MDP Tracking
% Copyright (c) 2015 CVGL Stanford
% Licensed under The MIT License [see LICENSE for details]
% Written by Yu Xiang
% --------------------------------------------------------
%
% offline training using ground truth annotations
function tracker = MDP_train_gt(seq_idx, tracker)

is_show = 0;
is_save = 1;
is_text = 0;
is_pause = 0;

opt = globals();
opt.is_show = is_show;

seq_name = opt.mot2d_train_seqs{seq_idx};
seq_num = opt.mot2d_train_nums(seq_idx);
seq_set = 'train';

if is_show
    close all;
end

% build the dres structure for images
filename = sprintf('%s/%s_dres_image.mat', opt.results, seq_name);
if exist(filename, 'file') ~= 0
    object = load(filename);
    dres_image = object.dres_image;
    fprintf('load images from file %s done\n', filename);
else
    dres_image = read_dres_image(opt, seq_set, seq_name, seq_num);
    fprintf('read images done\n');
    save(filename, 'dres_image', '-v7.3');
end

% read ground truth
filename = fullfile(opt.mot, opt.mot2d, seq_set, seq_name, 'gt', 'gt.txt');
dres_gt_all = read_mot2dres(filename);
dres_gt_all = fix_groundtruth(seq_name, dres_gt_all);

% generate training data
I = dres_image.Igray{1};
[dres_train, dres_det, labels] = generate_training_data(seq_idx, dres_image, opt);

% for debugging
% dres_train = {dres_train{6}};

% intialize tracker
if nargin < 2 || isempty(tracker) == 1
    fprintf('initialize tracker from scratch\n');
    tracker = MDP_initialize(I, dres_det, labels, opt);
else
    % continuous training
    fprintf('continuous training\n');
    tracker.image_width = size(I,2);
    tracker.image_height = size(I,1);
    tracker.max_width = max(dres_det.w);
    tracker.max_height = max(dres_det.h);
    tracker.max_score = max(dres_det.r);

    factive = MDP_feature_active(tracker, dres_det);
    index = labels ~= 0;
    tracker.factive = [tracker.factive; factive(index,:)];
    tracker.lactive = [tracker.lactive; labels(index)];
    tracker.w_active = svmtrain(tracker.lactive, tracker.factive, '-c 1 -q');
end

% for each training sequence
num_train = numel(dres_train);
for t = 1:num_train

    fprintf('tracking sequence %d\n', t);

    dres_gt = dres_train{t};

    % first frame
    fr = dres_gt.fr(1);
    id = dres_gt.id(1);

    % reset tracker
    tracker.prev_state = 1;
    tracker.state = 1;
    tracker.target_id = id;

    % start tracking
    while fr <= dres_gt.fr(end)
        if is_text
            fprintf('\nframe %d, state %d\n', fr, tracker.state);
        end

        % show results
        if is_show
            figure(1);
            % show ground truth
            subplot(2, 3, 1);
            show_dres(fr, dres_image.I{fr}, 'GT', dres_gt);

            % show detections
            subplot(2, 3, 2);
            show_dres(fr, dres_image.I{fr}, 'Detections', dres_det);
        end

        if tracker.state == 0
            break;

        elseif tracker.state == 1

            % initialize the LK tracker with gt
            tracker = LK_initialize(tracker, fr, id, dres_gt, 1, dres_image);
            tracker.state = 2;
            tracker.streak_occluded = 0;

            % build the dres structure
            dres_one = [];
            dres_one.fr = fr;
            dres_one.id = id;
            dres_one.x = dres_gt.x(1);
            dres_one.y = dres_gt.y(1);
            dres_one.w = dres_gt.w(1);
            dres_one.h = dres_gt.h(1);
            dres_one.r = 1;
            dres_one.state = 2;
            tracker.dres = dres_one;

        % tracking by association
        else

            % find the ground truths for association
            index_gt = find(dres_gt_all.fr == fr);
            dres = sub(dres_gt_all, index_gt);
            dres = MDP_crop_image_box(dres, dres_image.Igray{fr}, tracker);
            [dres, index_det, ctrack] = generate_association_index(tracker, fr, dres);

            % check if occluded or not
            index_gt = find(dres_gt.fr == fr);
            if isempty(index_gt) == 1 || dres_gt.covered(index_gt) ~= 0
                index_det = [];
            end

            if is_show
                figure(1);
                subplot(2, 3, 3);
                show_dres(fr, dres_image.I{fr}, 'Potential Associations', sub(dres, index_det));
                hold on;
                plot(ctrack(1), ctrack(2), 'ro', 'LineWidth', 2);
                hold off;
            end

            % compute features
            if isempty(index_det) == 0 && isempty(find(dres.id(index_det) == tracker.target_id, 1)) == 0
                % extract features with LK association
                dres_associate = sub(dres, index_det);
                features = MDP_feature_occluded(fr, dres_image, dres_associate, tracker);

                % compute labels
                m = size(features, 1);
                labels = -1 * ones(m, 1);
                ind = find(dres_associate.id == tracker.target_id);
                labels(ind) = 1;

                % update features
                tracker.f_occluded(end+1:end+m,:) = features;
                tracker.l_occluded(end+1:end+m) = labels;
                if is_text
                    fprintf('training examples in occluded state %d\n', size(tracker.f_occluded,1));
                end

                % update template
                dres_one = sub(dres_associate, ind);
                tracker = LK_associate(fr, dres_image, dres_one, tracker);

                tracker.prev_state = tracker.state;
                tracker.state = 2;
                % build the dres structure
                dres_one = [];
                dres_one.fr = fr;
                dres_one.id = tracker.target_id;
                dres_one.x = tracker.bb(1);
                dres_one.y = tracker.bb(2);
                dres_one.w = tracker.bb(3) - tracker.bb(1);
                dres_one.h = tracker.bb(4) - tracker.bb(2);
                dres_one.r = 1;
                dres_one.state = 2;

                if tracker.dres.fr(end) == fr
                    dres_tmp = tracker.dres;
                    index_tmp = 1:numel(dres_tmp.fr)-1;
                    tracker.dres = sub(dres_tmp, index_tmp);
                end
                tracker.dres = interpolate_dres(tracker.dres, dres_one);
                % update LK tracker
                tracker = LK_update(fr, tracker, dres_image.Igray{fr}, dres_associate, 1);

            else
                tracker.state = 3;
                dres_one = sub(tracker.dres, numel(tracker.dres.fr));
                dres_one.fr = fr;
                dres_one.id = tracker.target_id;
                dres_one.state = 3;

                if tracker.dres.fr(end) == fr
                    dres_tmp = tracker.dres;
                    index_tmp = 1:numel(dres_tmp.fr)-1;
                    tracker.dres = sub(dres_tmp, index_tmp);
                end
                tracker.dres = concatenate_dres(tracker.dres, dres_one);
            end

        end

        % show results
        if is_show
            figure(1);

            % show tracking results
            subplot(2, 3, 4);
            show_dres(fr, dres_image.I{fr}, 'Tracking', tracker.dres, 2);

            % show lost targets
            subplot(2, 3, 5);
            show_dres(fr, dres_image.I{fr}, 'Lost', tracker.dres, 3);

            subplot(2, 3, 6);
            show_templates(tracker, dres_image);

            fprintf('frame %d, state %d\n', fr, tracker.state);
            if is_pause
                pause();
            else
                pause(0.01);
            end

%             filename = sprintf('results/%s_%06d.png', seq_name, fr);
%             hgexport(h, filename, hgexport('factorystyle'), 'Format', 'png');
        end

        fr = fr + 1;
    end

end
tracker.w_occluded = svmtrain(tracker.l_occluded, tracker.f_occluded, '-c 1 -q -g 1 -b 1');
fprintf('Finish training %s\n', seq_name);

% save model
if is_save
    filename = sprintf('%s/%s_tracker_gt.mat', opt.results, seq_name);
    save(filename, 'tracker');
end

author: yuxng


text:
function varargout = spm_brainwarp(varargin)
% Part of old nonlinear spatial normalisation - a compiled routine
%_______________________________________________________________________
% [Alpha,Beta,Var] = spm_brainwarp(VG,VF,Affine,basX,basY,basZ,...
%                                   dbasX,dbasY,dbasZ,T,fwhm,VW,VW2)
% VG    - Template volume(s) (see spm_vol)
% VF    - Object volume
% Affine    - The affine transformation which maps between the object
%     and template.
% basX  - Basis vectors in X. # rows must eq. VG(1)
% basY  - Basis vectors in Y. # rows must eq. VG(2)
% basZ  - Basis vectors in Z. # rows must eq. VG(3)
% dbasX - Derivatives of basis vectors in X. # rows must eq. VG(1)
% dbasY - Derivatives of basis vectors in Y. # rows must eq. VG(2)
% dbasZ - Derivatives of basis vectors in Z. # rows must eq. VG(3)
% T - The current parameter estimates.
% fwhm  - The approximate smoothness of the images.
% VW    - an optional weighting volume for determining which voxels
%         should be weighted more heavily in the fitting process.
%         This volume should have the same dimensions and position
%         as the volumes in VG.
% VW2   - another optional weighting volume for determining which voxels
%         should be weighted more heavily in the fitting process.
%         This volume should have the same dimensions and position
%         as the volumes in VF.
% Without the weighting volumes, all voxels are assigned weights that
% are uniformly one.
%
% Alpha - A*A - where A is the design matrix
% Beta  - A*b - where f is the object image
% Var   - the approximate chi^2 (corrected for number of resels).
%_______________________________________________________________________
%
% The voxels of g1, g2.. are sampled according to the smoothness of the
% image (fwhm). The corresponding voxels of f are determined according
% to the current parameter estimates and the affine transform.  See
% "spm_write_sn.m" for more details about how this is done.
%
%
%-----------------------------------------------------------------------
%
% The design matrix A is generated internally from:
%
% diag(w)*[diag(df/dx)*B diag(df/dy)*B diag(df/dz)*B ...
%      diag(g1)*[1 x y z] ...
%      diag(g2)*[1 x y z] ...]
%
% where df/dx, df/dy & df/dz are column vectors containing the gradient
%   of image f with respect to displacements in x, y & z
%   (in the space of g).
%
%   B is generated from kron(basZ,kron(basY,BasX)). Each column of
%   B is a basis image.
%
%   g1, g2.. are template images.
%
%   x, y & z are simply the spatial coordinates of the voxels of f.
%
%   s1, s2.. are the current estimates for the required scaling
%   factors. These are derived from T(3*prod(VG(1:3))+1),
%   T(3*prod(VG(1:3))+2)...
%
%       w is an optional vector of weights, where w = 1/(1/w1 + 1/w2)
%       where w1 and w2 are derived from the optional weighting images.
%
% The vector b contains [diag(w)*(f - diag(g1)*s1 - diag(g1)*x*s2 - ...)].
%
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_brainwarp.m 4923 2012-09-13 15:01:01Z guillaume $

%-This is merely the help file for the compiled routine
error('spm_brainwarp.c not compiled - see Makefile')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% This is a testing environment for the files in the folder kernel/matrixfree
% 1. Based on data/contents, a list of required files is generated and it
%    is verified, that all files are present; additional files are listed.
% 2. All c-files are compiled.
% 3. All files are executed.
%==============================================================================

FAIRcheckFiles(mfilename);
testEnd;
%==============================================================================


text:
function [D,dir] = dimensions(points,varargin)

% Calculates the box-dimensions and dimension estimates of the point set
% "points". Returns also the corresponding direction vectors.

if nargin == 2
    P = varargin{1};
    points = P(points,:);
elseif nargin == 3
    P = varargin{1};
    Bal = varargin{2};
    I = vertcat(Bal{points});
    points = P(I,:);
end

if size(points,2) == 3
    X = cov(points);
    [U,S,~] = svd(X);

    dp1 = points*U(:,1);
    dp2 = points*U(:,2);
    dp3 = points*U(:,3);

    D = [max(dp1)-min(dp1) max(dp2)-min(dp2) max(dp3)-min(dp3) ...
        (S(1,1)-S(2,2))/S(1,1) (S(2,2)-S(3,3))/S(1,1) S(3,3)/S(1,1)];

    dir = [U(:,1)' U(:,2)' U(:,3)'];
else
    X = cov(points);
    [U,S,~] = svd(X);

    dp1 = points*U(:,1);
    dp2 = points*U(:,2);

    D = [max(dp1)-min(dp1) max(dp2)-min(dp2) ...
        (S(1,1)-S(2,2))/S(1,1) S(2,2)/S(1,1)];

    dir = [U(:,1)' U(:,2)'];
end

author: InverseTampere

repo: TreeQSM


text:
N1=2^16; % we are working with 2 seconds frames.

N2=2^12;
db_path='envsounds/solosDbdb/o2_scat/';
T=2^12; %128 ms for fs= 32 KHz
order=2;
Q1=1;
Q2=1;

%set the filters options
filt1_opt.filter_type = {'gabor_1d','morlet_1d'};
filt1_opt.Q = [Q1 Q2];
filt1_opt.J = T_to_J(T,filt1_opt);

% set the scattering options
sc1_opt.M=2;
sc1_opt.use_abs=1; %used in wavelet_1d, preparing for renormalization with conv(abs(x),phi)

%Create the renormalization operator
Wop1=normalize_wavelet_factory_1d(N1,filt1_opt,sc1_opt);

scatt_fun1 = @(x)(scat(x,Wop1));
feature_fun1={@(x)(format_scat(aggregate_scat(scatt_fun1(x),T)))};
[~,meta]=feature_fun1{1}(rand(N1,1));

src=instsnds_src_multi_obj('envsounds/sDbCoVl',N1);

prep_opt.parallel=1;
prep_opt.average=1;
prep_opt.sig_normalize=1;

db1=prepare_inst_db(src,feature_fun1,prep_opt);
db1.features=single(db1.features);

sigmas=threshold_parameters_factory(db1,meta,'median');
save([db_path 'sigmas_sDbCoVl_T12_Q11.mat'],'sigmas');

%%%sigmas=load('sigmas_sDbCoVl_T12_Q11.mat'); sigmas=sigmas.sigmas;
clear db1;

thresholds=[0 1e-5 1e-4 1e-3 1e-2 1e-1 1 2 3];

lgth=length(thresholds);

all_recog_rates={};
all_recog=zeros(size(thresholds));

for l=1:lgth

    sc1_opt.threshold=thresholds(l);
    sc1_opt.sigmas=sigmas;
    sc1_

    %Prepare the normalizing standardizing and thresholding operator
    Wop2=normalize_standardize_threshold_wavelet_factory_1d(N1,filt1_opt,sc1_opt);

    scatt_fun2=@(X)(format_scat(scat(X,Wop2)));

    feature_fun={@(X)(scatt_fun2(X))};

    [~ ,meta]=feature_fun{1}(rand(N1,1));

    db=prepare_inst_db(src,feature_fun,prep_opt);
    db.features=single(db.features);

    %%Remove silent frames;
%     if ~isfield(db.src,'rejects')
%         db.src=filter_src(db.src);
%     end

    db=poor_whiten2(db);

    db=svm_calc_kernel(db,'gaussian','square',1:size(db.features,2));

    partitions = load('sDbCoVl_fparts.mat');

    files_train=partitions.files_train;
    files_test=partitions.files_test;

    [train_set ,test_set]= create_filt_feats_parts(db.src, files_train,files_test,[],[]);

    optt.kernel_type = 'gaussian';

    optt.C=2^16;
    optt.gamma=2^(-8);

    % The data is unbalanced, use weights
    optt.w=1;

    labels={};
    feature_lbls={};
    votes={};

    C=optt.C;
    gamma=optt.gamma;

    optt1 = optt;

    model = svm_train(db,train_set,optt1,db_weights);
    [labels votes feature_lbls] = svm_test(db,model,test_set);

    [recog recog_rates]=classif_recog(labels,test_set,db.src)
    err = 1-recog;

    all_recog_rates{l}=recog_rates;
    all_recog(l)=recog;
end

save('all_recog_rates_sDbCoVl.mat','all_recog_rates','all_recog','thresholds');
sendmail('michelkapoko@gmail.com','take this best cf','check',{'all_recog_rates_sDbCoVl.mat'});

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
%--- help for mdd.global_options ---
%
%mdd.global_options is a function.
%    opts = global_options(opts0)
%

author: jmaih


text:
%% load an RGBD sequence
data = loadSUN3Dv2('2014-05-01_19-32-43_260595134347');

%% demo to show how to use a frame
frameID = 11;

%% load the frame data
image = readImage(data,frameID);
depth = readDepth(data,frameID); % use this to read the denoised depth map that is closet to the image

IR = readIR(data,frameID);

%% get 3D point cloud from the depth map
XYZcamera(:,:,1)=data.camera.D.X .* depth;
XYZcamera(:,:,2)=data.camera.D.Y .* depth;
XYZcamera(:,:,3)=depth .* (~isnan(data.camera.D.X) & ~isnan(data.camera.D.Y));
XYZcamera(:,:,4)=depth>0 & ~isnan(data.camera.D.X) & ~isnan(data.camera.D.Y);
valid = logical(XYZcamera(:,:,4));
valid = valid(:)';
XYZ = reshape(XYZcamera,[],4)';
XYZ = XYZ(1:3,valid);

points2ply('demo_all_points.ply', XYZ);

%% get the color for each 3D point by tranforming and projecting each 3D point
XYZrgb = transformPointCloud(XYZ,data.camera.D2RGB.Rt);
xyRGB = 1+project_points2(XYZrgb,data.camera.RGB.fc,data.camera.RGB.cc,data.camera.RGB.kc,data.camera.RGB.alpha_c);
uvRGB = round(xyRGB);
valid = 1<= uvRGB(1,:) & uvRGB(1,:) <= data.camera.RGB.width & 1<= uvRGB(2,:) & uvRGB(2,:) <= data.camera.RGB.height;
XYZ = XYZ(1:3,valid);
ind = sub2ind([data.camera.RGB.height data.camera.RGB.width],uvRGB(2,valid),uvRGB(1,valid));
RGB =[image(ind); image(ind+data.camera.RGB.height*data.camera.RGB.width); image(ind+data.camera.RGB.height*data.camera.RGB.width*2)];
points2ply('demo_points_with_color.ply', XYZ, RGB);

%% project the 3D points to the undistort color image
xyzRGBnormal = data.camera.RGB.K2render * XYZrgb;
xyzRGBnormal(1,:) = 1+ xyzRGBnormal(1,:) ./ xyzRGBnormal(3,:);
xyzRGBnormal(2,:) = 1+ xyzRGBnormal(2,:) ./ xyzRGBnormal(3,:);

%% undistort the RGB image into an pin hole camera with data.camera.RGB.K2render as intrinsics
undistortImage = undistort_image(im2double(image), data.camera.RGB.KK, data.camera.RGB.kc, data.camera.RGB.width, data.camera.RGB.height, data.camera.RGB.K2render);

%% undistort the depth map, transform the depth, and project it using OpenGL to the undistorted image with data.camera.RGB.K2render as intrincis
% manually find one
%[~,imageDepth] = WarpDepthMatlab(XYZcamera,data.camera.RGB.K2render, data.camera.D2RGB.Rt, data.camera.RGB.width, data.camera.RGB.height);
% or better to use this one
imageDepth = depth4RGB(data, frameID); % use this to read a depth map that is a combination of the closest two depth map
% another useful function to do the same thing is
% frames = getRGBDframe(sequenceName,frameID);

%% undistort the IR image. Typicaly, you don't need this
undistortIR = undistort_image(im2double(IR), data.camera.D.KK, data.camera.D.kc, data.camera.D.width, data.camera.D.height);

%% undistore the depth map to a pixel camera, and do NOT align it with the RGB image
[~,undistortDepth] = WarpDepthMatlab(XYZcamera,data.camera.D.KK, [eye(3) zeros(3,1)], data.camera.D.width, data.camera.D.height);

%% visualization
figure
subplot(2,3,1);
imagesc(depth);
axis equal; axis tight; axis off;
title('raw depth');

subplot(2,3,2);
imshow(image);
hold on;
plot(xyRGB(1,:),xyRGB(2,:),'.')
axis equal;
axis tight
title('raw image');

subplot(2,3,3);
imshow(undistortImage);
hold on;
plot(xyzRGBnormal(1,:),xyzRGBnormal(2,:),'.')
axis equal;
axis tight
title('undistorted image');

subplot(2,3,4);
imagesc(undistortDepth);
axis equal; axis tight; axis off;
title('undistorted depth');

subplot(2,3,5);
imshow(undistortIR);
title('undistorted IR');

subplot(2,3,6);
imagesc(imageDepth);
axis equal; axis tight; axis off;
title('depth on undistorted image');

author: jianxiongxiao

repo: ProfXkit


text:
function pass = test_empty( pref )
% For empty chebfun2v objects, does each command deal with them
% appropriately?
% Alex Townsend, March 2013.

%%
% Check that the chebfun2v commands work for empty Chebfun2v objects.

pass = 1;
F = chebfun2v;
try
    conj(F);
    cross(F,F);
    F.';
    curl(F);
    divergence(F);
    dot(F,F);
    F(1,1);
    imag(F);
    isempty(F);
    laplacian(F);
    F-F;
    F.*F;
    norm(F);
    F+F;
    F.^2;
    real(F);
    roots(F);
    F*F;
    F';
catch
   pass = 0;
end

end

author: chebfun

repo: chebfun


text:
classdef KnEA < ALGORITHM
% <many> <real/integer/label/binary/permutation> <constrained/none>
% Knee point driven evolutionary algorithm
% rate --- 0.5 --- Rate of knee points in the population

%------------------------------- Reference --------------------------------
% X. Zhang, Y. Tian, and Y. Jin, A knee point-driven evolutionary algorithm
% for many-objective optimization, IEEE Transactions on Evolutionary
% Computation, 2015, 19(6): 761-776.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            rate = Algorithm.ParameterSet(0.5);

            %% Generate random population
            Population = Problem.Initialization();
            FrontNo    = NDSort(Population.objs,Population.cons,inf);
            KneePoints = zeros(1,Problem.N);     % Set of knee points
            r          = -ones(1,2*Problem.N);	% Ratio of size of neighorhood
            t          = -ones(1,2*Problem.N);	% Ratio of knee points

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = MatingSelection(Population.objs,FrontNo,KneePoints);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                Population = [Population,Offspring];
                [FrontNo,MaxFNo]                = NDSort(Population.objs,Population.cons,Problem.N);
                [KneePoints,Distance,r,t]       = FindKneePoints(Population.objs,FrontNo,MaxFNo,r,t,rate);
                [Population,FrontNo,KneePoints] = EnvironmentalSelection(Population,FrontNo,MaxFNo,KneePoints,Distance,Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%  Figure 10.08      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% fig10_08.m is a script to generate Fig. 10.8  the
% frequency response of the satellite with low-gain PD compensation

np =[0.0360    0.9100];

dp =[1.0000    0.0396    1.0010    0.0000    0.0000];

nc2=0.001*[30 1];

nol2=conv(nc2,np);
dol2=dp;

hold off ; clf
w=logspace(-2,.2);
w(46)=1;
[magol2, phol2]= bode(nol2,dol2,w);
subplot(211) ; loglog(w,magol2); grid; hold on;
xlabel('\omega (rad/sec)');
ylabel('Magnitude |D_2(s)G(s)|');
loglog(w,ones(size(magol2)),'g');
title('Fig. 10.8 Frequency response of low-gain satellite PD design')
phol2a=[phol2, -180*ones(size(phol2))];
subplot(212);  semilogx(w, phol2a); grid; hold on;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');

author: Sable