text:
function blastStructure=getBlastFromExcel(models,blastFile,organismId)
% getBlastFromExcel
%   Retrieves gene homology information from Excel files. Used as
%   input to getModelFromHomology.
%
%   Input:
%   models          a cell array of model structures
%   blastFile       Excel file with homology information
%   organismId      the id of the organism of interest (as described in the
%                   Excel file)
%
%   Output:
%   blastStructure  structure containing the information in the Excel
%                   sheets.
%
%   The Excel file should contain a number of spreadsheets which in turn
%   contain the bidirectional homology measurements between the genes in the
%   organisms. The first and second column headers in each sheet is the
%   "to" and "from" model ids (as defined in models or for the new organism).
%   The entries should correspond to the gene names in those models. The third,
%   fourth, fifth, sixth and seventh columns represent the E-value, alignment
%   length, identity, bitscore and percentage of positive-scoring matches for
%   each measurement (captions should be "E-value", "Alignment length",
%   "Identity", "Bitscore" and "PPOS").
%
%   Usage: blastStructure=getBlastFromExcel(models,blastFile,organismId)

if ~isfile(blastFile)
    error('BLAST result file %s cannot be found',string(blastFile));
end
organismId=char(organismId);

blastStructure=[];

%Get a list of model IDs
organisms=cell(numel(models)+1,1);
organisms{1}=organismId;
for i=1:numel(models)
    organisms{i+1}=models{i}.id;
end

%Get all the spreadsheets in the file
[type, sheets]=xlsfinfo(blastFile);

%Check if the file is a Microsoft Excel Spreadsheet
if ~any(regexp(type,'Excel Spreadsheet'))
    EM='The file is not a Microsoft Excel Spreadsheet';
    dispEM(EM);
end

for i=1:numel(sheets)
    %Check if the sheet has the right header and deal with organisms that
    %are in "models"
    [values,dataSheet]=xlsread(blastFile,i);
    labels=dataSheet(1,:);
    if strcmpi(labels{3},'E-value') && strcmpi(labels{4},'Alignment length') ...
            && strcmpi(labels{5},'Identity') && strcmpi(labels{6},'Bitscore') ...
            && strcmpi(labels{7},'PPOS')
        %At least one of the organisms must have a model
        fromID=find(strcmpi(labels{1},organisms));
        toID=find(strcmpi(labels{2},organisms));
        %Check that the organism ids exist and that one of them is the
        %organism of interest
        if any(fromID) && any(toID) && (toID==1 || fromID==1)
            %Check that no gene ids are empty. This could for example be
            %the case if the gene names are wrongly formatted as numbers
            %instead of strings
            emptyNames=cellfun(@isempty,dataSheet(2:end,1)) | cellfun(@isempty,dataSheet(2:end,2));
            if any(emptyNames)
                if all(emptyNames)
                    EM=['Only empty gene names in sheet from ' organisms{fromID} ' to ' organisms{toID}];
                    dispEM(EM);
                else
                    EM=['Empty gene names in sheet from ' organisms{fromID} ' to ' organisms{toID} '. Ignoring genes with empty names'];
                    dispEM(EM,false);
                end
            end
            blastStructure(numel(blastStructure)+1).toId=organisms{toID};
            blastStructure(numel(blastStructure)).fromId=organisms{fromID};
            blastStructure(numel(blastStructure)).fromGenes=dataSheet(2:end,1);
            blastStructure(numel(blastStructure)).toGenes=dataSheet(2:end,2);
            blastStructure(numel(blastStructure)).evalue=values(:,1);
            blastStructure(numel(blastStructure)).aligLen=values(:,2);
            blastStructure(numel(blastStructure)).identity=values(:,3);
            blastStructure(numel(blastStructure)).bitscore=values(:,4);
            blastStructure(numel(blastStructure)).ppos=values(:,5);

            %Remove matches where any of the values is NaN. This would have
            %been done anyways in getModelFromHomology, but it's neater to
            %do it here
            I=isnan(blastStructure(end).evalue) | isnan(blastStructure(end).aligLen) | isnan(blastStructure(end).identity);
            blastStructure(end).fromGenes(I)=[];
            blastStructure(end).toGenes(I)=[];
            blastStructure(end).evalue(I)=[];
            blastStructure(end).aligLen(I)=[];
            blastStructure(end).identity(I)=[];
            blastStructure(end).bitscore(I)=[];
            blastStructure(end).ppos(I)=[];
        else
            if isempty(toID) || isempty(fromID)
                EM=['The data in sheet ' sheets{i} ' has no corresponding model. Ignoring sheet'];
                dispEM(EM,false);
            else
                EM=['The data in sheet ' sheets{i} ' does not involve the organism of interest. Ignoring sheet'];
                dispEM(EM,false);
            end
        end
    else
        EM=['The data in sheet ' sheets{i} ' is not correctly formatted. Ignoring sheet'];
        dispEM(EM,false);
    end
end

end


text:
% Mathematics Q2199546
% https://math.stackexchange.com/questions/2199546
% Minimizing the Sum of Quadratic Form with Equality Constraint
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     29/07/2017
%   *   First release.

%% General Parameters

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = OFF;

%% Simulation Parameters

numRows = 4;
numCols = 3; %<! Number of Vectors - i (K in the question)

numIterations   = 1000;
stepSize        = 0.05;

%% Generate Data

vE          = ones([numRows, 1]);
tA          = randn([numRows, numRows, numCols]);
mC          = randn([numRows, numCols]);
paramLambda = 1;

for ii = 1:numCols
    tA(:, :, ii) = tA(:, :, ii).' * tA(:, :, ii);
end

hProjEquality = @(mX) mX - ((sum(mX, 2) - vE) / numCols);

%% Solution by CVX

cvx_begin('quiet')
    cvx_precision('best');
    variable mXCvx(numRows, numCols)
    objVal = quad_form( mXCvx(:, 1), tA(:, :, 1) ) + (paramLambda * (mC(:, 1).' * mXCvx(:, 1)));
    for ii = 2:numCols
        objVal = objVal + quad_form( mXCvx(:, ii), tA(:, :, ii) ) + (paramLambda * (mC(:, ii).' * mXCvx(:, ii)));
    end
    minimize( objVal )
    subject to
        mXCvx * ones([numCols, 1]) == vE; %<! vE to be Eigen Vector for Eigen Value 1
cvx_end

disp([' ']);
disp(['CVX Solution Summary']);
disp(['The CVX Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(cvx_optval)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(mXCvx(:).'), ' ]']);
disp([' ']);

%% Solution by Projected Gradient Descent

mX = zeros([numRows, numCols]);

for ii = 1:numIterations
    for jj = 1:numCols
        mX(:, jj) = mX(:, jj) - (stepSize * ((2 * tA(:, :, jj) * mX(:, jj)) + (paramLambda * mC(:, jj))));
    end
    mX = hProjEquality(mX);
end

objVal = 0;
for ii = 1:numCols
    objVal = objVal + (mX(:, ii).' * tA(:, :, ii) * mX(:, ii)) + (paramLambda * mC(:, ii).' * mX(:, ii));
end

disp([' ']);
disp(['Projected Gradient Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(objVal)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(mX(:).'), ' ]']);
disp([' ']);

%% Restore Defaults

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

author: RoyiAvital

repo: StackExchangeCodes


text:
function [ num_int, pint ] = plane_imp_triangle_int_add_3d ( p1, p2, dist1, ...
  dist2, num_int, pint )

%*****************************************************************************80
%
%% PLANE_IMP_TRIANGLE_INT_ADD_3D is a utility for plane/triangle intersections.
%
%  Discussion:
%
%    This routine is called to consider the value of the signed distance
%    from a plane of two nodes of a triangle.  If the two values
%    have opposite signs, then there is a point of intersection between
%    them.  The routine computes this point and adds it to the list.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(3), P2(3), the coordinates of two vertices
%    of a triangle.
%
%    Input, real DIST1, DIST2, the signed distances of the
%    two vertices from a plane.
%
%    Input, integer NUM_INT, the number of intersection points.
%
%    Input, real PINT(3,NUM_INT), the intersection points.
%
%    Output, integer NUM_INT, the updated number of intersection points.
%
%    Output, real PINT(3,NUM_INT), the updated intersection points.
%
  dim_num = 3;

  if ( dist1 == 0.0 )
    num_int = num_int + 1;
    pint(1:dim_num,num_int) = p1(1:dim_num);
  elseif ( dist2 == 0.0 )
    num_int = num_int + 1;
    pint(1:dim_num,num_int) = p2(1:dim_num);
  elseif ( dist1 * dist2 < 0.0 )
    alpha = dist2 / ( dist2 - dist1 );
    num_int = num_int + 1;
    pint(1:dim_num,num_int) = alpha * p1(1:dim_num) + ( 1.0 - alpha ) * p2(1:dim_num);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
EEGDATA.buffer(:,1+mod(EEGDATA.smax:EEGDATA.smax+size(EEGDATA_chunk,2)-1,EEGDATA.buffer_len)) = EEGDATA_chunk;
EEGDATA.smax = EEGDATA.smax + size(EEGDATA_chunk,2);

repo: REST


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

% Fast Fourier Transform of the sequence x=[1 2 3],n=0,1,2

x=[1 2 3];

Xk1=fft(x)

Xk2=dft(x)

author: Sable

repo: mcbench-benchmarks


text:
function I = vl_test_pattern(n)
% VL_TEST_PATTERN  Generate test pattern
%   I=VL_TEST_PATTERN(N) returns the N-th test pattern.

ur    = linspace(-1,1,128) ;
vr    = linspace(-1,1,128) ;
[u,v] = meshgrid(ur,vr);

switch n
  case 1
    %I     = u.^2 + v.^2 > (1/4).^2 ;
    I     = abs(u) + abs(v) > (1/4) ;
    I     = 255 * I ;
    I(1:64,:) = 0 ;

  case 2
    I = zeros(100,100) ;
    I(20:100-20+1,20:100-20+1) = 128 ;
    I(30:100-30+1,30:100-30+1) = 200 ;
    I(50,50)                   = 255 ;
    I(50,55)                   = 250 ;
    I(50,45)                   = 245 ;
    I = 255 - I ;

  case 3
    I = 255 * vl_imsmooth(checkerboard(10,10),1) ;

  case 4
    I = 255 * rand(32,32) ;

  case 101
    I = 255 * vl_imreadbw(fullfile(vlfeat_root,'data','a.jpg')) ;

  case 102
	 I = 255 * vl_imreadbw(fullfile(vlfeat_root,'data','box.pgm')) ;

 case 'cone'
	I = sqrt(u.^2+v.^2) ;

end

author: jianxiongxiao

repo: ProfXkit


text:
%% Execute this file to collect validation data for model identification

x0 = x(end,:);
tspanV =[tspan(end):dt:20];

options = odeset('RelTol',1e-10,'AbsTol',1e-10*ones(1,n));

switch InputSignalType
    case 'sine2'
        A = 2;
        forcing = @(x,t) [(A*(sin(1*t)+sin(.1*t))).^2];
        [t_valid,x_valid]=ode45(@(t,x) LorenzSys(t,x,forcing(x,t),p),tspanV,x0,options);
        u_valid = forcing(0,tspanV);

    case 'chirp'
        A = 10;
        forcing = @(x,t) A*chirp(t,[],max(tspanV),1.).^2;
        [t_valid,x_valid]=ode45(@(t,x) LorenzSys(t,x,forcing(x,t),p),tspanV,x0,options);
        u_valid = forcing(0,tspanV);

    case 'noise'
        vareps = 0.01;
        Diff = @(t,x) [0; vareps];
        SDE = sde(@(t,x) LorenzSys(t,x,forcing(x,t),p),Diff,'StartState',x0');
        rng(1,'twister')
        [x_valid, t_valid, u_valid] = simByEuler(SDE, length(tspanV), 'DeltaTime', dt);
        u_valid = u_valid';
        x_valid = x_valid(1:end-1,:); t_valid = t_valid(1:end-1);
    case 'prbs'
        A = 1;
        taulim = [0.1 5];
        states = [-1 1];
        Nswitch = 300;
        forcing = @(x,t) A*prbs(taulim, Nswitch, states, t,0);

        [t_valid,x_valid]=ode45(@(t,x) LorenzSys(t,x,forcing(x,t),p),tspanV,x0,options);

        u_valid = zeros(size(tspanV));
        for i = 1:length(tspanV)
            u_valid(i) = forcing(0,tspanV(i));
        end
        figure,plot(tspanV,u_valid)

    case 'sphs'
        Pf = 1; % Fundamental period
        K = 8;
        A = 10;
        forcing = @(x,t) A*sphs(Pf,K,t);
        [t_valid,x_valid]=ode45(@(t,x) LorenzSys(t,x,forcing(x,t),p),tspanV,x0,options);
        u_valid = forcing(0,tspanV);

    case 'mixed'
        A = [4 5 4 3];
        Pf = 8; % Fundamental period
        K = 16;
        taulim = [0.1 3];
        states = [-1 1];
        Nswitch = 40;

        tspan1 = [0:dt:100];
        forcing = @(x,t) SI_Input(t,A,Pf,K,tspan1,taulim, Nswitch, states);
        [t_valid,x_valid] = ode45(@(t,x) LorenzSys(t,x,forcing(x,t),p),tspan1,x0,options);

        u_valid = zeros(size(tspan1));
        for i = 1:length(tspan1)
            u_valid(i) = forcing(0,tspan1(i));
        end
        figure,plot(tspan1,u_valid)

        A = 2;
        tspanv = [100:dt:200];
        forcingV = @(x,t) [(A*(sin(1*t)+sin(.1*t))).^2];
        [tv,xv]=ode45(@(t,x) LorenzSys(t,x,forcing(x,t),p),tspanv,x_valid(end,:),options);
        uv = forcingV(0,tspanv);

        u_valid = [u_valid, uv(2:end)];
        x_valid = [x_valid; xv(2:end,:)];
        t_valid = tspanV';
end

T = length(tspanV);
N = length(tspanV);

%% Show data
figure;
plot(t_valid,x_valid,'LineWidth',1.5)
xlabel('Time')
ylabel('xi')
legend('x','y','z')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')

author: eurika-kaiser

repo: SINDY-MPC


text:
function cc_io_test02 ( )

%*****************************************************************************80
%
%% CC_IO_TEST02 tests CC_HEADER_READ and CC_DATA_READ.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 July 2014
%
%  Author:
%
%    John Burkardt
%
  prefix = 'simple';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_IO_TEST02\n' );
  fprintf ( 1, '  Read a sparse matrix in CC format from 3 files.\n' );
%
%  Read the header.
%
  [ ncc, n ] = cc_header_read ( prefix );
%
%  Read the matrix data.
%
  [ icc, ccc, acc ] = cc_data_read ( prefix, ncc, n );
%
%  Print the CC matrix.
%
  m = n;
  cc_print ( m, n, ncc, icc, ccc, acc, '  The matrix in 1-based CC format:' );

  return
end


text:
function out = FDSInit(startPos, endPos, map, scalling)

%% generate shape pattern
    radius = 10/scalling;
    mat = zeros(radius*2,radius*2);
    dista = @(a,b) sqrt(a*a+b*b);

    shapePattern = [];

    for x = 1:1:2*radius
        for y =1:1:2*radius
            if dista(x-radius-0.5,y-radius-0.5) > radius -1 && dista(x-radius-0.5,y-radius-0.5) < radius
                mat(x,y) = 1;
                shapePattern = [shapePattern; floor(x-radius) floor(y-radius) 0 0];
            end
        end
    end
    neighbours = [
        0 1 0 0;
        -1 0 0 0;
        0 -1 0 0;
        1 1 0 0;
        -1 -1 0 0;
        1 0 0 0;
        1 -1 0 0;
        -1 1 0 0;
    ]';

%% prepare all data

	out.map = map;
    out.startPos = [startPos(2:-1:1); 0; 0];
    out.endPos = [endPos(2:-1:1); 0; 0];
    startPos = out.startPos;

    out.sacalling = scalling;
    out.pattern = shapePattern';
    out.ucc = neighbours;
    out.height = ceil(length(out.map(:,1)));
    out.width = ceil(length(out.map(1,:)));
    out.graph = zeros(out.height, out.width,7);
    out.graph(:,:,1:2) = inf;
    out.graph(:,:,3) = -1;
    out.graph(:,:,5) = -1;
    out.graph(:,:,6) = 0;
    out.graph(:,:,7) = 0;
    out.kM = 0;
    SQRT2 = sqrt(2)-1;
    out.comparator = DStarcmp;
    out.stack = java.util.PriorityQueue(180247, out.comparator);

    setk(out.endPos, 0);
    seth(out.endPos, 0);
    setQ(out.endPos);
    out.endPos(3:4) = [heur(out.endPos); 0];
    add(out.stack, out.endPos);

    %-----------------------------------------------------------
    function setQ(s)
         out.graph(s(1), s(2), 3) = 1;
    end
    function seth(s, val)
        out.graph(s(1), s(2),1) = val;
    end
    function setk(s, val)
        out.graph(s(1), s(2),2) = val;
    end
    function out = heur(s)
        k = abs(startPos - s);
        out = SQRT2*min(k(1:2)) + max(k(1:2));
    end

end

author: LazyFalcon


text:
function cx = legendre_associated_normalized ( n, m, x )

%*****************************************************************************80
%
%% LEGENDRE_ASSOCIATED: normalized associated Legendre functions.
%
%  Discussion:
%
%    The unnormalized associated Legendre functions P_N^M(X) have
%    the property that
%
%      Integral ( -1 <= X <= 1 ) ( P_N^M(X) )^2 dX
%      = 2 * ( N + M )! / ( ( 2 * N + 1 ) * ( N - M )! )
%
%    By dividing the function by the square root of this term,
%    the normalized associated Legendre functions have norm 1.
%
%    However, we plan to use these functions to build spherical
%    harmonics, so we use a slightly different normalization factor of
%
%      sqrt ( ( ( 2 * N + 1 ) * ( N - M )! ) / ( 4 * pi * ( N + M )! ) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%  Parameters:
%
%    Input, integer N, the maximum first index of the Legendre
%    function, which must be at least 0.
%
%    Input, integer M, the second index of the Legendre function,
%    which must be at least 0, and no greater than N.
%
%    Input, real X, the point at which the function is to be
%    evaluated.  X must satisfy -1 <= X <= 1.
%
%    Output, real CX(1:N+1), the values of the first N+1 function.
%
  if ( m < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEGENDRE_ASSOCIATED_NORMALIZED - Fatal error!\n' );
    fprintf ( 1, '  Input value of M is %d\n', m );
    fprintf ( 1, '  but M must be nonnegative.\n' );
    error ( 'LEGENDRE_ASSOCIATED_NORMALIZED - Fatal error!' );
  end

  if ( n < m )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEGENDRE_ASSOCIATED_NORMALIZED - Fatal error!\n' );
    fprintf ( 1, '  Input value of M = %d\n', m );
    fprintf ( 1, '  Input value of N = %d\n', n );
    fprintf ( 1, '  but M must be less than or equal to N.\n' );
    error ( 'LEGENDRE_ASSOCIATED_NORMALIZED - Fatal error!' );
  end

  if ( x < -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEGENDRE_ASSOCIATED_NORMALIZED - Fatal error!\n' );
    fprintf ( 1, '  Input value of X = %f\n', x );
    fprintf ( 1, '  but X must be no less than -1.\n' );
    error ( 'LEGENDRE_ASSOCIATED_NORMALIZED - Fatal error!' );
  end

  if ( 1.0 < x )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEGENDRE_ASSOCIATED_NORMALIZED - Fatal error!\n' );
    fprintf ( 1, '  Input value of X = %f\n', x );
    fprintf ( 1, '  but X must be no more than 1.\n' );
    error ( 'LEGENDRE_ASSOCIATED_NORMALIZED - Fatal error!' );
  end

  cx(1:m) = 0.0;

  cx(m+1) = 1.0;
  somx2 = sqrt ( 1.0 - x * x );

  fact = 1.0;
  for i = 1 : m
    cx(m+1) = -cx(m+1) * fact * somx2;
    fact = fact + 2.0;
  end

  if ( m < n )
    cx(m+2) = x * ( 2 * m + 1 ) * cx(m+1);
  end

  for i = m+2 : n
    cx(i+1) = ( ( 2 * i     - 1 ) * x * cx(i) ...
              + (   - i - m + 1 ) *     cx(i-1) ) ...
              / (     i - m     );
  end
%
%  Normalization.
%
  for mm = m : n
    factor = sqrt ( ( ( 2 * mm + 1 ) * r8_factorial ( mm - m ) ) ...
      / ( 4.0 * pi * r8_factorial ( mm + m ) ) );
    cx(mm+1) = cx(mm+1) * factor;
  end

  return
end


text:
function tests = test_ft_preproc_rereference

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preproc_rereference

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  fn = localfunctions;
  for i=1:numel(fn)
    feval(fn{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testOptions(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nchan   = 8;
nsample = 1000;
dat     = randn(nchan, nsample) + 1;

handlenan = false;
leadfield = randn(8, 1000);

result = {};
result{end+1} = ft_preproc_rereference(dat, 1,   'avg', handlenan);
result{end+1} = ft_preproc_rereference(dat, 1:4, 'avg', handlenan);
result{end+1} = ft_preproc_rereference(dat, 1:8, 'avg', handlenan);
result{end+1} = ft_preproc_rereference(dat, 1:4, 'median', handlenan);
result{end+1} = ft_preproc_rereference(dat, 1:8, 'median', handlenan);
result{end+1} = ft_preproc_rereference(dat, 1,   'rest', handlenan, leadfield);
result{end+1} = ft_preproc_rereference(dat, 1:4, 'rest', handlenan, leadfield);
result{end+1} = ft_preproc_rereference(dat, 1:8, 'rest', handlenan, leadfield);

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequal(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function ok = validate(key)
% Checks that an object of type Key obeys certain rules that
% must always be true.
% Inputs:
%   key: the object to be checked.
% Outputs:
%   ok: a boolean value indicating whether the object is valid.

assert(nargin==1)
assert(isa(key,'Key'))

ok = iscell(key.modelset);
ok = ok && iscell(key.segset);

nummods = length(key.modelset);
numsegs = length(key.segset);

ok = ok && (size(key.tar,1)==nummods);
ok = ok && (size(key.tar,2)==numsegs);

ok = ok && (size(key.non,1)==nummods);
ok = ok && (size(key.non,2)==numsegs);

author: nesl


text:
function D = spm_eeg_average(S)
% Average each channel over trials or trial types
% FORMAT D = spm_eeg_average(S)
%
% S        - optional input struct
%    fields of S:
% D        - MEEG object or filename of M/EEG mat-file with epoched data
% S.robust      - (optional) - use robust averaging
%                 .savew  - save the weights in an additional dataset
%                 .bycondition - compute the weights by condition (1,
%                                default) or from all trials (0)
%                 .ks     - offset of the weighting function (default: 3)
% S.prefix     - prefix for the output file (default - 'm')
%
% Output:
% D        - MEEG object (also written on disk)
%__________________________________________________________________________
% Copyright (C) 2008-2017 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: spm_eeg_average.m 7125 2017-06-23 09:49:29Z guillaume $

SVNrev = '$Rev: 7125 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','M/EEG averaging'); spm('Pointer','Watch');

%-Get input parameters
%--------------------------------------------------------------------------
if ~isfield(S, 'prefix'),       S.prefix = 'm';           end
if ~isfield(S, 'robust'),       S.robust = 0;             end

%-Configure robust averaging
%--------------------------------------------------------------------------
if isstruct(S.robust)
    if ~isfield(S.robust, 'savew'),        S.robust.savew =  0;        end
    if ~isfield(S.robust, 'bycondition'),  S.robust.bycondition = 0;   end
    if ~isfield(S.robust, 'ks'),           S.robust.ks =  3;           end
    if ~isfield(S.robust, 'removebad'),    S.robust.removebad =  0;    end

    robust      = 1;
    savew       = S.robust.savew;
    bycondition = S.robust.bycondition;
    ks          = S.robust.ks;
    removebad   = S.robust.removebad;
else
    robust = 0;
end

%-Get MEEG object
%--------------------------------------------------------------------------
D = spm_eeg_load(S.D);

montage_ind = D.montage('getindex');
D           = D.montage('switch', 0);

%-Check that there is any good data available
%--------------------------------------------------------------------------
if ntrials(D)==0 || isempty(indtrial(D, D.condlist, 'GOOD'))
    warning('No good trials for averaging were found. Nothing to do.');
    return;
end

%-Redirect to Time-Frequency averaging if necessary
%--------------------------------------------------------------------------
if strncmpi(D.transformtype,'TF',2) % TF and TFphase
    D = spm_eeg_average_TF(S);
    return
end

%-Generate new MEEG object with new files
%--------------------------------------------------------------------------
Dnew = clone(D, [S.prefix fname(D)], [D.nchannels D.nsamples D.nconditions]);
Dnew = type(Dnew, 'evoked');

if robust && savew
    Dw = clone(D, ['W' fname(D)]);
end

%-Do the averaging
%--------------------------------------------------------------------------
cl   = D.condlist;

ni = zeros(1, D.nconditions);
for i = 1:D.nconditions
    w = indtrial(D, deblank(cl{i}), 'GOOD');
    ni(i) = length(w);
    if ni(i) == 0
        warning('%s: No trials for trial type %d', D.fname, cl{i});
    end
end

goodtrials  = indtrial(D, cl, 'GOOD');
chanind     = D.indchantype({'MEEG', 'LFP'});

if prod(size(D))*8 < spm('memory')

    %-Average everything at once if there is enough memory
    %======================================================================
    spm_progress_bar('Init', D.nconditions, 'Averages done');
    if D.nconditions > 100, Ibar = floor(linspace(1, D.nconditions, 100));
    else Ibar = 1:D.nconditions; end

    if robust && ~bycondition
        W1      = ones(D.nchannels, D.nsamples, length(goodtrials));
        Y       = D(chanind, :, goodtrials);
        if removebad
            bad    = badsamples(D, chanind, ':', goodtrials);
            Y(bad) = NaN;
        end
        [Y, W2] = spm_robust_average(Y, 3, ks);
        W1(chanind, :, :) = W2;
        if savew
            Dw(:, :, goodtrials) = W1;
        end
        W = zeros(size(D));
        W(:, :, goodtrials) = W1;
    end

    for i = 1:D.nconditions

        w = indtrial(D, deblank(cl{i}), 'GOOD');

        if isempty(w)
            continue;
        end

        if ~robust
            Dnew(:, :, i) = mean(D(:, :, w), 3);
        else
            if bycondition
                W       = ones(D.nchannels, D.nsamples, length(w));
                Y       = D(chanind, :, w);
                if removebad
                    bad     = badsamples(D, chanind, ':', w);
                    Y(bad)  = NaN;
                end
                [Y, W1] = spm_robust_average(Y, 3, ks);
                W(chanind, :, :)    = W1;
                Dnew(chanind, :, i) = Y;
                if length(chanind)<D.nchannels
                    Dnew(setdiff(1:D.nchannels, chanind), :, i) = mean(D(setdiff(1:D.nchannels, chanind), :, w), 3);
                end
                if savew
                    Dw(:, :, w)   = W;
                end
            else
                X = D(:, :, w);
                X(isnan(X))      = 0;
                Dnew(:, :, i) = ...
                    sum(W(:, :, w).*X, 3)./sum(W(:, :, w), 3);
            end
        end

        if ismember(i, Ibar), spm_progress_bar('Set', i); end

    end  % for i = 1:D.nconditions
else
    %-Averaging by channel
    %======================================================================
    spm_progress_bar('Init', D.nchannels, 'Channels completed');
    if D.nchannels > 100, Ibar = floor(linspace(1, D.nchannels, 100));
    else Ibar = [1:D.nchannels]; end
    for j = 1:D.nchannels
        if robust && ~bycondition
            if ismember(j, chanind)
                Y       = D(j, :, goodtrials);
                if removebad
                    bad     = badsamples(D, j, ':', goodtrials);
                    Y(bad)  = NaN;
                end
                [Y, W1] = spm_robust_average(Y, 3, ks);
                W = zeros([1 D.nsamples D.ntrials]);
                W(1, :, goodtrials) = W1;
            else
                W1 = ones(1, D.nsamples, length(goodtrials));
            end

            if savew
                Dw(j, :, goodtrials) = W1;
            end
        end

        for i = 1:D.nconditions

            w = indtrial(D, deblank(cl{i}), 'GOOD');

            if isempty(w)
                continue;
            end

            if ~robust || ~ismember(j, chanind)
                Dnew(j, :, i) = mean(D(j, :, w), 3);
            else
                if bycondition
                    Y       = D(j, :, w);
                    if removebad
                        bad     = badsamples(D, j, ':', w);
                        Y(bad)  = NaN;
                    end
                    [Y, W] = spm_robust_average(Y, 3, ks);
                    Dnew(j, :, i) = Y;
                    if savew
                        Dw(j, :, w)   = W;
                    end
                else
                    X = D(j, :, w);
                    X(isnan(X))      = 0;
                    Dnew(j, :, i) = ...
                        sum(W(1, :, w).*X, 3)./sum(W(1, :, w), 3);
                end
            end

        end  % for i = 1:D.nconditions
        if ismember(j, Ibar), spm_progress_bar('Set', j); end
    end
end

spm_progress_bar('Clear');

%-Update some header information
%--------------------------------------------------------------------------
Dnew = conditions(Dnew, ':', cl);
Dnew = repl(Dnew, ':', ni);
Dnew = montage(Dnew, 'switch', montage_ind);

%-Display averaging statistics
%--------------------------------------------------------------------------
fprintf('%s: Number of replications per contrast:\n', Dnew.fname);      %-#
for i = 1:D.nconditions
    fprintf('  average %s: %d trials\n', cl{i}, ni(i));                 %-#
end

if robust
    disp('Robust averaging might have introduced high frequencies in the data. It is advised to re-apply low-pass filter.');
end

%-Save new evoked M/EEG dataset
%--------------------------------------------------------------------------
Dnew = Dnew.history(mfilename, S);
save(Dnew);

if robust && savew
    Dw = Dw.history(mfilename, S);
    save(Dw);
end

D = Dnew;

%-Cleanup
%--------------------------------------------------------------------------
fprintf('%-40s: %30s\n','Completed',spm('time'));                       %-#
spm('FigName','M/EEG averaging: done'); spm('Pointer','Arrow');

author: spm


text:
% Generate analysis time instants, optionaly using a given f0 curve
%
% Octave compatible
%
% Inputs
%  wav            : Wavform vector
%  fs             : Sampling frequency
%  tmargin        : Margin at start and end
%  usesampletimes : Time instants are thoses of samples.
%                   Time instants can't be between two sample times,
%                   such as window centers correspond to samples times
%  method
%     0: Regular instants.
%        varargin{1}: step size [s]
%     1: Keep f0sin time instants and fill with interpolated values.
%        varargin{1}: An f0 curve
%        varargin{2}: number of instants per period
%     2: Generate new instants from the given f0 curve.
%        varargin{1}: An f0 curve
%        varargin{2}: number of instants per period
%        Skip nan, inf and zero values
%
% Outputs
%  f0sout : f0 feature with the generated times
%
% Copyright (c) 2011 University of Crete - Computer Science Department (UOC-CSD)
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function at = gen_analysis_times(wav, fs, tmargin, usesampletimes, method, varargin)

    if nargin<5;    usesampletimes = true; end
    if nargin<6;    method = 3; end

    wavdur = (length(wav)-1)/fs;

    if method==0
        step = varargin{1};
        disp('    Generate regular analysis time');
        if usesampletimes
            at = round(tmargin*fs)/fs:round(step*fs)/fs:round((wavdur-tmargin)*fs)/fs;
        else
            at = tmargin:step:(wavdur-tmargin);
        end
        at = at';

    elseif method==1
        f0sin = varargin{1};
        nbperper = varargin{2};
        if nbperper>1
            disp(['    Add extra ' num2str(nbperper) ' analysis instants between each f0 value']);
            f0sin = f0sin(find(~isnan(f0sin(:,2)) & ~isinf(f0sin(:,2)) & f0sin(:,2)~=0),:);
            at = [];
            for n=1:size(f0sin,1)-1
                for m=0:nbperper-1
                    t = f0sin(n,1) + (m/nbperper)*(f0sin(n+1,1)-f0sin(n,1));
                    at = [at; t];
                end
            end
        else
            at = f0sin(:,1);
        end

    elseif method==2
        f0sin = varargin{1};
        nbperper = varargin{2};
        disp(['    Adapt analysis instants to f0 curve with ' num2str(nbperper) ' analysis intants per period']);
        f0sin = f0sin(find(~isnan(f0sin(:,2)) & ~isinf(f0sin(:,2)) & f0sin(:,2)~=0),:);
        at = [];
        nt = tmargin;
        while nt < wavdur - tmargin
            at = [at; nt];
            nt = nt + (1/nbperper)/interp1td(f0sin, nt);
        end

    end
    idx = find(at>tmargin & at<wavdur-tmargin);
    at = at(idx);

    if usesampletimes; at = round(fs*at)/fs; end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
% Flush out the MATLAB.
close all;
clc;
clear all;

% Read the desired image file.
ImageData=imread('auto.pnm');

% Display the original image.
figure,imshow(ImageData);
title(' Original Image: ');

% Take the input of the sigma from user.
SigmaValue = input('Enter the Sigma = ');

% Gaussian Filter implementation.
OutputImage = imgaussfilt(ImageData,SigmaValue);

% Display the output image.
figure, imshow(OutputImage);
title(' Final Image: ');

% Add noise to oriiginal image.
ImageDataNoise = imnoise(ImageData,'Salt & Pepper', 0.04);

% Display the original image with noise.
figure,imshow(ImageDataNoise);
title(' Original Image with noise: ');

% Gaussian Filter implementation.
OutputImageWithNoise = imgaussfilt(ImageDataNoise,SigmaValue);

% Display the output image.
figure,imshow(OutputImageWithNoise);
title(' Final Image(with noise): ');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


% initialize the MH by minimizing the - logposterior kernel
T = length(y);
vSwitch = ones(2, 1);
% % The parameters might be bounded, but the minimization routine will not
% % care, thus we need to make sure that the routine can go anywhere and
% % that we are still able to stay within the bounds.
% x0 = boundsINV(guess);
x0 = guess;
% posterior maximization
[fh, xh, gh, H, itct, fcount, retcodeh] = csminwel('logpostkernel',x0,.01*eye(length(x0)),[],10e-5,1000,centered,y, vSwitch, iPriorScale);
% processing the output of the maximization
postmode = xh;
JJ       = jacob(xh);
HH       = JJ * H * JJ';


text:
function adj = tet_mesh_order4_adj_set ( node_num, tetra_num, tetra_node, ...
  adj_num, adj_row )

%*****************************************************************************80
%
%% TET_MESH_ORDER4_ADJ_SET sets the nodal adjacency matrix.
%
%  Discussion:
%
%    A compressed format is used for the nodal adjacency matrix.
%
%    It is assumed that we know ADJ_NUM, the number of adjacency entries
%    and the ADJ_ROW array, which keeps track of the list of slots
%    in ADJ where we can store adjacency information for each row.
%
%    We essentially repeat the work of TET_MESH_ORDER4_ADJ_COUNT, but
%    now we have a place to store the adjacency information.
%
%    A copy of the ADJ_ROW array is useful, as we can use it to keep track
%    of the next available entry in ADJ for adjacencies associated with
%    a given row.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer TETRA_NUM, the number of tetrahedrons.
%
%    Input, integer TETRA_NODE(4,TETRA_NUM), the nodes that make up the
%    tetrahedrons.
%
%    Input, integer ADJ_NUM, the total number of adjacency relationships,
%
%    Input, integer ADJ_ROW(NODE_NUM+1), the ADJ pointer array.
%
%    Output, integer ADJ(ADJ_NUM), the adjacency information.
%

%
%  Each order 4 tetrahedron defines 6 adjacency pairs.
%
  pair(1,            1:  tetra_num) = tetra_node(1,1:tetra_num);
  pair(2,            1:  tetra_num) = tetra_node(2,1:tetra_num);

  pair(1,  tetra_num+1:2*tetra_num) = tetra_node(1,1:tetra_num);
  pair(2,  tetra_num+1:2*tetra_num) = tetra_node(3,1:tetra_num);

  pair(1,2*tetra_num+1:3*tetra_num) = tetra_node(1,1:tetra_num);
  pair(2,2*tetra_num+1:3*tetra_num) = tetra_node(4,1:tetra_num);

  pair(1,3*tetra_num+1:4*tetra_num) = tetra_node(2,1:tetra_num);
  pair(2,3*tetra_num+1:4*tetra_num) = tetra_node(3,1:tetra_num);

  pair(1,4*tetra_num+1:5*tetra_num) = tetra_node(2,1:tetra_num);
  pair(2,4*tetra_num+1:5*tetra_num) = tetra_node(4,1:tetra_num);

  pair(1,5*tetra_num+1:6*tetra_num) = tetra_node(3,1:tetra_num);
  pair(2,5*tetra_num+1:6*tetra_num) = tetra_node(4,1:tetra_num);

  pair_num = 6 * tetra_num;
%
%  Force the nodes of each pair to be listed in ascending order.
%
  pair = i4col_sort2_a ( 2, pair_num, pair );
%
%  Rearrange the columns in ascending order.
%
  pair = i4col_sort_a ( 2, pair_num, pair );
%
%  Mark all entries of ADJ so we will know later if we missed one.
%
  adj(1:adj_num) = -1;
%
%  Copy the ADJ_ROW array and use it to keep track of the next
%  free entry for each row.
%
  adj_row_copy(1:node_num) = adj_row(1:node_num);
%
%  Now set up the ADJ_ROW counts.
%
  for k = 1 : pair_num

    if ( 1 < k )
      if ( pair(1,k-1) == pair(1,k) && pair(2,k-1) == pair(2,k) )
        continue
      end
    end

    i = pair(1,k);
    j = pair(2,k);

    adj(adj_row_copy(i)) = j;
    adj_row_copy(i) = adj_row_copy(i) + 1;
    adj(adj_row_copy(j)) = i;
    adj_row_copy(j) = adj_row_copy(j) + 1;

  end

  return
end


text:
% OZDEMO
% This program creates the Ornstein-Zernike example in Chapter 3.
% [H,C]=OZDEMO returns the solution on a grid with a mesh
% spacing of 1/256.
%
function [h,c]=ozdemo
global L U rho
n=257;
epsilon=.1; sigma=2.0; rho=.2; beta=10; L=9;
dx=L/(n-1); r=0:dx:L; r=r';
%
% Compute the potential and store it in a global variable.
%
U=elj(r,sigma,epsilon,beta);
%
tol=[1.d-8,1.d-8];
x=zeros(2*n,1);
parms=[40,80,-.1];
[sol, it_hist, ierr] = nsoli(x,'oz',tol);
%
% Unpack h and c.
%
h=sol(1:n); c=sol(n+1:2*n);
%
% Plot the solution.
%
figure(1)
subplot(1,2,1)
plot(r,h,'-');
ylabel('h','Rotation',1);
xlabel('r');
subplot(1,2,2)
plot(r,c,'-');
ylabel('c','Rotation',1);
xlabel('r');
%
% Do a second solve with constant forcing terms.
% Are you getting the same results each time?
%
fb=it_hist(1,1);
parms=[40,80,-.1]; x=zeros(2*n,1);
[sola, it_hist1, ierr] = nsoli(x,'oz',tol,parms);
norm(sola-sol)
%
% plot residual vs iteration counter
%
figure(2)
ni=length(it_hist(:,1));
n1=length(it_hist1(:,1));
semilogy(0:ni-1,it_hist(:,1)/fb,'-',...
0:n1-1,it_hist1(:,1)/fb,'--');
legend('default','.1');
xlabel('Nonlinear iterations');
ylabel('Relative residual');
%
% plot residual vs function counter
%
figure(3)
semilogy(it_hist(:,2),it_hist(:,1)/fb,'-',...
it_hist1(:,2),it_hist1(:,1)/fb,'--');
xlabel('Function evaluations');
ylabel('Relative residual');
legend('default','.1');

%
function u=elj(r,sigma,epsilon,beta)
n2=length(r);
ra=r(2:n2);
r12=(sigma./ra).^12; r6=(sigma./ra).^6;
ua=exp(-4*beta*epsilon*(r12-r6));
u=[0,ua']';

repo: Aero-matlab


text:
function y = var( varargin )

%STD    Internal cvx version.

y = square_pos( std( varargin{:} ) );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


function [BVAR] = bvar(y,lags,options)

if nargin<3
    BVAR = bvar_(y,lags);
else
    BVAR = bvar_(y,lags,options);
end

end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 'bvar' generates draws from the paramters of VAR model with error
% % distributed as a multivariate normal
%
% % Core Inputs:
% % - y, data columns variables
% % - lags, lag order of the VAR
%
% % Additonal Inputs collected options:
% % - options are not mandatory. if nothing is specified then a Jeffrey prior is
% %   assumed
% % - options.K, number of draws from the posterior distribution
% % - options.hor, horizon to compute the impulse response
% % - options.fhor, horizon to compute the out-of-sample forecast
% % - options.priors, a string with he priors for the autoregressive paramters and
% %   for the scaling matrix.
% % (...) see below
%
% % Output: Draws from the conditional distribution of Phi, Sigma and
% % Omega, impulse response with the cholesky decompotion and long run
% % restricions, forecast and marginal likelihood.
%
% % Filippo Ferroni, 6/1/2015
% % Revised, 2/15/2017
% % Revised, 3/21/2018
% % Revised, 9/11/2019
% % Revised, 27/02/2020
% % Revised, 27/04/2020
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% if nargin < 2
%     error('the BVAR toolbox needs at least two inputs: data and number of lags');
% end
% if lags < 1
%     error('lags cannot be zero or negative');
% end
% % number of observable variables
% ny                  = size(y, 2);
%
% %********************************************************
% %* DEFAULT SETTINGS
% %********************************************************
% % Control random number generator
% if isOctave == 0
%     isMatlab = 1;
%     rng('default');
%     rng(999);
% else
%     isMatlab = 0;
%     % pkg load optim
%     randn('state',999);
%     rand('state',999);
% end
%
%
% % Default Settings (they can all be changed in 'options' see below)
% K                   = 5000;         % number of draws from the posterior
% hor                 = 24;           % horizon for the IRF
% fhor                = 12;           % horizon for the forecasts
% firstobs            = lags+1;       % first observation
% presample           = 0;            % using a presample for setting the hyper-parameter of the Minnesosta prior
% noconstant          = 0;            % when 0, includes a constatn in the VAR
% timetrend           = 0;            % when 1, includes a time trend in the VAR
% minn_prior_tau      = 3;            % Minnesota prior Hyper-Param: Overall Tightness
% minn_prior_decay    = 0.5;          % Minnesota prior Hyper-Param: Tighness on lags>1
% minn_prior_lambda   = 5;            % Minnesota prior Hyper-Param: Sum-of-Coefficient
% minn_prior_mu       = 2;            % Minnesota prior Hyper-Param: Co-Persistence
% minn_prior_omega    = 2;            % Minnesota prior Hyper-Param: Shocks Variance
% long_run_irf        = 0;            % when 0, it does not compute long run IRF
% irf_1STD            = 1;            % when 1, IRF are computed as 1SD increase. Else, IRF are compued as unitary increase in the shock
% cfrcst_yes          = 0;            % no conditional forecast unless defined in options
% non_explosive_      = 0;            %
% heterosked          = 0;
%
% signs_irf           = 0;
% narrative_signs_irf = 0;
% zeros_signs_irf     = 0;
% proxy_irf           = 0;
% noprint             = 0;
% nexogenous          = 0;
% exogenous           = [];
%
% % for mixed frequecy / irregurerly sampled data.
% % Interpolate the missing values of each times series.
% if any(any(isnan(y)))==1
%     warning('Activating the Mixed Frequency BVAR')
%
%     mixed_freq_on = 1;
%     index_nan_var = find(sum(isnan(y),1) ~= 0);
%     yoriginal     = y;
%     % interporlate the variables that have nans
%     T = 1:1:length(y);
%     for kk  = 1 : length(index_nan_var)
%         v                       = y(isfinite(y(:,index_nan_var(kk))),index_nan_var(kk));
%         x                       = find(isfinite(y(:,index_nan_var(kk))));
%         if isMatlab == 1
%             y(:,index_nan_var(kk))  = interp1(x,v,T','spline');
%         else
%             y(:,index_nan_var(kk))  = spline(x,v,T');
%         end
%         yinterpol               = y;
%     end
% else
%     mixed_freq_on = 0;
% end
% if mixed_freq_on == 1 && nargin < 3
%     warning(['You did not specified the aggregation of the mixed freq. '...
%         'I will treat them as stocks.']);
%     index = zeros(size(y,2),1);
% end
%
% % Priors declaration: default Jeffrey prior
% dummy    = 0;
% flat     = 1;
% priors   = priors_( );
%
% % declaring the names for the observable variables
% for v = 1 : ny
%     eval(['varnames{'   num2str(v) '} =  ''Var' num2str(v) ''';'])
% end
%
%
% %********************************************************
% %* CUSTOMIZED SETTINGS
% %********************************************************
% if nargin > 2
%     if isfield(options,'vnames')==1
%         varnames = options.vnames;
%     end
%     %======================================================================
%     % Inference options
%     %======================================================================
%     if isfield(options,'K')==1
%         K = options.K;
%     end
%     if isfield(options,'firstobs')==1
%         firstobs = options.firstobs;
%         if firstobs < lags + 1
%             error('firstobs need to be larger than lags +1')
%         end
%     end
%     if isfield(options,'presample')==1
%         presample = options.presample;
%     end
%     if isfield(options,'noconstant')==1
%         noconstant = options.noconstant;
%     end
%     if isfield(options,'timetrend')==1
%         timetrend  = options.timetrend;
%         noconstant = 0;
%     end
%     if isfield(options,'non_explosive_')==1
%         non_explosive_  = options.non_explosive_;
%     end
%     if isfield(options,'heterosked_weights')==1
%         ww  = options.heterosked_weights;
%         heterosked = 1;
%     end
%     %======================================================================
%     % Exogenous Variables options
%     %======================================================================
%     if isfield(options,'exogenous')==1 || isfield(options,'controls')==1
%         if isfield(options,'controls')==1
%             exogenous = options.controls;
%         else
%             exogenous = options.exogenous;
%         end
%         nexogenous = size(exogenous,2);
%         if size(exogenous,1) ~= size(y,1)+ fhor && size(exogenous,1) ~= size(y,1)
%             error('Size Mismatch between endogenous and exogenos variables; exo must be either T or T+fhor');
%         end
%     end
%     %======================================================================
%     % Minnesota prior options
%     %======================================================================
%     if (isfield(options,'priors')==1 && strcmp(options.priors.name,'Minnesota')==1) || (isfield(options,'priors')==1 && strcmp(options.priors.name,'minnesota')==1) || ...
%        (isfield(options,'prior')==1 && strcmp(options.prior.name,'Minnesota')==1) || (isfield(options,'prior')==1 && strcmp(options.prior.name,'minnesota')==1)
%         %  MINNESOTA PRIOR
%         dummy = 1;
%         flat  = 0;
%         timetrend = 0;
%         priors.name= 'Minnesota';
%     end
%     if isfield(options,'minn_prior_tau')==1 || isfield(options,'bvar_prior_tau')==1
%         dummy = 1;
%         flat  = 0;
%         priors.name= 'Minnesota';
%         %  MINNESOTA PRIOR: tightness
%         if isfield(options,'bvar_prior_mu')==1
%             minn_prior_tau = options.bvar_prior_tau;
%         else
%             minn_prior_tau = options.minn_prior_tau;
%         end
%     end
%     if isfield(options,'minn_prior_decay')==1 || isfield(options,'bvar_prior_decay')==1
%         dummy = 1;
%         flat  = 0;
%         priors.name= 'Minnesota';
%         %  MINNESOTA PRIOR: decay
%         if isfield(options,'bvar_prior_mu')==1
%             minn_prior_decay = options.bvar_prior_decay;
%         else
%             minn_prior_decay = options.minn_prior_decay;
%         end
%     end
%     if isfield(options,'minn_prior_lambda')==1 || isfield(options,'bvar_prior_lambda')==1
%         dummy = 1;
%         flat  = 0;
%         priors.name= 'Minnesota';
%         %  MINNESOTA PRIOR: sum-of-coeff
%         if isfield(options,'bvar_prior_mu')==1
%             minn_prior_lambda = options.bvar_prior_lambda;
%         else
%             minn_prior_lambda = options.minn_prior_lambda;
%         end
%     end
%     if isfield(options,'minn_prior_mu')==1 || isfield(options,'bvar_prior_mu')==1
%         dummy = 1;
%         flat  = 0;
%         priors.name= 'Minnesota';
%         %  MINNESOTA PRIOR: co-persistence
%         if isfield(options,'bvar_prior_mu')==1
%             minn_prior_mu = options.bvar_prior_mu;
%         else
%             minn_prior_mu = options.minn_prior_mu;
%         end
%     end
%     if isfield(options,'minn_prior_omega')==1 || isfield(options,'bvar_prior_omega')==1
%         dummy = 1;
%         flat  = 0;
%         priors.name= 'Minnesota';
%         %  MINNESOTA PRIOR: variance
%         if isfield(options,'bvar_prior_omega')==1
%             minn_prior_omega = options.bvar_prior_omega;
%         else
%             minn_prior_omega = options.minn_prior_omega;
%         end
%     end
%     if isfield(options,'max_minn_hyper')==1 && options.max_minn_hyper ==1 && mixed_freq_on ==0
%         % maximize the hyper parameters of the minnesota prior
%         hyperpara(1) = minn_prior_tau;
%         hyperpara(2) = minn_prior_decay;
%         hyperpara(3) = minn_prior_lambda;
%         hyperpara(4) = minn_prior_mu;
%         hyperpara(5) = minn_prior_omega;
%
%         dummy = 1;
%         flat  = 0;
%         priors.name= 'Minnesota';
%         try
%             [postmode,lm,~] = bvar_max_hyper(hyperpara,y,lags,options);
%
%             minn_prior_tau      = postmode(1);
%             minn_prior_decay    = postmode(2);
%             minn_prior_lambda   = postmode(3);
%             minn_prior_mu       = postmode(4);
%             minn_prior_omega    = postmode(5);
%             disp('                       ')
%             disp('Maximization Successful: I will use the the mode values.')
%
%         catch
%             warning('Maximization NOT Successful')
%             disp('Using hyper parameter default values')
%         end
%     end
%     %======================================================================
%     % Conjugate/Hierachical prior options
%     %======================================================================
%     if (isfield(options,'priors')==1 && strcmp(options.priors.name,'Conjugate')==1) || (isfield(options,'priors')==1 && strcmp(options.priors.name,'conjugate')==1) || ...
%        (isfield(options,'prior')==1 && strcmp(options.prior.name,'Conjugate')==1) || (isfield(options,'prior')==1 && strcmp(options.prior.name,'conjugate')==1)
%
%         if isfield(options,'prior')==1
%             options.priors = options.prior;
%         end
%         if dummy == 1
%             warning('You have set both the Conjugate and Minnesota (perhaps via options.minn_prior_XX)');
%             warning('I will consider the Conjugate prior only');
%         end
%         dummy = 2;
%         %  warning('The Conjugate prior is still under construction ... ');
%         flat  = 0;
%         priors.name= 'Conjugate';
%         % Priors for the AR parameters
%         if isfield(options.priors,'Phi') == 1
%             % mean
%             if isfield(options.priors.Phi,'mean') == 1
%                 prior.Phi.mean  = options.priors.Phi.mean;
%                 if max(size(prior.Phi.mean) ~= [ny*lags+(1-noconstant)+timetrend+nexogenous   ny]) ~= 0
%                     error('Size mismatch')
%                 end
%             else
%                 warning(['You did not provide a prior mean for the AR coeff. ' ...
%                     'Assume zeros everywhere.'])
%                 prior.Phi.mean  = zeros(ny*lags+(1-noconstant)+timetrend+nexogenous , ny);
%             end
%             % variance
%             if isfield(options.priors.Phi,'cov') == 1
%                 prior.Phi.cov   = options.priors.Phi.cov;
% %                 if length(prior.Phi.cov) ~= (ny*lags+(1-noconstant) +timetrend ) * ny
%                 if length(prior.Phi.cov) ~= (ny*lags+(1-noconstant) +timetrend+nexogenous ) || size(prior.Phi.cov,1 )~=size(prior.Phi.cov,2)
%                     error('Size mismatch: Covariance Phi should be square, e.g. size(Phi.mean,1)x size(Phi.mean,1)x')
%                 end
%             else
%                 warning(['You did not provide a Covariance for the AR coeff. ' ...
%                     'Assume 10 times Identity Matrix.'])
% %                 prior.Phi.cov  = 10 * eye((ny*lags+(1-noconstant) + timetrend) * ny);
%                 prior.Phi.cov  = 10 * eye((ny*lags+(1-noconstant) + timetrend+nexogenous));
%             end
%         else
%             warning(['You did not provide prior mean and covariance for the AR coeff ' ...
%                 'Assume zeros everywhere with covariance 10 times Identity Matrix.'])
% %             prior.Phi.cov   = 10 * eye((ny*lags+(1-noconstant) +timetrend ) * ny);
%             prior.Phi.cov   = 10 * eye((ny*lags+(1-noconstant) + timetrend+nexogenous ));
%             prior.Phi.mean  = zeros(ny*lags+(1-noconstant) + timetrend+nexogenous ,  ny);
%         end
%         % Priors for the Residual Covariance
%         if isfield(options.priors,'Sigma') == 1
%             % scale
%             if isfield(options.priors.Sigma,'scale') == 1
%                 prior.Sigma.scale = options.priors.Sigma.scale;
%                 if size(prior.Sigma.scale) ~= [ny ny]
%                     error('Size mismatch')
%                 end
%             else
%                 warning(['You did not provide a prior mean for the Residual Covariance. ' ...
%                     'Assume identity matrix.'])
%                 prior.Sigma.scale = eye(ny);
%             end
%             % degrees of freedom
%             if isfield(options.priors.Sigma,'df') == 1
%                 prior.Sigma.df = options.priors.Sigma.df;
%                 if length(prior.Sigma.df) ~= 1
%                     error('Size mismatch')
%                 end
%                 if prior.Sigma.df/2 <= ny-1
%                     error('Too few degrees of freedom - Increase prior df')
%                 end
%             else
%                 warning(['You did not provide the degrees of freedom for the Residual Covariance. ' ...
%                     'Assume ny+1.'])
%                 prior.Sigma.df = ny + nexogenous + timetrend + 1;
%                 while prior.Sigma.df/2 <= ny-1 % too few df
%                     prior.Sigma.df =  prior.Sigma.df +1;
%                 end
%             end
%         else
%             warning(['You did not provide prior mean and variance for the Residual Covariance. ' ...
%                 'Assume an identity matrix matrix with N+1 degrees of freedom.'])
%             prior.Sigma.scale = eye(ny);
%             prior.Sigma.df    = ny + nexogenous + timetrend + 1;
%             while prior.Sigma.df/2 <= ny-1 % too few df
%                 prior.Sigma.df =  prior.Sigma.df +1;
%             end
%         end
%     end
%     %======================================================================
%     % IRF options
%     %======================================================================
%     if isfield(options,'hor') ==1
%         hor = options.hor;
%     end
%     if isfield(options,'long_run_irf')==1
%         % Activating Long run IRF
%         long_run_irf = options.long_run_irf;
%     end
%     if isfield(options,'signs')==1
%         % Activating IRF with sign restrictions (mulitple horizons allowed)
%         signs_irf       = 1;
%         signs           = options.signs;
%         if iscellstr(signs) == 0
%             error(['options.signs should be a cell array. Each cell must contain a string with the format'...
%                 '\n''y(a,b,c)<0'' or ''y(a,b,c)>0'' where a, b and c are integers.',...
%                 '\na = index of the variable',...
%                 '\nb = horizon',...
%                 '\nc = index of the shock'],class(zeros))
%         end
%     end
%     if isfield(options,'narrative')==1
%         if signs_irf  == 0
%             warning('You did not provide any sign restrictions.')
%             signs{1} = 'isempty(y(1,1,1))==0';
%         end
%         narrative_signs_irf = 1;
%         narrative           = options.narrative ;
%     end
%     if isfield(options,'zeros_signs')==1
%         if signs_irf  == 1
%             signs_irf       = 0;  % disactivating signs
%         end
%         if narrative_signs_irf  == 1
%             narrative_signs_irf       = 0;  % disactivating narrative
%         end
%         % Activating IRF with zeros and sign restrictions (mulitple horizons NOT allowed)
%         zeros_signs_irf  = 1;
%         zeros_signs      = options.zeros_signs;
%         if iscellstr(zeros_signs) == 0
%             error(['options.zeros_signs should be a cell array.'...
%                 '\nEach cell must contain a string with the following format'...
%                 '\nFor sign restrictions ''y(a,b)=1'' or ''y(a,b)=-1'',',...
%                 '\nFor short run zero restriction ''ys(a,b)=0'',',...
%                 '\nFor long run restriction ''yr(a,1,b)=0'' where a and b are integers.',...
%                 '\na = index of the variable',...
%                 '\nb = index of the shock'],class(zeros_signs))
%         end
%         [f,sr] = sign2matrix(zeros_signs,ny);
%         if isfield(options,'var_pos')==0
%             var_pos = ones(1,ny);
%         else
%             var_pos = options.var_pos;
%         end
%     end
%     if isfield(options,'proxy')==1
%         % Activating IRF with provy
%         proxy_irf   = 1;
%         in.proxies  = options.proxy;
%         in.vars     = y;
%         in.p        = lags;
%         if isfield(options,'proxy_end') == 1
%             in.T_m_end  = options.proxy_end;
%         else
%             in.T_m_end  = 0;  %if the times series of the instrument ends when VAR data ends
%         end
%         in.irhor    = hor;
%         if isnumeric(in.proxies) == 0
%             error(['options.proxy should be a numeric array (nans or inf not allowed)'],class(in.proxies))
%         end
%     end
%     if isfield(options,'irf_1STD')==1
%         % Activating of unitary IRF, i.e. a unitary increase in the shocks
%         % (instead of 1 STD)
%         irf_1STD = options.irf_1STD;
%     end
%     %======================================================================
%     % (Un)Conditional Forecasts options
%     %======================================================================
%     if isfield(options,'fhor')==1
%         fhor = options.fhor;
%         if fhor < 1
%             error('Forecast horizon must be positive')
%         end
%     end
%     if isfield(options,'endo_index')==1
%         % Forecast conditional on the path of an endogenous var
%         cfrcst_yes      = 1;
%         if isfield(options,'endo_path')== 0
%             error('You need to provide the p[ath for the endogenou variable')
%         end
%         % rows forecasts, column variables
%         endo_path       = options.endo_path;
%         endo_path_index = options.endo_index;
%         if length(endo_path_index) ~= size(endo_path)
%             error(['Mismatch beween the number of endogenous paths and the number of conditioned variables'...
%                 '\nE.g. the # of conditioned variables must coincide with the # of column in ''options.endo_path'''],class(endo_path_index));
%         end
%         if isfield(options,'exo_index')==1
%             % Forecast conditional on the path of an endo var using only a
%             % subset of shocks. notice that the # of endo and # exo must coincide
%             cfrcst_yes  = 2;
%             exo_index   = options.exo_index;
%             Omega       = eye(ny);
%             %             if isfield(options,'Omegaf')==1
%             %                 Omegaf = options.Omegaf;
%             %             end
%             if length(exo_index) ~= length(endo_path_index)
%                 error('the # of conditioned endogenous and # exogenous shocks used must coincide');
%             end
%         end
%
%     end
%     %======================================================================
%     % Missing Values options
%     %======================================================================
%     if mixed_freq_on == 1 && isfield(options,'mixed_freq_index')==1
%         index = options.mixed_freq_index;
%         if length(index) ~= size(y,2)
%             error(['You have to specify as many index as observables.'...'
%                 '\nIf no missing values for var j index(j)=0'...
%                 '\nIf missing values for var j, and var j is a stock index(j)=0'...
%                 '\nIf missing values for var j, and var j is a real flow index(j)=2'],class(index))
%         end
%     elseif mixed_freq_on == 1 && isfield(options,'mf_varindex')== 1
%         index = zeros(ny,1);
%         index(options.mf_varindex) = 2;
%     elseif (mixed_freq_on == 1 && isfield(options,'mixed_freq_index')== 0) || (mixed_freq_on == 1 && isfield(options,'mf_varindex')== 0)
%         warning(['You did not specified the aggregation of the mixed freq.'...
%             '\nI will treat them as stocks.']);
%         index = zeros(size(y,2),1);
%     end
%     if isfield(options,'noprint')==1
%         noprint = options.noprint;
%     end
% end
%
% %********************************************************
% %* Consistency Checks
% %********************************************************
% nobs = size(y,1)-firstobs+1;
% if (firstobs+ nobs-1)> size(y,1)
%     fprintf('Incorrect or missing specification of the number of observations. nobs can be at most %4u\n',size(y,1)-firstobs+1);
%     error('Inconsistent number of observations.')
% end
% if firstobs + presample + lags >= nobs
%     error('presample too large')
% end
% if firstobs + presample  <= lags
%     error('firstobs+presample should be > # lags (for initializating the VAR)')
% end
% if dummy == 1 && nexogenous > 0
%     warning('I will not use exogenous variables with Minnesota');
%     nexogenous = 0;
% end
% if mixed_freq_on == 1 && nexogenous > 0
%     warning('I will not use exogenous variables with missing observations');
%     nexogenous = 0;
% end
% if cfrcst_yes ~= 0 && nexogenous > 0
%     warning('I will not use exogenous variables with conditional forecasts');
%     nexogenous = 0;
% end
% if size(exogenous,1) == size(y,1) && nexogenous > 0
%     warning('For forecast purposes, I will assume that exo are zero out-of sample.')
%     fprintf('To change this, include the exogenous forecasts in options.exogenous.\n')
%     exogenous = [exogenous; zeros(fhor,nexogenous)];
% end
%
% %********************************************************
% %* Priors and Posterior Distributions
% %********************************************************
%
% idx = firstobs+presample-lags:firstobs+nobs-1;
% nx  = 1;
% if noconstant
%     nx = 0;
% end
%
% % organize data as  yy = XX B + E
% [yy,XX] = YXB_(y(idx, :),lags,[nx timetrend]);
%
% ydum    = [];
% xdum    = [];
% pbreaks = 0;
% lambda  = 0;
% mu      = 0;
%
% ydata   = y(idx, :);
% T       = size(ydata, 1);
% if T-lags < lags*ny + nx %+ flat*(ny+1)
%     error('Less observations than regressors: increase the # of obs or decrease the # of lags.')
% end
% xdata   = ones(T,nx);
% if timetrend == 1
%     % xdata = [xdata [1:T]'];
%     xdata = [xdata [1-lags : T-lags]'];
% end
% if nexogenous > 0
%     xdata = [xdata exogenous(idx,:)];
%     XX    = [XX exogenous(idx(1)+lags : idx(end),:)];
% end
%
% % OLS estimate [NO DUMMY]:
% if heterosked == 0
%     varols  = rfvar3(ydata, lags, xdata, [T; T], 0, 0);
% else
%     varols  = rfvar3(ydata, lags, xdata, [T; T], 0, 0, ww);
% end
%
%
% % specify the prior
% if dummy == 1
%     % MINNESOTA PRIOR:
%     mnprior.tight         = minn_prior_tau;
%     mnprior.decay         = minn_prior_decay;
%     % Use only initializations lags for the variance prior
%     % vprior.sig            = std(y(firstobs+presample-lags : firstobs+presample,:))';
%     vprior.sig            = std(y(firstobs-lags : firstobs+presample-1,:))';
%     vprior.w              = minn_prior_omega;
%     lambda                = minn_prior_lambda;
%     mu                    = minn_prior_mu;
%     [ydum, xdum, pbreaks] = varprior(ny, nx, lags, mnprior, vprior);
%     % Prior density
%     Tp = presample + lags;
%     if nx
%         xdata = xdata(1:Tp, :);
%     else
%         xdata = [];
%     end
% %     varp            = rfvar3([ydata(1:Tp, :); ydum], lags, [xdata; xdum], [Tp; Tp + pbreaks], lambda, mu);
%     varp           = rfvar3([y(firstobs-lags : firstobs+presample-1, :); ydum], lags, [xdata; xdum], [Tp; Tp + pbreaks], lambda, mu);
%     Tup            = size(varp.u, 1);
%     prior.df       = Tup - ny*lags - nx - flat*(ny+1);
%     prior.S        = varp.u' * varp.u;
%     prior.XXi      = varp.xxi;
%     prior.PhiHat   = varp.B;
%     prior.YYdum    = varp.y;
%     prior.XXdum    = varp.X;
%     if prior.df < ny
%         error('Too few degrees of freedom in the Inverse-Wishart part of prior distribution. You should increase training sample size.')
%     end
%     prior.minn_prior_tau    = minn_prior_tau;
%     prior.minn_prior_decay  = minn_prior_decay;
%     prior.minn_prior_lambda = minn_prior_lambda;
%     prior.minn_prior_mu     = minn_prior_mu;
%     prior.minn_prior_omega  = minn_prior_omega;
% elseif dummy == 0
%     % JEFFREY OR UNIFORMATIVE PRIOR:
%     % [priors] = jeffrey(y,lags);
%     prior.name  = 'Jeffrey';
%
% end
%
% % specify the posterior (the varols agin on actual+dummy)
% [posterior,var] = posterior_(y);
%
%
% %*********************************************************
% %* Compute the log marginal data density for the VAR model
% %*********************************************************
%
% try
%     posterior_int = matrictint(posterior.S, posterior.df, posterior.XXi);
% catch
%     warning('I could not compute the likelihood');
%     posterior_int = nan;
% end
%
% if dummy == 1 % only for minnesota dummy
%     prior_int = matrictint(prior.S, prior.df, prior.XXi);
%     lik_nobs  = posterior.df - prior.df;
%     log_dnsty = posterior_int - prior_int - 0.5*ny*lik_nobs*log(2*pi);
%
% elseif dummy == 0 % jeffrey prior
%     lik_nobs  = posterior.df;
%     log_dnsty = posterior_int - 0.5*ny*lik_nobs*log(2*pi);
%
% elseif dummy == 2 % conjugate MN-IW prior
%     log_dnsty = mniw_log_dnsty(prior,posterior,varols);
%
% end
%
%
% %**************************************************
% %* Generating draws form the Posterior Distribution
% %**************************************************
%
% % Preallocation of memory
% % Matrices for collecting draws from Posterior Density
% % Last dimension corresponds to a specific draw
% Phi_draws     = zeros(ny*lags+nx+timetrend + nexogenous, ny, K);   % Autoregressive Parameters
% Sigma_draws   = zeros(ny,ny,K);                     % Shocks Covariance
% ir_draws      = zeros(ny,hor,ny,K);                 % variable, horizon, shock and draws - Cholesky IRF
% irlr_draws    = zeros(ny,hor,ny,K);                 % variable, horizon, shock and draws - Long Run IRF
% Qlr_draws     = zeros(ny,ny,K);                     % long run impact matrix
% e_draws       = zeros(size(yy,1), ny,K);                  % residuals
% yhatfut_no_shocks         = NaN(fhor, ny, K);   % forecasts with shocks
% yhatfut_with_shocks       = NaN(fhor, ny, K);   % forecast without the shocks
% yhatfut_cfrcst            = NaN(fhor, ny, K);   % forecast conditional on endogenous path
% logL                      = NaN(K,1);
%
% if signs_irf == 1
%     irsign_draws = ir_draws;
%     Omega_draws    = Sigma_draws;
% end
% if nexogenous > 0
%     irx_draws = zeros(ny,hor,nexogenous,K);
%     % Ox_draws  = Sigma_draws;
% end
% if narrative_signs_irf == 1
%     irnarrsign_draws = ir_draws;
%     Omegan_draws     = Sigma_draws;
% end
% if zeros_signs_irf == 1
%     irzerosign_draws   = ir_draws;
%     Omegaz_draws       = Sigma_draws;
% end
% if proxy_irf == 1
%     irproxy_draws = ir_draws;
% end
%
% % Settings for the forecasts
% forecast_data.xdata       = ones(fhor, nx);
% if timetrend
%     forecast_data.xdata = [forecast_data.xdata (T-lags+1 : T-lags+fhor)'];
% end
% if nexogenous>0
%     forecast_data.xdata = [forecast_data.xdata exogenous(T-lags+1 : T-lags+fhor,:)];
% end
% forecast_data.initval     = ydata(end-lags+1:end, :);
%
% try
%     S_inv_upper_chol    = chol(inv(posterior.S));
% catch
%     warning('POSTERIOR MEAN of SIGMA IS ILL-BEHAVED (NON-POSITIVE DEFINITE)')
%     try
%         warning('I will try with the LDL decomposition')
%         iS               = inv(posterior.S);
%         [L, D, P]        = ldl(iS);
%         S_inv_upper_chol = sqrt(D) * L' * P';%chol()
%     catch
%         warning('I will add +1e-05 to the diagonal')
%         S_inv_upper_chol    = chol(inv(posterior.S + 1e-05*eye(ny)));
%     end
% end
%
%
% XXi_lower_chol      = chol(posterior.XXi)';
%
% nk                  = ny*lags+nx+timetrend + nexogenous;
%
% % Declaration of the companion matrix
% Companion_matrix = diag(ones(ny*(lags-1),1),-ny);
% p  = 0;
% dd = 0;
%
% waitbar_yes = 0;
% if K > 99
%     waitbar_yes = 1;
%     wb = waitbar(0, 'Generating draws from the Posterior Distribution');
% end
%
% for  d =  1 : K
%
%     while dd == 0
%         %======================================================================
%         % Inferece: Drawing from the posterior distribution
%         % Step 1: draw from the Covariance
%         Sigma = rand_inverse_wishart(ny, posterior.df, S_inv_upper_chol);
%
%         % Step 2: given the Covariance Matrix, draw from the AR parameters
%         Sigma_lower_chol = chol(Sigma)';
%         Phi1 = randn(nk * ny, 1);
%         Phi2 = kron(Sigma_lower_chol , XXi_lower_chol) * Phi1;
%         Phi3 = reshape(Phi2, nk, ny);
%         Phi  = Phi3 + posterior.PhiHat;
%
%         Companion_matrix(1:ny,:) = Phi(1:ny*lags,:)';
%         test = (abs(eig(Companion_matrix)));
%         if non_explosive_ == 1
%             % Checking the eigenvalues of the companion matrix (on or inside the
%             % unit circle)
%             if any(test>1.01) %any(test>1.0000000000001)
%                 p = p+1;
%             else
%                 dd = 1;
%             end
%         else
%             dd = 1;
%         end
%     end
%
%     % store the draws
%     Phi_draws(:,:,d)   = Phi;
%     Sigma_draws(:,:,d) = Sigma;
%     errors             = yy - XX * Phi;
%     e_draws(:,:,d)     = errors;
%
%
%     %======================================================================
%     % IRF
%     % Compute the impulse response functions
%     % with cholesky
%     if irf_1STD == 1
%         % one STD increase
%         ir_draws(:,:,:,d)      = iresponse(Phi,Sigma,hor,eye(ny));
%     else
%         % one percent increase
%         ir_draws(:,:,:,d)      = iresponse(Phi,Sigma,hor,eye(ny),0);
%     end
%     if nexogenous > 0
%         Phi1         = Phi(1 : ny*lags + nx + timetrend, 1 : ny);
%         Exo1         = zeros(ny,nexogenous);
%         Exo1(:,1:nx) = Phi(ny*lags + nx + timetrend +1 ....
%                             : ny*lags + nx + timetrend + nexogenous, 1 : ny)';
%         Sig1         = eye(ny);
%         % unitary increase
%         irx_draws(:,:,:,d) = iresponse(Phi1,Sig1,hor,Exo1);
%     end
%
%     % with long run restrictions
%     if long_run_irf == 1
%         [irlr,Qlr]             = iresponse_longrun(Phi,Sigma,hor,lags);
%         irlr_draws(:,:,:,d)    = irlr;
%         Qlr_draws(:,:,d)       = Qlr;
%     end
%     % with sign restrictions
%     if signs_irf == 1
%         [irsign,Omega]         = iresponse_sign(Phi,Sigma,hor,signs);
%         irsign_draws(:,:,:,d)  = irsign;
%         Omega_draws(:,:,d)     = Omega;
%     end
%     % with narrative and sign restrictions
%     if narrative_signs_irf == 1
%         [irnarrsign,Omega]         = iresponse_sign_narrative(errors,Phi,Sigma,hor,signs,narrative);
%         irnarrsign_draws(:,:,:,d)  = irnarrsign;
%         Omegan_draws(:,:,d)        = Omega;
%     end
%     % with zeros and sign restrictions
%     if zeros_signs_irf == 1         %= iresponse_zeros_signs( Phi,Sigma,bvar1.hor,lags,var_pos,f,sr);
%         [irzerosign,Omega]          = iresponse_zeros_signs(Phi,Sigma,hor,lags,var_pos,f,sr);
%         irzerosign_draws(:,:,:,d)   = irzerosign;
%         Omegaz_draws(:,:,d)         = Omega;
%     end
%     % with proxy
%     if proxy_irf == 1
%         in.res                  = e_draws(:,:,d);
%         in.Phi                  = Phi_draws(:,:,d)  ;
%         in.Sigma                = Sigma;
%         tmp_                    = iresponse_proxy(in);
%         irproxy_draws(:,:,1,d)  = tmp_.irs';
%         clear tmp_
%     end
%
%     %======================================================================
%     % Forecasts
%     % compute the out of sample forecast (unconditional)
%     [frcst_no_shock,frcsts_with_shocks] = forecasts(forecast_data,Phi,Sigma,fhor,lags);
%     yhatfut_no_shocks(:,:,d)            = frcst_no_shock;
%     yhatfut_with_shocks(:,:,d)          = frcsts_with_shocks;
%
%     if cfrcst_yes == 1
%         % Forecast conditional on the path of an endo var using all shocks
%         [sims_with_endopath,EPS(:,:,d)] = ...
%             cforecasts(endo_path,endo_path_index,forecast_data,Phi,Sigma);
%         yhatfut_cfrcst(:,:,d) = sims_with_endopath;
%
%     elseif cfrcst_yes == 2
%         % Forecast conditional on the path of an endo var using only a
%         % subset of shocks. notice that the # of endo and # exo must coincide
%         % Omega is the structural orthonormal matrix
%         [sims_with_endopath,EPS(:,:,d)] = ...
%             cforecasts2(endo_path,endo_path_index,exo_index,forecast_data,Phi,Sigma,Omega);
%         yhatfut_cfrcst(:,:,d) = sims_with_endopath;
%     end
%
%
%     %======================================================================
%     % Mixed Frequency
%     if mixed_freq_on
%         % Checking the eigenvalues of the companion matrix (on or inside the
%         % unit circle). Needed for the initial of KF
%         Companion_matrix(1:ny,:) = Phi(1:ny*lags,:)';
%         %test = (abs(eig(Companion_matrix)));
%         if any(test>1.0000000000001)
%             KFoptions.initialCond = 1;
%         end
%         KFoptions.index   = index;
%         KFoptions.noprint = noprint;
%         % Forward Kalman Filter and Smoother
%         [logL,KFout]          = kfilternan(Phi,Sigma,yoriginal,KFoptions);
%         yfill(:,:,d)     = KFout.smoothSt_plus_ss(:,KFout.index_var);
%         yfilt(:,:,d)     = KFout.filteredSt_plus_ss(:,KFout.index_var);
%         % recompute the posterior with smoothed data
%         [posterior1]     = posterior_(yfill(:,:,d));
%         S_inv_upper_chol = chol(inv(posterior1.S));
%         XXi_lower_chol   = chol(posterior1.XXi)';
%         posterior.PhiHat = posterior1.PhiHat;
%         forecast_data.initval = yfill(end-lags+1:end, :, d);
%         logL(d) = KFout.logL;
%     end
%     if waitbar_yes, waitbar(d/K, wb); end
%     dd = 0; % reset
% end
% if waitbar_yes, close(wb); end
%
%
% %********************************************************
% %* Storing the resutls
% %*******************************************************
%
% % classical inference: OLS estimator
% BVAR.Phi_ols    = varols.B;
% BVAR.e_ols      = varols.u;
% BVAR.Sigma_ols  = 1/(nobs-nk)*varols.u'*varols.u;
% [BVAR.InfoCrit.AIC, BVAR.InfoCrit.HQIC, BVAR.InfoCrit.BIC] = IC(BVAR, nobs, nk);
% % the model with the lowest IC is preferred
%
% % OLS irf
% BVAR.ir_ols      = iresponse(BVAR.Phi_ols,BVAR.Sigma_ols,hor,eye(ny));
% if long_run_irf == 1
%     [irlr,Qlr]              = iresponse_longrun(BVAR.Phi_ols,BVAR.Sigma_ols,hor,lags);
%     BVAR.irlr_ols           = irlr;
%     BVAR.Qlr_ols(:,:)       = Qlr;
% end
%
% % test the normality of the ols VAR residuals (matlab stat toolbox needed)
% if  exist('kstest') ==2
%     for gg = 1 : ny
%         [H,Pv] = kstest(BVAR.e_ols(:,gg)/sqrt(BVAR.Sigma_ols(gg,gg)));
%         BVAR.HP(gg,:) = [H,Pv];
%     end
% else
%     BVAR.HP = [];
% end
%
% % bayesian inference:
% % the last dimension of these objects corresponds to a draw from the posterior
%
% % inference and IRFs
% BVAR.Phi_draws    = Phi_draws;          % draws from the autoregressive part
% BVAR.Sigma_draws  = Sigma_draws;        % draws from the covarance matrix
% BVAR.alpha_draws  = Phi_draws;          % Older name: draws from the autoregressive part
% BVAR.sigma_draws  = Sigma_draws;        % Older name: draws from the covarance matrix
%
% BVAR.ir_draws     = ir_draws;           % draws from the IRF with cholesky
% BVAR.irlr_draws   = irlr_draws;         % draws from the IRF with Long Run
% BVAR.Qlr_draws    = Qlr_draws;          % Long Run Rotation matrix
% BVAR.lags         = lags;               % lags
% BVAR.N            = ny;                 % number of variables
% BVAR.e_draws      = e_draws;            % residuals
% BVAR.e            = e_draws;            % backward compatible with earlier versions
%
% BVAR.posterior    = posterior;
% BVAR.prior        = prior;             % priors used
% BVAR.logmlike     = log_dnsty;
% BVAR.X            = var.X;              % regressors (including dummy if Minnesota)
% BVAR.y            = var.y;              % dependent  (including dummy if Minnesota)
% BVAR.XX           = XX;                 % regressors (no dummy)
% BVAR.yy           = yy;                 % dependent  (no dummy)
%
% % prediction
% BVAR.fhor         = fhor;               % forecast horizon
% BVAR.hor          = hor;                % IRF horizon
% BVAR.forecasts.no_shocks      = yhatfut_no_shocks;         % trajectories of forecasts without shocks
% BVAR.forecasts.with_shocks    = yhatfut_with_shocks;       % trajectories of forecasts with shocks
% BVAR.forecasts.conditional    = [];            % trajectories of conditional forecasts
% BVAR.forecasts.EPScond        = [];            % shocks of conditional forecasts
% if cfrcst_yes ~= 0
%     BVAR.forecasts.conditional    = yhatfut_cfrcst;            % trajectories of forecasts
%     BVAR.forecasts.EPScond        = EPS;                       % shocks of forecasts
% end
% BVAR.forecast_data            = forecast_data;
%
% %
% BVAR.varnames     = varnames;
% BVAR.ndraws       = K;
%
% if signs_irf == 1 && narrative_signs_irf == 0
%     BVAR.irsign_draws = irsign_draws;
%     BVAR.Omegas       = Omega_draws;
% else
%     BVAR.irsign_draws = [];
%     BVAR.Omegas       = [];
% end
% if narrative_signs_irf == 1
%     BVAR.irnarrsign_draws = irnarrsign_draws;
%     BVAR.Omegan           = Omegan_draws;
% else
%     BVAR.irnarrsign_draws = [];
%     BVAR.Omegan           = [];
% end
% if zeros_signs_irf == 1
%     BVAR.irzerosign_draws   = irzerosign_draws;
%     BVAR.Omegaz             = Omegaz_draws;
% else
%     BVAR.irzerosign_draws = [];
%     BVAR.Omegaz           = [];
% end
% if proxy_irf == 1
%     BVAR.irproxy_draws = irproxy_draws;
% else
%     BVAR.irproxy_draws= [];
% end
% if nexogenous > 0
%     BVAR.irx_draws = irx_draws;
% else
%     BVAR.irx_draws = [];
% end
%
% % missing observations
% BVAR.data         = y;                  % raw data
% if mixed_freq_on
%     BVAR.yfill = yfill;
%     BVAR.yfilt = yfilt;
%     BVAR.yinterpol = yinterpol;
%     BVAR.logL  = logL;
% end
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % end of bvar.m
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% %********************************************************
% %********************************************************
%     function [posterior,var] = posterior_(y)
%
%         % This part is needed in the case of missing values. The posterior
%         % needs to be revaluated given a new value of the Kalman smoothed
%         % observables
%         %===========================================
%         ydata = y(idx, :);
%         %===========================================
%         T     =  size(ydata, 1);
%         xdata = ones(T,nx);
%         if timetrend ==1
%             xdata = [xdata [1-lags:T-lags]'];
%             % xdata = [xdata [1:T]'];
%         end
%         if nexogenous >0
%             xdata = [xdata exogenous(idx,:)];
%         end
%         % posterior density
%         if heterosked == 0
%             var = rfvar3([ydata; ydum], lags, [xdata; xdum], [T; T+pbreaks], lambda, mu);
%         else
%             var = rfvar3([ydata; ydum], lags, [xdata; xdum], [T; T+pbreaks], lambda, mu, ww);
%         end
%         Tu = size(var.u, 1);
%
%         if dummy ==  1
%             %********************************************************
%             % Minnesota Prior
%             %********************************************************
%             % Prior density
%             Tp = presample + lags;
%             if nx
%                 xdata = xdata(1:Tp, :);
%             else
%                 xdata = [];
%             end
%             % varp            = rfvar3([ydata(1:Tp, :); ydum], lags, [xdata; xdum], [Tp; Tp + pbreaks], lambda, mu);
%             varp            = rfvar3([y(firstobs-lags : firstobs+presample-1, :); ydum], lags, [xdata; xdum], [Tp; Tp + pbreaks], lambda, mu);
%             Tup             = size(varp.u, 1);
%             prior.df        = Tup - ny*lags - nx - flat*(ny+1);
%             prior.S         = varp.u' * varp.u;
%             prior.XXi       = varp.xxi;
%             prior.PhiHat    = varp.B;
%             priors.YYdum    = varp.y;
%             priors.XXdum    = varp.X;
%             if prior.df < ny
%                 error('Too few degrees of freedom in the Inverse-Wishart part of prior distribution. You should increase training sample size.')
%             end
%             posterior.df    = Tu - ny*lags - nx - flat*(ny+1);
%             posterior.S     = var.u' * var.u;
%             posterior.XXi   = var.xxi;
%             posterior.PhiHat = var.B;
%
%         elseif dummy == 2
%             %********************************************************
%             % Conjugate Prior
%             %********************************************************
%             Ai              = inv(prior.Phi.cov);
%             posterior.df    = Tu - ny*lags - nx - nexogenous - prior.Sigma.df;
%             posterior.XXi   = inv(var.X'*var.X + Ai);
%             posterior.PhiHat = posterior.XXi * (var.X' * var.y + Ai * prior.Phi.mean);
%             %posterior.S     = var.u' * var.u + prior.Sigma.scale ;
%             %posterior.S     = var.u' * var.u + prior.Sigma.scale  + (var.B - prior.Phi.mean)' * Ai * (var.B - prior.Phi.mean);
%             %posterior.S     = var.u' * var.u + prior.Sigma.scale  + ...
%             %    (posterior.PhiHat - prior.Phi.mean)' * Ai * (var.B - prior.Phi.mean);
%             % check
%             posterior.S = var.u' * var.u + prior.Sigma.scale + ...
%                 prior.Phi.mean' * Ai * prior.Phi.mean + ...
%                 var.B' * (var.X'*var.X) * var.B ...
%                 - posterior.PhiHat' * (var.X'*var.X + Ai) * posterior.PhiHat;
%
%         else
%             %********************************************************
%             % Flat Jeffrey Prior
%             %********************************************************
%             posterior.df    = Tu - ny*lags - nx + flat*(ny+1) - nexogenous;
%             posterior.S     = var.u' * var.u;
%             posterior.XXi   = var.xxi;
%             posterior.PhiHat = var.B;
%         end
%
%     end
% %********************************************************
% %********************************************************
%     function [priors] = priors_( )
%
%         priors.name  = 'N/A';
%
%     end
% %********************************************************
% %********************************************************
%
% end


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

%%
clear; close all; clc;

%% Syntax
% |[M]=unpadeven(M_pad,L_even_dim);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef TestTwoRankSequentialLaminate < handle

    properties (Access = public)
        tol = 1e-12;
    end

    properties (Access = private)
      C0
      C1
      mu
      lambda
      lambda2D
      FractionVolume
      d1
      d2
      m1
      m2
      FirstCheckedRank2Ch
      WebRank2Ch
      Rank2Ch
    end

    methods (Access = public)

        function obj = TestTwoRankSequentialLaminate()
            obj.createTensors();
            obj.createParameters();
            obj.loadHomogenizedCheckedConstitutiveTensor();
            obj.computeTwoRankSequentialLaminate();
            obj.computeTwoRankSequentialLaminateFromAllaireTestedWebPageCode();
        end

        function error = computeError(obj)
            InitCh  = double(obj.FirstCheckedRank2Ch);
            R2Ch    = double(obj.Rank2Ch);
            R2ChWeb = double(obj.WebRank2Ch);
            err1  = norm(InitCh - R2Ch)/norm(R2Ch);
            err2  = norm(R2ChWeb- R2Ch)/norm(R2Ch);
            error = max(err1,err2);
        end

    end

    methods (Access = private)

        function createTensors(obj)
            obj.createStiffTensor();
            obj.createWeakTensor();
            obj.storeLambda2D();
            obj.makeTensorsVoigtPlaneStress();
        end

        function createStiffTensor(obj)
            obj.lambda = 0.7500;
            obj.mu     = 0.3750;
            obj.C1 = IsotropicConstitutiveTensor.createWithLambdaAndMu(obj.lambda,obj.mu);
        end

        function createWeakTensor(obj)
            E  = obj.C1.getYoung();
            nu = obj.C1.getPoisson();
            epsil = 1.0000e-03;
            E0 = epsil*E;
            nu0 = nu;
            obj.C0 = IsotropicConstitutiveTensor(E0,nu0);
        end

        function storeLambda2D(obj)
            obj.lambda2D = obj.C1.getLambda2D();
        end

        function makeTensorsVoigtPlaneStress(obj)
            obj.C0 = obj.makeTensorVoigtPlaneStress(obj.C0);
            obj.C1 = obj.makeTensorVoigtPlaneStress(obj.C1);
        end

        function createParameters(obj)
            obj.loadDirections()
            obj.loadLaminationParameters()
            obj.loadFractionVolume()
        end

        function loadDirections(obj)
            obj.loadFirstDirection()
            obj.loadSecondDirection()
        end

        function loadFirstDirection(obj)
            dir = [1     0     0];
            obj.d1 = Vector3D;
            obj.d1.setValue(dir);
            obj.d1.normalize()
        end

        function loadSecondDirection(obj)
            dir = [1     3     2];
            obj.d2 = Vector3D;
            obj.d2.setValue(dir);
            obj.d2.normalize()
        end

        function loadLaminationParameters(obj)
            obj.m1  = 0.8;
            obj.m2  = 0.2;
        end

        function loadFractionVolume(obj)
            obj.FractionVolume = 0.8000;
        end

        function loadHomogenizedCheckedConstitutiveTensor(obj)
            obj.FirstCheckedRank2Ch = [    0.326824930177257   0.028535636891760  -0.155777642710044
                                           0.028535636891760   0.746824571833621  -0.074187675387998
                                          -0.155777642710044  -0.074187675387998   0.032148130248805];
        end

        function computeTwoRankSequentialLaminate(obj)
            dir{1} = obj.d1;
            dir{2} = obj.d2;
            params = [obj.m1;obj.m2];
            theta = obj.FractionVolume;
            c1    = obj.C1.getValue();
            c0    = obj.C0.getValue();
            lam2D = obj.lambda2D;
            muV   = obj.mu;
            Homogenizer = RankTwoLaminateHomogenizer(c1,c0,dir,params,...
                                                     theta,lam2D,muV);
            obj.Rank2Ch = Homogenizer.getTensor();
        end

        function computeTwoRankSequentialLaminateFromAllaireTestedWebPageCode(obj)
            %Lambda   = obj.lambda2D;
            Lambda2D = obj.lambda2D;
            muV = obj.mu();
            theta = obj.FractionVolume;
            epsil = 1.0000e-03;
            dir1 = obj.d1.getValue();
            dir2 = obj.d2.getValue();
            obj.WebRank2Ch = RankTwoLaminateHomogenizerFromAllaireTestedWebPageCode.ChinvOld(...
                             Lambda2D,muV,dir1,dir2,obj.m1,obj.m2,theta,epsil);
        end

    end

    methods (Access = private, Static)

        function CVoigtPS = makeTensorVoigtPlaneStress(C)
            CVoigt = Tensor2VoigtConverter.convert(C);
            CVoigtPS = PlaneStressTransformer.transform(CVoigt);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function letitsnow(figure_hdl, sf_color)
% LETITSNOW - Make the figure snow!
%
% Usage:
% ------
% letitsnow              -  Create snow in the CURRENT figure. If no figure
%                           exists, an empty figure will be created.
% letitsnow(figure_hdl)  -  Create snow in the figure represented by
%                           figure_hdl
% letitsnow(figure_hdl, color) - Create snow with the specified color
%                                'color' is a scalar ranged 0~1
%
% How to stop the snow
% -----
% Press 's' on the figure or directly close the figure
%
% Mingjing Zhang @ Simon Fraser University, Canada
% Copyright(C)Stellari Studio, 2011

%% Update History
% -------------------------------------------------------------------------
% Date                  Updater                 Modification
% -------------------------------------------------------------------------
% Dec 19, 2010          M. Zhang                wrote it
%%

%% TODO:
% The frost effect

%% Error Check
if ~exist('figure_hdl', 'var')||isempty(figure_hdl)
    figure_hdl = gcf;
%     axes
end
if ~exist('sf_color', 'var') || isempty(sf_color)
    sf_color = 0.9;
end

if strcmp(get(figure_hdl, 'UserData'), 'sn0w1ng')
    return;
end
%% Constant Declaration
fps = 25;
% - For snowflakes -
num_snowflakes = 120;
snowflakes_size = 11;
snowflakes_color = sf_color.*[1 1 1];
height_range = 3;

snow_timer = timer('TimerFcn', @updateSnow_Callback_nest,...
    'ExecutionMode','fixedRate',...
    'Period', 1./fps, 'BusyMode', 'queue');

% Position of snowflakes
snowflakes_x = rand(1, num_snowflakes);

snowflakes_y = rand(1, num_snowflakes);
snowflakes_y = (snowflakes_y .* height_range.^3).^(1/3)+1;
snowflakes_pos = [snowflakes_x; snowflakes_y];

% Velocity of snowflakes
snowflakes_v = rand(1, num_snowflakes).*0.1 + 0.1;  % (distance per second)
snowflakes_theta = (rand(1, num_snowflakes).*60 + 240).*pi/180;
snowflakes_vel = [snowflakes_v.*cos(snowflakes_theta); snowflakes_v.*sin(snowflakes_theta)];

% parameters for the frost image (not implemented for this version)
% frost_image_res = [300 400];
% frost_image_data = ones(frost_image_res(1), frost_image_res(2), 3);
% frost_alpha_data = zeros(frost_image_res); %rand(frost_image_res).*0.4+0.6;
% frost_snowpile_top = ones(1, frost_image_res(2)).*frost_image_res(1);

% direction_vectors = [1 0;1 -1;0 -1;-1 -1;-1 0;-1 1;0 1;1 1];

%% Create two axes (above all the existing axes, if any)
% Modify the figure (just a little bit)
original_figure_title = get(figure_hdl, 'Name');
set(figure_hdl, 'Name', [original_figure_title ' - Happy Holidays!'],...
                'Renderer', 'opengl',...
                'UserData','sn0w1ng',...
                'KeyReleaseFcn', @figureKeyUp_Callback, ...
                'CloseRequestFcn', @figureCloseReq_Callback);

% The axes for the snowing effect
snow_axes_hdl = axes('Parent', figure_hdl, ...
                     'Units', 'normalized', ...
                     'Position', [0 0 1 1], ...
                     'Color', [0 0 0], ...
                     'XLim', [0 1], ...
                     'YLim', [0 1], ...
                     'NextPlot','add', ...
                     'Visible', 'off', ...
                     'XTick', [], ...
                     'YTick', []);

% The axes for the frost effect (not implemented for this version)
frost_axes_hdl = axes('Parent', figure_hdl, ...
                     'Units', 'normalized', ...
                     'Position', [0 0 1 1], ...
                     'Color', [0 0 0], ...
                     'XLim', [0 1], ...
                     'YLim', [0 1], ...
                     'NextPlot','add', ...
                     'Visible', 'off', ...
                     'XTick', [], ...
                     'YTick', []);

%% Create the image objects
snow_plot_hdl = zeros(1,4);

% One snowflake is virtually four different markers overlapped
snow_plot_hdl(1) = plot(snowflakes_x, snowflakes_y, 'h', 'MarkerSize',snowflakes_size);
snow_plot_hdl(2) = plot(snowflakes_x, snowflakes_y, 'x', 'MarkerSize',snowflakes_size);
snow_plot_hdl(3) = plot(snowflakes_x, snowflakes_y, '+', 'MarkerSize',snowflakes_size.*3/4);
snow_plot_hdl(4) = plot(snowflakes_x, snowflakes_y, 's', 'MarkerSize',snowflakes_size.*3/5);

set(snow_plot_hdl, 'Color', snowflakes_color);
% set(snow_plot_hdl, 'Color', snowflakes_color, 'XDataSource', 'snowflakes_x', 'YDataSource', 'snowflakes_y');

% frost_image_hdl = image([0 1], [0 1], frost_image_data, ...
%                         'Parent', frost_axes_hdl,...
%                         'AlphaData', frost_alpha_data, ...
%                         'Visible', 'on');
% frost_alpha_data(150, 200) = 1;
start_time = tic;
last_time = toc(start_time);
uistack(snow_axes_hdl, 'top');
uistack(frost_axes_hdl, 'top');

% Start the falling snow!
start(snow_timer);

%% TimerFcn Callback function for updating the falling snow
    function updateSnow_Callback_nest(obj, event)
        current_time = toc(start_time);
        elapsed_time = current_time - last_time;

        % Update the position
        snowflakes_pos = snowflakes_pos + snowflakes_vel.*elapsed_time;

        % Find the snowflakes falling on the ground, and let them reappear
        % from the top of the sky.
        ind_snowflakes_outofbound = snowflakes_pos(2,:)<-0.05;
        snowflakes_pos(2, ind_snowflakes_outofbound) = snowflakes_pos(2, ind_snowflakes_outofbound) + 1.1;

        % Update the x/y data source
        snowflakes_x = snowflakes_pos(1,:);
        snowflakes_y = snowflakes_pos(2,:);
        %         refreshdata(snow_plot_hdl, 'caller'); % BROKEN IN R2012b
        set(snow_plot_hdl, 'XData', snowflakes_x, 'YData', snowflakes_y);
        drawnow;
        last_time = current_time;

    end
%% KeyReleaseFcn: stop the snow when 's' is pressed
    function figureKeyUp_Callback(obj, event)
        key = get(figure_hdl, 'CurrentKey');
        if strcmp(key, 's')
            stopsnow;
            set(figure_hdl, 'KeyReleaseFcn', []);
        end
    end

%% CloseRequestFcn: stop the snow when the figure window is closed
    function figureCloseReq_Callback(obj, event)
        if strcmp(get(figure_hdl, 'UserData'), 'sn0w1ng')
            stopsnow;
        end
        closereq;
    end

%% Clean up function
    function stopsnow()
        if isvalid(snow_timer)
            stop(snow_timer);
            delete(snow_timer);
            delete(snow_axes_hdl);
            delete(frost_axes_hdl);
        end

        set(figure_hdl, 'UserData', []);
        set(figure_hdl, 'Name', original_figure_title);
    end
end