text:
function shepard_interp_2d_test01 ( prob, g, p )

%*****************************************************************************80
%
%% SHEPARD_INTERP_2D_TEST01 tests SHEPARD_INTERP_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the problem number.
%
%    Input, integer G, the grid number.
%
%    Input, real P, the power used in the distance weighting.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SHEPARD_INTERP_2D_TEST01:\n' );
  fprintf ( 1, '  Interpolate data from TEST_INTERP_2D problem #%d\n', prob );
  fprintf ( 1, '  using grid #%d\n', g );
  fprintf ( 1, '  using Shepard interpolation with P = %g\n', p );

  nd = g00_size ( g );
  fprintf ( 1, '  Number of data points = %d\n', nd );

  [ xd, yd ] = g00_xy ( g, nd );

  zd = f00_f0 ( prob, nd, xd, yd );

  if ( p == 1.0 )
    r8vec3_print ( nd, xd, yd, zd, '  X, Y, Z data:' );
  end
%
%  #1:  Does interpolant match function at interpolation points?
%
  ni = nd;
  xi = xd;
  yi = yd;
  zi = shepard_interp_2d ( nd, xd, yd, zd, p, ni, xi, yi );

  int_error = norm ( zi - zd ) / ni;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  L2 interpolation error averaged per interpolant node = %g\n', int_error );
%
%  #2: Plot the linear interpolant.
%
  if ( p == 1.0 )
    figure ( 1 );
    clf
    hold on
    tri = delaunay ( xd, yd );
    h = trisurf ( tri, xd, yd, zd );
    ft = f00_title ( prob );
    title ( ft );
    grid on
    xlabel ( '<--- X --->' );
    ylabel ( '<--- Y --->' );
    zlabel ( '<---Z(X,Y)--->' );
    l = light ( 'Position', [ 0, 0, 25 ] );
    set ( gca, 'CameraPosition', [ -10, -10, 20 ] );
    lighting phong
    shading interp
    hold off
    filename = sprintf ( 'p%02d_data.png', prob );
    print ( '-dpng', filename );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Created plot file "%s".\n', filename );
  end
%
%  #3: Plot the Shepard interpolant.
%
  figure ( 2 );
  clf
  nix = 21;
  niy = 21;
  xmin = 0.0;
  xmax = 1.0;
  ymin = 0.0;
  ymax = 1.0;
  xi = linspace ( xmin, xmax, nix );
  yi = linspace ( ymin, ymax, niy );
  [XI,YI] = meshgrid ( xi, yi );
  NI = nix * niy;
  XI = reshape ( XI, NI, 1 );
  YI = reshape ( YI, NI, 1 );
  ZI = shepard_interp_2d ( nd, xd, yd, zd, p, NI, XI, YI );
  XI = reshape ( XI, nix, niy );
  YI = reshape ( YI, nix, niy );
  ZI = reshape ( ZI, nix, niy );
  surf ( XI, YI, ZI );
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y --->' );
  zlabel ( '<---Z(X,Y)--->' );
  title ( sprintf ( 'Shepard interpolant for problem %d with P = %g', prob, p ) )
  grid on
  hold off
  filename = sprintf ( 'p%02d_power%g.png', prob, p );
  print ( '-dpng', filename );
  fprintf ( 1, '  Created plot file "%s".\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function result=RS_PF_Rationality(m,forecasterror,x,pdate)
%%   INPUT:
%%   m = size of rolling window
%%	forecastserror = actual - forecast
%%	x is the variable you want to test coefficients equal to zero
%%   pdate is the vector of dates

T=length(forecasterror);
if m>=T
    err_dlg = errordlg('The window size of the Rossi-Sekhposyan (JAE,2016) test cannot be equal to or larger than the forecasted sample.');
    waitfor(err_dlg);
elseif m==0
    err_dlg = errordlg('The window size of the Rossi-Sekhposyan (JAE,2016) test cannot be equal to zero.');
    waitfor(err_dlg);
end

cvtable = [
      NaN     0.1000    0.2000    0.3000    0.4000    0.5000    0.6000    0.7000    0.8000    0.9000
    1.0000   10.5066    9.0503    8.0245    7.1035    6.3957    5.6112    5.1113    4.6141    3.9748
    2.0000   21.2392   18.0544   15.8290   13.9122   13.0720   11.1526   10.4549    9.0570    7.8723
    3.0000   31.4497   26.8866   23.7832   21.4577   19.6097   17.4180   15.3225   13.5010   11.4381
    4.0000   43.5150   36.9028   32.8187   28.4075   25.1774   23.3645   20.5785   17.6700   15.5384
    5.0000   52.4148   45.7998   39.6896   35.7848   32.0200   28.4850   26.2204   23.1738   19.1090
    6.0000   62.6771   54.3749   47.4711   42.4503   38.4920   34.9394   30.4063   27.9807   23.8787
    7.0000   74.8406   62.3659   56.2449   49.0721   44.4213   39.6189   36.4280   33.0852   26.9654
    8.0000   84.5728   72.8813   63.2267   56.8973   51.5069   45.7856   41.3975   36.7853   31.2008
    9.0000  109.6177   95.3818   87.2701   77.9691   72.5004   63.5162   60.0533   51.9429   47.0975
   10.0000  122.4825  107.9759   94.2844   88.3139   80.2026   71.6698   67.4220   58.3127   54.7817 ];

P     = size(forecasterror,1);
mu    = m/P;
cvcol = round(mu*10)+1;

if cvcol==1;
    cvcol = 2;
end;

if cvcol>10;
    cvcol = 10;
end;

nreg = size(x,2);
cv   = cvtable(nreg+1,cvcol);

resultt=[];
for t=m:P
    % Calculate OLS Wald-Test
    out     = bear.RS_PF_OLS_Wald(forecasterror(t-m+1:t,:),x(t-m+1:t,:));
    resultt = [resultt;out];
end;

ptruncdate = pdate(m:end,1);

result.MZ         = max(resultt);
result.cv         = cv;
result.ptruncdate = ptruncdate;
result.resultt    = resultt(:,1);
result.cvones     = cv*ones(size(resultt,1),1);

author: european-central-bank

repo: BEAR-toolbox


text:
function [dat] = read_plexon_ddt(filename, begsample, endsample)

% READ_PLEXON_DDT reads header or data from a Plexon *.ddt file,
% which is a Plexon continuous data file optimized for continuous
% (streaming) recording where every channel is continuously recorded
% without gaps and the recording includes any dead time between spikes.
%
% Use as
%   [hdr] = read_plexon_ddt(filename)
%   [dat] = read_plexon_ddt(filename, begsample, endsample)
%
% samples start counting at 1
% returned values are in mV

% Copyright (C) 2005-2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Version 100: Samples are assumed to be 12 bits.  All channels have the same NIDAQ gain, and preamp
% gain is not saved.
%     int     Version;         // =100
%     int     DataOffset;      // Offset into the file where the data starts
%     double  Freq;            // Digitization frequency
%     int     NChannels;       // Number of channels
%     int     Year;            // Time/date when the data was acquired
%     int     Month;
%     int     Day;
%     int     Hour;
%     int     Minute;
%     int     Second;
%     int     Gain;            // NIDAQ gain (the same gain for all channels
%     char    Comment[128];    // User-supplied comment
%     unsigned char Padding[256];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dat = [];
fid = fopen_or_error(filename, 'rb', 'ieee-le');

Version    = fread(fid, 1, 'int'    );
dat.Version    = Version;

% this is common for all subformats
dat.DataOffset = fread(fid, 1, 'int'    );
dat.Freq       = fread(fid, 1, 'double' );
dat.NChannels  = fread(fid, 1, 'int'    );
dat.Year       = fread(fid, 1, 'int'    );
dat.Month      = fread(fid, 1, 'int'    );
dat.Day        = fread(fid, 1, 'int'    );
dat.Hour       = fread(fid, 1, 'int'    );
dat.Minute     = fread(fid, 1, 'int'    );
dat.Second     = fread(fid, 1, 'int'    );
dat.Gain       = fread(fid, 1, 'int'    );
dat.Comment        = char(fread(fid, 128, 'char')');

if Version==100
  % the remainder of the header contains bytes for padding
  dat.Padding    = fread(fid, 256, 'char' );
elseif Version==101
  % Version 101: A field was added to indicate bits-per-sample (12 or 16).  All channels have the same NIDAQ
  % gain, and preamp gain is not saved.
  dat.BitsPerSample = fread(fid,   1, 'char');
  dat.Padding       = fread(fid, 255, 'char');
elseif Version==102
  % Version 102: A byte array of per-channel NIDAQ gains was added.  The Gain field now indicates the
  % preamp gain,   a value of 1 usually indicating that no preamp gain was specified in the application that
  % recorded the DDT.
  dat.BitsPerSample = fread(fid,   1, 'char');
  dat.ChannelGain   = fread(fid,  64, 'char');
  dat.Padding       = fread(fid, 191, 'char');
elseif Version==103
  % Version 103: A field was added to specify ADC maximal input voltage. The value is integer number in
  % millivolts.
  dat.BitsPerSample  = fread(fid,   1, 'char');
  dat.ChannelGain    = fread(fid,  64, 'char');
  dat.MaxMagnitudeMV = fread(fid,   1, 'short');
  dat.Padding        = fread(fid, 189, 'char');
else
  ft_error('unsupported version of ddt file');
end

% determine the number of samples by looking at the length of the datafile
offset = ftell(fid);
fseek(fid, 0, 'eof');
dat.NSamples = (ftell(fid) - offset) / (2 * dat.NChannels);

if nargin==1
  % do not read any data, only return the header
  fclose(fid);
  return
end

% A DDT file contains a file header followed by an array of A/D samples stored as 16-bit integers regardless
% of whether values are 12 or 16-bits. The samples are multiplexed.

% read the desired samples of the data
offset = dat.DataOffset + (begsample-1)*dat.NChannels;
nsampl = endsample - begsample + 1;
fseek(fid, offset, 'bof');
dat.data = fread(fid, [dat.NChannels nsampl], 'int16');

fclose(fid);

% apply the calibration to obtain the signal in physical units
if Version==100
  PreAmpGain = 1000;
  dat.data = dat.data * 5000 / (2048 * dat.Gain * PreAmpGain);
elseif Version==101
  PreAmpGain = 1000;
  dat.data = dat.data * 5000 / (0.5 * (2^dat.BitsPerSample) * dat.Gain * PreAmpGain);
elseif Version==102
  for i=1:dat.NChannels
    dat.data(i,:) = dat.data(i,:) * 5000 ./ (0.5 * (2^dat.BitsPerSample) * dat.ChannelGain(i));
  end
elseif Version==103
  % I am not sure whether the calibration like this is correct, since the
  % Plexon documentation does not explicitely specify how to do it for the
  % 104 file format. I presume that it is identical to the 103 format.
  for i=1:dat.NChannels
    dat.data(i,:) = dat.data(i,:) * 5000 ./ (0.5 * (2^dat.BitsPerSample) * dat.ChannelGain(i));
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


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

%*****************************************************************************80
%
%% HYPERBALL01_SAMPLE uniformly samples the unit hyperball in M dimensions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Russell Cheng,
%    Random Variate Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley, 1998, pages 168.
%
%    Reuven Rubinstein,
%    Monte Carlo Optimization, Simulation, and Sensitivity
%    of Queueing Networks,
%    Krieger, 1992,
%    ISBN: 0894647644,
%    LC: QA298.R79.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, real X(M,N), the points.
%
  x = randn ( m, n );
  norm = ones ( 1, m ) * ( x.^2 );
  norm = sqrt ( norm );
  for i = 1 : m
    x(i,1:n) = x(i,1:n) ./ norm(1:n);
  end

  for j = 1 : n
    r = rand ( 1, 1 );
    x(1:m,j) = r ^ ( 1.0 / m ) * x(1:m,j);
  end

  return
end


text:
function model = treeTrain(X, Y, opts)
% Train a random tree
% X is NxD, each D-dimensional row is a data point
% Y is Nx1 discrete labels of classes
% returned model is to be directly plugged into treeTest

d= 5; % max depth of the tree

if nargin < 3, opts= struct; end
if isfield(opts, 'depth'), d= opts.depth; end

u= unique(Y);
[N, D]= size(X);
nd= 2^d - 1;
numInternals = (nd+1)/2 - 1;
numLeafs= (nd+1)/2;

weakModels= cell(1, numInternals);
% if we can afford to store as non-sparse (100MB array, say), it is
% slightly faster.
if storage([N nd]) < 100
    dataix= zeros(N, nd); % boolean indicator of data at each node
else
    dataix= sparse(N, nd);
end

leafdist= zeros(numLeafs, length(u)); % leaf distribution

% Propagate data down the tree while training weak classifiers at each node
for n = 1: numInternals

    % get relevant data at this node
    if n==1
        reld = ones(N, 1)==1;
        Xrel= X;
        Yrel= Y;
    else
        reld = dataix(:, n)==1;
        Xrel = X(reld, :);
        Yrel = Y(reld);
    end

    % train weak model
    weakModels{n}= weakTrain(Xrel, Yrel, opts);

    % split data to child nodes
    yhat= weakTest(weakModels{n}, Xrel, opts);

    dataix(reld, 2*n)= yhat;
    dataix(reld, 2*n+1)= 1 - yhat; % since yhat is in {0,1} and double
end

% Go over leaf nodes and assign class statistics
for n= (nd+1)/2 : nd
    reld= dataix(:, n);
    hc = histc(Y(reld==1), u);
    hc = hc + 1; % Dirichlet prior
    leafdist(n - (nd+1)/2 + 1, :)= hc / sum(hc);
end

model.leafdist= leafdist;
model.depth= d;
model.classes= u;
model.weakModels= weakModels;
end

author: karpathy

repo: Random-Forest-Matlab


text:
function my_imshow(img_data)

try
    imshow(img_data);
catch err_info
    disp('image cannot display!');
end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [catdim] = checkinput(x)

if ~iscell(x)
  error('input should be a cell-array');
end

siz    = size2(x, [], 'cell');
ndim   = size(siz, 2);

if ndim>2
  error('number of dimensions within a cell > 2 is currently not supported');
end

if length(x)<numel(x)
  error('currently only a 1xN or Nx1 cell-array is supported in the input');
end

same   = false(1,ndim);
for k = 1:ndim
  same(k) = all(siz(1,k)==siz(:,k));
end

if sum(same)==1
  % data need to be concatenated across the other dimension
  catdim = ~same;
elseif sum(same)==0
  % data cannot be concatenated across any of the dimensions
  error('the cells in the input array cannot be concatenated across any of the dimensions in the cells');
elseif sum(same)>1
  % more than one possibility
  catdim = same;
end


text:
function [fx,dF_dX,dF_dTheta] = f_lin2D(Xt,Theta,ut,inF)
% dummy 2D linear evolution function

deltat = inF.deltat;

if ~ isfield (inF, 'a')
    inF.a = 1;
end
if ~ isfield (inF, 'b')
    inF.b = 1e-1;
end

inF.a = inF.a * exp (Theta(1));

A = [- inF.b, - inF.a;
     1, - inF.b];

fx = Xt + deltat * (A * Xt + ut);
dF_dX = eye (size (Xt, 1)) + deltat * A';

dF_dTheta = deltat * [- inF.a * Xt(2), 0];

author: MBB-team

repo: VBA-toolbox


text:
function [msi] = read_bti_m4d(filename);

% READ_BTI_M4D
%
% Use as
%   msi = read_bti_m4d(filename)

% Copyright (C) 2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: read_bti_m4d.m 945 2010-04-21 17:41:20Z roboos $

[p, f, x] = fileparts(filename);
if ~strcmp(x, '.m4d')
  % add the extension of the header
  filename = [filename '.m4d'];
end

fid = fopen(filename, 'r');
if fid==-1
  error(sprintf('could not open file %s', filename));
end

% start with an empty header structure
msi = struct;

% these header elements contain strings and should be converted in a cell-array
strlist = {
  'MSI.ChannelOrder'
  };

% these header elements contain numbers and should be converted in a numeric array
%   'MSI.ChannelScale'
%   'MSI.ChannelGain'
%   'MSI.FileType'
%   'MSI.TotalChannels'
%   'MSI.TotalEpochs'
%   'MSI.SamplePeriod'
%   'MSI.SampleFrequency'
%   'MSI.FirstLatency'
%   'MSI.SlicesPerEpoch'
% the conversion to numeric arrays is implemented in a general fashion
% and all the fields above are automatically converted
numlist = {};

line = '';

msi.grad.label = {};
msi.grad.pnt   = zeros(0,3);
msi.grad.ori   = zeros(0,3);
while ischar(line)
  line = cleanline(fgetl(fid));
  if isempty(line) || (length(line)==1 && all(line==-1))
    continue
  end

  sep = strfind(line, ':');
  if length(sep)==1
    key = line(1:(sep-1));
    val = line((sep+1):end);
  elseif length(sep)>1
    % assume that the first separator is the relevant one, and that the
    % next ones are part of the value string (e.g. a channel with a ':' in
    % its name
    sep = sep(1);
    key = line(1:(sep-1));
    val = line((sep+1):end);
  elseif length(sep)<1
    % this is not what I would expect
    error('unexpected content in m4d file');
  end

  if ~isempty(strfind(line, 'Begin'))
    sep = strfind(key, '.');
    sep = sep(end);
    key = key(1:(sep-1));

    % if the key ends with begin and there is no value, then there is a block
    % of numbers following that relates to the magnetometer/gradiometer information.
    % All lines in that Begin-End block should be treated seperately
    val = {};
    lab = {};
    num = {};
    ind = 0;
    while isempty(strfind(line, 'End'))
      line = cleanline(fgetl(fid));
      if isempty(line) || (length(line)==1 && all(line==-1)) || ~isempty(strfind(line, 'End'))
        continue
      end
      ind = ind+1;
      % remember the line itself, and also cut it into pieces
      val{ind} = line;
      % the line is tab-separated and looks like this
      % A68 0.0873437   -0.075789   0.0891512   0.471135    -0.815532   0.336098
      sep = find(line==9); % the ascii value of a tab is 9
      sep = sep(1);
      lab{ind} = line(1:(sep-1));
      num{ind} = str2num(line((sep+1):end));

    end % parsing Begin-End block
    val = val(:);
    lab = lab(:);
    num = num(:);
    num = cell2mat(num);
    % the following is FieldTrip specific
    if size(num,2)==6
      msi.grad.label = [msi.grad.label; lab(:)];
      % the numbers represent position and orientation of each magnetometer coil
      msi.grad.pnt   = [msi.grad.pnt; num(:,1:3)];
      msi.grad.ori   = [msi.grad.ori; num(:,4:6)];
    else
      error('unknown gradiometer design')
    end
  end

  % the key looks like 'MSI.fieldname.subfieldname'
  fieldname = key(5:end);

  % remove spaces from the begin and end of the string
  val = strtrim(val);

  % try to convert the value string into something more usefull
  if ~iscell(val)
    % the value can contain a variety of elements, only some of which are decoded here
    if ~isempty(strfind(key, 'Index')) || ~isempty(strfind(key, 'Count')) || any(strcmp(key, numlist))
      % this contains a single number or a comma-separated list of numbers
      val = str2num(val);
    elseif ~isempty(strfind(key, 'Names')) || any(strcmp(key, strlist))
      % this contains a comma-separated list of strings
      val = tokenize(val, ',');
    else
      tmp = str2num(val);
      if ~isempty(tmp)
        val = tmp;
      end
    end
  end

  % assign this header element to the structure
  msi = setsubfield(msi, fieldname, val);

end % while ischar(line)

fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION to remove spaces from the begin and end
% and to remove comments from the lines
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function line = cleanline(line)
if isempty(line) || (length(line)==1 && all(line==-1))
  return
end
comment = findstr(line, '//');
if ~isempty(comment)
  line(min(comment):end) = ' ';
end
line = strtrim(line);

author: PatternRecognition

repo: OpenBMI


text:
function set = exponential( varargin )

%EXPONENTIAL   Exponential cone.
%   EXPONENTIAL, called with no arguments, creates three scalar variables X,
%   Y, and Z and constraints them to lie in an exponetial cone. That is,
%   given the declaration
%       variables x y z
%       {x,y,z} == exp_cone
%   constraints the variables to satisfy
%       y*exp(x/y) <= z
%       y > 0
%   The inequality form does not obey the disciplined convex programming
%   ruleset, but a function EXP_P has been created to represent this
%   computation; so the set declaration above is equivalent to
%       EXP_P(X,Y) <= Z
%   EXP_CONE(SX), where SX is a size vector, creates three array variables
%   X, Y, and Z, each of size SX, which are constrained elementwise to
%   satisfy EXP_P(X,Y) <= Z. If SX is empty, then SX=[1,1] is assumed.

cvx_expert_check( 'exponential' );
[ sx, gp ] = cvx_get_dimlist( varargin ); %#ok
gp = ~isempty(gp) && gp;
if gp,
    cvx_begin gp set
        variables ex( sx ) ey( sx ) ez( sx )
        x = log( ex ); y = log( ey ); z = log( ez );
else
    cvx_begin set
        variable x( sx )
        variable y( sx ) nonnegative_
        variable z( sx ) nonnegative_
end
    cvx_pushcone( true, 'exponential', [ vec(x)' ; vec(y)' ; vec(z)' ] );
cvx_end
if gp,
    set = cvxtuple( struct( 'x', ex, 'y', ey, 'z', ez ) );
else
    set = cvxtuple( struct( 'x', x, 'y', y, 'z', z ) );
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% Combined wire sizing and spacing
% Section 5.5,  L. Vandenberghe, S. Boyd, and A. El Gamal
% "Optimizing dominant time constant in RC circuits"
% Original by Lieven Vandenberghe
% Adapted for CVX by Joelle Skaf - 11/27/05
% Modified by Michael Grant - 3/8/06
%
% The problem is to determine the optimal sizes of interconnect wires and
% the optimal distances between them. We will consider an example with 3
% wires, each consisting of 5 segments (see paper, fig.21). The variables
% are the widths wij , and the distances s1 and s2 between the wires.
% The difference with the models used in other scripts is that we include a
% parasitic capacitance between the wires.
% The objective is to minimize the total width s1+s2.
% The problem can be formulated with the following SDP:
%               mimimize    s1 + s2
%                    s.t.   Tmax*G(w11,..,w35)-C(w11,..,w35,t11,..,t23) >=0
%                           1/t1j <= s1 - w1j - 0.5*w2j     , j = 1,..,5
%                           1/t2j <= s2 - w3j - 0.5*w2j     , j = 1,..,5
%                           0 <= tij <= 1         , i = 1,2 , j = 1,..,5
%                           t1 >=0, t2 >= 0, s1 >=smin, s2>=smin
%                           0 <= wij <= wmax
% the 2nd and 3rd constraints are nonlinear convex constraints that can be
% cast as 3 x 3-LMIs (Please refer to the paper for more details).

%
% Circuit parameters
%

n = 6;           % number of nodes per wire
N = 3*n;         % total number of nodes
m = n-1;         % number of segments per wire
alpha = 1;       % conductance per segment is is alpha*size
beta = 0.5;      % capacitance per segment is twice beta*size
gamma = 2;       % coupling capacitance is twice gamma*distance
G0 = 100;        % source conductance
C0 = [10,20,30]; % loads of first, second, third wires
wmin = 0.1;      % minimum width
wmax = 2.0;      % maximum width
smin = 1.0;      % minimum distance between wires
smax = 50;       % upper bound on s1 and s2  (meant to be inactive)

%
% Construct the capacitance and conductance matrices
%   C(x) = C0 + w11 * C1 + w21 * C2 + ...
%   G(x) = G0 + w11 * G1 + w21 * G2 + ...
% and we assemble the coefficient matrices together as follows:
%   CC = [ C0(:) C1(:) C2(:) ... ]
%   GG = [ G0(:) G1(:) G2(:) ... ]
%

CC = zeros(N,N,5*m+1);
GG = zeros(N,N,3*m+1);
for w = 0 : 2,
    % Constant terms
    CC(w*n+n,w*n+n,1) = C0(w+1);
    GG(w*n+1,w*n+1,1) = G0;
    for i = 1 : m,
        % capacitances to ground
        CC(w*n+[i,i+1],w*n+[i,i+1],w*m+i+1) = beta*[1,0;0,1];
        if w < 2,
            % coupling capacitors
            CC(w*n+[i,  n+i  ],w*n+[i,  n+i  ],(w+3)*m+i+1) = gamma*[1,-1;-1,1];
            CC(w*n+[i+1,n+i+1],w*n+[i+1,n+i+1],(w+3)*m+i+1) = gamma*[1,-1;-1,1];
        end
        % segment conductances
        GG(w*n+[i,i+1],w*n+[i,i+1],w*m+i+1) = alpha*[1,-1;-1,1];
    end
end
% Reshape for Matlab use
CC = reshape(CC,N*N,5*m+1);
GG = reshape(GG,N*N,3*m+1);

%
% Compute points the tradeoff curve and the two desired points
%

npts    = 50;
delays  = linspace( 85, 200, npts );
xdelays = [ 130, 90 ];
xnpts   = length(xdelays);
areas   = zeros(1,npts);
xareas  = zeros(1,xnpts);
for j = 1 : npts + xnpts,

    if j > npts,
        xj = j - npts;
        delay = xdelays(xj);
        disp( sprintf( 'Particular solution %d of %d (Tmax = %g)', xj, xnpts, delay ) );
    else,
        delay = delays(j);
        disp( sprintf( 'Point %d of %d on the tradeoff curve (Tmax = %g)', j, npts, delay ) );
    end

    %
    % Construct and solve the convex model
    %

    cvx_begin sdp quiet
        variables w(m,3) t(m,2) s(1,2)
        variable G(N,N) symmetric
        variable C(N,N) symmetric
        minimize( sum(s) )
        subject to
            G == reshape( GG * [ 1 ; w(:) ], N, N );
            C == reshape( CC * [ 1 ; w(:) ; t(:) ], N, N );
            delay * G - C >= 0;
            wmin <= w(:) <= wmax;
            t( : ) <= 1 / smin;
            s( : ) <= smax;
            inv_pos( t(:,1) ) <= s(1) - w(:,1) - 0.5 * w(:,2);
            inv_pos( t(:,2) ) <= s(2) - w(:,3) - 0.5 * w(:,2);
    cvx_end
    ss = cvx_optval;

    if j <= npts,
        areas(j) = ss;
    else,
        xareas(xj) = ss;

        %
        % Draw the wires
        %

        figure(4*xj-2);
        m2 = 2 * m;
        x1 = reshape( [ 1 : m ; 1 : m ], 1, m2 );
        x2 = x1( 1, end : -1 : 1 );
        y  = [ ss*ones(m2,1), s(2) + 0.5*w(x1,2), zeros(m2,1) ; ...
               ss-w(x2,1),    s(2) - 0.5*w(x2,2), w(x2,3)     ; ...
               ss,            s(2) + 0.5*w(1,2),  0           ];
        x1 = reshape( [ 0 : m - 1 ; 1 : m ], m2, 1 );
        x2 = x1( end : -1 : 1, 1 );
        x  = [ x1 ; x2 ; 0 ];
        hold off;
        fill( x, y, 0.9 * ones(size(y)) );
        hold on
        plot( x, y, '-' );
        axis( [-0.1, m+0.1,-0.1, ss+0.1]);
        colormap(gray);
        caxis([-1,1])
        title(sprintf('Solution (%d), Tmax = %g',xj,delay));

        %
        % Build the state space models and plot step responses
        %

        A = -inv(C)*G;
        T = linspace(0,2*delay,1000);
        B = -A * kron( eye(3), ones(n,1) );
        for inp = 1 : 3,
            figure(4*xj-2+inp);
            Y1 = simple_step(A,B(:,inp),T(2),length(T));
            hold off;
            plot(T,Y1([n,2*n,3*n],:),'-');
            hold on;
            text(T(1000),Y1(  n,1000),'v1');
            text(T(1000),Y1(2*n,1000),'v2');
            text(T(1000),Y1(3*n,1000),'v3');
            axis([0 2*delay -0.1 1.1]);
            % show dominant time constant
            plot(delay*[1;1], [-0.1;1.1], '--');
            title(sprintf('Solution (%d), Tmax = %g, step applied to wire %d',xj,delay,inp));
        end

    end

end

%
% Plot the tradeoff curve
%

figure(1);
ind = isfinite(areas);
plot(areas(ind), delays(ind));
xlabel('total width s_1 + s_2');
ylabel('dominant time constant');
title('Width-delay tradeoff curve')
hold on;
for k = 1 : xnpts,
    text( xareas(k), xdelays(k), sprintf( '(%d)', k ) );
end

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function varargout = wtps(varargin)
% VL_WTPS  Thin-plate spline warping
%   [XP1,XP2]=VL_WTPS(PHI,YP) computes the thin-plate spline (TPS)
%   specified by the basis PHI and the warped control point Yp.
%
%   Yp is a 2xK matrix with one column per control point and the basis
%   PHI is calculated by means of the VL_TPS function.
%
%   The thin-palte spline is defined on a domain X1,X2 and specified
%   by a set of points Y and their warp YP. The spline passes
%   interpolates exaclty the control points.
%
%   The parameters X1,X2 and Y are used to compute the basis PHI. This
%   operation is fairily slow, but computing the spline for a given Yp
%   is then very quick, as the operation is just a linear combination
%   of the basis.
%
%   Example::
%     To calculate the warped grid [X1,X2] by moving the control points Y to
%     the control points YP use:
%       [xp1,xp2]=VL_WTPS(VL_TPS(x1,x2,Y),Yp).
%
%   See also: VL_TPS(), VL_HELP().
[varargout{1:nargout}] = vl_wtps(varargin{:});

author: yihui-he

repo: panorama

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


text:
function rms_x=rms(x)

%Calculates RMS for x, NaN's are ignored
x=x(~isnan(x)); %Removing nans
rms_x=sqrt( sum(x.^2)./numel(x) );

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function c = correlation_damped_cosine ( n, rho, rho0 )

%*****************************************************************************80
%
%% CORRELATION_DAMPED_COSINE evaluates the damped cosine correlation function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Petter Abrahamsen,
%    A Review of Gaussian Random Fields and Correlation Functions,
%    Norwegian Computing Center, 1997.
%
%  Parameters:
%
%    Input, integer N, the number of arguments.
%
%    Input, real RHO(N,1), the arguments.
%
%    Input, real RHO0, the correlation length.
%
%    Output, real C(N,1), the correlations.
%
  rho = rho ( : );

  rhohat = abs ( rho ) / rho0;

  c = exp ( - rhohat ) .* cos ( rhohat );

  return
end


text:
%% OCR of hand-written digits using HoG and SVM
%
% In this tutorial, we will build an SVM classifer to recognize
% hand-written digits (0 to 9), using Histogram of Oriented Gradients (HOG)
% as feature vectors.
%
% For additional resources, see this
% <https://www.youtube.com/watch?v=0Zib1YEE4LU HoG video lecture>.
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/dd/d3b/tutorial_py_svm_opencv.html>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/digits.py>
% * <https://www.mathworks.com/help/vision/examples/digit-classification-using-hog-features.html>
%

function svm_hog_ocr_digits_demo()
    %% Data
    % Load MNIST handwritten digits (5000 images)
    % (params in this demo are hardcoded assuming 20x20 images)
    [imgs, labels, img] = load_mnist_digits();
    assert(numel(imgs)==5000 && isequal(size(imgs{1}),[20 20]));

    %%
    % Show a portion of the big image
    imshow(img(1:350,1:450))
    %imtool(img)

    %%
    % Show a sample of an image before and after deskew is applied
    figure
    subplot(121), imshow(imgs{1}), title('original')
    subplot(122), imshow(deskew(imgs{1})), title('deskewed')

    %% HoG
    % choose an HoG implementation to use: OpenCV, MATLAB, or a simple custom one
    HOG_ALG = 'opencv';
    HOG_ALG = validatestring(HOG_ALG, {'custom', 'matlab', 'opencv'});

    %%
    % HoG features are stored in a matrix, one row per image
    switch HOG_ALG
        case 'custom'
            nbins = 16;          % number of bins to quantize gradient angles
            nfeats = 4 * nbins;  % feature length
        case 'matlab'
            opts = {'CellSize',[10 10], 'BlockSize',[2 2], ...
                'NumBins',9, 'UseSignedOrientation',true};
            nfeats = numel(extractHOGFeatures(imgs{1}, opts{:}));
        case 'opencv'
            % create HOG object
            opts = {'WinSize',[20 20], 'BlockSize',[10 10], ...
                'BlockStride',[5 5], 'CellSize',[10 10], ...
                'NBins',9, 'SignedGradient',true};
            %opts = [opts, 'GammaCorrection',false];
            hog = cv.HOGDescriptor(opts{:});
            nfeats = numel(hog.compute(imgs{1}));
            display(hog)
    end
    descs = zeros(numel(imgs), nfeats, 'single');

    %%
    % Process images and compute HOG descriptors for each
    tic
    disp('Extracting...')
    for i=1:numel(imgs)
        % deskew image
        img = deskew(imgs{i});
        %img = cv.threshold(img, 'Otsu');

        % compute HOG descriptors flattened as a vector
        switch HOG_ALG
            case 'custom'
                descs(i,:) = my_hog(img, nbins);
            case 'matlab'
                descs(i,:) = extractHOGFeatures(img, opts{:});
            case 'opencv'
                d = hog.compute(img).'; % returns a matrix, one row per window
                descs(i,:) = d(:);
        end
    end
    toc

    %% Classification
    % Create multi-class linear SVM classifier
    svm = cv.SVM();
    svm.Type = 'C_SVC';
    svm.KernelType = 'Linear';
    svm.C = 0.5;  % or do a trainAuto grid search to find optimal C
    display(svm)

    %%
    % split data samples and corresponding labels
    labels = int32(labels(:));
    N = 2500;  % fix(numel(labels)/2)
    whos descs labels
    if ~mexopencv.isOctave() && mexopencv.require('stats')
        % train/test split
        tabulate(labels(1:N))
        tabulate(labels(N+1:end))
    end

    %%
    % use half as training data
    tic
    disp('Training...')
    svm.train(descs(1:N,:), labels(1:N));
    toc

    %%
    % the other half for testing
    tic
    disp('Testing...')
    yhat = svm.predict(descs(N+1:end,:));
    toc
    %TODO: https://savannah.gnu.org/bugs/?50716

    %%
    % performance on test set
    acc = nnz(yhat == labels(N+1:end));
    confmat = accumarray([labels(N+1:end), yhat]+1, 1);
    display(confmat)
    fprintf('Accuracy = %f%%\n', (acc/N)*100);

    %%
    % Show some of the wrong predictions
    idx = (yhat ~= labels(N+1:end));
    fprintf('Number of wrong predictions = %d / %d\n', nnz(idx), N);
    figure
    idx = find(idx, 3*3, 'first');
    for i=1:numel(idx)
        ii = idx(i);
        subplot(3,3,i), imshow(imgs{ii+N})
        title(sprintf('%d \\rightarrow %d', labels(ii+N), yhat(ii)))
    end
end

%% Helper functions

function [imgs, labels, img] = load_mnist_digits()
    % Load MNIST handwritten digits, one big image
    fname = fullfile(mexopencv.root(), 'test', 'digits.png');
    if exist(fname, 'file') ~= 2
        disp('Downloading image...')
        url = 'https://cdn.rawgit.com/opencv/opencv/3.2.0/samples/data/digits.png';
        urlwrite(url, fname);
    end
    img = cv.imread(fname, 'Grayscale',true);

    % split it into 5000 small images (500 from each 0:9 digits),
    % each image is 20x20 pixels
    [h,w] = size(img);
    imgs = mat2cell(img, ones(1,50)*h/50, ones(1,100)*w/100);
    imgs = reshape(imgs', 500, [])';  % cell array of size 10x500
    labels = repmat((0:9)', 1, 500);  % each row is one digit
end

function img = deskew(img)
    % deskew image using its second order moments
    sz = 20;
    m = cv.moments(img);
    if abs(m.mu02) > 0.2  % 0.01
        skew = m.mu11 / m.mu02;
        M = [1 skew -0.5*sz*skew; 0 1 0];
        img = cv.warpAffine(img, M, 'DSize',[sz sz], 'WarpInverse',true);
    end
end

function histo = my_hog(img, nbins)
    % this is a simplified implementation of HoG descriptors
    % (fixed 10x10 cells, 1 block, no overlapping, no block normalization)

    % number of orientation histogram bins
    if nargin < 2, nbins = 16; end

    % compute image gradient (its magnitude and direction)
    gx = cv.Sobel(img, 'DDepth','single', 'XOrder',1, 'YOrder',0);
    gy = cv.Sobel(img, 'DDepth','single', 'XOrder',0, 'YOrder',1);
    mag = hypot(gx, gy);
    ang = atan2(gy, gx) + pi;  % directions in [0,2*pi] (signed gradient)

    % quantize directions into 16 bins
    bins = fix((nbins-1) * ang/(2*pi)) + 1; % bin indices in [1,16] inclusive

    % divide 20x20 image into 4 10x10 cells, and in each,
    % accumulate gradient magnitudes inside each direction bin
    % (binned orientations each weighted by its magnitude contribution)
    bb = mat2cell(bins, [10 10], [10 10]);
    mm = mat2cell(mag, [10 10], [10 10]);
    histo = cell(size(bb));
    for i=1:numel(bb)
        H = accumarray(bb{i}(:), mm{i}(:), [nbins 1]);
        if false && any(H)
            % transform histogram to space with Hellinger metric
            H = sqrt(H ./ sum(H));
            H = H ./ norm(H);
        end
        histo{i} = H;
    end

    % stack features from cells into one tall vector of length 4*16
    histo = cat(1, histo{:});
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function test_suite=test_moxunit_util_platform_version()
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_platform_version
    % get start of version
    v=moxunit_util_platform_version();
    assert(isnumeric(v));
    assert(numel(v)>=2);
    version_start=sprintf('%d.%d',v(1:2));

    % now compare with output from 'ver'
    if moxunit_util_platform_is_octave()
        name='Octave';
    else
        name='MATLAB';
    end

    version_struct=ver(name);
    version_s=version_struct.Version;
    n=numel(version_start);

    assertEqual(version_s(1:n),version_start);

function v=dots_to_vec(s)
    parts=regexp(s,'\.','split');
    v=cellfun(@str2num,parts);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%% Simulation Data
simu = simulationClass();               % Initialize Simulation Class
simu.simMechanicsFile = 'RM3.slx';      % Specify Simulink Model File
simu.mode = 'normal';                   % Specify Simulation Mode ('normal','accelerator','rapid-accelerator')
simu.explorer = 'on';                   % Turn SimMechanics Explorer (on/off)
simu.startTime = 0;                     % Simulation Start Time [s]
simu.rampTime = 100;                    % Wave Ramp Time [s]
simu.endTime = 400;                     % Simulation End Time [s]
simu.solver = 'ode4';                   % simu.solver = 'ode4' for fixed step & simu.solver = 'ode45' for variable step
simu.dt = 0.1; 							% Simulation time-step [s]

%% Wave Information
% % noWaveCIC, no waves with radiation CIC
% waves = waveClass('noWaveCIC');       % Initialize Wave Class and Specify Type

% Regular Waves
waves = waveClass('regular');           % Initialize Wave Class and Specify Type
waves.height = 2.5;                     % Wave Height [m]
waves.period = 8;                       % Wave Period [s]

% % Regular Waves with CIC
% waves = waveClass('regularCIC');          % Initialize Wave Class and Specify Type
% waves.height = 2.5;                       % Wave Height [m]
% waves.period = 8;                         % Wave Period [s]

% % Irregular Waves using PM Spectrum
% waves = waveClass('irregular');           % Initialize Wave Class and Specify Type
% waves.height = 2.5;                       % Significant Wave Height [m]
% waves.period = 8;                         % Peak Period [s]
% waves.spectrumType = 'PM';                % Specify Wave Spectrum Type

% % Irregular Waves using JS Spectrum with Equal Energy and Seeded Phase
% waves = waveClass('irregular');           % Initialize Wave Class and Specify Type
% waves.height = 2.5;                       % Significant Wave Height [m]
% waves.period = 8;                         % Peak Period [s]
% waves.spectrumType = 'JS';                % Specify Wave Spectrum Type
% waves.bem.option = 'EqualEnergy';         % Uses 'EqualEnergy' bins (default)
% waves.phaseSeed = 1;                      % Phase is seeded so eta is the same

% % Irregular Waves using PM Spectrum with Traditional and State Space
% waves = waveClass('irregular');           % Initialize Wave Class and Specify Type
% waves.height = 2.5;                       % Significant Wave Height [m]
% waves.period = 8;                         % Peak Period [s]
% waves.spectrumType = 'PM';                % Specify Wave Spectrum Type
% simu.stateSpace = 1;                      % Turn on State Space
% waves.bem.option = 'Traditional';         % Uses 1000 frequnecies

% % Irregular Waves with imported spectrum
% waves = waveClass('spectrumImport');      % Create the Wave Variable and Specify Type
% waves.spectrumFile = 'spectrumData.mat';  % Name of User-Defined Spectrum File [:,2] = [f, Sf]

% % Waves with imported wave elevation time-history
% waves = waveClass('elevationImport');          % Create the Wave Variable and Specify Type
% waves.elevationFile = 'elevationData.mat';     % Name of User-Defined Time-Series File [:,2] = [time, eta]

%% Body Data
% Float
body(1) = bodyClass('hydroData/rm3.h5');
    % Create the body(1) Variable, Set Location of Hydrodynamic Data File
    % and Body Number Within this File.
body(1).geometryFile = 'geometry/float.stl';    % Location of Geomtry File
body(1).mass = 'equilibrium';
    % Body Mass. The 'equilibrium' Option Sets it to the Displaced Water
    % Weight.
body(1).inertia = [20907301 21306090.66 37085481.11];  % Moment of Inertia [kg*m^2]

% Spar/Plate
body(2) = bodyClass('hydroData/rm3.h5');
body(2).geometryFile = 'geometry/plate.stl';
body(2).mass = 'equilibrium';
body(2).inertia = [94419614.57 94407091.24 28542224.82];

%% PTO and Constraint Parameters
% Floating (3DOF) Joint
constraint(1) = constraintClass('Constraint1'); % Initialize Constraint Class for Constraint1
constraint(1).location = [0 0 0];               % Constraint Location [m]

% Translational PTO
pto(1) = ptoClass('PTO1');                      % Initialize PTO Class for PTO1
pto(1).stiffness = 0;                           % PTO Stiffness [N/m]
pto(1).damping = 1200000;                       % PTO Damping [N/(m/s)]
pto(1).location = [0 0 0];                      % PTO Location [m]

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function [gKern, gVarmeans, gVarcovars, gInd] = whiteVardistPsi1Gradient(whitekern, vardist, Z, covGrad, learnInducing)

% WHITEVARDISTPSI1GRADIENT Compute gradient of white variational Psi1.
% FORMAT
% DESC description here.
% RETURN gKern :
% RETURN gVarmeans :
% RETURN gVarcovars :
% RETURN gInd :
% ARG whitekern : the kernel structure associated with the white kernel.
% ARG vardist :
% ARG Z :
% ARG covGrad :
%
% SEEALSO :
%
% COPYRIGHT : Michalis K. Titsias, 2009
%

% VARGPLVM

gKern = 0;
gVarmeans = zeros(1,prod(size(vardist.means)));
gInd = zeros(1,prod(size(Z)));
gVarcovars = zeros(1,prod(size(vardist.covars)));

author: SheffieldML

repo: vargplvm


text:
function [valFeature] = batch_peridoc_pitch_count_fast(testWave, Fs, numFrame, nSamplePerFrame, nSampleForward)
x_start   = 1;
x_end     = nSamplePerFrame;
act_frame_x = zeros(nSamplePerFrame*2,numFrame);
for j=1:numFrame
    act_frame_x(:,j) =  [testWave(x_start:x_end); zeros(x_end-x_start+1,1)];
    x_start = x_start + nSampleForward;
    x_end   = x_end   + nSampleForward;
end

[frSize, nFr] = size(act_frame_x);
binWidth = Fs/frSize;
hammWinFrame  = hamming(frSize);
HighPassFreqUpper = 800;
numHighPassFreqUpperBin = ceil(HighPassFreqUpper/binWidth);
mapX=[1:numHighPassFreqUpperBin];
mapX = 1/numHighPassFreqUpperBin.*mapX;
mapX = mapX';

frame_x = zeros(nSamplePerFrame*2,numFrame);
for i=1:nFr
    frame_x(:,i) = act_frame_x(:,i)/norm(act_frame_x(:,i));
% so that signal level is NOT at play here!!!
end

frame_x_hamm  = bsxfun(@times, frame_x, hammWinFrame);
abs_fft_x = abs(fft(frame_x_hamm));
abs_fft_x(1:numHighPassFreqUpperBin,:) = bsxfun(@times, abs_fft_x(1:numHighPassFreqUpperBin,:), mapX);
abs_fft_x = bsxfun(@times, abs_fft_x, 1./max(abs_fft_x(1:2*HighPassFreqUpper/binWidth,:)));

% The above operation with mapX reduces low frequency power, basically
% frequency filtering!!!

% extract noise level from 2K-3K freq range
abs_fft_x = abs_fft_x.^1.2;

L = floor(frSize/2);
binWidth = (Fs/(L*2));
PitchRangeLower = 80;
PitchRangeUpper = 250;

PitchIdxRange = ceil(PitchRangeUpper - PitchRangeLower)/binWidth;
val_histPeak = [];
val_histTrough = [];

SumPeakVal = [];
SumTroughVal = [];
beamWidth =2;

for (i=1:PitchIdxRange)
    pitchVal  = floor( ((PitchRangeLower/binWidth) +i-1)*binWidth);

    % number of harmonics, skipping the first harmonics bcos very noisy!!!
    for (j=2:7)
        idx          = floor((j*pitchVal)/binWidth);
        idx_trough   = floor(idx + (pitchVal/(2*binWidth)));

        Pi  = sum(abs_fft_x(idx-beamWidth:idx+beamWidth,:));
        Ti =  sum(abs_fft_x(idx_trough-beamWidth:idx_trough+beamWidth,:));

        SumPeakVal(:,i,j-1)   =  Pi;
        SumTroughVal(:,i,j-1) = Ti;
    end

end
SP = sum(SumPeakVal,3);
ST = sum(SumTroughVal,3);
val_histTrough  = ST.^2;

for i=1:PitchIdxRange
    currSumPeakVal = squeeze(SumPeakVal(:,i,:))';
    currSumTroughVal = squeeze(SumTroughVal(:,i,:))';
    val_histPeak(:,i)   = SP(:,i).^2 - var(currSumPeakVal)' - var(currSumTroughVal)';
end

[V1 iV1] = max(val_histPeak');
for i=1:nFr
    V2(i) = abs(val_histTrough(i,iV1(i)));
end
valFeature      = V1./V2;

author: singaxiong

repo: SignalGraph


text:
function []=panel1fdisp(N,n,T,Units,endo,Y,stringdates2,decimaldates2,Fstartlocation,Fendlocation,forecast_estimates,pref)

% preliminary task: gather in a cell the values to be plotted
% initiate the cell
plotdata={};
% because forecasts have to be computed for each unit, loop over units
for ii=1:N
% each cell entry is a matrix of actual and forecast values
% this matrix comprises 4 rows: the first row is actual data, while the three other rows are the estimates (point estimates and confidence bands) for the forecasts
% also, this matrix has a number of rows equal to the dimension of decimaldates2, which comprises the total period sample+forecasts
   % loop over variables
   for jj=1:n
      plotdata{jj,1,ii}=nan(4,size(decimaldates2,1));
      % record actual sample values
      plotdata{jj,1,ii}(1,1:T)=Y(:,jj,ii)';
      % copy the last point of the actual sample for the forecast part of the matrix (required to have a clean plot)
      plotdata{jj,1,ii}(:,Fstartlocation-1)=repmat(Y(Fstartlocation-1,jj,ii),4,1);
      % record forecast, lower bound
      plotdata{jj,1,ii}(2,Fstartlocation:Fendlocation)=forecast_estimates{jj,1,ii}(1,:);
      % record forecast, point estimate
      plotdata{jj,1,ii}(3,Fstartlocation:Fendlocation)=forecast_estimates{jj,1,ii}(2,:);
      % record forecast, upper bound
      plotdata{jj,1,ii}(4,Fstartlocation:Fendlocation)=forecast_estimates{jj,1,ii}(3,:);
   end
end
if pref.plot
% then plot the figure
forecast=figure('Tag','BEARresults');
set(forecast,'Color',[0.9 0.9 0.9]);
set(forecast,'name','unconditional forecasts');
% initiate the count
count=0;
% loop over units
for ii=1:N
   % loop over endogenous variables
   for jj=1:n
   % increment count
   count=count+1;
   % then plot
   subplot(N,n,count)
   hold on
   Xpatch=[decimaldates2(Fstartlocation-1:Fendlocation,1)' fliplr((decimaldates2(Fstartlocation-1:Fendlocation,1))')];
   Ypatch=[plotdata{jj,1,ii}(2,Fstartlocation-1:Fendlocation) fliplr(plotdata{jj,1,ii}(4,Fstartlocation-1:Fendlocation))];
   Fpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
   set(Fpatch,'facealpha',0.6);
   set(Fpatch,'edgecolor','none');
   plot(decimaldates2,plotdata{jj,1,ii}(3,:),'Color',[0.4 0.4 1],'LineWidth',2);
   plot(decimaldates2,plotdata{jj,1,ii}(1,:),'Color',[0 0 0],'LineWidth',2);
   hold off
   set(gca,'XLim',[decimaldates2(1,1) decimaldates2(end,1)],'FontName','Times New Roman');
   set(gca,'XGrid','on');
   set(gca,'YGrid','on');
      % top labels
      if count<=n
      title(endo{count,1},'FontWeight','normal');
      end
      % side labels
      if jj==1
      ylabel(Units{ii,1},'FontWeight','normal');
      end
   end
end
end

% save on Excel
% create the cell that will be saved on excel
forecastcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},size(stringdates2,1)+3,1);
horzspace=repmat({''},3,6*n);
% loop over units
for ii=1:N
% initiate the cell of results
unitcell={};
   % loop over endogenous variables (horizontal dimension)
   for jj=1:n
   % create a header
   header=[{[Units{ii,1} ': ' endo{jj,1}]} {''} {''} {''} {''};{''} {''} {''} {''} {''};{''} {'actual'} {'lower bound'} {'median'} {'upper bound'}];
   % complete the cell
   endocell=[[header;stringdates2 num2cell((plotdata{jj,1,ii})')] vertspace];
   % concatenate to the previous parts of unitcell
   unitcell=[unitcell endocell];
   end
% concatenate to the previous parts of afcell
forecastcell=[forecastcell;horzspace;unitcell];
end
% trim
forecastcell=forecastcell(4:end,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),forecastcell,'forecasts','B2');
end


text:
function [ o, c, f ] = lp_coefficients ( n )

%*****************************************************************************80
%
%% LP_COEFFICIENTS: coefficients of Legendre polynomials P(n,x).
%
%  Discussion:
%
%    The Legendre polynomial with index N will have O = 1 + (N/2) terms.
%    The monomials of orders N, N-2, N-2, ... will have nonzero coefficients.
%
%  First terms:
%
%     1
%     0     1
%    -1/2   0      3/2
%     0    -3/2    0     5/2
%     3/8   0    -30/8   0     35/8
%     0    15/8    0   -70/8    0     63/8
%    -5/16  0    105/16  0   -315/16   0    231/16
%     0   -35/16   0   315/16   0   -693/16   0    429/16
%
%     1.00000
%     0.00000  1.00000
%    -0.50000  0.00000  1.50000
%     0.00000 -1.50000  0.00000  2.5000
%     0.37500  0.00000 -3.75000  0.00000  4.37500
%     0.00000  1.87500  0.00000 -8.75000  0.00000  7.87500
%    -0.31250  0.00000  6.56250  0.00000 -19.6875  0.00000  14.4375
%     0.00000 -2.1875   0.00000  19.6875  0.00000 -43.3215  0.00000  26.8125
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Daniel Zwillinger, editor,
%    CRC Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996.
%
%  Parameters:
%
%    Input, integer N, the degree of the polynomial.
%
%    Output, integer O, the number of coefficients.
%
%    Output, real C(1:O,1), the coefficients of the Legendre polynomial
%    of degree N.
%
%    Output, integer F(1:O,1), the exponents.
%

%
%  Compute the table.
%
  ctable(1:n+1,1:n+1) = 0.0;

  ctable(1,1) = 1.0;

  if ( 1 <= n)

    ctable(2,2) = 1.0;

    for i = 2 : n
      ctable(i+1,1:i-1) =                     (   - i + 1 ) * ctable(i-1,1:i-1) / i;
      ctable(i+1,2:i+1) = ctable(i+1,2:i+1) + ( i + i - 1 ) * ctable(i  ,1:i  ) / i;
    end

  end
%
%  Extract the nonzero data from the alternating columns of the last row.
%
  o = floor ( ( n + 2 ) / 2 );

  c = zeros ( o, 1 );
  f = zeros ( o, 1 );

  k = o;
  for j = n + 1 : -2 : 1
    c(k) = ctable(n+1,j);
    f(k) = j - 1;
    k = k - 1;
  end

  return
end


text:
% Snell's Law And Critical Angle Reflection Example
%
% This example illustrates Snell's law by steering a tone burst from a
% linear array transducer within a layered heterogeneous medium. It builds
% on the Simulating Transducer Field Patterns Example.
%
% author: Bradley Treeby
% date: 10th December 2009
% last update: 24th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: 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.
%
% k-Wave 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.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Nx = 128;           % number of grid points in the x (row) direction
Ny = Nx;            % number of grid points in the y (column) direction
dx = 50e-3/Nx;    	% grid point spacing in the x direction [m]
dy = dx;            % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of a layered propagation medium
medium.alpha_power = 1.5;   % [dB/(MHz^y cm)]
medium.alpha_coeff = 0.75;  % [dB/(MHz^y cm)]
medium.density = 1000;      % [kg/m^3]
c0 = 1500;                  % [m/s]
medium_mulp = 2;
medium.sound_speed = c0*ones(Nx, Ny);
medium.sound_speed(Nx/2:end, :) = medium_mulp*c0;

% create the time array
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);
kgrid.t_array = 0:dt:1000*dt;

% define a source mask for a linear element transducer with an odd number
% of elements
num_elements = 61;
source.p_mask = zeros(Nx, Ny);
x_offset = 25;
y_offset = 20;
start_index = Ny/2 - round(num_elements/2) + 1 - y_offset;
source.p_mask(x_offset, start_index:start_index + num_elements - 1) = 1;

% define the properties of the tone burst
sampling_freq = 1/dt;   % [Hz]
steering_angle = 20;    % [deg]
element_spacing = dx;   % [m]
tone_burst_freq = 1e6;  % [Hz]
tone_burst_cycles = 8;

% create an element index relative to the centre element of the transducer
element_index = -(num_elements - 1)/2:(num_elements - 1)/2;

% use geometric beam forming to calculate the tone burst offsets for each
% transducer element based on the element index
tone_burst_offset = 200 + element_spacing*element_index*sin(steering_angle*pi/180)/(c0*dt);

% create the tone burst signals
source.p = toneBurst(sampling_freq, tone_burst_freq, tone_burst_cycles, 'SignalOffset', tone_burst_offset);

% assign the input options
input_args = {'DisplayMask', source.p_mask};

% run the simulation
kspaceFirstOrder2D(kgrid, medium, source, [], input_args{:});

repo: Aero-matlab


text:
function [post nlZ dnlZ] = infExactWarp(hyp, mean, cov, lik, x, y)

% Exact inference for a GP with Gaussian likelihood. Compute a parametrization
% of the posterior, the negative log marginal likelihood and its derivatives
% w.r.t. the hyperparameters. See also "help infMethods".
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2015-07-13.
%                                      File automatically generated using noweb.
%
% See also INFMETHODS.M.

if iscell(lik), likstr = lik{1}; else likstr = lik; end
if ~ischar(likstr), likstr = func2str(likstr); end
if ~strcmp(likstr,'likGaussWarpExact')               % NOTE: no explicit call to likGaussWarpExact
  error('Exact inference only possible with warped exact Gaussian likelihood');
end
warp = lik{2}; if ischar(warp); warp = str2func(warp); end

ng = feval(warp{:});               % number of hyperparameters for the warping function
nhyp = feval(lik{:});              % number of hyperparameters
% if nargin<4, varargout = {nhyp}; return, end       % report number of parameters
nhyp = eval(nhyp);
if nhyp>length(hyp.lik), error('not enough hyperparameters'), end

[gy,lgpy] = feval(warp{:},y,hyp.lik(1:ng));                     % evaluate warping function

% figure(100); scatter(y,gy); drawnow;

[n, D] = size(x);
K = feval(cov{:}, hyp.cov, x);                      % evaluate covariance matrix
m = feval(mean{:}, hyp.mean, x);                    % evaluate mean vector
[gm,lgpm] = feval(warp{:},m,hyp.lik(1:ng));         % evaluate warped mean

sn2 = exp(2*hyp.lik(end));                          % noise variance of likGauss
if sn2<1e-6                        % very tiny sn2 can lead to numerical trouble
  L = chol(K+sn2*eye(n)); sl =   1;   % Cholesky factor of covariance with noise
  pL = -solve_chol(L,eye(n));                            % L = -inv(K+inv(sW^2))
else
  L = chol(K/sn2+eye(n)); sl = sn2;                       % Cholesky factor of B
  pL = L;                                           % L = chol(eye(n)+sW*sW'.*K)
end
alpha = solve_chol(L,gy-gm)/sl;

post.alpha = alpha;                            % return the posterior parameters
post.sW = ones(n,1)/sqrt(sn2);                  % sqrt of noise precision vector
post.L = pL;

if nargout>1                               % do we want the marginal likelihood?
  nlZ = (gy-gm)'*alpha/2 + sum(log(diag(L))) + n*log(2*pi*sl)/2 - sum(lgpy);   % -log marg lik
  if nargout>2                                         % do we want derivatives?
    dnlZ = hyp;                                 % allocate space for derivatives
    Q = solve_chol(L,eye(n))/sl - alpha*alpha';     % precompute for convenience
    for i = 1:numel(hyp.cov)
      dnlZ.cov(i) = sum(sum(Q.*feval(cov{:}, hyp.cov, x, [], i)))/2;
    end
    for i = 1:ng
        [dgy,dlgpy] = feval(warp{:},y,hyp.lik(1:ng),i);
        dgm = feval(warp{:},m,hyp.lik(1:ng),i);
        dnlZ.lik(i) = -sum(dlgpy) + (dgy - dgm)'*alpha;
    end
    dnlZ.lik(ng+1) = sn2*trace(Q);
    for i = 1:numel(hyp.mean)
      dnlZ.mean(i) = -(exp(lgpm)'.*feval(mean{:}, hyp.mean, x, i)')*alpha;
    end
  end

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [acc_t,predictions_hard,predictions_soft,genplot] = standard_classifier_test(model,X,Y,T,options,binsize)
% Fit an already trained classifier given by the structure model to the
% test data given in X with labels Y.
%
% INPUT
% model: a classification model outputted from the function
%           standard_classifier_train
% X: Brain data, (total time by regions) or (time by trials by regions)
% Y: Stimulus, (total time by q). This variable can either be entered as a
%       single vector of class assignments (ie q=1 and y takes values
%       1,2,...N=total number of classes); or can be entered as a matrix of
%       labels where each column contains 1 or 0 denoting whether that
%       class is active.
% T: Length of series
%       options             structure with the following subfields:
%           classifier      String indicating calssifer type;
%                   logistic    Logistic Regression classifier (default)
%                   SVM         Linear support vector machine classifier
%                   LDA         Linear Discriminant Analysis classifier
%           regularisation  String indicating type of regularisation;
%                   L1          Lasso regularisation
%                   L2          Ridge regularisation (default)
%                   ARD         Automatic relevance determination
%           lambda          Vector of regularisation penalty values to be
%                           tested
% binsize: NOT YET IMPLMENTED
%           how many consecutive time points will be used for each estiamtion.
%           By default, this is 1, such that one decoding model
%               is estimated using 1 time point
%
% OUTPUT
%
%
% Author: Cam Higgins, OHBA, University of Oxford (2019)

if nargin < 5 || isempty(options), options = struct(); end
if nargin < 6, binsize = 1; end

if ~all(T==T(1)), error('All elements of T must be equal for time-aligned classification'); end
N = length(T); ttrial = T(1);

if size(Y,1) < size(Y,2); Y = Y'; end
q = size(Y,2);
if size(Y,1) == length(T) % one value per trial
    responses = Y;
elseif length(Y(:)) ~= (ttrial*N*q)
    error('Incorrect dimensions in Y')
else
    responses = reshape(Y,[ttrial N q]);
    Ystar = reshape(repmat(responses(1,:,:),[ttrial,1,1]),[ttrial*N q]);
    responses = permute(responses(1,:,:),[2 3 1]); % N x q
    if ~all(Y(:)==Ystar(:))
        error('For time-aligned classification, the same stimulus must be presented for the entire trial');
    end
end

Ycopy = Y;
if size(Ycopy,1) == N
    Ycopy = repmat(reshape(Ycopy,[1 N q]),[ttrial 1 1]);
else
    Ycopy = reshape(Ycopy,[ttrial N q]);
end
% no demeaning by default if this is a classification problem
if ~isfield(options,'demeanstim'), options.demeanstim = 0; end

options.Nfeatures = 0;
options.K = 1;
%[X,Y,T] = preproc4hmm(X,Y,T,options); % this demeans Y
p = size(X,2);

%reshape to 3d vectors:
X = reshape(X,[ttrial N p]);
Y = reshape(Y,[ttrial N q]);

if (binsize/2)==0
    warning(['binsize must be an odd number, setting to ' num2str(binsize+1)])
    binsize = binsize + 1;
end

if any(T~=ttrial), error('All trials must have the same length'); end
if nargin<5, binsize = 1; end

if ~isfield(model,'classifier')
    error('Model inputted does not conform to requirements; it should be trained using the function standard_classifier_train.m');
else
    classifier = model.classifier;
    if ~(strcmp(classifier,'logistic') || strcmp(classifier,'SVM') ||...
            strcmp(classifier,'LDA') || strcmp(classifier,'KNN') || ...
            strcmp(classifier,'decisiontree') || strcmp(classifier,'SVM_rbf'));
        error('model.classifier can only take values "logistic", "SVM" or "LDA"');
    end
end
if strcmp(classifier,'logistic')
    Q_star = size(model.betas,3); %check if multinomial regression
else
    Q_star=q;
end
if ~isfield(options,'generalisationplot'),options.generalisationplot=false;end
if isfield(model,'lambda')
    L=length(model.lambda);
else
    L=1;
end

halfbin = floor(binsize/2);
%nwin = round(ttrial / binsize);
%binsize = floor(ttrial / nwin);
%RidgePen = lambda * eye(p);

if strcmp(classifier,'logistic')
    Y(Y==-1)=0;
    Y_pred = zeros(ttrial,N,Q_star,L);
    Y_pred_genplot=zeros(ttrial,ttrial,N,Q_star,L);
    for iLambda=1:L
        for t=1:ttrial
            if p>1
                Y_pred(t,:,:,iLambda)=squeeze(X(t,:,:))*permute(model.betas(t,:,:,iLambda),[2,3,1,4]) + repmat(model.intercepts(t,:,iLambda),N,1);
            else
                Y_pred(t,:,:,iLambda)=X(t,:)'*permute(model.betas(t,:,:,iLambda),[2,3,1,4]) + repmat(model.intercepts(t,:,iLambda),N,1);
            end
            if options.generalisationplot
                for t2=1:ttrial
                    Y_pred_genplot(t,t2,:,:,iLambda)= ...
                        squeeze(X(t2,:,:))*permute(model.betas(t,:,:,iLambda),[2,3,1,4]) + repmat(model.intercepts(t,:,iLambda),N,1);
                end
            end
        end
    end

    for iL=1:L
        if Q_star~=q
            pred_temp=multinomLogRegPred(Y_pred(:,:,:,iL));
            predictions_soft(:,:,iL)=reshape(pred_temp,[ttrial*N,q]);
            predictions_hard(:,:,iL)=hardmax(predictions_soft(:,:,iL));
        else
            Y_pred = reshape(Y_pred,[ttrial*N,q,L]);
            predictions_soft(:,:,iL)=log_sigmoid(Y_pred(:,:,iL));
            if Q_star>1
                predictions_hard(:,:,iL)=hardmax(Y_pred(:,:,iL));
            else
                predictions_hard(:,:,iL) = Y_pred(:,:,iL)>0;
            end
        end
    end
    if options.generalisationplot
        if Q_star==q
            Y_pred_genplot = reshape(Y_pred_genplot,[ttrial, ttrial*N,Q_star,L]);
            for iL=1:iL
                for t=1:ttrial
                    Y_pred_genplot(t,:,:,iL) = hardmax(squeeze(Y_pred_genplot(t,:,:,iL)));
                end
            end
        else %multinomial:
            Y_pred_genplot_q = zeros([ttrial, ttrial,N,q,L]);
            for iL=1:iL;
                for t=1:ttrial
                    pred_temp = multinomLogRegPred(squeeze(Y_pred_genplot(t,:,:,:,iL)));
                    pred_temp = reshape(pred_temp,[ttrial*N,q]);
                    Y_pred_genplot_q(t,:,:,:,iL) = reshape(hardmax(pred_temp),[ttrial,N,q]);
                end
            end
            Y_pred_genplot=Y_pred_genplot_q;
        end
    end
elseif strcmp(classifier,'SVM') || strcmp(classifier,'SVM_rbf')
     % note this uses distance from hyperplane as soft score to allow
     % multiclass categorisation
     Y_pred = zeros(ttrial,N,q);
     Y_pred_genplot=zeros(ttrial,ttrial,N,q);
     for t=1:ttrial
        for iStim=1:q
            %[~,sc] = predict(model.SVM{t,iStim},squeeze(X(t,:,:)));
            [~,sc] = model.SVM{t,iStim}.predict(squeeze(X(t,:,:)));
            Y_pred(t,:,iStim)=sc(:,2);
            if options.generalisationplot
                for t2=1:ttrial
                    %[~,sc] = predict(model.SVM{t,iStim},squeeze(X(t2,:,:)));
                    [~,sc] = model.SVM{t,iStim}.predict(squeeze(X(t2,:,:)));
                    Y_pred_genplot(t,t2,:,iStim)= sc(:,2);
                end
            end
        end
     end
     Y_pred = reshape(Y_pred,[ttrial*N,q,L]);
     predictions_soft=Y_pred;
     predictions_hard=hardmax(Y_pred);
     if options.generalisationplot
        Y_pred_genplot = reshape(Y_pred_genplot,[ttrial, ttrial*N,q]);
        for t=1:ttrial
            if q>1
                Y_pred_genplot(t,:,:) = hardmax(squeeze(Y_pred_genplot(t,:,:)));
            else
                Y_pred_genplot(t,:,:) = squeeze(Y_pred_genplot(t,:,:))>0;
            end
        end
    end
elseif strcmp(classifier,'LDA')
    X = reshape(X,[ttrial*N,p]);
    [predictions_hard, predictions_soft] = LDApredict(model,repmat(eye(ttrial),N,1),X);
    if q>1
        predictions_soft = exp(predictions_soft - repmat(max(predictions_soft,[],2),1,q));
        predictions_soft = rdiv(predictions_soft,sum(predictions_soft,2));
    end
elseif strcmp(classifier,'KNN')
    if ~isfield(model,'K');model.K=1;end
    predictions_hard = zeros(ttrial,N,q);
    for t=1:ttrial
        candidates = squeeze(model.X_train(t,:,:));
        D=dist(squeeze(X(t,:,:)),candidates');
        for n=1:N
            [dists,inds]=sort(D(n,:));
            if options.K==1
                M=model.Y_train(t,inds(1:options.K),:);
                ind_winner = find(squeeze(M));
                predictions_hard(t,n,ind_winner) = 1;
            else
                Y_pred=squeeze(model.Y_train(t,inds(1:options.K),:));
                M = sum(Y_pred,1);
                [M,inds2] = sort(M,'descend');
                if M(1)>M(2)
                    predictions_hard(t,n,inds2(1)) = 1;
                else %tie break - just take closest point
                    ind_winner = find(squeeze(model.Y_train(t,inds(1),:)));
                    predictions_hard(t,n,inds(1)) = 1;
                end
            end
        end
    end
    predictions_hard = reshape(predictions_hard,[ttrial*N,q]);
    predictions_soft = predictions_hard;
elseif strcmp(classifier,'decisiontree');
    predictions_hard = zeros(ttrial,N,q);
    for t=1:ttrial
        temp = model.decisiontree{t}.predict(squeeze(X(t,:,:)));
        for n=1:N
            predictions_hard(t,n,temp(n))=1;
        end
    end
    predictions_hard = reshape(predictions_hard,[ttrial*N,q]);
    predictions_soft = predictions_hard;
end

%and compute accuracy metrics using hard classification output:
Y = reshape(logical(Y==1),[ttrial*N q]);
predictions_hard=logical(predictions_hard);
true_preds = all(~xor(predictions_hard,Y),2);
acc = mean(true_preds);
acc_t = mean(reshape(true_preds,[ttrial,N]),2);
if options.generalisationplot
    Y_true_genplot = repmat(permute(Ycopy,[4,1,2,3]),[ttrial,1,1,1]);
    Y_pred_genplot = reshape(Y_pred_genplot,[ttrial, ttrial,N,q]);
    if q>1
        accplot = sum(logical(Y_true_genplot) & logical(Y_pred_genplot),4);
    else
        accplot = ~xor(logical(Y_true_genplot),logical(Y_pred_genplot));
    end
    genplot = squeeze(mean(accplot,3));
else
    genplot=[];
end
end

function preds_hard = hardmax(Y_pred)
% assuming multiple binomial only; Y_pred of dimension [NT x q]
if ndims(Y_pred)>2
    error('Wrong dimensions entered for computing predictions');
end
if size(Y_pred,2)>1
    [~,a] = max(Y_pred,[],2);
else
    a = Y_pred>0;
end
preds_hard = zeros(size(Y_pred));
for i=1:length(preds_hard)
    preds_hard(i,a(i))=1;
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function [stepsize, newx, newkey, lsstats] = ...
  linesearch_adaptive(problem, x, d, f0, df0, options, storedb, key)
% Adaptive line search algorithm (step size selection) for descent methods.
%
% function [stepsize, newx, newkey, lsstats] =
%        linesearch_adaptive(problem, x, d, f0, df0, options, storedb, key)
%
% Adaptive linesearch algorithm for descent methods, based on a simple
% backtracking method. Contrary to linesearch.m, this function is not
% invariant under rescaling of the search direction d. These two line
% search methods vary mainly in their strategy to pick the initial step
% size.
%
% Below, the step is constructed as alpha*d, and the step size is the norm
% of that vector, thus: stepsize = alpha*norm_d. The step is executed by
% retracting the vector alpha*d from the current point x, giving newx.
%
% This line-search may create and maintain a structure called lsmem inside
% storedb.internal. This gives the linesearch the opportunity to remember
% what happened in the previous calls. This is typically used to make a
% first guess at the step size, based on previous events.
%
% Inputs/Outputs : see help for linesearch
%
% See also: steepestdescent conjugategradients linesearch

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec. 30, 2012.
% Contributors: Nicolas Boumal
% Change log:
%
%   Sept. 13, 2013 (NB) :
%       The automatic direction reversal feature was removed (it triggered
%       when df0 > 0). Direction reversal is a decision that needs to be
%       made by the solver, so it can know about it.
%
%	Nov. 7, 2013 (NB) :
%       The whole function has been recoded to mimick more closely the new
%       version of linesearch.m. The parameters are available through the
%       options structure passed to the solver and have the same names and
%       same meaning as for the base linesearch. The information is logged
%       more reliably.
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.
%
%   April 8, 2015 (NB):
%       Got rid of lsmem input/output: now maintained in storedb.internal.

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey();
    end

    % Backtracking default parameters. These can be overwritten in the
    % options structure which is passed to the solver.
    default_options.ls_contraction_factor = .5;
    default_options.ls_suff_decr = .5;
    default_options.ls_max_steps = 10;
    default_options.ls_initial_stepsize = 1;

    if ~exist('options', 'var') || isempty(options)
        options = struct();
    end
    options = mergeOptions(default_options, options);

    contraction_factor = options.ls_contraction_factor;
    suff_decr = options.ls_suff_decr;
    max_ls_steps = options.ls_max_steps;
    initial_stepsize = options.ls_initial_stepsize;

    % Compute the norm of the search direction.
    norm_d = problem.M.norm(x, d);

    % If this is not the first iteration, then lsmem should have been
    % filled with a suggestion for the initial step.
    if isfield(storedb.internal, 'lsmem')
        lsmem = storedb.internal.lsmem;
        if isfield(lsmem, 'init_alpha')
            % Pick initial step size based on where we were last time,
            alpha = lsmem.init_alpha;
        end
    % Otherwise, fall back to a user supplied suggestion.
    else
        alpha = initial_stepsize / norm_d;
    end

    % Make the chosen step and compute the cost there.
    newx = problem.M.retr(x, d, alpha);
    newkey = storedb.getNewKey();
    newf = getCost(problem, newx, storedb, newkey);
    cost_evaluations = 1;

    % Backtrack while the Armijo criterion is not satisfied
    while newf > f0 + suff_decr*alpha*df0

        % Reduce the step size,
        alpha = contraction_factor * alpha;

        % and look closer down the line
        newx = problem.M.retr(x, d, alpha);
        newkey = storedb.getNewKey();
        newf = getCost(problem, newx, storedb, newkey);
        cost_evaluations = cost_evaluations + 1;

        % Make sure we don't run out of budget
        if cost_evaluations >= max_ls_steps
            break;
        end

    end

    % If we got here without obtaining a decrease, we reject the step.
    if newf > f0
        alpha = 0;
        newx = x;
        newkey = key;
        newf = f0; %#ok<NASGU>
    end

    % As seen outside this function, stepsize is the size of the vector we
    % retract to make the step from x to newx. Since the step is alpha*d:
    stepsize = alpha * norm_d;

    % Fill lsmem with a suggestion for what the next initial step size
    % trial should be. On average we intend to do only one extra cost
    % evaluation. Notice how the suggestion is not about stepsize but about
    % alpha. This is the reason why this line search is not invariant under
    % rescaling of the search direction d.
    switch cost_evaluations
        case 1
            % If things go very well, push your luck.
            init_alpha = 2 * alpha;
        case 2
            % If things go reasonably well, try to keep pace.
            init_alpha = alpha;
        otherwise
            % If we backtracked a lot, the new stepsize is probably quite
            % small: try to recover.
            init_alpha = 2 * alpha;
    end
    storedb.internal.lsmem.init_alpha = init_alpha;

    % Return some statistics also, for possible analysis.
    lsstats.costevals = cost_evaluations;
    lsstats.stepsize = stepsize;
    lsstats.alpha = alpha;

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function x = subsasgn( x, S, y )
error( 'Subscripted assignment not allowed for these objects.' );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function xxdecode = i4_to_xxdecode ( i )

%*****************************************************************************80
%
%% I4_TO_XXDECODE returns the I-th character in the XXDECODE encoding.
%
%  Example:
%
%    I  I4_TO_UUDECODE
%
%    1  '+'
%    2  '-'
%    3  '0'
%   ..
%   64  'z'
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the index of the character.
%    1 <= I <= 64.
%
%    Output, character XXDECODE, the requested character.
%
  string = '+-0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz';

  if ( 1 <= i && i <= 64 )
    xxdecode = string(i);
  else
    xxdecode = ' ';
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [f,k,s0,gof] = relaxFitMt(data,delta,r1,s0,tr,flipAngle,fitMethod,useParfor)
%
% [f,k,gof] = relaxFitMt(data,delta,t1,s0,tr,flipAngle,[fitMethod='f'],[useParfor=false])
%
% Computes a nonlinear fit of the bound-pool (f) map.
%
% fitMethod = 'f' for fminsearch, 'l' for lsqnonlin.
%
% tr should be in seconds, flipAngle in radians. r1 is 1/t1.
%
% Returns:
%
% SEE ALSO:
%
% relaxFitT1.m to fit the t1 and pd maps used by this function.
%
% HISTORY:
% 2008.02.26 RFD: wrote it.
% 2009.12.04 RFD: changed starting values to be a closer to the typical
% white matter values. I also tightened the range of allowable values for k
% to reflect the published values for brain tissue. Both of these changes
% have the effect of reducing crazy f-values, mostly in the csf. They do
% not seem to affect resulting f-values in the white matter, except in
% regions of high noise where the fits are unstable.

if(~exist('fitMethod','var')||isempty(fitMethod))
    fitMethod = 'f';
end
if(~exist('useParfor','var')||isempty(useParfor))
    useParfor = false;
end

% RFD: empirically, [1 .08] is the median for brain tissue. But, we mostly
% care about white matter, so we'll use the typical values for white matter
x0 = [3 .10];
%x0 = [2.4, .10]'; %this is our initial guess, bese [3.4, .15]'

% lb = [.1 .03]; % [k f]
% ub = [5 .28];
% 2009.12.04 RFD: tighten the reigns on k
lb = [0.9 .03]; % [k f]
ub = [4.5 .28];

t_m = 8e-3; %bese 8e-3
t_s = 5e-3; %bese 5e-3
t_r = 19e-3; %bese 19e-3
T2_B = 11e-6;
% The following (omega_1rms) is derived from the MT RF pulse shape. We
% hard-code a value here, but we should check it's validity.
% The angular freq of the MT pulse- the equivalent flip angle of the MT
% pulse (expressed in radians/sec).
%w1rms = 2400;
w1rms = 2400; % 382Hz * 2pi = 2400

for ii = 1:length(delta)
    W_B(ii) = pi*(w1rms^2)*lorentzian(delta(ii), T2_B);
end;

if(fitMethod=='l')
   options = optimset('LargeScale','on','LevenbergMarquardt','on', 'Display', 'off', 'MaxIter', 50);
else
   options = optimset('LargeScale','off','LevenbergMarquardt','on', 'Display', 'off', 'MaxIter', 50);
end

sz = size(data);
if(sz(1)~=numel(delta))
   error('size(data,1) must = numel(delta)!');
end

f = zeros(1,sz(2));
k = zeros(1,sz(2));
gof = zeros(1,sz(2));

% What we compute here is actually W_F./R1_F. R1_F is computed in
% the fit function, but we precompute the rest out here to save a
% few cpu cycles in the loop below.
W_F = (w1rms./(2*pi*delta)).^2/.055;
warning('off','all');
if(useParfor)
    parfor(ii=1:sz(2))
        % Some voxels produce a "Input to EIG must not contain NaN
        % or Inf" error in lsqnonl in. Tweaking the bounds or
        % starting estimate can fix it sometimes, but they are
        % probably junk voxels anyway, so we'll catch and skip them.
        try
            if(fitMethod=='l')
                [x, resnorm, residual, exitflag, output] = lsqnonlin(@(x) relaxMtFitFunc(x, data(:,ii), W_B, W_F, T2_B, r1(ii), s0(ii), t_m, t_s, t_r), x0, lb, ub, options); %bese j-12;
            else
                [x, resnorm, exitflag] = fminsearch(@(x) relaxMtFitFuncLs(x, data(:,ii), W_B, W_F, T2_B, r1(ii), s0(ii), t_m, t_s, t_r), x0, options);
            end
            %disp([ii x])
            if(exitflag>0)
                k(ii) = x(1);
                f(ii) = x(2);
                gof(ii) = resnorm;
            else
                gof(ii) = NaN;
            end
        catch
            % Leave the fit values at zero.
        end
    end
else
    for(ii=1:sz(2))
        try
            if(fitMethod=='l')
                [x, resnorm, residual, exitflag, output] = lsqnonlin(@(x) relaxMtFitFunc(x, data(:,ii), W_B, W_F, T2_B, r1(ii), s0(ii), t_m, t_s, t_r), x0, lb, ub, options); %bese j-12;
            else
                [x, resnorm, exitflag] = fminsearch(@(x) relaxMtFitFuncLs(x, data(:,ii), W_B, W_F, T2_B, r1(ii), s0(ii), t_m, t_s, t_r), x0, options);
            end
            if(exitflag>0)
                k(ii) = x(1);
                f(ii) = x(2);
                gof(ii) = resnorm;
            else
                gof(ii) = NaN;
            end
        catch
           % Leave the fit values at zero.
        end
    end
end
warning on;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
tic

hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);
i_fish = 8;
[cIX_0,gIX_0,M] = LoadSingleFishDefault(i_fish,hfig,[2,1]);

% [cIX_0,gIX_0] = LoadCluster_Direct(i_fish,2,1);

cIX = cIX_0;
gIX = gIX_0;

%% PCA
[coeff,score,latent,tsquared,explained,mu] = pca(M);
% figure;subplot(121);imagesc(coeff);subplot(122);imagesc(score);title('M');

% data to pass on from PCA:
k = 100;
Z0 = score(:,1:k);
% figure;imagesc(Z0);
t1 = toc

%% ICA
tic
numK = 139;%20;
Z = Z0';
[Zica, W, T, mu] = fastICA(Z,numK);

t2 = toc
%%
n = size(M',2);
Zr = T \ W' * Zica + repmat(mu,1,n);

figure;
subplot(131)
imagesc(Zica);
subplot(132)
imagesc(W);
subplot(133)
imagesc(Zr);

%% clustering
gIX = kmeans(Zr',numK,'distance','correlation');

UpdateIndices_Manual(hfig,cIX,gIX);

figure('Position',[50,100,1400,800]);
% isCentroid,isPlotLines,isPlotBehavior,isPlotRegWithTS
subplot(121)
setappdata(hfig,'isPlotBehavior',1);
setappdata(hfig,'isStimAvr',0);
setappdata(hfig,'isPlotLines',0);
UpdateTimeIndex(hfig);
DrawTimeSeries(hfig,cIX,gIX);

% right plot
subplot(122)
I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
DrawCellsOnAnat(I);

% kmeansPlot(Zr,gIX)
%% bad? pick top 100 cells for each component
% gIX = zeros(size(gIX));
% for i = 1:numK
%     [B,IX] = sort(Zica(i,:),'descend');
%     IX_keep = IX(1:100); % this is too crude ~~~~~~
%
%     gIX(IX_keep) = i;
% end
% IX = (gIX~=0);
% cIX = cIX(IX);
% gIX = gIX(IX);

%% pick top ? cells for each component % still crude, no
cutoff_std = 3;

Z2 = Zica;Z2(Z2<cutoff_std)=nan;
% manual code to not assign values to NaN columns (facepalm)
[~,gIX_1] = max(vertcat(nan(1,size(Z2,2)),Z2));
% gIX_1 = gIX_1';
gIX = gIX_1-1;
IX = (gIX~=0);
cIX = cIX_0(IX);
gIX = gIX_1(IX);

%%
% [cIX,gIX] = SelectClusterRange(cIX,gIX,[3,4]);

UpdateIndices_Manual(hfig,cIX,gIX);
%%
figure('Position',[50,100,1400,800]);
% isCentroid,isPlotLines,isPlotBehavior,isPlotRegWithTS
subplot(121)
setappdata(hfig,'isPlotBehavior',1);
setappdata(hfig,'isStimAvr',0);
setappdata(hfig,'isPlotLines',0);
UpdateTimeIndex(hfig);
DrawTimeSeries(hfig,cIX,gIX);

% right plot
subplot(122)
I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
DrawCellsOnAnat(I);

author: xiuyechen

repo: FishExplorer


text:
classdef TestComponent < handle
%TestComponent Abstract base class for TestCase and TestSuite
%
%   TestComponent methods:
%       run          - Run all test cases in test component
%       print        - Display summary of test component to Command Window
%       numTestCases - Number of test cases in test component
%       setUp        - Initialize test fixture
%       tearDown     - Clean up text fixture
%
%   TestComponent properties:
%       Name - Name of test component
%       Location - Directory where test component is defined
%
%   See TestCase, TestSuite

%   Steven L. Eddins
%   Copyright 2008-2009 The MathWorks, Inc.

    properties
       Name = '';
       Location = '';
    end

    properties (Access = 'protected')
        PrintIndentationSize = 4
    end

    methods (Abstract)
       print()
       %print Display summary of test component to Command Window
       %   obj.print() displays information about the test component to the
       %   Command Window.

       run()
       %run Execute test cases
       %   obj.run() executes all the test cases in the test component

       numTestCases()
       %numTestCases Number of test cases in test component
    end

    methods
        function setUp(self)
            %setUp Set up test fixture
            %   test_component.setUp() is called at the beginning of the run()
            %   method.  Test writers can override setUp if necessary to
            %   initialize a test fixture.
        end

        function tearDown(self)
            %tearDown Tear down test fixture
            %   test_component.tearDown() is at the end of the method.  Test
            %   writers can override tearDown if necessary to clean up a test
            %   fixture.
        end

    end
end

author: ricedsp


text:
classdef Sigmoid < dagnn.ElementWise
  methods
    function outputs = forward(obj, inputs, params)
      outputs{1} = vl_nnsigmoid(inputs{1}) ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs{1} = vl_nnsigmoid(inputs{1}, derOutputs{1}) ;
      derParams = {} ;
    end
  end
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
% SVARGPLVMPREDICTIONS2 This script implements the predictions for a classification task
%
% SEEALSO : demClassificationSvargplvm, demClassificationSvargplvm3, demClassificationGeneralSvargplvm, svargplvmPredictions, svargplvmPredictions3
%
% COPYRIGHT : Andreas C. Damianou, 2011

% VARGPLVM

%%%--- Script: svargplvmPredictions2 %%%%-
%%% This script is like svargplvmPredictions with a few additions
if ~exist('testOnTraining')
    testOnTraining=0;
end
if ~exist('globalOpt') || ~isfield(globalOpt, 'predictionsOptionX')
    if ~exist('predictionsOptionX')
        predictionsOptionX = 1;
    end
end
if ~exist('globalOpt') || ~isfield(globalOpt, 'predictionsGenerationOption')
    predictionsGenerationOption = 1;
end
if ~exist('testDisplayIters')
    testDisplayIters=0; %%%
end
if ~exist('predictionVariance')
    predictionVariance = false;
end
if ~exist('predictionXVariance')
    predictionXVariance = true;
end

numberTestPoints = size(Yts{1},1); % 10;
if testOnTraining
   % perm = randperm(model.N);
   % testInd = perm(1:size(Yall{1},1));
   testInd = 1:size(Yall{1},1);
else
    perm = randperm(size(Yts{obsMod},1));
    %testInd = perm(1:numberTestPoints);
    if ~exist('testInd')
        testInd = 1:size(Yts{1},1); %%%%
    end
end

ZpredMuAll = [];
XpredAll = zeros(length(testInd), model.q);  %%%
varZpredAll = zeros(length(testInd), model.q);  %%%
pb = myProgressBar(length(testInd));
for i=1:length(testInd)
    pb = myProgressBar(pb,i);
    curInd = testInd(i);
   % fprintf('# Testing indice number %d ', curInd);
    if testOnTraining
        y_star = model.comp{obsMod}.y(curInd,:);
        x_star = model.comp{obsMod}.vardist.means(curInd,:);
        varx_star = model.comp{obsMod}.vardist.covars(curInd,:);
    else
        y_star = Yts{obsMod}(curInd,:);
        z_star = Yts{infMod}(curInd,:);
        dst = dist2(y_star, model.comp{obsMod}.y);
        [mind, mini(i)] = min(dst);

        model.comp{obsMod}.mini = mini(i);
        Init(i,:) = model.vardist.means(mini(i),:);
        vardistx = vardistCreate(model.comp{obsMod}.vardist.means(mini(i),:), model.q, 'gaussian');
        vardistx.covars = model.comp{obsMod}.vardist.covars(mini(i),:);
        model.comp{obsMod}.vardistx = vardistx;
        iters = globalOpt.reconstrIters;
        % Find p(X_* | Y_*) which is approximated by q(X_*)
        [x_star, varx_star, modelUpdated] = vargplvmOptimisePoint(model.comp{obsMod}, vardistx, y_star, testDisplayIters, iters);%%%
    end
    numberOfNN = 1;
    % Now we selected a datapoint X_* by taking into account only the
    % private dimensions for Y. Now, based on the shared dimensions of
    % that, we select the closest (in a NN manner) X from the training data.
   % fprintf('# Finding the %d NN of X_* with the training X based only on the shared dims.\n', numberOfNN);
    [ind, distInd] = nn_class(model.X(:,sharedDims), x_star(:,sharedDims), numberOfNN, 'euclidean');

    ZpredMu = zeros(length(ind), size(model.comp{infMod}.y,2));
    varZstar = zeros(size(ZpredMu));
    ZpredSigma = zeros(length(ind), size(model.comp{infMod}.y,2));

    % Find p(y_*|x_*) for every x_* found from the NN
   % fprintf('# Predicting from the NN of X_* ');
    for k=1:numberOfNN
        x_cur = svargplvmPredictX(x_star, model.X, ind(k), sharedDims, predictionsOptionX);
        %[ZpredMu(k,:), ZpredSigma(k,:)] =
        %vargplvmPosteriorMeanVar(model.comp{infMod}, model.X(ind(k),:));
        if predictionsGenerationOption
            if ~predictionXVariance
                varx_star = varx_star * 0;
            end
            % NEW: varx_star
            if ~predictionVariance
                ZpredMu(k,:) = vargplvmPosteriorMeanVar(model.comp{infMod},x_cur, varx_star);
            else
                [ZpredMu(k,:) varZstar(k,:)] = vargplvmPosteriorMeanVar(model.comp{infMod},x_cur, varx_star);
            end
        else
            [ind2, distInd] = nn_class(model.X, x_cur, 1,'euclidean');
            ZpredMu(k,:) = model.comp{infMod}.y(ind2,:);
        end
    end
    ZpredMuAll = [ZpredMuAll; ZpredMu(1,:)];
    if predictionVariance
        varZpredAll = [varZpredAll; varZstar(1,:)];
    end
    XpredAll(i,:) = x_cur;
   % fprintf('\n\n');
end
fprintf(1, '\n');

ZpredMuAllOrig = ZpredMuAll;
defThresh = 0.5;
if min(min(Yts{2} == -1))
    defThresh = 0;
end
ZpredMuAll(ZpredMuAll >= defThresh) = 1;
ZpredMuAll(ZpredMuAll < defThresh) = 0;


text:
function this = MrAffineTransformation_constructor(this, testVariants)
% Unit test for MrAffineTransformation Constructor
%
%   Y = MrUnitTest()
%   run(Y, 'MrAffineTransformation_constructor')
%
% This is a method of class MrUnitTest.
%
% IN
%
% OUT
%
% EXAMPLE
%   MrAffineTransformation_constructor
%
%   See also MrUnitTest

% Author:   Saskia Bollmann
% Created:  2017-10-19
% Copyright (C) 2017 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% Unit test for MrAffineTransformation Constructor

switch testVariants
    case 'propVal' % test prop/val syntax

        % actual solution
        actSolution = this.make_affineTransformation_reference(0);

        % expected solution
        % get classes path
        classesPath = tapas_uniqc_get_path('classes');
        % make full filename
        solutionFileName = fullfile(classesPath, '@MrUnitTest' , 'affineTrafo.mat');
        expSolution = load(solutionFileName);
        expSolution = expSolution.affineTrafo;

        % verify equality of expected and actual solution
        % import matlab.unittests to apply tolerances for objects
        import matlab.unittest.TestCase
        import matlab.unittest.constraints.IsEqualTo
        import matlab.unittest.constraints.AbsoluteTolerance
        import matlab.unittest.constraints.PublicPropertyComparator

        this.verifyThat(actSolution, IsEqualTo(expSolution,...
            'Within', AbsoluteTolerance(10e-7),...
            'Using', PublicPropertyComparator.supportingAllValues));

    case 'matrix' % test affine transformation as input
        % expected solution
        expSolution = this.make_affineTransformation_reference(0);
        expSolution = expSolution.affineMatrix;
        % actual solution
        % make actual solution from affine matrix of expected solution
        actSolution = MrAffineTransformation(expSolution);
        actSolution = actSolution.affineMatrix;

        % verify equality of expected and actual solution
        this.verifyEqual(actSolution, expSolution, 'absTol', 10e-7);
end
end

author: translationalneuromodeling

repo: tapas


text:
function[L,ia,ib,num]=blocklen(x,delta)
%BLOCKLEN Counts the lengths of 'blocks' in an array.
%
%   Suppose X is a column vector which contains blocks of identical
%   values, e.g. X=[0 0 0 1 1 3 2 2 2]';
%
%   L=BLOCKLEN(X) counts the lengths of contiguous blocks containing
%   identical values of X, and returns an array L of size SIZE(X).
%   Each elements of L specifies the length of the block to which the
%   corresponding element of X belongs.
%
%   In the above example, L=[3 3 3 2 2 1 3 3 3]';
%
%   [L,IA,IB,NUM]=BLOCKLEN(X) optionally returns indices IA and IB
%   into the first and last elements, respectively, of each block, as
%   as well as the block number NUM. NUM is the same size as L.
%
%   [...]=BLOCKLEN(X,D) defines the junction between two blocks as
%   locations where ABS(DIFF(X))>D.  Thus D=1 is a 'rate of change'
%   definition, and X=[1 2 3 5 6 10 16 17 18]'; will yield the same
%   result for L as in the previous example.
%
%   See also BLOCKNUM.
%
%   Usage: [L,ia,ib]=blocklen(x);
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2000--2014 J.M. Lilly --- type 'help jlab_license' for details

% 06.09.04 JML fixed incorrect IA, IB output length

if strcmpi(x,'--t')
  blocklentest;return;
end

if nargin==1
  delta=0;
end

ii=(1:length(x))';

index=find(diff(x)~=delta);
ia=[1;index+1];
ib=[index;length(ii)];

L=zeros(size(ii));
L(ia)=ib-ia+1;
L=cumsum(L);

L0=zeros(size(ii));
ib2=ib(1:end-1);
ia2=ia(1:end-1);
L0(ib2+1)=ib2-ia2+1;
L0=cumsum(L0);
L=L-L0;

if nargout ==4
    num=zeros(size(x));
    num(ia)=1;
    num=cumsum(num);
end

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

function[]=blocklentest
x=   [0 0 0 1 1 3 2 2 2]';
y=   [3 3 3 2 2 1 3 3 3]';
num= [1 1 1 2 2 3 4 4 4]';

[y2,ia2,ib2,num2]=blocklen(x);
reporttest('BLOCKLEN with D==0',all(y==y2)&&all(num==num2))

x=[1 2 3 5 6 10 16 17 18]';
[y2,ia2,ib2,num2]=blocklen(x,1);
reporttest('BLOCKLEN with D==1',all(y==y2)&&all(num==num2))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
%| denoise_threshold_test.m
%| Test denoising based on l_1 type penalty, cf thresholding
%| min_x 1/2 |y - x|^2 + \beta pot(x)
%|
%| Copyright 2005-4-22, Jeff Fessler, University of Michigan

yi = linspace(-10,10,101)';
mask = true(size(yi));
A = Gdiag(ones(size(mask))); % identity "matrix"

if 1
	f.l2b_q = 1;
%	f.l2b_n = 5;
	f.l2b_n = 1;
	f.cut = 5; % cutoff point
%	f.delta = 0.2;
	f.delta = f.cut / (1 + 2^f.l2b_n); % for broken parabola
	f.tik = (1 + 2^f.l2b_n) * f.delta;
	f.niter = 40;
%	f.type = 'cauchy';
	f.type = 'broken';

	Rq = Reg1(mask, 'type_denom', 'matlab', ...
		'offsets', 0, ... % trick for identity
		'beta', 2^f.l2b_q);
%		'pot_arg', {'quad'}, 'beta', 2^f.l2b_q);

	xq = pwls_sps_os(0*yi(:), yi(:), [], A, Rq, ...
			2, [-inf inf], [], [], 1);
end

if 1
	Rc = Reg1(mask, 'type_denom', 'matlab', ...
		'offsets', 0, ... % trick for identity
		'pot_arg', {f.type, f.delta}, 'beta', 2^f.l2b_n);

	xc = pwls_sps_os(0*yi(:), yi(:), [], A, Rc, ...
		f.niter, [-inf inf], [], [], 1);

	Rh = Reg1(mask, 'type_denom', 'matlab', ...
		'offsets', 0, ... % trick for identity
		'pot_arg', {'hyper3', f.delta}, 'beta', 2^f.l2b_n);

	xh = pwls_sps_os(yi(:), yi(:), [], A, Rh, ...
		f.niter, [-inf inf], [], [], 1);
end

if im
	clf, plot(yi, yi, ':', ...
		yi, xh(:,end), '-', ...
		yi, xc(:,end), '-.', ...
		yi, xq(:,end), '--')
%	axis equal
	axis square
	tik = [min(yi), -f.tik 0 f.tik max(yi)];
%%%%%%%%	xtick(tik), ytick(tik)
	grid
	ir_legend({'I', 'hyper', f.type, 'quad'})
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function ROTY = ROTY(theta)

ROTY = [cos(theta) 0 sin(theta);
        0 1 0;
        -sin(theta) 0 cos(theta)];
end

author: yorgoon

repo: minimum-snap-geometric-control


text:
function [ok,catalog] = ZmapImportManager(fun, funArguments, varargin)
    % ZMAPIMPORTMANAGER - handles the import of zmap catalogs.
    %
    % OK = ZMAPIMPORTMANAGER( FUN ) provide the import function as a parameter. Catalog will replace
    % the primary catalog in ZmapGlobal.Data.  OK indicates whether a catalog was successfully
    % loaded
    %
    % [OK, CATALOG] = ZMAPIMPORTMANAGER( FUN ) return the catalog as output INSTEAD of updating the
    % primary catalog.
    %
    % ZMAPIMPORTMANAGER( FUN , ARGS ) ARGS is a cell of arguments that will be
    % passed to the import function FUN
    %
    % The ZMAPIMPORTMANAGER exists to do cleanup and shutdown as necessary of other catalogs in
    % memory. this includes the creation of maps and timeplots.
    %
    % catalogs imported throught he ZmapImportManager will be sorted in ascending Date order

    % yes, the argument fun expects the return arguments in reverse order from the ZmapImportManager
    % this is purely practical.
    assert(nargout(fun)==2,'import function must have two output arguments : [catalog, ok]');
    if exist('funArguments','var')
        assert(iscell(funArguments),...
            'Second argument should be a CELL containing arguments to be passed to fun');
        [catalog,ok] = fun(funArguments{:});
    else
        [catalog,ok] = fun();
    end
    if ok
        sort(catalog,'Date','ascend')
    end

    returnTheCatalog = nargout==2;
    if returnTheCatalog
        % do not assume we are modifying the primary catalog. This might be some other load
        if ~ok
            catalog = ZmapCatalog();
        end
    else
        % assume we replace the primary catalog
        if ok
            post_load();
        end
    end

    % save this catalog
    saveFile = fullfile(ZmapGlobal.Data.Directories.working, ZmapGlobal.Data.CatalogOpts.LastCatalogFilename);
    [pathstr,name,ext] = fileparts(saveFile);
    if isempty(ext)
        ext='.mat';
    end
    saveFile=fullfile(pathstr,[name ext]);
    if ~exist(pathstr,'dir')
        mkdir(pathstr);
    end
    try
        save(saveFile, 'catalog');
    catch ME
        warning('ZMAP:unableToSaveCatalog','unable to save the catalog');
        warning(ME.message);
    end

    function post_load()
        disp('post load')
        ZG = ZmapGlobal.Data;
        ZG.primeCatalog = catalog; % since a catalog is a handle, they point to same thing

        % ZG.mainmap_plotby='depth';

        setDefaultValues(catalog);
        ZG.maepi=catalog.subset(catalog.Magnitude > ZG.CatalogOpts.BigEvents.MinMag);

        % OPTIONALLY CLEAR SHAPE
        if ~isempty(ShapeGeneral.ShapeStash)
            % ask whether to keep shape
            switch questdlg('Keep current shape?','SHAPE','yes','no','no')
                case 'no'
                    ShapeGenreal.ShapeStash(ShapeGeneral)
                case 'yes'
                    % do nothing
            end
        end

        % OPTIONALLY CLEAR GRID
        if ~isempty(ZG.Grid)
            switch questdlg('Keep curent grid?','GRID','yes','no','no')
                case 'no'
                    ZG.Grid.delete()
                case 'yes'
                    % do nothing
            end
        end

        uimemorize_catalog();

    end

end

function setDefaultValues(A)
    % SETDEFAULTVALUES sets certain Zmap Global values based on catalog.
    ZG=ZmapGlobal.Data; % get zmap globals

    %  default values
    [t0b, teb] = bounds(A.Date) ;
    ttdif = days(teb - t0b);
    if ~exist('bin_dur','var')
        ZG.bin_dur = days(ceil(ttdif/100));
    elseif ttdif<=10  &&  ttdif>1
        ZG.bin_dur = days(0.1);
    elseif ttdif<=1
        ZG.bin_dur = days(0.01);
    end
    ZG.CatalogOpts.BigEvents.MinMag = max(A.Magnitude) -0.2;
end

author: CelsoReyes