text:
function [ o, x, w ] = en_r2_11_1 ( n, option )

%*****************************************************************************80
%
%% EN_R2_11_1 implements the Stroud rule 11.1 for region EN_R2.
%
%  Discussion:
%
%    The rule has order
%
%      O = ( 4 * N^5 - 20 * N^4 + 140 * N^3 - 130 * N^2 + 96 * N + 15 ) / 15.
%
%    The rule has precision P = 11.
%
%    EN_R2 is the entire N-dimensional space with weight function
%
%      w(x) = exp ( - x1^2 - x2^2 ... - xn^2 )
%
%    There are two versions of each rule, chosen by setting the
%    OPTION variable to 1 or 2.
%
%    The rule as tabulated by Stenger is available for N = 2 through 20.
%    This function accepts N = 3 through 5.
%
%     N    O
%    __  ___
%     3  151
%     4  417
%     5  983
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971,
%    ISBN: 0130438936,
%    LC: QA311.S85.
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%    3 <= N <= 5.
%
%    Input, integer OPTION, chooses rule option 1 or 2.
%
%    Output, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  if ( n < 3 | 5 < n )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'EN_R2_11_1 - Fatal error!\n' );
    fprintf ( 1, '  3 <= N <= 5 required.\n' );
    error ( 'EN_R2_11_1 - Fatal error!' )
  end

  if ( nargin < 2 )
    option = 1;
  end

  if ( option < 1 | 2 < option )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'EN_R2_11_1 - Fatal error!\n' );
    fprintf ( 1, '  1 <= OPTION <= 2 required.\n' );
    error ( 'EN_R2_11_1 - Fatal error!' )
  end

  o = ( 4 * n^5 - 20 * n^4 + 140 * n^3 - 130 * n^2 + 96 * n + 15 ) / 15;
  volume = sqrt ( pi^n );

  if ( n == 3 & option == 1 )
    u =     0.235060497367449E+01;
    v =     0.436077411927617E+00;
    w2 =    0.133584907401370E+01;
    b0 =  - 0.881591029957858E+01;
    b1 =  - 0.751996143360650E-01;
    b2 =    0.621743189471515E+01;
    b3 =    0.241426451456494E+00;
    b4 =  - 0.120709739276065E-02;
    b5 =  - 0.427751221210138E+01;
    b6 =    0.550169924840163E-01;
    b7 =    0.237084999634707E-01;
    b8 =  - 0.169791992887741E-02;
    b9 =  - 0.252266276123350E-04;
    b10 =   0.326777873717691E+01;
    b11 =   0.968469949206802E-02;
    b12 =   0.789754514877422E-03;
    b13 =   0.000000000000000E+00;
    b14 =   0.000000000000000E+00;
    b15 =   0.000000000000000E+00;
  elseif ( n == 3 & option == 2 )
    u =     0.235060497367449E+01;
    v =     0.133584907401370E+01;
    w2 =    0.436077411927617E+00;
    b0 =  - 0.141214037032900E+02;
    b1 =  - 0.803730274707282E-01;
    b2 =    0.235546545595906E+00;
    b3 =    0.888123191556611E+01;
    b4 =    0.142467131155533E-03;
    b5 =    0.582993124006494E-01;
    b6 =  - 0.561099173155661E+01;
    b7 =  - 0.204028691521686E-02;
    b8 =    0.252880089932256E-01;
    b9 =  - 0.814378678627283E-04;
    b10 =   0.804353953375146E-02;
    b11 =   0.393451849690453E+01;
    b12 =   0.171183493169724E-03;
    b13 =   0.000000000000000E+00;
    b14 =   0.000000000000000E+00;
    b15 =   0.000000000000000E+00;
  elseif ( n == 4 & option == 1 )
    u =     0.235060497367449E+01;
    v =     0.436077411927617E+00;
    w2 =    0.133584907401370E+01;
    b0 =    0.241502736147339E+03;
    b1 =  - 0.196095938531478E+00;
    b2 =  - 0.128675737999280E+03;
    b3 =    0.307568784278696E+00;
    b4 =  - 0.480908422319460E-02;
    b5 =    0.698087019367085E+02;
    b6 =    0.631837143743771E-01;
    b7 =    0.392226151971179E-01;
    b8 =  - 0.300948471646799E-02;
    b9 =  - 0.650235306755170E-04;
    b10 = - 0.386951974646715E+02;
    b11 =   0.171656829095787E-01;
    b12 =   0.139980343116450E-02;
    b13 =   0.101552487093372E-04;
    b14 =   0.222435922356439E+02;
    b15 =   0.000000000000000E+00;
  elseif ( n == 4 & option == 2 )
    u =     0.235060497367449E+01;
    v =     0.133584907401370E+01;
    w2 =    0.436077411927617E+00;
    b0 =  - 0.151944464736584E+03;
    b1 =  - 0.223498438689039E+00;
    b2 =    0.243574919068010E+00;
    b3 =    0.634373877008693E+02;
    b4 =  - 0.782065187814018E-04;
    b5 =    0.911833754536616E-01;
    b6 =  - 0.238927288245914E+02;
    b7 =  - 0.422314408318853E-02;
    b8 =    0.448218289217760E-01;
    b9 =  - 0.138053374667391E-03;
    b10 =   0.607473265800655E-02;
    b11 =   0.697375246129742E+01;
    b12 =   0.303414841680135E-03;
    b13 = - 0.314574391771792E-05;
    b14 =   0.409103498175100E-02;
    b15 =   0.000000000000000E+00;
  elseif ( n == 5 & option == 1 )
    u =     0.235060497367449E+01;
    v =     0.436077411927617E+00;
    w2 =    0.133584907401370E+01;
    b0 =    0.255885269311763E+04;
    b1 =  - 0.439598677491526E+00;
    b2 =  - 0.106541406144610E+04;
    b3 =    0.453540909054264E+00;
    b4 =  - 0.132100905623778E-01;
    b5 =    0.418606568954203E+03;
    b6 =    0.511394563043680E-01;
    b7 =    0.645581013845604E-01;
    b8 =  - 0.533417277494500E-02;
    b9 =  - 0.137981626254496E-03;
    b10 = - 0.147436933189884E+03;
    b11 =   0.304253807765057E-01;
    b12 =   0.248108698207828E-02;
    b13 =   0.113652094546015E-04;
    b14 =   0.394257407160391E+02;
    b15 =   0.331725011358320E-05;
  elseif ( n == 5 & option == 2 )
    u =     0.235060497367449E+01;
    v =     0.133584907401370E+01;
    w2 =    0.436077411927617E+00;
    b0 =  - 0.761305347548192E+03;
    b1 =  - 0.536360805019297E+00;
    b2 =    0.110669832078736E+00;
    b3 =    0.246421088923968E+03;
    b4 =  - 0.773649327968607E-03;
    b5 =    0.169088641205970E+00;
    b6 =  - 0.670700680243651E+02;
    b7 =  - 0.856090560229205E-02;
    b8 =    0.794446232770302E-01;
    b9 =  - 0.220272863263544E-03;
    b10 = - 0.373515812228225E-02;
    b11 =   0.123606544052884E+02;
    b12 =   0.537788804557843E-03;
    b13 = - 0.122101861480881E-04;
    b14 =   0.725117070759373E-02;
    b15 =   0.331725011358320E-05;
  end

  x = zeros(n,o);
  w = zeros(o,1);

  k = 0;
%
%  1 point.
%
  k = k + 1;
% x(1:n,k) = 0.0;
  w(k) = b0;
%
%  2 * N points.
%
  for i = 1 : n
    k = k + 1;
    x(i,k) = - u;
    w(k) = b1;
    k = k + 1;
    x(i,k) = + u;
    w(k) = b1;
  end
%
%  2 * N points.
%
  for i = 1 : n
    k = k + 1;
    x(i,k) = - v;
    w(k) = b2;
    k = k + 1;
    x(i,k) = + v;
    w(k) = b2;
  end
%
%  2 * N points.
%
  for i = 1 : n
    k = k + 1;
    x(i,k) = - w2;
    w(k) = b3;
    k = k + 1;
    x(i,k) = + w2;
    w(k) = b3;
  end
%
%  4 * ( N * ( N - 1 ) / 2 ) points.
%
  for i = 1 : n - 1
    for j = i + 1 : n
      k = k + 1;
      x(i,k) = - u;
      x(j,k) = - u;
      w(k) = b4;
      k = k + 1;
      x(i,k) = - u;
      x(j,k) = + u;
      w(k) = b4;
      k = k + 1;
      x(i,k) = + u;
      x(j,k) = - u;
      w(k) = b4;
      k = k + 1;
      x(i,k) = + u;
      x(j,k) = + u;
      w(k) = b4;
    end
  end
%
%  4 * ( N * ( N - 1 ) / 2 ) points.
%
  for i = 1 : n - 1
    for j = i + 1 : n
      k = k + 1;
      x(i,k) = - v;
      x(j,k) = - v;
      w(k) = b5;
      k = k + 1;
      x(i,k) = - v;
      x(j,k) = + v;
      w(k) = b5;
      k = k + 1;
      x(i,k) = + v;
      x(j,k) = - v;
      w(k) = b5;
      k = k + 1;
      x(i,k) = + v;
      x(j,k) = + v;
      w(k) = b5;
    end
  end
%
%  4 * ( N * ( N - 1 ) / 2 ) points.
%
  for i = 1 : n - 1
    for j = i + 1 : n
      k = k + 1;
      x(i,k) = - w2;
      x(j,k) = - w2;
      w(k) = b6;
      k = k + 1;
      x(i,k) = - w2;
      x(j,k) = + w2;
      w(k) = b6;
      k = k + 1;
      x(i,k) = + w2;
      x(j,k) = - w2;
      w(k) = b6;
      k = k + 1;
      x(i,k) = + w2;
      x(j,k) = + w2;
      w(k) = b6;
    end
  end
%
%  4 * ( N * ( N - 1 ) ) points.
%
  for i = 1 : n - 1
    for j = i + 1 : n
      k = k + 1;
      x(i,k) = - u;
      x(j,k) = - v;
      w(k) = b7;
      k = k + 1;
      x(i,k) = - u;
      x(j,k) = + v;
      w(k) = b7;
      k = k + 1;
      x(i,k) = + u;
      x(j,k) = - v;
      w(k) = b7;
      k = k + 1;
      x(i,k) = + u;
      x(j,k) = + v;
      w(k) = b7;
      k = k + 1;
      x(i,k) = - v;
      x(j,k) = - u;
      w(k) = b7;
      k = k + 1;
      x(i,k) = - v;
      x(j,k) = + u;
      w(k) = b7;
      k = k + 1;
      x(i,k) = + v;
      x(j,k) = - u;
      w(k) = b7;
      k = k + 1;
      x(i,k) = + v;
      x(j,k) = + u;
      w(k) = b7;
    end
  end
%
%  4 * ( N * ( N - 1 ) ) points.
%
  for i = 1 : n - 1
    for j = i + 1 : n
      k = k + 1;
      x(i,k) = - u;
      x(j,k) = - w2;
      w(k) = b8;
      k = k + 1;
      x(i,k) = - u;
      x(j,k) = + w2;
      w(k) = b8;
      k = k + 1;
      x(i,k) = + u;
      x(j,k) = - w2;
      w(k) = b8;
      k = k + 1;
      x(i,k) = + u;
      x(j,k) = + w2;
      w(k) = b8;
      k = k + 1;
      x(i,k) = - w2;
      x(j,k) = - u;
      w(k) = b8;
      k = k + 1;
      x(i,k) = - w2;
      x(j,k) = + u;
      w(k) = b8;
      k = k + 1;
      x(i,k) = + w2;
      x(j,k) = - u;
      w(k) = b8;
      k = k + 1;
      x(i,k) = + w2;
      x(j,k) = + u;
      w(k) = b8;
    end
  end
%
%  8 * ( N * ( N - 1 ) * ( N - 2 ) / 6 ) points.
%
  for i = 1 : n - 2
    for j = i + 1 : n - 1
      for l = j + 1 : n
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = - u;
        x(l,k) = - u;
        w(k) = b9;
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = - u;
        x(l,k) = + u;
        w(k) = b9;
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = + u;
        x(l,k) = - u;
        w(k) = b9;
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = + u;
        x(l,k) = + u;
        w(k) = b9;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = - u;
        x(l,k) = - u;
        w(k) = b9;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = - u;
        x(l,k) = + u;
        w(k) = b9;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = + u;
        x(l,k) = - u;
        w(k) = b9;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = + u;
        x(l,k) = + u;
        w(k) = b9;
      end
    end
  end
%
%  8 * ( N * ( N - 1 ) * ( N - 2 ) / 6 ) points.
%
  for i = 1 : n - 2
    for j = i + 1 : n - 1
      for l = j + 1 : n
        k = k + 1;
        x(i,k) = - v;
        x(j,k) = - v;
        x(l,k) = - v;
        w(k) = b10;
        k = k + 1;
        x(i,k) = - v;
        x(j,k) = - v;
        x(l,k) = + v;
        w(k) = b10;
        k = k + 1;
        x(i,k) = - v;
        x(j,k) = + v;
        x(l,k) = - v;
        w(k) = b10;
        k = k + 1;
        x(i,k) = - v;
        x(j,k) = + v;
        x(l,k) = + v;
        w(k) = b10;
        k = k + 1;
        x(i,k) = + v;
        x(j,k) = - v;
        x(l,k) = - v;
        w(k) = b10;
        k = k + 1;
        x(i,k) = + v;
        x(j,k) = - v;
        x(l,k) = + v;
        w(k) = b10;
        k = k + 1;
        x(i,k) = + v;
        x(j,k) = + v;
        x(l,k) = - v;
        w(k) = b10;
        k = k + 1;
        x(i,k) = + v;
        x(j,k) = + v;
        x(l,k) = + v;
        w(k) = b10;
      end
    end
  end
%
%  8 * ( N * ( N - 1 ) * ( N - 2 ) / 6 ) points.
%
  for i = 1 : n - 2
    for j = i + 1 : n - 1
      for l = j + 1 : n
        k = k + 1;
        x(i,k) = - w2;
        x(j,k) = - w2;
        x(l,k) = - w2;
        w(k) = b11;
        k = k + 1;
        x(i,k) = - w2;
        x(j,k) = - w2;
        x(l,k) = + w2;
        w(k) = b11;
        k = k + 1;
        x(i,k) = - w2;
        x(j,k) = + w2;
        x(l,k) = - w2;
        w(k) = b11;
        k = k + 1;
        x(i,k) = - w2;
        x(j,k) = + w2;
        x(l,k) = + w2;
        w(k) = b11;
        k = k + 1;
        x(i,k) = + w2;
        x(j,k) = - w2;
        x(l,k) = - w2;
        w(k) = b11;
        k = k + 1;
        x(i,k) = + w2;
        x(j,k) = - w2;
        x(l,k) = + w2;
        w(k) = b11;
        k = k + 1;
        x(i,k) = + w2;
        x(j,k) = + w2;
        x(l,k) = - w2;
        w(k) = b11;
        k = k + 1;
        x(i,k) = + w2;
        x(j,k) = + w2;
        x(l,k) = + w2;
        w(k) = b11;
      end
    end
  end
%
%  8 * ( N * ( N - 1 ) * ( N - 2 ) / 2 ) points.
%
  for i = 1 : n - 2
    for j = i + 1 : n - 1
      for l = j + 1 : n
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = - u;
        x(l,k) = - v;
        w(k) = b12;
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = - u;
        x(l,k) = + v;
        w(k) = b12;
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = + u;
        x(l,k) = - v;
        w(k) = b12;
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = + u;
        x(l,k) = + v;
        w(k) = b12;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = - u;
        x(l,k) = - v;
        w(k) = b12;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = - u;
        x(l,k) = + v;
        w(k) = b12;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = + u;
        x(l,k) = - v;
        w(k) = b12;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = + u;
        x(l,k) = + v;
        w(k) = b12;
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = - v;
        x(l,k) = - u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = - v;
        x(l,k) = + u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = + v;
        x(l,k) = - u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = + v;
        x(l,k) = + u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = - v;
        x(l,k) = - u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = - v;
        x(l,k) = + u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = + v;
        x(l,k) = - u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = + v;
        x(l,k) = + u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = - v;
        x(j,k) = - u;
        x(l,k) = - u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = - v;
        x(j,k) = - u;
        x(l,k) = + u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = - v;
        x(j,k) = + u;
        x(l,k) = - u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = - v;
        x(j,k) = + u;
        x(l,k) = + u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = + v;
        x(j,k) = - u;
        x(l,k) = - u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = + v;
        x(j,k) = - u;
        x(l,k) = + u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = + v;
        x(j,k) = + u;
        x(l,k) = - u;
        w(k) = b12;
        k = k + 1;
        x(i,k) = + v;
        x(j,k) = + u;
        x(l,k) = + u;
        w(k) = b12;
      end
    end
  end
%
%  16 * ( N * ( N - 1 ) * ( N - 2 ) * ( N - 3 ) / 24 ) points.
%
  for i = 1 : n - 3
    for j = i + 1 : n - 2
      for l = j + 1 : n - 1
        for m = l + 1 : n
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = - u;
          x(l,k) = - u;
          x(m,k) = - u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = - u;
          x(l,k) = - u;
          x(m,k) = + u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = - u;
          x(l,k) = + u;
          x(m,k) = - u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = - u;
          x(l,k) = + u;
          x(m,k) = + u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = + u;
          x(l,k) = - u;
          x(m,k) = - u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = + u;
          x(l,k) = - u;
          x(m,k) = + u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = + u;
          x(l,k) = + u;
          x(m,k) = - u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = + u;
          x(l,k) = + u;
          x(m,k) = + u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = - u;
          x(l,k) = - u;
          x(m,k) = - u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = - u;
          x(l,k) = - u;
          x(m,k) = + u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = - u;
          x(l,k) = + u;
          x(m,k) = - u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = - u;
          x(l,k) = + u;
          x(m,k) = + u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = + u;
          x(l,k) = - u;
          x(m,k) = - u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = + u;
          x(l,k) = - u;
          x(m,k) = + u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = + u;
          x(l,k) = + u;
          x(m,k) = - u;
          w(k) = b13;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = + u;
          x(l,k) = + u;
          x(m,k) = + u;
          w(k) = b13;
        end
      end
    end
  end
%
%  16 * ( N * ( N - 1 ) * ( N - 2 ) * ( N - 3 ) / 24 ) points.
%
  for i = 1 : n - 3
    for j = i + 1 : n - 2
      for l = j + 1 : n - 1
        for m = l + 1 : n
          k = k + 1;
          x(i,k) = - v;
          x(j,k) = - v;
          x(l,k) = - v;
          x(m,k) = - v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = - v;
          x(j,k) = - v;
          x(l,k) = - v;
          x(m,k) = + v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = - v;
          x(j,k) = - v;
          x(l,k) = + v;
          x(m,k) = - v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = - v;
          x(j,k) = - v;
          x(l,k) = + v;
          x(m,k) = + v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = - v;
          x(j,k) = + v;
          x(l,k) = - v;
          x(m,k) = - v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = - v;
          x(j,k) = + v;
          x(l,k) = - v;
          x(m,k) = + v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = - v;
          x(j,k) = + v;
          x(l,k) = + v;
          x(m,k) = - v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = - v;
          x(j,k) = + v;
          x(l,k) = + v;
          x(m,k) = + v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = + v;
          x(j,k) = - v;
          x(l,k) = - v;
          x(m,k) = - v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = + v;
          x(j,k) = - v;
          x(l,k) = - v;
          x(m,k) = + v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = + v;
          x(j,k) = - v;
          x(l,k) = + v;
          x(m,k) = - v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = + v;
          x(j,k) = - v;
          x(l,k) = + v;
          x(m,k) = + v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = + v;
          x(j,k) = + v;
          x(l,k) = - v;
          x(m,k) = - v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = + v;
          x(j,k) = + v;
          x(l,k) = - v;
          x(m,k) = + v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = + v;
          x(j,k) = + v;
          x(l,k) = + v;
          x(m,k) = - v;
          w(k) = b14;
          k = k + 1;
          x(i,k) = + v;
          x(j,k) = + v;
          x(l,k) = + v;
          x(m,k) = + v;
          w(k) = b14;
        end
      end
    end
  end
%
%  All quintuples UUUUU with 32 sign combinations.
%
  for i1 = 1 : n - 4
    for i2 = i1 + 1 : n - 3
      for i3 = i2 + 1 : n - 2
        for i4 = i3 + 1 : n - 1
          for i5 = i4 + 1 : n
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = - u;
            x(i3,k) = - u;
            x(i4,k) = - u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = - u;
            x(i3,k) = - u;
            x(i4,k) = - u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = - u;
            x(i3,k) = - u;
            x(i4,k) = + u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = - u;
            x(i3,k) = - u;
            x(i4,k) = + u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = - u;
            x(i3,k) = + u;
            x(i4,k) = - u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = - u;
            x(i3,k) = + u;
            x(i4,k) = - u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = - u;
            x(i3,k) = + u;
            x(i4,k) = + u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = - u;
            x(i3,k) = + u;
            x(i4,k) = + u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = + u;
            x(i3,k) = - u;
            x(i4,k) = - u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = + u;
            x(i3,k) = - u;
            x(i4,k) = - u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = + u;
            x(i3,k) = - u;
            x(i4,k) = + u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = + u;
            x(i3,k) = - u;
            x(i4,k) = + u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = + u;
            x(i3,k) = + u;
            x(i4,k) = - u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = + u;
            x(i3,k) = + u;
            x(i4,k) = - u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = + u;
            x(i3,k) = + u;
            x(i4,k) = + u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = - u;
            x(i2,k) = + u;
            x(i3,k) = + u;
            x(i4,k) = + u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = - u;
            x(i3,k) = - u;
            x(i4,k) = - u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = - u;
            x(i3,k) = - u;
            x(i4,k) = - u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = - u;
            x(i3,k) = - u;
            x(i4,k) = + u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = - u;
            x(i3,k) = - u;
            x(i4,k) = + u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = - u;
            x(i3,k) = + u;
            x(i4,k) = - u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = - u;
            x(i3,k) = + u;
            x(i4,k) = - u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = - u;
            x(i3,k) = + u;
            x(i4,k) = + u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = - u;
            x(i3,k) = + u;
            x(i4,k) = + u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = + u;
            x(i3,k) = - u;
            x(i4,k) = - u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = + u;
            x(i3,k) = - u;
            x(i4,k) = - u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = + u;
            x(i3,k) = - u;
            x(i4,k) = + u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = + u;
            x(i3,k) = - u;
            x(i4,k) = + u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = + u;
            x(i3,k) = + u;
            x(i4,k) = - u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = + u;
            x(i3,k) = + u;
            x(i4,k) = - u;
            x(i5,k) = + u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = + u;
            x(i3,k) = + u;
            x(i4,k) = + u;
            x(i5,k) = - u;
            w(k) = b15;
            k = k + 1;
            x(i1,k) = + u;
            x(i2,k) = + u;
            x(i3,k) = + u;
            x(i4,k) = + u;
            x(i5,k) = + u;
            w(k) = b15;
          end
        end
      end
    end
  end

  return
end


text:
%[2010]-"Firefly algorithm,stochastic test functions and design
%optimization"

% (9/12/2020)

function FA = jFireflyAlgorithm(feat,label,opts)
% Parameters
lb    = 0;
ub    = 1;
thres = 0.5;
alpha = 1;      % constant
beta0 = 1;      % light amplitude
gamma = 1;      % absorbtion coefficient
theta = 0.97;   % control alpha

if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'beta0'), beta0 = opts.beta0; end
if isfield(opts,'gamma'), gamma = opts.gamma; end
if isfield(opts,'alpha'), alpha = opts.alpha; end
if isfield(opts,'theta'), theta = opts.theta; end
if isfield(opts,'thres'), thres = opts.thres; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial
X   = zeros(N,dim);
for i = 1:N
  for d = 1:dim
    X(i,d) = lb + (ub - lb) * rand();
  end
end
% Fitness
fit  = zeros(1,N);
fitG = inf;
for i = 1:N
  fit(i) = fun(feat,label,(X(i,:) > thres),opts);
  % Best solution
  if fit(i) < fitG
    fitG = fit(i);
    Xgb  = X(i,:);
  end
end
% Pre
curve = zeros(1,max_Iter);
curve(1) = fitG;
t = 2;
% Generation
while t <= max_Iter
  % Alpha update
  alpha      = alpha * theta;
  % Rank firefly based on their light intensity
  [fit, idx] = sort(fit,'ascend');
  X          = X(idx,:);
	for i = 1:N
    % The attractiveness parameter
    for j = 1:N
      % Update moves if firefly j brighter than firefly i
      if fit(i) > fit(j)
        % Compute Euclidean distance
        r    = sqrt(sum((X(i,:) - X(j,:)) .^ 2));
        % Beta (2)
        beta = beta0 * exp(-gamma * r ^ 2);
        for d = 1:dim
          % Update position (3)
          eps    = rand() - 0.5;
          X(i,d) = X(i,d) + beta * (X(j,d) - X(i,d)) + alpha * eps;
        end
        % Boundary
        XB = X(i,:); XB(XB > ub) = ub; XB(XB < lb) = lb;
        X(i,:) = XB;
        % Fitness
        fit(i) = fun(feat,label,(X(i,:) > thres),opts);
        % Update global best firefly
        if fit(i) < fitG
          fitG = fit(i);
          Xgb  = X(i,:);
        end
      end
    end
  end
  curve(t) = fitG;
  fprintf('\nGeneration %d Best (FA)= %f',t,curve(t))
  t = t + 1;
end
% Select features
Pos   = 1:dim;
Sf    = Pos((Xgb > thres) == 1);
sFeat = feat(:,Sf);
% Store results
FA.sf = Sf;
FA.ff = sFeat;
FA.nf = length(Sf);
FA.c  = curve;
FA.f  = feat;
FA.l  = label;
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
function [lstmStates, attnData, attnInfos] = rnnLayerForward(W_rnn, W_emb, prevState, input, masks, params, rnnFlags, attnData, model)
% Running Multi-layer RNN for one time step.
% Input:
%   W_rnn: recurrent connections of multiple layers, e.g., W_rnn{ll}.
%   prevState: previous hidden state, e.g., for LSTM, prevState.c{ll}, prevState.h{ll}.
%   input: indices for the current batch
%   isTest: 1 -- don't store intermediate results in each state
%   isAttn: for attention, require attnData to be non-empty, has
%     attnData.srcHidVecsOrig and attnData.srcLens.
%   isDecoder: 0 -- encoder, 1 -- decoder
% Output:
%   nextState
%
% Thang Luong @ 2015, <lmthang@stanford.edu>

T = size(input, 2);
lstmStates = cell(T, 1);

% attention
attnInfos = cell(T, 1);
if rnnFlags.attn && rnnFlags.decode == 0 % encoder
  assert(T <= params.numSrcHidVecs);
  attnData.srcHidVecsOrig = zeroMatrix([params.lstmSize, params.curBatchSize, params.numSrcHidVecs], params.isGPU, params.dataType);
end

for tt=1:T % time
  % local monotonic alignment
  if params.attnLocalMono
    attnData.tgtPos = tt;
  end

  % multi-layer RNN
  [prevState, attnInfos{tt}] = rnnStepLayerForward(W_rnn, W_emb(:, input(:, tt)), prevState, masks(:, tt), params, rnnFlags, attnData, model);

  % encoder, attention
  if rnnFlags.attn && rnnFlags.decode == 0
    attnData.srcHidVecsOrig(:, :, tt) = prevState{end}.h_t;
  end

  % store all states
  lstmStates{tt} = prevState;
end

author: lmthang


text:
classdef PropertyChangedData < event.EventData & dynamicprops
    % Event data for widget property value changes
    %
    % Syntax:
    %           obj = wt.eventdata.PropertyChangedData(propName,newValue,oldValue)
    %           obj = wt.eventdata.PropertyChangedData(propName,valueChangedEvent)
    %           obj = wt.eventdata.PropertyChangedData(...,'p1',v1,...)
    %

    % Copyright 2020-2021 The MathWorks, Inc.

    %% Properties
    properties (SetAccess = protected)
        Value
        PreviousValue
        Property char
    end %properties

    %% Constructor / destructor
    methods
        function obj = PropertyChangedData(propName,newValue,varargin)

            % Set the changed property name
            obj.Property = propName;

            % Is input a MATLAB eventdata?
            if isa(newValue,'matlab.ui.eventdata.ValueChangedData')
                obj.Value = newValue.Value;
                obj.PreviousValue = newValue.PreviousValue;
            elseif isa(newValue,'matlab.ui.eventdata.ValueChangingData')
                obj.Value = newValue.Value;
            else
                % No - use the value directly
                obj.Value = newValue;
            end

            % Was a previous value provided?
            if mod(numel(varargin),2)
                obj.PreviousValue = varargin{1};
                varargin(1) = [];
            end

            % Any remaining varargin are dynamic property-value pairs
            for idx=1:numel(varargin)
                thisProp = varargin{idx};
                thisValue = remArgs.(thisProp);
                obj.addprop(thisProp);
                obj.(thisProp) = thisValue;
            end

        end %constructor
    end %methods

end % classdef

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function varargout = process_import_channel( varargin )
% PROCESS_IMPORT_CHANNEL: Import a raw file in the database.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2012-2017

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % ===== EEG DEFAULTS =====
    strList = {''};
    % Get registered Brainstorm EEG defaults
    bstDefaults = bst_get('EegDefaults');
    % Build a list of strings representing all the defaults
    for iGroup = 1:length(bstDefaults)
        for iDef = 1:length(bstDefaults(iGroup).contents)
            strList{end+1} = [bstDefaults(iGroup).name ': ' bstDefaults(iGroup).contents(iDef).name];
        end
    end

    % ===== PROCESS =====
    % Description the process
    sProcess.Comment     = 'Set channel file';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Import', 'Channel file'};
    sProcess.Index       = 30;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Epilepsy#Prepare_the_channel_file';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'raw', 'matrix'};
    sProcess.OutputTypes = {'data', 'raw', 'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    sProcess.isSeparator = 1;
    % Option: File to import
    sProcess.options.channelfile.Comment = 'File to import:';
    sProcess.options.channelfile.Type    = 'filename';
    sProcess.options.channelfile.Value   = {...
        '', ...                                % Filename
        '', ...                                % FileFormat
        'open', ...                            % Dialog type: {open,save}
        'Import channel file', ...             % Window title
        'ImportChannel', ...                   % LastUsedDir: {ImportData,ImportChannel,ImportAnat,ExportChannel,ExportData,ExportAnat,ExportProtocol,ExportImage,ExportScript}
        'single', ...                          % Selection mode: {single,multiple}
        'files_and_dirs', ...                  % Selection mode: {files,dirs,files_and_dirs}
        bst_get('FileFilters', 'channel'), ... % Get all the available file formats
        'ChannelIn'};                          % DefaultFormats
    % Option: Default channel files
    sProcess.options.usedefault.Comment = 'Or use default:';
    sProcess.options.usedefault.Type    = 'combobox';
    sProcess.options.usedefault.Value   = {1, strList};
    % Separator
    sProcess.options.separator.Type = 'separator';
    sProcess.options.separator.Comment = ' ';
    % Align sensors
    sProcess.options.channelalign.Comment = 'Align sensors using headpoints';
    sProcess.options.channelalign.Type    = 'checkbox';
    sProcess.options.channelalign.Value   = 1;
    % Fix units
    sProcess.options.fixunits.Comment = 'Fix distance units automatically';
    sProcess.options.fixunits.Type    = 'checkbox';
    sProcess.options.fixunits.Value   = 1;
    % Fix units
    sProcess.options.vox2ras.Comment = 'Apply voxel=>subject transformation from the MRI';
    sProcess.options.vox2ras.Type    = 'checkbox';
    sProcess.options.vox2ras.Value   = 1;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    OutputFiles = {};

    % ===== GET FILE =====
    % Get filename to import
    ChannelFile = sProcess.options.channelfile.Value{1};
    FileFormat  = sProcess.options.channelfile.Value{2};

    % ===== GET DEFAULT =====
    if isempty(ChannelFile)
        % Get registered Brainstorm EEG defaults
        bstDefaults = bst_get('EegDefaults');
        % Get default channel file
        iSel   = sProcess.options.usedefault.Value{1};
        strDef = sProcess.options.usedefault.Value{2}{iSel};
        % If there is something selected
        if ~isempty(strDef)
            % Format: "group: name"
            cDef   = strtrim(str_split(strDef, ':'));
            % Find the selected group in the list
            iGroup = find(strcmpi(cDef{1}, {bstDefaults.name}));
            % If group is found
            if ~isempty(iGroup)
                % Find the selected default in the list
                iDef = find(strcmpi(cDef{2}, {bstDefaults(iGroup).contents.name}));
                % If default was found
                if ~isempty(iDef)
                    ChannelFile = bstDefaults(iGroup).contents(iDef).fullpath;
                end
            end
        end
        FileFormat = 'BST';
    end
    % If nothing was selected (file or default)
    if isempty(ChannelFile)
        bst_report('Error', sProcess, [], 'No channel file selected.');
        return
    end

    % ===== SET CHANNEL FILE =====
    % Get channel studies
    [tmp, iChanStudies] = bst_get('ChannelForStudy', [sInputs.iStudy]);
    iChanStudies = unique(iChanStudies);
    % Channel align
    ChannelAlign = 2 * double(sProcess.options.channelalign.Value);
    % Other options
    if isfield(sProcess.options, 'fixunits') && isfield(sProcess.options.fixunits, 'Value')
        isFixUnits = sProcess.options.fixunits.Value;
    else
        isFixUnits = 1;
    end
    if isfield(sProcess.options, 'vox2ras') && isfield(sProcess.options.vox2ras, 'Value')
        isApplyVox2ras = sProcess.options.vox2ras.Value;
    else
        isApplyVox2ras = 1;
    end
    % Import channel files
    ChannelReplace = 2;
    isSave = 1;
    import_channel(iChanStudies, ChannelFile, FileFormat, ChannelReplace, ChannelAlign, isSave, isFixUnits, isApplyVox2ras);
    % Return all the files in input
    OutputFiles = {sInputs.FileName};
end

author: brainstorm-tools


text:
function Show_EWT2D_Tensor(ewt2d)

%=========================================================================
%
% function Show_Tensor_EWT2D(ewt2d)
%
% Show each subband image in a subplot of a single figure.
% The images are renormalized with respect to the min and max of all
% images.
%
% Author: Jerome Gilles
% Institution: UCLA - Department of Mathematics
% Year: 2013
% Version: 1.0
% ========================================================================

minI=255.0;
maxI=0.0;

[nr,nc]=size(ewt2d);

% Find the min and max over all images
for r=1:nr;
    for c=1:nc;
        minI=min(minI,min(ewt2d{r}{c}(:)));
        maxI=max(maxI,max(ewt2d{r}{c}(:)));
    end
end

% Plot each subband with the same normalization
figure;
for c=1:nc;
    for r=1:nr;
       n=(c-1)*nr+r;
       subplot(nc,nr,n);imshow(ewt2d{r}{c},[minI maxI]);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
%% QUAD_BATCH_LOCAL uses the BATCH command to run the QUAD code locally.
%
%  Discussion:
%
%    The QUAD code is a function, so first we must write a script
%    called QUAD_SCRIPT that runs the function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2010
%
%  Author:
%
%    John Burkardt
%
  clear

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUAD_BATCH_LOCAL\n' );
  fprintf ( 1, '  Run QUAD_SCRIPT locally.\n' );
%
%  BATCH sends the script for execution.
%
  job = batch ( 'quad_script', ...
    'Configuration', 'local', ...
    'CaptureDiary', true, ...
    'FileDependencies', { 'quad_fun' }, ...
    'matlabpool', 7 );
%
%  WAIT pauses MATLAB session til the job completes.
%
  wait ( job );
%
%  DIARY displays any messages printed during execution.
%
  diary ( job );
%
%  LOAD makes the script's workspace available.
%
%  value = estimated integral.
%
  load ( job );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimated integral = %f\n', value );
  fprintf ( 1, '  Exact integral    is %f\n', pi );
  fprintf ( 1, '  Error is             %e\n', abs ( value - pi ) );
%
%  DESTROY cleans up data about the job we no longer need.
%
  destroy ( job );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUAD_BATCH_LOCAL\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

author: johannesgerer

repo: jburkardt-m


text:
function [xmin,fmin]=Trust_Region_Method(x0)

eta1=0.01;
eta2=0.75;
gamma1=0.5;
gamma2=2.0;
Delta_bar=3;
Delta0=1;
epsilon=1e-6;
k=0;
xk=x0;
Deltak=Delta0;

while 1
    % stop criteria
    gk=gradient(xk);
    Bk=Hessian(xk);
    if norm(gk)<=epsilon
        xmin=xk;
        fmin=objective_fun(xmin);
        k
        return;
    end

    % solve subproblem
    % Cauchy displacement
    if gk'*Bk*gk<=0
        dk=-(Deltak/norm(gk))*gk;
    else
        dk=-min(norm(gk)^2/(gk'*Bk*gk),Deltak/norm(gk))*gk;
    end

    % Dogleg method
%     dku=-((norm(gk)^2)/(gk'*Bk*gk))*gk;
%     dkN=-inv(Bk)*gk;
%     if norm(dku)>=Deltak
%         dk=-(Deltak/norm(gk))*gk;
%     elseif norm(dku)<Deltak && norm(dkN)<=Deltak
%         dk=dkN;
%     else
%         temp1=inf;
%         for tau=0:0.01:1
%             if norm((1-tau)*dku+tau*dkN)-Deltak<temp1
%                 dk=(1-tau)*dku+tau*dkN;
%             end
%         end
%     end

    % ratio
    rhok=(objective_fun(xk)-objective_fun(xk+dk))/(mk(zeros(length(xk),1),xk)-mk(dk,xk));

    % adjust trust region radius
    if rhok<=eta1
        Deltak=gamma1*Deltak;
    elseif rhok>=eta2 && (norm(dk)<=Deltak+epsilon && norm(dk)>=Deltak-epsilon)
        Deltak=min(gamma2*Deltak,Delta_bar);
    end

    % update iteration point
    if rhok>eta1
        xk=xk+dk;
    end
    k=k+1;
end

repo: Optimization-Theory-and-Algorithm


text:
function spm_dartel_kernel(job)
% Generate Fisher kernel from flow fields
% FORMAT spm_dartel_kernel(job)
% job.flowfields - Flow-fields
% job.rform      - Form of L
% job.rparam     - Parameters of L
% job.dotprod    - Part of filename for results
%
% k(x_1,x_2) = <x_1,L x_2> = <L x_1, x_2>
%
% This is very slow, and is not in a form that would be
% suited to weighting according to location in the image.
% For this, the "square root" of L would need to be used
% in order to convert the flow fields into (e.g.) their
% Jacobian tensor fields.  For linear elasticity, this
% field would be decomposed by J = (J+J')/2 + (J-J')/2.
% The elements of the symetric part (along with its trace)
% would then be used to generate the kernel.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_dartel_kernel.m 4492 2011-09-16 12:11:09Z guillaume $

P      = strvcat(job.flowfields);
[pth,nam,ext] = fileparts(job.dotprod);
ofname = fullfile(pwd,['dp_' nam '.mat']);
form   = job.rform;
param  = job.rparam;

prm = [form 1 1 1 param];

N   = nifti(P);
dm  = size(N(1).dat);
n   = numel(N);
Phi = zeros(n,n);
spm_progress_bar('Init',n*n,'Generating kernel','Elements done');
for i=1:n,
    x1 = single(squeeze(N(i).dat(:,:,:,end,:)));
    x1 = dartel3('vel2mom',x1,prm);
    for j=i:n,
        x2       =squeeze(N(j).dat(:,:,:,end,:,:));
        d        = x1(:)'*x2(:);
        Phi(i,j) = d;
        Phi(j,i) = d;

        if ~rem(j,8)
            spm_progress_bar('Set',...
                n*n-(n+1-i)*(n+1-i)+(j-i)*2+1);
        end
    end
    input = job;
    typ   = 'flow1';
    save(ofname,'Phi','input','typ', spm_get_defaults('mat.format'));
end
spm_progress_bar('Clear');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ r, seed ] = r4_uniform_01 ( seed )

%*****************************************************************************80
%
%% R4_UNIFORM_01 is a uniform random number generator.
%
%  Discussion:
%
%    This routine implements the recursion
%
%      seed = 16807 * seed mod ( 2**31 - 1 )
%      r4_uniform_01 = seed / ( 2**31 - 1 )
%
%    The integer arithmetic never requires more than 32 bits,
%    including a sign bit.
%
%    If the initial seed is 12345, then the first three computations are
%
%      Input     Output      R4_UNIFORM_01
%      SEED      SEED
%
%         12345   207482415  0.096616
%     207482415  1790989824  0.833995
%    1790989824  2035175616  0.947702
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 July 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley Interscience, page 95, 1998.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer SEED, the integer "seed" used to generate
%    the output random number.  SEED should not be 0.
%
%    Output, real R, a random value between 0 and 1.
%
%    Output, integer SEED, the updated seed.  This would
%    normally be used as the input seed on the next call.
%
  seed = floor ( seed );

  seed = mod ( seed, 2147483647 );

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

  k = floor ( seed / 127773 );

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

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

  r = seed * 4.656612875E-10;

  return
end


text:
function ll = svargplvmLogLikelihood(model)

% SVARGPLVMLOGLIKELIHOOD Log-likelihood for a shared variational GP-LVM.
% FORMAT
% DESC returns the log likelihood for a given SVARGP-LVM model.
% ARG model : the model for which the log likelihood is to be
% computed. The model contains the data for which the likelihood is
% being computed in the 'y' component of the structure.
% RETURN ll : the log likelihood of the data given the model.
%
% SEEALSO:  svargplvmLogLikeGradients
%
% COPYRIGHT : Andreas C. Damianou, 2011

% VARGPLVM

% Functions f1 and f2 should be equivalent in the static case, but f1 might
% (?) be faster.
if ~isfield(model, 'dynamics') || isempty(model.dynamics)
    ll = f1(model);
    %ll = f2(model);
else
    ll = f2(model);
    for i = 1:model.numModels
        if isfield(model.comp{i}, 'KLweight')
            assert(model.comp{i}.KLweight == 0.5); % not implemented yet for dynamics
        end
    end
end

function f = f1(model)

% model.numModels=1; %%%%%%%%TEMP

% This works only when there are NOT any dynamics:
varmeans = sum(sum(model.vardist.means.*model.vardist.means));
varcovs = sum(sum(model.vardist.covars - log(model.vardist.covars)));
KLdiv = -0.5*(varmeans + varcovs) + 0.5*model.q*model.N;

model = svargplvmPropagateField(model, 'onlyLikelihood', 1);
ll=0;
for i=1:model.numModels
    ll = ll + vargplvmLogLikelihood(model.comp{i});
end
f = (ll + KLdiv);

% Here we take advantage of the fact that in the bound, the likelihood and
% the KL part break into ll+KL. So, in the shared case, we only have
% KL + ll1 + ll2 + ...
function f = f2(model)

% We want to count for the KL part only once, because it is shared. Then, we
% will add all likelihood parts from all sub-models. The above can be
% achieved by calculating the bound for the first sub-model as usual, and
% then calculate only the likelihood part for the rest of the sub-models
% (because the KL part would be the same for all models since we called
% expandParam previously).
f = vargplvmLogLikelihood(model.comp{1}); % ll_1 + KL_1 (where KL_1 == KL_i for all i)
% The following call will not affect the model after returning from this
% function, since the model is not returned.
model = svargplvmPropagateField(model, 'onlyLikelihood', 1, true);
ll=0;
for i=2:model.numModels
    ll = ll + vargplvmLogLikelihood(model.comp{i}); % ll_i
end

f = f + ll;

author: SheffieldML

repo: vargplvm


text:
classdef prtPreProcPca < prtPreProc
    % prtPreProcPca   Principle Component Analysis
    %
    %   PCA = prtPreProcPca creates a Principle Component Analysis selfect.
    %
    %   PCA = prtPreProcPca('nComponents',N) constructs a
    %   prtPreProcPCP selfect PCA with nComponents set to the value N.
    %
    %   A prtPreProcPca selfect has the following properites:
    %
    %   nComponents    - The number of principle componenets
    %
    %   A prtPreProcPca selfect also inherits all properties and functions from
    %   the prtAction class
    %
    %   Example:
    %
    %   dataSet = prtDataGenFeatureSelection;    % Load a data set
    %   pca = prtPreProcPca;            % Create a prtPreProcPca selfect
    %
    %   pca = pca.train(dataSet);       % Train the prtPreProcPca selfect
    %   dataSetNew = pca.run(dataSet);  % Run
    %
    %   % Plot
    %   plot(dataSetNew);
    %   title('PCA Projected Data');
    %
    %   See Also: prtPreProc, prtPreProcPca, prtPreProcPls,
    %   prtPreProcHistEq, prtPreProcZeroMeanColumns, prtPreProcLda,
    %   prtPreProcZeroMeanRows, prtPreProcLogDisc, prtPreProcZmuv,
    %   prtPreProcMinMaxRows

    % Copyright (c) 2013 New Folder Consulting
    %
    % Permission is hereby granted, free of charge, to any person obtaining a
    % copy of this software and associated documentation files (the
    % "Software"), to deal in the Software without restriction, including
    % without limitation the rights to use, copy, modify, merge, publish,
    % distribute, sublicense, and/or sell copies of the Software, and to permit
    % persons to whom the Software is furnished to do so, subject to the
    % following conditions:
    %
    % The above copyright notice and this permission notice shall be included
    % in all copies or substantial portions of the Software.
    %
    % THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
    % OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
    % MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
    % NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
    % DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
    % OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
    % USE OR OTHER DEALINGS IN THE SOFTWARE.

    properties (SetAccess=private)
        name = 'Principal Component Analysis' % Principal Component Analysis
        nameAbbreviation = 'PCA'  % PCA
    end

    properties
        nComponents = 3;   % The number of principle components
    end
    properties

        means = [];           % A vector of the means
        pcaVectors = [];      % The PCA vectors.

        trainingTotalVariance = []; % The total variance contained in the
        % training data
        totalVariance = []; % The variance contained in the reduced
        % dimension data.
        totalVarianceCumulative = []; % The variance contained in the
        % reduced dimension data as a
        % function of the number of
        % components
        totalPercentVarianceCumulative = []; %The perceont of the total training variance explained in totalVarianceCumulative
    end
    properties (Hidden)
        removeMean = true;
    end

    methods

        % Allow for string, value pairs
        function self = prtPreProcPca(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end
    end

    methods

        function dataSet = approximate(self,dataSet)
            % dataSetApp = approximate(pca,dataSet)
            %  Generate the PCA basis approximation to the data in dataSet.
            %
            dataSetScores = self.run(dataSet);
            dataSet = reconstruct(self,dataSetScores);
        end

        function dataSet = reconstruct(self,dataSetScores)
            % dataSetApp = reconstruct(pca,dataSetScores)
            %  Generate the PCA basis approximation using the scores in dataSetScores
            %

            xOut = repmat(self.means,dataSetScores.nObservations,1);
            if self.nComponents > 0
                xOut = xOut + (dataSetScores.X*self.pcaVectors');
            end
            dataSet = dataSetScores;
            dataSet.X = xOut;
        end

        function self = set.nComponents(self,nComp)
            if ~isnumeric(nComp) || ~isscalar(nComp) || nComp <= 0 || round(nComp) ~= nComp
                error('prt:prtPreProcPca','nComponents (%s) must be a positive scalar integer',mat2str(nComp));
            end
            self.nComponents = nComp;
        end

    end

    methods (Hidden = true)
        function featureNameModificationFunction = getFeatureNameModificationFunction(self) %#ok<MANU>
            featureNameModificationFunction = prtUtilFeatureNameModificationFunctionHandleCreator('PC Score #index#');
        end
    end

    methods (Access = protected, Hidden = true)
        function self = trainAction(self,dataSet)

            if self.removeMean
                self.means = prtUtilNanMean(dataSet.getObservations(),1);
            else
                self.means = zeros(1,dataSet.nFeatures);
            end

            if self.nComponents == 0
                return
            end
            x = bsxfun(@minus,dataSet.getObservations(),self.means);

            maxComponents = min(size(x));

            if self.nComponents > maxComponents
                warning('prt:prtPreProcPca','User specified # PCA components (%d) is > maximum number of PCA allowed (min(size(dataSet.data)) = %d)',self.nComponents,maxComponents);
                self.nComponents = maxComponents;
            end

            [s,u,v] = svds(x,self.nComponents); %#ok<ASGLU>

            self.pcaVectors = v;

            self.trainingTotalVariance = sum(var(x));
            pcaVariance = cumsum(var(x*v));

            self.totalVarianceCumulative = pcaVariance;
            self.totalVariance = self.totalVarianceCumulative(end);
            self.totalPercentVarianceCumulative = self.totalVarianceCumulative./self.trainingTotalVariance;
        end

        function dataSet = runAction(self,dataSet)
            dataSet.X = self.runActionFast(dataSet.getObservations);
        end

        function xOut = runActionFast(self,xIn)

            if self.nComponents == 0
                xOut = nan(size(xIn,1),1);
                return;
            end
            if self.removeMean
                xOut = bsxfun(@minus,xIn,self.means);
                if self.nComponents > 0
                    xOut = xOut*self.pcaVectors;
                end
            else
                xOut = xIn*self.pcaVectors;
            end
        end

    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function M = oned_bilinear ( kernel, phi, test, w_g )

%*****************************************************************************80
%
%% ONED_BILINEAR integrates a kernel times a basis times a test function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 December 2012
%
%  Author:
%
%    Jeff Borggaard
%
%  Parameters:
%
%    Input, real KERNEL(N_GAUSS), the kernel function evaluated
%    at the Gauss points.
%
%    Input, real PHI(N_GAUSS,N_DOF), the element test functions evaluated
%    at the Gauss points.
%
%    Input, real TEST(N_GAUSS,N_TEST), the test functions evaluated
%    at the Gauss points.
%
%    Input, real W_G(N_GAUSS), the Gauss weights.
%
%    Output, real M(N_TEST,N_DOF), the integral values.
%
  M = test' * diag ( kernel .* w_g ) * phi;

  return
end


text:
%	OVERVIEW:
%       This demonstration analyzes a segment of data collected in the
%       intensive care unit (ICU) which contains ECG, ABP, and PPG signals
%
%   OUTPUT:
%       HRV Metrics exported to .cvs files
%
%   DEPENDENCIES & LIBRARIES:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   REFERENCE:
%       Vest et al. "An Open Source Benchmarked HRV Toolbox for Cardiovascular
%       Waveform and Interval Analysis" Physiological Measurement (In Press), 2018.
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Giulia Da Poian
%	COPYRIGHT (C) 2018
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; clc; close all;

run(['..' filesep 'startup.m'])

% Remove old files generated by this demo
OldFolder = [pwd,filesep, 'OutputData', filesep, 'ResultsICU'];
if exist(OldFolder, 'dir')
    rmdir(OldFolder, 's');
    fprintf('Old Demo Folder deleted \n');
end

HRVparams = InitializeHRVparams('demoICU'); % include the project name
HRVparams.poincare.on = 0; % Pinocare analysis off for this demo
HRVparams.DFA.on = 0; % DFA analysis off for this demo
HRVparams.MSE.on = 0; % MSE analysis off for this demo
HRVparams.HRT.on = 0; % HRT analysis off for this demo

[subjectIDs,filesTBA] = GenerateListOfFilesTBA(HRVparams.ext,HRVparams.readdata,0);

idx = find(strcmp(subjectIDs,'TestICUdata'));

i_patient = idx;

% 1. Load Raw Patient Data (ECG Waveform)
load(filesTBA{i_patient});

% 2. Analyze data using HRV PhysioNet Cardiovascular Signal Toolbox
[~, resFilename] = Main_HRV_Analysis(signal(:,1),[],'ECGWaveform',...
    HRVparams,subjectIDs(i_patient),[],[],signal(:,5),'PPG',signal(:,3),'ABP');

% 3. Load annotations ans SQI for plot
AnnFile = strcat(HRVparams.writedata, filesep, 'Annotation', filesep, ...
                                                    subjectIDs{i_patient});
jqrs_ann = read_ann(AnnFile,'jqrs');
[~,~,sqijw] = read_ann(AnnFile,'sqijw'); % SQI is x100 to be stored

ppg_ann = read_ann(AnnFile,'ppg');
qppg(signal(:,5),HRVparams.Fs);
[~,ppgsqi,ppgNumSqi] = read_ann(AnnFile,'sqippg');

abpann = read_ann(AnnFile,'abpm');
features =  abpfeature(signal(:,3), abpann, HRVparams.Fs);
[BeatQ, goodbeats] = jSQI(features, abpann, signal(:,3));

% 5. Plotting
HRVparams.gen_figs = 1;
if HRVparams.gen_figs
    % ECG
    Plot_SignalDetection_SQI(time, signal(:,1), jqrs_ann, sqijw'./100,'ECG')

    % ABP
    Plot_SignalDetection_SQI(time, signal(:,3), abpann, double(~BeatQ(:,1))','ABP')

    % RESP
    figure(88)
    plot(time,signal(:,6));
    title('Respiration')

    % PLETH
    Plot_SignalDetection_SQI(time, signal(:,5), ppg_ann, ppgNumSqi'./100,'PPG')

end

% 6. Pulse Transit Time
ptt = pulsetransit(jqrs_ann, abpann);

% 7. Plot BP vs PTT
syst = features(1:length(ptt),2);
if HRVparams.gen_figs
    figure(22)
    plot(syst,ptt(:,3)./HRVparams.Fs,'o');
    xlabel('BP (mmHg)'); ylabel('PTT (s)');
    title('Pulse Transit Time - BP vs PTT (ABP - QRS)')

end

% 8. Compare generated output file with the reference one

currentFile = strcat(HRVparams.writedata, filesep, resFilename.HRV, '.csv');
referenceFile = strcat('ReferenceOutput', filesep, 'ICU_HRV_allwindows.csv');
testHRV = CompareOutput(currentFile,referenceFile);

if testHRV
    fprintf('** DemoRawDataICU: TEST SUCCEEDED ** \n ')
     fprintf('A file named %s.csv \n has been saved in %s \n', ...
    resFilename.HRV, HRVparams.writedata);
else
    fprintf('** DemoRawDataICU: TEST FAILED ** \n')
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [ cc_est, r_min, r_max ] = cube3d_grid_centralize ( ng, oc )

%*****************************************************************************80
%
%% CUBE3D_GRID_CENTRALIZE centralizes grid data from the surface of a 3D cube.
%
%  Discussion:
%
%    We generate a hypersphere grid of 1D index NG, with origin at OC,
%    and project these points onto the 3D cube of center [0,0,0]
%    and radius 1.
%
%    The cube is aligned with the coordinate axes.
%
%    We seek to estimate the value of the center.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NG, the grid index.  The hypersphere grid will have
%    order N = (2*NG+1) * (NG+1)^(M-2).
%
%    Input, real OC(M,1), the observation point, the original of the
%    hypersphere coordinate system.  OC should be inside the cube.
%
%    Output, real CC_EST(M,1), the estimated cube center.
%
%    Output, real R_MIN, R_MAX, the minimum and maximum distances
%    between the observation point and sample points on the cube surface.
%
  m = 3;
  n = ( 2 * ng + 1 ) * ( ng + 1 ) ^ ( m - 2 );
%
%  Destroy all row vectors!
%
  oc = oc(:);
%
%  Compute sample points on the cube surface.
%
  [ x1, x2, x3 ] = cube3d_grid ( ng, oc );
%
%  "Flatten" the data into vectors.
%
  x1 = reshape ( x1, n, 1 );
  x2 = reshape ( x2, n, 1 );
  x3 = reshape ( x3, n, 1 );
%
%  The center estimate is simply the average.
%
  cc_est = zeros ( 3, 1 );
  cc_est(1,1) = sum ( x1 ) / n;
  cc_est(2,1) = sum ( x2 ) / n;
  cc_est(3,1) = sum ( x3 ) / n;
%
%  For each point on the cube surface, the "radius" ROC is the norm of X - OC.
%
  roc = [ x1, x2, x3 ]' - repmat ( oc, 1, n );
  roc = roc.^2;
  roc = sum ( roc, 1 );
  roc = sqrt ( roc );
%
%  Return the minimum and maximum radiuses.
%
  r_min = min ( roc );
  r_max = max ( roc );

  return
end


text:
function createVolAnatSKERI(imgPathName,outputVAnatName,outputUnfoldName, clipVals)
% createVolAnat([imgPathName],[outputVAnatName],[outputUnfoldName], ...
% [clipVals=get from UI])
%
% imgPathName is the path to one of the I-files or to an analyze-format file.
% A file browser comes up if none is given.
%
% If the optional clipVals (in the format of [minVal maxVal]) is provided,
% will clip the vAnatomy before converting to uint8 (0-255, a big
% annoyance resulting from older software constraints) for the .dat file.
% Otherwise, pops up some example slices and a histogram, and lets you
% select the cutoff.
%
% 2000.01.20 RFD (bob@white.stanford.edu): wrote it.
% 2001.02.21 RFD: renamed from createVanatomy and added improved
%            the clipping UI by checking for obvious cases, like
%            when no clipping is necessary.
% 2002.05.22 RFD: fixed a bug that would apply an incorrect upper clip
%            when the lower clip was not zero. THis probably didn't affect
%            anyone with the old (crappy) anatomy data, since the lower clip
%            was always best left at zero anyway.
%            I also now get the mmPerPix from the image header.
% 2002.07.30 RFD: now allow Analyze format files as input. Reorients standard
%            axial Analyze data to match the mrGray sagittal convention.
% 2002.08.05 RFD: fixed crash when imgPathName is passed that was introduced on 7/30.
% 2003.06.25 DHB (brainard@psych.upenn.edu): Allow passing of output names.
% 2004.01.24 Junjie: change GE_readHeader to readIfileHeader
% Get the input filename
if ~exist('imgPathName','var') | isempty(imgPathName)
    [fname, fullPath] = uigetfile({'*.img;*.nii;*.nii.gz','MRI images';'*.img','analyze';'*.nii','NIFTI';'*.nii.gz','NIFTI gz';'*.*','all files'}, 'Select analyze/NIFTI I-files...');
    if(isnumeric(fname))
        error('user canceled.');
    end
    imgPathName = fullfile(fullPath, fname);
end
if ~exist('outputVAnatName','var')
    outputVAnatName = [];
end
if ~exist('outputUnfoldName','var') | isempty(outputUnfoldName)
    outputUnfoldName = 'UnfoldParams';
end
[fullPath,fname,ext] = fileparts(imgPathName);
if(strcmp(ext,'.hdr') | strcmp(ext,'.img'))
    [img,mmPerPix] = loadAnalyze(fullfile(fullPath,[fname '.hdr']));
    % To make analyze data the same orientation as vAnatomy, we reorient the data.
    % Note that we assume Analyze orientation code '0'- transverse, unflipped.
    % If that is the case, then the following will make the data look correct in mrGray.
    % It will even get left/right correct (left will be left and right right).
    img = permute(img,[3,2,1]);
    mmPerPix = [mmPerPix(3),mmPerPix(2),mmPerPix(1)];
    % flip each slice ud (ie. flip along matlab's first dimension, which is our x-axis)
    for jj=1:size(img,3)
        img(:,:,jj) = flipud(squeeze(img(:,:,jj)));
    end
    % flip each slice lr(ie. flip along matlab's second dimension, which is our y-axis)
    for jj=1:size(img,3)
        img(:,:,jj) = fliplr(squeeze(img(:,:,jj)));
    end
elseif any(strcmpi(ext,{'.nii','.gz'}))
% 	nii = load_nii(fullfile(fullPath,[fname,'.nii']));		% old_RGB flag only matters if datatype=128 & bitpix=24
% 	if isfield(nii.hdr.hist,'rot_orient')		% already RAS?
% 		img = nii.img;
% 		mmPerPix = nii.hdr.dime.pixdim(2:4);
% 	else
% 		error('unknown orientation')
% % 		k = [? ? ?];
% % 		img = permute(nii.img,k);
% % 		mmPerPix = nii.hdr.dime.pixdim(k+1);
% 	end
% 	% now orient to IPR for writeVolAnat to yield PIR vAnatomy
% 	img = permute(img,[3 2 1]);
% 	mmPerPix =mmPerPix([3 2 1]);
% 	img = flipdim(img,1);
% 	img = flipdim(img,2);

	nii = niftiRead(imgPathName);
	mmPerPix = nii.pixdim;
	img = double(nii.data);

	% get nifti rotation matrix
	qb = nii.quatern_b;
	qc = nii.quatern_c;
	qd = nii.quatern_d;
	qa = sqrt(1-qb^2-qc^2-qd^2);
	R = [ qa^2+qb^2-qc^2-qd^2, 2*(qb*qc+qa*qd), 2*(qb*qd-qa*qc);...
			2*(qb*qc-qa*qd), qa^2-qb^2+qc^2-qd^2, 2*(qc*qd+qa*qb);...
			2*(qb*qd+qa*qc), 2*(qc*qd-qa*qb), qa^2-qb^2-qc^2+qd^2];
	R = R';	% something about the handedness??? only if positive determinant of R???
				% rows of R now [R;A;S]
	R(:,3) = nii.qfac * R(:,3);
	disp(R)

	% now orient to IPR for writeVolAnat to yield PIR vAnatomy
	IPR = [0 0 0];
	[junk,IPR(3)] = max(abs(R(1,:)));	% L-R dimension
	[junk,IPR(2)] = max(abs(R(2,:)));	% P-A dimension
	[junk,IPR(1)] = max(abs(R(3,:)));	% I-S dimension
	if numel(unique(IPR))~=3
		error('bad rotation matrix')
	end
	if any(abs([R(1,IPR(3)),R(2,IPR(2)),R(3,IPR(1))]) < 0.8)
		disp('nifti quaternion inconsistent with 90deg rotations')
	end
	img = permute(img,IPR);
	if R(3,IPR(1)) >= 0
		img = flipdim(img,1);
	end
	if R(2,IPR(2)) >= 0
		img = flipdim(img,2);
	end
	if R(1,IPR(3)) < 0
		img = flipdim(img,3);
	end

% elseif(strcmp(ext,'.gz') | strcmp(ext,'.nii'))
elseif strcmp(ext,'.gz')
    ni = niftiRead(imgPathName);
    img = permute(double(ni.data),[3,2,1]);
    mmPerPix = [ni.pixdim(3),ni.pixdim(2),ni.pixdim(1)];
    for(jj=1:size(img,3))
        img(:,:,jj) = flipud(squeeze(img(:,:,jj)));
    end
    for(jj=1:size(img,3))
        img(:,:,jj) = fliplr(squeeze(img(:,:,jj)));
    end
else
    % get some info from the I-file header
    [su_hdr,ex_hdr,se_hdr,im_hdr,pix_hdr,im_offset] = readIfileHeader(imgPathName);
    mmPerPix = [im_hdr.pixsize_Y, im_hdr.pixsize_X, im_hdr.slthick];
%    nSlices = se_hdr.se_numimages;
% Now, the DICOM format of i-files does not provide numimages info, so you
% need to count all i-files in the folder.
    nSlices = length(getIfileNames(imgPathName));
    imDimXY = [im_hdr.dim_X, im_hdr.dim_Y];

    disp('Image header says:');
    disp(['  ',num2str(nSlices),' slices']);
    disp(['  ',num2str(imDimXY),' pixels inplane (x,y)']);
    disp(['  ',num2str(mmPerPix),' mm per pixel']);
    % load ifiles
    disp('Loading I-files...');
    img = makeCubeIfiles(fullfile(fullPath,fname), imDimXY, [1:nSlices]);
end
% Get image information
volume_pix_size = 1./mmPerPix;
disp('Finding optimal clip values...');
minVal = min(img(:));
maxVal = max(img(:));

% Set the minumum to zero
if(minVal~=0)
    disp(['Image minimum is not zero (',num2str(minVal),')- adjusting...']);
    img = img-minVal;
	 maxVal = maxVal - minVal;
	 minVal = 0;
end

% clip and scale, if necessary
if (maxVal <= 255)
    disp('Image values are within 8-bit range- no clipping or scaling is necessary.');
	 lowerClip = minVal;
	 upperClip = maxVal;
elseif exist('clipVals', 'var') && ~isempty(clipVals)
    % Clip and scale image values to be 0-255
    disp('Clipping intensities and scaling to 0-255...');
    img(img > clipVals(2)) = clipVals(2);
    img(img < clipVals(1)) = clipVals(1);
    img = round((img-clipVals(1))./(clipVals(2)-clipVals(1)) * 255);
	 lowerClip = clipVals(1);
	 upperClip = clipVals(2);
else
    figure(99);
    disp('Computing histogram...');
    [clippedImg, suggestedClipVals] = mrAnatHistogramClip(img, 0.5, 0.98);
    lowerClip = suggestedClipVals(1);
    upperClip = suggestedClipVals(2);
    [n,x] = hist(img(:),100);
    histHandle = subplot(2,2,1);

%     bar(x,n,1);
%     lcLine = line([lowerClip,lowerClip], [0,max(n)]);
%     ucLine = line([upperClip,upperClip], [0,max(n)]);

% set yscale so you can see more than the largest histogram peak
	bar(x,n,1,'b');
	set(histHandle,'ylim',[-0.05 1.05]*2e5,'xlim',[-0.05 1.05]*max(x))
	cLineY = [-0.05,1.05]*max(n);
	lcLine = line([lowerClip,lowerClip], cLineY, 'color','r');
	ucLine = line([upperClip,upperClip], cLineY, 'color','r');

    % display some slices
    slice{1} = squeeze(img(round(end/2),:,:));
    slice{2} = squeeze(img(:,round(end/2),:));
    slice{3} = squeeze(img(:,:,round(end/2)));
    for(i=[1:3])
        slice{i}(slice{i}<lowerClip) = lowerClip;
        slice{i}(slice{i}>upperClip) = upperClip;
        slice{i} = slice{i}-lowerClip;
        slice{i} = round(slice{i}./upperClip*255);
        subplot(2,2,2);
        image(slice{1});colormap(gray(256));axis image;axis off;
        subplot(2,2,3);
        image(slice{2});colormap(gray(256));axis image;axis off;
        subplot(2,2,4);
        image(slice{3});colormap(gray(256));axis image;axis off;
    end
    % Loop until user says it looks good
    ok = 0;
    while(~ok)
%         answer = inputdlg({'lower clip value (min=0): ',...
%                 ['upper clip value (max=',num2str(maxVal),'):']}, ...
%                 'Set intensity clip values', 1, ...
%                 {num2str(lowerClip),num2str(upperClip)}, 'on');
%         if ~isempty(answer)
%             lowerClip = str2num(answer{1});
%             upperClip = str2num(answer{2});
%         end

        disp(['intensity clip values are: ' num2str(lowerClip) ', ' num2str(upperClip)]);
        infoStr = 'Click twice on the histogram to set clip values. Define the smallest range that keeps the two right peaks.';
        uiwait(helpdlg(infoStr, 'clip info'))
        axes(histHandle);
        [x,y] = ginput(2);
        lowerClip = min(x);
        upperClip = max(x);

%         delete(lcLine); delete(ucLine);
%         lcLine = line([lowerClip,lowerClip], [0,max(n)]);
%         ucLine = line([upperClip,upperClip], [0,max(n)]);
			set(lcLine,'xdata',[lowerClip,lowerClip])
			set(ucLine,'xdata',[upperClip,upperClip])

        % Display some slices
        slice{1} = squeeze(img(round(end/2),:,:));
        slice{2} = squeeze(img(:,round(end/2),:));
        slice{3} = squeeze(img(:,:,round(end/2)));
        for(i=[1:3])
            slice{i} = rescale2(slice{i}, [lowerClip, upperClip]);
%             slice{i}(slice{i}<lowerClip) = lowerClip;
%             slice{i}(slice{i}>upperClip) = upperClip;
%             slice{i} = slice{i}-lowerClip;
%             slice{i} = round(slice{i}./upperClip*255);
            subplot(2,2,2);
            image(slice{1});colormap(gray(256));axis image;axis off;
            subplot(2,2,3);
            image(slice{2});colormap(gray(256));axis image;axis off;
            subplot(2,2,4);
            image(slice{3});colormap(gray(256));axis image;axis off;
        end

        bn = questdlg('Is this OK?','Confirm clip values','Yes','No','Cancel','Yes');
        if(strcmp(bn,'Cancel'))
            error('Cancelled.');
        else
            ok = strcmp(bn,'Yes');
        end
    end

    % Clip and scale image values to be 0-255
    %
    % I haven't been fully satisfied by just automatically clipping off the top
    % few % of intensities (this is what mrInitRet currently does to the inplanes).
    % So, lately I've been looking at the histogram and picking intensityClip
    % by hand, then viewing the images to see how well it worked.
    disp('Clipping intensities and scaling to 0-255...');
    img(img>upperClip) = upperClip;
    img(img<lowerClip) = lowerClip;
    img = round((img-lowerClip)./(upperClip-lowerClip) * 255);
    %img = rescale2(img, [lowerClip, upperClip], [0, 255]);
end
subplot(2,2,2);
image(squeeze(img(round(end/2),:,:)));colormap(gray(256));axis image; set(gca,'xtick',[],'ytick',[]),	xlabel('right \rightarrow'),ylabel('\leftarrow posterior') %axis off;
subplot(2,2,3);
image(squeeze(img(:,round(end/2),:)));colormap(gray(256));axis image; set(gca,'xtick',[],'ytick',[]), xlabel('right \rightarrow'),ylabel('\leftarrow inferior') %axis off;
subplot(2,2,4);
image(squeeze(img(:,:,round(end/2))));colormap(gray(256));axis image; set(gca,'xtick',[],'ytick',[]), xlabel('posterior \rightarrow'),ylabel('\leftarrow inferior') %axis off;
% Crop the image cube
% (you can skip this if you want, and just do the crop in mrGray)
%figure(1);image(squeeze(img2(round(end/2),:,:))./clip.*255);colormap(gray(256));axis image;
%figure(2);image(squeeze(img2(:,round(end/2),:))./clip.*255);colormap(gray(256));axis image;
%figure(3);image(squeeze(img2(:,:,round(end/2)))./clip.*255);colormap(gray(256));axis image;
%img2 = img2(20:200,30:240,1:124);
%figure(1);image(squeeze(img(round(end/2),:,:))./clip.*255);colormap(gray(256));axis image;
%figure(2);image(squeeze(img(:,round(end/2),:))./clip.*255);colormap(gray(256));axis image;
%figure(3);image(squeeze(img(:,:,round(end/2)))./clip.*255);colormap(gray(256));axis image;
% Save
disp('Saving volume anatomy and unfolding params ...');
path = writeVolAnat(img,mmPerPix,outputVAnatName);
save([path outputUnfoldName], 'volume_pix_size', 'lowerClip', 'upperClip');
disp(['vAnatomy and UnfoldParams saved to ' path]);
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [imageSeg]= getPlanSeg_complete(XYZworldframeTest,SpaceTest,imgSize,imageNum)
% input XYZworldframeTest is alinged in Z
onPlaneThreshold =0.055;
sizethr =750;
normalAgreeThreshold =0.8;

removeNaN = find(~sum(isnan(XYZworldframeTest),2));
XYZworldframeTest = XYZworldframeTest(removeNaN,:);
[zhist,zc]=hist(XYZworldframeTest(:,3),round(20*(SpaceTest.Rz(2)-SpaceTest.Rz(1))));
[~,zpks] =findpeaks([0,zhist,0],'minpeakdistance',3);
zpks = zpks-1;
% figure,
% bar(bincenter,count);
% hold on;
% plot(bincenter(pks),count(pks),'rx');
% remove
range =[max(min(XYZworldframeTest(:,3)),zc(max(1,zpks-1))'),...
        min(zc(min(size(zc,2),zpks+1)),max(XYZworldframeTest(:,3)))';
        zc(end-2),max(XYZworldframeTest(:,3));
        min(XYZworldframeTest(:,3)),zc(3)];
removepts = inrange(XYZworldframeTest(:,3),range);
XYZworldframeTestRemove=XYZworldframeTest;
XYZworldframeTestRemove(removepts>0,:) =[];

% hough lines
[H,T,R] = votelines(XYZworldframeTestRemove(:,1),XYZworldframeTestRemove(:,2));
% blur it
allangleVote = max(H);
allangleVote = [allangleVote,allangleVote];
shift= 0;
clear rightAngle
for i =1:length(shift)
    rightAngle(i,:) = allangleVote(91+shift(i):180+shift(i))+allangleVote(1:90);
end
[~,AnglepickedLin]= max(rightAngle(:));
[s,Anglepicked]=ind2sub(size(rightAngle),AnglepickedLin);
[~,Rpicked1] = max(H(:,Anglepicked));
[~,Rpicked2] = max(H(:,Anglepicked+90+shift(s)));
P = [R(Rpicked1),T(Anglepicked);R(Rpicked2),T(Anglepicked)+90+shift(s)];

rotation = min(P(:,2));
Rot = getRotationMatrix('z',-1*rotation/180*pi);
XYZworldframeTestNew = [[Rot(1:2,1:2)*XYZworldframeTest(:,[1,2])']', XYZworldframeTest(:,3)];
normals = points2normals(XYZworldframeTestNew);

% find  corner
[xhist,xc]=hist(XYZworldframeTestNew(:,1),round((SpaceTest.Rx(2)-SpaceTest.Rx(1))/0.05));
[~,xpks] =findpeaks([0,xhist,0],'MINPEAKHEIGHT',50);
xpks = xpks-1;
[yhist,yc]=hist(XYZworldframeTestNew(:,2),round((SpaceTest.Rx(2)-SpaceTest.Rx(1))/0.05));
[~,ypks] =findpeaks([0,yhist,0],'MINPEAKHEIGHT',50);
ypks = ypks-1;
[zhist,zc]=hist(XYZworldframeTestNew(:,3),round((SpaceTest.Rx(2)-SpaceTest.Rx(1))/0.05));
[~,zpks] =findpeaks([0,zhist,0],'MINPEAKHEIGHT',50);
zpks = zpks-1;

%% project points onto this plan and caculate conected component
gid =1;
Segment ={};
SegmentSize = [];
imageSegFill = zeros(imgSize);
for p =1:size(xpks,2)
    inLine =find(abs(XYZworldframeTestNew(:,1)-xc(xpks(p)))<onPlaneThreshold&abs(normals(1,:)')>normalAgreeThreshold);
    conthisSegMask = zeros(imgSize);
    conthisSegMask(removeNaN(inLine)) =1;
    label = bwlabel(conthisSegMask,8);
    unique_label = unique(label);
    unique_label(unique_label==0)=[];
    for i =1:length(unique_label),
         if sum(label(:)==unique_label(i))>sizethr,
            Segment(gid) = {label==unique_label(i)};
            SegmentSize(gid) = sum(label(:)==unique_label(i));
            imageSegFill(label==unique_label(i))=gid;
            gid = gid+1;
         end
    end
end

for p =1:size(ypks,2)
    inLine =find(abs(XYZworldframeTestNew(:,2)-yc(ypks(p)))<onPlaneThreshold&abs(normals(2,:)')>normalAgreeThreshold);
    conthisSegMask = zeros(imgSize);
    conthisSegMask(removeNaN(inLine)) =1;
    label = bwlabel(conthisSegMask,8);
    unique_label = unique(label);
    unique_label(unique_label==0)=[];
    for i =1:length(unique_label),
         if sum(label(:)==unique_label(i))>sizethr,

            Segment(gid) = {label==unique_label(i)};
            SegmentSize(gid) = sum(label(:)==unique_label(i));
            imageSegFill(label==unique_label(i))=gid;
            gid = gid+1;
         end
    end
end

floorSize = 0;
floorMask = zeros(imgSize);
Lowest = nanmin(XYZworldframeTestNew(:,3));
for p =1:size(zpks,2)
    inLine =find(abs(XYZworldframeTestNew(:,3)-zc(zpks(p)))<onPlaneThreshold&abs(normals(3,:)')>normalAgreeThreshold);
    conthisSegMask = zeros(imgSize);
    conthisSegMask(removeNaN(inLine)) =1;
    label = bwlabel(conthisSegMask,8);
    unique_label = unique(label);
    unique_label(unique_label==0)=[];
    for i =1:length(unique_label),
         if sum(label(:)==unique_label(i))>sizethr,
            Segment(gid) = {label==unique_label(i)};
            SegmentSize(gid) = sum(label(:)==unique_label(i));
            if p<2&&SegmentSize(gid)>floorSize&&zc(zpks(p))-Lowest<0.25,
                floorSize = SegmentSize(gid);
                floorMask = label==unique_label(i);
            end
            imageSegFill(label==unique_label(i))=gid;
            gid = gid+1;
         end
    end
end
%%

[SegmentSize,segid] = sort(SegmentSize,'descend');
Segment = Segment(segid);
cnt =1;
imageSeg = zeros(imgSize);
for i =1:length(Segment)
    idx = Segment{i}&imageSeg==0;
    se = strel('disk',3);
    idx2 = imclose(idx,se);
    idx(imageSegFill==0&(idx2&~idx))=1;
    if sum(idx(:))>sizethr
       imageSeg(idx)=cnt;
       cnt =cnt+1;
    end
end

while true,
    % remaining points
    remaining_points =imageSeg==0;
    cntPre = cnt;
    label = bwlabel(remaining_points,8);
    unique_label = unique(label);
    unique_label(unique_label==0)=[];
    directions = icosahedron2sphere(0);
    directions = directions(directions(:,2) <= 0,:);
    for i =1:length(unique_label)
        remaining_points_conn = label==unique_label(i);
        remaining_points_conn_ptsidx = find(remaining_points_conn(removeNaN));
        XYZworldframeTestRemaining_conn = XYZworldframeTest(remaining_points_conn_ptsidx,:);
        normals_conn = normals(:,remaining_points_conn_ptsidx);
        if size(XYZworldframeTestRemaining_conn,1)>sizethr
            [~, idx] = max(abs(directions * normals_conn),[],1);
            majorDirIdx = mode(idx);
            normalAgreeId = find(idx ==majorDirIdx);
            [B,~,inliers] = ransacfitplane(XYZworldframeTestRemaining_conn(normalAgreeId,:)', onPlaneThreshold);
            if sum(inliers)>sizethr
                thisSeg = zeros(imgSize);
                thisSeg(removeNaN(remaining_points_conn_ptsidx(normalAgreeId(inliers)))) =1;
                thisSeg_label = bwlabel(thisSeg,8);
                thisSeg_unique_label = unique(label);
                thisSeg_unique_label(thisSeg_unique_label==0)=[];
                for j = 1:length(thisSeg_unique_label),
                    if sum(thisSeg_label(:)==thisSeg_unique_label(j))>sizethr,
                        idxSeg = thisSeg_label==thisSeg_unique_label(j);
                        idxSeg = imclose(idxSeg,se);
                        imageSeg(idxSeg) =cnt;
                        cnt=cnt+1;
                    end
                end

            end
        end
    end
    if cntPre==cnt,
        break;
    end
end
map= [0,randperm(max(imageSeg(:)+1))];
map(map==1) =[];
imageSeg = map(imageSeg+1);
if sum(floorMask(:))>0,
    floorMask = imclose(floorMask,se);
    imageSeg(floorMask) = 1;
end
figure(1),imagesc(imageSeg);
figure(2),imagesc(floorMask);
%Boudary = [D1 D2 D3 D4];
if imageNum> 0,
    im = getImagesc(imageSeg);
    mkdir(segpath);
    imwrite(im,sprintf('%s/%04d.jpg',segpath,imageNum));

    im = getImagesc(floorMask);
    imwrite(im,sprintf('%s/%04d_floor.jpg',segpath,imageNum));

    save(sprintf('%s/%04d.mat',segpath,imageNum),'imageSeg','floorMask')

end
end

function [hough_transform,T,R] = votelines(X,Y)
         thetaResolution = 1;
         rhoResolution = 0.1;
         T = [1:180];
         theta = T/180*pi;
         rho = X(:)*cos(theta)+ Y(:)*sin(theta);
         % quantize the rho
         rhoNorm = max(1,round((rho-min(rho(:)))/rhoResolution));
         R =[1:max(rhoNorm(:))]*rhoResolution+min(rho(:));
         hough_transform = zeros(max(rhoNorm(:)),size(T,2));
         TT = repmat(T,[size(rhoNorm,1),1]);
         hough_transform = accumarray([rhoNorm(:),TT(:)],1,size(hough_transform));
end

author: thusiyuan


text:
classdef X2kKernel < Kernel
    % Object of type kernel but optimized for log so that the computation of the
    % radial quadrature goes faster
    properties (Access = public)
        C = 1;
        k = 1; % such that G(r) =C*r^(2*k)
    end

    methods
      	function[kernel] = X2kKernel(kk,CC)
            if nargin == 0
                kk = 1;
            end
            if nargin <= 1
                CC = 1;
            end
            kernel@Kernel(@(x)(CC*x.^(2*kk)),@(x)(2*kk*CC*x.^(2*kk - 1)))
            kernel.scalFunc = @(a,b,rho)(CC*X2kSP([a,b],rho,kk));
            kernel.normFunc = @(a,b)(abs(CC)*sqrt(2*pi*4*kk^2/(4*kk - 1)*(b^(4*kk - 1) - a^(4*kk- 1))));
            kernel.gamma_est = @X2kKernel.gamma_est;
            kernel.k = kk;
            kernel.C = CC;
            kernel.singular = false;
        end
    end
    methods (Access = public)
        function[out] = dilatation(this,lambda)
            out = X2kKernel(this.k,lambda^(2*this.k));
            % No changes in gamma_est, scalFunc nor normFunc
        end
        function[out] = mtimes(this,mu)
            if isa(this,'Kernel')
                assert(and(isa(mu,'double'),isscalar(mu)));
                out = X2kKernel(this.k,mu*this.C);
            else
                out = mtimes(mu,this);
            end
        end
        function[out] = radialQuadKernel(this,a,tol,varargin)
            Cmem = this.C;
            this = (1/Cmem)*this;
            rq = RadialQuadrature(a,this,tol/abs(Cmem),varargin{:});
            out = Cmem*rq;
        end
    end
    methods (Access = protected)

    end
    methods (Static, Access = protected)
        function[low,up] = gamma_est(~)
            % Helps the radial quadrature to guess the number of components
            up = 4;
            low= 0;
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [scores model binVals bins] = loglossboostLearnRandomFeatures(...
  data,labels,numIters,initWt,binVals,bins,params,obj,str)

if nargin<8,
  obj = [];
end

if nargin<9,
  str = '';
end

numEx = size(data,1);
wt = initWt;
model = struct('dim',{},'error',{},'dir',{},'tr',{},'alpha',{});

% Initialize with dummy models.
dummywk = struct('dim',1,'error',0.5','dir',1,'tr',0,'alpha',0);
for ndx = 1:numIters
  model(ndx) = dummywk;
end
scores = zeros(numEx,1);

clk = tic;
etimehist = [];
etimehistlength = 25;
nittoutput = 3;
origbins = bins; origdata = data; origbinVals = binVals;

numFeatures = size(origbins,1);
if numFeatures<10000,
  featureSetRatio = 1.1;
  refreshRatio = 100000; % never refresh
elseif numFeatures < 20000
  featureSetRatio = 0.7;
  refreshRatio = 1/10;
else
  featureSetRatio = 0.4;
  refreshRatio = 1/20;
end

% params.numSample = min(params.numSample,round(numEx/2));
for itt = 1:numIters
  if (ceil(itt/refreshRatio/numIters)-ceil( (itt-1)/refreshRatio/numIters) > 0.5)
    % every 1/20th iterations select features to train on.
    selFeatures = rand(1,size(origbins,1))>featureSetRatio;
    feature_map = find(~selFeatures);
    % MAYANK NOV 4 2015. copying and deleting as below is slow
%     bins = origbins;
%     data = origdata;
%     binVals = origbinVals;
%     bins(selFeatures,:) = [];
%     data(:,selFeatures) = [];
%     binVals(:,selFeatures) = [];
    bins = origbins(~selFeatures,:);
    data = origdata(:,~selFeatures);
    binVals = origbinVals(:,~selFeatures);
  end
  count = 0;
  while(count<1)
    [wkRule,wksel] = findWeakRuleSamples(data,labels,wt,binVals,bins,params);
    sel = false(size(labels));
    sel(wksel) = true;
    wkRule.dim = feature_map(wkRule.dim);
    tr = wkRule.tr;
    dir = wkRule.dir;
    dim = wkRule.dim;
    if dir>0,
      tt = ((origdata(:,dim)> tr)-0.5)*2;
    else
      tt = ((origdata(:,dim)<= tr)-0.5)*2;
    end
    curError = sum( (tt.*labels(:)).*wt(:));
    if curError>0,
      break;
    else
      curError = - curError;
      wkRule.dir = -wkRule.dir;
    end
    count = count + 1;
  end
%   if count == 11,
%     fprintf('Too much training\n');
%     break;
%   end

  wkRule.error = 0.5-curError/2;
  wkRule.alpha = 1-2*wkRule.error;
  model(itt) = wkRule;
  scores = myBoostClassify(origdata,model(1:itt));
  tt = scores.*labels;
  wt = initWt./(1+exp(tt));
  wt = wt./sum(wt);

  if( nargin>7 && mod(itt,nittoutput)==0)
    etime = toc(clk);
    etimehist(end+1) = etime;
    if numel(etimehist) > etimehistlength,
      etimehist = etimehist(end-etimehistlength+1:end);
    end
    if ~isempty(obj)
      obj.SetStatus('%s %d%% training done. Time Remaining:%ds ',str,...
        round(itt/numIters*100),round((numIters-itt)/nittoutput*mean(etimehist)));
      drawnow();
    end
    clk = tic;
  end

end
dump = toc(clk);
end

author: kristinbranson

repo: JAABA


text:
function lagrange_interp_1d_test ( )

%*****************************************************************************80
%
%% LAGRANGE_INTERP_1D_TEST tests the LAGRANGE_INTERP_1D library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 August 2012
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../r8lib' )
  addpath ( '../test_interp_1d' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGRANGE_INTERP_1D_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the LAGRANGE_INTERP_1D library.\n' );
  fprintf ( 1, '  The R8LIB library is needed.\n' );
  fprintf ( 1, '  The test needs the TEST_INTERP_1D library.\n' );

  prob_num = p00_prob_num ( );
  for prob = 1 : prob_num
    for nd = [ 4, 8, 16, 32, 64, 256 ]
      lagrange_interp_1d_test02 ( prob, nd );
    end
  end

  prob_num = p00_prob_num ( );
  for prob = 1 : prob_num
    for nd = [ 4, 8, 16, 32, 64, 256 ]
      lagrange_interp_1d_test03 ( prob, nd );
    end
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGRANGE_INTERP_1D_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../r8lib' )
  rmpath ( '../test_interp_1d' )

  return
end


text:
% UPDATE_POSTERIOR updates GPML posterior given a new observation.
%
% This function provides a fast update for a GPML posterior structure
% given a single new observation, avoiding O(N^3) retraining time.
% This implementation assumes exact inference (inf = @infExact) and a
% Gaussian likelihood (lik = @likGauss). It will not work correctly
% for other combinations of likelihood and inference!
%
% Usage:
%
%   new_posterior = update_posterior(hyperparameters, mean_function, ...
%           covariance_function, x, posterior, x_star, y_star)
%
% Inputs:
%
%       hyperparameters: a GPML hyperparameter struct
%         mean_function: a GPML mean function
%   covariance_function: a GPML covariance function
%                     x: the previous observations (N x D)
%             posterior: a GPML posterior struct
%                x_star: the new observation location (1 x D)
%                y_star: the new observation value
%
% Output:
%
%   new_posterior: the updated GPML posterior struct for the GP
%                  conditioned on ([x; x_star], [y; y_star]).
%
% See also INFMETHODS.

% Copyright (c) 2013--2014 Roman Garnett.

function new_posterior = update_posterior(hyperparameters, mean_function, ...
          covariance_function, x, posterior, x_star, y_star, y_sd)

  % check input
  if (numel(y_star) > 1)
    error('gpml_extensions:not_supported', ...
          'update_posterior only supports rank-one updates!');
  end

  % Need to decide how to deal with noise
  if nargin < 8; y_sd = []; end

  noise_variance = exp(2 * hyperparameters.lik);

  k      = feval(covariance_function{:}, hyperparameters.cov, x_star);
  k_star = feval(covariance_function{:}, hyperparameters.cov, x, x_star);

  [m_star, v_star] = gp(hyperparameters, [], mean_function, ...
                        covariance_function, [], x, posterior, x_star);

  % update posterior.L; handle both high- and low-noise
  % parameterizations of posterior
  if (is_chol(posterior.L))
    % high-noise parameterization: posterior.L contains chol(K / sigma^2 + I)

    alpha_update = solve_chol(posterior.L, k_star) * (1 / noise_variance);

    new_L_column = linsolve(posterior.L, k_star, ...
                            struct('UT', true, 'TRANSA', true)) * (1 / noise_variance);
    new_posterior.L = [posterior.L, new_L_column; ...
                       zeros(1, size(posterior.L, 1)), ...
                       sqrt(1 + k / noise_variance - new_L_column' * new_L_column)];
  else
    % low-noise parameterization: posterior.L contains -inv(K + \sigma^2 I)

    alpha_update = -posterior.L * k_star;

    v = -alpha_update / v_star;
    new_posterior.L = [posterior.L + v * alpha_update', -v; -v', -1 / v_star];
  end

  % alpha_update now contains (K + \sigma^2 I) \ k*
  new_posterior.alpha = ...
      [posterior.alpha; 0] + ...
      (m_star - y_star) / v_star * [alpha_update; -1];

  % noise vector is constant; just add one more entry
  new_posterior.sW = posterior.sW([1; (1:end)']);

end

%--------------------------------------------------------------------------
function result = is_chol(L)

  result = (ismatrix(L) && ...
            (size(L, 1) == size(L, 2)) && ...
            isreal(diag(L))  && ...
            all(diag(L) > 0) && ...
            isequal(L, triu(L)));

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


function dataout=modeDRAEC_BSS(nummics, numrefs, datain)
%
% Perform dr and aec together, then bss.
% nummics:              no. of mic channels
% numrefs:              no. of reference channels
% datain:               input data
% dataout:              output data
%

%% perform stft
addpath('stft2');
config;

M=nummics;
R=numrefs;
N=M;

Xtf=cell(M+R, 1);
for m=1:M+R
    Xtf{m}=stft(datain(:, m), stftshift, fftsize, false);
end
[K, T]=size(Xtf{1});

Ytf=cell(M, 1);
for m=1:M
    Ytf{m}=zeros(K, T);
end

%% params go to config

%% space for dr and aec
% current mic data backup
Micbuffer=zeros(K, M*(DR_DELAY+1));
% ref and delayed mic data
draecfsize=R*AEC_FLEN+M*DR_FLEN;
Refmicdelay=zeros(K, draecfsize);

% mic-ref correlation
Cxr=cell(K, 1);
for k=1:K
    Cxr{k}=zeros(M, draecfsize);
end

% reference auto correlation
Crr=cell(K, 1);
for k=1:K
    Crr{k}=zeros(draecfsize, draecfsize);
end

% reverb and echo path
REPath=cell(K, 1);
for k=1:K
    REPath{k}=zeros(M, draecfsize);
end

%% space for bss
% the weighted correlation matrices
C1=cell(K, 1);
C2=cell(K, 1);
for k=1:K
    C1{k}=STABLE_EPS*eye(M, M);
    C2{k}=STABLE_EPS*eye(M, M);
end

% demixing matrices
Demix=cell(K, 1);
for k=1:K
    Demix{k}=eye(N, M);
end

%% perform iteration
for t=1:T
    %% perform dr and aec together
    % direct nearend and early reverberation
    Early=zeros(K, M);

    %
    % shift in new data
    %
    Micbuffer=circshift(Micbuffer, M, 2);
    for m=1:M
        Micbuffer(:, m)=Xtf{m}(:, t);
    end

    % shift in reference data
    Refmicdelay(:, 1:R*AEC_FLEN)=circshift(Refmicdelay(:, 1:R*AEC_FLEN), R, 2);
    for r=1:R
        Refmicdelay(:, r)=Xtf{M+r}(:, t);
    end

    % delayed mic data
    Refmicdelay(:, R*AEC_FLEN+1:end)=circshift(Refmicdelay(:, R*AEC_FLEN+1:end), M, 2);
    Refmicdelay(:, R*AEC_FLEN+1:R*AEC_FLEN+M)=Micbuffer(:, end-M+1:end);

    for k=1:K
        mic=permute(Micbuffer(k, 1:M), [2,1]);
        ref=permute(Refmicdelay(k, :), [2,1]);

        % calculate late reverberation and echo
        late=REPath{k}*ref;

        % direct nearend and early reverberation
        early=mic-late;
        % output data
        Early(k, :)=early;

        %
        % calculate nonlinearity
        %
        xsq=abs(mic).^2;
        ysq=abs(early).^2;

        phi=sum(ysq(ysq<xsq)) + sum(xsq(ysq>=xsq));
        phi=(1-DRAEC_FORGET)*(phi+VAR_BIAS)^((GAMMA-2)/2);

        % update mic ref correlation
        Cxr{k}=DRAEC_FORGET*Cxr{k}+phi*(mic*ref');

        % update ref auto-correlation
        Crr{k}=DRAEC_FORGET*Crr{k}+phi*(ref*ref');

        % update echo and reverb path
        REPath{k}=Cxr{k}/(Crr{k}+DRAEC_DIAGLOAD*eye(draecfsize, draecfsize));
    end

    %% perform bss
    Bssout=zeros(K, M);

    %
    % calculate nonlinearity
    %
    phi1=0;
    phi2=0;

    for k=1:K
        x=Early(k, :).';
        y=Demix{k}*x;
        % output data
        Bssout(k, :)=y.';

        phi1=phi1+abs(y(1))^2;
        phi2=phi2+abs(y(2))^2;
    end

    phi1=(1-BF_FORGET)*(phi1+VAR_BIAS)^((GAMMA-2)/2);
    phi2=(1-BF_FORGET)*(phi2+VAR_BIAS)^((GAMMA-2)/2);

    % update the demixing matrices
    for k=1:K
        %
        % accumulate the weighted correlation
        %
        x=Early(k, :).';
        C1{k}=BF_FORGET*C1{k}+phi1*(x*x');
        C2{k}=BF_FORGET*C2{k}+phi2*(x*x');

        % solve gev problem
        D=heig2(BF_DIAGLOAD, C2{k}, C1{k});
        Demix{k}=D;
    end

    for m=1:M
        Ytf{m}(:, t)=Bssout(:, m);
    end
end

%% perform istft and output signal
dataout=zeros(dataLength(T, stftshift, fftsize ), N);
for n=1:N
    dataout(:, n)=istft(Ytf{n}, stftshift, false);
end

end


text:
% ManhNMF: Manhattan NMF (Guan et al. 2013)
% process_video('NMF', 'ManhNMF', 'dataset/demo.avi', 'output/demo_ManhNMF.avi');
rank = 1;
[W,H] = ManhNMF(M,rank);
L = W' * H;
S = M - L;

author: andrewssobral

repo: lrslibrary


text:
function mask=getMask(l,h,w,batchS,depth,posTempX,posTempY)
mask=abs(posTempX-repmat(l.mu_x,[h,w,1]));
mask=mask+abs(posTempY-repmat(l.mu_y,[h,w,1]));
mask=max(1-mask.*repmat(reshape(l.weights{3},[1,1,depth]),[h,w,1,batchS]),-1);
mask(:,:,l.filter~=1,:)=1;
end

repo: interpretableCNN


text:
function[overlap] = computeBlobOverlapAny(propBlobsA, propBlobsB, imageSize)
% [overlap] = computeBlobOverlapAny(propBlobsA, propBlobsB, imageSize)
%
% Take Sel. Search regions and reconstruct whether blobs overlap at all.
% (result is a matrix).
%
% Copyright by Holger Caesar, 2015

if numel(imageSize) == 3,
    imageSize = imageSize(1:2);
end;
propCountA = numel(propBlobsA);
propCountB = numel(propBlobsB);

% Precompute blob inds
blobsIndsA = cell(propCountA, 1);
blobsIndsB = cell(propCountB, 1);
for i = 1 : propCountA,
    blob = propBlobsA(i);
    blobsIndsA{i} = int32(blobToImageInds(blob, imageSize));
end;
if isequal(propBlobsA, propBlobsB),
    blobsIndsB = blobsIndsA;
else
    for i = 1 : propCountB,
        blob = propBlobsB(i);
        blobsIndsB{i} = int32(blobToImageInds(blob, imageSize));
    end;
end;

overlap = computeBlobOverlapAnyPair(blobsIndsA, blobsIndsB); % Mex-version is 6 times faster for Barcelona dataset

% function[overlap] = computeBlobOverlapAnyPair(blobsIndsA, blobsIndsB) %#ok<DEFNU>
% % [overlap] = computeBlobOverlapAnyPair(blobsIndsA, blobsIndsB)
% %
% % Compute if two blobs overlap at all (boolean).
% % To be replaced by a mex function.
%
% % Initialize
% propCountA = numel(blobsIndsA);
% propCountB = numel(blobsIndsB);
% overlap = zeros(propCountA, propCountB); %Dense is fine
%
% for i = 1 : propCountA,
%     for j = 1 : propCountB,
%         overlap(i, j) = any(builtin('_ismemberhelper', blobsIndsA{i}, blobsIndsB{j}));
%     end;
% end;

author: nightrome

repo: matconvnet-calvin