text:
function x = line_ccvt_lloyd ( n, a, b, it_num, header, x )

%*****************************************************************************80
%
%% LINE_CCVT_LLOYD carries out the constrained Lloyd algorithm.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of generators.
%
%    Input, real A, B, the left and right endpoints.
%
%    Input, integer IT_NUM, the number of iterations to take.
%
%    Input, string HEADER, an identifying string.
%
%    Input, real X(N), the initial point locations.
%
%    Output, real X(N), the final point locations.
%
  x = x(:);
%
%  Print the initial generators.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial generators:\n' );
  fprintf ( 1, '\n' );
  for k = 1 : n
    fprintf ( 1, '  %4d  %f\n', k, x(k,1) );
  end
%
%  Initialize the plotting arrays.
%
  step = 1 : it_num;
  e = nan ( it_num, 1 );
  xm = nan ( it_num, 1 );

  for it = 1 : it_num

    x_plot ( 1:n, it ) = x(1:n,1);

    x_new = line_ccvt_lloyd_step ( n, a, b, x );

    e(it) = line_cvt_energy ( n, a, b, x );
    e(it) = max ( e(it), eps );
%
%  Display the energy.
%
    figure ( 1 )
    plot ( step, log ( e ), 'm-*' )
    title ( 'Log (Energy)' )
    xlabel ( 'Step' )
    ylabel ( 'Energy' )
    grid
%
%  Compute the generator motion.
%
    xm(it,1) = sum ( ( x_new(:) - x(:) ).^2 ) / n;
%
%  Display the generator motion.
%
    figure ( 2 )
    plot ( step, log ( xm ), 'm-*' )
    title ( 'Log (Average generator motion)' )
    xlabel ( 'Step' )
    ylabel ( 'Motion' )
    grid
%
%  Update the generators.
%
    x(1:n,1) = x_new(1:n,1);

  end

  x_plot(1:n,it_num+1) = x(1:n,1);
%
%  Print the current generators.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Current generators:\n' );
  fprintf ( 1, '\n' );
  for k = 1 : n
    fprintf ( 1, '  %4d  %f\n', k, x(k,1) );
  end
%
%  Plot the evolution of the locations of the generators.
%
  figure ( 3 )

  y = ( 0 : it_num );
  for k = 1 : n
    plot ( x_plot(k,1:it_num+1), y, 'LineWidth', 1 )
    hold on;
  end
  grid on
  hold off;

  title ( 'Generator evolution.' );
  xlabel ( 'Generator positions' );
  ylabel ( 'Iterations' );
%
%  Save the plots.
%
  figure ( 1 )
  filename = strcat ( header, '_energy.png' );
  print ( '-dpng', filename );
  figure ( 2 )
  filename = strcat ( header, '_motion.png' );
  print ( '-dpng', filename );
  figure ( 3 )
  filename = strcat ( header, '_evolution.png' );
  print ( '-dpng', filename );

  return
end


function Spath = Simulate_General_paths_func( N_sim, M, mult, T, S_0, r, q, model, modelParams, jumpModel, jumpParams )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
M_mult = M*mult;  %Total number of Sampled Stock prices per path
modelType = modelParams.modelType;

if modelType == 1 % Stochastic Volatility (with Jumps)
    Spath = Simulate_StochVol_Jumps_func( N_sim, M_mult+1, T, S_0, r, q, model, modelParams, jumpModel, jumpParams);

elseif modelType == 2 % Jump Diffusion (and regular diffusion, e.g. Black Scholes)
    sigma = modelParams.sigma;
    Spath = Simulate_Jump_Diffusion_func( N_sim, M_mult+1, T, S_0, r, q, sigma, jumpModel, jumpParams);

elseif modelType == 3  % Stochastic Local Volatility (SLV)
    Spath = Simulate_SLV_func( N_sim, M_mult+1, T, S_0, r, q, model, modelParams);
end

end


text:
function value = stirling2_value ( n, m )

%*****************************************************************************80
%
%% STIRLING2_VALUE computes a Stirling number of the second kind.
%
%  Discussion:
%
%    S2(N,M) represents the number of distinct partitions of N elements
%    into M nonempty sets.  For a fixed N, the sum of the Stirling
%    numbers S2(N,M) is represented by B(N), called "Bell's number",
%    and represents the number of distinct partitions of N elements.
%
%    For example, with 4 objects, there are:
%
%    1 partition into 1 set:
%
%      (A,B,C,D)
%
%    7 partitions into 2 sets:
%
%      (A,B,C) (D)
%      (A,B,D) (C)
%      (A,C,D) (B)
%      (A) (B,C,D)
%      (A,B) (C,D)
%      (A,C) (B,D)
%      (A,D) (B,C)
%
%    6 partitions into 3 sets:
%
%      (A,B) (C) (D)
%      (A) (B,C) (D)
%      (A) (B) (C,D)
%      (A,C) (B) (D)
%      (A,D) (B) (C)
%      (A) (B,D) (C)
%
%    1 partition into 4 sets:
%
%      (A) (B) (C) (D)
%
%    So S2(4,1) = 1, S2(4,2) = 7, S2(4,3) = 6, S2(4,4) = 1, and B(4) = 15.
%
%
%  First terms:
%
%    N/M: 1    2    3    4    5    6    7    8
%
%    1    1    0    0    0    0    0    0    0
%    2    1    1    0    0    0    0    0    0
%    3    1    3    1    0    0    0    0    0
%    4    1    7    6    1    0    0    0    0
%    5    1   15   25   10    1    0    0    0
%    6    1   31   90   65   15    1    0    0
%    7    1   63  301  350  140   21    1    0
%    8    1  127  966 1701 1050  266   28    1
%
%  Recursion:
%
%    S2(N,1) = 1 for all N.
%    S2(I,I) = 1 for all I.
%    S2(I,J) = 0 if I < J.
%
%    S2(N,M) = M * S2(N-1,M) + S2(N-1,M-1)
%
%  Properties:
%
%    sum ( 1 <= K <= M ) S2(I,K) * S1(K,J) = Delta(I,J)
%
%    X**N = sum ( 0 <= K <= N ) S2(N,K) X_K
%    where X_K is the falling factorial function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of rows of the table.
%
%    Input, integer M, the number of columns of the table.
%
%    Output, integer VALUE, the value of S2(N,M).
%
  if ( n <= 0 )
    value = 0;
    return
  end

  if ( m <= 0 )
    value = 0;
    return
  end

  s2(1,1) = 1;
  s2(1,2:m) = 0;

  for i = 2 : n

    s2(i,1) = 1;

    for j = 2 : m
      s2(i,j) = j * s2(i-1,j) + s2(i-1,j-1);
    end

  end

  value = s2(n,m);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef prtPreProcLabelVoting < prtPreProc
    % Switches labels of the training data to be the most popular label amoungst k neihbors
    % This isn't necessarily a great idea but it has been proposed by several people.
    %
    % Example:
    %
    % ds = prtDataGenXor;
    %
    % eta = 0.2;
    % flipLabel = rand(ds.nObservations,1) < eta;
    %
    % newY = ds.Y;
    % newY((ds.Y == 0) & flipLabel) = 1;
    % newY((ds.Y == 1) & flipLabel) = 0;
    %
    % dsFlippedY = ds;
    % dsFlippedY.Y = newY;
    %
    % class = train(prtClassMap('rvs',prtRvGmm('nComponents',2)), dsFlippedY);
    %
    % subplot(1,2,1)
    % plot(class)
    %
    % class = train(prtPreProcLabelVoting('k',10) + prtClassMap('rvs',prtRvGmm('nComponents',2)), dsFlippedY);
    %
    % subplot(1,2,2)
    % plot(class.actionCell{2})

    properties (SetAccess=private)
        name = 'Label Voting'
        nameAbbreviation = 'LVNN'
    end
    properties
        k = 3;
    end

    methods
        function self = prtPreProcLabelVoting(varargin)
            self = prtUtilAssignStringValuePairs(self, varargin{:});
        end
    end
    methods (Access=protected,Hidden=true)
        function self = preTrainProcessing(self,DataSet)
            if ~self.verboseStorage
                warning('prtPreProcLabelVoting:verboseStorage:false','prtPreProcLabelVoting requires verboseStorage to be true; overriding manual settings');
            end
            self.verboseStorage = true;
            self = preTrainProcessing@prtPreProc(self,DataSet);
        end
        function self = trainAction(self,~)
            %Do nothing; we've already specified "verboseStorage = true",
            %so the ".dataSet" field will be set when it comes time to test
        end

        function ds = runActionOnTrainingData(self, ds)

            distanceMat = prtDistanceEuclidean(ds.X,ds.X);
            [~,I] = sort(distanceMat,1,'ascend');
            I = I(1:self.k+1,:);
            L = ds.Y(I);

            if self.k ~= 1
                L = L';
            end
            L = L(:,2:end); % First one is yourself...

            uClasses = ds.uniqueClasses;
            confOut = zeros(size(L,1), ds.nClasses);
            for iY = 1:ds.nClasses
                confOut(:,iY) = sum(L == uClasses(iY),2);
            end
            [~, classInd] = max(confOut,[],2);

            ds.Y = uClasses(classInd);

        end
        function ds = runAction(~, ds)
            % Nadda
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function	[Wout,Weff] = make_test_weight(parm)
% Make output weigtht for test
% -Input
% parm.Ydim   =  Output dim
% parm.Xdim   =  Input dim
% parm.Meff   =  Effective Input dim
% parm.Tau    =  Lag time steps
% parm.Dtau   =  Number of embedding dimension
% parm.WXmax = maximum weight of X-dim
% parm.WYmax = bias weight for Y-dim
% parm.Wtau  = effctive time delay length of weight
% -Output
% Wout  : Weight for time embedded input [Ydim, Xdim*Dtau]
% Weff  : nonzero weight for [Ydim, Meff*Dtau]
%         effective input dimension : [Meff x Wtau]

% Input dim
if isfield(parm,'Xdim')
	Xdim = parm.Xdim ;
else
	Xdim = 100;
end
% Output dim
if isfield(parm,'Ydim')
	Ydim = parm.Ydim ;
else
	Ydim = 1;
end
% Effective Input dim
if isfield(parm,'Meff')
	Meff = parm.Meff ;
else
	Meff = 10; ;
end
% Number of embedding dimension
if isfield(parm,'Dtau')
	Dtau  = parm.Dtau  ;
else
	Dtau  = 1;
end

if Meff > Xdim, Meff = Xdim; end;
if Dtau==0, Dtau=1; end;

% effctive time delay length of weight
if isfield(parm,'Wtau')
	tau  = parm.Wtau  ;
else
	tau  = Dtau;
end
% maximum weight of X-dim
if isfield(parm,'WXmax')
	WXmax = parm.WXmax;
else
	WXmax = 100;
end
% maximum bias weight for Y-dim
if isfield(parm,'WYmax')
	WYmax = parm.WYmax;
else
	WYmax = 10;
end

%
% --- Make Output Weight matrix
%
% Output weight

% Spatial weight
Wbase = rand(Ydim,Meff)*WXmax + WYmax;

% Temporal filter
t = (0:tau-1)/max((tau-1)/2,1);
Wtau  = exp( - 2* t.^2);

ix_eff = 1:Meff;
Weff   = zeros(Ydim,Meff*Dtau);
Wout   = zeros(Ydim,Xdim*Dtau);

for n = 1:tau
	Weff(:,ix_eff + Meff*(n-1)) = Wtau(n) * Wbase;
	Wout(:,ix_eff + Xdim*(n-1)) = Wtau(n) * Wbase;
end

author: KamitaniLab

repo: GenericObjectDecoding


text:
function c = eval_get_color_new(cidxs)

% http://colorbrewer2.org/

color = cell(8,4);

% qualitative
color{1}{1} = [55,126,184];
color{2}{1} = [255,127,0];
color{3}{1} = [255,255,51];%[178,223,138];
color{4}{1} = [0,0,0];
color{5}{1} = [77,175,74];
color{6}{1} = [228,26,28];
% color{6}{1} = [255,255,51];
color{7}{1} = [152,78,163];
color{8}{1} = [247,129,191];

% sequential
color{1}{2} = [253,141,60];
color{1}{3} = [254,204,92];
color{1}{4} = [255,255,178];
color{1}{5} = [254,240,217];

color{2}{2} = [107,174,214];
color{2}{3} = [189,215,231];
color{2}{4} = [117,107,177];

color{3}{2} = [194,230,153];

color{4}{2} = [140,150,198];
color{4}{3} = [179,205,227];
color{4}{4} = [237,248,251];

c = color{cidxs(1)}{cidxs(2)};

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function distance = Distance(positions,reference,varargin)

%Distance - Compute instantaneous distance to a reference point.
%
%  USAGE
%
%    distance = Distance(positions,reference,<options>)
%
%    positions      a list of position samples
%    reference      reference point (same coordinates as positions)
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'type'        'l' if X is linear (default), 'c' if X is circular (for
%                    1D positions only, which should be in [0..1])
%    =========================================================================
%
%  SEE
%
%    See also Diff.

% Copyright (C) 2004-2012 by Michaël Zugaro
%
% 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.

% Default values
type = 'linear';

if nargin < 2 | mod(length(varargin),2) ~= 0,
	error('Incorrect number of parameters (type ''help Distance'' for details).');
end
if ~issamples(positions,'#1') && ~issamples(positions,'#2'),
	error('Incorrect positions - should be a vector (type ''help Distance'' for details).');
end
if ~isdvector(reference,'#1') && ~isdvector(reference,'#2'),
	error('Incorrect reference - should be a vector (type ''help Distance'' for details).');
end
if length(reference) ~= size(positions,2)-1,
	error('Positions and reference have incompatible sizes (type ''help Distance'' for details).');
end

% Parse parameter list
for j = 1:2:length(varargin),
	if ~ischar(varargin{j}),
		error(['Parameter ' num2str(j+2) ' is not a property (type ''help Distance'' for details).']);
	end
	switch(lower(varargin{j})),
		case 'type',
			type = varargin{j+1};
			if ~isstring_FMAT(type,'l','c'),
				error('Incorrect value for property ''type'' (type ''help Distance'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{j}) ''' (type ''help Distance'' for details).']);
	end
end

distance = [];
if isempty(positions), return; end

distance = positions(:,1);

if issamples(positions,'#1'),
	if type(1) == 'l',
		distance(:,2) = abs(positions(:,2)-reference(1));
	else
		% Make sure X is normalized
		if max(positions(:,2)) > 1 || min(positions(:,2)) < 0,
			positions(:,2) = ZeroToOne(positions(:,2));
			warning('Positions should contain values in [0 1]. The data will now be transformed accordingly.');
		end
		distance1 = abs(positions(:,2)-reference(1));
		distance2 = 1-abs(positions(:,2)-reference(1));
		distance(:,2) = min([distance1 distance2],[],2);
	end
else
	distance(:,2) = sqrt((positions(:,2)-reference(1)).^2+(positions(:,3)-reference(2)).^2);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function scrollsubplot(dx,x,varargin);

% Cette fonction permet l'ajout d'un curseur dynamique sur les figures
% graphique tout en permettant l'utilisation de sous-graphiques subplot.
% Cette fonction est une variation de la fonction scrollplot de Wolfgang
% Stiegmaier et de scrollplotdemo de Steven Lord, slord@mathworks.com.
%
% Syntaxe
%   scrollsubplot(dx,x,h1,h2,...hn);
%
% Description
% dx: détermine la largeur de la fenetre
% x: vecteur temps
% h1...hn: les références aux sous-graphiques
%
% exemple
%
% temps = linspace(0,10*pi,100001);
% Y1 = sin(temps);
% Y2 = cos(temps);
% Y3 = sin(temps).^2;
%
% h1=subplot(3,1,1);
% plot(temps,Y1),grid;
% h2=subplot(3,1,2);
% plot(temps,Y2),grid;
% h3=subplot(3,1,3);
% plot(temps,Y3),grid;
% scrollsubplot(pi,temps,h1,h2,h3);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Benoit Cantin, 7 mars 2005. %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Original message
% :Created by Steven Lord, slord@mathworks.com
% :Uploaded to MATLAB Central
% :http://www.mathworks.com/matlabcentral
% :7 May 2002
% :Permission is granted to adapt this code for your own use.
% :However, if it is reposted this message must be intact.

% get current axes
a=gca;
% This avoids flickering when updating the axis
set(gcf,'doublebuffer','on');
% Set appropriate axis limits and settings
set(a,'xlim',[0 dx]);
%set(a,'ylim',[min(y) max(y)]);

% Generate constants for use in uicontrol initialization
pos=get(a,'position');
xmax=max(x);
xmin=min(x);

% This will create a slider which is just underneath the axis
% but still leaves room for the axis labels above the slider
Newpos=[pos(1) pos(2)-0.1 pos(3) 0.05];

for i = 1:length(varargin);
    %initialize postion of plot
    set(varargin{i},'xlim',[xmin xmin+dx]);

    % Création des variables h1...hn
    eval(['h',num2str(i), ' = varargin{i};']);

    % Setting up callback string to modify XLim of axis (gca)
    % based on the position of the slider (gcbo)
    if i == 1;
        S=['set(h', num2str(i), ' ,''xlim'',get(gcbo,''value'')+[0 ' num2str(dx) '])'];
    else
        S=[S ', set(h', num2str(i), ' ,''xlim'',get(gcbo,''value'')+[0 ' num2str(dx) '])'];
    end
end

% Creating Uicontrol with initial value of the minimum of x
h=uicontrol('style','slider',...
    'units','normalized','position',Newpos,...
    'callback',S,'min',xmin,'max',xmax-dx,'value',xmin);

author: Sable

repo: mcbench-benchmarks


text:
classdef SuperEllipseRhoBoundsComputer < handle

    properties (Access = private)
       q
       txi
       mxMax
       myMax
       mxMin
       myMin
       superEllipse
    end

    methods (Access = public)

        function obj = SuperEllipseRhoBoundsComputer(cParams)
            obj.init(cParams);
        end

        function [rhoMin,rhoMax] = compute(obj)
            rhoMin = obj.computeRhoMin();
            rhoMax = obj.computeRhoMax();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.q = cParams.q;
            obj.txi = cParams.txi;
            obj.mxMin = cParams.mxMin;
            obj.myMin = cParams.myMin;
            obj.mxMax = cParams.mxMax;
            obj.myMax = cParams.myMax;
            obj.superEllipse = SuperEllipseParamsRelator;
        end

        function rhoMin = computeRhoMin(obj)
            rhoMinMx = obj.rhoFromMx(obj.mxMax);
            rhoMinMy = obj.rhoFromMy(obj.myMax);
            rhoMin = max(rhoMinMx,rhoMinMy);
        end

        function rhoMax = computeRhoMax(obj)
            rhoMaxMx = obj.rhoFromMx(obj.mxMin);
            rhoMaxMy = obj.rhoFromMy(obj.myMin);
            rhoMax = min(rhoMaxMx,rhoMaxMy);
        end

        function rho = rhoFromMx(obj,mx)
            rho = obj.superEllipse.rhoFromMxAndTxi(mx,obj.txi,obj.q);
        end

        function rho = rhoFromMy(obj,my)
            rho = obj.superEllipse.rhoFromMyAndTxi(my,obj.txi,obj.q);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%% Test Problems for LEVMAR Interface
% J.Currie 2012
clc
clear all

%% NLS2 [0.3578]
clc
%Function
i = (1:40)';
fun = @(x) x(1)*exp(-x(2)*i) + x(3);
%Fitting Data
ydata=[5.8728, 5.4948, 5.0081, 4.5929, 4.3574, 4.1198, 3.6843, 3.3642, 2.9742, 3.0237, 2.7002, 2.8781,...
       2.5144, 2.4432, 2.2894, 2.0938, 1.9265, 2.1271, 1.8387, 1.7791, 1.6686, 1.6232, 1.571, 1.6057,...
       1.3825, 1.5087, 1.3624, 1.4206, 1.2097, 1.3129, 1.131, 1.306, 1.2008, 1.3469, 1.1837, 1.2102,...
       0.96518, 1.2129, 1.2003, 1.0743];
% %Setup Options
opts = optiset('solver','levmar');
% %Build & Solve
x0=[1.0; 1.0; 1.0];
Opt = opti('fun',fun,'ydata',ydata,'ndec',3,'options',opts)
[x,fval,exitflag,info] = solve(Opt,x0)

%%
clc
grad = @(x) mklJac(fun,x);

opts = [];
opts.display = 2;
opts.maxiter = 50;
lb = [];[0;0.5;2];
ub = [];[10;10;10];

[x,f,e,i,ii] = levmar(fun,[],x0,ydata,lb,ub,[],[],[],[],opts)

%% NLS2 [0.3578]
clc
%Function
i = (1:40)';
fun = @(x) exp(-x(1)*i);
%Fitting Data
ydata=[5.8728, 5.4948, 5.0081, 4.5929, 4.3574, 4.1198, 3.6843, 3.3642, 2.9742, 3.0237, 2.7002, 2.8781,...
       2.5144, 2.4432, 2.2894, 2.0938, 1.9265, 2.1271, 1.8387, 1.7791, 1.6686, 1.6232, 1.571, 1.6057,...
       1.3825, 1.5087, 1.3624, 1.4206, 1.2097, 1.3129, 1.131, 1.306, 1.2008, 1.3469, 1.1837, 1.2102,...
       0.96518, 1.2129, 1.2003, 1.0743];
% %Setup Options
opts = optiset('solver','levmar');
% %Build & Solve
x0=[1.0];
Opt = opti('fun',fun,'ydata',ydata,'ndec',1,'options',opts)
[x,fval,exitflag,info] = solve(Opt,x0)

%%
grad = @(x) mklJac(fun,x);

opts = [];
opts.display = 2;
opts.maxiter = 50;
lb = [0];
ub = [10];

[x,f,e,i,ii] = nl2sol(fun,grad,x0,ydata,[],[],opts)

%%
clc
prob = nls_prob(20);
opts = optiset('solver','nl2sol','display','iter');
Opt = opti(prob,opts)

[x,f,e,i] = solve(Opt)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function slshow()
    % slshow slide show of images in a directory
    % now essentially self-contained.
    % - CGR

    report_this_filefun();
    setup_slideshow();
end

function setup_slideshow()
    hodi = ZmapGlobal.Data.hodi;
    file_extensions = {'jpg','png','bmp','gif'};
    slides = dir('this will be an empty anwser');
    for n=1:numel(file_extensions)
        slides = [slides; dir(fullfile(hodi, 'resrc','slides', ['*.', file_extensions{n}]))];  %#ok<AGROW>
    end

    fn = fullfile(slides(1).folder,{slides.name});
    slides = sort(fn);

    % set up slide viewer

    fsl = figure_w_normalized_uicontrolunits('Menubar','none','NumberTitle','off');
    viewport = axes('pos',[0 0 1 1]);axis off; axis ij

    userdata.slides = slides;
    userdata.idx = 0;
    userdata.ga = viewport;

    uicontrol(fsl,...
        'Units','normalized',...
        'Callback',@(s,e)close(fsl),...
        'String','Close',...
        'BackgroundColor' ,[0.8 0.8 0.8],...
        'Position',[0.90 0.01 0.10 0.06]);

    next_button = uicontrol(fsl,...
        'Units','normalized',...
        'Callback',@next_slide,...
        'String','Next ',...
        'BackgroundColor' ,[0.8 0.8 0.8],...
        'Position',[0.90 0.12 0.10 0.06],...
        'UserData', userdata);
    next_slide(next_button,[]);
end

function next_slide(s, ~)
    s.UserData.idx = s.UserData.idx + 1;
    if s.UserData.idx > numel(s.UserData.slides)
        s.UserData.idx = 1;
    end
    this_slide_info = s.UserData.slides{s.UserData.idx};
    try
        [x,~] = imread(this_slide_info);
    catch ME
        error_handler(ME, @do_nothing)
    end

    image(s.UserData.ga, x);
    drawnow
end

author: CelsoReyes


text:
function tno_sparse_grid_write ( header, l, m, n, x, w )

%*****************************************************************************80
%
%% TNO_SPARSE_GRID_WRITE writes a Truncated Normal Odd Sparse Grid to X and W files.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string HEADER, a header for the file names.
%
%    Input, integer L, the level of the rule.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real X(M,N), the abscissas.
%
%    Input, real W(N), the weights.
%
  w_filename = sprintf ( '%s_d%d_level%d_w.txt', header, m, l );
  x_filename = sprintf ( '%s_d%d_level%d_x.txt', header, m, l );

  r8vec_write ( w_filename, n, w );
  r8mat_write ( x_filename, m, n, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  W data written to "%s".\n', w_filename );
  fprintf ( 1, '  X data written to "%s",\n', x_filename );

  return
end


text:
function test_suite = test_mergeBoxes
%testMergeBoxes  One-line description here, please.
%
%   output = testMergeBoxes(input)
%
%   Example
%   testMergeBoxes
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-08-06,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testBasic(testCase) %#ok<*DEFNU>
box1 = [5 20 10 25];
box2 = [0 15 15 20];
exp  = [0 20 10 25];
box = mergeBoxes(box1, box2);
testCase.assertEqual(exp, box, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% unique_bc - unique backward compatible with Matlab versions prior to 2013a

% Copyright (C) 2013 Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [C,IA,IB] = unique_bc(A,varargin)

errorFlag = error_bc;

v = version;
indp = find(v == '.');
v = str2num(v(1:indp(2)-1));
if v > 7.19, v = floor(v) + rem(v,1)/10; end

if nargin > 2
    ind = strmatch('legacy', varargin);
    if ~isempty(ind)
        varargin(ind) = [];
    end
end

if v >= 7.14
    [C,IA,IB] = unique(A,varargin{:},'legacy');
    if errorFlag
        [C2,IA2] = unique(A,varargin{:});
        if ~isequal(C, C2) || ~isequal(IA, IA2) || ~isequal(IB, IB2)
            warning('backward compatibility issue with call to unique function');
        end
    end
else
    [C,IA,IB] = unique(A,varargin{:});
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function triangle_symq_rule_test05 ( degree, numnodes, vert1, vert2, vert3 )

%*****************************************************************************80
%
%% TRIANGLE_SYMQ_RULE_TEST05 calls TRIASYMQ for a quadrature rule of given order and region.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    27 June 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer DEGREE, the desired total polynomial degree
%    exactness of the quadrature rule.  0 <= DEGREE <= 50.
%
%    Input, integer NUMNODES, the number of nodes to be used by the rule.
%
%    Input, real VERT1(2), VERT2(2), VERT3(2), the
%    vertices of the triangle.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGLE_SYMQ_RULE_TEST05\n' );
  fprintf ( 1, '  Compute a quadrature rule for a triangle.\n' );
  fprintf ( 1, '  Check it by integrating orthonormal polynomials.\n' );
  fprintf ( 1, '  Polynomial exactness degree DEGREE = %d\n', degree );

  area = triangle_area ( vert1, vert2, vert3 );
%
%  Retrieve a symmetric quadrature rule.
%
  [ rnodes, weights ] = triasymq ( degree, vert1, vert2, vert3,  numnodes );
%
%  Construct the matrix of values of the orthogonal polynomials
%  at the user-provided nodes
%
  npols = ( ( degree + 1 ) * ( degree + 2 ) ) / 2;

  rints = zeros ( npols, 1 );

  for i = 1 : numnodes
    z(1) = rnodes(1,i);
    z(2) = rnodes(2,i);
    r = triangle_to_ref ( vert1, vert2, vert3, z );
    pols = ortho2eva ( degree, r );
    rints(1:npols) = rints(1:npols) + weights(i) * pols(1:npols);
  end

  scale = sqrt ( sqrt ( 3.0 ) ) / sqrt ( area );
  rints(1:npols) = rints(1:npols) * scale;

  d = ( rints(1) - sqrt ( area ) )^2 + sum ( rints(2:npols).^2 );
  d = sqrt ( d ) / npols;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  RMS integration error = %g\n', d );

  return
end


text:
function z_callBack_virtualTeachPendant_J7_plus()
%% Soft realtime joints control Joint X button callback
% Copyright: Mohammad SAFEEA, 27th-July-2018
global virtual_smartPad_w_joints_command;
disp('Moving J7 in positive direction');
virtual_smartPad_w_joints_command(7)=1;
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function D = dtiSelfDiffusionOfWater(temperature, pressure)
%
% D = dtiSelfDiffusionOfWater([temperature=37], [pressure=101.325])
%
% Returns the self-diffusion coefficient of water
% (in micrometers^2/millisecond) at the given temperature (in
% Centigrade) and pressure (in kPa). (These default to body temp
% and standard atmospheric pressure.)
%
% The function used to compute D is from Krynicki, Green & Sawyer
% (1978). Pressure and temperature dependence of self-diffusion in
% water. Faraday Discuss. Chem. Soc., 66, 199 - 208.
%
% The results also agree pretty well with the values presented in Mills
% (1973). Self-Diffusion in Normal and Heavy Water. JPhysChem 77(5), pg.
% 685 - 688. Mills estimates the self-diffusion of normal water at 35 deg C
% to be 2.919 micrometers^2/ms. This function returns 2.862 at that temp, a
% difference of about 2%.
%
% E.g.:
% x=[20:40];figure;plot(x,dtiSelfDiffusionOfWater(x));xlabel('T (C)');ylabel('D (\mum^2/ms)')
%
% See also:
%  http://www.lsbu.ac.uk/water/explan5.html
%
% Note that the relationship between diffusivity and temperature is nearly
% linear over the biologically useful range (e.g. 20-40 deg C) and is very
% well fit by a 3rd order polynomial. Thus, to find t (in Centigrade) given
% d and standard pressure, use:
%
% p=[0.6056 -6.6855 39.2757 -36.8468]; t = polyval(p,d);
%
%
% HISTORY:
% 2006.10.06 Bob Dougherty (RFD): Wrote it.

if(~exist('pressure','var')||isempty(pressure))
  P = 101.325;
else
  P = pressure;
end
if(~exist('temperature','var')||isempty(temperature))
  T = 273.15+37;
else
  T = temperature+273.15;
end

% D is in micrometers^2/millisec
% T is in Kelvins (C + 273.15)
% P is in kPa (100 kPa = 1 bar)
D = 12.5 .* exp(P*-5.22*10^-6) .* sqrt(T) .* exp(-925.*exp(P.*-2.6.*10^-6)./(T-(95+P.*2.61.*10^-4)));
return;

fprintf('D = %0.2f micrometers^2/msec (@ %0.1f C & %0.1f kPa).\n',D,T-273.15,P);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [G, cliques, fill_ins] = triangulate(G, order)
% TRIANGULATE Ensure G is triangulated (chordal), i.e., every cycle of length > 3 has a chord.
% [G, cliques, fill_ins, cliques_containing_node] = triangulate(G, order)
%
% cliques{i} is the i'th maximal complete subgraph of the triangulated graph.
% fill_ins(i,j) = 1 iff we add a fill-in arc between i and j.
%
% To find the maximal cliques, we save each induced cluster (created by adding connecting
% neighbors) that is not a subset of any previously saved cluster. (A cluster is a complete,
% but not necessarily maximal, set of nodes.)

MG = G;
n = length(G);
eliminated = zeros(1,n);
cliques = {};
for i=1:n
  u = order(i);
  U = find(~eliminated); % uneliminated
  nodes = myintersect(neighbors(G,u), U); % look up neighbors in the partially filled-in graph
  nodes = myunion(nodes, u); % the clique will always contain at least u
  G(nodes,nodes) = 1; % make them all connected to each other
  G = setdiag(G,0);
  eliminated(u) = 1;

  exclude = 0;
  for c=1:length(cliques)
    if mysubset(nodes,cliques{c}) % not maximal
      exclude = 1;
      break;
    end
  end
  if ~exclude
    cnum = length(cliques)+1;
    cliques{cnum} = nodes;
  end
end

fill_ins = sparse(triu(max(0, G - MG), 1));

%assert(check_triangulated(G)); % takes 72% of the time!

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%  Copyright (c) 2014, Karen Simonyan
%  All rights reserved.
%  This code is made available under the terms of the BSD license (see COPYING file).

% Modified by Aruni Roy Chowdhury

function [ res, extra ] = evalBestThresh( scores, gt )
%EVALBESTTHRESH
% finds an optimal threshold - the threshold which maximises the accuracy

    % threshold scores and get the sign

    res = -Inf;
    extra.bestThresh = [];

    % thresh loop
    for i=1:numel(scores)

        curThresh = scores(i);
        class = 2 * (scores >= curThresh) - 1;

        % class-n accuracy
        acc = mean(class == gt);

        if acc > res

            res = acc;
            extra.bestThresh = curThresh;
        end
    end

    res = res * 100;

end

author: suhangpro

repo: mvcnn

save_path: github-repos/MATLAB/suhangpro-mvcnn


text:
function [sv,stat0]=ephpos(time,teph,sat,nav,iode,sv) %#ok

global glc gls
stat0=1; tt=1e-3; sv1=gls.sv;
[sys,~]=satsys(sat);

if sys==glc.SYS_GPS||sys==glc.SYS_GAL||sys==glc.SYS_BDS||sys==glc.SYS_QZS
%     [eph,stat]=searcheph(teph,sat,iode,nav);
    [eph,stat]=searcheph_h(teph,sat,-1,nav.eph);
    if stat==0,stat0=0;return;end
    sv=eph2pos(time,eph,sv);
    time=timeadd(time,tt);
    sv1=eph2pos(time,eph,sv1);
    sv.svh=eph.svh;
elseif sys==glc.SYS_GLO
%     [geph,stat]=searchgeph(teph,sat,iode,nav);
    [geph,stat]=searchgeph_h(teph,sat,-1,nav.geph);
    if stat==0,stat0=0;return;end
    sv=geph2pos(time,geph,sv);
    time=timeadd(time,tt);
    sv1=geph2pos(time,geph,sv1);
    sv.svh=geph.svh;
else
    stat0=0; return;
end

sv.vel=(sv1.pos-sv.pos)/tt;
sv.dtsd=(sv1.dts-sv.dts)/tt;

return

repo: GINav


text:
function varargout = symamd(X,variablesonly)
%SYMAMD (overloaded)

if nargin < 2
    variablesonly = 0;
end
 if isa(X,'blkvar')
    X = sdpvar(X);
 end

if variablesonly
    Z = X.basis(:,2:end);
else
    Z = X.basis;
end

Z = reshape(sum(abs(Z),2),X.dim(1),X.dim(2));
Z = Z~=0;
varargout{1}=symamd(Z);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [val] = optGeneFitness(rxn_vector_matrix, model, targetRxn, rxnListInput, isGeneList)
% The fitness function
%
% USAGE:
%
%    [val] = optGeneFitness(rxn_vector_matrix, model, targetRxn, rxnListInput, isGeneList)
%
% INPUTS:
%    rxn_vector_matrix:    reactions vectors in a matrix
%    model:                model structure
%    targetRxn:            target reactions
%    rxnListInput:         list of reactions
%    isGeneList:           bolean checking if it is a gene list
%
% OUTPUT:
%    val:                  fitness value

global MaxKnockOuts
%size(rxn_vector_matrix)

popsize = size(rxn_vector_matrix,1);
val = zeros(1,popsize);

%rxnGeneMat is a required field for this function if rxnList is a gene List, so if it does not exist,
%build it.
if isGeneList && ~isfield(model,'rxnGeneMat')
    model = buildRxnGeneMat(model);
end

for i = 1:popsize
    rxn_vector = rxn_vector_matrix(i,:);
    rxnList = rxnListInput(logical(rxn_vector));

    %see if we've done this before
    val_temp = memoize(rxn_vector);
    if ~ isempty(val_temp)
        val(i) = val_temp;
        continue;
    end

    % check to see if mutations is above the max number allowed
    nummutations = sum(rxn_vector);
    if nummutations > MaxKnockOuts
        continue;
    end

	% generate knockout.
    if isGeneList
        modelKO = deleteModelGenes(model, rxnList);
    else % is reaction list
        [isValidRxn,removeInd] = ismember(rxnList,model.rxns);
        removeInd = removeInd(isValidRxn);
        modelKO = model;
        modelKO.ub(removeInd) = 0;
        modelKO.lb(removeInd) = 0;
    end

    % find growthrate;
%     slnKO = optimizeCbModel(modelKO);
%     growthrate1 = slnKO.f; %max growth rate.
    if exist('LPBasis', 'var')
        modelKO.LPBasis = LPBasis;
    end

    [slnKO, LPOUT] = solveCobraLPCPLEX(modelKO, 0,1);
    LPBasis = LPOUT.LPBasis;
    growthrate = slnKO.obj;

    % check to ensure that GR is above a certain value
    if growthrate < .10
        continue;
    end

%    display('second optimization');
    % find the lowesest possible production rate (a hopefully high number)
    % at the max growth rate minus some set factor gamma (a growth rate slightly
    % smaller than the max). A positive value will eliminate solutions where the
    % production envelope has a vertical line at the max GR, a "non-unique"
    % solution. Set value to zero if "non-unique" solutions are not an issue.
    gamma = 0.01; % proportional to Grwoth Rate (hr-1), a value around 0.5 max.

    %find indicies of important vectors
    indBOF = find(modelKO.c);
    indTar = findRxnIDs(modelKO, targetRxn);
    % generate a model with a fixed max KO growth rate
    modelKOsetGR = modelKO;
    modelKOsetGR.lb(indBOF) = growthrate - gamma; % this growth rate is required as lb.
    modelKOsetGR.c = zeros(size(modelKO.c));
    modelKOsetGR.c(indTar) = -1; % minimize for this variable b/c we want to look at the very minimum production.

    % find the minimum production rate for the targeted reaction.

%     slnKOsetGR = optimizeCbModel(modelKOsetGR);
%     minProdAtSetGR1 = -slnKOsetGR.f;  % This should be a negative value b/c of the minimization setup, so -1 is necessary.

    if exist('LPBasis2', 'var')
        modelKOsetGR.LPBasis = LPBasis2;
    end

    [slnKOsetGR, LPOUT2] = solveCobraLPCPLEX(modelKOsetGR, 0,1);
    LPBasis2 = LPOUT2.LPBasis;
    minProdAtSetGR = -slnKOsetGR.obj;

    % objective function for optGene algorithm = val (needs to be a negative value, since it is
    % a minimization)
    val(i) = -minProdAtSetGR;
    % penalty for a greater number of mutations
    % val(i) = -minProdAtSetGR * (.98^nummutations);
    % select best substrate-specific productivity
    % val(i) = -minProdAtSetGR * (.98^nummutations) * growthrate;

    % check to prevent very small values from being considerered improvments
    if val(i) > -1e-3
        val(i) = 0;
    end

    memoize(rxn_vector, val(i));
end

return;

%% Memoize
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% MEMOIZE %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% internal function to speed things up.

function [value] =  memoize(gene_vector, value)
global HTABLE
hashkey = num2str(gene_vector);
hashkey = strrep(hashkey,' ',''); % cut out white space from string (more space efficient).

if nargin == 1
    value = HTABLE.get(hashkey);
    return;
else
    if HTABLE.size() > 10000
        HTABLE = java.util.Hashtable;  %reset the hashtable if more than 10,000 entries.
    end
    HTABLE.put(hashkey, value);
    value = [];
    return;
end
return

author: opencobra

repo: cobratoolbox


text:
function seed = i4_seed_advance ( seed )

%*****************************************************************************80
%
%% I4_SEED_ADVANCE "advances" the seed.
%
%  Discussion:
%
%    This routine implements one step of the recursion
%
%      SEED = ( 16807 * SEED ) mod ( 2^31 - 1 )
%
%    This version of the routine does not check whether the input value of
%    SEED is zero.  If the input value is zero, the output value will be zero.
%
%    If we repeatedly use the output of SEED_ADVANCE as the next input,
%    and we start with SEED = 12345, then the first few iterates are:
%
%         Input      Output
%          SEED        SEED
%
%         12345   207482415
%     207482415  1790989824
%    1790989824  2035175616
%    2035175616    77048696
%      77048696    24794531
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Second Edition,
%    Springer, 1987,
%    ISBN: 0387964673,
%    LC: QA76.9.C65.B73.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, December 1986, pages 362-376.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley, 1998,
%    ISBN: 0471134031,
%    LC: T57.62.H37.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, Number 2, 1969, pages 136-143.
%
%  Parameters:
%
%    Input, integer SEED, the seed value.
%
%    Output, integer SEED, the "next" seed.
%
  i4_huge = 2147483647;

  seed = floor ( seed );

  seed = mod ( seed, i4_huge );

  if ( seed < 0 )
    seed = seed + i4_huge;
  end

  k = floor ( seed / 127773 );

  seed = 16807 * ( seed - k * 127773 ) - k * 2836;

  if ( seed < 0 )
    seed = seed + i4_huge;
  end

  return
end

repo: jburkardt-m


text:
function [cd, score] = distance_metrics(model_gt, model, theta)

if nargin < 3
    theta = 0.1;
end

if ~isfield(model_gt, 'vtx')
    [model_gt.vtx] = model_gt.vertices;
    model_gt = rmfield(model_gt,'vertices');
    [model_gt.mesh] = model_gt.faces;
    model_gt = rmfield(model_gt,'faces');
end

% compare cad model
% normalize both models
mask = (model_gt.anchor~=0)&(model.anchor~=0);
mean_est = mean(model.vtx(model.anchor(mask), :), 1);
if isequal(size(model_gt.vtx,1), 36) && isequal(numel(model_gt.anchor), 8) % TODO
    model_gt.anchor(2) = 36;
    mean_gtr = mean(model_gt.vtx(model_gt.anchor(mask), :), 1);
else
    mean_gtr = mean(model_gt.vtx(model_gt.anchor(mask), :), 1);
end
model.vtx = bsxfun(@minus, model.vtx, mean_est);
model_gt.vtx = bsxfun(@minus, model_gt.vtx, mean_gtr);

std_est = mean(std(model.vtx(model.anchor(mask), :), 1, 1));
std_gtr = mean(std(model_gt.vtx(model_gt.anchor(mask), :), 1, 1));
model.vtx = model.vtx/std_est;
model_gt.vtx = model_gt.vtx/std_gtr;

R = align_models(model, model_gt);
model.vtx = model.vtx*R';

model = computeMeshInfo(model);
model_gt = computeMeshInfo(model_gt);

% compute distance from model to target
kdtree = KDTreeSearcher(model_gt.vtx);
[U,~] = surfProjection(model, model_gt, kdtree);
dist_model_target = sum(sum((model.vtx - U).^2, 2)) / size(model.vtx, 1);

dist_model_target_theta = sum(sum((model.vtx - U).^2, 2) > theta)...
    / size(model.vtx, 1);

% compute distance from target to model
kdtree = KDTreeSearcher(model.vtx);
[U,~] = surfProjection(model_gt, model, kdtree);
dist_target_model = sum(sum((model_gt.vtx - U).^2, 2)) / size(model_gt.vtx, 1);

dist_target_model_theta = sum(sum((model_gt.vtx - U).^2, 2) > theta) ...
    / size(model_gt.vtx, 1);

cd = dist_model_target + dist_target_model;
score = dist_model_target_theta + dist_target_model_theta;

author: jhonykaesemodel


text:
function h = plotConfEllipse(thetaOpt,confStats)
%PLOTCONFELLIPSE  Plot Confidence Ellipse/Ellipsoid

scov = confStats.Cov;
nparam = length(confStats.ConfInt);
ndata = size(confStats.ConfBnds.bnds,1);
clim = confStats.Conf;
ci = confStats.ConfInt;

%Find confidence ellipse/ellipsoid
[Ve,D] = eig(scov); e = diag(D); dfe = ndata-nparam; %note update ndata
scale = nparam*(ndata-1)/dfe * rmathlib('qf',clim,nparam,dfe);

if(length(thetaOpt)==2)
    t = linspace(0,2*pi,5e2)';
    uxc = sqrt(e(1)*scale)*cos(t); uyc = sqrt(e(2)*scale)*sin(t);
    c = ((Ve')\[uxc, uyc]')';
    xc = c(:,1) + thetaOpt(1); yc = c(:,2) + thetaOpt(2);
    modestr = 'Ellipse';
    h(1) = plot(thetaOpt(1),thetaOpt(2),'r*');
    hold on
    h(2) = plot(xc,yc,'b');
    h(3) = plot([-ci(1) -ci(1) ci(1) ci(1) -ci(1)]+thetaOpt(1),[-ci(2) ci(2) ci(2) -ci(2) -ci(2)]+thetaOpt(2),'k');
    hold off
    str = sprintf('%g%%',clim*100);
    legend('Least Squares Optimal Solution',[str ' Confidence ',modestr],['Rectangular ' str ' Confidence Interval'],'location','nw');
else
    optiwarn('opti:plotEllipse','3D Not implemented yet');
%     u = linspace(0,2*pi,100)'; v = linspace(0,pi,length(u))';
%     [U,V] = meshgrid(u,v);
%     uxc = sqrt(e(1)*scale).*cos(U).*sin(V); uyc = sqrt(e(2)*scale).*sin(U).*cos(V); uzc = sqrt(e(2)*scale)*cos(V);
%
%     c = ((Ve')\[uxc(:) uyc(:) uzc(:)]')';
%     xc = c(:,1) + thetaOpt(1); yc = c(:,2) + thetaOpt(2); zc = c(:,3) + thetaOpt(3);
%
%     n = length(u);
%     xc = reshape(xc,n,n);
%     yc = reshape(yc,n,n);
%     zc = reshape(zc,n,n);
%
%     h(2) = mesh(xc,yc,zc);
    modestr = 'Ellipsoid';
%     hold on
%         h(1) = plot3(thetaOpt(1),thetaOpt(2),thetaOpt(3),'r*');
%     hold off;
end
title(sprintf('Optimal Parameter %g%% Confidence %s',clim*100,modestr));
xlabel('\theta_1'); ylabel('\theta_2'); zlabel('\theta_3');


text:
function [J,J_od,J_id,J_bl] = jdegree(CIJ)
%JDEGREE    Joint degree distribution
%
%   [J,J_od,J_id,J_bl] = jdegree(CIJ);
%
%   This function returns a matrix in which the value of each element (u,v)
%   corresponds to the number of nodes that have u outgoing connections
%   and v incoming connections.
%
%   Input:      CIJ,    directed (weighted/binary) connection matrix
%
%   Outputs:    J,      joint degree distribution matrix (shifted by one)
%               J_od,   number of vertices with od>id.
%               J_id,   number of vertices with id>od.
%               J_bl,   number of vertices with id=od.
%
%   Note: Weights are discarded.
%
%
%   Olaf Sporns, Indiana University, 2002/2006/2008

% ensure CIJ is binary...
CIJ = double(CIJ~=0);

N = size(CIJ,1);

id = sum(CIJ,1);    % indegree = column sum of CIJ
od = sum(CIJ,2)';   % outdegree = row sum of CIJ

% Create the joint degree distribution matrix
% Note:  the matrix is shifted by one, to accomodate zero id and od in the first row/column.
% Upper triangular part of the matrix has vertices with an excess of
%    outgoing edges (od>id)
% Lower triangular part of the matrix has vertices with an excess of
%    outgoing edges (id>od)
% Main diagonal has units with id=od

szJ = max(max(id,od))+1;
J = zeros(szJ);

for i=1:N
   J(id(i)+1,od(i)+1) = J(id(i)+1,od(i)+1) + 1;
end;

J_od = sum(sum(triu(J,1)));
J_id = sum(sum(tril(J,-1)));
J_bl = sum(diag(J));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef ConstraintManager < handle

    properties (SetAccess = private)
        C_equality;
        C_inequality;
    end

    methods (Access = public)

        function obj = ConstraintManager()
            obj.C_equality = containers.Map();
            obj.C_inequality = containers.Map();
        end

        function add_eq_constraint(obj, key, C_eq1, C_eq2)
            obj.C_equality(key) =  {C_eq1, C_eq2};
        end

        function add_ineq_constraint(obj, key, C_ineq1, C_ineq2)
            obj.C_inequality(key) =  {C_ineq1, C_ineq2};
        end

        function remove_eq_constraint(obj, key)
            remove(obj.C_equality, key);
        end

        function remove_ineq_constraint(obj, key)
            remove(obj.C_inequality, key);
        end

        function [C_eq1, C_eq2] = combine_all_eq_constraints(obj)
            C_eq_array = values(obj.C_equality);
            C_eq1 = [];
            C_eq2 = [];
            for i = 1:numel(C_eq_array)
                C_eq1 = [C_eq1; C_eq_array{i}{1}];
                C_eq2 = [C_eq2; C_eq_array{i}{2}];
            end
        end

        function [C_ineq1, C_ineq2] = combine_all_ineq_constraints(obj)
            C_ineq_array = values(obj.C_inequality);
            C_ineq1 = [];
            C_ineq2 = [];
            for i = 1:numel(C_ineq_array)
                C_ineq1 = [C_ineq1; C_ineq_array{i}{1}];
                C_ineq2 = [C_ineq2; C_ineq_array{i}{2}];
            end
        end
    end
end

author: HiroIshida

repo: robust-tube-mpc


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 = Stunorm4(x)
    y = 1/sqrt(2)*InverseCDF4(x);
end


text:
function [err, Itrans] = AFNI_IndexChange (Info, Iorig, Direction, DispOrient)
%
%   [err, Itrans] = AFNI_IndexChange (Info, Iorig, Direction, [DispOrient])
%
%Purpose:
%   Change the AFNI Index system between the display's and AFNI's coordinate system
%
%
%Input Parameters:
%   Info: The data structure output from BrikInfo
%   Iorig : anNx3 matrix containing the Ix Iy Iz indices of N points
%   Direction : A string : either 'A2D' or 'D2A' meaning Afni to Display or vice versa
%   DispOrient : The orientation string (or vector) for AFNI's display. This parameter is optional
%    and the defualt is 'RAI", it could be 'LAI' if AFNI's using Left is Left option.
%
%Output Parameters:
%   err : 0 No Problem
%       : 1 Mucho Problems
%   Itrans: if Iorig is specified, Itrans is Iorig in the new coordinate system
%
%
%
%Key Terms:
%
%More Info :
%   Test_AFNI_IndexChange
%
%   see also AFNI_CoordChange
%
%
%     Author : Ziad Saad
%     Date : Fri Sep 8 12:21:01 PDT 2000
%     LBC/NIMH/ National Institutes of Health, Bethesda Maryland

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

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;
Itrans= [];
maplocation = [0 0 0];
mapsign = [0 0 0];

if (nargin == 3),
	DispOrient = [0 3 4]; %that's RAI
end

%Assume we're going from Display to Afni ('D2A)
if (ischar(DispOrient)),
	[err, OrCode] = AFNI_OrientCode (DispOrient);
else
	OrCode = DispOrient;
end

	[err, TrCode] = AFNI_OrientCode (Info.Orientation(:,1)');

%check if that's what is required
if (strmatch(Direction, 'A2D')),
	%we're going from Afni to Display coordinates
		tmp = TrCode;
		TrCode = OrCode;
		OrCode = tmp;
else
	if (~strmatch(Direction, 'D2A')),
		err = ErrEval(FuncName,'Err_Bad Direction string');
		return;
	end
end

%get maplocation and mapsign, automatically from AFNI_CoordChange (no need to rewrite the code here)
	[err,maplocation, mapsign] = AFNI_CoordChange (OrCode, TrCode);

Itrans = Iorig;
if (strmatch(Direction, 'A2D')),
		fprintf ('Doing A2D\n');
		%(1-mapsign(i))./2 is 0 when mapsign(i) = 1 and 1 when mapsign(i) = -1;  This way, the if condition for mapsign(i) can be done away with
		%I left the simple method for the second loop for clarity. i don't think there's much efficiency difference between the two.
	for (i=1:1:3),
		Itrans(:,i) = ( (1-mapsign(i))./2 .* (Info.DATASET_DIMENSIONS(maplocation(i)) - 1) ) + ( mapsign(i) .* Iorig(:,maplocation(i)) );
	end
else
		fprintf ('Doing D2A: \n');
	for (i=1:1:3),
		Itrans(:,i) = Iorig(:,maplocation(i));
		if (mapsign(i) < 0),
			Itrans(:,i) = Info.DATASET_DIMENSIONS((i)) -1 - Itrans(:,i);
		end
	end
end

err = 0;
return;


text:
function [J, grad] = costFunction(theta, X, y)
%COSTFUNCTION Compute cost and gradient for logistic regression
%   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
%   parameter for logistic regression and the gradient of the cost
%   w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

h = sigmoid(X*theta);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Note: grad should have the same dimensions as theta
%

J = (-y'*log(h)-(1-y)'*log(1-h))/m;

grad = X'*(h-y)/m;

% =============================================================

end

repo: MachineLearningEx


text:
function r = cot(a)
%COT          Taylor cotangent  cot(a)
%
%Thanks to George Corliss for providing the Taylor expansion
%

% written  06/03/09     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                   % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  K = getappdata(0,'INTLAB_TAYLOR_ORDER');

  r = a;
  ct = a.t;
  N = size(a.t,2);
  r.t(1,:) = cot(a.t(1,:));
  ct(1,:) = 1+r.t(1,:).^2;     % 1+a^2
  r.t(2,:) = - ct(1,:) .* a.t(2,:);
  for i=2:K
    ct(i,:) = sum( r.t(1:i,:).*r.t(i:-1:1,:) , 1 );
    r.t(i+1,:) = - sum( ct(1:i,:).*a.t(i+1:-1:2,:).*repmat((i:-1:1)',1,N) , 1 )/i;
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


function [testdata, target_interf_echo] = simulateDataImageRIR(sirdb, serdb, rt60, repeat, seed)
%
% Generate test data with simulated rir.
% sirdb:                specified SIR (dB)
% serdb:                specified SER (dB)
% rt60:                 reverberation time (s)
% testdata:             the generated data
% target_interf_echo:   target, interference, and echo signals
%

addpath('RIR-Generator-master');
% to reproduce experiments
rng(seed);

%% base directories
base_target='target/*.wav';
base_interf='interference/*.wav';
base_ref='reference/*.wav';

%% load sources
%
% target
%
sub=dir(base_target);
idx=randperm(size(sub, 1), 1);
target=audioread([sub(idx).folder, '/', sub(idx).name]);
siglen=length(target);
% set target signal activity
target=[zeros(siglen, 1);target;target;zeros(siglen, 1)];

%
% interference
%
sub=dir(base_interf);
idx=randperm(size(sub, 1), 1);
interf=audioread([sub(idx).folder, '/', sub(idx).name]);
% set interference activity
begin=randi(length(interf)-4*siglen);
interf=interf(begin+1:begin+4*siglen);

%
% reference
%
sub=dir(base_ref);
idx=randperm(size(sub, 1), 1);
ref=audioread([sub(idx).folder, '/', sub(idx).name]);
% set reference activity
begin=randi(length(ref)-2*siglen);
ref=ref(begin+1:begin+2*siglen);
ref=[zeros(length(ref), 1);ref];

% repeat the same utterance
if repeat > 1
    target_bk=target; interf_bk=interf; ref_bk=ref;
    for i=1:repeat-1
        target = [target;target_bk];
        interf = [interf;interf_bk];
        ref = [ref;ref_bk];
    end
end

%% generate simulated rirs
room_bound_x=[3, 6];
room_bound_y=[4, 8];
room_bound_z=[2.5, 4];
% sample the room size
room_x=round(rand*(room_bound_x(2)-room_bound_x(1))+room_bound_x(1), 2);
room_y=round(rand*(room_bound_y(2)-room_bound_y(1))+room_bound_y(1), 2);
room_z=round(rand*(room_bound_z(2)-room_bound_z(1))+room_bound_z(1), 2);

c = 340;                    % Sound velocity (m/s)
fs = 16000;                 % Sample frequency (samples/s)
L = [room_x, room_y, room_z];  % Room dimensions [x y z] (m)
beta = rt60;                % Reverberation time (s)
n = 4096;                   % Number of samples
mtype = 'omnidirectional';  % Type of microphone
order = -1;                 % -1 equals maximum reflection order!
dim = 3;                    % Room dimension
orientation = 0;            % Microphone orientation (rad)
hp_filter = 1;              % Enable high-pass filter
early_reverb = fs * 0.05;   % define early reverberation

%
% generate mic position
%
d = 0.1;                    % mic spacing (m)
theta=rand(1, 1)*2*pi;
R=[cos(theta), -sin(theta), 0; sin(theta), cos(theta), 0; 0, 0, 1];
m1=R*[-d/2; 0; 0];
m2=R*[d/2; 0; 0];

space=0.5;
r0=space*ones(1, 3)+rand(1, 3).*(L-2*space*ones(1, 3));
r = [r0+m1'; r0+m2'];       % Receiver positions [x_1 y_1 z_1 ; x_2 y_2 z_2] (m)

%
% generate source position
%
dist_bound = [0.2, 3];  % allowed source-receiver distance
while 1
    elevation=asin(2*rand - 1);
    azimuth=2*pi*rand;
    % favor large distance
    radii=dist_bound(1)+(dist_bound(2)-dist_bound(1))*(rand.^(1/3));
    [offset_x, offset_y, offset_z] = sph2cart(azimuth, elevation, radii);
    offset_xyz = [offset_x offset_y offset_z];
    s1 = offset_xyz + r0;
    % check if the source position is within the correct range, otherwise resample
    if s1 <= (L-space) & s1 >= space
        break
    end
end
while 1
    elevation=asin(2*rand - 1);
    azimuth=2*pi*rand;
    % favor large distance
    radii=dist_bound(1)+(dist_bound(2)-dist_bound(1))*(rand.^(1/3));
    [offset_x, offset_y, offset_z] = sph2cart(azimuth, elevation, radii);
    offset_xyz = [offset_x offset_y offset_z];
    s2 = offset_xyz + r0;
    % check if the source position is within the correct range, otherwise resample
    if s2 <= (L-space) & s2 >= space
        break
    end
end

% % show mic-source positions
% plot(r(:,1),r(:,2),'o'); xlim([0,L(1)]); ylim([0,L(2)]);
% hold on; plot(s1(1),s1(2), 'x'); plot(s2(1),s2(2), 'v'); hold off;

% loudspeaker is 15 cm below mic
s3=r0;
s3(3)=s3(3)-0.15;

%
% generate rir
%
h=rir_generator(c, fs, r, s1, L, beta, n, mtype, order, dim, orientation, hp_filter);
a11=h(1, :);
a21=h(2, :);
% tt=0:length(a11)-1;
% plot(tt, a11, tt, a21);

h=rir_generator(c, fs, r, s2, L, beta, n, mtype, order, dim, orientation, hp_filter);
a12=h(1, :);
a22=h(2, :);

h=rir_generator(c, fs, r, s3, L, beta, n, mtype, order, dim, orientation, hp_filter);
b1=h(1, :);
b2=h(2, :);
% tt=0:length(b1)-1;
% plot(tt, b1, tt, b2);

%% mix data
%
% generate data image
%
target_img=[conv(target, a11), conv(target, a21)];
target_img(length(target)+1:end, :)=[];
target_early_reverb=conv(target, a11(1:early_reverb));
target_early_reverb(length(target)+1:end, :)=[];

interf_img=[conv(interf, a12), conv(interf, a22)];
interf_img(length(interf)+1:end, :)=[];
interf_early_reverb=conv(interf, a12(1:early_reverb));
interf_early_reverb(length(interf)+1:end, :)=[];

echo=[conv(ref, b1), conv(ref, b2)];
echo(length(ref)+1:end, :)=[];

%
% adjust snr
%
[target_img, gain]=adjustSNR(target, target_img, 0);
target_early_reverb=target_early_reverb*gain;
[interf_img, gain]=adjustSNR(target_img, interf_img, sirdb);
interf_early_reverb=interf_early_reverb*gain;
[echo, gain]=adjustSNR(target_img, echo, serdb);

% generate the mixed signals
testdata=[target_img+interf_img+echo, ref];
target_interf_echo=[target_early_reverb, interf_early_reverb, echo(:, 1)];

end


text:
function [ mapX, mapY ] = WarpRegion( xmin, ymin, mask, triX, coeffs )
%WARPREGION Summary of this function goes here
%   Detailed explanation goes here

    %%
    [h, w] = size(mask);
    mapX = zeros(size(mask));
    mapY = zeros(size(mask));

    ys = [1:h]' * ones(1, w) + ymin - 1;
    xs = ([1:w]' * ones(1, h))' + xmin - 1;

    for t=0:size(coeffs,1)-1

        trimap = triX == t;

        a = coeffs(t+1,:);

        xo = a(1) + a(2) * xs + a(3) * ys;

        mapX(trimap) = xo(trimap);

        yo = a(4) + a(5) * xs + a(6) * ys;
        mapY(trimap) = yo(trimap);

    end

    mapX(~mask) = -1;
    mapY(~mask) = -1;

    %%
%     [h, w] = size(mask);
%     mapX_2 = zeros(size(mask));
%     mapY_2 = zeros(size(mask));
%
%     ys = [1:h]' * ones(1, w) + ymin - 1;
%     xs = ([1:w]' * ones(1, h))' + xmin - 1;
%
%     ys = ys(:);
%     xs = xs(:);
%
%     xos = coeffs(1,:) + bsxfun(@times, coeffs(2,:), xs) + bsxfun(@times, coeffs(3,:), ys);
%     yos = coeffs(4,:) + bsxfun(@times, coeffs(5,:), xs) + bsxfun(@times, coeffs(6,:), ys);
%
%     maps = repmat(trimap(:),1, size(coeffs,1));
%     maps = repmat

end

author: TadasBaltrusaitis

repo: OpenFace


text:
 function y = proj_subset(ob, x, ia_start, ia_inc)
%function y = proj_subset(ob, x, ia_start, ia_inc)
%	y = G(:,ia, :,:) * x
%	caution: ia_start is from 1, not 0!

%if ob.is.empty, error empty, end
if ob.apower ~= 1, error 'power not done', end
if ob.is_transpose, error 'transpose not done', end

if ob.is_masked
	x = embed(x, ob.mask);
end

y = wtfmex('dsc,proj', ob.arg', single(x), ob.mask, ...
	int32(ia_start-1), int32(ia_inc), int32(ob.nthread), ob.chat);

y = y(:,ia_start:ia_inc:ob.na);

y = ob.scale * double(y(:));

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function J = computeCost(X, y, theta)
%COMPUTECOST Compute cost for linear regression
%   J = COMPUTECOST(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta
%               You should set J to the cost.

J = sum((X*theta - y).^2) / (2 * m);

% =========================================================================

end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
function [all_opts,use_libigl_static_library] = gptoolbox_mexopts(varargin)

  use_libigl_static_library = [];
  debug = false;
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Static','Debug'}, ...
    {'use_libigl_static_library','debug'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  verbose = false;
  MEXOPTS={'-largeArrayDims','-DMEX'};
  if verbose
    MEXOPTS={'-v',MEXOPTS{:}};
  end
  MSSE42='CXXFLAGS=$CXXFLAGS -msse4.2';
  STDCPP11='CXXFLAGS=$CXXFLAGS -std=c++11';
  try
    ELTOPO_INC= sprintf('-I%s/',path_to_eltopo);
    ELTOPO_LIB= strsplit(sprintf('-L%s/eltopo3d -leltopo_release',path_to_eltopo));
  catch e
    ELTOPO_INC='-DELTOPO_NOT_FOUND';
    ELTOPO_LIB={'-DELTOPO_NOT_FOUND'};
  end
  CLANG={'CXX=/usr/bin/clang++','LD=/usr/bin/clang++'};
  FRAMEWORK_LDFLAGS='LDFLAGS=\$LDFLAGS -framework Foundation -framework AppKit -framework Accelerate';
  NOOPT_LDOPTIMFLAGS='LDOPTIMFLAGS="-O "';
  if debug
    MEXOPTS={'-g',MEXOPTS{:}};
  end

  % See libigl documentation. In short, Libigl is a header-only library by
  % default: no compilation needed (like Eigen). There's an advanced **option**
  % to precompile libigl as a static library. This cuts down on compilation time.
  % It is optional and more difficult to set up. Set this to true only if you
  % know what you're doing.
  %use_libigl_static_library  = false;
  if isempty(use_libigl_static_library)
    use_libigl_static_library = exist([path_to_libigl '/lib/libigl.a'],'file')~=0;
  end
  if use_libigl_static_library
    r = -1;
    if strcmp(char(java.lang.System.getProperty('user.name')),'ajx')
      cmd = sprintf('cd "%s"/lib && /usr/local/bin/cmake ../optional/ -DLIBIGL_USE_STATIC_LIBRARY=ON -DLIBIGL_WITH_ANTTWEAKBAR=ON -DLIBIGL_WITH_BBW=ON -DLIBIGL_WITH_CORK=ON -DLIBIGL_WITH_CGAL=ON -DLIBIGL_WITH_COMISO=ON -DLIBIGL_WITH_EMBREE=ON -DLIBIGL_WITH_LIM=ON -DLIBIGL_WITH_MATLAB=ON -DLIBIGL_WITH_MOSEK=ON -DLIBIGL_WITH_NANOGUI=OFF -DLIBIGL_WITH_OPENGL=ON -DLIBIGL_WITH_PNG=ON -DLIBIGL_WITH_TETGEN=ON -DLIBIGL_WITH_TRIANGLE=ON -DLIBIGL_WITH_VIEWER=ON -DLIBIGL_WITH_XML=ON -DCMAKE_BUILD_TYPE=Release && ../scripts/make.sh -j3',path_to_libigl);
    else
      cmd = sprintf('cd "%s"/lib && /usr/local/bin/cmake ../optional/ -DLIBIGL_USE_STATIC_LIBRARY=ON -DLIBIGL_WITH_ANTTWEAKBAR=ON -DLIBIGL_WITH_BBW=ON -DLIBIGL_WITH_CORK=ON -DLIBIGL_WITH_CGAL=ON -DLIBIGL_WITH_COMISO=ON -DLIBIGL_WITH_EMBREE=ON -DLIBIGL_WITH_LIM=ON -DLIBIGL_WITH_MATLAB=ON -DLIBIGL_WITH_MOSEK=ON -DLIBIGL_WITH_NANOGUI=OFF -DLIBIGL_WITH_OPENGL=ON -DLIBIGL_WITH_PNG=ON -DLIBIGL_WITH_TETGEN=ON -DLIBIGL_WITH_TRIANGLE=ON -DLIBIGL_WITH_VIEWER=ON -DLIBIGL_WITH_XML=ON -DCMAKE_BUILD_TYPE=Release && make -j3',path_to_libigl);
    end
    if verbose
      fprintf('%s\n',cmd);
    end
    [r,c] = system(cmd);
    if r ~= 0
      warning('libigl Make error:\n%s',c);
      use_libigl_static_library = false;
    end
  end

  %use_libigl_static_library= false;
  LIBIGL_INC=sprintf('-I%s/include',path_to_libigl);
  if use_libigl_static_library
    LIBIGL_FLAGS='-DIGL_STATIC_LIBRARY';
    LIBIGL_LIB=strsplit(sprintf('-L%s/lib -ligl',path_to_libigl));
    LIBIGL_LIBEMBREE='-ligl_embree';
    LIBIGL_LIBMATLAB='-ligl_matlab';
    LIBIGL_LIBCGAL='-ligl_cgal';
    LIBIGL_LIBCORK='-ligl_cork';
  else
    % `mex` has a silly requirement that arguments be non-empty, hence the NOOP
    % defines
    LIBIGL_FLAGS='-DIGL_SKIP';
    LIBIGL_LIB={'-DIGL_SKIP'};
    LIBIGL_LIBMATLAB='-DIGL_SKIP';
    LIBIGL_LIBEMBREE='-DIGL_SKIP';
    LIBIGL_LIBCGAL='-DIGL_SKIP';
    LIBIGL_LIBCORK='-DIGL_SKIP';
  end
  LIBIGL_BASE={LIBIGL_INC,LIBIGL_FLAGS,LIBIGL_LIB{:},LIBIGL_LIBMATLAB};

  if exist(sprintf('%s/external/nanogui/ext/eigen/',path_to_libigl))
    EIGEN_INC=sprintf('-I%s/external/nanogui/ext/eigen/',path_to_libigl);
  elseif exist('/usr/local/include/eigen3')
    EIGEN_INC='-I/usr/local/include/eigen3';
  elseif exist('/opt/local/include/eigen3')
    EIGEN_INC='-I/opt/local/include/eigen3';
  end

  EMBREE=[path_to_libigl '/external/embree'];
  EMBREE_INC=strsplit(sprintf('-DWITH_EMBREE -I%s -I%s/include/',EMBREE,EMBREE));
  %EMBREE_LIB=strsplit(sprintf('-ltbb -L%s -lembree_avx -lembree_avx2 -lembree_sse42 -limage -llexers -llights -llights_ispc -lnoise -lnoise_ispc -lscenegraph -lsimd -lsys -ltasking -ltexture -ltexture_ispc -ltutorial -ltutorial_device -ltutorial_device_ispc ',[path_to_libigl '/external/embree/build/']));
  EMBREE_LIB=strsplit(sprintf('-L%s -lembree',[path_to_libigl '/external/embree/build/']));

  CORK_INC=sprintf('-I%s/src',[path_to_libigl '/external/cork']);
  CORK_LIB=strsplit(sprintf('-L%s -lcork',[path_to_libigl '/lib']));

  TINYXML2=[path_to_libigl '/external/tinyxml2'];
  TINYXML2_INC=sprintf('-I%s/',TINYXML2);
  TINYXML2_LIB=strsplit(sprintf('-L%s/ -ltinyxml2',TINYXML2));

  if exist('/usr/local/include/CGAL')
    CGAL='/usr/local/';
  elseif exist('/opt/local/include/CGAL')
    CGAL='/opt/local/';
  end
  CGAL_INC=sprintf('-I%s/include',CGAL);
  CGAL_LIB=strsplit(sprintf('-L%s/lib -lCGAL -lCGAL_Core -lgmp -lmpfr',CGAL));
  CGAL_FLAGS='CXXFLAGS=\$CXXFLAGS -frounding-math';

  BOOST='/usr/local/';
  BOOST_INC=sprintf('-I%s/include',BOOST);
  BOOST_LIB=strsplit(sprintf('-L%s/lib -lboost_thread-mt -lboost_system-mt',BOOST));

  all_opts = { ...
         MEXOPTS{:},...
         MSSE42, ...
         STDCPP11, ...
         EIGEN_INC, ...
         ELTOPO_INC, ...
         ELTOPO_LIB{:}, ...
         FRAMEWORK_LDFLAGS, ...
         LIBIGL_BASE{:}, ...
         LIBIGL_LIBEMBREE, ...
         LIBIGL_LIBCGAL, ...
         LIBIGL_LIBCORK, ...
         EMBREE_INC{:}, EMBREE_LIB{:}, ...
         CORK_INC,CORK_LIB{:}, ...
         TINYXML2_INC,TINYXML2_LIB{:}, ...
         CGAL_INC,CGAL_LIB{:},CGAL_FLAGS, ...
         BOOST_INC,BOOST_LIB{:}, ...
         };
end

author: alecjacobson

repo: gptoolbox


text:
function [vacc, hacc, vcm, hcm, pg] = crfTestCV3(imsegs, pvSP, phSP, adjlist, pE, crfw, param, edgelen, priors)

ncv = numel(crfw);
nimages = numel(imsegs);

for f = 1:nimages
    disp(num2str(f))
    c = ceil(f/nimages*ncv);
    pg{f} = crfTestImage(pvSP{f}, phSP{f}, adjlist{f}, pE{f}, crfw{c}, param{c}, edgelen{f}, priors{c});
    [tmpvacc2, tmphacc2] = mcmcProcessResult(imsegs(f), pg(f));
    [tmpvacc1, tmphacc1] = mcmcProcessResult(imsegs(f), {[pvSP{f}(:, 1) repmat(pvSP{f}(:, 2), 1, 5).*phSP{f} pvSP{f}(:, 3)]});
    disp([num2str(tmpvacc1) '-->' num2str(tmpvacc2) '   ' num2str(tmphacc1) '-->' num2str(tmphacc2)])
    im = rgb2gray(im2double(imread(['../images/all_images/' imsegs(f).imname])));
    figure(1), displayLabels(imsegs(f), pvSP{f}, im);
    figure(2), displayLabels(imsegs(f), pg{f}, im);
    drawnow;
    pause(1)
end
[vacc, hacc, vcm, hcm] = mcmcProcessResult(imsegs, pg);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function displayLabels(imsegs, pg, im)
if size(pg, 2)==7
    v2 = pg(:, 1);
    v1 = sum(pg(:, 2:6), 2);
    v3 = pg(:, 7);
else
    v2 = pg(:, 1);
    v1 = pg(:, 2);
    v3 = pg(:, 3);
end
lim = cat(3, v1(imsegs.segimage), v2(imsegs.segimage), v3(imsegs.segimage));
imagesc(lim*0.5 + repmat(im, [1 1 3])*0.5), axis image

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

function pg = crfTestImage(pvSP, phSP, adjlist, pE, crfw, param, edgelen, crfpriors)

% get f1(y1, y2) = log(P(y1|x)) + log(P(y2|x))
% and f2(y1, y2) = I(y1=y2)*log(P(y1=y2|x)) + I(y1~=y2)*log(P(y1~=y2|x))

npair = size(adjlist, 1);

f1 = zeros(npair, 7, 7);
f2 = zeros(npair, 7, 7);
f3 = zeros(npair, 7, 7);

pg = [pvSP(:, 1)  repmat(pvSP(:, 2), 1, 5).*phSP  pvSP(:, 3)];
s1 = adjlist(:, 1);
s2 = adjlist(:, 2);

nb = zeros(size(pg, 1), 1);
edgeperc = zeros(numel(s1), 2);
nb = zeros(size(pg, 1), 1);
for k = 1:numel(s1)
    nb(s1(k)) = nb(s1(k))+1;
    nb(s2(k)) = nb(s2(k))+1;
    edgeperc(k, 1) = edgelen(s1(k), s2(k)) / (sum(edgelen(s1(k), :))+sum(edgelen(:, s1(k))));
    edgeperc(k, 2) = edgelen(s1(k), s2(k)) / (sum(edgelen(s2(k), :))+sum(edgelen(:, s2(k))));
end

for k1 = 1:7
    for k2 = 1:7
        kp = param(k1 + (k2-1)*7, :);
        f1(:, k1, k2) = log(pg(s1, k1))./nb(s1) + log(pg(s2, k2))./nb(s2);
        if k1==k2
            f2(:, k1, k2) = ...
                log(kp(3)./(1+exp(-kp(1)-kp(2)*(log(pE)-log(1-pE))))).*...
                (edgeperc(:, 1) + edgeperc(:, 2));
            %f2(:, k1, k2) = log(1);
        else
            f3(:, k1, k2) = ...
                log(kp(3)./(1+exp(-kp(1)-kp(2)*(log(pE)-log(1-pE))))).*...
                (edgeperc(:, 1) + edgeperc(:, 2));
            %f2(:, k1, k2) = log(0.000001);
        end
    end
end

nsp = size(pvSP, 1);
adjmat = zeros(nsp, nsp);
potentials = cell(nsp, nsp);

evidences = 1/7*ones(7, nsp); % exp(crfw(1)*log(pg'));
%evidences = exp(log(pg'));

%f2 = f2./nb(:, 1).^crfw(end) + f2./nb(:, 2).^crfw(end);
%f3 = f3./nb(:, 1).^crfw(end) + f3./nb(:, 2).^crfw(end);

for k = 1:numel(s1)
    adjmat(s1(k), s2(k)) = 1;
    adjmat(s2(k), s1(k)) = 1;
    potentials{s1(k), s2(k)} = 0.5*squeeze(exp(crfw(1)*f1(k, :, :) + crfw(2)*f2(k, :, :) + crfw(3)*f3(k, :, :) + reshape(crfpriors, [1 7 7])));
    potentials{s2(k), s1(k)} = potentials{s1(k), s2(k)}';
end

% mrf2 = mk_mrf2(adjmat, potentials);
% engine = belprop_mrf2_inf_engine(mrf2);
% engine = enter_soft_evidence(engine, evidences);
%
% pg = zeros(nsp, 7);
% for s = 1:nsp
%     disp(marginal_nodes(engine, s))
%     pg(s, :) = marginal_nodes(engine, s)';
% end

[pg, niter] = bp_mrf2_old(adjmat, potentials, evidences, 'max_iter', 50);
pg = pg';

author: Cloud-CV

repo: object-proposals


text:
function value = f2 ( x )

%*****************************************************************************80
%
%% F2 evaluates a function of a scalar used in defining P2(Q).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ian Sloan, Stephen Joe,
%    Lattice Methods for Multiple Integration,
%    Oxford, 1994,
%    ISBN: 0198534728,
%    LC: QA311.S56
%
%  Parameters:
%
%    Input, real X, the value of the argument.
%
%    Output, real VALUE, the value of F2(X).
%
  value = 1.0 + 2.0 * pi * pi * ( x * x - x + 1.0 / 6.0 );

  return
end

author: johannesgerer


text:
function [Fm,Vm]=cap_patchcylinder(F1,V1,F2,V2,nr,nz)

% function [Fm,Vm]=cap_patchcylinder(F1,V1,F2,V2,nr,nz)
% ------------------------------------------------------------------------
% This function assumes the inputs F1, V1 and F2, V2 define the faces and
% vertices of two concentrid cylinders and closes to top and bottom faces
% by connecting the cylinders together. The input nr defines the number of
% radial points to add for the connection. The input nz defines the number
% of steps used in the z-direction for the cylinders.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% Change log:
% 2014/09/25
% 2019/03/29 Updated help documentation
%
%------------------------------------------------------------------------

%%

Fm=[F1; (F2+size(V1,1))];
Vm=[V1; V2];

%Coordinates created from nr*nz matrix

%Top edge indices
It=ones(nr-1,1); Jt=(1:1:nr-1)';
INDt=sub2ind([nz,nr-1],It,Jt);
% Vt1=V1(INDt,:); Vt2=V2(INDt,:); Vt=[Vt1; Vt2];

%Bottom edge indices
Ib=nz.*ones(nr-1,1); Jb=(1:1:nr-1)';
INDb=sub2ind([nz,nr-1],Ib,Jb);
% Vb1=V1(INDb,:); Vb2=V2(INDb,:); Vb=[Vb1; Vb2];

Fp=(ones(nr-1,1)*[1 2 2 1])+(((1:1:nr-1*ones)'-1)*ones(1,4)); Fp(Fp==nr)=1; %Quad order

Fq=[INDt(Fp(:,1:2)) INDt(Fp(:,3:4))+size(V1,1)];
Ftt=[Fq(:,[1 3 2]); Fq(:,[1 4 3])]; %Tri order

Fq=[INDb(Fp(:,1:2)) INDb(Fp(:,3:4))+size(V1,1)];
Ftb=[Fq(:,[1 3 2]); Fq(:,[1 4 3])]; %Tri order

Fm=[Fm; Ftt; Ftb];

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% Name: tanner_LDPC.m
% Description: This script displays different parameters of an LDPC code
% given his A-list format file.
% Copyright (c) 2009. Robert Morelos-Zaragoza. All rights reserved.

% clear all
%
% name = input('enter file name: ');
name = 'NR_ldpc';
[H] = alist2sparse('NR_ldpc.alist')

figure(1)
spy(H)
title( strcat('H matrix of LDPC code with A-list file name: ', name))

n1 = size(H,1); n2 = size(H,2);

% Adjacency matrix as defined in Matlab
A = [ zeros(n2) H'; H zeros(n1) ];

xy = [];
for i=1:n2
    xy = [ xy; i 2 ];                 % Locations of code nodes
end
for i=1:n1
    xy = [ xy; i 1 ];                 % Locations of check nodes
end

figure(2)
gplot(A,xy)
axis([0 n2+1 0.9 2.1])
title( strcat('Tanner graph of LDPC code with A-list file name: ', name))

figure(3)
for i=1:n1
    cnd(i) = sum(H(i,:));
end
stem(cnd)
title('Check node degree distribution')

figure(4)
for j=1:n2
    vnd(j) = sum(H(:,j));
end
stem(vnd)
title('Variable node degree distribution')

fprintf('Density = %e\n',nnz(H)/(n1*n2));

author: xiaoshaoning


text:
%% Author: epokh
%% Website: www.epokh.org/drupy
%% This software is under GPL

%%This function plan a trajectory passing in N
%% points using N-1 3th order polynomial:
%% the user must choose: intial,final and intermediate points
%% timePoints is a matrix [3,N] where first row is the time,second row
%% is the position points, third row is the velocities
%% intermediate velocities (initial and final ones are imposed to be 0)
%% Output: a matrix [N-1,4] whose columns are the polynomail coefficients

function [TJ]=trajectory3Velconstraint(timePointsVels)
N=length(timePointsVels);

if(N<2)
fprintf('Error: you defined only 1 point');
return;
end
TJ=[];
%%set the initial and final velocities to 0
timePointsVels(3,1)=0;
timePointsVels(3,end)=0;
for k=1:N-1
    tk=timePointsVels(1,k);
    tk1=timePointsVels(1,k+1);
    qk=timePointsVels(2,k);
    qk1=timePointsVels(2,k+1);
    velk=timePointsVels(3,k);
    velk1=timePointsVels(3,k+1);
    coeffs=pol3interpol(tk,tk1,qk,qk1,velk,velk1);
    TJ=[TJ;coeffs'];
end

end

author: Sable

repo: mcbench-benchmarks


text:
function plotquadtfdist(p, varargin);
%PLOTQUADTFDIST Plot quadratic time-frequency distribution
%   Usage: plotquadtfdist(p);
%
%   'plotquadtfdist(p)' plots the quadratic time-frequency distribution
%   on the time-frequency plane. The quadratic time-frequency distribution
%   should be a square matrix.
%
%   PLOTQUADTFDIST takes the following additional arguments:
%
%     'dynrange',r
%              Limit the dynamical range to r by using a colormap in
%              the interval [chigh-r,chigh], where chigh is the highest
%              value in the plot. The default value of [] means to not
%              limit the dynamical range.
%
%     'db'     Apply 20*log_{10} to the coefficients. This makes
%              it possible to see very weak phenomena, but it might show
%              too much noise. A logarithmic scale is more adapted to
%              perception of sound. This is the default.
%
%     'dbsq'   Apply 10*log_{10} to the coefficients. Same as the
%              'db' option, but assume that the input is already squared.
%
%     'lin'    Show the coefficients on a linear scale. This will
%              display the raw input without any modifications. Only works for
%              real-valued input.
%
%     'linsq'  Show the square of the coefficients on a linear scale.
%
%     'linabs'  Show the absolute value of the coefficients on a linear scale.
%
%     'tc'     Time centring. Move the beginning of the signal to the
%              middle of the plot.
%
%     'clim',clim   Use a colormap ranging from clim(1) to clim(2). These
%                   values are passed to imagesc. See the help on imagesc.
%
%     'image'       Use imagesc to display the plot. This is the default.
%
%     'contour'     Do a contour plot.
%
%     'surf'        Do a surf plot.
%
%     'colorbar'    Display the colorbar. This is the default.
%
%     'nocolorbar'  Do not display the colorbar.
%
%     'display'     Display the figure. This is the default.
%

% AUTHOR: Jordy van Velthoven

complainif_notenoughargs(nargin, 1, 'PLOTQUADTFDIST');

[N1, N2] = size(p);

if N1==N2
 yr = [0, 1-2/N1];
else
 error('%s: The input should be a square matrix.', upper(mfilename));
end;

tfplot(p, 1, yr, varargin{:});

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [msh,lights,tenseMsh] = mrmConvertEMSEMesh(fileName,mmPerVox, host, id, varargin);
%
%  [msh,lights,tenseMsh] = mrmConvertEMSEMesh(fileName,mmPerVox, host, id, varargin);
%
%
% Author: ARW (based on RFD's mrmBuildMesh)
% Purpose:
%  Take a mesh in EMSE's .wfr format and convert it into a mrLoadRet /
%  mrMesh-type mesh. Allows you to do the usual mrMesh type things like
%  relaxation...
%  Additional processing options can be set as well.  These
%  options include:
%     * 'RelaxIterations'- the next value specifies how many extra smoothing iterations.
%     * 'SavetenseMsh'- if you specify any RelaxIterations and add this option, then
% you will get two meshes- the relaxed mesh and the unrelaxed (tense) mesh.
%
% Returns:
% The mesh structure and a lights structure.
% You can also get the unrelaxed mesh out (tenseMsh).
%
% See Also
%  mrmMapVerticesToGray
%
% Examples:
%   mrmConvertEMSEMesh('x:\anatomy\norcia\EMSE\test.wfr', view.mmPerVox, 'Background', [0.3,0.4,0.5]);
%   mrmBuildMesh('x:\anatomy\norcia\EMSE\test.wfr', view.mmPerVox, 'localhost', -1);
%   mrmBuildMesh('x:\anatomy\norcia\EMSE\test.wfr', view.mmPerVox, 'localhost', -1, ...
%                     'RelaxIterations', relaxIterations, 'SavetenseMsh');
%
% Notes:
%  2003.09.17 RFD: vertex-to-volume mapping is now done in a separate
%  function (mrmMapVerticesToGray), and we don't call it. So, the calling
%  function will need to compute that mapping and add the appropriate
%  fields to the mesh struct.

% transparency is off by default because it is slow.
mrGlobals;

summaryParams = 1;
meshName = '';
QueryFlag = 1;
relaxIter = 1;
saveTense = 0;
backColor = [0.2 0 0];  % CHanged from Gray - partly just to keep track of EMSE meshes...

if ieNotDefined('fileName'), error('EMSE format file is required.'); end
if ieNotDefined('mmPerVox'), error('mmPerVox is required.'); end
if ieNotDefined('host'), host = 'localhost'; end
if ieNotDefined('id'), id = 1; end

if(nargout>2),  saveTense = 1; end

% Parse the varargin values
for(ii=1:length(varargin))
    if    (strcmpi(varargin{ii}, 'RelaxIterations')), relaxIter = varargin{ii+1};
    elseif(strcmpi(varargin{ii}, 'QueryFlag')),       QueryFlag = varargin{ii+1};
    elseif(strcmpi(varargin{ii}, 'MeshName')),        meshName = varargin{ii+1};
    elseif (strcmpi(varargin{ii},'Background')),      backColor = varargin{ii+1};
    end
end

% Set initial parameters for the mesh.
msh = meshDefault(host,id,mmPerVox,relaxIter,meshName);

if QueryFlag,  msh = meshQuery(msh,summaryParams); end

% If the window is already open, no harm is done.
msh = mrmInitHostWindow(msh)

[msh, lights] = mrmInitMesh(msh,backColor);

disp('Importing mesh data from EMSE')
%

% ***************************************************************
% We need to construct a basic msh structure from the file data.
[vertex,face,edge,mesh]=mesh_emse2matlab2(fileName);

% Vertices come back rotated in an odd way. The entire cortex is rotated 90
% degrees about the L/R axis (i.e. the sag view is rotated 90degrees
% anticlockwise).
v=vertex(1:2,:);
v=v-128;

rotMat=[0 -1;1 0];
v=v'*rotMat;
vertex(1:2,:)=v'+128;
LR=vertex(3,:);

LR=(256-LR);
vertex(3,:)=LR;

face=face([3 2 1],:);

nVerts=length(vertex);

%vertex=vertex*10000;
p.vertices=(vertex);
p.triangles=face;
p.triangles = p.triangles - 1
p.class='mesh';

% Sometimes we pass in the whole array of meshes.  Mostly, just one,
% though.
host = meshGet(msh,'host');
windowID = meshGet(msh,'windowid');

if isempty(host), host = 'localhost'; end

if isempty(windowID), error('Mesh must specify a window'); end
%[id, status, result] = mrMesh ('localhost', windowID, 'add_actor', p)
p.scale = meshGet(msh,'mmPerVox');
p = setSmooth(p,msh,1);
p = setDecimate(p,msh,1);
p.actor = meshGet(msh,'actor');

p.colors=ones(4,length(p.vertices))*255;

[a,b,c]=mrmesh('localhost',id,'set_mesh',p)

        p.scale = meshGet(msh,'mmPerVox');
        p = setSmooth(p,msh,1);
        p = setDecimate(p,msh,1);
        %p.actor = meshGet(msh,'actor');
        p.normals=mrmGet(msh,'normals');
        p.actor = meshGet(msh,'actor');

%disp('Building smoothed and decimated mesh for display...');
%mrmSet(msh,'buildMeshAndDecimateAndSmooth',voxels);
%[msh] = mrmSet(msh,'smooth')

% This is a little 'center object' routine.  We could put this into mrmSet,
% really.
vertices = mrmGet(msh,'vertices');
mrmSet(msh,'origin',-mean(vertices'));

% Save these unsmoothed data.
unSmoothedData = mrmGet(msh,'data');
msh.initVertices=vertices;
msh.grayLayers=3;
view=getSelectedVolume;
    disp('Finding vertex to gray map');
    % You have to do this before smoothing...

v2gMap = mrmMapVerticesToGray(vertices,view.nodes,[1 1 1]);

msh.vertexGrayMap=v2gMap;
msh.grayToVertexMap = mrmMapGrayToVertices(view.nodes,vertices, [1 1 1]);

% Now also try to find the other mapping: the mapping from the mesh to all
% the gray nodes: So that we can ask: for any arbitray gray node, which
% mesh point is it closest to?

% If we smooth, the mesh, it is done here.
if(relaxIter>0)

    mrmSet(msh,'smooth');

    % We get the curvature colors from the uninflated mesh
    % We should figure out the correct value to use as our threshold. The mean
    % isn't ideal, since it is moved around by the large areas with arbitrary
    % curvature (eg. corpus callosum). What we really want is the value that
    % corresponds to zero curvature.
    % maybe make the specific curvature map colors adjustable?
    % We now use the actual curvature values, so we know that zero is zeros
    % curvature.
    disp('Setting up the curvature colors');

    % Attach curvature data to the mesh.  We turn on the color later, I think.
    msh=mrmSet(msh,'curvature');
    p.colors=[repmat(msh.curvature,3,1); ones(1,length(msh.curvature))]*255;
    %[a,b,c]=mrmesh('localhost',id,'set_mesh',p)

    curvColorIntensity = 128*meshGet(msh,'curvatureModDepth'); % mesh.curvature_mod_depth;

    monochrome = uint8(round((double(msh.curvature>0)*2-1)*curvColorIntensity+127.5));

    msh = mrmSet(msh,'colors',monochrome);

    disp('Storing the smoothed mesh data computed by mrMesh...');
    data = mrmGet(msh,'data');
    msh = meshSet(msh,'data',data);
    msh = meshSet(msh,'connectionMatrix',1);

else
    % We don't smooth the data with mrMesh.  We just assign it.
    msh = meshSet(msh,'data',unSmoothedData);
    msh = meshSet(msh,'connectionMatrix',1);
end

% In either case, we return a version of the data without smoothing.  This
% mesh is used to register with the gray coordinates.  This is a little
% disorganized.  If we don't smooth, tenseMesh is identically msh.
tenseMsh = msh;
tenseMsh = meshSet(tenseMsh,'data',unSmoothedData);

return;

%---------------------------------------
function p = setSmooth(p,mesh,val)
if val
    p.do_smooth = 1;
    p.smooth_iterations = meshGet(mesh,'smoothiterations');
    p.smooth_relaxation = meshGet(mesh,'smoothrelaxation');
    p.smooth_sinc_method = meshGet(mesh,'smoothmethod');
    p.do_smooth_pre = meshGet(mesh,'smooth_pre');
else
    p.do_smooth = 0;
    p.do_smooth_pre = 0;
end

return;

%----------------------------
function p = setDecimate(p,mesh,val)

if val
    p.do_decimate = 1;
    p.decimate_reduction = meshGet(mesh,'decimatereduction');
    p.decimate_iterations = meshGet(mesh,'decimateiterations');
else
    p.do_decimate = 0;
end

return;

author: vistalab

repo: vistasoft


text:
close all;
clearvars;
clc;
[LoD,HiD,LoR,HiR] = wfilters('db4');

load noisdopp;

% perform 4 level decomposition
[coefficients, levels] = wavedec(noisdopp,4, LoD,HiD);

% plot approximation and detail components
figure;
plot(coefficients); title('Coefficients');

export_fig images/noisdopp_db4_l4_decomposition.png -r120 -nocrop;

% reconstruct the signal
reconstructed = waverec(coefficients, levels, LoR, HiR);

% measure the maximum difference
max_abs_diff = max(abs(noisdopp-reconstructed))

figure;
for level=0:4
    level_app_coeffs = appcoef(coefficients, levels, LoR, HiR, level);
    subplot(511+level);
    plot(level_app_coeffs);
    title(sprintf('Approximation coefficients @ level-%d', level));
end

export_fig images/noisdopp_db4_l4_appcoeffs.png -r120 -nocrop;

author: indigits

repo: sparse-plex


text:
function rgb = gray2rgb(im,map,n)
  % GRAY2RGB
  %
  % rgb = gray2rgb(im,map,n)
  %
  % Use a colormap, map, to convert a grayscale/intensity image with n colors
  % to pseudocolor rgb image. The n colors in the original image are mapped
  % evenly to colors in the given colormap, not according to value.
  %
  % rgb = gray2rgb(im,map)
  %
  % Same as above but compute number of colors in im
  %
  % rgb = gray2rgb(im)
  %
  % Use default colormap and compute number of colors in im
  %
  % Inputs:
  %  im  grayscale image, 2d array (h,w)
  %  map  colormap
  %  n  (optional) number of colors in im, default is to compute unique colors,
  %    this is usually preferred but it's time consuming so if the value is
  %    already known then it can be given.
  % Outputs:
  %  rgb  color Red, Green, Blue image, 3d array (h,w,3)
  %
  % See also: rgb2gray, hsv2rgb, rgb2hsv, imshow
  %

  % number of unique colors
  if(~exist('n','var'))
    n = size(unique(im),1);
  end

  if(~exist('map','var'))
    if(~exist('gray2ind'))
      error( ...
        'Image processing toolbox seems not installed, try again without map')
    end
    if(isempty(get(0,'CurrentFigure')))
      map = jet(n);
    else
      map = colormap;
    end
    %map = jet(n);
    m = size(map,1);
  else
    m = size(map,1);
  end

  if(exist('gray2ind')&&n<=65536)

    in = gray2ind(im,n);

    if(strcmp(class(in),'uint8') && (n/(n/m)) > 256)
      in = uint16(in);
    end

    if( m <= 256 )
      rgb = label2rgb(idivide(in,n/m,'fix')+1,map);
    else
      rgb = ind2rgb(idivide(in,n/m,'fix'),map);
    end

    % if it's not going to change the result it's faster to convert to and
    % from uint8 and use label2rgb
    if(strcmp(class(im),'double'))
      rgb = im2double(rgb);
    elseif(strcmp(class(im),'uint16'))
      rgb = im2uint16(rgb);
    end
  else
    % slope
    m = 1/(1/4);
    % precomputation
    mim = m*im;
    nmim = -mim;
    x1 = im<1/8;
    x2 = im>=1/8 & im<3/8;
    x3 = im>=3/8 & im<5/8;
    x4 = im>=5/8 & im<7/8;
    x5 = im>=7/8;
    rgb = cat( 3, ...
      x3.*(mim - (3/8)*m)   + x4.*1 + x5.*(nmim + (9/8)*m), ...
      x2.*(mim - (1/8)*m)   + x3.*1 + x4.*(nmim + (7/8)*m), ...
      x1.*(mim + 0.5)       + x2.*1 + x3.*(nmim + (5/8)*m));
  end
end


text:
function pass = test_curl( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e2*pref.techPrefs.chebfuneps;

% Example 1 : (0,x,xy)
Vx = ballfun(@(x,y,z)0);
Vy = ballfun(@(r,lam,th)r.*sin(th).*cos(lam), 'spherical');
Vz = ballfun(@(r,lam,th)r.*sin(th).*cos(lam).*r.*sin(th).*sin(lam), 'spherical');
V = ballfunv(Vx,Vy,Vz);
W = curl(V);
Exactx = ballfun(@(r,lam,th)r.*sin(th).*cos(lam), 'spherical');
Exacty = ballfun(@(r,lam,th)-r.*sin(th).*sin(lam), 'spherical');
Exactz = ballfun(@(r,lam,th)1, 'spherical');
Exact = ballfunv(Exactx,Exacty,Exactz);
pass(1) = norm(W-Exact)<tol;

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function result = blurOnGraph(signal,laplacian,time,steps)

h = time/steps;
nVertices = size(laplacian,1);

timeStep = @(x) (speye(nVertices) - h*laplacian) \ x;

result = signal;
for i=1:steps
    result = timeStep(result);

    % help fix numerical issues
    result(result<0) = 0;
    result = result / sum(result) * sum(signal);
end

author: gpeyre


text:
function imgOut = WardHistAdjTMO(img, nBin, LdMin, LdMax, bPlotHistogram, bDownsampling)
%
%        imgOut = WardHistAdjTMO(img, nBin, LdMin, LdMax, bPlotHistogram, bDownsampling)
%
%
%        Input:
%           -img: input HDR image
%           -nBin: number of bins for calculating the histogram (1,+Inf)
%           -LdMin: minimum luminance value of the dispLay
%           -LdMax: maximum luminance value of the dispLay
%           -bPlotHistogram:
%
%        Output:
%           -imgOut: tone mapped image
%
%     Copyright (C) 2010-21  Francesco Banterle
%
%     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/>.
%
%     The paper describing this technique is:
%     "A Visibility Matching Tone Reproduction Operator for High Dynamic Range Scenes"
% 	  by Gregory Ward Larson, Holly Rushmeier, Christine Piatko
%     in IEEE Transactions on Visualization and Computer Graphics 1997
%

%is it a gray/three color channels image?
check13Color(img);

checkNegative(img);

if(~exist('nBin', 'var'))
    nBin = 100;
end

if(nBin < 1)
    nBin = 100;
end

if(~exist('LdMin', 'var'))
    LdMin = 1; %cd/m^2
end

if(LdMin < 0.0)
    LdMin = 1;
end

if(~exist('LdMax', 'var'))
    LdMax = 100; %cd/m^2
end

if(LdMax <= 0.0)
    LdMax = 100;
end

if(LdMax < LdMin)
    tmp = LdMin;
    LdMin = LdMax;
    LdMax = tmp;
end

if(~exist('bPlotHistogram', 'var'))
    bPlotHistogram = 0;
end

if(~exist('bDownsampling', 'var'))
    bDownsampling = 0;
end

epsilon = 1e-6;

%compute luminance channel
L = lum(img);

%downsample according to fovea...
if(bDownsampling)
    L2 = WardDownsampling(L + epsilon);
else
    L2 = L + 1e-6;
end

%compute stastistics
LMin = min(L2(:));
LMax = max(L2(:));

Llog  = log(L2);

LlMin = log(LMin);
LlMax = log(LMax);

LldMin = log(LdMin + epsilon);
LldMax = log(LdMax + epsilon);

%compute the histogram H
H = zeros(nBin, 1);
delta = (LlMax - LlMin) / nBin;

for i=1:nBin
    indx = find(Llog > (delta * (i - 1) + LlMin) & Llog <= (delta * i + LlMin));
    H(i) = numel(indx);
end

%apply the histogram ceiling
maxH = max(H);
x_vis = LlMin:((LlMax - LlMin) / (nBin -1)):LlMax;

if(bPlotHistogram)
    bar(x_vis, H/maxH);
    hold on;
end

H = histogram_ceiling(H, delta / (LldMax - LldMin));

if(bPlotHistogram)
    bar(x_vis, H / maxH);
    hold off;
end

%compute P(x)
P = cumsum(H);
P = P / max(P);

%calculate tone mapped luminance
L(L > LMax) = LMax;
x = (LlMin:((LlMax - LlMin) / (nBin - 1)):LlMax)';
P_L = interp1(x , P , real(log(L)), 'linear');
Ld  = exp(LldMin + (LldMax - LldMin) * P_L);
%normalize in [0,1]
Ld  = (Ld - LdMin) / (LdMax - LdMin);

%change luminance
imgOut = ChangeLuminance(img, L, Ld);
end

author: banterle


text:
function [d, pre, post, height, cycle, pred] = dfs(adj_mat, start, directed)
% DFS Perform a depth-first search of the graph starting from 'start'.
% [d, pre, post, height, cycle, pred] = dfs(adj_mat, start, directed)
%
% d(i) is the time at which node i is first discovered.
% pre is a listing of the nodes in the order in which they are first encountered (opened).
% post is a listing of the nodes in the order in which they are last encountered (closed).
% A node is last encountered once we have explored all of its neighbors.
% If the graph is directed, i's neighbors are its children.
% If the graph is a tree, preorder is parents before children, and
% postorder is children before parents.
% For a DAG, topological order = reverse(postorder).
% height(i) is the height (distance) of node i from the start.
% 'cycle' is true iff a (directed) cycle is found.
% pred(i) is the parent of i in the dfs tree rooted at start.
% See Cormen, Leiserson and Rivest, "An intro. to algorithms" 1994, p478.

% We can detect undirected cycles by checking if we are about to visit a node n which we have
% already visited. To detect *directed* cycles, we need to know if n has been closed or is still open.
% For example (where arcs are directed down)
%   1    2
%   \   /
%     3
% Assume we visit 1, 3 and then 2 in order. The fact that a child of 2 (namely, 3) has
% already been visited is okay, because 3 has been closed.
% The algorithms in Aho, Hopcroft and Ullman, and Sedgewick, do not detect directed cycles.

n = length(adj_mat);

global white gray black
white = 0; gray = 1; black = 2;

color = white*ones(1,n);
d = zeros(1,n);
height = zeros(1,n);
pred = zeros(1,n);
pre = [];
post = [];
cycle = 0;
global count
count = 0;
h = 0;
[d, pre, post, height, cycle, color, pred] = ...
    dfs2(adj_mat, start, directed, h, d, pre, post, height, cycle, color, pred);

%%%%%%%%%%

function [d, pre, post, height, cycle, color, pred] = ...
    dfs2(adj_mat, i, directed, h, d, pre, post, height, cycle, color, pred)

global count
global white gray black

color(i) = gray;
count = count + 1;
d(i) = count;
pre = [pre i];
height(i) = h;
if directed
  ns = children(adj_mat, i);
else
  ns = neighbors(adj_mat, i);
end
for j=1:length(ns)
  n=ns(j);
  if ~directed & n==pred(i) % don't go back up the edge you just came down
    % continue
  else
    if color(n) == gray % going back to a non-closed vertex via a new edge
      %fprintf(1, 'cycle from %d to %d\n', i, n);
      cycle = 1;
    end
    if color(n) == white % not visited n before
      pred(n)=i;
      [d, pre, post, height, cycle, color, pred] = ...
	  dfs2(adj_mat, n, directed, h+1, d, pre, post, height, cycle, color, pred);
    end
  end
end
color(i) = black;
post = [post i];

author: bayesnet

repo: bnt


text:
% NMF-PG: NMF solved by Projected Gradient
% process_video('NMF', 'NMF-PG', 'dataset/demo.avi', 'output/demo_NMF-PG.avi');
alg_path_aux = fullfile(lrs_conf.nmf_path,'NMF-DTU-Toolbox');
addpath(genpath(alg_path_aux));

M = sparse(M);

% cjlin: Alternative non-negative least squares using projected gradients.
[W, H] = nmf(M,2,'cjlin');

L = W * H;
S = M - L;

rmpath(genpath(alg_path_aux));

author: andrewssobral

repo: lrslibrary


text:
function Rlabel=RegionMerging(varargin)
% post-processing over segment Image based RAG region merging
% input:
% I by M x N gray Image
% label: initial segment result
% Treshold: interation temination condion default is 1000
% RegionNum: the Region Number after Mereging defalut is 2(forground and
% background
% output:
% Rlable: merging result
% reference:
% writed by radishgiant

[I,Rlabel,Treshold,RegionNumber] = parse_inputs(varargin{:});
minsim=Treshold-1;
K=length(unique(Rlabel));
[matrixSize(1),matrixSize(2),channel]=size(I);
X=reshape(I,matrixSize(1)*matrixSize(2),channel);
Vset=cell(K,1);
Pset=cell(K,1);
delnode=[];
labelstart=min(Rlabel(:));
if labelstart==0
    Rlabel=Rlabel+1;
elseif labelstart>1
    Rlabel=Rlabel-labelstart+1;
end
for i=1:K
    [r,c]=find(Rlabel==i);
    Vset{i}=[r,c];
    Pset{i}=X(sub2ind(matrixSize,r,c),:);
end

while length(unique(Rlabel))>RegionNumber || minsim(1)<=Treshold
Adj=imRAG(Rlabel);
sim=partionsimliar(Pset,Adj,delnode);
[minsim,simInd]=sort(sim,1,'ascend');
node1=Adj(simInd(1),1);
node2=Adj(simInd(1),2);
node=min(node1,node2);
nodem=max(node1,node2);
Rlabel(Rlabel==node1|Rlabel==node2)=node;
Pset{node}=[Pset{node1};Pset{node2}];
Pset{nodem}=[];
Vset{node}=[Vset{node1};Vset{node2}];
Vset{nodem}=[];
delnode=[delnode;nodem];
end

end

function [I,label,Treshold,RegionNumber] = parse_inputs(varargin)

narginchk(2,4);
I = varargin{1};
label=varargin{2};
validateattributes(I,{'numeric'}, {'2d' '3d'}, ...
    mfilename, 'I', 1);
validateattributes(label,{'numeric'}, {'2d'}, ...
    mfilename, 'label', 2);
if nargin<3
    Treshold=1000;
    RegionNumber=2;
elseif nargin<4
    Treshold=varargin{3};
    RegionNumber=2;
else
    Treshold=varargin{3};
    RegionNumber=varargin{4};
validateattributes(Treshold,{'numeric'}, {'scalar','real'}, ...
    mfilename, 'Treshold', 3);
validateattributes(RegionNumber,{'numeric'}, {'scalar','real','>=',2}, ...
    mfilename, 'RegionNumber', 4);
end

end

author: radishgiant

repo: ThresholdAndSegment