text:
function [cl,cd] = bpolar(BladeState,AC)
% To make your own bpolar:
% BladeState is the state of the blade element with matrix components:
% BladeState =
%            M: Mach number
%        alpha: angle of attack (radians)
%            r: radius location (nondimensional)
%          psi: azimuth location (radians)
%     alphadot: rate of change of angle of attack (rad/s)
%           Re: Reynolds number
% all components on BladeState inputs are m x n, where n is the number of
% blade elements, and m is the uniformly distributed number of azimuth
% elements.
%
% NB: this function should run VERY fast, so avoid loops and extensive
% interpolation.
% see also lininterp1f
%
% AC is a structure with AC data (see documentation)
%
% cl and cd are the m x n arrays containing each blade element's lift and
% drag coefficients.
% note that output coefficients will be non-dimensionalized by a constant
% blade chord, typically the average blade chord.
%
% The blade polar function is responsible for capturing all important rotor
% phenomena, including reverse flow, retreating blade stall, and advancing
% tip critical mach number.
%
% to capture reverse flow effects, the drag polar should be able to handle
% alphas +/- 360 deg
%
% critical mach drag rise and stall should also be modeled.

% This file is an example of using wind tunnel data (NACA 0012), very fast
% interpolation, and a critical mach drag rise approximation.

% Drag Coefficient data
%from http://www.cyberiad.net
% -------------------------------- REYNOLDS NUMBER -----------------------

% RE = [5000000] ;

DragData = [    0.0064
    0.0064
    0.0066
    0.0068
    0.0072
    0.0076
    0.0081
    0.0086
    0.0092
    0.0098
    0.0106
    0.0118
    0.0130
    0.0143
    0.0159
    0.0177
    0.0198
    0.0229
    0.1480
    0.2740
    0.2970
    0.3200
    0.3440
    0.3690
    0.3940
    0.4200
    0.4460
    0.4730
    0.5700
    0.7450
    0.9200
    1.0750
    1.2150
    1.3450
    1.4700
    1.5750
    1.6650
    1.7350
    1.7800
    1.8000
    1.8000
    1.7800
    1.7500
    1.7000
    1.6350
    1.5550
    1.4650
    1.3500
    1.2250
    1.0850
    0.9250
    0.7550
    0.5750
    0.4200
    0.3200
    0.2300
    0.1400
    0.0550
    0.0250
    0.0550
    0.1400
    0.2300
    0.3200
    0.4200
    0.5750
    0.7550
    0.9250
    1.0850
    1.2250
    1.3500
    1.4650
    1.5550
    1.6350
    1.7000
    1.7500
    1.7800
    1.8000
    1.8000
    1.7800
    1.7350
    1.6650
    1.5750
    1.4700
    1.3450
    1.2150
    1.0750
    0.9200
    0.7450
    0.5700
    0.4730
    0.4460
    0.4200
    0.3940
    0.3690
    0.3440
    0.3200
    0.2970
    0.2740
    0.1480
    0.0229
    0.0198
    0.0177
    0.0159
    0.0143
    0.0130
    0.0118
    0.0106
    0.0098
    0.0092
    0.0086
    0.0081
    0.0076
    0.0072
    0.0068
    0.0066
    0.0064
    0.0064];
LiftData =[          0
    0.1100
    0.2200
    0.3300
    0.4400
    0.5500
    0.6600
    0.7700
    0.8800
    0.9900
    1.1000
    1.1842
    1.2673
    1.3242
    1.3423
    1.3093
    1.2195
    1.0365
    0.9054
    0.8412
    0.8233
    0.8327
    0.8563
    0.8903
    0.9295
    0.9718
    1.0193
    1.0680
    0.9150
    1.0200
    1.0750
    1.0850
    1.0400
    0.9650
    0.8750
    0.7650
    0.6500
    0.5150
    0.3700
    0.2200
    0.0700
    -0.0700
    -0.2200
    -0.3700
    -0.5100
    -0.6250
    -0.7350
    -0.8400
    -0.9100
    -0.9450
    -0.9450
    -0.9100
    -0.8500
    -0.7400
    -0.6600
    -0.6750
    -0.8500
    -0.6900
    0
    0.6900
    0.8500
    0.6750
    0.6600
    0.7400
    0.8500
    0.9100
    0.9450
    0.9450
    0.9100
    0.8400
    0.7350
    0.6250
    0.5100
    0.3700
    0.2200
    0.0700
    -0.0700
    -0.2200
    -0.3700
    -0.5150
    -0.6500
    -0.7650
    -0.8750
    -0.9650
    -1.0400
    -1.0850
    -1.0750
    -1.0200
    -0.9150
    -1.0680
    -1.0193
    -0.9718
    -0.9295
    -0.8903
    -0.8563
    -0.8327
    -0.8233
    -0.8412
    -0.9054
    -1.0365
    -1.2195
    -1.3093
    -1.3423
    -1.3242
    -1.2673
    -1.1842
    -1.1000
    -0.9900
    -0.8800
    -0.7700
    -0.6600
    -0.5500
    -0.4400
    -0.3300
    -0.2200
    -0.1100
    0];
Alpha =[0
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    30
    35
    40
    45
    50
    55
    60
    65
    70
    75
    80
    85
    90
    95
    100
    105
    110
    115
    120
    125
    130
    135
    140
    145
    150
    155
    160
    165
    170
    175
    180
    185
    190
    195
    200
    205
    210
    215
    220
    225
    230
    235
    240
    245
    250
    255
    260
    265
    270
    275
    280
    285
    290
    295
    300
    305
    310
    315
    320
    325
    330
    333
    334
    335
    336
    337
    338
    339
    340
    341
    342
    343
    344
    345
    346
    347
    348
    349
    350
    351
    352
    353
    354
    355
    356
    357
    358
    359
    360]*pi/180;
modinpcount = 0;
while any(BladeState.alpha(:)<0)
    ind = BladeState.alpha<0;
    BladeState.alpha(ind) = BladeState.alpha(ind)+2*pi;
    modinpcount=modinpcount+1;
end
if modinpcount>2
    disp('crazy stuff going on in bpolar - very very negative alphas')
end
modinpcount = 0;
while any(BladeState.alpha(:)>2*pi)
    ind = BladeState.alpha>2*pi;
    BladeState.alpha(ind) = BladeState.alpha(ind)-2*pi;
    modinpcount=modinpcount+1;
end
if modinpcount>1
    disp('crazy stuff going on in bpolar - very very positive alphas')
end
cd = mylerp1(Alpha,DragData,BladeState.alpha);%,'linear');
cl = mylerp1(Alpha,LiftData,BladeState.alpha);%,'linear');
if any(isnan([cd(:);cl(:)]));
    disp('    check this out')
end
% Data from http://adg.stanford.edu/aa241/drag/dragrise.html
% at Cl of .3
% M-Mcrit:
Mdiff0 =[ -1
    -0.15
    -0.125
    -0.1
    -0.075
    -0.05
    -0.025
    0
    0.025
    0.05
    0.075
    0.1
    0.125
    0.15
    0.175
    0.2
    0.225
    10000.225];

Cdinvisc =[0
    0
    0
    0.000025
    0.00005
    0.00015
    0.0005
    0.002
    0.0052
    0.0103
    0.0173
    0.026
    0.037
    0.048
    0.068
    0.1
    0.148
    20000];

Mcrit = max(.55,-.02002*abs(BladeState.alpha)+.765);
% from ESDU AERO W.00.03.01
% Critical Mach number for high speed aerofoil sections.

cd = cd + mylerp1(Mdiff0,Cdinvisc,abs(BladeState.M)-Mcrit);

%tip loss
ind = BladeState.r>.97;
cl(ind)=cl(ind).*(1-(BladeState.r(ind)-.97)/.03);

%% root cutout
cl(BladeState.r<.15) = 0;
cd(BladeState.r<.15) = .1;

if any([isnan(cl(:)) isnan(cd(:))])
    disp('NaN in bpolar function')
    cd(isnan(cd)) = 1;
    cl(isnan(cl)) = 1;
end

end

function yi = mylerp1(x,y,xi)
% Y = interp1(x,y,X,'linear');

toreshape = size(xi);
xi = xi(:);
yi = lininterp1f(x,y,xi,NaN); %mex .dll function
yi = reshape(yi,toreshape);

end


text:
function [ a, rank ] = npart_sf_lex_successor ( n, npart, a, rank )

%*****************************************************************************80
%
%% NPART_SF_LEX_SUCCESSOR computes SF NPART partition.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the integer to be partitioned.
%    N must be positive.
%
%    Input, integer NPART, the number of parts of the partition.
%    1 <= NPART <= N.
%
%    Input/output, integer A(NPART), contains the partition.
%    A(1) through A(NPART) contain the nonzero integers which
%    sum to N.  The values in A must be in DESCENDING order.
%
%    Input/output, integer RANK, the rank.
%    If RANK = -1 on input, then the routine understands that this is
%    the first call, and that the user wishes the routine to supply
%    the first element in the ordering, which has RANK = 0.
%    In general, the input value of RANK is increased by 1 for output,
%    unless the very last element of the ordering was input, in which
%    case the output value of RANK is 0.
%

%
%  Return the first element.
%
  if ( rank == -1 )
    a = i4vec_part2 ( n, npart );
    rank = 0;
    return
  end
%
%  Check.
%
  ierror = part_sf_check ( n, npart, a );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'NPART_SF_LEX_SUCCESSOR - Fatal error!\n' );
    fprintf ( 1, '  The input array is illegal.\n' );
    fprintf ( 1, '  IERROR = %d\n', ierror );
    error ( 'NPART_SF_LEX_SUCCESSOR - Fatal error!' );
  end
%
%  Find the last entry that is 2 or more.
%
  for i = npart : -1 : 1
    if ( 1 < a(i) )
      indx = i;
      break
    end
  end
%
%  As long as the last nonunit occurs after the first position,
%  have it donate 1 to the left.
%
  if ( 1 < indx )

    a(indx) = a(indx) - 1;
    a(indx-1) = a(indx-1) + 1;
    indx = indx - 1;

    while ( 1 )

      if ( indx <= 1 )
        break
      end

      if ( a(indx) <= a(indx-1) )
        break
      end

      temp      = a(indx);
      a(indx)   = a(indx-1);
      a(indx-1) = temp;

      indx = indx - 1;

    end
%
%  Sum the tail.
%
    temp = sum ( a(indx+1:npart) );
%
%  Partition the tail sum equally over the tail.
%
    a(indx+1:npart) = i4vec_part2 ( temp, npart - indx );

    rank = rank + 1;
%
%  If A(2) through A(NPART) are 1, then this is the last element.
%  Return the first one.
%
  else

    a = i4vec_part2 ( n, npart );
    rank = 0;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = sgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy )

%*****************************************************************************80
%
%% SGEMV computes y := alpha * A * x + beta * y for general matrix A.
%
%  Discussion:
%
%    SGEMV performs one of the matrix-vector operations
%      y := alpha*A *x + beta*y
%    or
%      y := alpha*A'*x + beta*y,
%    where alpha and beta are scalars, x and y are vectors and A is an
%    m by n matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 February 2014
%
%  Author:
%
%    Jack Dongarra, Jeremy Du Croz, Sven Hammarling,  Richard Hanson.
%
%  Parameters:
%
%    Input, character TRANS, specifies the operation to be performed:
%    'n' or 'N'   y := alpha*A *x + beta*y.
%    't' or 'T'   y := alpha*A'*x + beta*y.
%    'c' or 'C'   y := alpha*A'*x + beta*y.
%
%    Input, integer M, the number of rows of the matrix A.
%    0 <= M.
%
%    Input, integer N, the number of columns of the matrix A.
%    0 <= N.
%
%    Input, real ALPHA, the scalar multiplier for A * x.
%
%    Input, real A(LDA,N).  The M x N subarray contains
%    the matrix A.
%
%    Input, integer LDA, the the first dimension of A as declared
%    in the calling routine.  max ( 1, M ) <= LDA.
%
%    Input, real X(*), an array containing the vector to be
%    multiplied by the matrix A.
%    If TRANS = 'N' or 'n', then X must contain N entries, stored in INCX
%    increments in a space of at least ( 1 + ( N - 1 ) * abs ( INCX ) )
%    locations.
%    Otherwise, X must contain M entries, store in INCX increments
%    in a space of at least ( 1 + ( M - 1 ) * abs ( INCX ) ) locations.
%
%    Input, integer INCX, the increment for the elements of
%    X.  INCX must not be zero.
%
%    Input, real BETA, the scalar multiplier for Y.
%
%    Input/output, real Y(*), an array containing the vector to
%    be scaled and incremented by A*X.
%    If TRANS = 'N' or 'n', then Y must contain M entries, stored in INCY
%    increments in a space of at least ( 1 + ( M - 1 ) * abs ( INCY ) )
%    locations.
%    Otherwise, Y must contain N entries, store in INCY increments
%    in a space of at least ( 1 + ( N - 1 ) * abs ( INCY ) ) locations.
%
%    Input, integer INCY, the increment for the elements of
%    Y.  INCY must not be zero.
%

%
%  Test the input parameters.
%
  info = 0;
  if ( ~ ( trans == 'N' || trans == 'n' || ...
           trans == 'T' || trans == 't' || ...
           trans == 'C' || trans == 'c' ) )
    info = 1;
    fprintf ( 1, 'SGEMV rejects input argument number %d\n', info );
    error ( 'SGEMV - Fatal error!' );
  end

  if ( m < 0 )
    info = 2;
    fprintf ( 1, 'SGEMV rejects input argument number %d\n', info );
    error ( 'SGEMV - Fatal error!' );
  end

  if ( n < 0 )
    info = 3;
    fprintf ( 1, 'SGEMV rejects input argument number %d\n', info );
    error ( 'SGEMV - Fatal error!' );
  end

  if ( lda < max ( 1 : m ) )
    info = 6;
    fprintf ( 1, 'SGEMV rejects input argument number %d\n', info );
    error ( 'SGEMV - Fatal error!' );
  end

  if ( incx == 0 )
    info = 8;
    fprintf ( 1, 'SGEMV rejects input argument number %d\n', info );
    error ( 'SGEMV - Fatal error!' );
  end

  if ( incy == 0 )
    info = 11;
    fprintf ( 1, 'SGEMV rejects input argument number %d\n', info );
    error ( 'SGEMV - Fatal error!' );
  end
%
%  Quick return if possible.
%
  if ( ( m == 0 ) || ...
       ( n == 0 ) || ...
       ( alpha == 0.0 && beta == 1.0 ) )
   return
  end
%
%  Set LENX and LENY, the lengths of the vectors x and y, and set
%  up the start points in X and Y.
%
  if ( trans == 'N' || trans == 'n' )
    lenx = n;
    leny = m;
  else
    lenx = m;
    leny = n;
  end

  if ( 0 < incx )
    kx = 1;
  else
    kx = 1 - ( lenx - 1 ) * incx;
  end

  if ( 0 < incy )
    ky = 1;
  else
    ky = 1 - ( leny - 1 ) * incy;
  end
%
%  Start the operations. In this version the elements of A are
%  accessed sequentially with one pass through A.
%
%  First form  y := beta*y.
%
  if ( beta ~= 1.0 )
    if ( incy == 1 )
      if ( beta == 0.0 )
        y(1:leny) = 0.0;
      else
        y(1:leny) = beta * y(1:leny);
      end
    else
      iy = ky;
      if ( beta == 0.0 )
        for i = 1 : leny
          y(iy) = 0.0;
          iy = iy + incy;
        end
      else
        for i = 1 : leny
          y(iy) = beta * y(iy);
          iy = iy + incy;
        end
      end
    end
  end

  if ( alpha == 0.0 )
    return
  end
%
%  Form y := alpha*A*x + y.
%
  if ( trans == 'N' || trans == 'n' )
    jx = kx;
    if ( incy == 1 )
      for j = 1 : n
        if ( x(jx) ~= 0.0 )
          temp = alpha * x(jx);
          for i = 1 : m
            y(i) = y(i) + temp * a(i,j);
          end
        end
        jx = jx + incx;
      end
    else
      for j = 1 : n
        if ( x(jx) ~= 0.0 )
          temp = alpha * x(jx);
          iy = ky;
          for i = 1 : m
            y(iy) = y(iy) + temp * a(i,j);
            iy = iy + incy;
          end
        end
        jx = jx + incx;
      end
    end
%
%  Form y := alpha*A'*x + y.
%
  else
    jy = ky;
    if ( incx == 1 )
      for j = 1 : n
        temp = 0.0;
        for i = 1 : m
          temp = temp + a(i,j) * x(i);
        end
        y(jy) = y(jy) + alpha * temp;
        jy = jy + incy;
      end
    else
      for j = 1 : n
        temp = 0.0;
        ix = kx;
        for i = 1 : m
          temp = temp + a(i,j) * x(ix);
          ix = ix + incx;
        end
        y(jy) = y(jy) + alpha * temp;
        jy = jy + incy;
      end
    end
  end

  return
end


text:
function hintonDiagram(w,max_m,min_m)
% Display a matrix W where the square's area represents W(i,j)
% and red = negative, green = positive
%
%  Examples
%
%    W = randn(4,5); hintonDiagram(W)

% This file is from pmtk3.googlecode.com

% Based on hintonw function from Mathworks neural net toolbox

if nargin < 1,error('Not enough input arguments.');end
if nargin < 2, max_m = max(max(abs(w))); end
if nargin < 3, min_m = max_m / 100; end
if max_m == min_m, max_m = 1; min_m = 0; end

% DEFINE BOX EDGES
xn1 = [-1 -1 +1]*0.5;
xn2 = [+1 +1 -1]*0.5;
yn1 = [+1 -1 -1]*0.5;
yn2 = [-1 +1 +1]*0.5;

% DEFINE POSITIVE BOX
xn = [-1 -1 +1 +1 -1]*0.5;
yn = [-1 +1 +1 -1 -1]*0.5;

% DEFINE POSITIVE BOX
xp = [xn [-1 +1 +1 +0 +0]*0.5];
yp = [yn [+0 +0 +1 +1 -1]*0.5];

[S,R] = size(w);

cla reset
hold on
set(gca,'xlim',[0 R]+0.5);
set(gca,'ylim',[0 S]+0.5);
set(gca,'xlimmode','manual');
set(gca,'ylimmode','manual');
xticks = get(gca,'xtick');
set(gca,'xtick',xticks(find(xticks == floor(xticks))))
yticks = get(gca,'ytick');
set(gca,'ytick',yticks(find(yticks == floor(yticks))))
set(gca,'ydir','reverse');
if get(0,'screendepth') > 1
  %set(gca,'color',[1 1 1]*.5);
  %set(gcf,'color',[1 1 1]*.3);
end

for i=1:S
  for j=1:R
    m = sqrt((abs(w(i,j))-min_m)/max_m);
    m = min(m,max_m)*0.95;
    if real(m)
      if w(i,j) >= 0
        fill(xn*m+j,yn*m+i,[0 0.8 0])
        plot(xn1*m+j,yn1*m+i,'w',xn2*m+j,yn2*m+i,'k')
      elseif w(i,j) < 0
        fill(xn*m+j,yn*m+i,[0.8 0 0]);
        plot(xn1*m+j,yn1*m+i,'k',xn2*m+j,yn2*m+i,'w');
      end
    end
  end
end

plot([0 R R 0 0]+0.5,[0 0 S S 0]+0.5,'w');
%xlabel('Input');
%ylabel('Neuron');
grid on

end

author: andrewssobral


text:
function [c] = APGLASSOup(b,A,para)

%% Object: c = argmin (1/2)\|y-Dx\|_2^2+lambda\|x\|_1+\mu\|x_I\|_2^2
%%                s.t. x_T >0
%% input arguments:
%%         b -------  D'*y transformed object vector
%%         A -------  D'*D transformed dictionary
%%         para ------ Lambda: Sparsity Level
%%                     (lambda1: template set; lambda2:trivial template; lambda3:\mu)
%%                     Lip: Lipschitz Constant for F(x)
%%                     Maxit: Maximal Iteration number
%%                     nT: number of templates
%% output arguments:
%%         c ------  output Coefficient vetor

%  Initialization
ColDim = size(A,1);
xPrev = zeros(ColDim,1);
x = zeros(ColDim,1);
tPrev = 1;
t = 1;
lambda = para.Lambda;
Lip = para.Lip;
maxit = para.Maxit;
nT = para.nT;

temp_lambda = zeros(ColDim,1);
temp_lambda(1:nT) = lambda(1);

%% main loop
for iter =1:maxit
    tem_t = (tPrev-1)/t;
    tem_y = (1+tem_t)*x - tem_t*xPrev;
    temp_lambda(nT+1:end) = lambda(3)*tem_y(nT+1:end);
    tem_y = tem_y - (A*tem_y-b+temp_lambda)/Lip; % update residual
    xPrev = x;
    x(1:nT) = max(tem_y(1:nT),0);
    x(nT+1:end) = softthres(tem_y(nT+1:end),lambda(2)/Lip);
    tPrev = t;
    t = (1+sqrt(1+4*t^2))/2;
end
c = x;

%% soft thresholding operator
function y = softthres(x,lambda)
y = max(x-lambda,0)-max(-x-lambda,0);

author: yuxng


text:
function [pval, k, K] = circ_kuipertest(alpha1, alpha2, res, vis_on)

% [pval, k, K] = circ_kuipertest(sample1, sample2, res, vis_on)
%
%   The Kuiper two-sample test tests whether the two samples differ
%   significantly.The difference can be in any property, such as mean
%   location and dispersion. It is a circular analogue of the
%   Kolmogorov-Smirnov test.
%
%   H0: The two distributions are identical.
%   HA: The two distributions are different.
%
% Input:
%   alpha1    fist sample (in radians)
%   alpha2    second sample (in radians)
%   res       resolution at which the cdf is evaluated
%   vis_on    display graph
%
% Output:
%   pval        p-value; the smallest of .10, .05, .02, .01, .005, .002,
%               .001, for which the test statistic is still higher
%               than the respective critical value. this is due to
%               the use of tabulated values. if p>.1, pval is set to 1.
%   k           test statistic
%   K           critical value
%
% References:
%   Batschelet, 1980, p. 112
%
% Circular Statistics Toolbox for Matlab

% Update 2012
% By Marc J. Velasco and Philipp Berens, 2009
% velasco@ccs.fau.edu

if nargin < 3
    res = 100;
end
if nargin < 4
    vis_on = 0;
end

n = length(alpha1(:));
m = length(alpha2(:));

% create cdfs of both samples
[phis1 cdf1 phiplot1 cdfplot1] = circ_samplecdf(alpha1, res);
[foo, cdf2 phiplot2 cdfplot2] = circ_samplecdf(alpha2, res); %#ok<ASGLU>

% maximal difference between sample cdfs
[dplus, gdpi] = max([0 cdf1-cdf2]);
[dminus, gdmi] = max([0 cdf2-cdf1]);

% calculate k-statistic
k = n * m * (dplus + dminus);

% find p-value
[pval K] = kuiperlookup(min(n,m),k/sqrt(n*m*(n+m)));
K = K * sqrt(n*m*(n+m));

% visualize
if vis_on
    figure
    plot(phiplot1, cdfplot1, 'b', phiplot2, cdfplot2, 'r');
    hold on
    plot([phis1(gdpi-1), phis1(gdpi-1)], [cdf1(gdpi-1) cdf2(gdpi-1)], 'o:g');
    plot([phis1(gdmi-1), phis1(gdmi-1)], [cdf1(gdmi-1) cdf2(gdmi-1)], 'o:g');
    hold off
    set(gca, 'XLim', [0, 2*pi]);
    set(gca, 'YLim', [0, 1.1]);
    xlabel('Circular Location')
    ylabel('Sample CDF')
    title('CircStat: Kuiper test')
    h = legend('Sample 1', 'Sample 2', 'Location', 'Southeast');
    set(h,'box','off')
    set(gca, 'XTick', pi*(0:.25:2))
    set(gca, 'XTickLabel', {'0', '', '', '', 'pi', '', '', '', '2pi'})
end

end

function [p K] = kuiperlookup(n, k)

load kuipertable.mat;
alpha = [.10, .05, .02, .01, .005, .002, .001];
nn = ktable(:,1);  %#ok<NODEF>

% find correct row of the table
[easy row] = ismember(n, nn);
if ~easy
   % find closest value if no entry is present)
   row = length(nn) - sum(n<nn);
   if row == 0
       error('N too small.');
   else
      warning('N=%d not found in table, using closest N=%d present.',n,nn(row)) %#ok<WNTAG>
   end
end

% find minimal p-value and test-statistic
idx = find(ktable(row,2:end)<k,1,'last');
if ~isempty(idx)
  p = alpha(idx);
else
  p = 1;
end
K = ktable(row,idx+1);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%% levelset2isosurface
% Below is a demonstration of the features of the |levelset2isosurface| function

%%
clear; close all; clc;

%% Syntax
% |[F,V]=levelset2isosurface(K,controlPar);|

%% Description
% This function computes the isosurface for an input levelset image at a
% desired intensity level.
%
% The input consists of:
%%
%
% * A 3D levelset image |K|
% * A structure with control parameters |controlPar| with the fields:
% contourLevel (the level for the isosurface), voxelSize (specifying the
% voxel size), capOpt (specifying how to cap the surface if needed), and
% nSub (specifying if the surface should be drawn at the native resolution
% nSub=[1 1 1] or if it should be courser, e.g. by skippin every 2nd voxels
% if nSub=[2 2 2]).

%% Examples

%% Import image data for this demo

defaultFolder = fileparts(fileparts(mfilename('fullpath'))); %Set main folder
pathNameImageData=fullfile(defaultFolder,'data','DICOM','0001_human_calf');
loadNameImageData=fullfile(pathNameImageData,'IMDAT','IMDAT.mat');
IMDAT_struct=load(loadNameImageData); %The image data structure
G = IMDAT_struct.G; %Geometric/spatial information
v=G.v; %The voxel size
M= IMDAT_struct.type_1; %The image data

pathName=fullfile(defaultFolder,'data','imseg'); %Folder name for contours
contourName='imseg_calf_tibia';

%% Control parameters
pointSpacing=2;
numSmoothSteps=10;

%% Compute levelset

loadName=fullfile(pathName,contourName);
load(loadName); %Load segmentation structure
Vcs=saveStruct.ContourSet; %Access the contour data

[K]=contour2levelset(M,v,Vcs,2);

%%
% Visualize levelset and contours together

%Visualize levelset
vizStruct.colormap=warmcold(250); %Set colormap for levelset visualization
vizStruct.clim=[-abs(min(K(:))) abs(min(K(:)))]; %Set color limits
[~,indMin]=min(K(:)); %assume centre of shape is at lowest value
[i,j,k]=ind2sub(size(M),indMin); %Convert index to subscript indices
vizStruct.sliceIndices=[i,j,k]; %Set indices as slices to plot

hf2=sv3(K,v,vizStruct); %Open slice viewer for levelset

%Visualize contours
optionStruct.Color='k';
plotContours({Vcs},optionStruct);  %Plot contours
drawnow;

%% Compute surface

controlPar.contourLevel=0;
controlPar.voxelSize=v;
controlPar.capOpt=1;
controlPar.nSub=[1 1 1]; %ceil(pointSpacing./v);
[Fi,Vi]=levelset2isosurface(K,controlPar);
%
% %% Smoothen
%
% controlPar_smooth.Method='HC';
% controlPar_smooth.n=numSmoothSteps;
% [Vi]=patchSmooth(Fi,Vi,[],controlPar_smooth);
%
% %% Remesh evenly
%
% controlPar_remesh.pointSpacing=pointSpacing; %Set desired point spacing
% controlPar_remesh.disp_on=1; % Turn off command window text display
% [Fi,Vi]=ggremesh(Fi,Vi,controlPar_remesh);

%% Visualise surface and contours on image

sv3(M,v);
gpatch(Fi,Vi,'bw','none',0.5);
optionStruct.pathName=pathName;
optionStruct.Color='w';
plotContours(contourName,optionStruct);  %Plot contours
axisGeom;
% camlight headlight;
drawnow;

%%
cFigure;
subplot(1,2,1); hold on;
gpatch(Fi,Vi,'bw','k',1);
axisGeom;
camlight headlight;

subplot(1,2,2); hold on;
gpatch(Fi,Vi,'bw','none',0.5);
optionStruct.pathName=pathName;
optionStruct.Color='g';
plotContours(contourName,optionStruct);  %Plot contours
axisGeom;
camlight headlight;

drawnow;

%%

%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*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 test_ft_appendsens

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY

% regular
elec1.elecpos = [1 1 1; 2 2 2; 3 3 3];
elec1.chanpos = elec1.elecpos;
elec1.label = {'1';'2';'3'};
elec1.unit  = 'mm';
elec1.coordsys  = 'acpc';

elec2.elecpos = [4 4 4; 5 5 5];
elec2.chanpos = elec2.elecpos;
elec2.label = {'4';'5'};
elec2.unit  = 'mm';
elec2.coordsys  = 'acpc';

elec = ft_appendsens([], elec1, elec2);

% duplicate chanpos, elecpos, and channel
elec1.elecpos = [1 1 1; 2 2 2; 3 3 3];
elec1.chanpos = elec1.elecpos;
elec1.label = {'1';'2';'3'};
elec1.unit  = 'mm';
elec1.coordsys  = 'acpc';
elec1.tra  = [1 0 0; 0 1 0; 0 0 1];
elec2.elecpos = [2 2 2; 0 0 0]; % 2 2 2 is a duplicate
elec2.chanpos = elec2.elecpos;
elec2.label = {'2';'7'}; % 2 is a duplicate
elec2.unit  = 'mm';
elec2.coordsys  = 'acpc';
elec2.tra  = [1 0; 0 1];
elec = ft_appendsens([], elec1, elec2);

% duplicate elecpos (eg same elec used for two bipolar derivations)
elec1.elecpos = [1 1 1; 2 2 2];
elec1.chanpos = [1.5 1.5 1.5];
elec1.label = {'1-2'};
elec1.unit  = 'mm';
elec1.coordsys  = 'acpc';
elec1.tra  = [1 -1];
elec2.elecpos = [2 2 2; 3 3 3]; % 2 2 2 is a duplicate
elec2.chanpos = [2.5 2.5 2.5];
elec2.label = {'2-3'}; % 2 is a duplicate
elec2.unit  = 'mm';
elec2.coordsys  = 'acpc';
elec2.tra  = [1 -1];

elec = ft_appendsens([], elec1, elec2);

% duplicate chanpos and elecpos, but no duplicate label (eg two electrodes
% localized to the same location)
elec1.elecpos = [1 1 1; 2 2 2; 3 3 3];
elec1.chanpos = elec1.elecpos;
elec1.label = {'1';'2';'3'};
elec1.unit  = 'mm';
elec1.coordsys  = 'acpc';
elec1.tra  = [1 0 0; 0 1 0; 0 0 1];
elec2.elecpos = [2 2 2; 0 0 0]; % 2 2 2 is a duplicate
elec2.chanpos = elec2.elecpos;
elec2.label = {'8';'7'}; % no duplicate label
elec2.unit  = 'mm';
elec2.coordsys  = 'acpc';
elec2.tra  = [1 0; 0 1];

try
  elec = ft_appendsens([], elec1, elec2);
  % this SHOULD throw an error because two labels for the same chanpos are not allowed
  error('an error was expected but not thrown')
catch
  fprintf('ft_appendsens threw an error as expected\n')
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function computeAssociations_Clinical_Old(pathUnivariate,cohorts,clinical,outcomes,timeToEvent)

startpath = pwd;
cd(pathUnivariate)

nameOutcomes = fieldnames(outcomes.(cohorts{1})); nOutcomes = numel(nameOutcomes); nCohort = numel(cohorts);
clinicVar = fieldnames(clinical.(cohorts{1})); nClinic = numel(clinicVar);
rsMat_Binary = zeros(nClinic,nOutcomes); pMat_Binary = zeros(nClinic,nOutcomes); stringCell_Binary = cell(nClinic,nOutcomes);
rsMat_Time = zeros(nClinic,nOutcomes); pMat_Time = zeros(nClinic,nOutcomes); stringCell_Time = cell(nClinic,nOutcomes);
for i = 1:nClinic
    for j = 1:nOutcomes
        var = []; outcome = []; time = [];
        for c = 1:nCohort
            cohort = cohorts{c};
            var = [var;clinical.(cohort).(clinicVar{i})];
            outcome = [outcome;outcomes.(cohort).(nameOutcomes{j})];
            time = [time;timeToEvent.(cohort).(nameOutcomes{j})];
        end
        [rsMat_Binary(i,j),pMat_Binary(i,j)] = corr(var,outcome,'type','Spearman','rows','pairwise');
        [rsMat_Time(i,j),pMat_Time(i,j)] = corr(var,time,'type','Spearman','rows','pairwise');
        if pMat_Binary(i,j) < 0.01
            stringCell_Binary{i,j} = ['rs = ',num2str(rsMat_Binary(i,j),'%.2f'),', p = ',num2str(pMat_Binary(i,j),'%.2i')];
        else
            stringCell_Binary{i,j} = ['rs = ',num2str(rsMat_Binary(i,j),'%.2f'),', p = ',num2str(pMat_Binary(i,j),'%.2f')];
        end
        if pMat_Time(i,j) < 0.01
            stringCell_Time{i,j} = ['rs = ',num2str(rsMat_Time(i,j),'%.2f'),', p = ',num2str(pMat_Time(i,j),'%.2i')];
        else
            stringCell_Time{i,j} = ['rs = ',num2str(rsMat_Time(i,j),'%.2f'),', p = ',num2str(pMat_Time(i,j),'%.2f')];
        end
    end
end
for j = 1:nOutcomes
    pVal = pMat_Binary(:,j);
    [significance] = benjamini_hochberg(pVal,0.10);
    for i = 1:nClinic
        if significance(i)
            stringCell_Binary{i,j} = ['*',stringCell_Binary{i,j},'*']; % Showing significance;
        end
    end
    pVal = pMat_Time(:,j);
    [significance] = benjamini_hochberg(pVal,0.10);
    for i = 1:nClinic
        if significance(i)
            stringCell_Time{i,j} = ['*',stringCell_Time{i,j},'*']; % Showing significance;
        end
    end
end
Locoregional = stringCell_Binary(:,1); Distant = stringCell_Binary(:,2); Death = stringCell_Binary(:,3);
results.tableCorr = table(Locoregional,Distant,Death,'RowNames',clinicVar);
results.rsMat = rsMat_Binary; results.pMat = pMat_Binary;
save('clinical_UniV_Binary','results'), clear results
Locoregional = stringCell_Time(:,1); Distant = stringCell_Time(:,2); Death = stringCell_Time(:,3);
results.tableCorr = table(Locoregional,Distant,Death,'RowNames',clinicVar);
results.rsMat = rsMat_Time; results.pMat = pMat_Time;
save('clinical_UniV_Time','results'), clear results

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
%
% Usage:   Y =mexSort(X);
%
% Name: mexSort
%
% Description: sort the elements of X using quicksort
%
% Inputs: X:  double vector of size n
%
% Output: Y: double  vector of size n
%
% Author: Julien Mairal, 2010

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function ef = cpf_qlim_event(cb_data, cx)
%CPF_QLIM_EVENT  Event function to detect gen reactive power limit violations
%   EF = CPF_QLIM_EVENT(CB_DATA, CX)
%
%   CPF event function to detect generator reactive power limit violations,
%   i.e. Qg <= Qmin or Qg >= Qmax.
%
%   Inputs:
%       CB_DATA : struct of data for callback functions
%       CX : struct containing info about current point (continuation soln)
%
%   Outputs:
%       EF : event function value

%   MATPOWER
%   Copyright (c) 2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   and Shrirang Abhyankar, Argonne National Laboratory
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% event function value is 2 ng x 1 vector equal to:
%%      [ Qg - Qmax ]
%%      [ Qmin - Qg ]

%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;

%% get updated MPC
d = cb_data;
mpc = cpf_current_mpc(d.mpc_base, d.mpc_target, ...
    d.Ybus, d.Yf, d.Yt, d.ref, d.pv, d.pq, cx.V, cx.lam, d.mpopt);

%% compute Qg violations for on-line gens, not at PQ buses
nb = size(mpc.bus, 1);
ng = size(mpc.gen, 1);
on = find(mpc.gen(:, GEN_STATUS) > 0 & ...  %% which generators are on?
          mpc.bus(mpc.gen(:, GEN_BUS), BUS_TYPE) ~= PQ);  %% ... and are not PQ buses
gbus = mpc.gen(on, GEN_BUS);                %% what buses are they at?
ngon = size(on, 1);

%% build connection matrix, element i, j is 1 if gen on(i) at bus j is ON
Cg = sparse((1:ngon)', gbus, ones(ngon, 1), ngon, nb);
C = Cg * Cg';

%% violations are based on total violation at bus, not individual violations
%% (see https://github.com/MATPOWER/matpower/issues/26)
v_Qmax = NaN(ng, 1);
v_Qmin = v_Qmax;
v_Qmax(on) = C * (mpc.gen(on, QG) - mpc.gen(on, QMAX));
v_Qmin(on) = C * (mpc.gen(on, QMIN) - mpc.gen(on, QG));

%% assemble event function value
ef = [v_Qmax; v_Qmin];

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%STATSDTC Stats Decision tree Classifier (Matlab Stats Toolbox)
%
%   W = STATSDTC(A,'PARAM1',val1,'PARAM2',val2,...)
%   W = A*STATSDTC([],'PARAM1',val1,'PARAM2',val2,...)
%   D = B*W
%
% INPUT
%   A          Dataset used for training
%   PARAM1     Optional parameter, see CLASSIFICATIONTREE.FIT
%   B          Dataset used for evaluation
%
% OUTPUT
%   W          Decision tree classifier
%   D          Classification matrix, dataset with posteriors
%
% DESCRIPTION
% This is the PRTools interface to the CLASSIFICATIONTREE of the Matlab
% Stats toolbox. See there for more information. It is assumed that objects
% labels, feature labels and class priors are included in the dataset A.
%
% The decision tree is stored in W and can be retrieved by T = +W or by
% T = getdata(W). The Stats toolbox command VIEW can be used to visualize
% it, either in the command window (default) or graphically setting the
% 'mode' options to 'graph'.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, DTC, TREEC, CLASSIFICATIONTREE, VIEW

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com

function W = statsdtc(varargin)

  name = 'Stats DecTree';

  if mapping_task(varargin,'definition')
    W = define_mapping(varargin,[],name);
  elseif mapping_task(varargin,'training')
    A       = varargin{1};
    data    = +A;
    labels  = getlabels(A);
    prior   = getprior(A);
    featlab = getfeatlab(A);
    if ischar(featlab)
      featlab = cellstr(featlab);
    end
    tree    = ClassificationTree.fit(data,labels,'prior',prior, ...
    'PredictorNames',featlab,varargin{2:end});
    W = trained_mapping(A,tree);
  else % evaluation
    [A,W]    = deal(varargin{:});
    tree     = getdata(W);
    [dummy,post] = predict(tree,+A);
    W        = setdat(A,post,W);
  end

return

author: marianux

repo: ecg-kit

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


text:
function test_tutorial_connectivity3(datadir)

% MEM 6gb
% WALLTIME 00:10:00
% DEPENDENCY ft_timelockanalysis ft_sourceanalysis ft_connectivityanalysis ft_prepare_sourcemodel headsurface

% This is the third section of the connectivity tutorial, which
% starts with the CMC dataset, extracts a virtual channel and performs
% connectivity analysis on the virtual channel time series.

if nargin==0
  % this is where the data should be located
  datadir = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/connectivity');
end

load(fullfile(datadir, 'source.mat'));

[maxval, maxindx] = max(source.avg.coh);
maxpos = source.pos(maxindx,:);

load(fullfile(datadir,'data.mat'));

%% compute the beamformer filter
cfg                   = [];
cfg.covariance        = 'yes';
cfg.channel           = 'MEG';
cfg.vartrllength      = 2;
cfg.covariancewindow  = 'all';
timelock              = ft_timelockanalysis(cfg, data);

cfg                  = [];
cfg.method           = 'lcmv';
cfg.hdmfile          = fullfile(datadir,'SubjectCMC.hdm');
cfg.sourcemodel.pos  = maxpos;
cfg.keepfilter       = 'yes';
source               = ft_sourceanalysis(cfg, timelock);

%% construct the 3-D virtual channel at the location of interest
beamformer = source.avg.filter{1};

chansel = ft_channelselection('MEG', data.label); % find the names
chansel = match_str(data.label, chansel);         % find the indices

sourcedata = [];
sourcedata.label = {'x', 'y', 'z'};
sourcedata.time = data.time;
for i=1:length(data.trial)
  sourcedata.trial{i} = beamformer * data.trial{i}(chansel,:);
end

cfg = [];
cfg.viewmode = 'vertical';  % you can also specify 'butterfly'
ft_databrowser(cfg, sourcedata);

%% construct a single virtual channel in the maximum power orientation
timeseries = cat(2, sourcedata.trial{:});

[u, s, v] = svd(timeseries, 'econ');

% whos u s v
%   Name           Size              Bytes  Class     Attributes
%
%   s              3x3                  72  double
%   u              3x3                  72  double
%   v         196800x3             4723200  double

% this is equal to the first column of matrix V, apart from the scaling with s(1,1)
timeseriesmaxproj = u(:,1)' * timeseries;

virtualchanneldata = [];
virtualchanneldata.label = {'cortex'};
virtualchanneldata.time = data.time;
for i=1:length(data.trial)
  virtualchanneldata.trial{i} = u(:,1)' * beamformer * data.trial{i}(chansel,:);
end

%% combine the virtual channel with the two EMG channels
cfg = [];
cfg.channel = 'EMG';
emgdata = ft_selectdata(cfg, data);

cfg = [];
combineddata = ft_appenddata(cfg, virtualchanneldata, emgdata);

% save combineddata combineddata

%% compute the spectral decomposition
cfg            = [];
cfg.output     = 'fourier';
cfg.method     = 'mtmfft';
cfg.foilim     = [5 100];
cfg.tapsmofrq  = 5;
cfg.keeptrials = 'yes';
cfg.channel    = {'cortex' 'EMGlft' 'EMGrgt'};
freq    = ft_freqanalysis(cfg, combineddata);

cfg = [];
cfg.method = 'coh';
coherence = ft_connectivityanalysis(cfg, freq);

cfg = [];
cfg.zlim = [0 0.2];
figure
ft_connectivityplot(cfg, coherence);
title('coherence')

figure
plot(coherence.freq, squeeze(coherence.cohspctrm(1,2,:)))
title(sprintf('connectivity between %s and %s', coherence.label{1}, coherence.label{2}));
xlabel('freq (Hz)')
ylabel('coherence')


text:
function X = value(X)
% VALUE (overloaded)

X = reshape(value(sdpvar(X)),X.dim);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef antelopesource < dataretrieval.spatiotemporal_database
   %antelopesource Summary of this class goes here
   %   Detailed explanation goes here

   properties
      dbname = ''; %name of antelope database
      isopen = false; % status of antelope database
      searchcriteria % collection of search parameter strings to send to antelope
      chaninfo % array of channelTags
      matlabstart % one or more startdates (in matlab format)
      matlabend   % one or more enddates in matlab format (matching # of startdate)
      dbpointer = []; % antelope database pointer
      filtereddb = []; % antelope database pointer (filtered)
      trpointer = { -1 }; % pointer to antelope trace
      applycalibrations = true; % apply calibration to tr before converting
      combinetraces = true;

   end

   properties(Dependent)
      nrecs  % number or records for current database ptr
      dbfields % field names for current database ptr
      epochstart % start time(s) as epoch number
      epochend % end time(s) as  epoch number
   end

   methods
      % DATE get/set
      function val = get.epochstart(obj)
         val = antelopsource.mat2epoch(obj.matlabstart);
      end
      function val = get.epochend(obj)
         val = antelopesource.mat2epoch(obj.matlabend);
      end
      function obj = set.epochstart(obj, val)
         obj.matlabstart = antelopesource.epoch2mat(val);
      end
      function obj = set.epochend(obj, val)
         obj.matlabend = antelopesource.epoch2mat(val);
      end
      function obj = set.matlabstart(obj, val)
         if ischar(val) || iscell(val)
            val = datnum(val);
         end
         obj.matlabstart = val;
      end
      function obj = set.matlabend(obj, val)
         if ischar(val) || iscell(val)
            val = datnum(val);
         end
         obj.matlabend = val;
      end

      function data = retrieve(obj, where_, from_ , until_)
         assert(numel(datenum(from_)) == numel(datenum(until_)),...
            'Waveform:load_antelope:startEndMismatch', 'Unequal number of start and end times');
         % fill antelope source with proper search criteria
         obj.matlabstart = from_;
         obj.matlabend = until_;
         obj.chaninfo = where_;
         % database to open needs to be determined from all the above
         % do once for each date? each channel?
         obj = obj.open('r');
         % the following needs to be repeated for each where_
         obj = obj.subsetbyLocation(1);
         obj = obj.subsetbyTime(1);
         try
            w = obj.getAllTraces();
         catch
            w = obj.getTracesOneAtATime();
         end

         % return if no rows
         if dbnrecs(obj.dbpointer)==0
            %no data returned
            return
         end
         fprintf('Got %d matching records\n',dbnrecs(db));
         [wfid,sta,chan,st,et]=dbgetv(db, 'wfid','sta','chan','time', 'endtime');
         sta = cellstr(sta);
         chan = cellstr(chan);
         for c=1:numel(sta)
            fprintf('%12d %s %s %f %f\n',wfid(c), sta{c}, chan{c}, st(c), et(c));
         end

         % if starttime and endtime blank, get from min/max times in wfdisc table
         if isempty(starttime) & isempty(endtime)
            starttime = min(st);
            endtime = max(et);
         end
      end

      function p = get.dbfields(obj)
         try
            p = dbquery(obj.dbpointer,'dbTABLE_FIELDS');
         catch
            disp('unable to get fields');
            p = {};
         end
      end

      function tf = isdbfield(obj, candidate)
         %isdbfield   returns true if a string matches a current dbfieldname
         tf = ~isempty(candidate) && ismember(candidate, obj.fields );
      end

      function n = get.nrecs(obj)
         %nrecs   return either the number of records or 0
         try
            n = dbnrecs(obj.dbpointer);
         catch
            disp('unable to get nrecs');
            n = 0;
         end
      end

      function obj = open(obj, attrib)
         %open   open an antelope database to the wfdisc
         %   obj = open(obj, 'r') open antelope database for reading and
         %   sets the talbe to wfdisc
         if obj.isopen()
            error('the database might already be open');
         end
         try
            obj.dbptr = dbopen(obj.dbname, attrib);
            obj.isopen = true;
         catch er
            disp(er)
            error('Unable to open database');
         end
         try
            mydb = dblookup_table(mydb,'wfdisc');
         catch er
            disp(er)
            error('Able to open database, but not to open the wfdisc table');
         end
      end

      function obj = close(obj)
         %close   close an antelope database
         if obj.isopen()
            dbclose(obj.dbptr);
         else
            error('the database might already be closed');
         end
         obj.isopen = false;
      end

      function tf = pointsToAntelopeDatabase(obj)
         %pointsToAntelopeDatabase makes sure pointer points to valid db
         tf = ~isempty(obj.dbpointer) && ...
            isa(obj.dbpointer,'struct') && ...
            any(strcmpi(fieldnames(obj.dbpointer),'database'));
      end

      function obj = subsetbyLocation(obj, n)
         searchstr = antelopesource.buildSearchCriterion(obj.chaninfo(n));
         obj.dbpointer = dbsubset(obj.dbpointer, searchstr);
      end
      function obj = subsetbyTime(obj, n)
         searchstr = sprintf('time <= %f && endtime >= %f',obj.epochend(n), obj.epochstart(n))
         obj.dbpointer = dbsubset(obj.dbpointer, searchstr);
      end

      outputWaveforms = load_antelope(request, specificDatabase)

      function outputWaveforms = RecursivelyLoadFromEachDatabase(request)
         TRY_MULTIDAY = false;
         database = getDatabase(request, TRY_MULTIDAY);
         database = unique(database, 'stable'); %  avoid changing requested order
         outputWaveforms = emptyWaveform();
         for n = 1 : numel(database)
            w = load_antelope(request, database(n));
            outputWaveforms = [outputWaveforms; w(:)]; %#ok<AGROW>
         end
      end

      function database = getDatabase(request, TRY_MULTIDAY)
         [ds, chanInfo, startTimes, endTimes] = unpackDataRequest(request);
         if TRY_MULTIDAY %good luck splicing
            dbDatesToCheck = subdivide_files_by_date(ds,startTimes, endTimes);
            database =  getfilename(ds,chanInfo, dbDatesToCheck);
         else
            database =  getfilename(ds,chanInfo, startTimes);
         end
      end

      function w = cycleThroughTraces(tr, COMBINE_WAVEFORMS)
         %cycleThroughTraces  converts each trace into one or more waveforms
         %   returns a Nx1 waveform
         w(numel(tr)).waves = waveform;
         for traceidx = 1:numel(tr)
            w(traceidx) = wavesfromtraces(tr, traceidx, COMBINE_WAVEFORMS);
         end
         w = [w.waves]; %horizontal cat
         w = w(:);
      end

      function w_scnl = wavesfromtraces(tr, traceidx, COMBINE_WAVEFORMS)
         if ~isstruct(tr{traceidx}) % marker for no data
            w_scnl.waves = emptyWaveform();
         else
            w_scnl.waves = traceToWaveform(tr{traceidx}); %create waveform list
            trdestroy(tr{traceidx});
         end
         if COMBINE_WAVEFORMS %combine all of this trace's records
            w_scnl.waves = combine(w_scnl.waves);
         end;
         w_scnl.waves = reshape(w_scnl.waves,  numel(w_scnl.waves), 1);
      end

      function obj = buildSearchCriteria(obj)
         %buildSearchCriteria creates searchstrings from chaninfo
         %
         % sample output for *.OKTU..EHZ :
         %         {'sta=~/OKTU/ && cha=~/EHZ/ && net=~/.*/'}

         for n=1:numel(obj.chaninfo)
            sc(n) = {buildSearchCriterion(obj.chaninfo(n))};
         end
         sc(cellfun(@isempty,sc))=[]; %remove empty values
         obj.searchcriteria = sc;

         function x = addto(x, fieldname, value)
            %addto   appends ' && fieldname=/value/' to x
            if ~isempty(value)
               if ~isempty(x); x = [x ' && ']; end
               value = replaceStars(value);
               x = [x, fieldname, '=~/', value,'/'];
            end
         end

         function value = replaceStars(value)
            %replaceStars   replace * with .*, but leave existing .* alone
            value = regexprep(value,'(?<!\.)\*','\.\*');
         end
      end

      [tr, rawDb, filteredDb] =  get_antelope_traces(obj, startdates, enddates, criteriaList, database)

   end
   methods(Static, Access=protected)
      [units, data_description] = segtype2units(unitCode)

      function sc = buildSearchCriterion(chaninfo)
         sc = addto('','sta',chaninfo.station);
         sc = addto(sc,'chan',chaninfo.channel);
         % the following aren't understood, and therefore aren't used
         %sc = addto(sc,'net',chaninfo.network);
         %sc = addto(sc,'loc',chaninfo.location);
      end
      function w = emptyWaveform()
         %emptyWaveform   create an empty waveform
         w = waveform;
         w = w([]);
      end

      function w = blankWaveformWithCalib()
         %blankWaveformWithCalib   get blank waveform with default calib
         persistent wBlank
         if numel(wBlank) == 0  % cannot simply test if wBlank is empty
            wBlank = waveform;
            wBlank =  addfield(wBlank,'calib', 0);
            wBlank = addfield(wBlank, 'calibration_applied', 'NO');
         end
         w = wBlank;
      end

      function X = makecell(X)
         %makecell   ensure value is a cell, if not make it
         if ~iscell(X)
            X = {X};
         end
      end

      function Ep = mat2epoch(M)
         %mat2epoch   convert matlab time to epoch time (millisecond precision)
         %   uses Antelope's conversion
         %
         %   if string is converted, then expected format is: 'yyyy-mm-dd HH:MM:SS.FFF'
         %
         %   See also str2epoch, datestr
         if isnumeric(M)
            M = datestr(M,'yyyy-mm-dd HH:MM:SS.FFF');
         end
         Ep = str2epoch(M);
      end
      function M = epoch2mat(Ep)
         %epoch2mat   convert epoch time to matlab time (millisecond precision)
         %   Uses Antelope's conversion
         %
         %  See also epoch2str, datenum
         M = datenum(epoch2str(Ep,'%G %T'));
      end

      function trs = trsTrimStarts(trs, desiredStartEpoch)
         %trsTrimStarts   trims data and adjusts starttime of trStructs
         %   trs = trsTrimStarts(trs, desiredStartEpoch)
         %   adjusts the trs.time, trs.strtime, and trs.data fields
         %
         %   if desiredStartEpoch < trs.time, this does nothing.
         %
         %   See also tr2struct
         for n=1:numel(trs)
            actualStartEpoch = trs(n).time;
            secondsToTrim = desiredStartEpoch - actualStartEpoch;
            if secondsToTrim < 0, continue; end
            samplesToTrim = fix(secondsToTrim * trs(n).samplerate);
            secondsActuallyTrimmed =  samplesToTrim / trs(n).samplerate;
            trs(n).time = actualStartEpoch - secondsActuallyTrimmed;
            trs(n).data(1:samplesToTrim) = []; % remove leading samples
            startString = trs(n).strtime;
            startMat = datenum(startString);
            matSecondsTrimmed = datenum([0,0,0,0,0,secondsActuallyTrimmed]);
            trs(n).strtime = datestr(startMat - matSecondsTrimmed,...
               'mm/dd/yyyy HH:MM:SS.FFF'); %notice date format!
         end

      end

      function trStructs = safe_trload(db, st, et)
         %safe_trload   load a trace (get the struct) and expand time range on failure
         %  trStructs = safe_trload(db, startEpoch, endEpoch)
         %  assumes that db is already subset
         %  returned trace structure may have more data than expected
         earliestStart = min([db.start]);
         tr = [];
         while st >= min([db.start])
            try
               tr = trload_css(mydb, st, et);
            catch
               st = max(st - 60, earliestStart); % back up 60 seconds & try again
            end
         end
         if isempty(tr)
            error('Unable to acess data by merely backing up in time');
         end
         trsplice(tr, 20)
         trStructs = tr2struct(tr);
         trfree(tr);
         trStructs = antelopesource.trsTrimStarts(trStructs, st);
      end

      function T = trstruct2SeismicTrace(trs)
         %trstruct2SeismicTrace  convert all trstructs to SeismicTraces
         %  trstruct2SeismicTrace(trs) where trs is the result of tr2struct,
         %  performed on a tr pointer.
         T = SeismicTrace();
         T = repmat(T,size(trs));
         for n=1:numel(trs)
            T(n).start = trs(n).strtime;
            T(n).data = trs(n).data;
            T(n).station = trs(n).sta;
            T(n).channel = trs(n).chan;
            T(n).network = trs(n).net;
            T(n).samplerate = trs(n).samprate;
            T(n).units = segtype2units(trs(n).segtype);
            T(n).calib.applied = false;
            T(n).calib.value = trs(n).calib;
         end
      end
   end

   end

   %{
   example of a trStuct
   ans =
          evid: -1
          arid: -1
           net: '-'
           sta: 'AKT'
          chan: 'HHZ'
          time: 1.3215e+09
       endtime: 1.3215e+09
     processid: -1
         nsamp: 93399
      samprate: 99.9985
       instype: '-'
         calib: 1.3962
        calper: -1
      response: 0
      datatype: 'sd'
       segtype: 'V'
      procflag: 0
          data: [93399x1 double]
           lat: -999
           lon: -999
          elev: -999
        refsta: '-'
        dnorth: 0
         deast: 0
        edepth: 0
          hang: 0
          vang: 0
       comment: 0
            m0: 0
            m1: 0
           mt0: -1.0000e+10
           mt1: -1.0000e+10
    bundletype: 0
         dbptr: [1x1 struct]
       strtime: '11/16/2011  16:27:50.006'
    strendtime: '11/16/2011  16:43:23.999'

   %}

author: geoscience-community-codes

repo: GISMO


text:
function [dag, best_score, cache] = learn_struct_gs2(data, nodesizes, seeddag, varargin)
%
% LEARN_STRUCT_GS2(data,seeddag) learns a structure of Bayesian net by Greedy Search.
% dag = learn_struct_gs(data, nodesizes, seeddag)
%
% dag: the final structure matrix
% Data : training data, data(i,m) is the m obsevation of node i
% Nodesizes: the size array of different nodes
% seeddag: given seed Dag for hill climbing, optional
% cache : data structure used to memorize local score computations
%   (cf. SCORE_INIT_CACHE function)
%
% by Gang Li @ Deakin University (gli73@hotmail.com)
% (use mk_nbrs_of_dag_topo, developped by Wei Hu, instead of mk_nbrs_of_dag)
% (Caching implementation : ofrancois.olivier.c.h@gmail.com, philippe.leray@univ-nantes.fr)
%

[N ncases] = size(data);
if (nargin < 3 )
    seeddag = zeros(N,N); % mk_rnd_dag(N); %call BNT function
elseif ~acyclic(seeddag)
    seeddag = mk_rnd_dag(N); %zeros(N,N);
end;

% set default params
scoring_fn = 'bic';
verbose  = 'yes';
cache=[];

% get params
args = varargin;
nargs = length(args);
if length(args) > 0
    if isstr(args{1})
    	for i = 1:2:nargs
    		switch args{i}
    		case 'scoring_fn', scoring_fn = args{i+1};
    		case 'verbose',  verbose  = strcmp(args{i+1},'yes');
    		case 'cache',  cache=args{i+1} ;
    		end;
    	end;
    end;
end;

done = 0;
[best_score cache] = score_dags(data,nodesizes, {seeddag},'scoring_fn',scoring_fn,'cache',cache);
while ~done
    [dags,op,nodes] = mk_nbrs_of_dag_topo(seeddag);
    nbrs = length(dags);
    [scores cache] = score_dags(data, nodesizes, dags,'scoring_fn',scoring_fn,'cache',cache);
    max_score = max(scores);
    new = find(scores == max_score );
    if ~isempty(new) & (max_score > best_score)
        p = sample_discrete(normalise(ones(1, length(new))));
        best_score = max_score;
        seeddag = dags{new(p)};
    else
        done = 1;
    end;
end;

dag = seeddag;

outcount = 0;
[best_score cache] = score_dags(data,nodesizes, {seeddag},'scoring_fn',scoring_fn,'cache',cache);
while outcount < 2
    innercount = 0;
    for i=1:N
        for j=1:N
           if i==j, continue;    end;
           if seeddag(i,j) == 0  % No edge i-->j, then try to add it
               tempdag = seeddag;
               tempdag(i,j) = 1;
               if acyclic(tempdag)
                    [temp_score cache] = score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn,'cache',cache);
                    if temp_score > best_score
                        seeddag = tempdag
                        best_score= temp_score;
                        innercount = innercount +1;
                    end;
               end
           else  % exists edge i--j, then try reverse it or remove it
               tempdag = seeddag;
               tempdag(i,j) = 0; tempdag(j,i) = 1;
               if acyclic(tempdag)
                   [temp_score cache] = score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn,'cache',cache);
                   if temp_score > best_score
                       seeddag = tempdag;
                       best_score = temp_score;
                       innercount = innercount +1;
                   else
                       tempdag = seeddag;
                       tempdag(i,j) = 0;
                       [temp_score cache] = score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn,'cache',cache);
                       if temp_score > best_score
                           seeddag = tempdag;
                           best_score= temp_score;
                           innercount = innercount +1;
                       end;
                   end;
               else
                   tempdag = seeddag;
                   tempdag(i,j)=0;
                   [temp_score cache] = score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn,'cache',cache);
                   if temp_score > best_score
                       seeddag = tempdag
                       best_score= temp_score;
                       innercount = innercount +1;
                   end;
               end;
           end;
        end; % end for j
    end; % end for i
    if innercount == 0
        outcount = outcount +1;
    end;
end;  % end while

dag = seeddag;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef DifferentialDrive < simiam.robot.dynamics.Dynamics

% Copyright (C) 2013, Georgia Tech Research Corporation
% see the LICENSE file included with this software

    properties
        wheel_radius
        wheel_base_length
    end

    methods
        function obj = DifferentialDrive(wheel_radius, wheel_base_length)
            obj = obj@simiam.robot.dynamics.Dynamics();
            obj.wheel_radius = wheel_radius;
            obj.wheel_base_length = wheel_base_length;
        end

        function pose_t_1 = apply_dynamics(obj, pose_t, dt, vel_r, vel_l)
            R = obj.wheel_radius;
            L = obj.wheel_base_length;

%             fprintf('(vel_r,vel_l) = (%0.6g,%0.6g)\n', vel_r,vel_l);

            v = R/2*(vel_r+vel_l);
            w = R/L*(vel_r-vel_l);

%             fprintf('Calculated velocities (v,w): (%0.3g,%0.3g)\n', v, w);

            [x_k, y_k, theta_k] = pose_t.unpack();

            options = odeset('RelTol',1e-8,'AbsTol',1e-8);
            [t,z] = ode23(@obj.dynamics, [0 dt], [x_k, y_k, theta_k, v, w], options);

%             x_k_1 = x_k + dt*(v*cos(theta_k));
%             y_k_1 = y_k + dt*(v*sin(theta_k));
%             theta_k_1 = theta_k + dt*w;

%             pose_t_1 = simiam.ui.Pose2D(x_k_1, y_k_1, theta_k_1);
            pose_t_1 = simiam.ui.Pose2D(z(end,1),z(end,2),z(end,3));
        end

        function dz = dynamics(obj, t, z)
            dz = zeros(5,1);
            dz(1:2) = z(4)*[cos(z(3));sin(z(3))];
            dz(3) = z(5);
        end

        function [vel_r,vel_l] = uni_to_diff(obj,v,w)
            R = obj.wheel_radius;
            L = obj.wheel_base_length;

            %% START CODE BLOCK %%
            vel_r = (2*v+w*L)/(2*R);
            vel_l = (2*v-w*L)/(2*R);
            %% END CODE BLOCK %%
        end

        function [v,w] = diff_to_uni(obj,r,l)
            R = obj.wheel_radius;
            L = obj.wheel_base_length;

            v = R/2*(r+l);
            w = R/L*(r-l);
        end
    end

end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
function [ y, ierror ] = y_check_bahai ( y )

%*****************************************************************************80
%
%% Y_CHECK_BAHAI checks a Bahai year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, the year, which must not be 0.
%
%    Output, integer IERROR, is 0 if Y is legal, and 1 otherwise.
%
  if ( y <= 0 )
    ierror = 1;
  else
    ierror = 0;
  end

  return
end


text:
function y = randomwalk(x0,p,nsteps)

y    = zeros(1,nsteps);
y(1) = x0;
for istep = 2:nsteps
      y(istep) = y(istep-1) + step(p);
end

function y = step(p)

y = rand;
if ( rand < 1-p )
      y = -y; % variable step
end

% %% Generate and plot a random walk
% x = randomwalk(0,0.5,10000);
% plot(x)

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
%
% function masks = cands2labels(candidates,ms)
%
% From a hierarchy (lp: leaves partition and ms: merging sequence) and a
% list of candidates as sets of labels (candidates), we get each candidate
% as a list of labels of the leaves partition lp.
% ------------------------------------------------------------------------
function labels = cands2labels(candidates,ms)
    labels = mex_cands2labels(double(ms),candidates);
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function cond = bodewig_condition ( )

%*****************************************************************************80
%
%% BODEWIG_CONDITION returns the L1 condition of the BODEWIG matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real COND, the L1 condition.
%
  cond = 10.436619718309862;

  return
end


text:
function K = covMaterniso(d, hyp, x, z, i)

% Matern covariance function with nu = d/2 and isotropic distance measure. For
% d=1 the function is also known as the exponential covariance function or the
% Ornstein-Uhlenbeck covariance in 1d. The covariance function is:
%
%   k(x^p,x^q) = sf^2 * f( sqrt(d)*r ) * exp(-sqrt(d)*r)
%
% with f(t)=1 for d=1, f(t)=1+t for d=3 and f(t)=1+t+t²/3 for d=5.
% Here r is the distance sqrt((x^p-x^q)'*inv(P)*(x^p-x^q)), P is ell times
% the unit matrix and sf2 is the signal variance. The hyperparameters are:
%
% hyp = [ log(ell)
%         log(sf) ]
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2010-09-10.
%
% See also COVFUNCTIONS.M.

if nargin<3, K = '2'; return; end                  % report number of parameters
if nargin<4, z = []; end                                   % make sure, z exists
xeqz = isempty(z); dg = strcmp(z,'diag');                       % determine mode

ell = exp(hyp(1));
sf2 = exp(2*hyp(2));
if all(d~=[1,3,5]), error('only 1, 3 and 5 allowed for d'), end         % degree

switch d
  case 1, f = @(t) 1;               df = @(t) 1;            % df(t) = f(t)-f'(t)
  case 3, f = @(t) 1 + t;           df = @(t) t;
  case 5, f = @(t) 1 + t.*(1+t/3);  df = @(t) t.*(1+t)/3;
end
          m = @(t,f) f(t).*exp(-t); dm = @(t,f) df(t).*exp(-t).*t;

% precompute distances
if dg                                                               % vector kxx
  K = zeros(size(x,1),1);
else
  if xeqz                                                 % symmetric matrix Kxx
    K = sq_dist(sqrt(d)/ell*x');
  else                                                   % cross covariances Kxz
    K = sq_dist(sqrt(d)/ell*x',sqrt(d)/ell*z');
  end
end

if nargin<5                                                        % covariances
  K = sf2*m(sqrt(K),f);
else                                                               % derivatives
  if i==1
    K = sf2*dm(sqrt(K),f);
  elseif i==2
    K = 2*sf2*m(sqrt(K),f);
  else
    error('Unknown hyperparameter')
  end
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function E = bitreslice(A,b_old, b_new, echo_on)
% function E = bitreslice(A,b_old, b_new, echo_on)
% Author: Damon Bradley
% Purpose: A is an [M by N] matrix of integers. Implicit is each element's
% representaiton in binary, which is b_old bits wide. Matrix E is output
% which is an [M by L] matrix of integers, with the bits of A concatenated
% column-wise and re-split into b_new bit-wide elements. The number of bits
% per row, b_old*N must be evenly divisible by b_new to work.
%
% Inputs:
% A         - An M by N matrix of unsigned integers
% b_old     - Number of bits to represent A, should be based on maximum value
% b_new     - New number of bits to slice each bit-concatenated row of A
% echo_on   - Choose 1 to display underlying bit representation. Don't specify otherwise

if nargin > 3,
    lookunderhood = true;
else
    lookunderhood = false;
end
[M N] = size(A);
TotalBitsPerRow = N*b_old;
L = TotalBitsPerRow/b_new;

%A = [3 4 7; 3 2 0; 10 5 5]
[MaxA pos] = max(A(:));
required_old_bits = ceil(log2(MaxA));
if b_old < required_old_bits,
    [i j] = ind2sub(size(A),pos);
    str = ['A(' num2str(i) ',' num2str(j) ')=' num2str(MaxA)];
    error(['D''oh! ' num2str(b_old) ...
        ' bits are not enough to represent largest element, (' str '), in the input matrix.' ])
end

if lookunderhood,
    B = dec2binMx(A,b_old)
    C = cell2mat(B);
    D = mat2cell(C,ones(1,M), b_new*ones(1,L))
    E = bin2decMx(D);
else
    B = dec2binMx(A,b_old);
    C = cell2mat(B);
    D = mat2cell(C,ones(1,M), b_new*ones(1,L));
    E = bin2decMx(D);
end

end

function B = dec2binMx(A, bits)
% function y = dec2binMx(x, bits)
% Author: Damon Bradley
% Purpose: perform dec2bin but on matrices. If size(A) = [M N], then B is a
% cell array of string representations of the b-bit representation of the entries of A, sized [M N].

M = size(A,1);
Z = cellstr(dec2bin(A(:),bits));
B = reshape(Z,M,[]);

end

function B = bin2decMx(A)
% function B = bin2decMx(A, bits)
% Author: Damon Bradley
% Purpose: perform bin2dec on cell array obinary entries.

B = cellfun(@bin2dec2,A);

end

function dec = bin2dec2(binstring)
    bitvec = double(binstring) == 49;
    dec = polyval(bitvec, 2);
end

author: Sable

repo: mcbench-benchmarks


text:
function init(obj, varargin)
% init(obj, varargin)
%
% Initialize the CalvinNN network with the default options.
% Also starts up GPU pool.
%
% Copyright by Holger Caesar, 2016

defnnOpts.expDir = fullfile('data', 'exp');
defnnOpts.continue = true;
defnnOpts.batchSize = 2;
defnnOpts.numSubBatches = 2;
defnnOpts.gpus = [];
defnnOpts.numEpochs = 16;
defnnOpts.learningRate = [repmat(1e-3, [1, 12]), repmat(1e-4, [1, 4])];
defnnOpts.weightDecay = 0.0005;
defnnOpts.momentum = 0.9;
defnnOpts.derOutputs = {'objective', 1};
defnnOpts.lossFnObjective = 'softmaxlog'; % Default is softmax
defnnOpts.extractStatsFn = @CalvinNN.extractStats;
defnnOpts.testFn = @(imdb, nnOpts, net, inputs, batchInds) error('Error: Test function not implemented'); % function used at test time to evaluate performance
defnnOpts.misc = struct(); % fields used by custom layers are stored here
defnnOpts.plotEval = true;
defnnOpts.plotAccuracy = true;

% Network options
defnnOpts.convertToTrain = true;

% Fast R-CNN options
defnnOpts.fastRcnn = true;
defnnOpts.fastRcnnParams = true; % learning rates and weight decay
defnnOpts.misc.roiPool.use = true;
defnnOpts.misc.roiPool.freeform.use = false;
defnnOpts.bboxRegress = true;

% Merge input settings with default settings
nnOpts = vl_argparse_old(defnnOpts, varargin, 'nonrecursive');

% Check settings
assert(numel(nnOpts.learningRate) == 1 || numel(nnOpts.learningRate) == nnOpts.numEpochs);

% Do not create directory in evaluation mode
if ~exist(nnOpts.expDir, 'dir') && ~isempty(nnOpts.expDir),
    mkdir(nnOpts.expDir);
end

% Setup GPU
numGpus = numel(nnOpts.gpus);
assert(numGpus <= 1);
if numGpus == 1,
    gpuDevice(nnOpts.gpus);
end

% Set new fields
obj.nnOpts = nnOpts;

author: nightrome

repo: matconvnet-calvin


text:
function hits = buffon_simulate ( a, l, trial_num )

%*****************************************************************************80
%
%% BUFFON_SIMULATE simulates a Buffon needle experiment.
%
%  Discussion:
%
%    In the Buffon needle experiment, we suppose that the plane has been
%    ruled by vertical lines with a spacing of A units, and that a
%    needle of length L is dropped "at random" onto this grid.
%
%    Because of the various symmetries, we may assume that this eye of
%    this needle lands in the first infinite strip, and we may further
%    assume that its Y coordinate is 0.  Thus, we have
%    the eye as (X1,Y1) with 0 <= X1 <= A and Y1 = 0.
%
%    ANGLE, the angle that the needle makes is taken to be uniformly random.
%    The point of the needle, (X2,Y2), therefore lies at
%
%      (X2,Y2) = ( X1+L*cos(ANGLE), Y1+L*sin(ANGLE) )
%
%    The needle will have crossed at least one grid line if any of the
%    following are true:
%
%      X2 <= 0, A <= X2.
%
%    The probability of a crossing on a single trial is
%
%      P(A,L) = ( 2 * L ) / ( PI * A )
%
%    and therefore, a record of the number of hits for a given number of
%    trials can be used as a very roundabout way of estimating PI.
%
%    Note that this routine will try to generate 4 * TRIAL_NUM random
%    double precision values at one time, using automatic arrays.
%    When I tried this with TRIAL_NUM = 1,000,000, the program failed,
%    because of internal system limits on such arrays.
%
%    Such a problem could be avoided by using a DO loop running through
%    each trial individually, but this tend to run much more slowly than
%    necessary.
%
%    Since this routine invokes the FORTRAN90 random number generator,
%    the user should initialize the random number generator, particularly
%    if it is desired to control whether the sequence is to be varied
%    or repeated.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, the horizontal spacing between the
%    vertical grid lines.  0 <= A.
%
%    Input, real L, the length of the needle.
%
%    Input, integer TRIAL_NUM, the number of times the needle is
%    to be dropped onto the grid.
%
%    Output, integer BUFFON_SIMULATE, the number of times the needle
%    crossed at least one line of the grid of cells.
%
%  Local Parameters:
%
%    Local, integer BATCH_SIZE, specifies the number of trials to be done
%    in a single batch.  Setting BATCH_SIZE to 1 will be very slow.
%    Replacing it by TRIAL_NUM would be fine except that your system
%    may have a limit on the size of automatic arrays.  We have set a default
%    value of 10,000 here which should be large enough to be efficient
%    but small enough not to annoy the system.
%
  batch_size = 10000;

  hits = 0;

  for batch = 1 : batch_size: trial_num

    n = min ( batch_size, trial_num + 1 - batch );
%
%  Randomly choose the location (X1,Y1) of the eye of the needle
%  in [0,0]x[A,0], and the angle the needle makes.
%
    x1(1:n) = a * rand ( 1, n );
    angle(1:n) = 2.0D+00 * pi * rand ( 1, n );
%
%  Compute the location of the point of the needle.
%  We only need to know the value of X2, not Y2!
%
    x2(1:n) = x1(1:n) + l * cos ( angle(1:n) );
%
%  Count the end locations that lie outside the cell.
%
    hits = hits + length ( find ( x2(1:n) <= 0.0 | a <= x2(1:n) ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function llsq_test ( )

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

  llsq_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LLSQ_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function gui_l (x, y)
%
% This file is needed by FASTICAG

% The load dialog for loading new data
% and new initial guess.

% @(#)$Id$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global variables

% Handle to the window
global hf_FastICA_Load;

% Handles to some of the controls in window
global he_FastICA_file;

% What is the load type of load dialog
global g_FastICA_loadType;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Configuration options
FIGURENAME = 'FastICA: Load';
FIGURETAG = 'f_FastICALoad';
FIGURESIZE = [x y 450 150];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check to see if this figure is already open - it should not!
% Can't have more than one copy - otherwise the global
% variables and handles can get mixed up.
if ~isempty(findobj('Tag',FIGURETAG))
  error('Error: load dialog already open!');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize some of the controls' values

% What are we loading - who is calling?
caller = get(gcf, 'CurrentObject');

switch get(caller, 'Tag')
 case 'b_LoadData'                             % Do we load new data...
  loadString = 'Load data from variable in Matlab.';
  g_FastICA_loadType = 'data';
  FIGURENAME = 'FastICA: Load data';

 case 'b_LoadGuess'                            % ... or new initial guess?
  loadString = 'Load initial guess for mixing matrix A from variable in Matlab.';
  g_FastICA_loadType = 'guess';
  FIGURENAME = 'FastICA: Load initial guess';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the figure
a = figure('Color',[0.8 0.8 0.8], ...
	   'PaperType','a4letter', ...
	   'Name', FIGURENAME, ...
	   'NumberTitle', 'off', ...
	   'Tag', FIGURETAG, ...
	   'Position', FIGURESIZE, ...
	   'MenuBar', 'none');
set (a, 'Resize', 'off');

hf_FastICA_Load = a;

set(hf_FastICA_Load, 'HandleVisibility', 'callback');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% From here on it get's ugly as I have not had time to clean it up

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the frames
pos_l=2;
pos_w=FIGURESIZE(3)-4;
pos_h=FIGURESIZE(4)-4;
pos_t=2;
h_f_load_background = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'Style','frame', ...
  'Tag','f_load_background');

pos_w=120;
pos_l=FIGURESIZE(3)-(pos_w+2+2);
pos_h=FIGURESIZE(4)-2*4;
pos_t=4;
h_f_load_side = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'Style','frame', ...
  'Tag','f_load_side');

pos_l=4;
pos_w=FIGURESIZE(3)-8-pos_w-2;
pos_h=FIGURESIZE(4)-8;
pos_t=4;
h_f_load = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'Style','frame', ...
  'Tag','f_load');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Controls in f_load
bgc = get(h_f_load, 'BackgroundColor');

pos_w=230;

pos_frame=get(h_f_load, 'Position');
pos_h = 40;
pos_t = pos_frame(2) + pos_frame(4) - pos_h - 6;
pos_l = pos_frame(1) + 6;

b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String',loadString, ...
  'Style','text', ...
  'Tag','t_93');

pos_h = 20;
pos_t = pos_t - pos_h - 10;
pos_l = pos_frame(1) + 6;

b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Name of the variable:', ...
  'Style','text', ...
  'Tag','t_92');

pos_w = 200;
pos_l = pos_l + 30;
pos_t = pos_t - pos_h;
he_FastICA_file = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','', ...
  'Style','edit', ...
  'Tag','e_file');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Controls in f_load_side
pos_vspace = 6;
pos_hspace = 10;
pos_frame = get(h_f_load_side, 'Position');
pos_w = 100;
pos_h = 30;
pos_l = pos_frame(1) + pos_hspace;
pos_t = pos_frame(2) + pos_frame(4) - pos_h - pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_lc Load', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Load', ...
  'Tag','b_lLoad');

pos_t=pos_t-pos_h-pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_lc Cancel', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Cancel', ...
  'Tag','b_lCancel');

pos_t = pos_frame(2) + pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_lc Help', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Help', ...
  'Tag','b_lHelp');


text:
function sub = block(blocks, block_sizes)
% BLOCK Return a vector of subscripts corresponding to the specified blocks.
% sub = block(blocks, block_sizes)
%
% e.g., block([2 5], [2 1 2 1 2]) = [3 7 8].

blocks = blocks(:)';
block_sizes = block_sizes(:)';
skip = [0 cumsum(block_sizes)];
start = skip(blocks)+1;
fin = start + block_sizes(blocks) - 1;
sub = [];
for j=1:length(blocks)
  sub = [sub start(j):fin(j)];
end


text:
function [AF, u, v, w] = arrayFactor(xPos, yPos, zPos, elementWeights, f, c, thetaScanAngles, phiScanAngles, thetaSteerAngle, phiSteerAngle)
%arrayFactor - Calculate array factor of 1D, 2D or 3D array
%
%This matlab function calculates the array factor of a 1D, 2D or 3D array based
%on the position of the elements/sensors and the weight associated with
%each sensor. If no angle is given as input, the scanning angle is theta
%from -90 to 90, and phi from 0 to 360 degrees with 1 degree resolution
%
%[AF, u, v, w] = arrayFactor(xPos, yPos, zPos, elementWeights, f, c, thetaScanAngles, phiScanAngles, thetaSteerAngle, phiSteerAngle)
%
%IN
%xPos            - 1xP vector of x-positions
%yPos            - 1xP vector of y-positions
%zPos            - 1xP vector of z-positions
%elementWeights  - 1xP vector of element weights
%f               - Wave frequency
%c               - Speed of sound
%thetaScanAngles - 1xM vector or MxN matrix of theta scanning angles in degrees (optional)
%phiScanAngles   - 1XN vector or MxN matrix of phi scanning angles in degrees (optional)
%thetaSteerAngle - Theta steering angle in degrees (optional)
%phiSteerAngle   - Phi steering angle in degrees (optional)
%
%OUT
%AF              - Calculated array factor
%u               - MxN matrix of u coordinates in UV space [sin(theta)*cos(phi)]
%v               - MxN matrix of v coordinates in UV space [sin(theta)*sin(phi)]
%w               - MxN matrix of w coordinates in UV space [cos(theta)]
%
%
%Created by J?rgen Grythe
%Last updated 2017-02-27

if ~isvector(xPos)
    error('X-positions of array elements must be a 1xP vector where P is number of elements')
end

if ~isvector(yPos)
    error('Y-positions of array elements must be a 1xP vector where P is number of elements')
end

if ~isvector(elementWeights)
    error('Weighting of array elements must be a 1xP vector where P is number of elements')
end

if ~isscalar(f)
    error('The input frequency must be a single value')
end

%theta is the elevation and is the normal incidence angle from -90 to 90
if ~exist('thetaScanAngles', 'var')
    thetaScanAngles = -pi/2:pi/180:pi/2;
else
    thetaScanAngles = thetaScanAngles*pi/180;
end

%phi is the azimuth, and is the angle in the XY-plane from 0 to 360
if ~exist('phiScanAngles', 'var')
    phiScanAngles = 0:pi/180:2*pi;
else
    phiScanAngles = phiScanAngles*pi/180;
end

%theta, phi steering angles
if ~exist('thetaSteerAngle', 'var')
    thetaSteerAngle = 0;
else
    thetaSteerAngle = thetaSteerAngle*pi/180;
end

if ~exist('phiSteerAngle', 'var')
    phiSteerAngle = 0;
else
    phiSteerAngle = phiSteerAngle*pi/180;
end

%Wavenumber
k = 2*pi*f/c;

%Number of elements/sensors in the array
P = length(xPos);

%Calculating wave vector in spherical coordinates
if isvector(thetaScanAngles)

    %Size of vectors containing theta and phi angles
    M = length(thetaScanAngles);
    N = length(phiScanAngles);

    %Calculate UV coordinates
    u = sin(thetaScanAngles)'*cos(phiScanAngles);
    v = sin(thetaScanAngles)'*sin(phiScanAngles);
    w = repmat(cos(thetaScanAngles)', 1, N);

    % Apply steering
    us = u - sin(thetaSteerAngle)*cos(phiSteerAngle);
    vs = v - sin(thetaSteerAngle)*sin(phiSteerAngle);
    ws = w - cos(thetaSteerAngle);
else

    %Size of matrix containing theta and phi angles
    [M, N] = size(thetaScanAngles);

    %Calculate UV coordinates
    u = sin(thetaScanAngles).*cos(phiScanAngles);
    v = sin(thetaScanAngles).*sin(phiScanAngles);
    w = cos(thetaScanAngles);

    % Apply steering
    us = u - sin(thetaSteerAngle).*cos(phiSteerAngle);
    vs = v - sin(thetaSteerAngle).*sin(phiSteerAngle);
    ws = w - cos(thetaSteerAngle);
end

%Calculate array factor
uu = bsxfun(@times, us, reshape(xPos, 1, 1, P));
vv = bsxfun(@times, vs, reshape(yPos, 1, 1, P));
ww = bsxfun(@times, ws, reshape(zPos, 1, 1, P));

g = repmat(reshape(elementWeights, 1, 1, P), M, N);

AF = sum(g.*exp(1j*k*(uu + vv + ww)), 3);

%Normalising
AF = abs(AF)./max(max(abs(AF)));

%
%                 N
%AF(theta, phi) = sum [ g_n * exp{jk(u*x_n + v*y_n + w*z_n)} ]
%                n=1
%
%u =
%|sin(theta_0)*cos(phi_0) sin(theta_0)*cos(phi_1) .. sin(theta_0)*cos(phi_N)|
%|sin(theta_1)*cos(phi_0) sin(theta_1)*cos(phi_1) .. sin(theta_1)*cos(phi_N)|
%|    .                           .                         .               |
%|sin(theta_M)*cos(phi_0) sin(theta_M)*cos(phi_1) .. sin(theta_M)*cos(phi_N)|

%v =
%|sin(theta_0)*sin(phi_0) sin(theta_0)*sin(phi_1) .. sin(theta_0)*sin(phi_N)|
%|sin(theta_1)*sin(phi_0) sin(theta_1)*sin(phi_1) .. sin(theta_1)*sin(phi_N)|
%|    .                           .                         .               |
%|sin(theta_M)*sin(phi_0) sin(theta_M)*sin(phi_1) .. sin(theta_M)*sin(phi_N)|

%w =
%|cos(theta_0) cos(theta_0) .. cos(theta_0)|
%|cos(theta_1) cos(theta_1) .. cos(theta_1)|
%|    .            .                .      |
%|cos(theta_N) cos(theta_N) .. cos(theta_N)|

%uu =
%   --------
%  /       /|
% / xPos  / |
%---------  | M (length theta)
%|       |  |
%|   u   |  /
%|       | / P (# elements)
%---------/
%   N (length phi)

%g =
%   --------
%  / g_P   /|
% /       / |
%---------  | M
%|g1   g1|  |
%|   g1  |  /
%|g1   g1| / P (# elements)
%---------/
%   N

author: jorgengrythe

repo: beamforming


text:
function b = is3DImage(img)
%IS3DIMAGE  Check if an image is 3D
%
%   B = isColorImage(IMG);
%   Returns TRUE if the image IMG is 3D. An image is assumed to be 3D if
%   all conditions are satisfied:
%   - number of dimensions is >= 3
%   - size of third dimension is not equal to 3, or all sizes equal 3
%
%   Example
%   % planar image
%   is3DImage(imread('cameraman.tif'))
%   ans =
%       0
%
%   % Three-D image
%   is3DImage(ones(5, 5, 5))
%   ans =
%       1
%
%   % planar color image
%   is3DImage(ones(5, 5, 3))
%   ans =
%       0
%
%   Three-D color image
%   is3DImage(ones(5, 5, 3, 4));
%   ans =
%       1
%
%   See also
%     isColorImage
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-05-20,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% check number of dimension
dim = size(img);
if length(dim) < 3
    b = false;
    return
end

% check all dimensions equal to 3
if sum(dim([1:2 4:end])~=3) == 0
    b = true;
    return
end

% check third dimension
if dim(3) == 3
    dim(3) = [];
end

% 3D if 3 dimensions...
b = length(dim) == 3;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function test_headmodel_dipoli_new_old

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_headmodel_dipoli mesh_sphere ft_prepare_vol_sens ft_compute_leadfield

% tests the old and new version of dipoli wrapper

% generate a unit sphere
[pnt, tri] = mesh_sphere(162);

% create the BEM geometries
geom = [];
geom.bnd(1).pnt = pnt * 100;
geom.bnd(1).tri = tri;
geom.bnd(2).pnt = pnt * 90;
geom.bnd(2).tri = tri;
geom.bnd(3).pnt = pnt * 80;
geom.bnd(3).tri = tri;

elec.chanpos = pnt * 100;
elec.elecpos = pnt * 100;
for i=1:size(pnt,1)
  elec.label{i} = sprintf('%d', i);
end

arg(1).name = 'conductivity';
arg(1).value = {[], [1 1/20 1], [0.33 0.125 0.33], [1 1 1], [0.1 0.1 0.1]};

arg(2).name = 'isolatedsource';
arg(2).value = {'yes' , 'no'};

optarg = constructalloptions(arg);

for i=1:size(optarg,1)

  vol = {};
  arg = optarg(i,:);

  % new way - low level:
  vol{1} = ft_headmodel_dipoli(geom.bnd,arg{:});

  % old way:
  tmpcfg = ft_keyval2cfg(arg{:});
  tmpcfg.method = 'dipoli';
  vol{2} = ft_prepare_bemmodel(tmpcfg,geom);
  vol{2} = rmfield(vol{2},'unit');

  % new way - high level:
  tmpcfg = ft_keyval2cfg(arg{:});
  tmpcfg.method = 'bem_dipoli';
  vol{3} = ft_prepare_headmodel(tmpcfg,geom.bnd);
  vol{3} = rmfield(vol{3},'unit');

  % compute the leadfields for a comparison
  [vol{1}, elec] = ft_prepare_vol_sens(vol{1}, elec);
  [vol{2}, elec] = ft_prepare_vol_sens(vol{2}, elec);
  [vol{3}, elec] = ft_prepare_vol_sens(vol{3}, elec);
  lf{1} = ft_compute_leadfield([0 10 60], elec, vol{1});
  lf{2} = ft_compute_leadfield([0 10 60], elec, vol{2});
  lf{3} = ft_compute_leadfield([0 10 60], elec, vol{3});

  % compare the leadfields in all possible combinations
  comb = nchoosek(1:numel(vol),2);
  for j=1:size(comb,1)
    chk = comb(j,:);
    err = norm(lf{chk(1)} - lf{chk(2)}) / norm(lf{chk(1)});
    if err>0.001
      error('combination %d %d not successful\n',chk(1),chk(2));
    end
  end

end

author: fieldtrip

repo: fieldtrip


text:
function [gx,gy]=gaussgradient(IM,sigma)
%GAUSSGRADIENT Gradient using first order derivative of Gaussian.
%  [gx,gy]=gaussgradient(IM,sigma) outputs the gradient image gx and gy of
%  image IM using a 2-D Gaussian kernel. Sigma is the standard deviation of
%  this kernel along both directions.
%
%  Contributed by Guanglei Xiong (xgl99@mails.tsinghua.edu.cn)
%  at Tsinghua University, Beijing, China.

%determine the appropriate size of kernel. The smaller epsilon, the larger
%size.
epsilon=1e-2;
halfsize=ceil(sigma*sqrt(-2*log(sqrt(2*pi)*sigma*epsilon)));
size=2*halfsize+1;
%generate a 2-D Gaussian kernel along x direction
for i=1:size
    for j=1:size
        u=[i-halfsize-1 j-halfsize-1];
        hx(i,j)=gauss(u(1),sigma)*dgauss(u(2),sigma);
    end
end
hx=hx/sqrt(sum(sum(abs(hx).*abs(hx))));
%generate a 2-D Gaussian kernel along y direction
hy=hx';
%2-D filtering
gx=imfilter(IM,hx,'replicate','conv');
gy=imfilter(IM,hy,'replicate','conv');

function y = gauss(x,sigma)
%Gaussian
y = exp(-x^2/(2*sigma^2)) / (sigma*sqrt(2*pi));

function y = dgauss(x,sigma)
%first order derivative of Gaussian
y = -x * gauss(x,sigma) / sigma^2;


text:
function signal_mat = focus(kgrid, input_signal, source_mask, focus_position, sound_speed)
%FOCUS    Create input signal based on source mask and focus position.
%
% DESCRIPTION:
%       focus takes a single input signal and a source mask and creates an
%       input signal matrix (with one input signal for each source point).
%       The appropriate time delays required to focus the signals at a
%       given position in Cartesian space are automatically added based on
%       the user inputs for focus_position and sound_speed.
%
% USAGE:
%       input_signal_mat = focus(kgrid, input_signal, source_mask, focus_position, sound_speed)
%
% INPUTS:
%       kgrid           - k-Wave grid structure returned by makeGrid
%       input_signal    - single time series input
%       source_mask     - matrix specifying the positions of the time
%                         varying source distribution (i.e., source.p_mask
%                         or source.u_mask)
%       focus_position  - position of the focus in Cartesian coordinates
%       sound_speed     - scalar sound speed
%
% OUTPUTS:
%       input_signal_mat - matrix of time series following the source
%                         points using MATLAB's column-wise linear index
%                         ordering
%
% ABOUT:
%       author          - Bradley Treeby
%       date            - 21st February 2012
%       last update     - 25th July 2013
%
% 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/>.

% check that kgrid.t_array is defined
if strcmp(kgrid.t_array, 'auto')
    error('kgrid.t_array must be defined');
end

% check that the input sound speed is scalar
if numel(sound_speed) ~= 1
    error('sound_speed input must be scalar');
end

% calculate the distance from every point in the source mask to the focus
% position
switch kgrid.dim
    case 1
        dist = abs(kgrid.x(source_mask == 1) - focus_position(1));
    case 2
        dist = sqrt( (kgrid.x(source_mask == 1) - focus_position(1)).^2 + ...
            (kgrid.y(source_mask == 1) - focus_position(2)).^2 );
    case 3
        dist = sqrt( (kgrid.x(source_mask == 1) - focus_position(1)).^2 + ...
            (kgrid.y(source_mask == 1) - focus_position(2)).^2 + ...
            (kgrid.z(source_mask == 1) - focus_position(3)).^2 );
end

% convert the distance to units of time points
dist = round(dist./(kgrid.dt*sound_speed));

% convert time points to relative delays
dist = -(dist - max(dist(:)));
max_delay = max(dist(:));

% create an input matrix
signal_mat = zeros(length(dist), length(input_signal) + max_delay);

% assign the input signal
for source_index = 1:length(dist)
    delay = dist(source_index);
    signal_mat(source_index, :) = [zeros(1, delay), input_signal, zeros(1, max_delay - delay)];
end

repo: Aero-matlab


text:
function [muX, varX, model]  = vargplvmDynamicsUpdateModelTestVar(model, x, varx, y)
% VARGPLVMDYNAMICSUPDATEMODELTESTVAR return the original variational means and
% variances for the test data and update the model to be prepared for
% prediction
% DESC in the dynamical vargplvm the original variational parameters are
% reparametrized so that the optimiser optimises the new free ones and then
% these are mapped back to the original ones. This function takes the new
% free parameters which are optimised and returns the original ones. Also,
% because there is coupling in the variational parameters (because the
% original ones are obtained via a nonlinear operation on the free ones
% also involving Kt) after introducing test points (which have their own
% latent positions) the whole variational distribution is changed.
% Therefore, the updated model is also returned.
% ARG model : The vargplvm model which contains dynamics.
% ARG x : the optimised variational means (the free parameters, not the ones) of the test points
% ARG varx : the optimised variational variances (the free parameters, not the
% ones) of the test points
% ARG y: the test points
% RETURN muX : the original,true variational means corresponding to x
% RETURN varX : the original,true variational variances corresponding to varx
% RETURN model : the updated model, prepared for prediction

%
% COPYRIGHT : Michalis K. Titsias, 2009-2011
% COPYRIGHT : Neil D. Lawrence, 2009-2011
%
% VARGPLVM

Nstar = size(x,1);
N = model.N;

if isfield(model.dynamics, 'seq') && ~isempty(model.dynamics.seq)
   mask = sum(isnan(y),2);
   indexObservedData = find(mask==0)';
   yOb = y(indexObservedData, :);

   modelTest = model;
   modelTest.dynamics.t = model.dynamics.t_star;
   modelTest.dynamics.vardist.means = x;
   modelTest.dynamics.vardist.covars = varx;
   modelTest.dynamics.N = size(modelTest.dynamics.vardist.means, 1);
   modelTest.vardist.numData = modelTest.dynamics.N;
   modelTest.vardist.nParams = 2*prod(size(modelTest.dynamics.vardist.means));
   Kt = kernCompute(model.dynamics.kern, model.dynamics.t_star);
   modelTest.dynamics.Kt = Kt;
   modelTest.dynamics.fixedKt = 1;
   % This enables the computation of the variational means and covars
   modelTest.dynamics.seq = [];
   modelTest = vargpTimeDynamicsUpdateStats(modelTest);

   % Variational means and covars for the test block data
   % latent variables
   muX = modelTest.vardist.means;
   varX = modelTest.vardist.covars;

   if ~isempty(yOb)
      % Since the test block also might have a fully observed subblock, you
      % need to update the model structure (this for the subseqeucent prediction
      % using the PostMeanVar function)
      vardistOb = modelTest.vardist;
      vardistOb.means = modelTest.vardist.means(indexObservedData, :);
      vardistOb.covars = modelTest.vardist.covars(indexObservedData, :);
      vardistOb.nParams = 2*prod(size(vardistOb.means));
      vardistOb.numData = size(vardistOb.means,1);
      % Psi statistics for the data of the new block/sequence which have a fully
      % observed features/dimensions
      obsPsi0 = kernVardistPsi0Compute(model.kern, vardistOb);
      obsPsi1 = kernVardistPsi1Compute(model.kern, vardistOb, model.X_u);
      obsPsi2 = kernVardistPsi2Compute(model.kern, vardistOb, model.X_u);
      model.N = model.N + size(yOb,1);
      myOb = yOb - repmat(model.bias,size(yOb,1),1);
      myOb = myOb./repmat(model.scale,size(yOb,1),1);
      model.m = [model.m; myOb];
      model.dynamics.t = [model.dynamics.t; model.dynamics.t_star(indexObservedData)];
      Kt = zeros(N+size(indexObservedData,2),N+size(indexObservedData,2));
      Kt(1:N,1:N) = model.dynamics.Kt;
      tmpKt = modelTest.dynamics.Kt(indexObservedData,:);
      tmpKt = tmpKt(:,indexObservedData);
      Kt(N+1:end, N+1:end) = tmpKt;
      model.dynamics.Kt = Kt;
      model.dynamics.seq = [model.dynamics.seq, (N+size(indexObservedData,2))];

      model.Psi1 = [model.Psi1; obsPsi1];
      model.Psi2 = model.Psi2 + obsPsi2;
      model.Psi0 = model.Psi0 + obsPsi0;

      % dynamics (reparametrized) variational parameters (barmu and lambdas)
      model.dynamics.vardist.means = [model.dynamics.vardist.means; x(indexObservedData,:)];
      model.dynamics.vardist.covars = [model.dynamics.vardist.covars; varx(indexObservedData,:)];
      model.dynamics.N = size(model.dynamics.vardist.means, 1);
      model.dynamics.vardist.transforms.index = (size(model.dynamics.vardist.means,2)*model.N+1):size(model.dynamics.vardist.means,2)*model.N*2;
      model.dynamics.vardist.numData =  size(model.dynamics.vardist.means, 1);
      model.dynamics.vardist.nParams = 2*prod(size(model.dynamics.vardist.means));

      % the actual variational distrbution
      model.vardist.means = [model.vardist.means; vardistOb.means];
      model.vardist.covars = [model.vardist.covars; vardistOb.covars];
      model.vardist.transforms.type = model.vardist.transforms.type;
      model.vardist.transforms.index = (size(model.vardist.means,2)*model.N+1):size(model.vardist.means,2)*model.N*2;
      model.vardist.numData =  size(model.vardist.means, 1);
      model.vardist.nParams = 2*prod(size(model.vardist.means));

      % update in the model structure useful for prediction
      model.C = model.invLm * model.Psi2 * model.invLmT;
      model.TrC = sum(diag(model.C)); % Tr(C)
      model.At = (1/model.beta) * eye(size(model.C,1)) + model.C; % At = beta^{-1} I + C
      model.Lat = jitChol(model.At)';
      model.invLat = model.Lat\eye(size(model.Lat,1));
      model.invLatT = model.invLat';
      model.logDetAt = 2*(sum(log(diag(model.Lat)))); % log |At|
      model.P1 = model.invLat * model.invLm; % M x M
      model.P = model.P1 * (model.Psi1' * model.m);
      model.TrPP = sum(sum(model.P .* model.P));
   end
else
   % Augment the time vector to include the timestamp of the new point
   model.dynamics.t = [model.dynamics.t; model.dynamics.t_star];
   % Augment the reparametrized variational parameters mubar and lambda
   model.dynamics.vardist.means = [model.dynamics.vardist.means; x];
   model.dynamics.vardist.covars = [model.dynamics.vardist.covars; varx];
   model.dynamics.N =  model.dynamics.N + Nstar;
   model.vardist.numData = model.dynamics.N;
   model.dynamics.vardist.numData = model.dynamics.N;
   model.vardist.nParams = 2*prod(size(model.dynamics.vardist.means));
   model.dynamics.vardist.nParams = 2*prod(size(model.dynamics.vardist.means));
   model.dynamics.seq = model.dynamics.N; %%%%%% Chec
   model = vargpTimeDynamicsUpdateStats(model);
   vardist2 = model.vardist;
   vardist2.means = model.vardist.means(1:end-Nstar,:);
   vardist2.covars = model.vardist.covars(1:end-Nstar,:);
   vardist2.nParams = 2*prod(size(vardist2.means));
   vardist2.numData = size(vardist2.means,1);
   model.Psi0 = kernVardistPsi0Compute(model.kern, vardist2);
   model.Psi1 = kernVardistPsi1Compute(model.kern, vardist2, model.X_u);
   model.Psi2 = kernVardistPsi2Compute(model.kern, vardist2, model.X_u);
   model.C = model.invLm * model.Psi2 * model.invLmT;
   model.TrC = sum(diag(model.C)); % Tr(C)
   model.At = (1/model.beta) * eye(size(model.C,1)) + model.C; % At = beta^{-1} I + C
   model.Lat = jitChol(model.At)';
   model.invLat = model.Lat\eye(size(model.Lat,1));
   model.invLatT = model.invLat';
   model.logDetAt = 2*(sum(log(diag(model.Lat)))); % log |At|
   model.P1 = model.invLat * model.invLm; % M x M
   model.P = model.P1 * (model.Psi1' * model.m);
   model.TrPP = sum(sum(model.P .* model.P));

   % Variational means and covars for the test block data
   % latent variables
   muX = model.vardist.means(N+1:end,:);
   varX = model.vardist.covars(N+1:end,:);

end

model.X = model.vardist.means;

author: SheffieldML

repo: vargplvm