text:
%%
% calculates the loglikelihood scores for random simulations
%%

global no1 bo1 inb1 inb2

compflag = 'on';
tmp_newt2 = newt2;

%%
% get the name of the p value grid file and load it
%%

prompt = {'Enter the name of the random simulation file',...
    'Enter the name of the real data file (p value estimation)',...
    'Enter the forecast start time (days from mainshock)',...
    'Enter the forecast end time (days from mainshock)'...
    'Enter the number of random simulations'};
title = 'Filename Input';
lines = 1;
def = {'','','10','15','100'};
answer = inputdlg(prompt,title,lines,def);
rand_in = answer{1,1};
real_in = answer{2,1};
fore_sday = str2double(answer{3,1});
fore_eday = str2double(answer{4,1});
numrand = str2double(answer{5,1});

real_in = 'test.mat';
load(real_in);

%%
% get a b c p and k from input file (in that order) for spatially varying models
%%
tvg(:,1) = bpvg(:,16);
tvg(:,2) = bpvg(:,18);
tvg(:,3) = bpvg(:,14);
tvg(:,4) = bpvg(:,11);
tvg(:,5) = bpvg(:,3);
tvg(:,6) = bpvg(:,4);

rand_in = '100rand.mat';
load(rand_in)
nodes = tmpgri;

%bpvg = rpvg;

%if(bpvg(1,15) == 1)
%         % this grid was not calculated with constant area so stop.
%         disp('%%%%ERROR%%%% This input grid was NOT calculated with constant radius please use one that is.')
%         return
%end

%%
% cut cat to only forecast time period
%%
fore_start = maepi(:,3) + fore_sday/365;
fore_end = maepi(:,3) + fore_eday/365;
ll = tmp_newt2.Date >= fore_start & tmp_newt2.Date <= fore_end;
fore_cat = tmp_newt2(ll,:);

dt = .1
calc_mag = [];
calc_cummag = [];

%%
% call distofault to get the distance from each grid node to the estimated fault
% gx and gy come from the input file
%%

faultdist

lt1 = fdkm < 1;
fdkm(lt1) = 1;

disp('values for sequence specific model');
pvalcat
ssa = rja;
ssb = rjb;
ssc = cv;
ssp = pv;

%%
% use the distance to taper the a value for the Generic calculation
% and the sequence specific calculation
%%

%totd = sum(1./(fdkm.^2));
totd = sum(fdkm.^2);
da = -1.67-log(totd);
dssa = ssa/totd;
numevents = length(fore_cat);
for i = 1:length(fdkm)
    id(i) = totd/(fdkm(i)^2);
end
sid =  sum(id);
sler = id/sid;
%%
% gevents is the number of events scaled with distance from fault
% this is not actually observed but is only the number of total
% obs events distributed over the entire grid with number of events
% taper as 1/r^2
%
% a is scaled by subtracting the log of the ratio of number of
% events at the node to the entire number of events from the
% overall a value.
%%
gevents = sler*numevents;
for i = 1:length(gevents)
    tapera(i) = -1.67-log(numevents/gevents(i));
    sstapera(i) = ssa-log(numevents/gevents(i));
end

%plot_tapera(tapera,s1x,s1y,s2x,s2y,tmpgri,xvect,yvect)
magco_fixed = 1.0
lvary_a = [];
lvary_ab = [];
lgca = [];

%%
% calculate the overall b-value
%%
[magml bo1 stanml,  avml] =  bmemag(newt2);

%%
% loop over all grid nodes and calculate various forecasts
%%

allcount = 0;
itotal = length(bpvg);
wait = waitbar(0,'Thinking...');
[lgl dum,  dum] = size(bpvg);
for gloop = 1:itotal
    allcount=allcount+1;

    %%
    %     Calculate the forecast for varying a & b and varying a w/constant b
    %
    %  select the events within the radius used for the parameter calc. (rd) from
    %  the forecast time period, to compare the calculated number of events with
    %%
    xg = bpvg(gloop,3);
    yg = bpvg(gloop,4);
    l = sqrt(((fore_cat(:,1)-xg)*cos(pi/180*yg)*111).^2 + ((fore_cat(:,2)-yg)*111).^2) ;
    [s,is] = sort(l);
    fore_cat = fore_cat(is(:,1),:) ;       % re-orders matrix to agree row-wise

    %%
    % get fore events  w/in the original radius
    %%
    l3 = l <= bpvg(1,5);
    obs_events = fore_cat(l3,:);      % obs_events is the events at the node w/in the radius of the original

    for rand_loop = 1:numrand+1
        if rand_loop == 1  % only calculate these the first time thru -- not for each rand loop

            %%
            % Calculate forecast for:
            % Sequence Specific model with a-values smoothed as 1/r^2 from fault and calculate the PDF
            %%

            mag_events = [];
            ct=1;
            for m = magco_fixed-0.1:0.1:7
                tct = 1;
                for t = fore_sday:dt:fore_eday
                    mag_events_ss(tct,ct) = 10^(sstapera(gloop)+ssb*(maepi(1,6)-m))*(t+ssc).^(-ssp);
                    tct = tct + 1;
                end
                ct = ct + 1;
            end

            %%
            % calculate the PDF and the log likelihood score for the ss model -- a tapered away from fault
            %%
            sum_mag = sum(mag_events_ss,1)';
            diff_mag = -diff(sum_mag,1);
            ss_pdf = poisspdf(obs_mag',diff_mag);

            lss(:,gloop) = log(ss_pdf);
            l = isinf(lss(:,gloop));
            lss(l,gloop) = -350;

            %%
            % get a b c p and k from input file (in that order) for spaitally varying models
            %%
            x(1) = tvg(gloop,1);
            x(2) = tvg(gloop,2);
            x(3) = tvg(gloop,3);
            x(4) = tvg(gloop,4);

        else %% if rand_loop ~= 1
            x(1) = rpvg(gloop, 1, rand_loop-1);
            x(2) = rpvg(gloop, 2, rand_loop-1);
            x(4) = rpvg(gloop, 3, rand_loop-1);
            x(3) = rpvg(gloop, 4, rand_loop-1);
        end

        %%
        % Calculate forecast for:
        % Generic CA model with a-values smoothed as 1/r^2 from fault and calculate the PDF
        %%

        num_nodes = length(nodes);
        gen_b = .91;
        gen_p =  1.08;
        gen_c = .05;

        mag_events = [];
        ct=1;
        for m = magco_fixed-0.1:0.1:7
            tct = 1;
            for t = fore_sday:dt:fore_eday
                mag_events_gca(tct,ct) = 10^(tapera(gloop)+gen_b*(maepi(1,6)-m))*(t+gen_c).^(-gen_p);
                tct = tct + 1;
            end
            ct = ct + 1;
        end

        %%
        % calculate the PDF and the log likelihood score for Generic model -- a tapered away from fault
        %%
        sum_mag = sum(mag_events_gca,1)';
        diff_mag = -diff(sum_mag,1);
        [obs_mag,gca_magbin] = hist(obs_events(:,6),magco_fixed-0.05:0.1:7);
        gca_pdf = poisspdf(obs_mag',diff_mag);

        lgca(:,gloop) = log(gca_pdf);
        l = isinf(lgca(:,gloop));
        lgca(l,gloop) = -350;

        %%
        % calculate the forecast for each magnitude bin for varying a and constant b
        % bo1 is the constant b value calculated above using max likelihood
        %%
        mag_events = [];
        ct=1;
        for m = magco_fixed-0.1:0.1:7
            tct = 1;
            for t = fore_sday:dt:fore_eday
                mag_events_ab(tct,ct) = 10^(x(1)+x(2)*(maepi(1,6)-m))*(t+x(3)).^(-x(4));
                mag_events_a(tct,ct) = 10^(x(1)+bo1*(maepi(1,6)-m))*(t+x(3)).^(-x(4));
                tct = tct + 1;
            end
            ct = ct + 1;
        end
        %%
        % calculate the PDF and the log likelihood score for varying a & b
        %%
        sum_mag = sum(mag_events_ab,1)';
        diff_mag = -diff(sum_mag,1);
        [obs_mag,magbin] = hist(obs_events(:,6),magco_fixed-0.05:0.1:7);
        vary_abpdf = poisspdf(obs_mag',diff_mag);

        lvary_ab(:,gloop,rand_loop) = log(vary_abpdf);
        l = isinf(lvary_ab(:,gloop,rand_loop));
        lvary_ab(l,gloop,rand_loop) = -350;

        %%
        % calculate the PDF and the log likelihood score for varying a & constant b
        %%
        sum_mag = sum(mag_events_a,1)';
        diff_mag = -diff(sum_mag,1);
        [obs_mag,vary_amagbin] = hist(obs_events(:,6),magco_fixed-0.05:0.1:7);
        vary_apdf = poisspdf(obs_mag',diff_mag);

        lvary_a(:,gloop,rand_loop) = log(vary_apdf);
        l = isinf(lvary_a(:,gloop,rand_loop));
        lvary_a(l,gloop,rand_loop) = -350;

        waitbar(allcount/itotal);
    end
end

close(wait);

slvary_ab=[];
slvary_a=[];
gdiffab=[];
gdiffa=[];

%%
% sum the scores
%%
slgca = sum(lgca(:,gloop),2);
slss = sum(lss(:,gloop),2);
[dum numm,  dum] = size(lvary_ab);
for rloop2 = 1:numrand+1
    %%
    % sum the scores at each node
    %%
    slgca(rloop2,1:numm) = sum(lgca(:,:,rloop2),1);
    slvary_ab(rloop2,1:numm) = sum(lvary_ab(:,:,rloop2),1);
    slvary_a(rloop2,1:numm) = sum(lvary_a(:,:,rloop2),1);
    if rloop2 >= 2
        gdiffab(rloop2,1:length(slvary_ab(1,:))) = slvary_ab(rloop2,:) - slgca.subset(rloop2);
        gdiffa(rloop2,1:length(slvary_a(1,:))) = slvary_a(rloop2,:) - slgca.subset(rloop2);
    end
end

meanab = mean(gdiffab,2);
meana = mean(gdiffa,2);

[xll, yll] =size(slvary_ab);
gcatot = sum(slgca);
sstot = sum(slss);
abtot = sum(slvary_ab(2:numrand+1,:));
atot = sum(slvary_a(2:numrand+1,:));

rdabtot = sum(slvary_ab(1,:));
rdatot = sum(slvary_a(1,:));

tdiff = gcatot - abtot;

compflag = 'of';

author: CelsoReyes


text:
clear all, close all, clc
n = 256;
w = exp(-i*2*pi/n);

% Slow
for i=1:n
    for j=1:n
        DFT(i,j) = w^((i-1)*(j-1));
    end
end

% Fast
[I,J] = meshgrid(1:n,1:n);
DFT = w.^((I-1).*(J-1));
imagesc(real(DFT))

author: dynamicslab


text:
% SRH is a robust pitch tracker.
%
% Octave compatible
%
% Description
%  The Summation of the Residual Harmonics (SRH) method is described in [1].
%  This algorithm exploits a criterion taking into the strength of the
%  harmonics and subharmonics of the residual excitation signal in order to
%  determine both voicing decision and F0 estimates. It is shown in [1] to
%  be particularly interesting in adverse conditions (low SNRs with various
%  types of additive noises).
%
%
% Inputs
%  wave            : [samples] [Nx1] input signal (speech signal)
%  fs              : [Hz]      [1x1] sampling frequency
%  f0min           : [Hz]      [1x1] minimum possible F0 value
%  f0max           : [Hz]      [1x1] maximum possible F0 value
%  hopsize         : [ms]      [1x1] time interval between two consecutive
%                    frames (i.e. defines the rate of feature extraction).
%
% Outputs
%  F0s             : vector containing the F0 estimates (values are
%                    provided even in unvoiced parts).
%  VUVDecisions    : vector containing the binary voicing decisions.
%  SRHVal          : vector containing the SRH values (according the
%                    harmonic criterion - voicing decision are derived from
%                    these values by simple thresholding).
%  time            : [s] Analysis instants of the features described above.
%
% Example
%  Please see the HOWTO_glottalsource.m example file.
%
% References
%  [1] T.Drugman, A.Alwan, "Joint Robust Voicing Detection and Pitch Estimation
%      Based on Residual Harmonics", Interspeech11, Firenze, Italy, 2011.
%      Publication available at the following link:
%      http://tcts.fpms.ac.be/~drugman/files/IS11-Pitch.pdf
%
% Copyright (c) 2011 University of Mons, FNRS
%
% License
%  This code is a part of the GLOAT toolbox with the following
%  licence:
%  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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Thomas Drugman thomas.drugman@umons.ac.be
%
% Modified
%  John Kane kanejo@tcd.ie September 27th 2014 - Bug fix and efficiency

function [F0s,VUVDecisions,SRHVal,time] = pitch_srh(wave,fs,f0min,f0max,hopsize)

if length(wave)/fs<0.1
    display('SRH error: the duration of your file should be at least 100ms long');
    pause(0.001)
    F0s=0;VUVDecisions=0;SRHVal=0;time=0;
end

if f0max<=f0min
    display('You look funny! Your f0min should be lower than f0max!!')
    pause(0.001)
end

if fs>16000
    display('Sample rate is bigger than 16kHz. Audio is resampled.')
    wave=resample(wave,16000,fs);
    fs=16000;
end

if nargin < 5
    hopsize=10;
end

%% Setings
nHarmonics = 5;
SRHiterThresh = 0.1;
SRHstdThresh = 0.05;
VoicingThresh = 0.07;
VoicingThresh2 = 0.085;
LPCorder=round(3/4*fs/1000);
Niter=2;

%% Compute LP residual
[res] = lpcresidual(wave,round(25/1000*fs),round(5/1000*fs), ...
                    LPCorder);

%% Create frame matrix
waveLen = length(wave);
clear wave;
frameDuration = round(100/1000*fs)-2; % Minus 2 to make equivalent
                                      % to original
frameDuration = round(frameDuration/2)*2; % Enforce evenness of the frame's length
shift = round(hopsize/1000*fs);
halfDur = round(frameDuration/2);
time = halfDur+1:shift:waveLen-halfDur;
N = length(time);
frameMat=zeros(frameDuration,N);
for n=1:N
    frameMat(:,n) = res(time(n)-halfDur:time(n)+halfDur-1);
end
clear res;

%% Create window matrix and apply to frames
win = blackman(frameDuration);
frameMat = bsxfun(@times, frameMat, win);

%% Do mean subtraction
frameMat = bsxfun(@minus, frameMat, mean(frameMat, 1));

%% Compute spectrogram matrix
specMat = zeros(floor(fs/2), size(frameMat, 2));
idx = 1:floor(fs/2);
for i = 1:size(frameMat,2)
    tmp = abs(fft(frameMat(:, i), fs));
    specMat(:, i) = tmp(idx);
end
clear frameMat tmp idx;
specMat = bsxfun(@rdivide, specMat, sqrt(sum(specMat.^2, 1)));

%% Estimate the pitch track in 2 iterations
no_pitch_range_adjustment = true;
for Iter=1:Niter

    [F0s,SRHVal] = SRH( specMat, nHarmonics, f0min, f0max );

    if max(SRHVal) > SRHiterThresh
        F0medEst = median( F0s( SRHVal > SRHiterThresh ) );

        % Only refine F0 limits if within the original limits
        if round(0.5*F0medEst) > f0min
            f0min=round(0.5*F0medEst);
            no_pitch_range_adjustment = false;
        end
        if round(2*F0medEst) < f0max
            f0max=round(2*F0medEst);
            no_pitch_range_adjustment = false;
        end
    end
    if no_pitch_range_adjustment
        break
    end
end
clear specMat;
time=time/fs;

%% Voiced-Unvoiced decisions are derived from the value of SRH (Summation of
%% Residual Harmonics)
VUVDecisions = zeros(1,N);

if std(SRHVal) > SRHstdThresh
   VoicingThresh = VoicingThresh2;
end

VUVDecisions( SRHVal > VoicingThresh ) = 1;

end

function [F0,SRHVal] = SRH( specMat, nHarmonics, f0min, f0max )

% Function to compute Summation of Residual harmonics function
% on a spectrogram matrix, with each column corresponding to one
% spectrum.

% Initial settings
N = size(specMat, 2);
SRHmat = zeros(f0max,N);

fSeq = f0min:f0max;
fLen = length(fSeq);

% Prepare harmonic indeces matrices.
plusIdx = repmat( (1:nHarmonics)',1,fLen) .* repmat(fSeq,nHarmonics,1);
subtrIdx = round( repmat( (1:nHarmonics-1)'+.5,1,fLen) .* ...
                  repmat(fSeq,nHarmonics-1,1) );

% avoid costly repmat operation by adjusting indices
plusIdx = mod(plusIdx-1,size(specMat,1))+1;
subtrIdx = mod(subtrIdx-1,size(specMat,1))+1;

% Do harmonic summation
for n=1:N
    specMatCur = specMat(:,n);
    SRHmat(fSeq,n) = ( sum( specMatCur(plusIdx), 1 ) - ...
                       sum( specMatCur(subtrIdx), 1 ) )';
end

% Retrieve f0 and SRH value
[SRHVal,F0] = max( SRHmat );

end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function disp(A)
%DISP   Display a CHEBOP by converting it to a pretty-printed string.
%
% See also DISPLAY.

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

loose = strcmp(get(0, 'FormatSpacing'), 'loose');

if ( isempty(A) )
    fprintf('   (empty chebop)\n');
    if ( loose )
        fprintf('\n')
    end

else
    % Store the name of the variable(s) to the operator
    args = [];

    % Check whether the chebop is linear in order to be able to display
    % linearity information.
    if ( isempty(A.op) )
        fprintf('   Empty operator');
    elseif ( islinear(A) )
        fprintf('   Linear operator');
    else
        fprintf('   Nonlinear operator');
    end
    if ( ~isempty(A.op) )
        fprintf(':\n')
        [str, args] = formatOperator(A.op);
        if ( ~isempty(A.opShow) )
            str = A.opShow;
        end
        fprintf('      %s\n', str);
        if ( loose )
            fprintf('\n')
        end
        fprintf('  ');
    end
    fprintf(' operating on chebfun objects defined on:\n')
    Adom = A.domain;
    dom = strtrim(sprintf('%g ', Adom));
    % Replace ' ' with commas in the domain for printing. If the numbers are not
    % all integers, we also want to keep whitespace between the numbers.
    if ( all(floor(Adom) == ceil(Adom)) )
        dom = strrep(dom, ' ', ',');
    else
        dom = strrep(dom, ' ', ', ');
    end
    % Print the domain
    fprintf('      [%s]\n', dom);
    if ( loose )
        fprintf('\n')
    end

    if ( ~isempty(A.lbc) || ~isempty(A.rbc) || ~isempty(A.bc) )
        fprintf('   with\n');
        if ( loose )
            fprintf('\n')
        end
    end

    if ( ~isempty(A.lbc) )
        fprintf('    left boundary condition(s):\n');
        printLRbc(A.lbc, A.lbcShow, args, loose);
    end

    if ( ~isempty(A.rbc) )
        fprintf('    right boundary condition(s):\n');
        printLRbc(A.rbc, A.rbcShow, args, loose);
    end

    if ( ~isempty(A.bc) )
        if ( strcmpi(A.bc, 'periodic') )
            fprintf('    periodic boundary conditions.\n');
        else
            fprintf('    constraints:\n');
            if ( ischar(A.lbc) )
                fprintf('      %s\n', A.bc);
            else
                str = stripHandle(func2str(A.bc));
                fprintf('      %s = 0\n', str);
            end
        end
        if ( loose )
            fprintf('\n')
        end
    end

    if ( ~isempty(A.maxnorm) )
        mn = A.maxnorm;
        mn = mn(:)'; % Ensure row vector for printing
        if ( length(mn) > 1 )
            maxnormString = ['[' num2str(mn) ']'];
        else
            maxnormString = num2str(mn);
        end
        fprintf('   enforcing a maximum norm of solution: %s\n', maxnormString);
    end

end

end  % End of function DISP().

function str = stripHandle(op)
idx = strfind(op, ')');
str = op(idx(1)+1:end);

end

function printLRbc(bc, bcShow, args, loose)
% Print left or right boundary condition(s)
if ( ischar(bcShow) )
    % We had an assignment such as N.lbc = 'neumann';
    fprintf('      %s\n', bcShow);
elseif ( isnumeric(bcShow) )
    % We had an assignment such as N.lbc = [2; 3];
    if ( ~isempty(args) )
        % How many conditions are we dealing with?
        numBC = length(bcShow);

        % Check if we're dealing with a scalar case (e.g. u=2, u'=3) or system
        % case (e.g. u=2, v=3):
        if ( isempty(strfind(args, ',')) )  % Scalar case

            % Extra whitespace to make things align
            if ( numBC < 4)
                extraWS = repmat(' ', 1, numBC - 1);
            else
                extraWS = '    ';
            end

            % We always need to print the first condition:
            fprintf('      %s%s = %s\n', args, extraWS, num2str(bcShow(1)));
            % If we got more conditions, print them as well:
            if ( numBC >= 2 )
                fprintf('      %s''%s = %s\n', args, ...
                    extraWS(1:end-1), num2str(bcShow(2)));
            end
            if ( numBC >= 3 )
                fprintf('      %s''''%s = %s\n', args, ...
                    extraWS(1:end-2), num2str(bcShow(3)));
            end
            % Print all remaining conditions:
            for bcCounter = 4:numBC
                fprintf('      %s^(%i) = %s\n', args, ...
                    bcCounter-1, num2str(bcShow(bcCounter)));
            end
        else
           % Vector case. Split args back into individual variables before so we
           % can print them one per line
           args = strsplit(args,',');

           for bcCounter = 1:numBC
               fprintf('      %s = %s\n', args{bcCounter}, num2str(bcShow(bcCounter)));
           end
        end
    else
        fprintf('      %s\n', num2str(bcShow));
    end
else
    str = stripHandle(func2str(bc));
    fprintf('      %s = 0\n', str);
end
if ( loose )
    fprintf('\n')
end
end

function [str, args] = formatOperator(op)
% What type of a function did we get passed?
opFunc = functions(op);

% Convert the function to a string:
op = func2str(op);

if ( strcmp(opFunc.type, 'simple') )
    % If the function passed was of the simple type, we simply return the
    % results of func2str above
    str = op;

else
    % We got passed an anonymous function. Return a string to print on a nice
    % form.
    idx = strfind(op, ')');
    args = op(3:idx(1)-1);
    % If we have any commas in args, the independent variable (e.g. x) must be
    % included. For nicer output, don't print it
    idxCommas = strfind(args, ',');
    if ( ~isempty(idxCommas) )
        args = args(idxCommas(1)+1:end);
    end

    % Isolate the actual operator part from the string
    op = op(idx(1)+1:end);

    % If we're dealing with a system, throw away the end [ and ]. Also, split up
    % components into lines.
    if ( op(1) == '[' && op(end) == ']' )
        op = op(2:end-1);

        % Replace ; with the ASCII character for new line. Add indentation to
        % make it look nice
        wSpace = repmat(' ', 1, length(args) + 12);
        op = strrep(op, ';', [sprintf('\n'), wSpace]);
    end

    % Output string
    str = [args, ' |--> ', op];
end

end

function s = mat2char(V)
% MAT2CHAR  Convert the varmat to pretty-print string (incl. BCs)

print_bc = true;       % Print BCs or not
defreal = 6;           % Size of matrix to display
dom = get(V,'domain');
numints = numel(dom.endsandbreaks);

if numints == 1 && isempty(V,'bcs')
    print_bc = false;
end

defreal = max(2,ceil(defreal/(numints-1)));
s1 = ['   with n = ',int2str(defreal),' realization:'];

% Evaluate the linop
if print_bc && ~isempty(V,'bcs')
    Vmat = feval(V,defreal,'bc',dom);
else
    Vmat = feval(V,defreal,'nobc',dom);
end

% Some pretty scaling for large/small numbers
M = max(max(abs(Vmat)));
if ( (M > 1e2 || M < 1) && M ~= 0 )
    M = 10^round(log10(M)-1);
    Vmat = Vmat/M;
    s1 = [s1 sprintf('\n    %2.1e * ',M)];
end

if isreal(Vmat)
    s2 = num2str(Vmat,'  %8.4f');
    for k = 1:size(s2,1)
        s2(k,:) = strrep(s2(k,:),' 0.0000','      0');
        if numel(s2(k,:)) > 5,
            s2(k,1:6) = strrep(s2(k,1:6),'0.0000','     0');
        end
        s2(k,:) = strrep(s2(k,:),'-0.0000','      0');
    end
else
    s2 = num2str(Vmat,'%8.2f');
    for k = 1:size(s2,1)
        s2(k,:) = strrep(s2(k,:),'0.00 ','   0 ');
    end
end

% Pad with appropriate spaces
s2 = [ repmat(' ',size(s2,1),5) s2 ];
s = char(s1,'',s2);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% Copyright (C) 2014-2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym erfc (@var{x})
%% Symbolic erfc function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = erfc (x)
%%   @result{} y = (sym) erfc(x)
%% @end group
%% @end example
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @end defmethod

function y = erfc(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('erfc', x);
end

%!error erfc (sym(1), 2)
%!assert (isequaln (erfc (sym(nan)), sym(nan)))

%!shared x, d
%! d = 1;
%! x = sym('1');

%!test
%! f1 = erfc(x);
%! f2 = erfc(d);
%! assert( abs(double(f1) - f2) < 1e-15 )

%!test
%! D = [d d; d d];
%! A = [x x; x x];
%! f1 = erfc(A);
%! f2 = erfc(D);
%! assert( all(all( abs(double(f1) - f2) < 1e-15 )))

%!test
%! % round trip
%! y = sym('y');
%! A = erfc (d);
%! f = erfc (y);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B, -eps)

repo: octsympy


text:
function lifetimes = getStateLifeTimes (Gamma,T,options,threshold,threshold_Gamma,do_concat)
% Computes the state life times for the state time courses, in number of time points
%
% Gamma can be the probabilistic state time courses (time by states),
%   which can contain the probability of all states or a subset of them,
%   or the Viterbi path (time by 1).
% In the first case, threshold_Gamma is used to define when a state is active.
%
% The parameter threshold is used to discard visits that are too
%   short (deemed to be spurious); these are discarded when they are shorter than
%   'threshold' time points
%
% The parameter 'options' must be the same than the one supplied to the
% hmmmar function for training.
%
% lifetimes is then a cell (subjects/segments/trials by states), where each element is a
%   vector of life times
%
% if do_concat is specified, the lifetimes are concatenated across
% segments of data, so that lifetimes is a single cell (1 by states)
%
% Diego Vidaurre, OHBA, University of Oxford (2016)

raise_warning=false;
if nargin<3, options = struct(); options.Fs = 1; options.downsample = 1; end
if ~isfield(options,'Fs'), options.Fs = 1; end
if ~isfield(options,'downsample'), options.downsample = options.Fs; end
if ~isfield(options, 'dropstates'), options.dropstates = 0; end
if nargin<4 || isempty(threshold), threshold = 0; end
if nargin<5 || isempty(threshold_Gamma), raise_warning=true; threshold_Gamma = (2/3); end
if nargin<6 || isempty(do_concat), do_concat = true; end

is_vpath = (size(Gamma,2)==1 && all(rem(Gamma,1)==0));
if iscell(T)
    if size(T,1) == 1, T = T'; end
    for i = 1:length(T)
        if size(T{i},1)==1, T{i} = T{i}'; end
    end
    Nsubj = length(T);
    trials2subjects = zeros(length(cell2mat(T)),1); ii = 1;
    for i = 1:length(T)
        Ntrials = length(T{i});
        trials2subjects(ii:ii+Ntrials-1) = i;
        ii = ii + Ntrials;
    end
    T = cell2mat(T);
else
    Nsubj = length(T);
    trials2subjects = 1:Nsubj;
end
N = length(T);

r = 1;
if isfield(options,'downsample') && options.downsample>0
    r = (options.downsample/options.Fs);
end

if ~isfield(options,'order') && ~isfield(options,'embeddedlags')
   options.order = (sum(T) - size(Gamma,1)) / length(T);
end

if isfield(options,'tuda') && options.tuda | ...
        isfield(options, 'order') && options.order == 0
    T = ceil(r * T);
elseif isfield(options,'order') && options.order > 0
    T = ceil(r * T);
    T = T - options.order;
elseif isfield(options,'embeddedlags') && length(options.embeddedlags) > 1
    d1 = -min(0,options.embeddedlags(1));
    d2 = max(0,options.embeddedlags(end));
    T = T - (d1+d2);
    T = ceil(r * T);
end

if is_vpath % viterbi path
    vpath = Gamma;
    if options.dropstates==1
        K = length(unique(vpath));
    else
        K = options.K;
    end
    Gamma = zeros(length(vpath),K);
    for k = 1:K
       Gamma(vpath==k,k) = 1;
    end
else
    if raise_warning
        warning(['Using the Viterbi path is here recommended instead of the state ' ...
            'probabilistic time courses (Gamma)'])
    end
    K = size(Gamma,2);
    Gamma = Gamma > threshold_Gamma;
end

if ~do_concat
    lifetimes = cell(Nsubj,K);
    for j = 1:N
        t0 = sum(T(1:j-1));
        ind = (1:T(j)) + t0;
        if length(ind)==1, continue; end
        jj = trials2subjects(j);
        for k = 1:K
            lifetimes{jj,k} = [lifetimes{jj,k} aux_k(Gamma(ind,k),threshold)];
        end
    end
else
    lifetimes = cell(1,K);
    for j = 1:N
        t0 = sum(T(1:j-1));
        ind = (1:T(j)) + t0;
        if length(ind)==1, continue; end
        for k = 1:K
            lifetimes{k} = [lifetimes{k} aux_k(Gamma(ind,k),threshold)];
        end
    end
end

end

function lifetimes = aux_k(g,threshold)

lifetimes = [];

while ~isempty(g)

    t = find(g==1,1);
    if isempty(t), break; end
    lt = find(g(t+1:end)==0,1);

    if isempty(lt) % end of trial
        if (length(g)-t+1)>threshold
            lifetimes = [lifetimes (length(g)-t+1)];
        else
            g(t:end) = 0; % too short visit to consider
        end
        break;
    end

    if lt<threshold % too short visit, resetting g and trying again
        g(t:t+lt-1) = 0;
        continue;
    end

    lifetimes = [lifetimes lt];
    tnext = t + lt;
    g = g(tnext:end);

end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function [beta, t, pvals, convals, con_t, con_pvals, sigma, Phi, df, stebeta, conste, F] = fit_gls(y, X, c, p, varargin)
% Fit a linear model using generalized least squares and an AR(p) model
%
% :Usage:
% ::
%
%     [beta, t, pvals, convals, con_t, con_pvals, sigma, Phi, df, stebeta, conste, F] = fit_gls(y,X,c,p,[PX, equal to pinv(X), for speed], [Weights])
%
% This program uses the Cochrane-Orcutt algorithm to iteratively find the
% GLS solution and estimate the noise parameters.
%
% Step 1: Find the OLS solution.
%
% Step 2: Use residuals from the previous fit to estimate the parameters in
% the AR(p) noise model.
%
% Step 3: Find the GLS solution using the covariance matrix corresponding
% to an AR(p) model with parameters estimated in Step 2 inserted.
%
% Step 4: Repeat steps 2-3 until convergence.
%
% :Inputs:
%
%   **y:**
%        fMRI time course (T x 1 vector)
%
%   **X:**
%        Design matrix (T x param matrix)
%
%   **c:**
%        contrast vector(s) (param x # contrasts matrix)
%
%   **p:**
%        order of AR model.
%
%   **PX:**
%        pinv(X), for speeded, repeated calculations with different y vectors
%
%        Note: if using weights, px = inv(X' * W * X) * X' * W;
%              where W = diag(Weights);
%
%   **Weights:**
%        Optional, vector of weights for each observation
%
%        Empty or missing: Unweighted analysis.
%
% Note that setting p=0 implies a white noise model.
%
% :Output:
%
%   **t:**
%        t-value for the contrast c'beta
%
%   **df:**
%        degrees of freedom using Satterthwaite approximation
%
%   **beta:**
%        beta vector
%
%   **Phi:**
%        vector of coefficients in AR(p) model
%
%   **sigma:**
%        standard deviation
%
%   **stebeta:**
%        standard error of betas
%
% ..
%    by Martin Lindquist
%
%    Last updated: 3/29/08, Tor Wager, added weighted least squares
%                        verified that beta and t-values are identical to
%                        glmfit.m in matlab7.5 with ar p = 0
%                        ***AR(p) with weighted least squares needs to be
%                        checked.  behaving reasonably.
%               4/1/08,  Tor : output stats for boht betas and contrasts
%                        Reorganized order of outputs
%               5/15/08  Tor : Weird things happening with single inputs; particulaly with aryule; force
%                         double
% ..

y = double(y);
X = double(X);

T = length(y);      % Length of time course

k = size(X, 2);     % predictors

if length(varargin) > 1
    w = double(varargin{2});        % weights for weighted least squares
else
    w = ones(T, 1);
end
W = diag(w);
sqrtW = sqrt(W);     % for weighted residuals

if ~isempty(varargin) && ~isempty(varargin{1})
    px = double(varargin{1});
else
    invxvx = inv(X' * W * X);   % we can re-use this later if p == 0
    px = invxvx * X' * W;
end

% Step 1: Find the OLS solution
beta = px*y;                                      % Beta values; weighted, if weights are used
resid = sqrtW * y - sqrtW * X * beta;             % Residuals (Weighted, if weights are used)
sigma = sqrt((1 / (T - k)) * resid' * resid);     %sum(resid.^2)));  % Estimate of Sigma

% Weighted residuals: Three equivalent ways
% We pick one that is compatible with AR estimation
% 1)
% beta = px*y;                       % Beta values
% resid = y - X*beta;                     % Residuals
% sigma = sqrt((1 / (T - k)) * resid' * W * resid); %sum(resid.^2)));  % Estimate of Sigma

% 2)
%r2 = sqrt(W) * resid; sigma2 = sqrt((1 / (T - k)) * r2' *  r2)

% 3)
% r2 = sqrt(W) * y - sqrt(W) * X * beta;
% sigma2 = sqrt((1 / (T - k)) * r2' *  r2)

% Stuff needed for future iterations
iV = W;  % for ar p = 0 case
A = W;   % for ar p = 0 case
Phi = 0;
betaold = 0;

% Steps 2-4: Find the GLS solution by iteration
% If p=0, skip this step. Appropriate solution already calculated above.
% Continue iteration until convergence or for at most 10 loops.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note to Jack and Tor: Keith Worsley uses a similar algorithm when fitting
% a GLM with an AR(p) noise model. However, he skips the iterative step
% and only goes through the loop one time. He claims that this is enough. I
% am not entirely convinced, therefore it is probably better to go through
% a few times if needed.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

i=1;                % Set counter

while i < 10 && p > 0 && (i == 1 || sum((beta - betaold).^2) > 0.001)
    % Do up to 10 iterations, if arp > 0 and either first iteration or
    % there's a difference from last iteration

    % resid = y - X*beta;             % Calculate residuals of current model fit

    resid = sqrtW * y - sqrtW * X * beta;             % Residuals (Weighted, if weights are used)

    % Estimate AR parameters using residuals
    [a,e] = aryule(resid, p);
    Phi =zeros(length(a)-1,1);
    Phi(1:p) = -a(2:(p+1));

    %sigma = sqrt(e); % swrtW*e??  % Moved later, because never used until
    %after iteration loop

    % Find the inverse of the covariance matrix
    A = sqrtW; % ***should be sqrt(W)?

    for j=1:p
        %A = A + diag(-Phi(j)*ones((T-j),1),-j);

        A = A + sqrtW * diag(-Phi(j)*ones((T-j),1),-j);
    end;

    %create_figure('A'); imagesc(A, [-.2 .2]); colorbar, drawnow, input(' ')

    iV = A*A';                      %  New weights, with AR estimates, The inverse of the covariance matrix

    betaold = beta;                 % Set old solution to be betaold
    beta = inv(X'*iV*X)*X'*iV*y;    % Calculate new solution

    i = i+1;                        % Add one to counter

    %create_figure('COV'); imagesc(iV, [-.02 .02]); colorbar, drawnow, input(' ')
end

if p > 0
    % re-calc sigma
    sigma = sqrt(e);

    invxvx = inv(X'*iV*X);

    % Should we use the kind of thing below? Seems like sqrt(e) is unweighted,
    % though it seems reasonable...
    % sqrtiv = sqrt(iV);
    % resid = sqrtiv * y - sqrtiv * X * beta;             % Residuals (Weighted, if weights are used)
    % sigma = sqrt((1 / (T - k)) * resid' * resid)
end

R = (eye(T) - X * invxvx * X' * iV);        % Residual inducing matrix

Wd = R * A * inv(iV) * A';                  % inv(iV) = Covariance matrix
df = (trace(Wd).^2)./trace(Wd*Wd);       % Satterthwaite approximation for degrees of freedom

% Should check on below: this creates difference in t-values from
% glmfit...but why wouldn't we need to re-calculate a (larger) sigma if we have reduced df?
% if df ~= (T - k)                    % Tor added 3/29, have to re-calculate sigma
%     sigma = sqrt((1 / df) * resid' * iV * resid); %sum(resid.^2)));  % Estimate of Sigma
% end

varbeta = sigma^2.* invxvx;            % Var(beta)
stebeta = diag(varbeta).^.5;
t = beta ./ stebeta;

convals = [];
conste = [];
con_t = [];
con_pvals = [];
F = [];
if ~isempty(c)
    % Contrast(s)
    convals = (c' * beta);
    conste = diag(c' * varbeta * c).^.5;                   % tor added as output
    con_t = convals ./ conste;                            %sqrt(c'*varbeta*c);             % t-value
end

% get rid of nuisance regressors that aren't in any contrast
% wh = stebeta > 0;
% c = c(wh,wh);
% beta = beta(wh);
% stebeta = stebeta(wh);
%
% if ~isempty(c)
%     beta = (c' * beta);
%     %beta = beta(wh);
%     t = beta ./ stebeta;   %sqrt(c'*varbeta*c);            % t-value
% else
%     t = beta(wh) ./ stebeta;
% end

%%%%%% Test H0: beta_1 = beta_2 = .... = beta_param = 0
if nargout > 6
    SSE = y'*y - beta'*X'*y;                 % Error sum of squares
    mSSE = SSE/df;
    J = ones(T);
    SST=y'*y - (1/T).*y'*J*y;           % Total sum of squares
    SSM=SST-SSE;                        % Model sum of squares

    dfSSM=length(c) - 1;                % degrees of freedom for model (param - 1)
    mSSM=SSM/dfSSM;

    F=mSSM/mSSE;         % F-statistic - compare with F-distruibution with (param-1, df) degrees of freedom
end

% % get contrast values if we need those
% if ~isempty(c)
%     beta = (beta' * c)';
% end

if nargout > 2
    pvals = 2 .* (1 - tcdf(abs(t), df)); % two-tailed

    % make sure p-values for valid results are not zero...
    pvals(pvals == 0 & beta ~= 0 & ~isnan(beta)) = 1000*eps;

    if ~isempty(c) && nargout > 5
        con_pvals = 2 .* (1 - tcdf(abs(con_t), df)); % two-tailed

        con_pvals(con_pvals == 0 & convals ~= 0 & ~isnan(convals)) = 1000*eps;
    end

end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [reg_corner,rho,eta,reg_param] = l_curve(U,sm,b,method,L,V)
%L_CURVE Plot the L-curve and find its "corner".
%
% [reg_corner,rho,eta,reg_param] =
%                  l_curve(U,s,b,method)
%                  l_curve(U,sm,b,method)  ,  sm = [sigma,mu]
%                  l_curve(U,s,b,method,L,V)
%
% Plots the L-shaped curve of eta, the solution norm || x || or
% semi-norm || L x ||, as a function of rho, the residual norm
% || A x - b ||, for the following methods:
%    method = 'Tikh'  : Tikhonov regularization   (solid line )
%    method = 'tsvd'  : truncated SVD or GSVD     (o markers  )
%    method = 'dsvd'  : damped SVD or GSVD        (dotted line)
%    method = 'mtsvd' : modified TSVD             (x markers  )
% The corresponding reg. parameters are returned in reg_param.  If no
% method is specified then 'Tikh' is default.  For other methods use plot_lc.
%
% Note that 'Tikh', 'tsvd' and 'dsvd' require either U and s (standard-
% form regularization) or U and sm (general-form regularization), while
% 'mtvsd' requires U and s as well as L and V.
%
% If any output arguments are specified, then the corner of the L-curve
% is identified and the corresponding reg. parameter reg_corner is
% returned.  Use routine l_corner if an upper bound on eta is required.

% Reference: P. C. Hansen & D. P. O'Leary, "The use of the L-curve in
% the regularization of discrete ill-posed problems",  SIAM J. Sci.
% Comput. 14 (1993), pp. 1487-1503.

% Per Christian Hansen, IMM, July 26, 2007.

% Set defaults.
if (nargin==3), method='Tikh'; end  % Tikhonov reg. is default.
npoints = 200;  % Number of points on the L-curve for Tikh and dsvd.
smin_ratio = 16*eps;  % Smallest regularization parameter.

% Initialization.
[m,n] = size(U); [p,ps] = size(sm);
if (nargout > 0), locate = 1; else locate = 0; end
beta = U'*b; beta2 = norm(b)^2 - norm(beta)^2;
if (ps==1)
  s = sm; beta = beta(1:p);
else
  s = sm(p:-1:1,1)./sm(p:-1:1,2); beta = beta(p:-1:1);
end
xi = beta(1:p)./s;

if (strncmp(method,'Tikh',4) | strncmp(method,'tikh',4))

  eta = zeros(npoints,1); rho = eta; reg_param = eta; s2 = s.^2;
  reg_param(npoints) = max([s(p),s(1)*smin_ratio]);
  ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
  for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end
  for i=1:npoints
    f = s2./(s2 + reg_param(i)^2);
    eta(i) = norm(f.*xi);
    rho(i) = norm((1-f).*beta(1:p));
  end
  if (m > n & beta2 > 0), rho = sqrt(rho.^2 + beta2); end
  marker = '-'; txt = 'Tikh.';

elseif (strncmp(method,'tsvd',4) | strncmp(method,'tgsv',4))

  eta = zeros(p,1); rho = eta;
  eta(1) = abs(xi(1))^2;
  for k=2:p, eta(k) = eta(k-1) + abs(xi(k))^2; end
  eta = sqrt(eta);
  if (m > n)
    if (beta2 > 0), rho(p) = beta2; else rho(p) = eps^2; end
  else
    rho(p) = eps^2;
  end
  for k=p-1:-1:1, rho(k) = rho(k+1) + abs(beta(k+1))^2; end
  rho = sqrt(rho);
  reg_param = (1:p)'; marker = 'o';
  if (ps==1)
    U = U(:,1:p); txt = 'TSVD';
  else
    U = U(:,1:p); txt = 'TGSVD';
  end

elseif (strncmp(method,'dsvd',4) | strncmp(method,'dgsv',4))

  eta = zeros(npoints,1); rho = eta; reg_param = eta;
  reg_param(npoints) = max([s(p),s(1)*smin_ratio]);
  ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
  for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end
  for i=1:npoints
    f = s./(s + reg_param(i));
    eta(i) = norm(f.*xi);
    rho(i) = norm((1-f).*beta(1:p));
  end
  if (m > n & beta2 > 0), rho = sqrt(rho.^2 + beta2); end
  marker = ':';
  if (ps==1), txt = 'DSVD'; else txt = 'DGSVD'; end

elseif (strncmp(method,'mtsv',4))

  if (nargin~=6)
    error('The matrices L and V must also be specified')
  end
  [p,n] = size(L); rho = zeros(p,1); eta = rho;
  [Q,R] = qr(L*V(:,n:-1:n-p),0);
  for i=1:p
    k = n-p+i;
    Lxk = L*V(:,1:k)*xi(1:k);
    zk = R(1:n-k,1:n-k)\(Q(:,1:n-k)'*Lxk); zk = zk(n-k:-1:1);
    eta(i) = norm(Q(:,n-k+1:p)'*Lxk);
    if (i < p)
      rho(i) = norm(beta(k+1:n) + s(k+1:n).*zk);
    else
      rho(i) = eps;
    end
  end
  if (m > n & beta2 > 0), rho = sqrt(rho.^2 + beta2); end
  reg_param = (n-p+1:n)'; txt = 'MTSVD';
  U = U(:,reg_param); sm = sm(reg_param);
  marker = 'x'; ps = 2;  % General form regularization.

else
  error('Illegal method')
end

% Locate the "corner" of the L-curve, if required.
if (locate)
  [reg_corner,rho_c,eta_c] = l_corner(rho,eta,reg_param,U,sm,b,method);
end

% Make plot.
plot_lc(rho,eta,marker,ps,reg_param);
if locate
  ax = axis;
  HoldState = ishold; hold on;
  loglog([min(rho)/100,rho_c],[eta_c,eta_c],':r',...
         [rho_c,rho_c],[min(eta)/100,eta_c],':r')
  title(['L-curve, ',txt,' corner at ',num2str(reg_corner)]);
  axis(ax)
  if (~HoldState), hold off; end
end

repo: REST


text:
% another simple convection example using method of lines, still a bit
% messy, wait for a blog post
% Using Matlab ODE solvers to solve a convective flux that
% moves a plume in diagonal direction
% see how terribly diffusive the solution is!
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc
% define a 2D mesh
H = 1;
W = 1;
Nx = 100;
Ny = 100;
mesh1 = createMesh2D(Nx, Ny, W, H);
% velocity field
u = createCellVariable(mesh1, 0.001*ones(Nx,Ny));
uf = arithmeticMean(u);
% uf.yvalue(:,:)=0;
% diffusion field
D = 1e-2*createCellVariable(mesh1, 1e-2);
Df = arithmeticMean(D);
% transient term coefficient
alfa = createCellVariable(mesh1, 1);
% define the boundaries
% dirichlet on all the boundaries
BC = createBC(mesh1); % all Neumann
BC.bottom.a(:) = 0; BC.bottom.b(:) = 1; BC.bottom.c(:) = 0;
BC.left.a(:) = 0; BC.left.b(:) = 1; BC.left.c(:) = 0;
% BC.right.a(:) = 0; BC.right.b(:) = 1; BC.right.c(:) = 0;
% BC.top.a(:) = 0; BC.top.b(:) = 1; BC.top.c(:) = 0;
% Initial values
phi_old = createCellVariable(mesh1, 0, BC);
phi_old.value(5:10, 5:10) = 1;
% define the convection term
Mconv = convectionUpwindTerm(uf);
Mdif = diffusionTerm(Df);
% define the BC term
[Mbc, RHSbc] = boundaryCondition(BC);
% solver
final_t = 500;
% define the transient term
[M, RHS] = combineBC2D(BC, Mconv, ...
    zeros((Nx+2)*(Ny+2),1));
% eq: dcdt = D d2c/dx2
dcdt = @(t,c)(-M*c-RHS);
[t_temp, c_temp] = ode45(dcdt, [0 final_t], internalCells(phi_old));
for i =1:length(t_temp)
    phi.value = reshape(c_temp(i,:), Nx, Ny);
    figure(1);pcolor(phi.value');title(["t = ", num2str(t_temp(i))]); colorbar;drawnow;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function inurbs = nrbdegelev(nurbs, ntimes)
%
% Function Name:
%
%   nrbdegelev - Elevate the degree of the NURBS curve or surface.
%
% Calling Sequence:
%
%   ecrv = nrbdegelev(crv,utimes);
%   esrf = nrbdegelev(srf,{utimes,vtimes});
%
% Parameters:
%
%   crv		: NURBS curve, see nrbmak.
%
%   srf		: NURBS surface, see nrbmak.
%
%   utimes	: Increase the degree along U direction utimes.
%
%   vtimes	: Increase the degree along V direction vtimes.
%
%   ecrv	: new NURBS structure for a curve with degree elevated.
%
%   esrf        : new NURBS structure for a surface with degree elevated.
%
%
% Description:
%
%   Degree elevates the NURBS curve or surface. This function uses the
%   B-Spline function bspdegelev, which interface to an internal 'C'
%   routine.
%
% Examples:
%
%   Increase the NURBS surface twice along the V direction.
%   esrf = nrbdegelev(srf, 0, 1);
%
% See:
%
%   bspdegelev

%  D.M. Spink
%  Copyright (c) 2000.

if nargin < 2
  error('Input argument must include the NURBS and degree increment.');
end

if ~isstruct(nurbs)
  error('NURBS representation is not structure!');
end

if ~strcmp(nurbs.form,'B-NURBS')
  error('Not a recognised NURBS representation');
end

degree = nurbs.order-1;

if iscell(nurbs.knots)
  % NURBS represents a surface
  [dim,num1,num2] = size(nurbs.coefs);

  % Degree elevate along the v direction
  if ntimes(2) == 0
    coefs = nurbs.coefs;
    knots{2} = nurbs.knots{2};
  else
    coefs = reshape(nurbs.coefs,4*num1,num2);
    [coefs,knots{2}] = bspdegelev(degree(2),coefs,nurbs.knots{2},ntimes(2));
    num2 = size(coefs,2);
    coefs = reshape(coefs,[4 num1 num2]);
  end

  % Degree elevate along the u direction
  if ntimes(1) == 0
    knots{1} = nurbs.knots{1};
  else
    coefs = permute(coefs,[1 3 2]);
    coefs = reshape(coefs,4*num2,num1);
    [coefs,knots{1}] = bspdegelev(degree(1),coefs,nurbs.knots{1},ntimes(1));
    coefs = reshape(coefs,[4 num2 size(coefs,2)]);
    coefs = permute(coefs,[1 3 2]);
  end

else

  % NURBS represents a curve
  if isempty(ntimes)
    coefs = nurbs.coefs;
    knots = nurbs.knots;
  else
    [coefs,knots] = bspdegelev(degree,nurbs.coefs,nurbs.knots,ntimes);
  end

end

% construct new NURBS
inurbs = nrbmak(coefs,knots);

author: Sable

repo: mcbench-benchmarks


text:
function X = positive_definite_intrinsic_mean(A)
% Computes an intrinsic mean of a collection of positive definite matrices.
%
% function X = positive_definite_intrinsic_mean(A)
%
% Input:  A 3D matrix A of size nxnxm such that each slice A(:, :, k) is a
%         positive definite matrix of size nxn.
%
% Output: A positive definite matrix X of size nxn which is an intrinsic mean
%         of the m matrices in A, that is, X minimizes the sum of squared
%         Riemannian distances to the matrices in A:
%            f(X) = sum_k=1^m .5*dist^2(X, A(:, :, k))
%         The distance is defined by the natural metric on the set of
%         positive definite matrices: see sympositivedefinitefactory.
%
% This simple example is not the best way to compute intrinsic means. Its
% purpose it to serve as base code to explore other algorithms. In
% particular, in the presence of large noise, this algorithm seems not to
% be able to reach points with a very small gradient norm. This may be
% caused by insufficient accuracy in the gradient computation.
%
% See also: sympositivedefinitefactory

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Nicolas Boumal, Sept. 3, 2013
% Contributors:
%
% Change log:
%     Sep. 15, 2022 (NB):
%         Changed name from positive_definite_karcher_mean, clarified
%         some comments.

    % Generate some random data to test the function if none is given.
    if ~exist('A', 'var') || isempty(A)
        n = 5;
        m = 50;
        A = zeros(n, n, m);
        ref = diag(max(.1, 1+.1*randn(n, 1)));
        for i = 1 : m
            noise = 0.01*randn(n);
            noise = (noise + noise')/2;
            [V, D] = eig(ref + noise);
            A(:, :, i) = V*diag(max(.01, diag(D)))*V';
        end
    end

    % Retrieve the size of the problem:
    % There are m matrices of size nxn to average.
    n = size(A, 1);
    m = size(A, 3);
    assert(n == size(A, 2), ...
           ['The slices of A must be square, i.e., the ' ...
	        'first and second dimensions of A must be equal.']);

    % Our search space is the set of positive definite matrices of size n.
    % Notice that this is the only place we specify on which manifold we
    % wish to compute Karcher means. Replacing this factory for another
    % geometry will yield code to compute Karcher means on that other
    % manifold, provided that manifold is equipped with a dist function and
    % a logarithmic map log.
    M = sympositivedefinitefactory(n);

    % Define a problem structure, specifying the manifold M, the cost
    % function and its gradient.
    problem.M = M;
    problem.cost = @cost;
    problem.grad = @grad;

    % Explicitly pick an approximate Hessian for the trust-region method.
    % (This is only to show an example of how it can be done; the solver
    % below, rlbfgs, does not use the approximate Hessian; trustregions
    % would, but it would figure it out automatically with default
    % stepsize if the line below is omitted.)
    problem.approxhess = approxhessianFD(problem, struct('stepsize', 1e-4));

    % The functions below make many redundant computations. This
    % performance hit can be alleviated by using the caching system. We go
    % for a simple implementation here, as a tutorial example.

    % Cost function
    function f = cost(X)
        f = 0;
        for k = 1 : m
            f = f + M.dist(X, A(:, :, k))^2;
        end
        f = f/(2*m);
    end

    % Riemannian gradient of the cost function
    function g = grad(X)
        g = M.zerovec(X);
        for k = 1 : m
            % Update g in a linear combination of the form
            % g = g - [something]/m.
            g = M.lincomb(X, 1, g, -1/m, M.log(X, A(:, :, k)));
        end
    end

    % Execute some checks on the derivatives for early debugging.
    % These things can be commented out of course.
    % The slopes should agree on part of the plot at least. In this case,
    % it is sometimes necessary to inspect the plot visually to make the
    % call, but it is indeed correct.
    % checkgradient(problem);
    % pause;

    % Execute this if you want to force using a proper parallel vector
    % transport. This is not necessary. If you omit this, the default
    % vector transport is the identity map, which is (of course) cheaper
    % and seems to perform well in practice.
    % M.transp = M.paralleltransp;

    % Issue a call to a solver. Default options are selected.
    % Our initial guess is the first data point. Most solvers work well
    % with this problem. Limited-memory BFGS is one good example:
    X = rlbfgs(problem, A(:, :, 1));

end

author: NicolasBoumal

repo: manopt


text:
function [cleanZ] = PatchDCTGG(Z,patchSize,noiseSD,imsize,W,invW,excludeList)
% a simple thresholding denoiser - approx. corresponds to a sprase prior
% over the marginals
if ~exist('excludeList','var')
    excludeList = [];
end

meanZ = mean(Z);
Z = bsxfun(@minus,Z,meanZ);
if (~isempty(excludeList))
    WZ = W*Z(:,excludeList);
else
    WZ = W*Z;
end

t = noiseSD*3;
WZ(abs(WZ)<t)=0;
cleanZ = Z;
if ~isempty(excludeList)
    cleanZ(:,excludeList) = invW*WZ;
else
    cleanZ = invW*WZ;
end
cleanZ = bsxfun(@plus,cleanZ,meanZ);

author: lbasek

repo: image-denoising-benchmark


text:
function prob_test098 ( )

%*****************************************************************************80
%
%% TEST098 tests LOGISTIC_MEAN, LOGISTIC_SAMPLE, LOGISTIC_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST098\n' );
  fprintf ( 1, '  For the Logistic PDF:\n' );
  fprintf ( 1, '  LOGISTIC_MEAN computes the mean;\n' );
  fprintf ( 1, '  LOGISTIC_SAMPLE samples;\n' );
  fprintf ( 1, '  LOGISTIC_VARIANCE computes the variance.\n' );

  a = 2.0;
  b = 3.0;

  check = logistic_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST098 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = logistic_mean ( a, b );
  variance = logistic_variance ( a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF parameter B =             %14f\n', b );
  fprintf ( 1, '  PDF mean =                    %14f\n', mean );
  fprintf ( 1, '  PDF variance =                %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = logistic_sample ( a, b, seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function LMprogressiveinstall(folderlist, HOMEIMAGES, HOMEANNOTATIONS);
%
% Downloads a set of folders from the LabelMe dataset.
%  LMprogressiveinstall(folderlist, HOMEIMAGES, HOMEANNOTATIONS);
%
% This function is useful if you do not want to download the entire
% dataset. You can first browse the dataset, and when you have decided for
% a set of folders that you want to use, you can use this function to
% download only those folders.
%
% Contribute to the dataset by labeling few images:
% http://people.csail.mit.edu/brussell/research/LabelMe/intro.html

Nfolders = length(folderlist);
webpageanno = 'http://people.csail.mit.edu/brussell/research/LabelMe/Annotations'
webpageimg = 'http://people.csail.mit.edu/brussell/research/LabelMe/Images'

% create folders:
disp('create folders');
for i = 1:Nfolders
    [SUCCESS,MESSAGE,MESSAGEID] = mkdir(HOMEIMAGES, folderlist{i});
    [SUCCESS,MESSAGE,MESSAGEID] = mkdir(HOMEANNOTATIONS, folderlist{i});
end

disp('progressively downloading images ..., # of folders:')
for f = 1:Nfolders
    disp(sprintf('%d/%d, %s', f, Nfolders, folderlist{f}))
    wpi = [webpageimg  '/' folderlist{f}];
    images = urldir(wpi, 'IMG');
    Nimages = length(images);
    for i = 1:Nimages
        % if the image already exist, then overwrite
        fp = fopen(fullfile(HOMEIMAGES,folderlist{f},images(i).name));
        if fp < 0
            [F,STATUS] = urlwrite([wpi '/' images(i).name], fullfile(HOMEIMAGES,folderlist{f},images(i).name));
%            disp(sprintf('downloading: %s', images(i).name))
        else
            fclose(fp);
%            disp(sprintf('no need to download: %s', images(i).name))
        end
    end
    end

disp('downloading all annotations..., # of folders:')
for f = 1:Nfolders
    disp(sprintf('%d/%d, %s', f, Nfolders, folderlist{f}))

    wpa = [webpageanno '/' folderlist{f}];
    annotations = urldir(wpa, 'TXT');
    Nanno= length(annotations);
    for i = 1:Nanno
        [F,STATUS] = urlwrite([wpa '/' annotations(i).name], fullfile(HOMEANNOTATIONS,folderlist{f},annotations(i).name));
    end
end

author: CSAILVision

repo: LabelMeToolbox


text:
function msfun_simiam_ui(block)

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

%%
%% The setup method is used to set up the basic attributes of the
%% S-function such as ports, parameters, etc. Do not add any other
%% calls to the main body of the function.
%%
app = [];
setup(block);

 %endfunction

%% Function: setup ===================================================
%% Abstract:
%%   Set up the basic characteristics of the S-function block such as:
%%   - Input ports
%%   - Output ports
%%   - Dialog parameters
%%   - Options
%%
%%   Required         : Yes
%%   C-Mex counterpart: mdlInitializeSizes
%%
function setup(block)

% Register number of ports
block.NumInputPorts  = 0;
block.NumOutputPorts = 0;

% Set up the states
block.NumContStates = 0;
% block.NumDworks = 0;

% Register parameters
block.NumDialogPrms     = 0;

% Register sample times
%  [0 offset]            : Continuous sample time
%  [positive_num offset] : Discrete sample time
%
%  [-1, 0]               : Inherited sample time
%  [-2, 0]               : Variable sample time
block.SampleTimes = [0.033 0];

% Specify the block simStateCompliance. The allowed values are:
%    'UnknownSimState', < The default setting; warn and assume DefaultSimState
%    'DefaultSimState', < Same sim state as a built-in block
%    'HasNoSimState',   < No sim state
%    'CustomSimState',  < Has GetSimState and SetSimState methods
%    'DisallowSimState' < Error out when saving or restoring the model sim state
block.SimStateCompliance = 'DefaultSimState';

%% -----------------------------------------------------------------
%% The MATLAB S-function uses an internal registry for all
%% block methods. You should register all relevant methods
%% (optional and required) as illustrated below. You may choose
%% any suitable name for the methods and implement these methods
%% as local functions within the same file. See comments
%% provided for each function for more information.
%% -----------------------------------------------------------------

block.RegBlockMethod('InitializeConditions', @InitializeConditions);
block.RegBlockMethod('Start', @Start);
block.RegBlockMethod('Outputs', @Outputs);     % Required
block.RegBlockMethod('Update', @Update);
block.RegBlockMethod('Derivatives', @Derivatives);
block.RegBlockMethod('Terminate', @Terminate); % Required

end %end setup

%%
%% InitializeConditions:
%%   Functionality    : Called at the start of simulation and if it is
%%                      present in an enabled subsystem configured to reset
%%                      states, it will be called when the enabled subsystem
%%                      restarts execution to reset the states.
%%   Required         : No
%%   C-MEX counterpart: mdlInitializeConditions
%%
function InitializeConditions(block)

end %end InitializeConditions

%%
%% Start:
%%   Functionality    : Called once at start of model execution. If you
%%                      have states that should be initialized once, this
%%                      is the place to do it.
%%   Required         : No
%%   C-MEX counterpart: mdlStart
%%
function Start(block)

root_path = fullfile(fileparts(mfilename('fullpath')), '../');
addpath(genpath(root_path));
root_path

app = simiam.ui.AppWindow(root_path, true);
app.load_ui();

% store info in custom data;
customData = containers.Map('UniformValues', false);
simulator = app.simulator_;
customData('simulatorHandle') = simulator;
set(block.BlockHandle, 'UserData', customData, 'UserDataPersistent', 'off');

end %endfunction

%%
%% Outputs:
%%   Functionality    : Called to generate block outputs in
%%                      simulation step
%%   Required         : Yes
%%   C-MEX counterpart: mdlOutputs
%%
function Outputs(block)

end %end Outputs

%%
%% Update:
%%   Functionality    : Called to update discrete states
%%                      during simulation step
%%   Required         : No
%%   C-MEX counterpart: mdlUpdate
%%
function Update(block)

customData = get(block.BlockHandle, 'UserData');
simulator = customData('simulatorHandle');

if(~app.is_state_crashed_)
    simulator.step([],[]);
else
    error('Collision detected!');
end

end %end Update

%%
%% Derivatives:
%%   Functionality    : Called to update derivatives of
%%                      continuous states during simulation step
%%   Required         : No
%%   C-MEX counterpart: mdlDerivatives
%%
function Derivatives(block)

end %end Derivatives

%%
%% Terminate:
%%   Functionality    : Called at the end of simulation for cleanup
%%   Required         : Yes
%%   C-MEX counterpart: mdlTerminate
%%
function Terminate(block)

end %end Terminate

end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
function r8r8vec_index_insert_unique_test ( )

%*****************************************************************************80
%
%% R8R8VEC_INDEX_INSERT_UNIQUE_TEST tests R8R8VEC_INDEX_INSERT_UNIQUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  n_max = 30;

  x_max = 4.0;
  x_min = 1.0;
  y_max = 3.0;
  y_min = 1.0;

  n = 0;
  x = [];
  y = [];
  indx = [];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8R8VEC_INDEX_INSERT_UNIQUE_TEST\n' );
  fprintf ( 1, '  R8R8VEC_INDEX_INSERT_UNIQUE inserts unique values into an\n' );
  fprintf ( 1, '    index sorted array.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Generate %d random values:\n', n_max );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Index    XVAL    YVAL\n' );
  fprintf ( 1, '\n' );

  seed = 123456789;

  for i = 1 : n_max

    [ xval, seed ] = r8_uniform_ab ( x_min, x_max, seed );
    xval = round ( xval );
    [ yval, seed ] = r8_uniform_ab ( y_min, y_max, seed );
    yval = round ( yval );

    [ n, x, y, indx, ival, ierror ] = r8r8vec_index_insert_unique ( ...
      n, x, y, indx, xval, yval );

    fprintf ( 1, '    %6d  %6f  %6f\n', ival, xval, yval );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Vector of unique X Y values:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I  X(I)   Y(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '    %6d  %6f  %6f\n', i, x(i), y(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  X, Y sorted by index\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I  INDX(I)  X(INDX(I))  Y(INDX(I))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '    %6d  %6f  %6f\n', i, indx(i), x(indx(i)), y(indx(i)) );
  end

  return
end


text:
function D=bwdistsc(bw,aspect)
% D=BWDISTSC(BW,ASPECT)
% BWDISTSC computes Euclidean distance transform of a binary 3D image BW.
% Distance transform assigns to each pixel in BW a number that is the
% distance from that pixel to the nearest nonzero pixel in BW. BWDISTSC
% can accept a regular 2D image, a 3D array, and a cell array of 2D slices.
% ASPECT is a 3-component vector defining the voxel-aspect-ratio for BW.
% If ASPECT is not given, [1 1 1] isotropic aspect ratio is assumed.
%
% BWDISTSC uses fast optimized scan algorithm and cell-arrays to
% represent internal data, and is less demanding to physical memory as
% well as in many cases up to 10 times faster than MATLAB's native bwdist.
%
% Example:
% bw=zeros(100,100,100);
% bw(40:60,40:60,40:60)=1;
% tic;D=bwdist(bw);toc
% tic;D=bwdistsc(bw);toc
%
% BWDISTSC tries to use MATLAB bwdist from image processing toolbox for 2D
% scans if possible, which is faster, otherwise BWDISTSC will use its own
% algorithm to also perform 2D scans. Own algorithm is also used if x- and
% y-anisotropy scales are not equal; therefore, if your data has only one
% axis that is anisotropic, it is always advantageous to feed it to
% BWDISTSC so that the anisotropic axis is z.
%
%(c) Yuriy Mishchenko HHMI JFRC Chklovskii Lab JUL 2007
% Updated Yuriy Mishchenko (Toros University) SEP 2013

% This implementation uses optimized forward-backward scan version of the
% algorithm of the original bwdistsc (2007), which substantially improves
% its speed and simplifies the code. The improvement is described in the
% part on the selection initial point in the SIVP paper below. The original
% implementation is still used in bwdistsc1, since forward-backward scan
% does not allow limiting computation to a fixed distance value MAXVAL.

% This code is free for use or modifications, just please give credit
% where appropriate. And if you modify code or fix bugs, please drop
% me a message at gmyuriy@hotmail.com.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Scan algorithms below use following Lema:                     %
% LEMA: let F(X,z) be lower envelope of a family of parabola:   %
% F(X,z)=min_{i} [G_i(X)+(z-k_i)^2];                            %
% and let H_k(X,z)=A(X)+(z-k)^2 be a parabola.                  %
% Then for H_k(X,z)==F(X,z) at each X there exist at most       %
% two solutions k1<k2 such that H_k12(X,z)=F(X,z), and          %
% H_k(X,z)<F(X,z) is restricted to at most k1<k2.               %
% Here X is any-dimensional coordinate.                         %
%                                                               %
% Thus, simply scan away from any z such that H_k(X,z)<F(X,z)   %
% in either direction as long as H_k(X,z)<F(X,z) and update     %
% F(X,z). Note that need to properly choose starting point;     %
% starting point is any z such that H_k(X,z)<F(X,z); z==k is    %
% usually, but not always the starting point!!!                 %
% usually, but not always the starting point!                   %
%                                                               %
% Citation:                                                     %
% Mishchenko Y. (2013) A function for fastcomputation of large  %
% discrete Euclidean distance transforms in three or more       %
% dimensions in Matlab. Signal, Image and Video Processing      %
% DOI: 10.1007/s11760-012-0419-9.                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% parse inputs
if(nargin<2 || isempty(aspect)) aspect=[1 1 1]; end

% determine geometry of the data
if(iscell(bw)) shape=[size(bw{1}),length(bw)]; else shape=size(bw); end

% correct this for 2D data
if(length(shape)==2) shape=[shape,1]; end
if(length(aspect)==2) aspect=[aspect,1]; end

% allocate internal memory
D=cell(1,shape(3)); for k=1:shape(3) D{k}=zeros(shape(1:2)); end

%%%%%%%%%%%%% scan along XY %%%%%%%%%%%%%%%%
for k=1:shape(3)
    if(iscell(bw)) bwXY=bw{k}; else bwXY=bw(:,:,k); end

    % initialize arrays
    DXY=zeros(shape(1:2));
    D1=zeros(shape(1:2));

    % if can, use 2D bwdist from image processing toolbox
    if(exist('bwdist') && aspect(1)==aspect(2))
        D1=aspect(1)^2*bwdist(bwXY).^2;
    else    % if not, use full XY-scan
        %%%%%%%%%%%%%%% X-SCAN %%%%%%%%%%%%%%%
        % reference for nearest "on"-pixel in bw in x direction down

        %  scan bottow-up (for all y), copy x-reference from previous row
        %  unless there is "on"-pixel in that point in current row, then
        %  that's the nearest pixel now
        xlower=repmat(Inf,shape(1:2));

        xlower(1,find(bwXY(1,:)))=1;    % fill in first row
        for i=2:shape(1)
            xlower(i,:)=xlower(i-1,:);  % copy previous row
            xlower(i,find(bwXY(i,:)))=i;% unless there is pixel
        end

        % reference for nearest "on"-pixel in bw in x direction up
        xupper=repmat(Inf,shape(1:2));

        xupper(end,find(bwXY(end,:)))=shape(1);
        for i=shape(1)-1:-1:1
            xupper(i,:)=xupper(i+1,:);
            xupper(i,find(bwXY(i,:)))=i;
        end

        % build (X,Y) for points for which distance needs to be calculated
        idx=find(~bwXY); [x,y]=ind2sub(shape(1:2),idx);

        % update distances as shortest to "on" pixels up/down in the above
        DXY(idx)=aspect(1)^2*min((x-xlower(idx)).^2,(x-xupper(idx)).^2);

        %%%%%%%%%%%%%%% Y-SCAN %%%%%%%%%%%%%%%
        % this will be the envelop of parabolas at different y
        D1=repmat(Inf,shape(1:2));

        p=shape(2);
        for i=1:shape(2)
            % some auxiliary datasets
            d0=DXY(:,i);

            % selecting starting point for x:
            % * if parabolas are incremented in increasing order of y,
            %   then all below-envelop intervals are necessarily right-
            %   open, which means starting point can always be chosen
            %   at the right end of y-axis
            % * if starting point exists it should be below existing
            %   current envelop at the right end of y-axis
            dtmp=d0+aspect(2)^2*(p-i)^2;
            L=D1(:,p)>dtmp;
            idx=find(L);
            D1(idx,p)=dtmp(L);

            % these will keep track along which X should
            % keep updating distances
            map_lower=L;
            idx_lower=idx;

            % scan from starting points down in increments of 1
            for ii=p-1:-1:1
                % new values for D
                dtmp=d0(idx_lower)+aspect(2)^2*(ii-i)^2;

                % these pixels are to be updated
                L=D1(idx_lower,ii)>dtmp;
                D1(idx_lower(L),ii)=dtmp(L);

                % other pixels are removed from scan
                map_lower(idx_lower)=L;
                idx_lower=idx_lower(L);

                if(isempty(idx_lower)) break; end
            end
        end
    end
    D{k}=D1;
end

%%%%%%%%%%%%% scan along Z %%%%%%%%%%%%%%%%
D1=cell(size(D));
for k=1:shape(3)
  D1{k}=repmat(Inf,shape(1:2));
end

% start building the envelope
p=shape(3);
for k=1:shape(3)
    % if there are no objects in this slice, nothing to do
    if(isinf(D{k}(1,1)))
      continue;
    end

    % selecting starting point for (x,y):
    % * if parabolas are incremented in increasing order of k, then all
    %   intersections are necessarily at the right end of the envelop,
    %   and so the starting point can be always chosen as the right end
    %   of the axis

    % check which points are valid starting points, & update the envelop
    dtmp=D{k}+aspect(3)^2*(p-k)^2;
    L=D1{p}>dtmp;
    D1{p}(L)=dtmp(L);

    % map_lower keeps track of which pixels can be yet updated with the
    % new distance, i.e. all such XY that had been under the envelop for
    % all Deltak up to now, for Deltak<0
    map_lower=L;

    % these are maintained to keep fast track of whether map is empty
    idx_lower=find(map_lower);

    % scan away from the starting points in increments of -1
    for kk=p-1:-1:1
        % new values for D
        dtmp=D{k}(idx_lower)+aspect(3)^2*(kk-k)^2;

        % these pixels are to be updated
        L=D1{kk}(idx_lower)>dtmp;
        map_lower(idx_lower)=L;
        D1{kk}(idx_lower(L))=dtmp(L);

        % other pixels are removed from scan
        idx_lower=idx_lower(L);

        if(isempty(idx_lower)) break; end
    end
end

% prepare the answer
if(iscell(bw))
    D=cell(size(bw));
    for k=1:shape(3) D{k}=sqrt(D1{k}); end
else
    D=zeros(shape);
    for k=1:shape(3) D(:,:,k)=sqrt(D1{k}); end
end

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
classdef dme_load_opf < mp.dme_load & mp.dme_shared_opf

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

%     properties
%     end     %% properties

%     methods
%     end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function value = gen_hermite_integral ( expon, alpha )

%*****************************************************************************80
%
%% GEN_HERMITE_INTEGRAL evaluates a monomial generalized Hermite integral.
%
%  Discussion:
%
%    The integral:
%
%      integral ( -oo < x < +oo ) x^n |x|^alpha exp(-x^2) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, int EXPON, the exponent of the monomial.
%
%    Input, real ALPHA, the exponent of |X| in the integral.
%    -1.0 < ALPHA.
%
%    Output, real VALUE, the value of the integral.
%
  if ( mod ( expon, 2 ) == 1 )

    value = 0.0;

  else

    a = alpha + expon;

    if ( a <= -1.0 )

      value = - r8_huge ( );

    else

      value = gamma ( ( a + 1.0 ) / 2.0 );

    end

  end

  return
end


text:
function superima = QLsuperimage(C, S, iscale)
%------------------------------------------------------------------------------
% Creates an image that goes with a multiresolutiondecomposition that has been
% made in the way of LISQ. The superimage, that is an image that comprises
% all levels, is constructed as a rectangular gridfunction.
%
% C        is a one-dimensional array that contains the coefficients of the
%          Quincunx Lifting Scheme (LISQ) decompositions.
% S        is the bookkeeping vector
% iscale   optional inputparameter
%          1 show superimage as is (default)
%          2 show superimage after removal of outliers
%          3 show superimage after enhancement of contrast using
%            histogram equalisation (if available)
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: November 24, 2003.
%  2003 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
 if nargin == 3
   iopt = iscale;
 elseif nargin == 2
   iopt = 1;
 else
   error(' QLsuperimage - number of arguments should be either 2 or 3 ');
 end
%
 [nS, mS] = size(S);
 if mS ~= 6
   error(' QLsuperimage - unexpected dimensions of bookkeeping vector ');
 end
levels = S(nS,1);
if mod(levels, 2) == 1
  error('  QLsuperimage - only an even number of levels is accepted ');
elseif levels < 1
  superima = [];
else
%
% Firstly, we determine the dimensions of the superimage to be and create the
% necessary space (see LISQ/storeR.m and LISQ/QLiftRec2Nevill.m)
%
  [Detail10, Detail01] = retrieveQ1001(1, 'd', C, S);
  sizeR = size(Detail10) + size(Detail01);
  n  = sizeR(1); m = sizeR(2);
  nH = round((n+m+1.999)/2)+round((m+0.999)/2);
  nm = round((n+m+1.999)/2)+1;
  nV = 1;
  levels = S(nS,1);
  for k = 1:2:levels
    nV = nV + nm;
    nm = round((nm+0.999)/2);
  end
  Approx = retrieveR(levels, 'a', C, S);
  [nl, ml] = size(Approx); nV = nV + nl;
  superima = ones(nV, nH);
%
% Secondly, we fill the superimage with all subsequent detail coefficients
% Put the details at the rotated quincunx grid (level 1)
  background=max(max(max(Detail10)), max(max(Detail01)));
  RotaDet=rota1001fill(Detail10, Detail01, background);
  clear Detail10 Detail01;
%-----intermezzo---begin--
  hgram=ones(1,64);
  hgram(64)=round(31.5*(n/m+m/n))+1;
% hgram is a histogram that is desired to be matched by the image after
% its values have been transformed by histeq, see the documentation of
% histeq() in the Matlab Image Processing Toolbox.
% Due to the rotation of the quincunx grid to make it visible, there
% is a substantial area for padding with one colour (white or nearly
% white). Here follows an overview of the quantities of pixels involved
% at the rotated quincunx grid (asymptotically):
%  Let n and m be the dimensions of the original image.
%
%                               2
%                      ( n + m )
%  Grand total:         -------
%                          4
%
%                         n m
%  Used pixels:           ---
%                          2
%
%                        2    2
%                       n  + m
%  Pixels for padding:  --------
%                          4
%
%                            2    2
%                  n m      n  + m                 n   m
%  Used:padding =  ---   :  -------  = 63 : 31.5*( - + - )
%                   2          4                   m   n
%
%-----intermezzo---end----
  [nl, ml] = size(RotaDet);
  OrigV=1; OrigH=1;
  superima(OrigV:(OrigV+nl-1), OrigH:(OrigH+ml-1)) = imequi(RotaDet,iopt,hgram);
  Detail11 = retrieveR(2, 'd', C, S);
  OrigV=1; OrigH=ml+1;
  [nr, mr] = size(Detail11);
  superima(OrigV:(OrigV+nr-1), OrigH:(OrigH+mr-1)) = imequi(Detail11,iopt);
% Now the other levels, if any
  OrigV=OrigV+nl; OrigH=1;
  for k = 3:2:levels
%   Put the details at the rotated quincunx grid
    [Detail10, Detail01] = retrieveQ1001(k, 'd', C, S);
    background=max([max(max(Detail10))  max(max(Detail01))]);
    RotaDet = rota1001fill(Detail10, Detail01, background);
    clear Detail10 Detail01;
    [nl, ml] = size(RotaDet);
    superima(OrigV:(OrigV+nl-1), OrigH:(OrigH+ml-1)) = imequi(RotaDet,iopt,hgram);
    Detail11 = retrieveR(k+1, 'd', C, S);
    OrigH=ml+1;
    [nr, mr] = size(Detail11);
    superima(OrigV:(OrigV+nr-1), OrigH:(OrigH+mr-1)) = imequi(Detail11,iopt);
    clear Detail11;
    OrigV=OrigV+nl; OrigH=1;
  end
% Thirdly, we fill the superimage with approximation coefficients.
  [nl, ml] = size(Approx);
  superima(OrigV:(OrigV+nl-1), OrigH:(OrigH+ml-1)) = imequi(Approx,iopt);
end
%------------------------------------------------------------------------------


text:
function func_LuminanceMasking
% Show luminance masking effect of HVS
% This function is part of the toolbox "Basic introduction to HVS"
%
% Command line
% ----------------------
% func_LuminanceMasking
% input:  None
%
% output: CSF figure
%
% More information can be found in
% S. E. Palmer, "Vision Science: From Photons to Phenomenology," MIT Press,
% Cambridge, MA, 1999.
%
%
% Jing Tian Apr.24 2004
% Contact me : scuteejtian@hotmail.com
% Homepage : http://ikanchi.yeah.net
% This program is written in Apr.2003 during my postgraduate in
% NTU, Singapore.
% ----------------------

freq = 1;
C = 0.05;
L = 100;

x = linspace(-1.5 * pi, 0.5 * pi, 100);
y = linspace(150, 50, 100);
[xx,yy] = meshgrid(x, y);

i = 1;
for L = 100 : 20 : 200;
	z = L .* C .* sin(2 .* pi .* freq .* xx) + L;
	imagesc(z);
	colormap gray;
	shading interp;
	ch = ['luminancemasking', num2str(i), '.bmp'];
	imwrite(z, gray(256), ch, 'bmp');
	i = i + 1;
end


text:
% CLASSIF_RECOG Calculates the average recognition rate
%
% Usage
%    [rr_mean,recog_rate] =  CLASSIF_RECOG(labels,test_set,truth)
%
% Input
%    labels (int): The labels attributed to the testing instances.
%    test_set (int): The object indices of the testing instances.
%    truth: the actual labels of the testing instances
%
% Output
%    rr_mean (real): The mean recognition rate
%    recog_rate(real): array containing the individual recognition rates of
%       each class.
%
% Description
%    This function computes the average recognition rate of each class.
%    It is a widely used measure of the performance of the classifier when the
%    data set is highly unbalanced.
%
% See also
%    SVM_TEST, CLASSIF_ERR, CLASSIF_RECOG

function [rg_mean,recog_rate]=classif_recog(labels,test_set,truth)
    if isstruct(truth)
        src=truth;
        truth = [src.objects.class];
    end

    % Normally, the test_set contains samples of all classes
    recog_rate=zeros(1,max(truth));
    gdTruth=truth(:,test_set);

    for k=1:max(truth)
        mask=k==gdTruth;
        good_elts=find(labels==k & mask);

        mask1=numel(find(mask==1));
        recog_rate(k)=numel(good_elts)/mask1;
    end

    rg_mean=mean(recog_rate);
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
% NDIMS - number of dimension of memory mapped underlying array
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Nov. 2008

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function res = ndims(obj)

    if length(obj.dimensions) <= 2 || all(obj.dimensions(3:end) == 1), res = 2;
    else                                                               res = length(obj.dimensions);
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [X, varX, model, grChek] = svargplvmOptimisePoint(model, vardistx, y, display, iters, varargin)

% SVARGPLVMOPTIMISEPOINT Optimise the postion of one or more latent points
% given observations from 1 or more modalities
% FORMAT
% DESC optimises the location of a group of points in latent space
% given an initialisation and the corresponding observed data point.
% ARG model : the MRD model for which the point will be optimised.
% ARG vardistx : the initialisation of the points in the latent space.
% ARG y : the observed data points for which the latent points are to
% be optimised.
% ARG display : whether or not to display the iterations of the
% optimisation (default: true)
% ARG iters : maximum number of iterations for the optimisation
% (default 2000).
% RETURN x : the optimised means in the latent space.
% RETURN varx : the optimised variances in the latent space.
% RETURN model: the model which is augmented to also include the new test
% points and the quantities that change because of these, as there is
% coupling in the dynamics case.

%
% COPYRIGHT :  Andreas Damianou, 2013
%
% SEEALSO : vargplvmOptimisePoint

% VARGPLVM

grChek = [];

if nargin < 5 || isempty(iters), iters = 2000; end
if nargin < 4 || isempty(display), display = true; end

options = optOptions;
if display
    options(1) = 1;
    % options(9) = 1; % gradchek
end
options(14) = iters;

if isfield(model, 'optimiser')
    optim = str2func(model.optimiser);
else
    optim = str2func('scg');
end

if isfield(model, 'dynamics') && ~isempty(model.dynamics)
    % TODO
    error('Not implemented for dynamical case')
else
    x = vardistExtractParam(vardistx);
end

if ~(isfield(model, 'dynamics') && ~isempty(model.dynamics))
    model.DgtN_test = true;
    for i = model.testModalities
        % Perform precomputations which will result in faster execution. This
        % happens in two stages:
        % a) Precomputing constants only once here, instead in every call of
        % the objective and grad
        % b) Replace model.m with a reduced rank representation, since it only
        % appears in the expression model.m * model.m'. This is the same tricks
        % used in vargplvmCreate for DgtN flag.
        % TODO: Do the same for dynamics case
        if isfield(model.comp{i}, 'DgtN') && model.comp{i}.DgtN
            mOrig{i} = model.comp{i}.m;
        end
        model.comp{i}.testPrecomp.indexMissing = find(isnan(y{i}(1,:)));
        indexPresent = setdiff(1:model.comp{i}.d, model.comp{i}.testPrecomp.indexMissing );
        y{i} = y{i}(:,indexPresent);
        P = model.comp{i}.P1 * (model.comp{i}.Psi1' * model.comp{i}.m(:,model.comp{i}.testPrecomp.indexMissing));
        model.comp{i}.testPrecomp.TrPP = sum(sum(P .* P));
        model.comp{i}.testPrecomp.TrYY = sum(sum(model.comp{i}.m(:,model.comp{i}.testPrecomp.indexMissing) .* model.comp{i}.m(:,model.comp{i}.testPrecomp.indexMissing)));
        y{i} = y{i} - repmat(model.comp{i}.bias(indexPresent),size(y{i},1),1);
        y{i} = y{i}./repmat(model.comp{i}.scale(indexPresent),size(y{i},1),1);
        mPres = model.comp{i}.m(:, indexPresent);
        mPres = [mPres; y{i}];
        YYT = mPres * mPres'; % NxN
        [U S V]=svd(YYT);
        if isfield(model.comp{i}, 'DgtN') && model.comp{i}.DgtN
            model.comp{i}.testPrecomp.mReduced=U*sqrt(abs(S));
        end
        model.comp{i}.testPrecomp.TrYY2 = sum(diag(YYT)); % scalar

        % For grads
        mPresGrad = [y{i}; model.comp{i}.m(:, indexPresent)];
        YYT = mPresGrad * mPresGrad'; % NxN
        [U S V]=svd(YYT);
        if isfield(model.comp{i}, 'DgtN') && model.comp{i}.DgtN
            model.comp{i}.testPrecomp.mReducedGrad=U*sqrt(abs(S));
        end
        model.comp{i}.testPrecomp.mY = mPresGrad*y{i}';
        if isfield(model.comp{i}, 'DgtN') && model.comp{i}.DgtN
            model.comp{i}.m = []; % Less overhead in passing arguments (pass by value)
        end
    end
else
    model.comp.DgtN_test = false;
end

if length(varargin) == 2
    if strcmp(varargin{1}, 'gradcheck')
        assert(islogical(varargin{2}));
        %options(9) = varargin{2};
        doGradchek = varargin{2};
        if doGradchek
            [gradient, delta] = feval('gradchek', vardistExtractParam(vardistx), @svargplvmPointObjective, @svargplvmPointGradient, model, y);
            deltaf = gradient - delta;
            d=norm(deltaf - gradient)/norm(gradient + deltaf); %%
            d1=norm(deltaf - gradient,1)/norm(gradient + deltaf,1); %%
            grRatio = sum(abs(gradient ./ deltaf)) / length(deltaf);
            fprintf(1,' Norm1 difference: %d\n Norm2 difference: %d\n Ratio: %d\n',d1,d, grRatio);
            grChek = {delta, gradient, deltaf, d, d1};
        else
            grChek = [];
        end
    end
end

if iters > 0
    if strcmp(func2str(optim), 'optimiMinimize')
        % Carl Rasmussen's minimize function
        x = optim('vargplvmPointObjectiveGradient', x, options, model, y);
    else
        % NETLAB style optimization.
        x = optim('svargplvmPointObjective', x,  options, ...
            'svargplvmPointGradient', model, y);
    end
end

if isfield(model, 'dynamics') && ~isempty(model.dynamics)
    if isfield(model.dynamics, 'reoptimise') && model.dynamics.reoptimise   %%% RE-OPT-CODE-NEW
        % Expand and update the model.                                      %%% RE-OPT-CODE-NEW
        % In this case, x is [mu_bar lambda theta_t X_u]                    %%% RE-OPT-CODE-NEW
        [vardistx, model] = vargplvmPartExpand(model, x, 1);                %%% RE-OPT-CODE-NEW
    else %%% RE-OPT-CODE-NEW
        % now separate the variational disribution into the training part and the
        % testing part and update the original training model (only with the new training
        % variational distribution) and the test variational distribution
        % this is doing the expand
        x = reshape(x, vardist.numData, model.dynamics.q*2);
        xtrain = x(1:model.N,:);
        xtest = x(model.N+1:end,:);
        model.dynamics.vardist = vardistExpandParam(model.dynamics.vardist, xtrain);
        vardistx = vardistExpandParam(model.vardistx, xtest);
    end %%% RE-OPT-CODE-NEW
else
    vardistx = vardistExpandParam(vardistx,x);
end

X = vardistx.means;
varX = vardistx.covars;

if ~(isfield(model, 'dynamics') && ~isempty(model.dynamics))
    model = rmfield(model, 'DgtN_test');
    for i=model.testModalities
        model.comp{i} = rmfield(model.comp{i}, 'testPrecomp');
        if isfield(model.comp{i}, 'DgtN') && model.comp{i}.DgtN
            model.comp{i}.m = mOrig{i};
        end
    end
end

author: SheffieldML

repo: vargplvm


text:
function FigFont(opt)
% =======================================================================
% Set font in a figure to desired font size and font style
% =======================================================================
% FigFont(opt)
% -----------------------------------------------------------------------
% OPTIONAL INPUT
%   - opt: output of FigFontOption(fsize)
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

%% CHECK INPUT
% =======================================================================
if ~exist('opt','var')
    opt = FigFontOption;
end

%% SET FONT STYLE
% =======================================================================
% AXES
aux = findobj(gcf,'Type','axes');
set(aux,'Fontsize',opt.axes_size,'FontWeight',opt.axes_weight,'FontName',opt.axes_name);

% TITLES
if max(size(aux))==1
    aux = get(aux, 'title');
    set(aux,'Fontsize',opt.title_size,'FontWeight',opt.title_weight,'FontName',opt.title_name);
else
    aux1 = get(aux(1), 'title');
    set(aux1, 'Fontsize', opt.title_size,'FontWeight',opt.title_weight,'FontName',opt.title_name);
    aux2 = get(aux(2), 'title');
    set(aux2, 'Fontsize',opt.title_size,'FontWeight',opt.title_weight,'FontName',opt.title_name);
end

% SUPTITLE
aux = findobj(gcf,'tag','suptitleText');
set(aux,'Fontsize',opt.suptitle_size,'FontWeight',opt.suptitle_weight,'FontName',opt.suptitle_name);

% LEGEND
aux = findobj(gcf,'Type','axes','Tag','legend');
set(aux,'Fontsize',opt.legend_size, 'FontWeight',opt.legend_weight,'FontName',opt.legend_name);

% Y LABEL
aux = get(gca,'YLabel');
if isempty(aux)~=1
    set(aux,'Fontsize',opt.ylabel_size, 'FontWeight',opt.ylabel_weight,'FontName',opt.ylabel_name);
end

% X LABEL
aux = get(gca,'XLabel');
if isempty(aux)~=1
    set(aux,'Fontsize',opt.xlabel_size, 'FontWeight',opt.xlabel_weight,'FontName',opt.xlabel_name);
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox