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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deftypemethod  @@sym {} min (@var{a})
%% @deftypemethodx @@sym {} min (@var{a}, @var{b})
%% @deftypemethodx @@sym {} min (@var{a}, [], @var{dim})
%% @deftypemethodx @@sym {[@var{r}, @var{I}] =} min (@dots{})
%% Return minimum value of a symbolic vector or vectors.
%%
%% Example:
%% @example
%% @group
%% min(sym(1), sym(2))
%%   @result{} (sym) 1
%% m = min([1 sym(pi)/4 6])
%%   @result{} m = (sym)
%%       π
%%       ─
%%       4
%% [m, I] = min([sym(1) 0 6])
%%   @result{} m = (sym) 0
%%   @result{} I = 2
%% @end group
%% @end example
%%
%% @seealso{@@sym/max}
%% @end deftypemethod

function [z, I] = min(A, B, dim)

  if (nargout <= 1)
    if (nargin == 1)
      if (isvector(A))
        z = pycall_sympy__ ('return Min(*_ins[0])', A);
      else
        z = min(A, [], 1);
      end
    elseif (nargin == 2)
      z = elementwise_op('Min', sym(A), sym(B));
    elseif (nargin == 3)
      assert (isempty (B))
      assert (logical(dim == 1) || logical(dim == 2))

      cmd = { '(A, dim) = _ins'
              'if not A.is_Matrix:'
              '    A = sp.Matrix([A])'
              'if dim == 0:'
              '    if A.rows == 0:'
              '        return A'
              '    return Matrix([[Min(*A.col(i)) for i in range(0, A.cols)]])'
              'elif dim == 1:'
              '    if A.cols == 0:'
              '        return A'
              '    return Matrix([Min(*A.row(i)) for i in range(0, A.rows)])' };
      z = pycall_sympy__ (cmd, A, dim - 1);
    else
      print_usage ();
    end
    return
  end

  % We have second output: need the index of the minimum, can't use "Min"
  if (nargin == 1)
    if (isvector(A))
      cmd = { 'A = _ins[0]'
              'val, idx = min((val, idx) for (idx, val) in enumerate(A))'
              'return val, idx+1' };
      [z, I] = pycall_sympy__ (cmd, A);
    else
      [z, I] = min(A, [], 1);
    end
  elseif (nargin == 3)
    assert (isempty (B))
    assert (logical(dim == 1) || logical(dim == 2))

    cmd = { '(A, dim) = _ins'
            'def myargmin(z):'
            '    return min((val, idx) for (idx, val) in enumerate(z))'
            'if not A.is_Matrix:'
            '    A = sp.Matrix([A])'
            'if dim == 0:'
            '    if (A.cols == 0):'
            '        return (Matrix(0, 0, []), Matrix(0, 0, []))'
            '    if (A.rows == 0):'
            '        return (A, A)'
            '    val_idx_pairs = [myargmin(A.col(i)) for i in range(0, A.cols)]'
            '    m, I = zip(*val_idx_pairs)'
            '    return (Matrix([m]), Matrix([I]))'
            'elif dim == 1:'
            '    if (A.rows == 0):'
            '        return (Matrix(0,0,[]), Matrix(0,0,[]))'
            '    if (A.cols == 0):'
            '        return (A, A)'
            '    val_idx_pairs = [myargmin(A.row(i)) for i in range(0, A.rows)]'
            '    m, I = zip(*val_idx_pairs)'
            '    return (Matrix(m), Matrix(I))' };

    [z, I] = pycall_sympy__ (cmd, A, dim - 1);
    I = double(I);
    if (~isempty(I))
      I = I + 1;
    end
  else
    print_usage ();
  end

end

%!test
%! % scalars with dim
%! a = min(sym(pi), [], 1);
%! b = sym(pi);
%! assert (isequal (a, b));
%! a = min(sym(pi), [], 2);
%! assert (isequal (a, b));
%! a = max(sym(pi), [], 1);
%! assert (isequal (a, b));
%! a = max(sym(pi), [], 2);
%! assert (isequal (a, b));

%!shared A, D
%! D = [0 1 2 3];
%! A = sym(D);

%!test
%! % row vectors
%! assert (isequal (min(A), sym(min(D))))
%! assert (isequal (min(A), sym(0)))
%! assert (isequal (max(A), sym(max(D))))
%! assert (isequal (max(A), sym(3)))

%!test
%! % row vectors with dim
%! assert (isequal (min(A, [], 1), sym(min(D, [], 1))))
%! assert (isequal (min(A, [], 2), sym(min(D, [], 2))))
%! assert (isequal (max(A, [], 1), sym(max(D, [], 1))))
%! assert (isequal (max(A, [], 2), sym(max(D, [], 2))))

%!test
%! % column vectors
%! A = A.';
%! D = D.';
%! assert (isequal (min(A), sym(min(D))))
%! assert (isequal (min(A), sym(0)))
%! assert (isequal (max(A), sym(max(D))))
%! assert (isequal (max(A), sym(3)))

%!test
%! % row vectors with dim
%! assert (isequal (min(A, [], 1), sym(min(D, [], 1))))
%! assert (isequal (min(A, [], 2), sym(min(D, [], 2))))
%! assert (isequal (max(A, [], 1), sym(max(D, [], 1))))
%! assert (isequal (max(A, [], 2), sym(max(D, [], 2))))

%!shared

%!test
%! % empty
%! a = min(sym([]));
%! assert(isempty(a))
%! a = max(sym([]));
%! assert(isempty(a))

%!test
%! % matrix
%! A = [1 4 6; 2 2 5];
%! A = sym(A);
%! assert (isequal (min(A), sym([1 2 5])))
%! assert (isequal (min(A, [], 1), sym([1 2 5])))
%! assert (isequal (min(A, [], 2), sym([1; 2])))
%! assert (isequal (max(A), sym([2 4 6])))
%! assert (isequal (max(A, [], 1), sym([2 4 6])))
%! assert (isequal (max(A, [], 2), sym([6; 5])))

%!test
%! % index output is double not sym
%! [m, I] = min(sym(2), [], 1);
%! assert (strcmp(class(I), 'double'))
%! [m, I] = max(sym(2), [], 1);
%! assert (strcmp(class(I), 'double'))

%!test
%! % empty rows/columns, I is double
%! A = sym(zeros(0, 4));
%! [m, I] =  min(A, [], 1);
%! assert (strcmp(class(I), 'double'))
%! [m, I] =  max(A, [], 1);
%! assert (strcmp(class(I), 'double'))
%! A = sym(zeros(3, 0));
%! [m, I] =  min(A, [], 2);
%! assert (strcmp(class(I), 'double'))
%! [m, I] =  max(A, [], 2);
%! assert (strcmp(class(I), 'double'))

%!test
%! % index output
%! A = [0 1 9; 10 7 4];
%! B = sym(A);
%! [m1, I1] = min(A);
%! [m2, I2] = min(B);
%! assert (isequal (I1, I2))
%! assert (isequal (m1, double(m2)))
%! [m1, I1] = max(A);
%! [m2, I2] = max(B);
%! assert (isequal (I1, I2))
%! assert (isequal (m1, double(m2)))

%!test
%! % index output, with dim
%! A = [0 1 9; 10 7 4];
%! B = sym(A);
%! [m1, I1] = min(A, [], 1);
%! [m2, I2] = min(B, [], 1);
%! assert (isequal (I1, I2))
%! assert (isequal (m1, double(m2)))
%! [m1, I1] = min(A, [], 2);
%! [m2, I2] = min(B, [], 2);
%! assert (isequal (I1, I2))
%! assert (isequal (m1, double(m2)))
%! [m1, I1] = max(A, [], 1);
%! [m2, I2] = max(B, [], 1);
%! assert (isequal (I1, I2))
%! assert (isequal (m1, double(m2)))
%! [m1, I1] = max(A, [], 2);
%! [m2, I2] = max(B, [], 2);
%! assert (isequal (I1, I2))
%! assert (isequal (m1, double(m2)))

%!test
%! % empty columns
%! A = sym(zeros(0, 4));
%! [m, I] =  min(A, [], 1);
%! assert (isequal (size(m), [0 4]))
%! assert (isequal (size(I), [0 4]))
%! [m, I] =  max(A, [], 1);
%! assert (isequal (size(m), [0 4]))
%! assert (isequal (size(I), [0 4]))

%!test
%! % empty rows
%! A = sym(zeros(3, 0));
%! [m, I] =  min(A, [], 2);
%! assert (isequal (size(m), [3 0]))
%! assert (isequal (size(I), [3 0]))
%! [m, I] =  max(A, [], 2);
%! assert (isequal (size(m), [3 0]))
%! assert (isequal (size(I), [3 0]))

%!test
%! % another empty case
%! % we differ slightly from double which gives 1x0/0x1
%! A = sym(zeros(3, 0));
%! [m, I] =  min(A, [], 1);
%! assert (isempty (m))
%! assert (isempty (I))
%! A = sym(zeros(0, 3));
%! [m, I] =  min(A, [], 2);
%! assert (isempty (m))
%! assert (isempty (I))

%!test
%! % empty without index output
%! A = sym(zeros(3, 0));
%! assert (isempty (min (A, [], 1)))
%! assert (isempty (max (A, [], 1)))
%! assert (isempty (min (A, [], 2)))
%! assert (isempty (max (A, [], 2)))
%! A = sym(zeros(0, 3));
%! assert (isempty (min (A, [], 1)))
%! assert (isempty (max (A, [], 1)))
%! assert (isempty (min (A, [], 2)))
%! assert (isempty (max (A, [], 2)))

%!test
%! % binary op form, one a scalar
%! A = sym([3 1 9]);
%! m = min(A, sym(2));
%! M = max(A, sym(2));
%! assert (isequal (m, sym([2 1 2])))
%! assert (isequal (M, sym([3 2 9])))
%! m = min(sym(2), A);
%! M = max(sym(2), A);
%! assert (isequal (m, sym([2 1 2])))
%! assert (isequal (M, sym([3 2 9])))

%!test
%! % binary op form, both scalar
%! m = min(sym(1), sym(2));
%! M = max(sym(2), sym(2));
%! assert (isequal (m, sym(1)))
%! assert (isequal (M, sym(2)))

%!test
%! syms x y
%! assert (isequal (children (min (x, y)), [x y]))

%!test
%! syms x y z
%! A = [x 1; y z];
%! assert (isequal (min (A, [], 1), [min(x, y)  min(1, z)]))
%! assert (isequal (max (A, [], 1), [max(x, y)  max(1, z)]))
%! assert (isequal (min (A, [], 2), [min(x, 1); min(y, z)]))
%! assert (isequal (max (A, [], 2), [max(x, 1); max(y, z)]))

%!test
%! syms x y positive
%! a = min([x 2 y -6]);
%! assert (isequal (a, -6))
%! a = max([x y -6]);
%! assert (isequal (a, max(x, y)))

%!test
%! syms x negative
%! a = min([x 6 10]);
%! assert (isequal (a, x))
%! a = max([x -2 6]);
%! assert (isequal (a, 6))


text:
%
%
% Copyright (c) 2012 University of Crete - Computer Science Department (UOC-CSD)
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file 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 Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function p = wrappednormpdf(x, m, s, N)

    % N=10 gives a log error smaller than 1e-12 for s=8
    % with s=8 distribution is "almost" uniform (std < 1e-12)
    if nargin<4; N=10; end

    p = 0;
    for k=-N:N
        p = p + normpdf(x+2.*pi.*k, m, s);
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [c,ia,ib] = intersect_sorted(a,b)
%INTERSECT_SORTED     Set intersection between sorted sets.
% INTERSECT_SORTED(A,B) when A and B are vectors returns the values common
% to both A and B.  A and B must be sorted and unique, and the result will be
% sorted and unique.
% [C,IA,IB]=INTERSECT_SORTED(A,B) also returns indices such that C=A(IA) and C=B(IB).

% Written by Tom Minka
% (c) Microsoft Corporation. All rights reserved.

if nargout>1
    [tf,ib]=ismember_sorted(a,b);
    ib=ib(tf);
    c = a(tf);
    [tf,ia]=ismember_sorted(b,a);
    ia=ia(tf);
else
    c = a(ismember_sorted(a,b));
end

author: tminka

repo: lightspeed

save_path: github-repos/MATLAB/tminka-lightspeed


text:
function H =  sndspd_mean(z,ssp,d)
% harmonic_mean
%
% Val Schmidt
% Center for Coastal and Ocean Mapping
% University of New Hampshire
% 2012
%
% A routine to calculate the harmonic mean sound speed.
%
%   z:         Vector of depth values. (m)
%   ssp:       Vector of sound speed values at the associated depths. (m/s)
%   d:         [Optional] maximum depth (m) (see below)
%
% Suppose one is given a sound speed profile (sound speeds [ssp] and depths [z]).
% One wishes to know how long it will take a sound pulse (say from a sonar)
% to travel down through the water column.
%
% Most people will realize that one cannot simply average the sound speeds
% and divide by the total depth. However many will incorrectly assume that
% a depth-weighted mean sound will produce the correct answer.
% Unfortunately it does not.
%
% One must instead divide the total depth by the sum of the time it takes
% to the sound to pass through each layer of constant sound speed. This is
% called the 'harmonic mean'. The harbmonic mean is calculated like so.
%
%              |    i=n                | -1
%              |  SUM    delta z(i)    |
%              |   i=1   ---------     |
%              |            c(i)       |
% C_harmonic = |   ------------------- |
%              |    i=n                |
%              |  SUM     delta z(i)   |
%              |     i=1               |
%              |                       |
%
% When the maximum depth, 'd', is specified, two possible things may result
% depending on whether d is greater or less than the maximum depth in z.
%
% Sometimes one has a complete sound speed profile, but only needs the
% harmonic sound speed to some max depth less than the maximum depth of the
% profile. One may then specify the maximum depth as the argument 'd', and
% the function will return a harmonic mean to only that depth.
%
% Alternatively, one frequently has a sound speed profile for the upper
% portion of the water column only, but desireds a harmonic sound speed all
% the way to the bottom. Below some max depth temperature and salnity are
% constant and sound speed increases nearly linearly (to about 1 part in
% 10^4) with increasing depth. In this case, one may specify a the bottom
% depth, d. The function will then approximate the the sound speed at the
% bottom linearly from the final value in the ssp vector and calculate a
% harmonic sound speed for the entire path.
%
% When the sound speed profile is a constant linear gradient rather than
% discrete values one can calculate the time required to traverse the
% layer directly rather than approximating the gradient in discrete steps.
% That time is given by
%
%       1      c(i+1)
% t =   - ln ( ------ )
%       g       c(i)
%
% This value is then added to the summation of travel times (with a
% corresponding addition to the summation of depth differences) in the
% haronic mean calculation above.

%%

% Initialize a variable. variables
t=0;

% Mackenzie Formula first-order term for sound speed as a function of
% depth.
D1=.0163;

% Check some details.
if ( ~isequal( size(z),size(ssp) ) )
    error('ERROR input vectors are not of same size');
end

% The algorithm below will produce incorrect results if the first depth
% value is not zero. So we force it to be zero, shifting the other depths
% appropriately.
top=z(1);
z=z-top;

% One can optionally specify a max bottom depth. This value may be less
% than, or deeper than, the lowest depth in the z vector. If it is less than
% than z(end), the harmonic mean to the depth, d is calculated. If it is more
% than z(end), the sound speed profile is extended to the depth d using the
% first-order depth dependent term from the Mackenzie equation and the
% results of the Harmonic mean calculation for a gradient.

if ( nargin==3 )

    % First we shift the d value as we did the z values to make sure
    % everything matches.
    d=d-top;

    % If d is less than the max depth in z truncate the z and ssp vectors
    % to that depth.
    if (d < max(z) )
        idx=find( z <= d);
        z=z(idx);
        ssp=ssp(idx);

    % Otherwise we calculate the extra time in this final layer.
    else

    % We need the sound speed at hte specified bottom.
    ss_bottom = ssp(end) + ( d - z(end) ) * D1;

    % This is the equation for the time spent in a layer having a gradient
    % sound speed, D1, and where the difference in sound speed from the top to
    % the bottom of the layer is delta_ss.
    t= 1/D1 *log( ss_bottom / ssp(end) );

    end

else
    % Give d the max depth regardless.
    d=z(end);
end

% Calculate the harmonic mean sound speed.
H = ( (sum(diff(z)./ssp(1:(length(ssp)-1))) + t ) / d )^(-1);

author: Sable

repo: mcbench-benchmarks


text:
function call_price=american_call_perpetual(S, K, r, q, sigma)

%--------------------------------------------------------------------------
%
% DESCRIPTION:
%
% Price for an american perpetual call option
%
%
% Reference:
%
% John Hull, "Options, Futures and other Derivative Securities",
% Prentice-Hall, second edition, 1993.
%
%--------------------------------------------------------------------------
%
% INPUTS:
%
%  S:       spot price
%  K:       exercice price
%  r:       interest rate
%  q:       dividend yield
%  sigma:   volatility
%
%--------------------------------------------------------------------------
%
% OUTPUT:
%
% call_price: price of a call option
%
%--------------------------------------------------------------------------
%
% Author:  Paolo Z., February 2012
%
%--------------------------------------------------------------------------

sigma_sqr=(sigma^2);
h1 = 0.5 - ((r-q)/sigma_sqr);
h1 = h1 + sqrt( (((r-q)/sigma_sqr-0.5)^2)+2.0*r/sigma_sqr );

call_price=(K/(h1-1.0))*( (((h1-1.0)/h1)*(S/K))^h1 );


text:
function vgc_plot(k,kstar,v,Polc,g, METHOD)
% plot the value function against k
% plot capital policy for next period against captial this period
% save variable
% input is the k vector, optimal-- kstar, value --v, capital policy - g
% METHOD=1 if deterministic noninterpolation
% Method=2----Deterministic, continuous state space model, With interpolation
gstar=ppval(spline(k,g),kstar);  % g value at optimal kstar
vstar=ppval(spline(k,v),kstar);  % v value at optimal kstar

T=char( ' No Interpolation, Deterministic',...
        ' With Interpolation, Deterministic');

figure(1);
plot(k,v); hold on;
plot(kstar,vstar,'ro');
text(kstar+.1,vstar-.1,['k*=' num2str(kstar)]);
title(['Value Function--' T(METHOD,:)]);
xlabel('capital k');  ylabel('value');
axis square; hold off;
print('-dpsc',['ValFun' num2str(METHOD)]);

figure(2);
subplot(2,1,1); plot(k,Polc);
title('Policy Function-Consumption');
xlabel('Capital Kt');   ylabel('Consumption');

subplot(2,1,2); plot(k,g,k,k,'-'); hold on;
plot(kstar,gstar,'ro');
text(kstar+.1,gstar-.1,['g*=' num2str(gstar) '| k*=' num2str(kstar)]);
title('Policy Function-Capital');
xlabel('Capital Kt');   ylabel('Capital Next Period');
hold off;

print('-dpsc', ['PolFun' num2str(METHOD)]);
save(['growth' num2str(METHOD)]);


text:
function[p] = elec_load_brainstorm(p)

% elec_load_brainstorm - Load brainstorm Channel data into eeg_toolbox
%
% Usage: [p] = elec_load_brainstorm(p)
%
% where:    p.elec.path & p.elec.file define the brainsform file
%
% The brainstorm file consists of Channel, which is an array of
% structures.  The fields are:
%
%     Loc     - a 3x2 matrix of electrode coordinates (x,y,z in rows).
%               BrainStorm (x,y,z meters) = 3Dspace (x,y,z cm) / 100.
%     Orient  - a corresponding matrix of sensor orientations (MEG);
%               all zero for EEG.
%     Weight  - a vector of relative or absolute weights (eg, amplification);
%               all ones for this routine.
%     Type    - a character string, 'EEG' in this case.
%     Name    - a charater string indicating the electrode name.
%     Comment - a charater string indicating the reference electrode.  Empty
%               for active electrodes and 'EEG REF' for reference electrode.
%
% See brainstorm website at http://neuroimage.usc.edu/, including a
% download pdf file describing the brainstorm database formats.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:55 $

% Licence:  GNU GPL, no express or implied warranties
% History:  11/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('p', 'var'), error('...no p struct input\n');
elseif isempty(p),     error('...no p struct input\n');
end

fprintf('\nELEC_LOAD_BRAINSTORM...\n'); tic;

fprintf('...Converting brainstorm to p structure.\n');

[path,file,ext] = fileparts(strcat(p.elec.path,filesep,p.elec.file));
bsfile = fullfile(path,[file ext]);

load(bsfile, 'Channel');

for i=1:length(Channel),

    if findstr(lower(Channel(i).Comment),'ref'),
        p.elec.data.ref = Channel(i).Loc(:,1)';
    else
        p.elec.data.x(i,1) = Channel(i).Loc(1,1);
        p.elec.data.y(i,1) = Channel(i).Loc(2,1);
        p.elec.data.z(i,1) = Channel(i).Loc(3,1);
        p.elec.data.label{i,1} = Channel(i).Name;
    end
end

p.elec.data.centroid = [0 0 0];
p.elec.data.nasion = [];
p.elec.data.lpa = [];
p.elec.data.rpa = [];

t = toc; fprintf('...done (%6.2f sec).\n\n',t);

return

author: PatternRecognition

repo: OpenBMI


text:
function varsigma = fgplvmPosteriorVar(model, X);

% FGPLVMPOSTERIORVAR Variances of the posterior at points given by X.
% FORMAT
% DESC returns the posterior mean and variance for a given set of
% points.
% ARG model : the model for which the posterior will be computed.
% ARG x : the input positions for which the posterior will be
% computed.
% RETURN sigma : the variances of the posterior distributions.
%
% SEEALSO : gpPosteriorVar, fgplvmCreate
%
% COPYRIGHT : Neil D. Lawrence, 2009

% FGPLVM

varsigma = gpPosteriorVar(model, X);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
I = imread('bh.png');
%I = imread('C:\Users\Mostwanted\Desktop\Wuhan_China.jpg');
subplot(1,2,1);
imshow(I);

I=double(I);
f=I(:,:,1);
ff=I(:,:,2);
fff=I(:,:,3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k1=4;
k2=5;
r=161;
alf=1458;
nn=floor((r+1)/2);
for i=1:r
    for j=1:r
        b(i,j) =exp(-((i-nn)^2+(j-nn)^2)/(k1*alf))/(k2*pi*alf*10000); % Gaussian 1
   end
end

k1=8;
k2=8;
r=161;
alf=1458;
nn=floor((r+1)/2);
for i=1:r
    for j=1:r
        bb(i,j) =exp(-((i-nn)^2+(j-nn)^2)/(k1*alf))/(k2*pi*alf*10000);     % Gaussian 2
   end
end

k1=0.5;
k2=0.5;
r=161;
alf=1458;
nn=floor((r+1)/2);
for i=1:r
    for j=1:r
        bbb(i,j) =exp(-((i-nn)^2+(j-nn)^2)/(k1*alf))/(k2*pi*alf*10000);  % Gaussian 2 3
    end
    end
%%%%%%%%%%% R component of treatment %%%%%%%%%%%%%
Img = double(f);
[m,n]=size(f);

aa=125;

for i=1:m
    for j=1:n
        C(i,j)=log(1+aa*(Img(i,j)/I(i,j)));
    end
end

K=imfilter(Img,b);
KK=imfilter(Img,bb);
KKK=imfilter(Img,bbb);

for i=1:m
    for j=1:n
       G(i,j)=1/3*(log(Img(i,j)+1)-log(K(i,j)+1));
        G(i,j)=1/3*(log(Img(i,j)+1)-log(KK(i,j)+1))+G(i,j);
         G(i,j)=C(i,j)*(1/3*(log(Img(i,j)+1)-log(KKK(i,j)+1))+G(i,j));
    end
end

mi=min(min(G));
ma=max(max(G));
       L=(G-mi)*255/(ma-mi);
%%%%%%%%%%%%%% G Processing Components %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Img = double(ff);
[m,n]=size(ff);

aa=125;
for i=1:m
    for j=1:n
        CC(i,j)=log(1+aa*(Img(i,j)/I(i,j)));
    end
end

K=imfilter(Img,b);
KK=imfilter(Img,bb);
KKK=imfilter(Img,bbb);
for i=1:m
    for j=1:n
       G(i,j)=1/3*(log(Img(i,j)+1)-log(K(i,j)+1));
        G(i,j)=1/3*(log(Img(i,j)+1)-log(KK(i,j)+1))+G(i,j);
         G(i,j)=CC(i,j)*(1/3*(log(Img(i,j)+1)-log(KKK(i,j)+1))+G(i,j));
    end
end

mi=min(min(G));
ma=max(max(G));
       LL=(G-mi)*255/(ma-mi);
%%%%%%%%%%%%% With the B component of the Department %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Img = double(fff);
[m,n]=size(fff);

aa=125;
for i=1:m
    for j=1:n
        CCC(i,j)=log(1+aa*(Img(i,j)/I(i,j)));
    end
end

K=imfilter(Img,b);
KK=imfilter(Img,bb);
KKK=imfilter(Img,bbb);

for i=1:m
    for j=1:n
       G(i,j)=1/3*(log(Img(i,j)+1)-log(K(i,j)+1));
        G(i,j)=1/3*(log(Img(i,j)+1)-log(KK(i,j)+1))+G(i,j);
         G(i,j)=CCC(i,j)*(1/3*(log(Img(i,j)+1)-log(KKK(i,j)+1))+G(i,j));
    end
end

mi=min(min(G));
ma=max(max(G));

       LLL=(G-mi)*255/(ma-mi);
%%%%%%%%%%%% Department of Integrated color image of science %%%%%%%%%%%%%%%
msrcr=cat(3,L,LL,LLL);
subplot(1,2,2);
imshow(uint8(msrcr));
%imwrite(uint8(msrcr),'C:\Users\Mostwanted\Desktop\Wuhan_China_outcr1.jpg');
%imwrite(uint8(msrcr),'C:\Users\Mostwanted\Desktop\Washington_DC_outcr1.jpg');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function varargout=lpcbwexp(varargin)
%
% For calling details please see v_lpcbwexp.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpcbwexp.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpcbwexp.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpcbwexp(varargin{:});
else
    v_lpcbwexp(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%	SurfBox-MATLAB (c)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%	Yue Lu and Minh N. Do
%%
%%	Department of Electrical and Computer Engineering
%%	Coordinated Science Laboratory
%%	University of Illinois at Urbana-Champaign
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%	IRCrecF.m
%%
%%	First created: 08-14-05
%%	Last modified: 04-13-06
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function subband =IRCrecF(subband, K, level, Hflts)

%   Iteratively Resampled Checkerboard Filter Bank Reconstruction (Fourier
%   Domain Implementation)
%
%   *****  INPUT:  *****
%
%   subband: N-dimensional (N >= 2) input signal in the FREQUENCY domain.
%
%   K: channel index, 1 <= K <= N
%
%   level: a row vector specifying the decomposition levels for each pair
%   of dimensions. Note: only level(m) for m ~= k will be used.
%
%   Hflts: A cell array containing two decomposition checkerboard filters
%   {H0, c0, H1, c1}
%
%   *****  OUTPUT:  *****
%
%   subband: N-D output signal in the FREQUENCY domain.
%
%   See also IRCdecF.m
%

%% dimension of the problem
N = ndims(subband);

%% check parameter validity
if (length(level) ~= N) | (K > length(level)) ...
        | (sum(level < 0) ~= 1) | (level(K) >= 0)
    error('Input parameter LEVEL is not valid.');
end

%% prepration
szS = size(subband);
szH0 = size(Hflts{1}); szH1 = size(Hflts{3});
szFlt = max(szH0, szH1);

%% Subscript cell array
subary = repmat({':'}, [2 N]);

%% preallocate the memory for a temp variable
subtmp = repmat(complex(0), szS);

%% start working
for m = N : -1 : 1  %% N is a small number like 2,3, or 4.

    %% We work on pairs of dimensions: (K,1), (K,2), ... (K, N)
    if m == K   %% except on dimensions (K, K)
        continue;
    end

    %% Size of the 2-D slice
    ftSpat = zeros((szS(K)-1) * 2, szS(m) / 2^(level(m)-1)); %% get the full size

    %% This one is ``thinner'' in the m dimension. This will reduce
    %% the number of filter resampling (shearing) operations.
    ftSpatSml = zeros((szS(K)-1) * 2, szFlt(2));

    indm = [1: szS(m)];

    %% This variable is what we actually use to multiply with the subband
    ftSlice = repmat(complex(0), [szS(K), szS(m)]);

    %% IRCdec on dimensions (K, m)
    for le = level(m) : -1 : 1 %% Another small number, like 0, 1, 2, 3

        if any(size(ftSpat) < szFlt)
            error('Image size is smaller than the filter size. Try using more compact filters.');
        end

        %% select the right column indices
        selcol = mod(floor( (indm - 1) / (szS(m) / 2^le) ), 2) == 0;
        subary{1,m} = indm(selcol);
        subary{2,m} = indm(~selcol);

        for chan = 0 : 1    %% Two channels

            %% Upsampling by 2 along dimension m
            subdouble = subband;
            subdouble(subary{2-chan,:}) = subband(subary{1+chan,:});

            H = Hflts{2*chan + 1};  %% Get the decomposition filter
            ctr = Hflts{2*chan + 2};
            szH = size(H);

            %% To get back to the original K-m dimension, since ftSlice
            %% might have been transposed.
            ftSlice = reshape(ftSlice, [szS(K), szS(m)]);

            if le == 1
                %% Get the FFT of filters
                ftSpat(:) = 0;
                ftSpat(1:szH(1), 1:szH(2)) = H;
                ftFreq = fft2(circshift(ftSpat, 1 - ctr));
                 %% Remove the complex conjugate symmetric part
                ftSlice = ftFreq(1:szS(K), :);
            else
                %% We will resample the checkerboard filters before filtering.
                %% This has the advantage of carrying out the resampling
                %% operation on 2-D filters instead of on N-D signals. Also by
                %% doing this, we avoid having to backsample at the end.
                %% Resampling in the frequency domain is generally MUCH
                %% harder than what one might think at first glance,
                %% especially when the data is of rectangular sizes, i.e.,
                %% when size(x, 1) ~= size(x, 2).

                nsubs = 2^(le - 1); %% number of subbands
                tmp_sub = [1 : szS(m) / nsubs]; %% a set of subscripts

                for n = 0 : nsubs - 1
                    %% Get the shearing factor (SF)
                    sf = nsubs - 1 - 2 * n;

                    %% After shearing, the filter H becomes ``taller''. We need
                    %% to make sure that we can still handle this filter
                    %% size.
                    if size(ftSpatSml, 1) < (szH(1) + abs(sf)* (szH(2)-1))
                        error('Image size is smaller than the filter size (from shearing). Try using more compact filters.');
                    end

                    %% load the filter
                    ftSpatSml(:) = 0;
                    ftSpat(:) = 0;
                    if sf > 0
                        ftSpatSml(1:szH(1), 1: szH(2)) = H;
                        newctr = [ctr(1) + sf * (ctr(2)-1), ctr(2)];
                    else
                        ftSpatSml(end-szH(1)+1:end, 1: szH(2)) = H;
                        newctr = [size(ftSpatSml, 1) + ctr(1) - szH(1) + sf * (ctr(2)-1), ctr(2)];
                    end

                    %% resample the filter in the spatial domain
                    ftSpatSml = resampc(ftSpatSml, (3 + sign(sf))/2, abs(sf), 'per');

                    ftSpat(:, 1:size(ftSpatSml, 2)) = ftSpatSml;
                    ftFreq = fft2(circshift(ftSpat, 1 - newctr));
                    ftSlice(:, tmp_sub) = ftFreq(1:szS(K), :);

                    tmp_sub = tmp_sub + szS(m) / nsubs;
                end

            end

            if K > m
                ftSlice = ftSlice .';
            end

            %% Pointwise multiplication of Filters and Input Data
            sz = ones(1, N);
            sz([K,m]) = szS([K,m]);

            ftSlice = reshape(ftSlice, sz);
            if chan == 0
                subtmp = subdouble .* repmat(ftSlice, szS ./ sz);
            else
                subband = subtmp + subdouble .* repmat(ftSlice, szS ./ sz);
            end
            clear subdouble;

        end  %% chan = 0:1

       %% Double the size of the 2-D slice along dimension m
       ftSpat = zeros(size(ftSpat) .* [1 2]);

    end %% le = 1 : level(m)

    %% We will work on the next dimension.
    subary{1, m} = ':'; subary{2,m} = ':';
end

%%	This software is provided "as-is", without any express or implied
%%	warranty. In no event will the authors be held liable for any
%%	damages arising from the use of this software.

author: thomaskuestner


text:
function [mapUpdate, robPoseMapFrame, laserEndPntsMapFrame] = inv_sensor_model(map, scan, robPose, gridSize, offset, probOcc, probFree)
% Compute the log odds values that should be added to the map based on the inverse sensor model
% of a laser range finder.

% map is the matrix containing the occupancy values (IN LOG ODDS) of each cell in the map.
% scan is a laser scan made at this time step. Contains the range readings of each laser beam.
% robPose is the robot pose in the world coordinates frame.
% gridSize is the size of each grid in meters.
% offset = [offsetX; offsetY] is the offset that needs to be subtracted from a point
% when converting to map coordinates.
% probOcc is the probability that a cell is occupied by an obstacle given that a
% laser beam endpoint hit that cell.
% probFree is the probability that a cell is occupied given that a laser beam passed through it.

% mapUpdate is a matrix of the same size as map. It has the log odds values that need to be added for the cells
% affected by the current laser scan. All unaffected cells should be zeros.
% robPoseMapFrame is the pose of the robot in the map coordinates frame.
% laserEndPntsMapFrame are map coordinates of the endpoints of each laser beam (also used for visualization purposes).

% Initialize mapUpdate.
size(map);
mapUpdate = zeros(size(map));

% Robot pose as a homogeneous transformation matrix.
robTrans = v2t(robPose);

% TODO: compute robPoseMapFrame. Use your world_to_map_coordinates implementation.
%%%%%%%%%%%%%%%
%robPose
%%%%%%%%%%%%%%%
robPoseMapFrame(1:2) = world_to_map_coordinates(robPose(1:2),gridSize,offset);
robPoseMapFrame(3) = robPose(3);
%%%%%%%%%%%%%%%
%robPoseMapFrame
%%%%%%%%%%%%%%%
% Compute the Cartesian coordinates of the laser beam endpoints.
% Set the third argument to 'true' to use only half the beams for speeding up the algorithm when debugging.
laserEndPnts = robotlaser_as_cartesian(scan, 30, false);

% Compute the endpoints of the laser beams in the world coordinates frame.
laserEndPnts = robTrans*laserEndPnts;

% TODO: compute laserEndPntsMapFrame from laserEndPnts. Use your world_to_map_coordinates implementation.
laserEndPntsMapFrame = world_to_map_coordinates(laserEndPnts(1:2,:),gridSize,offset);
%laserEndPntsMapFrame(3,:) = laserEndPnts(3,:);
%%%%%%%%%%%%
%laserEndPntsMapFrame
%%%%%%%%%%%%

% freeCells are the map coordinates of the cells through which the laser beams pass.
freeCells = [];

% Iterate over each laser beam and compute freeCells.
% Use the bresenham method available to you in tools for computing the X and Y
% coordinates of the points that lie on a line.
% Example use for a line between points p1 and p2:
% [X,Y] = bresenham(map,[p1_x, p1_y; p2_x, p2_y]);
% You only need the X and Y outputs of this function.
for sc=1:columns(laserEndPntsMapFrame)
%%%	sc
        %TODO: compute the XY map coordinates of the free cells along the laser beam ending in laserEndPntsMapFrame(:,sc)
	%%%%%%%%%%%%%%
	%robPoseMapFrame(1)
	%robPoseMapFrame(2)
	%%%%%%%%%%%%%%
%	[~,~,map,X,Y] = bresenham(map,[robPoseMapFrame(1), robPoseMapFrame(2); laserEndPntsMapFrame(:,sc)'],0);
	[X,Y] = bresenham2([robPoseMapFrame(1), robPoseMapFrame(2); laserEndPntsMapFrame(1,sc), laserEndPntsMapFrame(2,sc)]);
        %TODO: add them to freeCells
	freeCells = [freeCells, [X;Y]];
%	X
%	Y

endfor

%TODO: update the log odds values in mapUpdate for each free cell according to probFree.
for i=1:size(freeCells,2)
    tmpR = freeCells(1,i);
    tmpC = freeCells(2,i);
    mapUpdate(tmpR,tmpC) = prob_to_log_odds(probFree);
endfor

%TODO: update the log odds values in mapUpdate for each laser endpoint according to probOcc.
for i=1:size(laserEndPnts,2)
    mapUpdate(laserEndPntsMapFrame(1,i),laserEndPntsMapFrame(2,i)) = prob_to_log_odds(probOcc);
endfor

end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function BF_AnnotatePoints(xy,TimeSeries,annotateParams)
% BF_AnnotatePoints     Annotates time series/metadata to a plot
%
%---INPUTS:
% xy, a vector (or cell) of x-y co-ordinates of points on the plot
% TimeSeries, a table of time series making up the plot
% annotateParams, structure of custom plotting parameters

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

numTimeSeries = height(TimeSeries);

% ------------------------------------------------------------------------------
%% Set default plotting parameters:
% ------------------------------------------------------------------------------
if isfield(annotateParams,'n')
    numAnnotate = annotateParams.n;
else
    numAnnotate = 6;
end
if numAnnotate==0 % no need to go further
    return
end
% Cannot allow to annotate more points than time series
numAnnotate = min(numAnnotate,numTimeSeries);
if isfield(annotateParams,'maxL')
    maxL = annotateParams.maxL;
else
    maxL = 300; % length of annotated time series segments
end
if isfield(annotateParams,'userInput')
    userInput = annotateParams.userInput;
else
    userInput = true; % user input points rather than being randomly chosen
end
if isfield(annotateParams,'fdim')
    fdim = annotateParams.fdim;
else
    fdim = [0.30,0.08]; % width, height
end
% textAnnotation can be: 'Name' (annotate name),
% 'ID' (annotate ts_id), or 'none' (not annotation)
if isfield(annotateParams,'textAnnotation')
    if ismember(annotateParams.textAnnotation,{'Name','ID','none'})
        textAnnotation = annotateParams.textAnnotation;
    else
        warning('Unknown annotation type: %s',annotateParams.textAnnotation);
        textAnnotation = 'none'; % no annotations by default
    end
else
    textAnnotation = 'none'; % no annotations by default
end
if isfield(annotateParams,'theLineWidth')
    theLineWidth = annotateParams.theLineWidth;
else
    theLineWidth = 0.8; % % line width for time series
end
% Work through colors of the rainbow if only one class:
if isfield(annotateParams,'doRainbow')
    doRainbow = annotateParams.doRainbow;
else
    doRainbow = true;
end
plotCircle = true; % circle around annotated points

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

% Want on a common scale for finding neighbors:
% xy_std = std(xy);
% xy_mean = mean(xy);
% xy_zscore = zscore(xy);

% Initializing basic parameters/variables:
ax = gca;
pxlim = ax.XLim; % plot limits
pylim = ax.YLim; % plot limits
pWidth = diff(pxlim); % plot width
pHeight = diff(pylim); % plot height
alreadyPicked = zeros(numAnnotate,1); % record those already picked

% Groups:
[groupLabels,classLabels,groupLabelsInteger,numGroups] = TS_ExtractGroupLabels(TimeSeries);

% Colors:
if ~isfield(annotateParams,'groupColors')
    groupColors = GiveMeGroupColors(annotateParams,numGroups); % Set colors
else
    groupColors = annotateParams.groupColors;
end
if numGroups==1
    if doRainbow
        rainbowColors = [BF_GetColorMap('set1',5,1); BF_GetColorMap('dark2',6,1)];
    else
        % all black:
        rainbowColors = repmat([0,0,0],20,1);
        rainbowColors = mat2cell(rainbowColors,ones(20,1));
    end
end

%-------------------------------------------------------------------------------
% Don't use user input to select points to annotate: instead they are selected randomly
%-------------------------------------------------------------------------------
if ~userInput
    % Annotate a random set:
    alreadyPicked = randsample(numTimeSeries,numAnnotate);
end

% ------------------------------------------------------------------------------
% Go through and annotate selected points
%-------------------------------------------------------------------------------
fprintf(1,['Annotating time series segments to %u points in the plot, ' ...
                    'displaying %u samples from each...\n'],numAnnotate,maxL);

for j = 1:numAnnotate
    title(sprintf('%u points remaining to annotate',numAnnotate-j+1));

    % Get user to pick a point, then find closest datapoint to their input:
    if userInput % user input
        newPointPicked = false;
        numAttempts = 0;
        % Keep selecting until you get a new point (can't get stuck in infinite
        % loop because maximum annotations is number of datapoints). Still, you
        % get 3 attempts to pick a new point
        while ~newPointPicked && (numAttempts < 3)
            point = ginput(1);
            % point_z = (point - xy_mean)./xy_std;
            iPlot = BF_ClosestPoint_ginput(xy,point,true);
            if ~ismember(iPlot,alreadyPicked)
                alreadyPicked(j) = iPlot;
                newPointPicked = true;
            else
                beep
                numAttempts = numAttempts + 1;
            end
        end
        if numAttempts==3 % maxed out
            alreadyPicked(j) = NaN;
            continue
        end
    else
        % Use a previously-selected point (pre-selected at random)
        iPlot = alreadyPicked(j);
    end

    % (x,y) co-ords of the point to plot:
    plotPoint = xy(iPlot,:);

    % Get the group index of the selected time series:
    theGroupIndex = groupLabelsInteger(iPlot);

    % Crop the time series:
    if ~isempty(maxL)
        timeSeriesSegment = TimeSeries.Data{iPlot}(1:min(maxL,end));
    else
        timeSeriesSegment = [];
    end

    % When only one group, cycle through rainbow colors for each annotation:
    if numGroups==1
        % cycle through rainbow colors sequentially:
        groupColors{1} = rainbowColors{rem(j,length(rainbowColors)-1)+1};
    end

    % Plot a circle around the annotated point:
    if plotCircle
        plot(plotPoint(1),plotPoint(2),'o','MarkerEdgeColor',groupColors{theGroupIndex},...
                            'MarkerFaceColor',brighten(groupColors{theGroupIndex},0.5));
    end

    % Add text annotations:
    switch textAnnotation
    case 'Name'
        % Annotate text with names of datapoints:
        if ismember('Group',TimeSeries.Properties.VariableNames)
            text(plotPoint(1),plotPoint(2)-0.01*pHeight,sprintf('%s-%u',...
                    TimeSeries.Name{iPlot},TimeSeries.Group(iPlot)),...
                    'interpreter','none','FontSize',8,...
                    'color',brighten(groupColors{theGroupIndex},-0.6));
        else
            text(plotPoint(1),plotPoint(2)-0.01*pHeight,TimeSeries.Name{iPlot},...
                    'interpreter','none','FontSize',8,...
                    'color',brighten(groupColors{theGroupIndex},-0.6));
        end
    case 'ID'
        % Annotate text with ts_id:
        text(plotPoint(1),plotPoint(2)-0.01*pHeight,...
                num2str(TimeSeries.ID(iPlot)),...
                    'interpreter','none','FontSize',8,...
                    'color',brighten(groupColors{theGroupIndex},-0.6));
    case 'length'
        text(plotPoint(1),plotPoint(2)-0.01*pHeight,...
                num2str(length(TimeSeries.Data{iPlot})),...
                'interpreter','none','FontSize',8,...
                'color',brighten(groupColors{theGroupIndex},-0.6));
    end

    % Adjust if annotation goes off axis x-limits
    px = plotPoint(1)+[-fdim(1)*pWidth/2,+fdim(1)*pWidth/2];
    if px(1) < pxlim(1), px(1) = pxlim(1); end % can't plot off left side of plot
    if px(2) > pxlim(2), px(1) = pxlim(2)-fdim(1)*pWidth; end % can't plot off right side of plot

    % Adjust if annotation goes above maximum y-limits
    py = plotPoint(2)+[0,fdim(2)*pHeight];
    if py(2) > pylim(2)
        py(1) = pylim(2)-fdim(2)*pHeight;
    end

    % Annotate the time series
    if ~isempty(timeSeriesSegment)
        plot(px(1)+linspace(0,fdim(1)*pWidth,length(timeSeriesSegment)),...
                py(1)+fdim(2)*pHeight*(timeSeriesSegment-min(timeSeriesSegment))/(max(timeSeriesSegment)-min(timeSeriesSegment)),...
                    '-','color',groupColors{theGroupIndex},'LineWidth',theLineWidth);
    end

end

% Clear the title
title('')

%-------------------------------------------------------------------------------
function groupColors = GiveMeGroupColors(annotateParams,numGroups) % Set colors
    if isstruct(annotateParams) && isfield(annotateParams,'cmap')
        if ischar(annotateParams.cmap)
            groupColors = BF_GetColorMap(annotateParams.cmap,numGroups,1);
        else
            groupColors = annotateParams.cmap; % specify the cell itself
        end
    else
        if numGroups < 10
            groupColors = BF_GetColorMap('set1',numGroups,1);
        elseif numGroups <= 12
            groupColors = BF_GetColorMap('set3',numGroups,1);
        elseif numGroups <= 22
            groupColors = [BF_GetColorMap('set1',numGroups,1); ...
                        BF_GetColorMap('set3',numGroups,1)];
        elseif numGroups <= 50
            groupColors = mat2cell(jet(numGroups),ones(numGroups,1));
        else
            error('There aren''t enough colors in the rainbow to plot this many classes!')
        end
    end
    if (numGroups == 1)
        groupColors = {'k'}; % Just use black...
    end
end
%-------------------------------------------------------------------------------

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%
% Script to exemplify the use of the Matlab port of the liblbfgs-library
%

% These are all available options
options = ...
struct('m',6,...                % The number of corrections to approximate the inverse hessian matrix.
    'epsilon',1e-7,...          % A minimization terminates when ||g|| < epsilon * max(1, ||x||)
    'past',0,...                % Distance for delta-based convergence test
    'delta',1e-5,...            % Delta for convergence test.
    'MaxIter',500,...           %The maximum number of iterations.
    'linesearch','more_thuente',...  % The line search algorithm (other options are 'backtracking_armijo', 'backtracking', 'backtracking_wolfe' or 'backtracking_strong_wolfe')
    'max_linesearch',40,...     % The maximum number of trials for the line search.
    'min_step',1e-20,...        % The maximum step of the line search.
    'max_step',1e20,...         % The minimum step of the line search routine.
    'ftol',1e-4,...             % A parameter to control the accuracy of the line search routine.
    'wolfe',0.9,...             % A coefficient for the Wolfe condition.
    'gtol',0.9,...              % A parameter to control the accuracy of the line search routine.
    'xtol',1e-16,...            % The machine precision for floating-point values.
    'orthantwise_c',0,...   	% Coefficient for the L1 norm of variables.
    'orthantwise_start',0,...   % First index for the parameters subject to L1-penalty
    'orthantwise_end',-1,...    % Last index for the parameters subject to L1-penalty
    'DerivativeCheck','off',... % Derivative check using finite differences  ('on','off')
    'Display','iter');          % Available options are 'final','iter' or 'none'

N = 10000;           % number of variables
x0 = randn(N,1);      % initial guess
[x,fval,msg] = liblbfgs(@objective,x0,options);

% Example of syntax for passing additional arguments to the objective function
% [x,fval,msg] = lbfgs(@(x) objective(x,par),x0,options);

repo: REST


text:
function imResult = blendMode_Phoenix(A, B, offsetW, offsetH)
%% Phoenix blending mode: This subtracts the lighter pixel from the darker
%   pixel, and adds 1, giving a bright result.
%
% Input:
%       A       -       Base Image
%       B       -       Top Image
%   offsetW     -   move picture B horizontally in respect to the top-left
%                   corner of picture A. Default value = 1.
%   offsetH     -   move picture B vertically in respect to the top-left
%                   corner of picture A. Default value = 1.
%
% Output:
%       imResult    -   Result of the blending, having the same size of the
%                       Base Image A.
%

%% Check Input
a = size(A);
b = size(B);
blendMode_checkInput(nargin, a, b, func2str(@blendMode_Phoenix));

if nargin < 3
    offsetW = 1;
    offsetH = 1;
end

if nargin < 4
    offsetH = 1;
end

%% Implementation
imResult = A;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    [A, B] = blendMode_ResizeImages(A, B, a, b, offsetW, offsetH);
end

C = min(A, B) - max(A, B) + 1;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    imResult = blendMode_CreateResult(imResult, C, offsetW, offsetH);
else
    imResult = C;
end

end


text:
function [gKern, gVarmeans, gVarcovars] = rbfardjitVardistPsi0Gradient(rbfardjitKern, vardist, covGrad)

% RBFARDJITVARDISTPSI0GRADIENT Description

% VARGPLVM

[gKern, gVarmeans, gVarcovars] = rbfard2VardistPsi0Gradient(rbfardjitKern, vardist, covGrad);

repo: vargplvm


text:
function CPD = generic_CPD(clamped)
% GENERIC_CPD Virtual constructor for generic CPD
% CPD = discrete_CPD(clamped)

if nargin < 1, clamped = 0; end

CPD.clamped = clamped;
CPD = class(CPD, 'generic_CPD');

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% PLOT2DKM - For a 2-D binary classification problem, plot2dkm plots the data,
%            the margin and error vectors and contours of constant margin
%            for the SVM classifier in memory.
%
% Syntax: plot2dkm
%
% Version 3.22e -- Comments to diehl@alumni.cmu.edu
%

function plot2dkm

% flags for example state
MARGIN    = 1;
ERROR     = 2;
RESERVE   = 3;
UNLEARNED = 4;

% define global variables
global ind;   % cell array containing indices of margin, error, reserve and unlearned vectors
global X;     % matrix of margin, error, reserve and unlearned vectors stored columnwise
global y;     % column vector of class labels (-1/+1) for margin, error, reserve and unlearned vectors

% plot examples with label -1
figure;
indn1 = find(y == -1);
scatter(X(1,indn1),X(2,indn1),40,'b','filled');
hold on;

% plot examples with label +1
ind1 = find(y == 1);
scatter(X(1,ind1),X(2,ind1),40,'r');

% plot margin vectors
scatter(X(1,ind{MARGIN}),X(2,ind{MARGIN}),120,'k');
scatter(X(1,ind{MARGIN}),X(2,ind{MARGIN}),150,'k');
scatter(X(1,ind{MARGIN}),X(2,ind{MARGIN}),200,'k');

% plot error vectors
scatter(X(1,ind{ERROR}),X(2,ind{ERROR}),120,'k');
scatter(X(1,ind{ERROR}),X(2,ind{ERROR}),150,'k');
scatter(X(1,ind{ERROR}),X(2,ind{ERROR}),200,'k');
scatter(X(1,ind{ERROR}),X(2,ind{ERROR}),120,'k','x');
scatter(X(1,ind{ERROR}),X(2,ind{ERROR}),150,'k','x');
scatter(X(1,ind{ERROR}),X(2,ind{ERROR}),200,'k','x');

% draw margin band
xl = xlim;
yl = ylim;
pd = min(xl(2)-xl(1),yl(2)-yl(1))/100;
x_range = xl(1):pd:xl(2);
y_range = yl(1):pd:yl(2);
f = zeros(length(y_range),length(x_range));
i = 1;
for xp = x_range
   j = 1;
   for yp = y_range
      f(j,i) = svmeval([xp ; yp]);
      j = j + 1;
   end;
   i = i + 1;
end;
contour(x_range,y_range,f,[-1 1]);

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function options = svmRunOptions(varargin)
% options = svmRunOptions(varargin)
% SUPPORT VECTOR MACHINE - RUN OPTIONS
% ---------------------------------------------------------
% Produces structure with parameters describing how to perform a
% classification.
%
% INPUTS
%	OPTIONS
%       'Procedure' -
%           'KFolds' - Split trials into k approximately equal folds,
%               removing each fold one by one, training on the remainder,
%               and testing on the withheld fold.
%           'ScanFolds' - KFolds where the folds are scans/runs.
%           'LeaveOneOut' - KFolds where the folds are individual trials.
%       'Verbose' - Whether or not to print MATLAB generated output.
%           DEFAULT: false
%       'Options' - Pre-existing options structure.
%       'OptionsFile' - Path to file containing pre-existing options
%           structure.
%       'k' - K parameter for KFolds. DEFAULT: 10
%       's' | 'SVM_Type' - See libSVM help.
%       't' | 'Kernel_Type' - See libSVM help.
%       'd' | 'Degree' - See libSVM help.
%       'g' | 'Gamma' - See libSVM help.
%       'r' | 'Coef0' - See libSVM help.
%       'c' | 'Cost' - See libSVM help.
%       'n' | 'Nu' - See libSVM help.
%       'p' | 'Epsilon_Loss' - See libSVM help.
%       'm' | 'CacheSize' - See libSVM help.
%       'e' | 'Epsilon' - See libSVM help.
%       'h' | 'Shrinking' - See libSVM help.
%       'b' | 'Probability_Estimates' - See libSVM help.
%       'w' | 'Weight' - See libSVM help.
%       'q' | 'Quiet' - See libSVM help.
%
% OUTPUTS
%   options - Structure containing all parameters to be used in svmRun, or
%       elsewhere.
%
% USAGE
%   Setting up the parameters to run a classification with KFolds, k = 20.
%   This options structure is then passed to svmRun.
%       options = svmRunOptions('Procedure', 'KFolds', 'k', 20);
%
%   Often, this function is called from within another function, such as
%   svmRun, which will accept the same inputs.
%       svmRun(..., 'Procedure', 'KFolds', 'k', 20);
%
% See also SVMINIT, SVMRUN, SVMEXPORTMAP, SVMBOOTSTRAP, SVMRELABEL,
% SVMREMOVE, SLINIT.
%
% renobowen@gmail.com [2010]
%

    % Lib SVM options
    options.svm_type = [];
    options.kernel_type = 0;
    options.degree = [];
    options.gamma = [];
    options.coef0 = [];
    options.cost = [];
    options.nu = [];
    options.epsilon_loss = [];
    options.cachesize = [];
    options.epsilon = [];
    options.shrinking = [];
    options.probability_estimates = [];
    options.weight = [];
    options.quiet = 1;
    options.subsample = false;

    % MATLAB side options
    options.procedure = 'scanfolds';
    options.verbose = false;
    options.k = 10;

    i = 1;
    while (i <= length(varargin))
        if (isempty(varargin{i})), break; end
        switch (lower(varargin{i}))
            case {'options'} % careful with this one and the next - we assume you know what you're doing with the options struct
                options = varargin{i + 1};
            case {'optionsfile'}
                load(varargin{i + 1});
            case {'procedure'}
                options.procedure = varargin{i + 1};
            case {'k'}
                options.k = varargin{i + 1};
            case {'log2cvector'}
                options.log2cvector = varargin{i + 1};
            case {'log2gvector'}
                options.log2gvector = varargin{i + 1};
            case {'svm_type' 's'}
                options.svm_type = varargin{i + 1};
            case {'kernel_type' 't'}
                options.kernel_type = varargin{i + 1};
            case {'degree' 'd'}
                options.degree = varargin{i + 1};
            case {'gamma' 'g'}
                options.gamma = varargin{i + 1};
            case {'coef0' 'r'}
                options.coef0 = varargin{i + 1};
            case {'cost' 'c'}
                options.coef = varargin{i + 1};
            case {'nu' 'n'}
                options.nu = varargin{i + 1};
            case {'epsilon_loss' 'p'}
                options.epsilon_loss = varargin{i + 1};
            case {'cachesize' 'm'}
                options.cachesize = varargin{i + 1};
            case {'epsilon' 'e'}
                options.epsilon = varargin{i + 1};
            case {'shrinking' 'h'}
                options.shrinking = varargin{i + 1};
            case {'probability_estimates', 'b'}
                options.probability_estimates = varargin{i + 1};
            case {'weight' 'w'}
                options.weight = [varargin{i + 1} varargin{i + 2}];
                i = i + 1;
            case {'quiet' 'q'} % Silence LIBSVM side outputs
                options.quiet = 1;
                i = i - 1;
            case {'verbose'}
                options.verbose = true;
                i = i - 1;
            case {'silent'} % Silence MATLAB side outputs
                options.verbose = false;
                i = i - 1;
            case {'subsample'}
                subsample = varargin{i+1};
            otherwise
                fprintf(1, 'Unrecognized option: ''%s''\n', varargin{i});
                return;
        end
        i = i + 2;
    end

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% GRADIENT STEP function.
% Computes the cost and gradient of the neural network for unrolled theta parameters.
function [J gradients] = nn_gradient_step(nn_params_unrolled, layers, X, y, lambda)
    % Reshape nn_params back into the matrix parameters.
    nn_params = nn_params_roll(nn_params_unrolled, layers);

    % Compute the cost.
    J = nn_cost_function(nn_params, layers, X, y, lambda);

    % Do backpropagation.
    gradients = nn_backpropagation(nn_params, layers, X, y, lambda);
end

author: trekhleb

repo: machine-learning-octave


text:
function pass = test_firstOrderIntegralEqn(~)

%% Setup, testing a first order integral equation
d = [0 5];
L = chebop(@(u) diff(u) + 2*u + 5*cumsum(u), d);

%% Two different ways of specifying the LBC, should give the same results:
L.lbc = 0;
u1 = L\1;

L.lbc = [];
L.bc = @(x,u) u(0);
u2 = L\1;

%% Did we get the same solution?
pass = (norm(u1-u2) == 0);

%% sum instead of cumsum, rbc
L = chebop(@(x,u) diff(u) + 2*sum((5-x)*sin(x)*u),d);
x = chebfun('x', L.domain);
L.rbc = 1;
rhs = sin(pi*x);
u = L\rhs;
pass(2) = (norm(L*u - rhs,inf) < 1e-10) && (u(end) - 1 < 1e-10);

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [xResampleV,yResampleV,zResampleV] = ...
    getResampledGrid(resampResolutionV,xValsV,yValsV,zValsV,...
    originV,gridAlignMethod,varargin)
%getResampledGrid_v2
% ------------------------------------------------------------------------
% INPUTS
% resampResolutionV  : Output voxel spacing in cm [dx dy dz]
% xValsV             : x coordinates of voxel centers in original scan
% yValsV             : y coordinates of voxel centers in original scan
% zValsV             : z coordinates of voxel centers in original scan
%------------------------------------------------------------------------
%Ref: https://arxiv.org/pdf/1612.07003.pdf
%AI 12/01/22

if ~exist('method','var')
    gridAlignMethod = 'center';
end

% Set default perturbation offsets:
if nargin<7
    perturbX = 0;
    perturbY = 0;
    perturbZ = 0;
else
    perturbV = varargin{1};
    perturbX = perturbV(1);
    perturbY = perturbV(2);
    perturbZ = perturbV(3);
end

%% Define voxel spacing
dx = abs(median(diff(xValsV)));
dy = abs(median(diff(yValsV)));
dz = abs(median(diff(zValsV)));
origResolutionV = [dx dy dz];

%Correct to match DICOM convention
origResolutionV(3) = -origResolutionV(3);
resampResolutionV(3) = -resampResolutionV(3);

if length(resampResolutionV)==3 && ~isnan(resampResolutionV(3))
    resamp3DFlag = 1;
else
    %Resample in-plane
    resamp3DFlag = 0;
end

%% No. voxels
origSizeV = [length(xValsV) length(yValsV) length(zValsV)];
resampSizeV = ceil( origSizeV.*origResolutionV ./ resampResolutionV);

switch(gridAlignMethod)

    case 'center'

        %Get output grid origin
        % In world coordinates:
        resampOriginV = originV + (origResolutionV.*(origSizeV-1) - ...
            resampResolutionV.*(resampSizeV-1))/2;

        %In grid co-ords:
        %resampOriginV = 0.5* (origSizeV- 1 -
        %resampResolutionV.*(resampSizeV-1)/origResolutionV);

        %Generate output grid
        xResampleV = resampOriginV(1):resampResolutionV(1):...
            resampOriginV(1)+(resampSizeV(1)-1)*resampResolutionV(1);

        yResampleV = resampOriginV(2):resampResolutionV(2):...
            (resampOriginV(2)+(resampSizeV(2)-1)*resampResolutionV(2));

        if resamp3DFlag
            zResampleV = resampOriginV(3):resampResolutionV(3):...
                (resampOriginV(3)+(resampSizeV(3)-1)*resampResolutionV(3));
            zResampleV = flip(zResampleV);
        else
            zResampleV = zValsV;
        end

    otherwise
        error('Unsupported grid alignment method %s',gridAlignMethod)
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function varargout=frq2midi(varargin)
%
% For calling details please see v_frq2midi.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_frq2midi.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: frq2midi.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_frq2midi(varargin{:});
else
    v_frq2midi(varargin{:});
end


text:
function dtiCompileMex

oldDir = pwd;
rootdir=fileparts(fileparts(which('mrDiffusion')));

mexFiles = {};

cd(fileparts(which('dtiFiberTracker.cxx')));
disp('Compiling dtiFiberTracker.cxx ...');
try
    mex -O -I./jama dtiFiberTracker.cxx
    mexFiles{end+1} = ['dtiFiberTracker.' mexext];
catch
    disp('Compiling dtiFiberTracker.cxx FAILED...');
end

cd(fileparts(which('nearpoints.cxx')));
disp('Compiling nearpoints.cxx ...');
try
    mex -O nearpoints.cxx
    mexFiles{end+1} = ['nearpoints.' mexext];
catch
    disp('Compiling nearpoints.cxx FAILED...');
end

cd(fileparts(which('magicwand1.c')));
disp('Compiling magicwand1.c ...');
try
    mex -O magicwand1.c
    mexFiles{end+1} = ['magicwand1.' mexext];
catch
    disp('Compiling magicwand1.c FAILED...');
end

cd(fileparts(which('dtiJointHist.c')));
disp('Compiling dtiJointHist.c ...');
try
    mex -O dtiJointHist.c
    mexFiles{end+1} = ['dtiJointHist.' mexext];
catch
    disp('Compiling dtiJointHist.c FAILED...');
end

cd(fileparts(which('niftiRead.c')));
disp('Compiling read/writeFileNifti ...');
try
    mex niftiRead.c nifti1_io.c znzlib.c zlib/adler32.c zlib/compress.c zlib/crc32.c zlib/deflate.c zlib/gzio.c zlib/infback.c zlib/inffast.c zlib/inflate.c zlib/inftrees.c zlib/trees.c zlib/zutil.c
    mexFiles{end+1} = ['niftiRead.' mexext];
    mex writeFileNifti.c nifti1_io.c znzlib.c zlib/adler32.c zlib/compress.c zlib/crc32.c zlib/deflate.c zlib/gzio.c zlib/infback.c zlib/inffast.c zlib/inflate.c zlib/inftrees.c zlib/trees.c zlib/zutil.c
    mexFiles{end+1} = ['writeFileNifti.' mexext];
catch
    disp('Compiling read/writeFileNifti.c FAILED...');
end

%      disp('Compiling dtiTensorInterp_Pajevic.c ...');
%      if(strcmp(computer,'GLNX86'))
%          mex -O dtiTensorInterp_Pajevic.c ./bcadtlib.so
%          mexFiles{end+1} = ['dtiTensorInterp_Pajevic.' mexext];
%      else
%          warning('dtiTensorInterp_Pajevic.c can only be compiled under 32-bit linux! Skipping it...'); %#ok<WNTAG>
%      end

%Confusing, multiple copies of sources and/or mex files in a different
%directory
%(1) 'mrManDist'; %Csource: c source
%(2) 'dijkstra2';%Csource under ManifoldUtilities AND within dir mrFlatMesh
%(3) 'dijkstra';%Csource

%'curvature'; -- where is the source?
%'smooth_mesh'; -- where is the source?
%'build_mesh'; -- where is the source?

%List of files to compile with path relative from the root vistasoft dir
functionsToCompile={
    ['mrAlign' filesep 'regInplanes'];
    ['mrAlign' filesep 'regHistogram'];
    ['mrAnatomy' filesep 'MrGray' filesep 'GrowGrayMatter' filesep 'grow_gray']; %cpp
    ['mrAnatomy' filesep 'mrMesh' filesep 'tcpToolbox' filesep 'pnet'];
    ['mrAnatomy' filesep 'VolumeUtilities' filesep 'mrAnatFastInterp3'];
    ['mrAnatomy' filesep 'VolumeUtilities' filesep 'myCinterp3'];
    ['mrDiffusion' filesep 'src' filesep 'dtiFitTensor'];
    ['mrDiffusion' filesep 'src' filesep 'dtiSplitTensor'];
    ['mrDiffusion' filesep 'xform' filesep 'fastDeformation' filesep 'trilin'];
    ['mrDiffusion' filesep 'xform' filesep 'fastDeformation' filesep 'initInvDiagsInC'];
    ['mrLoadRet' filesep 'Analysis' filesep 'SignalProc' filesep 'sumOfNeighbors'];
    ['mrLoadRet' filesep 'Analysis' filesep 'MotionComp' filesep 'MI' filesep 'MI' filesep 'Joint_Histogram' filesep 'spm_hist2_weighted_MI'];
    ['mrLoadRet' filesep 'Analysis' filesep 'MotionComp' filesep 'MI' filesep 'MI' filesep 'Joint_Histogram' filesep 'spm_hist2_weighted'];
    ['mrDiffusion' filesep 'src' filesep 'ndfun'];
    };

for CfileI =1:size(functionsToCompile, 1)
    cext='.c';
    CfileName=functionsToCompile{CfileI};
    [srcDir, f] = fileparts(which(CfileName));
    if(~exist(srcDir, 'dir'))
        warning(['Something went wrong! Where is ' srcDir ' directory?']); %#ok<WNTAG>
    end
    cd(srcDir);
    if(~exist([f cext], 'file'))
        cext='.cpp';
    end

    try
        disp(['Compiling ' f '...']);
        cmd = ['mex -O ' f cext];
        if(ispc)
            compileWithLapackLibToWinCommmand(cmd);
        else
            eval([cmd ' -lmwlapack -lmwblas']);
        end
        mexFiles{end+1} = [f '.' mexext];

    catch
        warning('Something went wrong! Check mex messages above.'); %#ok<WNTAG>
    end

end

fprintf('%d files compiled successfully.\n', length(mexFiles));

if(~isempty(mexFiles))
    fprintf('\n\n');
    fprintf('The following files were compiled and placed into the directories where respective source files were located%s:\n\n');
    for ii=1:length(mexFiles); fprintf('  %s\n', mexFiles{ii}); end
    %fprintf('\n\nYou probably want to move them to the usual place for mex files\n');
    %fprintf('(maybe %s ?)\n\n', fullfile(fileparts(srcDir),'DLLXX'));
end
cd(oldDir);

return;

function compileWithLapackLibToWinCommmand(cmd)

if exist(fullfile(matlabroot, 'extern', 'lib', 'win32', 'microsoft', 'libmwlapack.lib'),'file')
    libFile = fullfile(matlabroot, 'extern', 'lib', 'win32', 'microsoft', 'libmwlapack.lib');
elseif exist(fullfile(matlabroot, 'extern', 'lib', 'win32', 'microsoft', 'msvc60', 'libmwlapack.lib'),'file')
    libFile = fullfile(matlabroot, 'extern', 'lib', 'win32', 'microsoft', 'msvc60', 'libmwlapack.lib');
elseif exist(fullfile(matlabroot, 'extern', 'lib', 'win64', 'microsoft', 'libmwlapack.lib'),'file')
    libFile = fullfile(matlabroot, 'extern', 'lib', 'win64', 'microsoft', 'libmwlapack.lib');
else
    error('Cannot find libmwlapack.lib in standarad matlab directories!');
end

if isempty(find(matlabroot) == ' ')
    cmd = [cmd ' ' libFile];
    eval(cmd);
else
    copyfile(libFile);
    cmd = [cmd ' libmwlapack.lib'];
    eval(cmd);
    delete('libmwlapack.lib');
end

return;

% /home/bob/svn/vistasoft/trunk/fileFilters/nifti/niftiRead.mexa64
% /home/bob/svn/vistasoft/trunk/fileFilters/nifti/matToQuat.mexa64
% /home/bob/svn/vistasoft/trunk/fileFilters/nifti/writeFileNifti.mexa64
% /home/bob/svn/vistasoft/trunk/mrAlign/regInplanes.mexa64
% /home/bob/svn/vistasoft/trunk/mrAlign/regHistogram.mexa64
% /home/bob/svn/vistasoft/trunk/mrAnatomy/MrGray/GrowGrayMatter/grow_gray.mexa64
%
% /home/bob/svn/vistasoft/trunk/mrAnatomy/mrMesh/curvature.mexa64
% /home/bob/svn/vistasoft/trunk/mrAnatomy/mrMesh/tcpToolbox/pnet.mexa64
% /home/bob/svn/vistasoft/trunk/mrAnatomy/mrMesh/smooth_mesh.mexa64
% /home/bob/svn/vistasoft/trunk/mrAnatomy/mrMesh/build_mesh.mexa64
%
% /home/bob/svn/vistasoft/trunk/mrAnatomy/VolumeUtilities/nearpoints.mexa64
% /home/bob/svn/vistasoft/trunk/mrAnatomy/VolumeUtilities/mrAnatFastInterp3.mexa64
% /home/bob/svn/vistasoft/trunk/mrAnatomy/VolumeUtilities/myCinterp3.mexa64
% /home/bob/svn/vistasoft/trunk/mrAnatomy/ManifoldUtilities/mrManDist.mexa64
% /home/bob/svn/vistasoft/trunk/mrAnatomy/ManifoldUtilities/Csource/dijkstra2.mexa64
% /home/bob/svn/vistasoft/trunk/mrAnatomy/ManifoldUtilities/ShortestPaths/dijkstra.mexa64
%X /home/bob/svn/vistasoft/trunk/mrDiffusion/src/dtiFitTensor.mexa64
%X /home/bob/svn/vistasoft/trunk/mrDiffusion/src/dtiSplitTensor.mexa64
%X /home/bob/svn/vistasoft/trunk/mrDiffusion/src/dtiFiberTracker.mexa64
%X/home/bob/svn/vistasoft/trunk/mrDiffusion/src/ndfun.mexa64
% X/home/bob/svn/vistasoft/trunk/mrDiffusion/src/magicwand1.mexa64
% X/home/bob/svn/vistasoft/trunk/mrDiffusion/src/dtiJointHist.mexa64
% /home/bob/svn/vistasoft/trunk/mrDiffusion/xform/fastDeformation/trilin.mexa64
% /home/bob/svn/vistasoft/trunk/mrDiffusion/xform/fastDeformation/initInvDiagsInC.mexa64
% /home/bob/svn/vistasoft/trunk/mrLoadRet/Analysis/SignalProc/sumOfNeighbors.mexa64
% /home/bob/svn/vistasoft/trunk/mrLoadRet/Analysis/MotionComp/MI/MI/Joint_Histogram/spm_hist2_weighted_MI.mexa64
% /home/bob/svn/vistasoft/trunk/mrLoadRet/Analysis/MotionComp/MI/MI/Joint_Histogram/spm_hist2_weighted.mexa64


text:
function u=proj(x,v,y)
%
% u=proj(x,v,y);
% Projects piecewise constant function v with
% discontinuities at y onto
% piecewise constant function u with discontinuities at x.
%
n=length(x);
l=length(y);
if (l==0),
	u=v(1)*ones(size(x)-1);
else
  u=zeros(1,n-1);
  j=1;
  h=x(1);
  ns=y(j);
  while (ns<h),
    j=j+1;
    ns=y(j);
  end;
  for i=1:n-1,
    dx=x(i+1)-h;
    while (ns<x(i+1)),
      u(i)=u(i)+v(j)*(y(j)-h)/dx;
      h=ns;
      j=j+1;
      if ((j<=l)&&(j>0)),
        ns=y(j);
      else
        ns=x(n)+1;
      end;
    end;
    u(i)=u(i)+v(j)*(x(i+1)-h)/dx;
    h=x(i+1);
  end;
end;

repo: Aero-matlab


text:
function x = mogSample(model, numSamples)

% MOGSAMPLE Sample from a mixture of Gaussians model.
% FORMAT
% DESC samples from a mixture of Gaussians.
% ARG model  : the model that you want to sample from.
% ARG numSamples : the number of samples required.
% RETURN x : the samples from the model.
%
% COPYRIGHT : Neil D. Lawrence, 2008
%
% SEEALSO : mogCreate

% MLTOOLS

p = rand(numSamples, 1);
bins = cumsum(model.prior);
compNo = sum(repmat(p, 1, model.m)<repmat(bins, numSamples, 1), 2);
x = zeros(numSamples, model.d);

for i = 1:model.m
  ind = find(compNo == i);
  x(ind, :) = repmat(model.mean(i, :), length(ind), 1);
  switch model.covtype
   case 'ppca'
    samps = randn(length(ind), model.q);
    samps = samps*model.W{i}';
    samps = samps + randn(length(ind), model.d)*sqrt(model.sigma2(i));
    x(ind, :) = x(ind, :) + samps;
   case 'spherical'
    x(ind, :) = x(ind, :) + randn(length(ind), model.d)* ...
        sqrt(model.sigma(i));
  end
end


text:
function OUT = MovCorr(X,window,i)
% =======================================================================
% Computes correlation between X_i and X_j for all j different
% of i, from window to the number of observations.  The
% first window rows are NaN.
% =======================================================================
% OUT = MovCorr(X, window, i)
% -----------------------------------------------------------------------
% INPUT
%   - X: panel time series T observations x N variables
%   - window : size of the moving window
%   - i: the variable X_i against which correlations should be returned
% -----------------------------------------------------------------------
% OUTPUT
%   - OUT: matrix with moving correlation T observations x N-1
%       variables. The first window-1 rows are NaN.
% =======================================================================
% EXAMPLE
% X = rand(100,5);
% OUT = MovCorr(X,20,1)
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015

[nobs,nvars] = size(X);
OUT = nan(nobs,nvars-1);

for tt = window:nobs
    C = corrcoef(X(tt-window+1:tt, :));
    idx = setdiff(1:nvars, [i]);
    OUT(tt, :) = C(i, idx);
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
clear all; close all; clc;
%% TO PERFORM SPEECH ENHANCEMENT:
%
%% 1. Place the noisy speech files in 'set/test_noisy_speech'.
%
%% 2. Run this script.
%
%% 3. The enhanced speech outputs will be in 'out'.
%
%% Note: the script will finish with 'Inference complete.'.
%

if isunix
	matlab_bin = [matlabroot, '/bin/glnxa64'];
	command = ['sudo mv ', matlab_bin, '/libexpat.so.1 ', matlab_bin, '/libexpat.so.1.NOFIND'];
	system(command);
end

gpu_flag = lower(input('Use a GPU? (y/n):', 's'));
if ismember(gpu_flag, {'y', 'yes'})
    gpu_flag = true;
elseif ismember(gpu_flag, {'n', 'no'})
    gpu_flag = false;
else
    error('Please enter either ''y'' or ''n''.')
end

if gpu_flag
    system('source ~/venv/DeepXi/bin/activate; ./run.sh VER=resnet-1.1n INFER=1 GAIN=mmse-lsa');
else
    system('source ~/venv/DeepXi/bin/activate; ./run.sh VER=resnet-1.1n INFER=1 GAIN=mmse-lsa GPU=9999');
end

author: anicolson

repo: DeepXi

save_path: github-repos/MATLAB/anicolson-DeepXi


text:
function Loglike = Loglike_MixHP(Seqs, model, alg)

% given responsibility, calcuate the expected number of sequence belonging
% to the k-th cluster
EX = model.R;

% update parameters of Hawkes processes (mu_k, A_k), k=1,...,N
% initialize
A = model.beta;%random('exp', model.beta);
mu = sqrt(pi/2)*model.b;%random('rayl', model.b);

tmp1 = A(:)./(model.beta(:));
tmp1(isnan(tmp1)) = 0;
tmp1(isinf(tmp1)) = 0;

tmp2 = mu(:).^2./(2*model.b(:).^2);
tmp2(isnan(tmp2)) = 0;
tmp2(isinf(tmp2)) = 0;

tmp3 = log(mu(:));
tmp3(isnan(tmp3)) = 0;
tmp3(isinf(tmp3)) = 0;

Loglike = sum(tmp1)+ sum(tmp2)-sum(tmp3);

% E-step: evaluate the responsibility using the current parameters
for c = 1:length(Seqs)

    Time = Seqs(c).Time;
    Event = Seqs(c).Mark;
    Tstart = Seqs(c).Start;

    if isempty(alg.Tmax)
        Tstop = Seqs(c).Stop;
    else
        Tstop = alg.Tmax;
        indt = Time < alg.Tmax;
        Time = Time(indt);
        Event = Event(indt);
    end

    N = length(Time);
    % calculate the integral decay function in the log-likelihood function
    G = Kernel_Integration(Tstop - Time, model);

    for i = 1:N

        ui = Event(i);
        ti = Time(i);

        lambdai = mu(ui,:)+eps;
        pii = lambdai;

        if i>1
            tj = Time(1:i-1);
            uj = Event(1:i-1);

            gij = Kernel(ti-tj, model);
            auiuj = A(uj, :, :, ui);
            pij = repmat(gij, [1,1,model.K,1]).* auiuj;

            tmp = sum(sum(pij,1),2);
            lambdai = lambdai + tmp(:)';

        end

        LL = LL+log(lambdai);

    end
    LL = LL - (Tstop-Tstart).*sum(mu);
    tmp = sum(sum(repmat(G, [1,1,model.K]).*sum(A(Event,:,:,:),4),2),1);
    LL = LL - tmp(:)';

    Loglike = Loglike - EX(c,:)*LL(:);

end

Loglike = -Loglike;

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function p = prior_sqinvsinvchi2(varargin)
%PRIOR_SQINVSINVCHI2  Scaled-Inv-Chi^2 prior structure
%
%  Description
%    P = PRIOR_SQINVSINVCHI2('PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    creates Scaled-Inv-Chi^2 prior structure for square inverse of
%    the parameter in which the named parameters have the specified
%    values. Any unspecified parameters are set to default values.
%
%    P = PRIOR_SQINVSINVCHI2(P, 'PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    modify a prior structure with the named parameters altered
%    with the specified values.
%
%    Parameterisation is done by Bayesian Data Analysis,
%    second edition, Gelman et.al 2004.
%
%    Parameters for Scaled-Inv-Chi^2 [default]
%      s2       - scale squared (variance) [1]
%      nu       - degrees of freedom [4]
%      s2_prior - prior for s2 [prior_fixed]
%      nu_prior - prior for nu [prior_fixed]
%
%  See also
%    PRIOR_*

% Copyright (c) 2000-2001,2010,2012 Aki Vehtari
% Copyright (c) 2010 Jaakko Riihim�ki

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

  ip=inputParser;
  ip.FunctionName = 'PRIOR_SQINVSINVCHI2';
  ip.addOptional('p', [], @isstruct);
  ip.addParamValue('s2',1, @(x) isscalar(x) && x>0);
  ip.addParamValue('s2_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('nu',4, @(x) isscalar(x) && x>0);
  ip.addParamValue('nu_prior',[], @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  p=ip.Results.p;

  if isempty(p)
    init=true;
    p.type = 'SqInv-S-Inv-Chi2';
  else
    if ~isfield(p,'type') && ~isequal(p.type,'SqInv-S-Inv-Chi2')
      error('First argument does not seem to be a valid prior structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('s2',ip.UsingDefaults)
    p.s2 = ip.Results.s2;
  end
  if init || ~ismember('nu',ip.UsingDefaults)
    p.nu = ip.Results.nu;
  end
  % Initialize prior structure
  if init
    p.p=[];
  end
  if init || ~ismember('s2_prior',ip.UsingDefaults)
    p.p.s2=ip.Results.s2_prior;
  end
  if init || ~ismember('nu_prior',ip.UsingDefaults)
    p.p.nu=ip.Results.nu_prior;
  end

  if init
    % set functions
    p.fh.pak = @prior_sqinvsinvchi2_pak;
    p.fh.unpak = @prior_sqinvsinvchi2_unpak;
    p.fh.lp = @prior_sqinvsinvchi2_lp;
    p.fh.lpg = @prior_sqinvsinvchi2_lpg;
    p.fh.recappend = @prior_sqinvsinvchi2_recappend;
  end

end

function [w, s] = prior_sqinvsinvchi2_pak(p)

  w=[];
  s={};
  if ~isempty(p.p.s2)
    w = log(p.s2);
    s=[s; 'log(SqInv-Sinvchi2.s2)'];
  end
  if ~isempty(p.p.nu)
    w = [w log(p.nu)];
    s=[s; 'log(SqInv-Sinvchi2.nu)'];
  end
end

function [p, w] = prior_sqinvsinvchi2_unpak(p, w)

  if ~isempty(p.p.s2)
    i1=1;
    p.s2 = exp(w(i1));
    w = w(i1+1:end);
  end
  if ~isempty(p.p.nu)
    i1=1;
    p.nu = exp(w(i1));
    w = w(i1+1:end);
  end
end

function lp = prior_sqinvsinvchi2_lp(x, p)

  lJ = -log(x)*3 + log(2);  % log(-2/x^3) log(|J|) of transformation
  xt = x.^-2;               % transformation
  lp = -sum((p.nu./2+1) .* log(xt) + (p.s2.*p.nu./2./xt) + (p.nu/2) .* log(2./(p.s2.*p.nu)) + gammaln(p.nu/2)) + sum(lJ);

  if ~isempty(p.p.s2)
    lp = lp + p.p.s2.fh.lp(p.s2, p.p.s2) + log(p.s2);
  end
  if ~isempty(p.p.nu)
    lp = lp + p.p.nu.fh.lp(p.nu, p.p.nu) + log(p.nu);
  end
end

function lpg = prior_sqinvsinvchi2_lpg(x, p)
  lJg = -3./x;              % gradient of log(|J|) of transformation
  xt  = x.^-2;              % transformation
  xtg = -2/x.^3;             % derivative of transformation
  lpg = xtg.*(-(p.nu/2+1)./xt +p.nu.*p.s2./(2*xt.^2)) + lJg;

  if ~isempty(p.p.s2)
    lpgs2 = (-sum(p.nu/2.*(1./xt-1./p.s2)) + p.p.s2.fh.lpg(p.s2, p.p.s2)).*p.s2 + 1;
    lpg = [lpg lpgs2];
  end
  if ~isempty(p.p.nu)
    lpgnu = (-sum(0.5*(log(xt) + p.s2./xt + log(2./p.s2./p.nu) - 1 + digamma1(p.nu/2))) + p.p.nu.fh.lpg(p.nu, p.p.nu)).*p.nu + 1;
    lpg = [lpg lpgnu];
  end
end

function rec = prior_sqinvsinvchi2_recappend(rec, ri, p)
% The parameters are not sampled in any case.
  rec = rec;
  if ~isempty(p.p.s2)
    rec.s2(ri,:) = p.s2;
  end
  if ~isempty(p.p.nu)
    rec.nu(ri,:) = p.nu;
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
load 'sdae_mnist_vis.mat';

colors = colormap;

figure;
hold on;
for c = 1:10
    x = H(X_labels == (c-1), 1);
    y = H(X_labels == (c-1), 2);
    rndidx = randperm(length(x));
    x = x(rndidx(1:500));
    y = y(rndidx(1:500));

    plot(x, y, 'x', 'Color', colors(ceil(c/10 * 64), :));
end
hold off;

author: kyunghyuncho

repo: deepmat


text:
function pde = oscdiffdata3

pde = struct('f',1,'d',@d);

    function K =  d(p)  % Diffusion constant
    x = p(:,1); y = p(:,2); z = p(:,3);
    K = 2*(2+sin(10*pi*x).*cos(10*pi*y)*sin(20*pi*z));
    end
end

repo: ifem


text:
function plotData(X, y)
%PLOTDATA Plots the data points X and y into a new figure
%   PLOTDATA(x,y) plots the data points with + for the positive examples
%   and o for the negative examples. X is assumed to be a Mx2 matrix.

% Create New Figure
figure; hold on;

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the positive and negative examples on a
%               2D plot, using the option 'k+' for the positive
%               examples and 'ko' for the negative examples.
%

ids = find(y==0);
Xo = X(ids,:);
plot(Xo(:,1), Xo(:,2), 'ko');

ids = find(y==1);
Xo = X(ids,:);
plot(Xo(:,1), Xo(:,2), 'k+');

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

hold off;

end

repo: machine-learning-coursera


text:
function kSpaceLow = updownsample( kSpaceHigh, lInputImage, out_size, lAAfilter)
%
% updownsample - up-sample or down-sample an input series using fourier domain
%                input series needs to be continuous of a high degree
%
%
% input:    in_m                - input matrix for up/down sampling, in fourier domain
%
%           out_size            - desired number of pixels in the output image

%
% output:   out_m               - up/down sampled image
%
% NOTE: it is important to specify if the image is REAL or COMPLEX, since
%       if the image is REAL -> we have to use ABS() on the inverse fourier
%       transform (because of roundoff errors of the transform).
%
% NOTE: since a desired amount of pixels is needed at the output, there is
%       no attempt to use matrices which are in size of power of 2. this
%       optimization can not be used in this case
%
% NOTE: input series needs to be CONTINUOUS OF A HIGH DEGREE, since the
%       upsampling is done in the frequency domain, which samples the output
%       grid with SINE-like (harmonic) functions
%
%
%
% Theory:   the upsampling is done by zero-padding in the input domain BETWEEN the samples,
%           then at the fourier domain, taking the single spectrum (out of the repetition of spectrums)
%           i.e. low pass with Fcutoff=PI/upsample_factor, zeroing the rest of the spectrum
%           and then doing ifft to the distribution.
%           since we have a zero padding operation in time, we need to multiply by the fourier gain.
%
%              +-----------+     +-------+     +---------+     +--------+     +--------+
%   y[n,m] --> | up-sample | --> |  FFT  | --> |   LPF   | --> | * Gain | --> |  IFFT  | --> interpolated
%              | factor M  |     +-------+     | Fc=PI/M |     +--------+     +--------+
%              +-----------+                   +---------+
%
%           this operation is the same as the following one (which has less operations):
%
%              +-------+     +--------+     +--------------+     +--------+
%   y[n,m] --> |  FFT  | --> | * Gain | --> | Zero Padding | --> |  IFFT  | --> interpolated
%              +-------+     +--------+     +--------------+     +--------+
%
%           NOTE THAT, the zero-padding must be such that the D.C. ferquency remains the D.C. frequency
%           and that the zero padding is applied to both positive and negative frequencies.
%           The zero padding actually condences the frequency -> which yields a longer series in the
%           image domain, but without any additional information, thus the operation must be an interpolation.

if(nargin < 4)
    lAAfilter = true;
end

% ==============================================
% get input image size, and calculate the gain
% ==============================================
out_x_sz = out_size(1); out_y_sz = out_size(2); out_z_sz = out_size(3);
[in_x_sz,in_y_sz, in_z_sz, nCha, nTime] = size( kSpaceHigh );
gain_x = out_x_sz/in_x_sz;
gain_y = out_y_sz/in_y_sz;
gain_z = out_z_sz/in_z_sz;

% upsample or downsample as needed
% ==================================
if(lInputImage)
    % image -> kSpace
    kSpaceHigh = fftnshift(kSpaceHigh,1:3);
end

% check if up/down sampling is needed at all
if(gain_x == 1 && gain_y == 1 && gain_z == 1)
    kSpaceLow = kSpaceHigh;
    dImgLow = ifftnshift(kSpaceLow,1:3);
    return;
end

% if downsampling -> apply anti-aliasing filter
if(gain_x < 1 && gain_y < 1 && gain_z < 1 && lAAfilter)
    kSpaceHigh = kSpaceHigh .* repmat(windowND(@hamming, [size(kSpaceHigh,1), size(kSpaceHigh,2), size(kSpaceHigh,3)]),[1 1 1 size(kSpaceHigh,4) size(kSpaceHigh,5)]);
end

% build grid vectors for the up/down sampling
% ============================================
% if the input is even & output is odd-> use floor for all
% if the output is even & input is odd -> use ceil for all
% other cases - don't care
% for downsampling -> the opposite
if (~mod( in_x_sz,2 ) && (out_x_sz>in_x_sz)) || (mod( in_x_sz,2 ) && (out_x_sz<in_x_sz))
    x_output_space  = max(floor((out_x_sz-in_x_sz)/2),0) + [1:min(in_x_sz,out_x_sz)];
    x_input_space   = max(floor((in_x_sz-out_x_sz)/2),0) + [1:min(in_x_sz,out_x_sz)];
else
    x_output_space  = max(ceil((out_x_sz-in_x_sz)/2),0) + [1:min(in_x_sz,out_x_sz)];
    x_input_space   = max(ceil((in_x_sz-out_x_sz)/2),0) + [1:min(in_x_sz,out_x_sz)];
end
if (~mod( in_y_sz,2 ) && (out_y_sz>in_y_sz)) || (mod( in_y_sz,2 ) && (out_y_sz<in_y_sz))
   y_output_space  = max(floor((out_y_sz-in_y_sz)/2),0) + [1:min(in_y_sz,out_y_sz)];
   y_input_space   = max(floor((in_y_sz-out_y_sz)/2),0) + [1:min(in_y_sz,out_y_sz)];
else
   y_output_space  = max(ceil((out_y_sz-in_y_sz)/2),0) + [1:min(in_y_sz,out_y_sz)];
   y_input_space   = max(ceil((in_y_sz-out_y_sz)/2),0) + [1:min(in_y_sz,out_y_sz)];
end
if (~mod( in_z_sz,2 ) && (out_z_sz>in_z_sz)) || (mod( in_z_sz,2 ) && (out_z_sz<in_z_sz))
   z_output_space  = max(floor((out_z_sz-in_z_sz)/2),0) + [1:min(in_z_sz,out_z_sz)];
   z_input_space   = max(floor((in_z_sz-out_z_sz)/2),0) + [1:min(in_z_sz,out_z_sz)];
else
   z_output_space  = max(ceil((out_z_sz-in_z_sz)/2),0) + [1:min(in_z_sz,out_z_sz)];
   z_input_space   = max(ceil((in_z_sz-out_z_sz)/2),0) + [1:min(in_z_sz,out_z_sz)];
end

% perform the up/down sampling
kSpaceLow  = zeros( out_x_sz, out_y_sz, out_z_sz, size(kSpaceHigh,4), size(kSpaceHigh,5) );
% kSpaceHigh = fftshift(fftshift(fftshift(kSpaceHigh,1),2),3); % already zero-centered
kSpaceLow(x_output_space,y_output_space, z_output_space, :, :) = kSpaceHigh(x_input_space,y_input_space, z_input_space, :, :);

end


text:
function errstring = checkstructfields(ds,varargin)

% CHECKSTRUCTFIELDS Checks if a data structure has the required fields.
%
%   errstring = checkstructfields(ds, field1, field2, ...)
%
%   This function checks if data structure 'ds' has the fields specified by the
%   character arrays (strings) 'field1', 'field2', ... e.g.
%
%        missing_fields = checkstructfields(ds, 'field1', 'field2', 'field3');
%
%   The names of the missing fields are returned in 'errstring'.
%
%   Copyright (c) Oregon Health & Science University (2006)
%
%   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for
%   academic use only (see included license file) and can be obtained from
%   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the
%   software should contact rebel@csee.ogi.edu for commercial licensing information.
%
%   See LICENSE (which should be part of the main toolkit distribution) for more
%   detail.

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

errstring = '';

%-- Check existence of required data structure fields
for j=1:length(varargin),
  if ~ischar(varargin{j})
    error(' [ checkstructfields ]   The field name arguments should be character arrays (strings).');
  end
  if ~isfield(ds,varargin{j}),
    errstring = [varargin{j} ' '];
  end
end

repo: ekfmonoslam


text:
% StackExchange Signal Processing Q59089
% https://dsp.stackexchange.com/questions/59089
% The Gradient of Least Squares of 2D Image Convolution
% References:
%   1.  A
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes
% - 1.0.000     24/06/2019
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = OFF;

SHAPE_MODE_FULL   = 'full';
SHAPE_MODE_VALID  = 'valid';

DIFF_MODE_FORWARD   = 1;
DIFF_MODE_BACKWARD  = 2;
DIFF_MODE_CENTRAL   = 3;
DIFF_MODE_COMPLEX   = 4;

%% Simulation Parameters

numRowsImage = 13;
numColsImage = 10;

numRowsKernel = 5;
numColsKernel = 3;

diffMode    = DIFF_MODE_COMPLEX;
epsVal      = 1e-6;

hCorr2D = @(mX, mH, convShape) conv2(mX, flip(flip(mH, 1), 2), convShape);

%% Generate Data

mX = randn(numRowsImage, numColsImage);
mH = randn(numRowsKernel, numColsKernel);
mY = randn(numRowsImage - numRowsKernel + 1, numColsImage - numColsKernel + 1);

hConvX = @(vX) 0.5 * sum((conv2(reshape(vX, numRowsImage, numColsImage), mH, SHAPE_MODE_VALID) - mY) .^ 2, 'all');
hConvH = @(vH) 0.5 * sum((conv2(mX, reshape(vH, numRowsKernel, numColsKernel), SHAPE_MODE_VALID) - mY) .^ 2, 'all');

%% Gradient of Convolution

% Gradient with respect to X
vGNumeric   = CalcFunGrad(mX(:), hConvX, diffMode, epsVal);
vGAnalytic  = conv2((conv2(mX, mH, SHAPE_MODE_VALID) - mY), mH(end:-1:1, end:-1:1), SHAPE_MODE_FULL);
% vGAnalytic = hCorr2D((conv2(mX, mH, SHAPE_MODE_VALID) - mY), mH, SHAPE_MODE_FULL);

max(abs(vGNumeric(:) - vGAnalytic(:)))

% Gradient with respect to H
vGNumeric   = CalcFunGrad(mH(:), hConvH, diffMode, epsVal);
vGAnalytic  = conv2(mX(end:-1:1, end:-1:1), (conv2(mX, mH, SHAPE_MODE_VALID) - mY), SHAPE_MODE_VALID);

max(abs(vGNumeric(:) - vGAnalytic(:)))

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function display_regions(bndinfo, regions, type, varargin)

sps = region_to_sp(bndinfo, regions, type, varargin{:});
mask = sps(bndinfo.wseg);

imagesc(mask);
axis image

author: Cloud-CV

repo: object-proposals


text:
function view = makeHemisphereButtons(view)
%
% view = makeHemisphereButtons(view)
%
% Installs push buttons that set curSlice.
%
% djh, 1/97

% Left button Callback:
%   selectButton(view.ui.sliceButtons,1);
%   view=refreshScreen(view);
% The hemisphere buttons pass a '2' in the recompute field of refreshView.
% This lets refreshView deal with the image rotation slider bar.
callbackStr = ...
    ['selectButton(',view.name,'.ui.sliceButtons,1); ',...
	view.name,'=refreshScreen(',view.name,',2);'];
view.ui.sliceButtons(1) = ...
    uicontrol('Style','radiobutton',...
    'String','Left',...
	'BackgroundColor', get(gcf, 'Color'), ...
    'Units','normalized',...
    'Position',[.1 .01 .1 .05],...
    'Callback', callbackStr);

% Right button Callback:
%   selectButton(view.ui.sliceButtons,1);
%   view=refreshScreen(view);
callbackStr = ...
    ['selectButton(',view.name,'.ui.sliceButtons,2); ',...
	view.name,'=refreshScreen(',view.name,',2);'];
view.ui.sliceButtons(2) = ...
    uicontrol('Style','radiobutton',...
    'String','Right',...
    'Units','normalized',...
	'BackgroundColor', get(gcf, 'Color'), ...
    'Position',[.25 .01 .1 .05],...
    'Callback', callbackStr);
return;


text:
function net = load_cnn(fparams, im_size)

	net = load(['networks/' fparams.nn_name]);
	net = vl_simplenn_tidy(net);
	net.layers = net.layers(1:max(fparams.output_layer));

	if strcmpi(fparams.input_size_mode, 'cnn_default')
	    base_input_sz = net.meta.normalization.imageSize(1:2);
	elseif strcmpi(fparams.input_size_mode, 'adaptive')
	    base_input_sz = im_size(1:2);
	else
	    error('Unknown input_size_mode');
	end

	net.meta.normalization.imageSize(1:2) = round(base_input_sz .* fparams.input_size_scale);
	net.meta.normalization.averageImageOrig = net.meta.normalization.averageImage;

	if isfield(net.meta,'inputSize')
	    net.meta.inputSize = base_input_sz;
	end

	if size(net.meta.normalization.averageImage,1) > 1 || size(net.meta.normalization.averageImage,2) > 1
	    net.meta.normalization.averageImage = imresize(single(net.meta.normalization.averageImage), net.meta.normalization.imageSize(1:2));
	end

	net.info = vl_simplenn_display(net);
end

repo: CFWCR


text:
% LIONSIMBA example script
% Multiple cells scenario: simulates three cells in serial connection.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

% Clear the workspace
clear

% Define the integration times.
t0 = 0;
tf = 10^4;
% Define the parameters structure.
param{1} = Parameters_init;
param{2} = Parameters_init;
param{3} = Parameters_init;

% Modify the initial SOC of cell #1
param{1}.cs_n_init = 0.95*param{1}.cs_n_init;

% Double the length of the positive electrode of cell #2
param{2}.len_p=2*param{2}.len_p;

% Simulate the battery of 3 cells
results = startSimulation(t0,tf,[],-30,param);

ncells = length(param);

v_tot = 0;
for i=1:ncells
    v_tot = v_tot + results.Phis{i}(:,1)-results.Phis{i}(:,end);
end
%% Plot the results

figure(1)
subplot(2,ncells,(1:ncells))
plot(results.time{1},v_tot,'LineWidth',6)
grid on
box on
xlabel('Time [s]')
ylabel('Pack Voltage [V]')
hold on
xlim([0 results.time{1}(end)])
ylim([2.5*length(param) 4.2*length(param)])

figure_index = ncells+1;
for i=1:ncells
    subplot(2,ncells,figure_index)
    plot(results.time{i},results.Voltage{i},'--','LineWidth',3)
    grid on
    box on
    xlabel('Time [s]')
    ylabel(['Cell #',num2str(i),' Voltage [V]'])
    figure_index = figure_index+1;
    ylim([2.5 4.2])
    xlim([0 results.time{1}(end)])
end

n_plots = 1;
figure_index = 1;
figure(2)
for i=1:ncells
    subplot(ncells,n_plots,figure_index)
    ce_indices = [1 results.parameters{i}.Np+1 results.parameters{i}.Np+results.parameters{i}.Ns results.parameters{i}.Np+results.parameters{i}.Ns+results.parameters{i}.Nn];
    plot(results.time{i},results.ce{i}(:,ce_indices),'LineWidth',3)
    grid on
    box on
    if(figure_index==ncells)
        xlabel('Time [s]')
    else
%         set(gca,'xtick',[])
        set(gca,'xticklabel',[])
    end
    ylabel(['Cell #',num2str(i),' - c_{e} [mol/m^3}'])
    xlim([0 results.time{1}(end)])
    ylim([0 2200])
    figure_index = figure_index+1;
end

n_plots = 1;
figure_index = 1;
figure(3)
for i=1:ncells
    subplot(ncells,n_plots,figure_index)
    phie_indices = [1 results.parameters{i}.Np+1 results.parameters{i}.Np+results.parameters{i}.Ns results.parameters{i}.Np+results.parameters{i}.Ns+results.parameters{i}.Nn];
    plot(results.time{i},results.Phie{i}(:,phie_indices),'LineWidth',3)
    grid on
    box on
    if(figure_index==ncells)
        xlabel('Time [s]')
    else
%         set(gca,'xtick',[])
        set(gca,'xticklabel',[])
    end
    ylabel(['Cell #',num2str(i)])
    xlim([0 results.time{1}(end)])
    ylim([-0.3 0])
    figure_index = figure_index+1;
end

n_plots = 1;
figure_index = 1;
figure(4)
for i=1:ncells
    subplot(ncells,n_plots,figure_index)
    cs_indices = [1 results.parameters{i}.Np results.parameters{i}.Np+1 results.parameters{i}.Np+results.parameters{i}.Nn];
    plot(results.time{i},results.cs_surface{i}(:,cs_indices),'LineWidth',3)
    grid on
    box on
    if(figure_index==ncells)
        xlabel('Time [s]')
    else
%         set(gca,'xtick',[])
        set(gca,'xticklabel',[])
    end
    ylabel(['Cell #',num2str(i)])
    xlim([0 results.time{1}(end)])
    ylim([100 51554])
    figure_index = figure_index+1;
end

n_plots = 1;
figure_index = 1;
figure(5)
for i=1:ncells
    subplot(ncells,n_plots,figure_index)
    cs_indices = [1 results.parameters{i}.Np results.parameters{i}.Np+1 results.parameters{i}.Np+results.parameters{i}.Nn];
    plot(results.time{i},results.Temperature{i}(:,end),'LineWidth',3)
    grid on
    box on
    if(figure_index==ncells)
        xlabel('Time [s]')
    else
%         set(gca,'xtick',[])
        set(gca,'xticklabel',[])
    end
    ylabel(['Cell #',num2str(i),' - Temp [K]'])
    xlim([0 results.time{1}(end)])
    ylim([296 315])
    figure_index = figure_index+1;
end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function [p,q,r] = btf (A)                                                  %#ok
%BTF permute a square sparse matrix into upper block triangular form
% with a zero-free diagonal, or with a maximum number of nonzeros along the
% diagonal if a zero-free permutation does not exist.
%
% Example:
%       [p,q,r] = btf (A) ;
%       [p,q,r] = btf (A,maxwork) ;
%
% If the matrix has structural full rank, this is essentially identical to
%
%       [p,q,r] = dmperm (A)
%
% except that p, q, and r will differ in trivial ways.  Both return an upper
% block triangular form with a zero-free diagonal, if the matrix is
% structurally non-singular.  The number and sizes of the blocks will be
% identical, but the order of the blocks, and the ordering within the blocks,
% can be different.
%
% If the matrix is structurally singular, the q from btf will contain negative
% entries.  The permuted matrix is C = A(p,abs(q)), and find(q<0) gives a list
% of indices of the diagonal of C that are equal to zero.  This differs from
% dmperm, which does not place the maximum matching along the main diagonal
% of C=A(p,q), but places it above the diagonal instead.
%
% The second input limits the maximum amount of work the function does to
% be maxwork*nnz(A), or no limit at all if maxwork <= 0.  If the function
% terminates early as a result, a maximum matching may not be found, and the
% diagonal of A(p,abs(q)) might not have the maximum number of nonzeros
% possible.  Also, the number of blocks (length(r)-1) may be larger than
% what btf(A) or dmperm(A) would compute.
%
% See also maxtrans, strongcomp, dmperm, sprank

% Copyright 2004-2007, Tim Davis, University of Florida
% with support from Sandia National Laboratories.  All Rights Reserved.

error ('btf mexFunction not found') ;