text:
function [t1,pd] = relaxFitT1_lsq(data,flipAngles,tr,b1Map)
%
% [t1,pd] = relaxFitT1(data,flipAngles,tr,b1Map)
%
% Computes a linear fit of the the T1 estimate for all voxels. The data can
% be passed as either a 4d array of X x Y x Z x nT1Measurements or an array of
% size nVoxels x nT1Measurements. THe b1Map can be either be a scalar (to
% correct for a constant bias across the whole image), or an image the same
% size as one of the t1 measurments (ie. either X x Y x Z or nVoxels x 1).
%
% We allow multiple b1 maps to correct for different sets of measurements.
% If there is more than 1 b1 map, then each row of the flip angle array
% will be corrected by the corresponding b1 map. E.g., flipAngles(1,:) will
% be corrected by b1Map(:,:,:,1), flipAngles(2,:) will be corrected by
% b1Map(:,:,:,2), etc. The flip angle corresponding to each t1 should be
% specified using the linear index. E.g., flipAngle(1) for t1(:,:,:,1),
% flipAngle(2) for t1(:,:,:,2), etc.
%
% Returns:
%   T1: T1 estimate (seconds)
%   PD: a map that includes spin-density (M0), scanner scaling constant
%       (G), and T2*: PD = M0 * G * exp(-TE / T2*).
%
% SEE ALSO:
%
% relaxMtFit.m to fit the f and k maps to the output of this function.
%
% HISTORY:
% 2008.02.26 RFD: wrote it.

if(~exist('b1Map','var')||isempty(b1Map))
  b1Map = 1;
end

theta = flipAngles*pi/180;
szT1 = size(data);
szB1 = size(b1Map);
if(numel(szT1)>2)
   nVox = prod(szT1(1:3));
   nT1 = szT1(4);
   data = reshape(data,nVox,nT1);
   if(numel(b1Map)>1)
       nB1 = szB1(4);
       b1Map = reshape(b1Map,nVox,nB1);
   else
       nB1 = 1;
   end
else
   nVox = size(data,1);
   nT1 = size(data,2);
   nB1 = size(b1Map,2);
end
if(nB1>1)
    if(size(theta,1)~=nB1)
        error('The number of B1 maps must match the number of rows in the flip angle array.');
    end
    b1Inds = repmat([1:nB1]',1,size(theta,2));
else
    b1Inds = repmat([1],1,numel(theta));
end
theta = theta(:)';
b1Inds = b1Inds(:)';

% The code below (esp ndfun) doesn't work when we have Inf or NaN, so:
b1Map(~isfinite(b1Map)) = 1;
data(~isfinite(data)) = 0;

%% NON-LINEAR T1 FIT
%

% Fit a line to the data in each voxel to estimate T1.
% We'll use an eigenvector formulation since we already have a
% vectorized eigenvector decompostion coded up.

% Build a matrix M where M = [x1-x0 y1-y0; x2-x0 y2-y0; ... xn-x0 yn-y0].
% To make it work with ndfun, we need to reshape things a bit.
M = zeros(nT1,2,nVox);
for(ii=1:nT1)
  correctedFlip = theta(ii).*b1Map(:,b1Inds(ii));
  M(ii,1,:) = abs(data(:,ii)./tan(correctedFlip));
  M(ii,2,:) = abs(data(:,ii)./sin(correctedFlip));
end
M0 = mean(M,1);
for(ii=1:size(M,1))
  M(ii,:,:) = M(ii,:,:) - M0;
end
% The best-fitting line is the eigenvector corresponding to
% the largest eigenvalue of eig(M'*M):
[vec,val] = ndfun('eig',ndfun('mult',permute(M,[2 1 3]),M));
% The slope (m) of the line is simply the ratio of y to x and the
% intercept (b) is y-m*x
% Note that the eigenvalues from ndfun are sorted in descending
% order, opposite from Matlab's 'eig'.
m = vec(2,1,:)./vec(1,1,:);
%b = M0(:,2,:)-m.*M0(:,1,:);

% *** CHECK THIS
%m(m<=0) = NaN;
m = abs(squeeze(m));
m(m<0.5) = 0.5;

t1 = (-tr/1000)./log(m);
% Clip to plausible values
t1(isnan(t1)|t1<0.1) = 0.1;
t1(t1>5) = 5;

%pd = b./(1-m);
%pd(pd>5e5) = 5e5;

pd = zeros(nVox,nT1);
% We could use m here, but we've already clipped t1 to plausible values, so recomputing
% the slope here provides a sane range of pd values.
t1tr = exp(-tr/1000./t1);
for(ii=1:nT1)
    correctedFlip = theta(ii).*b1Map(:,b1Inds(ii));
    pd(:,ii) = data(:,ii)./sin(correctedFlip).*((1-cos(correctedFlip).*t1tr)./(1-t1tr));
end
% Could do least squares for better PDmap?
pd = mean(pd,2);

if(numel(szT1)>ndims(data))
   t1 = reshape(t1,szT1(1:3));
   pd = reshape(pd,szT1(1:3));
end

return

author: vistalab

repo: vistasoft


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Eugenio Alcala Baselga
% Date: 02/06/2018
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ A, B ] = A_OL_KINEMATIC_3VARS_DISCRETE( vd, omega, theta_Err, Ts )

    A3by3 = [1        omega*Ts   0;
            -omega*Ts   1     vd*(sin(theta_Err)/theta_Err)*Ts;
             0          0     1];

    B3by2 = [-Ts 0;
              0 0;
              0 -Ts];
    A = A3by3;
    B = B3by2;
end


text:
function [flag] = VBA_isBinary (X)
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [flag] = VBA_isBinary (X)
% check if X only contains 0 or 1 values
%
% IN:
%   - X: matrix, cell array of matrices, or structure to be checked
% OUT:
%   - flag: - true if all elements of X are binary (0 or 1)
%           - false otherwise
%           - NaN if any element or field is not numeric
%
% /////////////////////////////////////////////////////////////////////////

switch class (X)
    case 'logical'
        flag = true;

    case 'double'
        flag = all (ismember (X(:), [0, 1]));

    case 'cell'
        flag = all (cellfun (@VBA_isBinary, X));

    case 'struct'
        flag = all (structfun (@VBA_isBinary, X));

    otherwise
        flag = NaN;
end

author: MBB-team

repo: VBA-toolbox


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function [lower, upper] = ...
    UpperBound1(S, g, df, B, beta, lower, Nr, NSim, NSSim, ...
    getpaths, payoff)
% method from Broadie see Chapter 8

v = g(:, end);          % option value at maturity
c = zeros(NSim,Nr-1);   % continuation value

stoppingtime = Nr * ones(NSim,1);

 for i=Nr-1:-1:1
         index = find(g(:,i)>0);
         s = S(index, i+1);
         v = v * df(i+1);

         Acell = B(s);
         A = cell2mat(Acell{:,:});
         c(index,i) = A*beta(:,i);                % cont value itm paths

         earlyexercise = g(index,i) >= c(index,i);% exercise
         stoppingtime(index(earlyexercise)) = i;
         v(index(earlyexercise)) = g(index(earlyexercise),i);
 end

% Computing the martingale pi

L = zeros(NSim,1);
expectation = zeros(NSim,Nr);
expectation(:,1) = lower * ones(NSim,1);

pi = zeros(NSim,Nr+1);                  % martingale
pi(:,1) = lower * ones(NSim,1);         % set pi(0)

% perform the subsimulation for each value of the path set
for i=1:1:Nr-1
    for j=1:NSim
        expectation(j,i+1) = subsimulation(S(j,i+1), df, B, NSSim, Nr-i,...
            beta(:,i:end), getpaths, payoff) * prod(df(1:i));
    end
end

for i=1:1:Nr
    i_exercise = stoppingtime == i;     % check if early exercised

    if i < Nr
        L(i_exercise) = g(i_exercise, i) * prod(df(1:i));
        L(~i_exercise) = expectation(~i_exercise,i+1);
    else
        L(i_exercise) = g(i_exercise, i) * prod(df(1:i));
    end

    pi(:,i+1) = pi(:,i) + L - expectation(:,i);
end

maximum = zeros(NSim,1);

for j=1:1:NSim
    maximum(j) = L(j) + max(g(j,:) - pi(j,2:end)); % compute the max
end
upper = mean(maximum);

end

function y = subsimulation(S0, df, B, NSim, Nr, beta, gp,payoff)
    S2 = gp(S0,NSim,Nr); S2 = S2(:,2:end);   % paths
    g2 = payoff(S2);                         % payoff
    iVec = 1:NSim;

    exercise = Nr * ones(NSim,1);            % exercise per path
    % determine exercise strategy for path set S2
    for i=1:1:Nr-1
        i_nexercised = exercise == Nr;
        I_nexercise = iVec(i_nexercised);

        s = S2(i_nexercised,i);
        Acell = B(s);
        A = cell2mat(Acell{:,:});
        c = A * beta(:,i);

        i_exercise = g2(i_nexercised,i) >= c & g2(i_nexercised,i) > 0;
        exercise(I_nexercise(i_exercise)) = i;
    end
    summe=0;
    for j=1:1:NSim
        summe = summe + g2(j,exercise(j)) * prod(df(1:exercise(j)));
    end
    y = summe / NSim;                         % MC value from subsimulation
end

author: Sable

repo: mcbench-benchmarks


text:
function [divTr, divTe]= sample_chronKfold(label, folds)
%SAMPLE_CHRONKFOLD - Sampling function: chronological consequtive folds
%
%Synopsis:
%  [PARTR, PARTE]= sample_chronKfold(LABEL, FOLDS)
%
% IN  LABEL   - class labels, but here only the length of the 2nd dim is
%               used (to determine the number of samples)
%     FOLDS   - DOUBLE nFolds: number of folds into which the samples are
%               divided. Or FOLDS can be [nShifts nFolds] in which case
%               all partitions will also be generated in shifted versions.
%
% OUT PARTR   - Partitions of the training set
%               PARTR{n}: cell array holding the training sets
%               folds for shuffle #n, more specificially
%               PARTR{n}{m} holds the indices of the training set of
%               the m-th fold of shuffle #n.
%     PARTE   - analogue to PARTR, for the test sets

% Benjamin Blankertz

misc_checkType(label, 'DOUBLE[- -]');
misc_checkType(folds, 'DOUBLE|DOUBLE[2]');

nSamples= size(label,2);

if length(folds)==1
  folds= [1 folds];
end

divTr= {cell(1,folds(2))};
divTe= {cell(1,folds(2))};
for nn= 1:folds(1)
  div= round(linspace(0+(nn-1), nSamples-(folds(1)-nn), folds(2)+1));
  for kk= 1:folds(2)
    divTe{nn}{kk}= div(kk)+1:div(kk+1);
    divTr{nn}{kk}= setdiff(nn:nSamples-(folds(1)-nn), divTe{nn}{kk});
    % check that all classes are inhabited
    if ~all(sum(label(:,divTr{nn}{kk}))),
      error('empty classes in training set');
    end
  end
end

author: bbci


text:
function r = sobol_dataset ( m, n, skip )

%*****************************************************************************80
%
%% SOBOL_DATASET generates a Sobol dataset and writes it to a file.
%
%  Usage:
%
%    r = sobol_dataset ( m, n, skip )
%
%    where
%
%    * M, the spatial dimension,
%    * N, the number of points to generate,
%    * SKIP, the number of initial values to skip.
%    * R is the M by N array created.
%
%    creates an M by N Sobol dataset and writes it to the
%    file "sobol_M_N.txt".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 December 2009
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SOBOL_DATASET\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Generate a Sobol dataset.\n' );
%
%  Get the spatial dimension.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    m = input ( '  Enter the spatial dimension M: ' );
  else
    m = str2num ( m );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension M = %d\n', m );
%
%  Get the number of points.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    n = input ( '  Enter the number of points N: ' );
  else
    n = str2num ( n );
  end

  fprintf ( 1, '  Number of points N = %d\n', n );
%
%  Get SKIP.
%
  if ( nargin < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  SKIP is the number of points to skip.\n' );
    fprintf ( 1, '  0 is a common choice.\n' );
    fprintf ( 1, '  Another is the least power of 2 greater than or equal to N.\n' );
    skip = input ( '  Enter SKIP: ' );
  else
    skip = str2num ( skip );
  end

  fprintf ( 1, '  SKIP = %d\n', skip );
%
%  Compute the data.
%
  r = i4_sobol_generate ( m, n, skip );
%
%  Write it to a file.
%
  output_filename = ...
    strcat ( 'sobol_', num2str ( m ), '_', num2str ( n ), '.txt' );

  r8mat_write ( output_filename, m, n, r );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data was written to the file "%s".\n', ...
    output_filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SOBOL_DATASET:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function bit = i4_bit_hi1 ( n )

%*****************************************************************************80
%
%% I4_BIT_HI1 returns the position of the high 1 bit base 2 in an integer.
%
%  Example:
%
%       N    Binary     BIT
%    ----    --------  ----
%       0           0     0
%       1           1     1
%       2          10     2
%       3          11     2
%       4         100     3
%       5         101     3
%       6         110     3
%       7         111     3
%       8        1000     4
%       9        1001     4
%      10        1010     4
%      11        1011     4
%      12        1100     4
%      13        1101     4
%      14        1110     4
%      15        1111     4
%      16       10000     5
%      17       10001     5
%    1023  1111111111    10
%    1024 10000000000    11
%    1025 10000000001    11
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the integer to be measured.
%    N should be nonnegative.  If N is nonpositive, the value will always be 0.
%
%    Output, integer BIT, the number of bits base 2.
%
  i = floor ( n );
  bit = 0;

  while ( 1 )

    if ( i <= 0 )
      break;
    end

    bit = bit + 1;
    i = floor ( i / 2 );

  end

  return
end

function bit = i4_bit_lo0 ( n )

%*****************************************************************************80
%
%% I4_BIT_LO0 returns the position of the low 0 bit base 2 in an integer.
%
%  Example:
%
%       N    Binary     BIT
%    ----    --------  ----
%       0           0     1
%       1           1     2
%       2          10     1
%       3          11     3
%       4         100     1
%       5         101     2
%       6         110     1
%       7         111     4
%       8        1000     1
%       9        1001     2
%      10        1010     1
%      11        1011     3
%      12        1100     1
%      13        1101     2
%      14        1110     1
%      15        1111     5
%      16       10000     1
%      17       10001     2
%    1023  1111111111     1
%    1024 10000000000     1
%    1025 10000000001     1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the integer to be measured.
%    N should be nonnegative.
%
%    Output, integer BIT, the position of the low 1 bit.
%
  bit = 0;
  i = floor ( n );

  while ( 1 )

    bit = bit + 1;
    i2 = floor ( i / 2 );

    if ( i == 2 * i2 )
      break;
    end

    i = i2;

  end

  return
end
function [ quasi, seed ] = i4_sobol ( dim_num, seed )

%*****************************************************************************80
%
%% I4_SOBOL generates a new quasirandom Sobol vector with each call.
%
%  Discussion:
%
%    The routine adapts the ideas of Antonov and Saleev.
%
%    Thanks to Francis Dalaudier for pointing out that the range of allowed
%    values of DIM_NUM should start at 1, not 2!  17 February 2009.
%
%    This function was modified to use PERSISTENT variables rather than
%    GLOBAL variables, 13 December 2009.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Bennett Fox.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Antonov, Saleev,
%    USSR Computational Mathematics and Mathematical Physics,
%    Volume 19, 1980, pages 252 - 256.
%
%    Paul Bratley, Bennett Fox,
%    Algorithm 659:
%    Implementing Sobol's Quasirandom Sequence Generator,
%    ACM Transactions on Mathematical Software,
%    Volume 14, Number 1, pages 88-100, 1988.
%
%    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.
%
%    Ilya Sobol,
%    USSR Computational Mathematics and Mathematical Physics,
%    Volume 16, pages 236-242, 1977.
%
%    Ilya Sobol, Levitan,
%    The Production of Points Uniformly Distributed in a Multidimensional
%    Cube (in Russian),
%    Preprint IPM Akad. Nauk SSSR,
%    Number 40, Moscow 1976.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the number of spatial dimensions.
%    DIM_NUM must satisfy 1 <= DIM_NUM <= 40.
%
%    Input/output, integer SEED, the "seed" for the sequence.
%    This is essentially the index in the sequence of the quasirandom
%    value to be generated.  On output, SEED has been set to the
%    appropriate next value, usually simply SEED+1.
%    If SEED is less than 0 on input, it is treated as though it were 0.
%    An input value of 0 requests the first (0-th) element of the sequence.
%
%    Output, real QUASI(DIM_NUM), the next quasirandom vector.
%
  persistent atmost;
  persistent dim_max;
  persistent dim_num_save;
  persistent initialized;
  persistent lastq;
  persistent log_max;
  persistent maxcol;
  persistent poly;
  persistent recipd;
  persistent seed_save;
  persistent v;

  if ( isempty ( initialized ) )
    initialized = 0;
    dim_num_save = -1;
  end

  if ( ~initialized | dim_num ~= dim_num_save )

    initialized = 1;

    dim_max = 40;
    dim_num_save = -1;
    log_max = 30;
    seed_save = -1;
%
%  Initialize (part of) V.
%
    v(1:dim_max,1:log_max) = zeros(dim_max,log_max);

    v(1:40,1) = [ ...
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ]';

    v(3:40,2) = [ ...
            1, 3, 1, 3, 1, 3, 3, 1, ...
      3, 1, 3, 1, 3, 1, 1, 3, 1, 3, ...
      1, 3, 1, 3, 3, 1, 3, 1, 3, 1, ...
      3, 1, 1, 3, 1, 3, 1, 3, 1, 3 ]';

    v(4:40,3) = [ ...
               7, 5, 1, 3, 3, 7, 5, ...
      5, 7, 7, 1, 3, 3, 7, 5, 1, 1, ...
      5, 3, 3, 1, 7, 5, 1, 3, 3, 7, ...
      5, 1, 1, 5, 7, 7, 5, 1, 3, 3 ]';

    v(6:40,4) = [ ...
                     1, 7, 9,13,11, ...
      1, 3, 7, 9, 5,13,13,11, 3,15, ...
      5, 3,15, 7, 9,13, 9, 1,11, 7, ...
      5,15, 1,15,11, 5, 3, 1, 7, 9 ]';

    v(8:40,5) = [ ...
                           9, 3,27, ...
     15,29,21,23,19,11,25, 7,13,17, ...
      1,25,29, 3,31,11, 5,23,27,19, ...
     21, 5, 1,17,13, 7,15, 9,31, 9 ]';

    v(14:40,6) = [ ...
              37,33, 7, 5,11,39,63, ...
     27,17,15,23,29, 3,21,13,31,25, ...
      9,49,33,19,29,11,19,27,15,25 ]';

    v(20:40,7) = [ ...
                                         13, ...
     33,115, 41, 79, 17, 29,119, 75, 73,105, ...
      7, 59, 65, 21,  3,113, 61, 89, 45,107 ]';

    v(38:40,8) = [ ...
                                7, 23, 39 ]';
%
%  Set POLY.
%
    poly(1:40)= [ ...
        1,   3,   7,  11,  13,  19,  25,  37,  59,  47, ...
       61,  55,  41,  67,  97,  91, 109, 103, 115, 131, ...
      193, 137, 145, 143, 241, 157, 185, 167, 229, 171, ...
      213, 191, 253, 203, 211, 239, 247, 285, 369, 299 ];

    atmost = 2^log_max - 1;
%
%  Find the number of bits in ATMOST.
%
    maxcol = i4_bit_hi1 ( atmost );
%
%  Initialize row 1 of V.
%
    v(1,1:maxcol) = 1;

  end
%
%  Things to do only if the dimension changed.
%
  if ( dim_num ~= dim_num_save )
%
%  Check parameters.
%
    if ( dim_num < 1 | dim_max < dim_num )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'I4_SOBOL - Fatal error!\n' );
      fprintf ( 1, '  The spatial dimension DIM_NUM should satisfy:\n' );
      fprintf ( 1, '    1 <= DIM_NUM <= %d\n', dim_max );
      fprintf ( 1, '  But this input value is DIM_NUM = %d\n', dim_num );
      return
    end

    dim_num_save = dim_num;
%
%  Initialize the remaining rows of V.
%
    for i = 2 : dim_num
%
%  The bits of the integer POLY(I) gives the form of polynomial I.
%
%  Find the degree of polynomial I from binary encoding.
%
      j = poly(i);
      m = 0;

      while ( 1 )

        j = floor ( j / 2 );

        if ( j <= 0 )
          break;
        end

        m = m + 1;

      end
%
%  Expand this bit pattern to separate components of the logical array INCLUD.
%
      j = poly(i);
      for k = m : -1 : 1
        j2 = floor ( j / 2 );
        includ(k) = ( j ~= 2 * j2 );
        j = j2;
      end
%
%  Calculate the remaining elements of row I as explained
%  in Bratley and Fox, section 2.
%
      for j = m + 1 : maxcol
        newv = v(i,j-m);
        l = 1;
        for k = 1 : m
          l = 2 * l;
          if ( includ(k) )
            newv = bitxor ( newv, l * v(i,j-k) );
          end
        end
        v(i,j) = newv;
      end
    end
%
%  Multiply columns of V by appropriate power of 2.
%
    l = 1;
    for j = maxcol-1 : -1 : 1
      l = 2 * l;
      v(1:dim_num,j) = v(1:dim_num,j) * l;
    end
%
%  RECIPD is 1/(common denominator of the elements in V).
%
    recipd = 1.0 / ( 2 * l );

    lastq(1:dim_num) = 0;

  end

  seed = floor ( seed );

  if ( seed < 0 )
    seed = 0;
  end

  if ( seed == 0 )

    l = 1;
    lastq(1:dim_num) = 0;

  elseif ( seed == seed_save + 1 )
%
%  Find the position of the right-hand zero in SEED.
%
    l = i4_bit_lo0 ( seed );

  elseif ( seed <= seed_save )

    seed_save = 0;
    l = 1;
    lastq(1:dim_num) = 0;

    for seed_temp = seed_save : seed - 1
      l = i4_bit_lo0 ( seed_temp );
      for i = 1 : dim_num
        lastq(i) = bitxor ( lastq(i), v(i,l) );
      end
    end

    l = i4_bit_lo0 ( seed );

  elseif ( seed_save + 1 < seed )

    for seed_temp = seed_save + 1 : seed - 1
      l = i4_bit_lo0 ( seed_temp );
      for i = 1 : dim_num
        lastq(i) = bitxor ( lastq(i), v(i,l) );
      end
    end

    l = i4_bit_lo0 ( seed );

  end
%
%  Check that the user is not calling too many times!
%
  if ( maxcol < l )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_SOBOL - Fatal error!\n' );
    fprintf ( 1, '  Too many calls!\n' );
    fprintf ( 1, '  MAXCOL = %d\n', maxcol );
    fprintf ( 1, '  L =      %d\n', l );
    return
  end
%
%  Calculate the new components of QUASI.
%
  for i = 1 : dim_num
    quasi(i) = lastq(i) * recipd;
    lastq(i) = bitxor ( lastq(i), v(i,l) );
  end

  seed_save = seed;
  seed = seed + 1;

  return
end
function r = i4_sobol_generate ( m, n, skip )

%*****************************************************************************80
%
%% I4_SOBOL_GENERATE generates a Sobol dataset.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 December 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points to generate.
%
%    Input, integer SKIP, the number of initial points to skip.
%
%    Output, real R(M,N), the points.
%
  seed = skip;

  for j = 1 : n
    [ r(1:m,j), seed ]  = i4_sobol ( m, seed );
  end

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
    return;
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function bbci= bbci_calibrate_reset(bbci, varargin)
%BBCI_CALIBRATE_RESET - Revert calibration settings to original values
%
%Synopsis:
%  BBCI= bbci_calibrate_reset(BBCI)
%  BBCI= bbci_calibrate_reset(BBCI, 'param1', 'param2', ...)
%
%Argument:
%  BBCI: Structure holding specific settings for calibration in the
%     field 'calibrate'.
%
%Output:
%  BBCI - Updated BBCI structure in which all subfields of BBCI.settings
%     which are specified in BBCI.calibrate.auto_reset are reverted to
%     the values in BBCI.calibrate.default_settings.
%     For BBCI.calibrate.auto_reset=1 all subfields are reverted.

% 11-2011 - Benjamin Blankertz

BC= bbci.calibrate;
if ~isfield(BC, 'default_settings'),
  error('no default settings stored in BBCI structure');
end
BCS= bbci.calibrate.default_settings;

if nargin==1,
  fld_list= fieldnames(BCS);
elseif length(varargin)==1 && iscell(varargin{1}),
  fld_list= varargin{1};
else
  fld_list= varargin;
end

for ii= 1:length(fld_list),
  fld= fld_list{ii};
  bbci.calibrate.settings.(fld)= bbci.calibrate.default_settings.(fld);
end


text:
function [centroid, area] = polygonCentroid(varargin)
%POLYGONCENTROID Compute the centroid (center of mass) of a polygon
%
%   CENTROID = polygonCentroid(POLY)
%   CENTROID = polygonCentroid(PTX, PTY)
%   Computes center of mass of a polygon defined by POLY. POLY is a N-by-2
%   array of double containing coordinates of vertices.
%
%   [CENTROID AREA] = polygonCentroid(POLY)
%   Also returns the (signed) area of the polygon.
%
%   Example
%     % Draws the centroid of a paper hen
%     x = [0 10 20  0 -10 -20 -10 -10  0];
%     y = [0  0 10 10  20  10  10  0 -10];
%     poly = [x' y'];
%     centro = polygonCentroid(poly);
%     drawPolygon(poly);
%     hold on; axis equal;
%     drawPoint(centro, 'bo');
%
%   References
%   algo adapted from P. Bourke web page
%
%   See also:
%   polygons2d, polygonArea, polygonSecondAreaMoments, drawPolygon
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 05/05/2004.
%

% Algorithme P. Bourke, vectorized version

% HISTORY
% 2012.02.24 vectorize code

% parse input arguments
if nargin == 1
    var = varargin{1};
    px = var(:,1);
    py = var(:,2);
elseif nargin == 2
    px = varargin{1};
    py = varargin{2};
end

% vertex indices
N = length(px);
iNext = [2:N 1];

% compute cross products
common = px .* py(iNext) - px(iNext) .* py;
sx = sum((px + px(iNext)) .* common);
sy = sum((py + py(iNext)) .* common);

% area and centroid
area = sum(common) / 2;
centroid = [sx sy] / 6 / area;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function F = vertcat( varargin )
%VERTCAT   Vertical concatenation of BALLFUN objects.
%
% K = VERTCAT(F, G, H) is the vertical concatenation of BALLFUN objects F,
% G, and H. The function K is a BALLFUNV object.
%
% [F ; G ; H] is equivalent to VERTCAT(F, G, H).
%
% VERTCAT(F, G) returns an error. BALLFUNV objects have three components.
%
% VERTCAT(F) returns the BALLFUN F.
%
% See also BALLFUNV.

% Copyright 2018 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

if ( nargin == 1 )
    % VERTCAT of one argument just returns the same thing back to the user:
    F = varargin{1};
elseif ( nargin == 3 )
    if all( cellfun(@(F) isa(F,'ballfun'), varargin) )
        F = ballfunv( varargin{:} );
    else
        error('BALLFUN:vertcat:tooManyComponents', ...
            'Only BALLFUN objects are valid to concatenate.');
    end
else
    error('BALLFUN:vertcat:tooManyInputs', ...
        'Can only vertically concatenate three BALLFUN objects.');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
BTB_memo= BTB;
BTB.RawDir= fullfile(BTB.DataDir, 'demoRaw');
BTB.MatDir= fullfile(BTB.DataDir, 'demoMat');

% add more to the list if you want to do it in a row
subdir_list= {'VPkg_08_08_07'};
% you could have more files also
basename_list= {'calibration_motorimagery', ...
                'feedback_motorimagery'};

% definition of classes based on markers
stimDef= {1, 2, 3;
          'left','right', 'foot'};

% load raw files (with filtering), define classes and montage,
% and save in matlab format
for k= 1:length(subdir_list);
 for ib= 1:length(basename_list),
  subdir= subdir_list{k};
  sbj= subdir(1:find(subdir=='_',1,'first')-1);
  file= fullfile(subdir, [basename_list{ib} sbj]);
  fprintf('converting %s\n', file)

  [cnt, mrk_orig, hdr] = file_readBV(file);

  % create mrk and mnt
  mrk= mrk_defineClasses(mrk_orig, stimDef);
  mrk.orig= mrk_orig;
  mnt= mnt_setElectrodePositions(cnt.clab);
  mnt= mnt_setGrid(mnt, 'M+EOG+EMG');

  % save in matlab format
  file_saveMatlab(file, cnt, mrk, mnt, 'Vars','hdr');
 end
end

BTB= BTB_memo;


text:
function [gp, varargout] = gp_optim(gp, x, y, varargin)
%GP_OPTIM  Optimize paramaters of a Gaussian process
%
%  Description
%    GP = GP_OPTIM(GP, X, Y, OPTIONS) optimises the parameters of a
%    GP structure given matrix X of training inputs and vector
%    Y of training targets.
%
%    [GP, OUTPUT1, OUTPUT2, ...] = GP_OPTIM(GP, X, Y, OPTIONS)
%    optionally returns outputs of the optimization function.
%
%    OPTIONS is optional parameter-value pair
%      z      - optional observed quantity in triplet (x_i,y_i,z_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, expected
%               value for ith case.
%      optimf - function handle for an optimization function, which is
%               assumed to have similar input and output arguments
%               as usual fmin*-functions. Default is @fminscg.
%      opt    - options structure for the minimization function.
%               Use optimset to set these options. By default options
%               'GradObj' is 'on', 'LargeScale' is 'off'.
%      loss   - 'e' to minimize the marginal posterior energy (default) or
%               'loo' to minimize the negative leave-one-out lpd
%               'kfcv' to minimize the negative k-fold-cv lpd
%               'waic' to minimize the WAIC loss
%               only 'e' and 'loo' with Gaussian likelihood have gradients
%      k      - number of folds in kfcv
%
%  See also
%    GP_SET, GP_E, GP_G, GP_EG, FMINSCG, FMINLBFGS, OPTIMSET, DEMO_REGRESSION*
%
% Copyright (c) 2010-2012 Aki Vehtari

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

ip=inputParser;
ip.FunctionName = 'GP_OPTIM';
ip.addRequired('gp',@(x) isstruct(x) || isempty(x));
ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('optimf', @fminscg, @(x) isa(x,'function_handle'))
ip.addParamValue('opt', [], @isstruct)
ip.addParamValue('loss', 'e', @(x) ismember(lower(x),{'e', 'loo', 'kfcv', 'waic' 'waic' 'waicv' 'waicg'}))
ip.addParamValue('k', 10, @(x) isreal(x) && isscalar(x) && isfinite(x) && x>0)
ip.parse(gp, x, y, varargin{:});
if isempty(gp)
  gp=gp_set();
end
if isempty(gp_pak(gp))
  % nothing to optimize
  return
end
z=ip.Results.z;
optimf=ip.Results.optimf;
opt=ip.Results.opt;
loss=ip.Results.loss;
k=ip.Results.k;

switch lower(loss)
  case 'e'
    fh_eg=@(ww) gp_eg(ww, gp, x, y, 'z', z);
    optdefault=struct('GradObj','on','LargeScale','off');
  case 'loo'
    fh_eg=@(ww) gp_looeg(ww, gp, x, y, 'z', z);
    if isfield(gp.lik.fh,'trcov') || isequal(gp.latent_method, 'EP')
      optdefault=struct('GradObj','on','LargeScale','off');
    else
      % Laplace-LOO does not have yet gradients
      optdefault=struct('Algorithm','interior-point');
      if ismember('optimf',ip.UsingDefaults)
        optimf=@fmincon;
      end
    end
  case 'kfcv'
    % kfcv does not have yet gradients
    fh_eg=@(ww) gp_kfcve(ww, gp, x, y, 'z', z, 'k', k);
    optdefault=struct('Algorithm','interior-point');
    if ismember('optimf',ip.UsingDefaults)
      optimf=@fmincon;
    end
  case {'waic' 'waicv'}
    % waic does not have yet gradients
    fh_eg=@(ww) -gp_waic(gp_unpak(gp,ww), x, y, 'z', z);
    optdefault=struct('Algorithm','interior-point');
    if ismember('optimf',ip.UsingDefaults)
      optimf=@fmincon;
    end
  case 'waicg'
    % waic does not have yet gradients
    fh_eg=@(ww) -gp_waic(gp_unpak(gp,ww), x, y, 'z', z, 'method', 'G');
    optdefault=struct('Algorithm','interior-point');
    if ismember('optimf',ip.UsingDefaults)
      optimf=@fmincon;
    end
end
opt=setOpt(optdefault,opt);
w=gp_pak(gp);
if isequal(lower(loss),'e') || (isequal(lower(loss),'loo')) && (isfield(gp.lik.fh,'trcov') || isequal(gp.latent_method, 'EP'))
  switch nargout
    case 6
      [w,fval,exitflag,output,grad,hessian] = optimf(fh_eg, w, opt);
      varargout={fval,exitflag,output,grad,hessian};
    case 5
      [w,fval,exitflag,output,grad] = optimf(fh_eg, w, opt);
      varargout={fval,exitflag,output,grad};
    case 4
      [w,fval,exitflag,output] = optimf(fh_eg, w, opt);
      varargout={fval,exitflag,output};
    case 3
      [w,fval,exitflag] = optimf(fh_eg, w, opt);
      varargout={fval,exitflag};
    case 2
      [w,fval] = optimf(fh_eg, w, opt);
      varargout={fval};
    case 1
      w = optimf(fh_eg, w, opt);
      varargout={};
  end
else
  lb=repmat(-8,size(w));
  ub=repmat(10,size(w));
  switch nargout
    case 6
      [w,fval,exitflag,output,grad,hessian] = optimf(fh_eg, w, [], [], [], [], lb, ub, [], opt);
      varargout={fval,exitflag,output,grad,hessian};
    case 5
      [w,fval,exitflag,output,grad] = optimf(fh_eg, w, [], [], [], [], lb, ub, [], opt);
      varargout={fval,exitflag,output,grad};
    case 4
      [w,fval,exitflag,output] = optimf(fh_eg, w, [], [], [], [], lb, ub, [], opt);
      varargout={fval,exitflag,output};
    case 3
      [w,fval,exitflag] = optimf(fh_eg, w, [], [], [], [], lb, ub, [], opt);
      varargout={fval,exitflag};
    case 2
      [w,fval] = optimf(fh_eg, w, [], [], [], [], lb, ub, [], opt);
      varargout={fval};
    case 1
      w = optimf(fh_eg, w, [], [], [], [], lb, ub, [], opt);
      varargout={};
  end
end
gp=gp_unpak(gp,w);
end

function opt=setOpt(optdefault, opt)
  % Set default options
  opttmp=optimset(optdefault,opt);

  % Set some additional options for @fminscg
  if isfield(opt,'lambda')
    opttmp.lambda=opt.lambda;
  end
  if isfield(opt,'lambdalim')
    opttmp.lambdalim=opt.lambdalim;
  end
  opt=opttmp;
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function plot_3d_hdc_sheet_gtvo(visualDataFile, expMapHistoryFile, odoMapHistoryFile, groundTruthFile,vtHistoryFile, emHistoryFile, varargin)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

    %% Initial some variables
    % Transform the angle from degree to radian
    global DEGREE_TO_RADIAN;
    DEGREE_TO_RADIAN = pi / 180;

    % Transform the angle from radian to degree
    global RADIAN_TO_DEGREE;
    RADIAN_TO_DEGREE = 180 / pi;

    % define the size of block in reading video
    global BLOCK_READ;

    % define the render rate for dyawing vo
    global RENDER_RATE;

    global GT_ODO_X_SCALING;
    global GT_ODO_Y_SCALING;
    global GT_ODO_Z_SCALING;
    global GT_EXP_X_SCALING;
    global GT_EXP_Y_SCALING;
    global GT_EXP_Z_SCALING;

    global ODO_MAP_X_SCALING;
    global ODO_MAP_Y_SCALING;
    global ODO_MAP_Z_SCALING;
    global EXP_MAP_X_SCALING;
    global EXP_MAP_Y_SCALING;
    global EXP_MAP_Z_SCALING;

    % Process the parameters
    for i=1:(nargin - 6)
        if ischar(varargin{i})
            switch varargin{i}
                case 'BLOCK_READ', BLOCK_READ = varargin{i+1};
                case 'RENDER_RATE', RENDER_RATE = varargin{i+1};
                case 'GT_ODO_X_SCALING', GT_ODO_X_SCALING = varargin{i+1};
                case 'GT_ODO_Y_SCALING', GT_ODO_Y_SCALING = varargin{i+1};
                case 'GT_ODO_Z_SCALING', GT_ODO_Z_SCALING = varargin{i+1};
                case 'GT_EXP_X_SCALING', GT_EXP_X_SCALING = varargin{i+1};
                case 'GT_EXP_Y_SCALING', GT_EXP_Y_SCALING = varargin{i+1};
                case 'GT_EXP_Z_SCALING', GT_EXP_Z_SCALING = varargin{i+1};
                case 'ODO_MAP_X_SCALING', ODO_MAP_X_SCALING = varargin{i+1};
                case 'ODO_MAP_Y_SCALING', ODO_MAP_Y_SCALING = varargin{i+1};
                case 'ODO_MAP_Z_SCALING', ODO_MAP_Z_SCALING = varargin{i+1};
                case 'EXP_MAP_X_SCALING', EXP_MAP_X_SCALING = varargin{i+1};
                case 'EXP_MAP_Y_SCALING', EXP_MAP_Y_SCALING = varargin{i+1};
                case 'EXP_MAP_Z_SCALING', EXP_MAP_Z_SCALING = varargin{i+1};
            end
        end
    end

    %% vt
    global PREV_VT_ID;
    global SUB_VT_IMG;
    global VT_HISTORY;
    global IMG_TYPE;

    %% vo
    % for drawing images for estimating pitch and yaw
    global SUB_HORI_TRANS_IMG;
    global SUB_ROT_IMG;
    global SUB_VERT_TRANS_IMG;
    global SUB_PITCH_IMG;

    % for limit the velocity of gtvo
    global MAX_HORI_TRANS_V_THRESHOLD;
    global MAX_ROT_V_THRESHOLD;
    global MAX_VERT_TRANS_V_THRESHOLD;
    global MAX_PITCH_V_THRESHOLD;

    %% yaw_pitch_hdc
    global YAW_PITCH_HDC;

    % The dimension of yaw in yaw_pitch_hdc network
    global YAW_PITCH_HDC_Y_DIM;

    % The dimension of pitch in yaw_pitch_hdc network
    global YAW_PITCH_HDC_P_DIM;

    % The yaw theta size of each unit in radian, 2*pi/ YAW_PITCH_HDC_Y_DIM
    % radian e.g. 2*pi/360 = 0.0175
    global YAW_PITCH_HDC_Y_TH_SIZE;
    global YAW_PITCH_HDC_P_TH_SIZE;
    global MAX_ACTIVE_YAW_PITCH_HIS_PATH;

    %% 3D gridcells
    global GRIDCELLS;

    global GC_X_DIM;
    global GC_Y_DIM;
    global GC_Z_DIM;

    % set the initial position in the grid cell network
    [gcX, gcY, gcZ] = get_gc_initial_pos();

    global MAX_ACTIVE_XYZ_PATH;

    %% 3D EM
    global EXPERIENCES;

    % Experience history
    global EXP_HISTORY;

    % The namber of total experiences
    global NUM_EXPS;

    %% read imgs as VT input

    % get the visual data info
    [subFoldersPathSet, numSubFolders] = get_images_data_info(visualDataFile);

    % initial some variables for drawing pitch and yaw
    odoPitchTheta = [0 0 0];
    Theta = [0 0 0 0 0 0];
    Rho = [0 0 0 0 0 0];
    startpoint =[0 0];
    endpoint = [0.8 0.8];

    odoYawTheta = [0 0 0];

    hdcPitchTheta = [0 0 0];

    hdcYawTheta = [0 0 0];

    hdcPitchTheta(1,3)= 0;
    hdcYawTheta(1,3)= 0;

    zAxisHDC = 1:2:36;

    % for drawing the activity packet of HDC
    zThirdDim = 36;
    tempYawPitchHDC = zeros(YAW_PITCH_HDC_Y_DIM, YAW_PITCH_HDC_P_DIM, zThirdDim);

    % for drawing the activity packet in toriodal ring
    R = 12;
    r = 6;
    xTorus = nan(36,36);
    yTorus = nan(36,36);
    zTorus = nan(36,36);
    for j = 1 : 36
        for i = 1 : 36
            phi = (i - 1) * 2 * pi / 35;
            theta = (j - 1) * 2 * pi / 35;
            xTorus(i,j)= cos(phi) * (R + cos(theta) * r);
            yTorus(i,j)= sin(phi) * (R + cos(theta) * r);
            zTorus(i,j)= sin(theta) * r;
        end
    end

    %% for drawing contour3 in 3D ring
    % define color map based on activity in 3D ring
%     coMap = zeros(YAW_PITCH_HDC_Y_DIM * YAW_PITCH_HDC_P_DIM, 3);

    %% processing two ring model for hd integration
    % x,y,z, yaw, pitch
    expTrajectory(1,1) = 0;
    expTrajectory(1,2) = 0;
    expTrajectory(1,3) = 0;
    expTrajectory(1,4) = 0;
    expTrajectory(1,5) = 0;

    odoMapTrajectory(1,1) = 0;
    odoMapTrajectory(1,2) = 0;
    odoMapTrajectory(1,3) = 0;
    odoMapTrajectory(1,4) = 0;
    odoMapTrajectory(1,5) = 0;

    transOdoTraj(1,1) = 0;
    transOdoTraj(1,2) = 0;
    transOdoTraj(1,3) = 0;

    % load ground truth data
    [frameId, gt_x, gt_y, gt_z, gt_rx, gt_ry, gt_rz] = load_ground_truth_data(groundTruthFile);
%     groundTruthData = [gt_x gt_y gt_z];

    horiTransVHis = [];

    preHoriTransV = 0;
    preVertTransV = 0;
    preYawRotV = 0;
    prePitchRotV = 0;

    % for all frames do
    replace = 0;

    % end libviso2

    %% process
    for iSubFolder = 1 : numSubFolders

        [curFolderPath, imgFilesPathList, numImgs] = get_cur_img_files_path_list(subFoldersPathSet, IMG_TYPE, iSubFolder);

        if numImgs > 0
            for curFrame = 1 : numImgs

                [curImg] = read_current_image(curFolderPath, imgFilesPathList, curFrame);

                % visual templates and visual odo uses intensity so convert to grayscale
                curGrayImg = rgb2gray(curImg);
                curGrayImg = im2double(curGrayImg);

                if curFrame > 1

                    tempHoriTransV = sqrt((gt_x(curFrame) - gt_x(curFrame - 1))^2 + ...
                        (gt_y(curFrame) - gt_y(curFrame - 1))^2 );
                    if abs(tempHoriTransV) < MAX_HORI_TRANS_V_THRESHOLD
                        horiTransV = tempHoriTransV;
                        preHoriTransV = horiTransV;
                    else
                        horiTransV = preHoriTransV;
                    end

                    tempYawRotV = get_signed_delta_radian(gt_rz(curFrame), gt_rz(curFrame - 1));

                    if abs(tempYawRotV) < MAX_ROT_V_THRESHOLD
                        yawRotV = tempYawRotV;
                        preYawRotV = yawRotV;
                    else
                        yawRotV = preYawRotV;
                    end

                    tempPitchRotV = get_signed_delta_radian(gt_rx(curFrame),gt_rx(curFrame - 1));

                    if abs(tempPitchRotV) < MAX_PITCH_V_THRESHOLD
                        pitchRotV = tempPitchRotV;
                        prePitchRotV = pitchRotV;
                    else
                        pitchRotV = prePitchRotV;
                    end

                    odoPitchTheta(curFrame, 3) = odoPitchTheta(curFrame - 1, 3) + pitchRotV;
                    odoYawTheta(curFrame, 3) = odoYawTheta(curFrame -1, 3) + yawRotV;

                    tempVertTransV = sqrt((gt_x(curFrame) - gt_x(curFrame - 1))^2 + ...
                        (gt_y(curFrame) - gt_y(curFrame - 1))^2 + ...
                        (gt_z(curFrame) - gt_z(curFrame - 1))^2) * sin(odoPitchTheta(curFrame, 3));

                    if abs(tempVertTransV) < MAX_VERT_TRANS_V_THRESHOLD
                        vertTransV = tempVertTransV;
                        preVertTransV = vertTransV;
                    else
                        vertTransV = preVertTransV;
                    end

                else
                    horiTransV = 0;
                    pitchRotV = 0;
                    vertTransV = 0;
                    yawRotV = 0;

                    odoPitchTheta(curFrame, 3) = pitchRotV;
                    odoYawTheta(curFrame, 3) = yawRotV;

                end

                %% step3: processing vt
                % get the most active view template
                [vt_id] = visual_template(curGrayImg, gcX, gcY,gcZ, yawRotV, pitchRotV);

                %% step4: processing the integration of yaw_pitch_hdc
                yaw_pitch_hdc_iteration(vt_id, yawRotV, pitchRotV);

                [curYawTheta, curPitchTheta] = get_current_yaw_pitch_theta();

                % transform to radian
                curYawThetaInRadian = curYawTheta * YAW_PITCH_HDC_Y_TH_SIZE;
                curPitchThetaInRadian = curPitchTheta * YAW_PITCH_HDC_P_TH_SIZE;

                MAX_ACTIVE_YAW_PITCH_HIS_PATH = [MAX_ACTIVE_YAW_PITCH_HIS_PATH; curYawTheta curPitchTheta];

                % 3d grid cells interation
                gc_iteration(vt_id, horiTransV, curYawThetaInRadian, vertTransV, curPitchThetaInRadian);

                [gcX, gcY, gcZ] = get_gc_xyz();

                MAX_ACTIVE_XYZ_PATH = [MAX_ACTIVE_XYZ_PATH; gcX gcY gcZ];

                % 3d experience map interation
                exp_map_iteration(vt_id, horiTransV, yawRotV, vertTransV, pitchRotV, gcX, gcY, gcZ, curYawTheta, curPitchTheta);

                % for drawing HDC pitch theta
                hdcPitchTheta(curFrame, 1) = cos(curPitchTheta * YAW_PITCH_HDC_P_TH_SIZE);
                hdcPitchTheta(curFrame, 2) = sin(curPitchTheta * YAW_PITCH_HDC_P_TH_SIZE);
                hdcPitchTheta(curFrame, 3) = curPitchTheta * YAW_PITCH_HDC_P_TH_SIZE;

                if curFrame > 1
                    % for drawing odometry pitch theta
                    odoPitchTheta(curFrame, 1) = cos(odoPitchTheta(curFrame, 3));
                    odoPitchTheta(curFrame + 1, 2) = sin(odoPitchTheta(curFrame, 3));
    %                 odoPitchTheta(curFrame + 1, 3) = odoPitchTheta(curFrame, 3);

                    % for drawing HDC yaw theta
                    hdcYawTheta(curFrame, 1) = cos(curYawTheta * YAW_PITCH_HDC_Y_TH_SIZE);
                    hdcYawTheta(curFrame, 2) = sin(curYawTheta * YAW_PITCH_HDC_Y_TH_SIZE);
                    hdcYawTheta(curFrame, 3) = curYawTheta * YAW_PITCH_HDC_Y_TH_SIZE;

                    % for drawing odometry yaw theta
                    odoYawTheta(curFrame, 1) = cos(odoYawTheta(curFrame, 3));
                    odoYawTheta(curFrame, 2) = sin(odoYawTheta(curFrame, 3));
    %                 odoYawTheta(curFrame + 1, 3) = odoYawTheta(curFrame, 3) + yawRotV;

                    odoMapTrajectory(curFrame,1) = odoMapTrajectory(curFrame - 1, 1) + horiTransV * cos(odoYawTheta(curFrame, 3));
                    odoMapTrajectory(curFrame,2) = odoMapTrajectory(curFrame - 1, 2) + horiTransV * sin(odoYawTheta(curFrame, 3));
                    odoMapTrajectory(curFrame,3) = odoMapTrajectory(curFrame - 1, 3) + vertTransV;
                    odoMapTrajectory(curFrame,4) = odoYawTheta(curFrame, 3);
                    odoMapTrajectory(curFrame,5) = odoPitchTheta(curFrame, 3);
    %                 transOdoTraj(curFrame + 1,1) = odoMapTrajectory(curFrame + 1,1);
    %                 transOdoTraj(curFrame + 1,2) = odoMapTrajectory(curFrame + 1,2);
    %                 transOdoTraj(curFrame + 1,3) = odoMapTrajectory(curFrame + 1,3);

                end

                for ind = 1 : YAW_PITCH_HDC_Y_DIM
                    tempYawPitchHDC(:,:,ind) = YAW_PITCH_HDC;
                end

                for ind = 1 : NUM_EXPS
                    expTrajectory(ind,1) = EXPERIENCES(ind).x_exp;
                    expTrajectory(ind,2) = EXPERIENCES(ind).y_exp;
                    expTrajectory(ind,3) = EXPERIENCES(ind).z_exp;
                    expTrajectory(ind,4) = EXPERIENCES(ind).yaw_exp_rad;
                    expTrajectory(ind,5) = EXPERIENCES(ind).pitch_exp_rad;
%                     transExpTraj(ind,1) =   expTrajectory(ind,1);
%                     transExpTraj(ind,2) =   expTrajectory(ind,2);
%                     transExpTraj(ind,3) =   expTrajectory(ind,3);
                end

                % draw the results
                if (mod(curFrame, RENDER_RATE) == 1)

%
                    % drawing the activity of yaw_pitch_hdc
                    figure('color',[1 1 1]);
                    subplot(1, 1, 1,'replace');
                    [x,y]= meshgrid(1:1:36);
                    z = YAW_PITCH_HDC ;
                    surf(x,y,z,'EdgeColor','b');
                    % surf(x,y,z,'FaceColor','red','EdgeColor','none');
%                     shading interp
                    axis([1 YAW_PITCH_HDC_Y_DIM 1 YAW_PITCH_HDC_P_DIM 0 0.02]);
                    % meshc(x,y,z);
                    view(0, 90);
                    xl = xlabel('Pitch');
                    yl = ylabel('Yaw');
                    zlabel('Acitivity');
%                     set(xl,'Rotation',15);
%                     set(yl,'Rotation',-30);
%                     title('Yaw Pitch HDCell Activity');
                    axis on
%

                    drawnow;

                end

                PREV_VT_ID = vt_id;
            end

            % save history data (x, y, z, yaw, pitch)
%             save_history_data(expMapHistoryFile, expTrajectory);
%             save_history_data(odoMapHistoryFile, odoMapTrajectory);
%             save_em_history_data(emHistoryFile, VT_HISTORY);
%             save_vt_history_data(vtHistoryFile, EXP_HISTORY);
        end
    end
end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function im_rec = tomo_recon_lsqr(W, p, siz, tol, n_it)
%TOMO_RECON_LSQR   Tomographic reconstruction using LSQR method.
%
%   Phymhan
%   09-Aug-2013 09:53:51

if nargin < 5
    n_it = 100;
    if nargin < 4
        tol = 1e-6;
    end
end
%siz = size(im);
%[W, p, ~, ~] = buildWeightMatrix(im,angles);
f = lsqr(W,p,tol,n_it);
im_rec = reshape(f,siz);
end


text:
function y = G2a(X)
% mca for rtd function
% reduced order A matrix - no opamps
% X = [R1 R2 R3 R4 R5 R6 R7 R8 R9 RT E1]
R1=X(1);R2=X(2);R3=X(3);R4=X(4);R5=X(5);R6=X(6);
R7=X(7);R8=X(8);R9=X(9);RT=X(10);E1=X(11);
%
A=[1/R1+1/R4+1/RT -1/RT -1/R4 0 -1/R1;
-1/RT 1/R5+1/R6+1/RT -1/R5 0 0;
-1/R4 0 1/R2+1/R3+1/R4 -1/R3 0;
0 -1/R5 1/R5+1/R7 0 0;
0 0 0 0 -1/R9];
%
B=[0;E1/R6;E1/R2;0;E1/R8];
C=A\B;y=C(4);


text:
% designsim_gui_script
% Last modified 3/21/01 by Tor Wager		fix HRF bug / various
name = input('Enter output file name: ');

% sets defaults, reads GUI values from optimizeGUI, and runs simulation

disp('Designsim_gui_script')
disp('===============================')
disp('  ...Setting up variables from GUI.')
% returns vector of t values as t in workspace.

if ~exist('M') | ~isstruct(M),disp('No GA results (M struct) detected - using only random designs.'),end

clear SIM
clear H;

format compact

% defaults
noise_var = 0.3;
beta = [.5 .5 0 0;1 1 0 0];	% 20
c = [1 1 -1 -1];					% 13
niterations = 500;
conditions = [1 2 3 4];
ISI = 1.2;
restevery = 9;						% 26
restlength = 2;					% 25
HPlength = 60;						% 23
xc = [1 0];							% 24

% read values from GUI for GA
H = findobj('Tag', 'E1');
conditions = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E4');
ISI = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E12');
noise_var = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E13');
c = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E14');
niterations = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E15');
respMean = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E16');
respSD = str2num(get(H(1), 'String'));

H = findobj('Tag', 'C1');
addResponseVariability = get(H(1), 'Value');

% additional values required for random sim.

H = findobj('Tag', 'E2');
freqConditions = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E3');
scanLength = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E5');
TR = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E6');
cbalColinPowerWeights = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E7');
numGenerations = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E8');
sizeGenerations = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E9');
lowerLimit = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E10');
maxOrder = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E11');
plotFlag = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E20');
beta = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E23');
HPlength = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E24');
eval(['load ' get(H(1), 'String')]);

H = findobj('Tag', 'E25');
restlength = str2num(get(H(1), 'String'));

H = findobj('Tag', 'E26');
restevery = str2num(get(H(1), 'String'));

HRF = spm_hrf(.1);
HRF = HRF/ max(HRF);

for myISI = 1:16
ISI = myISI;
disp(['Starting simulation with ISI ' num2str(myISI)]);

SIM.beta = beta;
SIM.contrasts = c;
SIM.noise_var = noise_var;
SIM.xc = xc;
SIM.ISI = ISI;
SIM.HPlength = HPlength;
SIM.restevery = restevery;
SIM.restlen = restlength;
SIM.niterations = niterations;
SIM.noise_var = noise_var;
SIM.freqConditions = freqConditions;

disp('  ...Generating design list.')
% for building RANDOM designs to test against GA design
fitnessMatrix = zeros(4,sizeGenerations);
numStim = ceil(scanLength / (ISI));
numRestStim = (ceil(numStim/restevery) - 1) * restlength;
numStimEachCond = ceil((numStim-numRestStim) * freqConditions);			            % row vector of stim in each cond
numsamps = ceil(numStim*ISI/TR);

SIM.numsamps = numsamps;
if exist('M') & isstruct(M), SIM.numsamps = size(M.modelatTR,1);,end
SIM.numStimEachCond = numStimEachCond;

unsortedList = [];
for i = 1:size(conditions,2)									% unsorted list of all stims in proper frequencies
   startIndex = size(unsortedList,1)+1;
   unsortedList(startIndex:(startIndex + numStimEachCond(i)-1),1) = conditions(i);
end

if addResponseVariability == 1
    % THIS ISN'T FINISHED, AND I'M NOT SURE I SEE THE POINT, REALLY.

else
    disp(['  ...loading smoothing matrix with hrf LP and HP filter length of ' num2str(HPlength)])
    [hpS] = use_spm_filter(TR,SIM.numsamps,'none','specify',HPlength);
	[fullS] = use_spm_filter(TR,SIM.numsamps,'hrf','specify',HPlength);

	if exist('M') & isstruct(M),
	        stimList = M.stimlist(1:numsamps+1,1);   % cut off to make more efficient
                stimList = sampleInSeconds(stimList,ISI,.1);
		dmodel = getPredictors(stimList, HRF);
		dmodel = resample(dmodel,1,TR*10);
		% dmodel = dmodel / max(max(dmodel));                     % scale so that max = 1;
		if size(dmodel,1) > SIM.numsamps
			dmodel = dmodel(1:SIM.numsamps,:);
		elseif size(dmodel,1) < SIM.numsamps
		        warning('stimlist for GA is too short for this ISI!  Results may not be accurate for this ISI!')
	                dmodel = [dmodel; zeros(SIM.numsamps - size(dmodel,1),size(dmodel,2))];
	        end

  		disp(['  ...simulating GA results with noisy data at variance ' num2str(noise_var)])
                SIM
	        % Test idealized GA - no response variability
		fprintf('		...original	model')
    	[dummy,SIM.ga.se,SIM.ga.t] = designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,0,xc,niterations);
 		fprintf('...high-pass filtered')
    	[dummy,SIM.ga.HPse,SIM.ga.HPt] = designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,hpS,xc,niterations);
		fprintf('...smoothed\n')
		[dummy,SIM.ga.HLse,SIM.ga.HLt] = designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,fullS,xc,niterations);
	end

    disp(['  ...simulating random designs with ' num2str(niterations) ' iterations...']);

SIM
whos unsortedList

	% to do:
% cut off model at numscans  done, i think
% rebuild ga at different isis  done
% cut off stimlist beforehand cause it's too long.   done
% warning if stimlist in GA is too short for selected isi - then pad with zeros.  done

	% computing random SIM model
	% ============================================================================================

	for i = 1:niterations
      % Test Random models

	  	% Prepare Stimulus List
      	stimList = getRandom(unsortedList);
      	% ====== insert rests, if needed, and trim list ===================
        stimList = insert_rests(stimList,restevery,restlength);
        if size(stimList,1) < numStim, stimList(end+1:numStim,:) = 0;
        elseif size(stimList,1) > numStim,stimList = stimList(1:numStim,:);
        end
		stimList = sampleInSeconds(stimList,ISI,.1);
		dmodel = getPredictors(stimList, HRF);
		dmodel = resample(dmodel,1,TR*10);
		% dmodel = dmodel / max(max(dmodel));                     % scale so that max = 1;
		if size(dmodel,1) > SIM.numsamps
			dmodel = dmodel(1:SIM.numsamps,:);
		elseif size(dmodel,1) < SIM.numsamps
			dmodel = [dmodel; zeros(SIM.numsamps - size(dmodel,1),size(dmodel,2))];
		end

		% this gets the average over [numnoise] simulations with different noise vectors
	  % =====================================================================================================
      % no smoothing or filtering
      [dummy,SIM.rnd.se(:,i),SIM.rnd.t(:,i)] = designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,0,xc);
	  % high-pass only
	  [dummy, SIM.rnd.HPse(:,i),SIM.rnd.HPt(:,i)] = designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,hpS,xc);
      % now with smoothing
      [dummy, SIM.rnd.HLse(:,i),SIM.rnd.HLt(:,i)] = designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,fullS,xc);
      fprintf('.')
	  if mod(i,30) == 0, fprintf(' %d ',i),end
	  if mod(i,90) == 0, fprintf('\n'),end

	  % this gets the average over [numnoise] simulations with different noise vectors
	  % =====================================================================================================
      % no smoothing or filtering
      %[SIM.rnd.t(:,i), SIM.rnd.se(:,i)] = designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,0,xc,100);
	  % high-pass only
	  %[SIM.rnd.HPt(:,i), SIM.rnd.HPse(:,i)] = designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,hpS,xc,100);
      % now with smoothing
      %[SIM.rnd.HLt(:,i), SIM.rnd.HLse(:,i)] = designsim('myxc',dmodel,HRF,ISI,TR,noise_var,c,beta,fullS,xc,100);
      %fprintf('.')
      %if mod(i,30) == 0, fprintf(' %d\n',i),end
      %if mod(i,50) == 0, disp(['done ' num2str(i)]),end;

    end

end

      MyISIMean.ga.t(myISI) = mean(SIM.ga.t);
      MyISIMean.rnd.t(myISI) = mean(SIM.rnd.t);
      MyISIMean.ga.HPt(myISI) = mean(SIM.ga.HPt);
      MyISIMean.rnd.HPt(myISI) = mean(SIM.rnd.HPt);
      MyISIMean.ga.HLt(myISI) = mean(SIM.ga.HLt);
      MyISIMean.rnd.HLt(myISI) = mean(SIM.rnd.HLt);
      MyISIMean.ISI(myISI) = ISI;
      MyISIMean.ga.se(myISI) = std(SIM.ga.t)/sqrt(size(SIM.ga.t,2));
      MyISIMean.rnd.se(myISI) = std(SIM.rnd.t)/sqrt(size(SIM.rnd.t,2));
      MyISIMean.ga.HPse(myISI) = std(SIM.ga.HPt)/sqrt(size(SIM.ga.HPt,2));
      MyISIMean.rnd.HPse(myISI) = std(SIM.rnd.HPt)/sqrt(size(SIM.rnd.HPt,2));
      MyISIMean.ga.HLse(myISI) = std(SIM.ga.HLt)/sqrt(size(SIM.ga.HLt,2));
      MyISIMean.rnd.HLse(myISI) = std(SIM.rnd.HLt)/sqrt(size(SIM.rnd.HLt,2));

end % end loop through simulations
      eval(['save ' name ' MyISIMean'])

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function degree = grNodeInnerDegree(node, edges)
%GRNODEINNERDEGREE Inner degree of a node in a graph
%
%   DEG = grNodeInnerDegree(NODE, EDGES);
%   Returns the inner degree of a node in the given edge list, i.e. the
%   number of edges arriving to it.
%   NODE is the index of the node, and EDGES is a liste of couples of
%   indices (origin and destination node).
%
%   Note: Also works when node is a vector of indices
%
%   See Also:
%   grNodeDegree, grNodeOuterDegree
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2006-01-17
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).
%

%   HISTORY
%   2008-08-07 pre-allocate memory, update doc

% allocate memory
N = size(node, 1);
degree = zeros(N, 1);

% compute inner degree of each vertex
for i=1:length(node)
    degree(i) = sum(edges(:,2)==node(i));
end


text:
function xlsheets(sheetnames,varargin)

%XLSHEETS creates or opens existing Excel file and names sheets
%
% xlsheets(sheetnames,filename)
% xlsheets(sheetnames)
%
% xlsheets  : Creates new excel file (or opens it if file exists)
%               and name the sheets as listed in (sheetnames)
%               and saves the workbook as (filename).
%
%       sheetnames:     List of sheet names (cell array).
%       filename:       Name of excel file.
%
% NOTE: Follow the following rules when naming your sheets:
%       1- Make sure the name you entered does not exceed 31 characters.
%       2- Make sure the name does not contain any of the following characters:  :  \  /  ?  *  [  or  ]
%       3- Make sure you did not leave the name blank.
%       4- Make sure each sheet name is a character string.
%       5- Make sure you do not have two sheets or more with the same name.
%
% Example:
%
%      sheetnames = {'Mama','Papa','Son','Daughter','Dog'};
%      filename = 'family.xls';          % can be named without '.xls'
%      xlsheets(sheetnames,filename);
%      xlsheets(sheetnames);            % Will leave file open
%

%   Copyright 2004 Fahad Al Mahmood
%   Version: 1.0 $  $Date: 12-Feb-2004
%   Version: 1.5 $  $Date: 16-Feb-2004  (Open exisiting file feature)
%   Version: 2.0 $  $Date: 26-Feb-2004  (Fixed [Group] problem + Making process invisible)
%   Version: 2.1 $  $Date: 27-Feb-2004  (Fixed replacing existing sheets problem)
%   Version: 2.5 $  $Date: 15-Mar-2004  (Fixed filename problem)
%   Version: 3.0 $  $Date: 04-Apr-2004  (Fixed Naming to an existing sheetnames problem + Fixed Opening Multiple Excel Programs Problem)
%   Version: 3.1 $  $Date: 10-Apr-2004  (Added more help about the rules of naming Excel sheets)
%   Version: 3.2 $  $Date: 10-Apr-2004  (Supporting Full or Partial Path)

% Making sure the names of the sheets are according to Excel rules.
for n=1:length(sheetnames)
%  (1) Making sure each sheetname entered does not exceed 31 characters.
    if length(sheetnames{n})>31
        error(['sheet (' sheetnames{n} ') exceeds 31 characters! (see xlsheets help)'])
    end
%  (2) Making sure each sheetname does not contain any illegal character.
    if any(ismember([':','\','/','?','*'],sheetnames{n})) | ismember('[',sheetnames{n}(1))
        error(['sheet (' sheetnames{n} ') contains an illegal character! (see xlsheets help)'])
    end
%  (3) Making sure each sheetname is not blank.
    if isempty(sheetnames{n})
        error(['sheet ' int2str(n) ' is empty! (see xlsheets help)'])
    end
%  (4) Making sure each sheetname is a character string.
    if ~ischar(sheetnames{n})
        error(['sheet (' int2str(n) ') is NOT a character string! (see xlsheets help)'])
    end
end

%  (5) Making sure two or more sheets do not have the same name.
if length(sheetnames)>length(unique(sheetnames))
    error('Two or more sheets have the same name!')

end

% Opening Excel
target_num_sheets = length(sheetnames);
Excel = actxserver('Excel.Application');
if nargin==2
    filename = varargin{1};
    [fpath,fname,fext] = fileparts(filename);
    if isempty(fpath)
        out_path = pwd;
    elseif fpath(1)=='.'
        out_path = [pwd filesep fpath];
    else
        out_path = fpath;
    end
    filename = [out_path filesep fname fext];
    if ~exist(filename,'file')
        % The following case if file does not exist (Creating New Workbook)
        Workbook = invoke(Excel.Workbooks,'Add');
        % getting the number of sheets in new workbook
        numsheets = get(Excel,'SheetsInNewWorkbook');
        new=1;
    else
        % The following case if file does exist (Opening Workbook)
        Workbook = invoke(Excel.Workbooks, 'open', filename);
        % getting the number of sheets in new workbook
        workSheets = Excel.sheets;
        for i = 1:workSheets.Count
            sheet = get(workSheets,'item',i);
            description{i} = sheet.Name;
            if ~isempty(sheet.UsedRange.value)
                indexes(i) = true;
            else
                indexes(i) = false;
            end
        end
        descr = description(indexes);
        numsheets = length(descr);
        new=0;
    end
    leave_file_open = 0;
else
    % The following case if file does not exist (Creating New Workbook)
    Workbook = invoke(Excel.Workbooks,'Add');
    % getting the number of sheets in new workbook
    numsheets = get(Excel,'SheetsInNewWorkbook');
    new=1;
    leave_file_open = 1;
end

% making Excel visible only if workbook name is not specified or new workbook is created.
if nargin==1
    set(Excel,'Visible', 1);
end

if target_num_sheets > numsheets

    % Activating Last sheet of new (filename)
    Sheets = Excel.ActiveWorkBook.Sheets;
    sheet = get(Sheets, 'Item', numsheets);
    invoke(sheet, 'Activate');

    % Adding sheets to match the number of (sheetnames) specified.
    for i=1:target_num_sheets-numsheets
        invoke(Excel.Sheets,'Add');
    end

elseif target_num_sheets < numsheets

    % Deleting sheets to match the number of (sheetnames) specified.
    for i=numsheets-target_num_sheets:-1:1
        sheet = get(Excel.ActiveWorkBook.Sheets, 'Item', i);
        invoke(sheet, 'Delete');
    end
end

% Renaming sheets to temporary names
for i=1:target_num_sheets
    Sheets = Excel.Worksheets;
    sheet = get(Sheets, 'Item', i);
    invoke(sheet, 'Activate');
    Activesheet = Excel.Activesheet;
    temp_name = ['temp_' int2str(i)];
    set(Activesheet,'Name',temp_name);
end

% Renaming sheets to the designated names
for i=1:target_num_sheets
    Sheets = Excel.Worksheets;
    sheet = get(Sheets, 'Item', i);
    invoke(sheet, 'Activate');
    Activesheet = Excel.Activesheet;
    set(Activesheet,'Name',char(sheetnames(i)));
end

if nargin>1
    if new invoke(Workbook, 'SaveAs', filename);
    else invoke(Workbook, 'Save'); end
end

if ~leave_file_open invoke(Excel, 'Quit'); end
delete(Excel);


text:
function [mesh, cfg] = ft_prepare_mesh(cfg, data)

% FT_PREPARE_MESH creates a triangulated surface mesh or tetrahedral/hexahedral
% volume mesh that can be used as geometrical description for a volume conduction
% model. The mesh can either be created manually from anatomical MRI data or can be
% generated starting from a segmented MRI. This function can also be used to create a
% cortex hull, i.e. the smoothed envelope around the pial surface created by
% freesurfer.
%
% Use as
%   mesh = ft_prepare_mesh(cfg)
%   mesh = ft_prepare_mesh(cfg, mri)
%   mesh = ft_prepare_mesh(cfg, seg)
%
% Configuration options:
%   cfg.method      = string, can be 'interactive', 'projectmesh', 'iso2mesh', 'isosurface',
%                     'headshape', 'hexahedral', 'tetrahedral','cortexhull', 'fittemplate'
%   cfg.tissue      = cell-array with strings representing the tissue types, or numeric vector with integer values
%   cfg.numvertices = numeric vector, should have same number of elements as the number of tissues
%
% When providing an anatomical MRI or a segmentation, you should specify
%   cfg.downsample  = integer number (default = 1, i.e. no downsampling), see FT_VOLUMEDOWNSAMPLE
%   cfg.spmversion  = string, 'spm2', 'spm8', 'spm12' (default = 'spm8')
%
% For method 'headshape' you should specify
%   cfg.headshape   = a filename containing headshape, a Nx3 matrix with surface
%                     points, or a structure with a single or multiple boundaries
%
% For method 'cortexhull' you should not give input data, but specify
%   cfg.headshape   = string, filename containing the pial surface computed by freesurfer recon-all
%
% For method 'fittemplate' you should specify
%   cfg.headshape   = a filename containing headshape
%   cfg.template    = a filename containing headshape
% With this method you are fitting the headshape from the configuration to the template;
% the resulting affine transformation is applied to the input mesh (or set of meshes),
% which is subsequently returned as output variable.
%
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% Example
%   mri             = ft_read_mri('Subject01.mri');
%
%   cfg             = [];
%   cfg.output      = {'scalp', 'skull', 'brain'};
%   segmentation    = ft_volumesegment(cfg, mri);
%
%   cfg             = [];
%   cfg.tissue      = {'scalp', 'skull', 'brain'};
%   cfg.numvertices = [800, 1600, 2400];
%   mesh            = ft_prepare_mesh(cfg, segmentation);
%
%   cfg             = [];
%   cfg.method      = 'cortexhull';
%   cfg.headshape   = '/path/to/surf/lh.pial';
%   cfg.fshome      = '/path/to/freesurfer dir';
%   cortex_hull     = ft_prepare_mesh(cfg);
%
% See also FT_VOLUMESEGMENT, FT_PREPARE_HEADMODEL, FT_PLOT_MESH

% Undocumented functionality: at this moment it allows for either
%   mesh = ft_prepare_mesh(cfg)             or
%   mesh = ft_prepare_mesh(cfg, headmodel)
% but more consistent would be to specify a volume conduction model with
%   cfg.headmodel     = structure with volume conduction model, see FT_PREPARE_HEADMODEL
%   cfg.headshape     = name of file containing the volume conduction model, see FT_READ_HEADMODEL
%
% Undocumented options, I have no clue why they exist
%   cfg.method = {'singlesphere' 'concentricspheres' 'localspheres'}

% Copyrights (C) 2009-2012, Robert Oostenveld & Cristiano Micheli
% Copyrights (C) 2012-2013, Robert Oostenveld & Lilla Magyari
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data
ft_preamble trackconfig

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% we cannot use nargin, because the data might have been loaded from cfg.inputfile
hasdata = exist('data', 'var');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden', {'numcompartments', 'outputfile', 'sourceunits', 'mriunits'});

% get the options
cfg.downsample  = ft_getopt(cfg, 'downsample', 1);      % default is no downsampling
cfg.numvertices = ft_getopt(cfg, 'numvertices', 3000);  % set the default
cfg.smooth      = ft_getopt(cfg, 'smooth');             % no default
cfg.spmversion  = ft_getopt(cfg, 'spmversion', 'spm8');

% Translate the input options in the appropriate default for cfg.method
if isfield(cfg, 'headshape') && ~isempty(cfg.headshape)
  cfg.method = ft_getopt(cfg, 'method', 'headshape');
elseif hasdata && ~strcmp(ft_headmodeltype(data), 'unknown')
  cfg.method = ft_getopt(cfg, 'method', ft_headmodeltype(data));
elseif hasdata
  cfg.method = ft_getopt(cfg, 'method', 'projectmesh');
else
  cfg.method = ft_getopt(cfg, 'method', []);
end

if hasdata && cfg.downsample~=1
  % optionally downsample the anatomical volume and/or tissue segmentations
  tmpcfg = keepfields(cfg, {'downsample', 'showcallinfo'});
  data = ft_volumedownsample(tmpcfg, data);
  % restore the provenance information and put back cfg.smooth
  tmpsmooth = cfg.smooth;
  [cfg, data] = rollback_provenance(cfg, data);
  cfg.smooth = tmpsmooth;
end

switch cfg.method
  case 'interactive'
    % this makes sense with a non-segmented MRI as input
    % call the corresponding helper function
    mesh = prepare_mesh_manual(cfg, data);

  case {'projectmesh', 'iso2mesh', 'isosurface'}
    % this makes sense with a segmented MRI as input
    % call the corresponding helper function
    mesh = prepare_mesh_segmentation(cfg, data);

  case 'headshape'
    % call the corresponding helper function
    mesh = prepare_mesh_headshape(cfg);

  case 'hexahedral'
    % the MRI is assumed to contain a segmentation
    % call the corresponding helper function
    mesh = prepare_mesh_hexahedral(cfg, data);

  case 'tetrahedral'
    % the input data is assumed to contain a segmentation
    % call the corresponding helper function
    mesh = prepare_mesh_tetrahedral(cfg, data);

  case {'singlesphere' 'concentricspheres'}
    % the input data is assumed to contain a spherical head model
    data = ft_datatype_headmodel(data);   % ensure that it is consistent and up-to-date
    data = ft_determine_units(data);      % ensure that it has units

    if length(data.r)>1 && numel(cfg.numvertices)==1
      % use the same number of vertices for each tissue
      cfg.numvertices = repmat(cfg.numvertices, length(data.r));
    end

    mesh = [];
    for i=1:length(data.r)
      ft_info('triangulating sphere %d in the volume conductor\n', i);
      [pos, tri] = mesh_sphere(cfg.numvertices(i));
      mesh(i).pos(:,1) = pos(:,1)*data.r(i) + data.o(1);
      mesh(i).pos(:,2) = pos(:,2)*data.r(i) + data.o(2);
      mesh(i).pos(:,3) = pos(:,3)*data.r(i) + data.o(3);
      mesh(i).tri = tri;
    end

  case 'cortexhull'
    mesh = prepare_mesh_cortexhull(cfg);

  case 'fittemplate'
    M   = prepare_mesh_fittemplate(cfg.headshape.pos, cfg.template.pos);
    orig.mri = data;
    orig = ft_transform_geometry(M, orig);
    mesh = orig.data;

  otherwise
    ft_error('unsupported cfg.method')
end

% copy the geometrical units from the input to the output
if ~isfield(mesh, 'unit') && hasdata && isfield(data, 'unit')
  for i=1:numel(mesh)
    mesh(i).unit = data.unit;
  end
elseif ~isfield(mesh, 'unit')
  mesh = ft_determine_units(mesh);
end

% copy the coordinate system from the input to the output
if ~isfield(mesh, 'coordsys') && hasdata && isfield(data, 'coordsys')
  for i=1:numel(mesh)
    mesh(i).coordsys = data.coordsys;
  end
end

% smooth the mesh
if ~isempty(cfg.smooth)
  cfg.headshape = mesh;
  cfg.numvertices = [];
  mesh = prepare_mesh_headshape(cfg);
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous   data
ft_postamble provenance mesh
ft_postamble history    mesh

author: spm


text:
function [dicomInfo] = SiemensCsaParse(dicomInfo)
% Read values from the 'SIEMENS CSA HEADER' private group in a
% DICOM file produced by a Siemens MR scanner.
%
% [dicomInfo] = SiemensCsaParse(dicomInfo)
%
% dicomInfo is a Siemens MR DICOM file loaded using Matlab's dicominfo()
%        OR a string specifying the filename of a DICOM file to load.
%
% output is the original structure with an extra field 'csa' that contains
% the decoded 'SIEMENS CSA HEADER' block

% Inspiration for this comes from the gdcm tool (http://gdcm.sf.net/),
% specifically gdcmCSAHeader.cxx lines 490-65.

% Possible TODO:
% 1. Migrate to C++ MEX file.
% 2. Investigate whether there's a Siemens product .dll file (e.g.
%    Z:\n4\pkg\MrServers\MrProtSrv) for parsing these headers.
% 3. Understand meaning of syngodt (Syngo Datatype ???) field.

% Copyright: Chris Rodgers (University of Oxford), 2008-11.
% All rights reserved.

% $Id: SiemensCsaParse.m 5022 2012-02-01 11:26:39Z will $

if ischar(dicomInfo)
    dicomInfo = dicominfo(dicomInfo);
end

%% Locate SIEMENS CSA HEADER
tagId = getDicomPrivateTag(dicomInfo,'0029','SIEMENS CSA HEADER');

if isempty(tagId)
    % Give warning and return with unmodified dicomInfo
    warning('RodgersSpectroTools:NoSiemensCsaHeader','This DICOM file does not contain a ''SIEMENS CSA HEADER''.')
    return
end

%% Process the two parameter fields
SiemensCsaParse_ReadDicomTag(['Private_0029_' tagId '10']);
SiemensCsaParse_ReadDicomTag(['Private_0029_' tagId '20']);

return

%% Main code
function SiemensCsaParse_ReadDicomTag(strTag)
    currdx=0;

    if ~strcmp(char(private_read(4)),'SV10') || ~all(private_read(4)==[4 3 2 1])
        error('Unsupported CSA block format');
    end

    % This parsing code is translated from gdcm (http://gdcm.sf.net/)
    numElements = double(private_readuint32(1));

    % Sanity check
    if private_readuint32(1)~=77
        error('Unsupported CSA block format');
    end

    for tagdx=1:numElements
        tagName = private_readstring(64);

        % Fix up tagName
        tagName(tagName == '-') = [];

        vm = private_readuint32(1);
        vr = private_readstring(4);
        syngodt = private_readuint32(1);
        nitems = double(private_readuint32(1));

        checkbit = private_readuint32(1);

        if checkbit ~= 77 && checkbit ~= 205
            error('Unsupported CSA block format');
        end

        data = {};
        for itemdx=1:nitems
            header = double(private_readuint32(4));

            if (header(3) ~= 77 && header(3) ~= 205) || ...
                    (header(1) ~= header(2)) || ...
                    (header(1) ~= header(4))
                error('Unsupported CSA block format');
            end

            data{itemdx} = private_readstring(header(1));

            % Dump junk up to DWORD boundary
            private_read(mod(mod(4-header(1),4),4));
        end

        % Store this in the csa structure
        switch vr
            case {'CS', 'LO', 'LT', 'SH', 'SS', 'UI', 'UT', 'UN'} % Strings and unknown byte string
                if numel(data) < vm
                    % Pad if necessary. Siemens CSA format omits null strings.
                    data{vm} = '';
                end

                if vm == 1
                    dicomInfo.csa.(tagName) = data{1};
                else
                    dicomInfo.csa.(tagName) = data(1:vm);
                end
            case {'DS', 'FD', 'FL', 'IS', 'SL', 'ST', 'UL', 'US'} % Numbers
                dataNumeric = arrayfun(@str2double,data);

                if numel(dataNumeric) < vm
                    % Zero pad if necessary. Siemens CSA format omits zeros.
                    dataNumeric(vm) = 0;
                end

                dicomInfo.csa.(tagName) = dataNumeric(1:vm);
            otherwise
                warning('RodgersSpectroTools:UnknownVrType','Unknown VR type: "%s".',vr)
        end
    end

    %% Helper functions to simulate file I/O
    function [out] = private_read(numBytes)
        out = dicomInfo.(strTag)(currdx+(1:numBytes)).';
        currdx=currdx+numBytes;
    end

    function [out] = private_readuint32(num)
        out=typecast(private_read(4*num),'uint32');
    end

    function [out] = private_readstring(maxchar)
        out = reshape(char(private_read(maxchar)),1,[]);
        terminator = find(out==0,1);
        if numel(terminator)>0
            out=out(1:(terminator-1));
        end
    end

end

end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function scr = merge(param1,param2)
% scr = merge(scr_obj_array)
% scr = merge(scr1,scr2)
% Merges Scores objects.  The function can be called in two ways: with
% two Scores objects or with an object array of Scores.
% Inputs:
%   Either:
%   scr1: A Scores object.
%   scr2: Another Scores object.
%   Or:
%   scr_obj_array: An array of Scores objects.
% Outputs:
%   scr: A Scores object that contains the information from the
%     input Scores objects.

assert(nargin<3)
if nargin==2
    scr1 = param1;
    scr2 = param2;
    assert(isa(scr1,'Scores'))
    assert(isa(scr2,'Scores'))
    assert(scr1.validate())
    assert(scr2.validate())

    scr_obj_arr = Scores.empty(2,0);
    scr_obj_arr(1) = scr1;
    scr_obj_arr(2) = scr2;
    scr = Scores.merge(scr_obj_arr);

else
    % the output scr must have all models and segment in the input
    % scrs (only once) and the union of all the scoremasks.
    % It is an error if two of the input Scores objects have a
    % score for the same trial.
    assert(nargin==1)
    obj_array = param1;
    assert(isa(obj_array,'Scores'))
    numscrs = length(obj_array);
    scr = Scores();
    for ii=1:numscrs
	scr_new = Scores();
	scr1 = scr;
	scr2 = obj_array(ii);

	% create new scr with empty matrices
	scr_new.modelset = union(scr1.modelset, scr2.modelset);
	scr_new.segset = union(scr1.segset, scr2.segset);

	% expand scr1 matrices
	scoremat_1 = zeros(length(scr_new.modelset), length(scr_new.segset));
	scoremask_1 = false(length(scr_new.modelset), length(scr_new.segset));
	[dummy,model_index_a,model_index_b] = intersect(scr_new.modelset,scr1.modelset);
	[dummy,seg_index_a,seg_index_b] = intersect(scr_new.segset,scr1.segset);
	scoremat_1(model_index_a,seg_index_a) = scr1.scoremat(model_index_b,seg_index_b);
	scoremask_1(model_index_a,seg_index_a) = scr1.scoremask(model_index_b,seg_index_b);

	% expand scr2 matrices
	scoremat_2 = zeros(length(scr_new.modelset), length(scr_new.segset));
	scoremask_2 = false(length(scr_new.modelset), length(scr_new.segset));
	[dummy,model_index_a,model_index_b] = intersect(scr_new.modelset,scr2.modelset);
	[dummy,seg_index_a,seg_index_b] = intersect(scr_new.segset,scr2.segset);
	scoremat_2(model_index_a,seg_index_a) = scr2.scoremat(model_index_b,seg_index_b);
	scoremask_2(model_index_a,seg_index_a) = scr2.scoremask(model_index_b,seg_index_b);

	% check for clashes
	assert(~any(scoremask_1(:)&scoremask_2(:)))

	% merge masks
	scoremat = scoremat_1 + scoremat_2;
	scoremask = scoremask_1 | scoremask_2;

	% build new scr
	scr_new.scoremat = scoremat;
	scr_new.scoremask = scoremask;
	scr = scr_new;
    end
end

assert(scr.validate())

author: nesl


text:
function [ti,fi]=midpoint(t1,f1,t2,f2,k)
%MIDPOINT Mid-point construction used in the interference diagram.
%	[TI,FI]=MIDPOINT(T1,F1,T2,F2,K) gives the coordinates in the
%	time-frequency plane of the interference-term corresponding to
%	the points (T1,F1) and (T2,F2), for a distribution in the
%	affine class perfectly localized on power-law group-delays of
%	the form : tx(nu)=t0+c nu^(K-1).
%
%	T1 : time-coordinate of the first point
%	F1 : frequency-coordinate of the first point (>0)
%	T2 : time-coordinate of the second point
%	F2 : frequency-coordinate of the second point (>0)
%	K  : power of the group-delay law
%	  K = 2    : Wigner-Ville
%	  K = 1/2  : D-Flandrin
%	  K = 0    : Bertrand (unitary)
%	  K = -1   : Unterberger (active)
%	  K = inf  : Margenau-Hill-Rihaczek
%	TI : time-coordinate of the interference term
%	FI : frequency-coordinate of the interference term
%
%	See also PLOTSID.

%	P. Flandrin, September 1995 - F. Auger, April 1996.
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if f1<=0 | f2<=0,
 error('F1 and F2 must be >0');
end
[rt1,ct1]=size(t1);
[rt2,ct2]=size(t2);
[rf1,cf1]=size(f1);
[rf2,cf2]=size(f2);
if (rt1~=rt2|rt1~=rf1|rt1~=rf2) | (ct1~=ct2|ct1~=cf1|ct1~=cf2),
 error('T1, T2, F1 and F2 must have the same size');
end
if rt1>ct1,
 error('T1 must be a row-vector');
elseif rt2>ct2,
 error('T2 must be a row-vector');
elseif rf2>cf2,
 error('F2 must be a row-vector');
elseif rf1>cf1,
 error('F1 must be a row-vector');
end

if (k==2),
 fi=(f1+f2)/2;
 ti=(t1+t2)/2;
elseif (k==inf),
 ti=[t1;t2];
 fi=[f2;f1];
else
 I=find(abs(f1-f2)>sqrt(eps));
 if length(I)~=0,
  if (k==1),
   fi(I)=exp( (f1(I).*(log(f1(I))-1)-f2(I).*(log(f2(I))-1)) ./ ...
 	(f1(I)-f2(I)));
   ti(I)=(t1(I).*f1(I)-t2(I).*f2(I)) ./ (f1(I)-f2(I)) - ...
 	(t1(I)-t2(I)) ./ (log(f1(I))-log(f2(I)));
  elseif (k==0),
   fi(I)=(f1(I)-f2(I))./(log(f1(I))-log(f2(I)));
   ti(I)=(t1(I).*f1(I)-t2(I).*f2(I)) ./ (f1(I)-f2(I)) + ...
         f1(I) .* f2(I) .* (t2(I)-t1(I)) .* ...
    	(log(f1(I))-log(f2(I))) ./ (f2(I)-f1(I)).^2;
  else
   t0(I)=(t1(I).*f2(I).^(k-1)-t2(I).*f1(I).^(k-1)) ./ ...
 	(f2(I).^(k-1)-f1(I).^(k-1));
   fi(I)=((f1(I).^k-f2(I).^k) ./ (f1(I)-f2(I))/k).^(1/(k-1));
   ti(I)=t0(I)+(t2(I)-t1(I)) ./ (f2(I).^(k-1)-f1(I).^(k-1)) .*fi(I).^(k-1);
  end
 end
end

author: HeLiangHIT


text:
function Rob = createRobots(Robot)

% CREATEROBOTS Create robots structure array.
%   Rob = CREATEROBOTS(Robot) creates the Rob() structure array to be used
%   as SLAM data. The input Robot{}  is a cell array of structures as
%   specified by the user in userData.m. There must be one Robot{} per each
%   robot considered in the simulation. See userData.m for details.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

for rob = 1:numel(Robot)

    Ri = Robot{rob}; % input robot structure

    % identification
    Ro.rob     = rob;
    Ro.id      = Ri.id;
    Ro.name    = Ri.name;
    Ro.type    = Ri.type;
    Ro.motion  = Ri.motion;

    Ro.sensors = [];

    % Robot frame in quaternion form
    ep = [Ri.position;deg2rad(Ri.orientationDegrees)];
    EP = diag([Ri.positionStd;deg2rad(Ri.orientationStd)].^2);
    [qp,QP] = propagateUncertainty(ep,EP,@epose2qpose); % frame and cov. in quaternion

    % control and rest of the state
    switch Ri.motion

        case {'constVel'}
            % control
            Ro.con.u    = [Ri.dv;deg2rad(Ri.dwDegrees)];
            Ro.con.uStd = [Ri.dvStd;deg2rad(Ri.dwStd)];
            Ro.con.U    = diag(Ro.con.uStd.^2);
            Ro.com.W    = 1/Ro.con.U;  % Information matrix

            % velocity states
            v = [Ri.velocity;deg2rad(Ri.angularVelDegrees)];
            V = diag([Ri.velStd;deg2rad(Ri.angVelStd)].^2);

            % state
            Ro.state.x    = [qp;v]; % state
            Ro.state.P    = blkdiag(QP,V);
            Ro.state.dx   = zeros(12,1);

        case {'odometry'}
            % control
            Ro.con.u    = [Ri.dx;deg2rad(Ri.daDegrees)];
            Ro.con.uStd = [Ri.dxStd;deg2rad(Ri.daStd)];
            Ro.con.U    = diag(Ro.con.uStd.^2);
            Ro.com.W    = Ro.con.U^-1;  % Information matrix

            % state
            Ro.state.x    = qp; % state
            Ro.state.P    = EP;
            Ro.state.dx   = zeros(6,1);

        otherwise
            error('Unknown motion model ''%s'' for robot %d.',Robot{rob}.motion,Robot{rob}.id);
    end

    Ro.state.r  = []; % Points  Map.x  into  state.dx
    Ro.state.size = numel(Ro.state.x);   % state size
    Ro.state.dsize = numel(Ro.state.dx);

    Ro.frame.r = [];
    Ro.frame.x = qp;
    Ro.frame.P = QP;
    Ro.frame = updateFrame(Ro.frame);

    Rob(rob) = Ro; % output robot structure

end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

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

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function mode = getmod

% this function returns the 'mode' value, which
% determines the method used for the computation of sidereal
% time and the terrestrial-to-celestial transformation, and the
% accuracy of nutation and related calculations.

%  mode = selection for method and accuracy (out)

%         mode = 0 means cio-based method, full accuracy
%         mode = 1 means cio-based method, reduced accuracy
%         mode = 2 means equinox-based method, full accuracy
%         mode = 3 means equinox-based method, reduced accuracy

global imode

mode = imode;

author: goGPS-Project


text:
function [h, cm] = sldrawmanipts(X, s, crspec)
%SLDRAWMANIPTS Draws the sample points on a manifold
%
% $ Syntax $
%   - sldrawmanipts(X)
%   - sldrawmanipts(X, s, crspec)
%   - h = sldrawmanipts(X, s, crspec)
%   - [h, cm] = sldrawmanipts(X, s, crspec)
%
% $ Arguments $
%   - X:        The sample matrix (2 x n or 3 x n)
%   - s:        The scalar indicating the marker size (default = 5)
%   - crspec:   The specification of the colors, can be in either of the
%               following forms:
%               - a string (plot symbol) indicating a basic color
%               - a 1 x n row vector as a color map
%               - a function handle to calculate the color map
%                   c = f(x, y) for 2D samples
%                   c = f(x, y, z) for 3D samples
%               default = 'b';
%   - h:        The handles to the drawn objects
%   - cm:       The color map used
%
% $ History $
%   - Created by Dahua Lin, on Sep 9, 2006
%

%% parse and verify arguments

if ndims(X) ~= 2
    error('sltoolbox:invalidarg', ...
        'The sample matrix X should be a 2D matrix');
end

[d, n] = size(X);
if d ~= 2 && d ~= 3
    error('sltoolbox:invalidarg', ...
        'The samples should be 2D or 3D');
end
if d == 2
    x = X(1, :);
    y = X(2, :);
else
    x = X(1, :);
    y = X(2, :);
    z = X(3, :);
end

if nargin < 2 || isempty(s)
    s = 5;
else
    if ~isscalar(s)
        error('sltoolbox:invalidarg', ...
            's should be a scalar');
    end
end

if nargin < 3 || isempty(crspec)
    cm = 'b';
else
    if ischar(crspec)
        cm = crspec;
    elseif isnumeric(crspec)
        if ~isequal(size(crspec), [1,n])
            error('sltoolbox:sizmismatch', ...
                'The size of the color map is not consistent with the sample number');
        end
        cm = crspec;
    elseif isa(crspec, 'function_handle')
        if d == 2
            cm = crspec(x, y);
        else
            cm = crspec(x, y, z);
        end
    else
        error('sltoolbox:invalidarg', ...
            'crspec should be a string, a row vector or a function handle');
    end
end

%% draw

if d == 2
    if nargout == 0
        scatter(x, y, s, cm);
    else
        h = scatter(x, y, s, cm);
    end
else
    if nargout == 0
        scatter3(x, y, z, s, cm);
    else
        h = scatter3(x, y, z, s, cm);
    end
end

author: lmthang


text:
function [hog_data, valid_inds, vid_id] = Read_HOG_files(users, vid_ids, hog_data_dir)

    hog_data = [];
    vid_id = {};

    feats_filled = 0;

    for i=1:numel(users)

        hog_file = [hog_data_dir, '/train/' users{i} '.hog'];
        if(~exist(hog_file, 'file'))
            hog_file = [hog_data_dir, '/devel/' users{i} '.hog'];
        end

        f = fopen(hog_file, 'r');

        num_cols = fread(f, 1, 'int32');
        if(isempty(num_cols))
            break;
        end

        num_rows = fread(f, 1, 'int32');
        num_chan = fread(f, 1, 'int32');
        num_feats = num_rows * num_cols * num_chan + 1;

        % go to the beginning
        fseek(f, 0, 'bof');

        % Read only the relevant bits

        % Skip to the right start element (1 indexed)
        fseek(f, 4*(4+num_rows*num_rows*num_chan)*(vid_ids(i,1)-1), 'bof');

        feature_vec = fread(f, [4 + num_rows * num_cols * num_chan, vid_ids(i,2) - vid_ids(i,1)], 'float32');
        fclose(f);

        curr_data = feature_vec(4:end,:)';
        curr_ind = size(curr_data,1);

        vid_id_curr = cell(size(curr_data,1),1);
        vid_id_curr(:) = users(i);

        vid_id = cat(1, vid_id, vid_id_curr);

        % Assume same number of frames per video
        if(i==1)
            hog_data = zeros(sum(vid_ids(:,2)-vid_ids(:,1)), num_feats);
        end

        hog_data(feats_filled+1:feats_filled+curr_ind,:) = curr_data;

        feats_filled = feats_filled + curr_ind;

    end
    valid_inds = hog_data(:,1) > 0;
    hog_data = hog_data(:,2:end);
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [X_poly] = polyFeatures(X, p)
%POLYFEATURES Maps X (1D vector) into the p-th power
%   [X_poly] = POLYFEATURES(X, p) takes a data matrix X (size m x 1) and
%   maps each example into its polynomial features where
%   X_poly(i, :) = [X(i) X(i).^2 X(i).^3 ...  X(i).^p];
%

% You need to return the following variables correctly.
X_poly = zeros(numel(X), p);

% ====================== YOUR CODE HERE ======================
% Instructions: Given a vector X, return a matrix X_poly where the p-th
%               column of X contains the values of X to the p-th power.

m = size(X, 1);

for i=1:m

    poly_feature = zeros(p, 1);

    for j=1:p
        poly_feature(j) =  X(i).^j;
    end

    X_poly(i, :) = poly_feature;
end

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

end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
function address = resolveip(input)
%resolveip Summary of this function goes here
%   Detailed explanation goes here
    try
        address = java.net.InetAddress.getByName(input);
    catch
        error(sprintf('Unknown host %s.', input));
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function obj = weight(obj, mode, varargin)
% Note: w must be always within the range 0 and 1! 0 <= w <= 1

% Input parsing ----------------------------------------------------------------

validMode = {'DeltaAngle' 'Roughness'};

p = inputParser;
p.addRequired('mode', @(x) any(strcmpi(x, validMode)));
p.parse(mode, varargin{:});
p = p.Results;
% Clear required inputs to avoid confusion
clear mode

% Start ------------------------------------------------------------------------

procHierarchy = {'CORRPOINTS' 'WEIGHT'};
msg('S', procHierarchy);
msg('I', procHierarchy, sprintf('Corr. points label = ''%s''', obj.label));
msg('I', procHierarchy, sprintf('IN: Mode = ''%s''', p.mode));

% Correspondences present? -----------------------------------------------------

if obj.noCP == 0
    msg('I', procHierarchy, 'no correspondences present!');
    msg('E', procHierarchy);
    return;
end

% Weight -----------------------------------------------------------------------

if strcmpi(p.mode, 'DeltaAngle')

    w = abs(cosg(obj.dAlpha));

elseif strcmpi(p.mode, 'Roughness')

    roughnessMax    = max([obj.A1.roughness  obj.A2.roughness], [], 2); % max for each row
    roughnessMaxAll = max([obj.A1.roughness; obj.A2.roughness]);        % max of all points

    w = 1 - roughnessMax./roughnessMaxAll;

end

% Update weights
obj.A.w = obj.A.w .* w;

% End --------------------------------------------------------------------------

msg('E', procHierarchy);

end

author: pglira


text:
function linplus_test589 ( )

%*****************************************************************************80
%
%% TEST589 tests R8TO_INDICATOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST589\n' );
  fprintf ( 1, '  R8TO_INDICATOR sets up a R8TO indicator matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  a = r8to_indicator ( n );

  r8to_print ( n, a, '  The R8TO indicator matrix:' );

  return
end

author: johannesgerer