text:
% testGraphCutsScript3

if 0
    imdir = '../images/all_images';
    outdir = '../data';
    load([outdir '/rand_indices.mat']);
    load([outdir '/allimsegs2.mat']);
    load([outdir '/mcmcSuperpixelClassifier.mat']);
    load([outdir '/mcmcEdgeData.mat']);
    load([outdir '/mcmcEdgeClassifier.mat']);
end

alpha = [1:1:3];

for k = 1:numel(alpha);

    disp(num2str(alpha(k)))

    for cvf = 1:numel(cv_images)

        disp([num2str(cvf) ': ' imsegs(cv_images(cvf)).imname])

        f = cv_images(cvf);
        c = ceil(cvf/50);

        im = im2double(imread([imdir '/' imsegs(f).imname]));
        wseg = segmentWatershed(im, 1);

        cimages = pg2confidenceImages(imsegs(f), pg(cvf));
        [labv, labh] = testImageGraphCuts3(im, wseg, cimages{1}, alpha(k));

        nsp = numel(labv);
        pg2{cvf} = zeros(nsp, 7);
        lab = (labv==1) + (labv==2).*(1+labh) + 7*(labv==3);
        % set most likely label
        pg2{cvf}((lab-1)*nsp + [1:nsp]') = 1;
        % set so that h accuracy can be determined
        pg2{cvf}((labh)*nsp + [1:nsp]') = 0.5;

        [pv1, ph1] = splitpg(pg{cvf});
        lim1 = APPgetLabeledImage2(im, imsegs(f), pv1, ph1);
        [pv2, ph2] = splitpg(pg2{cvf});
        lim2 = APPgetLabeledImage2(im, wseg, pv2, ph2);
        figure(1), hold off, imagesc(lim1), axis image;
        figure(2), hold off, imagesc(lim2), axis image;
        drawnow; pause(1)
    end

    [vaccgc(k), haccgc(k)] = mcmcProcessResult(imsegs(cv_images), pg2);
    disp(num2str([vaccgc ; haccgc]))
end

author: zouchuhang

repo: LayoutNet


text:
% Section 8.1.1: Separating a point from a polyhedron
% Boyd & Vandenberghe "Convex Optimization"
% Joelle Skaf - 10/09/05
%
% The goal is to produce a hyperplane separating x0 and the polyhedron
% defined as {x | Ax <= b}
%           minimize    mu'*x0 - b'*lambda
%                       A'*lambda = mu
%                       norm(mu)* <= 1
%                       lambda >= 0

% Input data
randn('seed',0);
n  = 10;
m  = 2*n;
x0 = randn(n,1);
A  = randn(m,n);
b  = rand(m,1);

% CVX solution
fprintf(1,'Finding a separating hyperplane between the 2 polyhedra...');

cvx_begin quiet
    variables muu(n) lambda(m)
    maximize ( muu'*x0 - b'*lambda )
    A'*lambda == muu; %#ok
    norm(muu) <= 1; %#ok
    lambda >= 0; %#ok
cvx_end

fprintf(1,'Done! \n');

% Verification
disp('------------------------------------------------------------------');
disp('Note that 0 is in {x | Ax <= b} by construction...' );
disp('Verifying that x0 is separated from {x | Ax <= b} i.e. mu^T*x0 > 0');
disp([' mu^T*x0 = ' num2str(muu'*x0) ]);

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [gKern1, gKern2, gVarmeans, gVarcovars, gInd] = linard2biasVardistPsi2Gradient(linardKern, biasKern, vardist, Z, covGrad, learnInducing)

% LINARD2BIASVARDISTPSI2GRADIENT description.

% VARGPLVM

if nargin < 6
    learnInducing = 1;
end

% variational means
N = size(vardist.means,1);
%  inducing variables
[M Q] = size(Z);

[Psi2, Pnobias, Psi1] = linard2biasVardistPsi2Compute(linardKern, biasKern, vardist, Z);

% inverse variances
A = linardKern.inputScales;

% gradient for the bias parameter
gKern2 = sum(sum(Pnobias.*covGrad));
Bnm = biasKern.variance*ones(size(Psi1));
BPsi1Covg = (Bnm*covGrad);

gInd = zeros(M, Q);

for q=1:vardist.latentDimension
   %
   gKern(q) = sum(sum((vardist.means(:,q)*(Z(:,q)')).*BPsi1Covg));

   gVarmeans(:,q) = A(q)*sum((ones(vardist.numData,1)*Z(:,q)').*BPsi1Covg,2);
   if learnInducing
        gInd(:,q) = A(q)*sum((vardist.means(:,q)*ones(1,size(Z,1))).*BPsi1Covg,1)';
   end
   %
end
%

gKern1 = 2*gKern(:)';
% gVarmeans is N x Q matrix (N:number of data, Q:latent dimension)
% this will unfold this matrix column-wise
gVarmeans = 2*gVarmeans(:)';

% gInd is M x Q matrix (M:number of inducing variables, Q:latent dimension)
% this will unfold this matrix column-wise
gInd = 2*gInd(:)';

gVarcovars = zeros(1,prod(size(vardist.covars)));

author: SheffieldML

repo: vargplvm


text:
%% Example of clustering with AIC and BIC computed for K-means with RSS
clear all; close all;

%% Generate Perfect circle

nb_samples  = 100;
phi         = linspace(0,2 * pi,nb_samples);
r           = 10;
X           = [r .* cos(phi(:)), r .* sin(phi(:))];

%% Plot circle points

plot_options                = [];
plot_options.title          = 'Circle';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1  = ml_plot_data(X,plot_options);
axis square;

%% Generate Random Clusters from GMM

clear all; close all;
num_samples         = 400;
num_classes         = 4;
dim                 = 2;
[X,labels,gmm]      = ml_clusters_data(num_samples,dim,num_classes);
options.title       = 'Random Gaussian';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1 = ml_plot_data(X,options);
axis equal;

%% K-means

cluster_options             = [];
cluster_options.method_name = 'kmeans';
cluster_options.K           = 2;

result1                     = ml_clustering(X,cluster_options,'Start','plus','Distance','sqeuclidean','MaxIter',500);

result1.lambda              = 2;
[rss,aic,bic]               = ml_clustering_eval(X,result1);

% Plot decision boundary
if exist('hd','var') && isvalid(hd), delete(hd);end
hd = ml_plot_class_boundary(X,result1);
%% Search for optimal \lambda value in AIC(RSS)
clc;
cluster_options             = [];
cluster_options.method_name = 'kmeans';

% Compute Maximum RSS values (when K=1,\alpha)
[N,D]       = size(X);
alpha       = N/4;
repeats     = 1;
K           = 1:alpha;
[mus, stds] = ml_clustering_optimise(X,K,repeats,cluster_options,'Start','plus','Distance','sqeuclidean','MaxIter',500);

% Compute lower/upper bounds of \lambda values according to limits derived
% in slide 9,12,13 of TP2-Metrics-Recap
rss_k1         = mus(1,1);
rss_k2         = mus(1,alpha);
upper_lambda   = (rss_k1 - rss_k2) / (alpha*D);
lower_lambda   = (rss_k1) / (D*N);
lambda_range   = [floor(lower_lambda):1:ceil(upper_lambda)+1];

% Plot AIC(RSS) with a constrained range of lambda's
clc;
if exist('h_aic','var')     && isvalid(h_aic),     delete(h_aic);    end
h_aic = figure;hold on;
leg_names = [];
aic_mins = [];
for l=1:length(lambda_range)
    cluster_options.lambda    = lambda_range(l); %Magic number
    repeats                   = 1;
    Ks                        = 1:N;
    [mus, stds]               = ml_clustering_optimise(X,Ks,repeats,cluster_options,'Start','plus','Distance','sqeuclidean','MaxIter',500);
    plot(mus(2,:),'-');
    hold on;
    [min_aic, min_k] = min(mus(2,:));
    leg_names = [leg_names; {strcat('\lambda= ',num2str(lambda_range(l)))}];
    aic_mins = [aic_mins; [min_aic, min_k]];
end
legend(leg_names);
grid on;
box on;
xlabel('K');
title('AIC(RSS) for K-means')

%% Plot BIC(RSS)

if exist('h_bic','var')     && isvalid(h_bic),     delete(h_bic);    end
h_bic = figure;hold on;
plot(mus(3,:),'-');
grid on;
box on;
xlabel('K');
title('BIC(RSS) for K-means')

author: epfl-lasa


text:
classdef  VideoMaker_TopOpt_levelSet3D < VideoMaker
  properties

  end

  methods (Access = public)

        function obj = VideoMaker_TopOpt_levelSet3D(cParams)
            obj.init(cParams);
        end

        function Make_video_design_variable(obj,output_video_name)
            post = Postprocess_TopOpt_levelSet();
            field2print = post.levelSet_name;
            componentfield = post.levelSet_name_component;
            obj.Make_video_characteristic_function(field2print,componentfield,output_video_name)
        end

  end

  methods (Access = private)
         function Make_video_characteristic_function(obj,field2print,componentfield,output_video_name_in)
            file_tcl_name = 'tcl_gid.tcl';
            file_list = obj.create_file_list(obj.iterations_to_print,obj.file_name,obj.files_folder);

            [output_video_name] = obj.replace_special_character(output_video_name_in);
            [output_photo]=obj.replace_special_character(strcat(obj.files_folder,'\',obj.file_name,'.png'));
            [file_list] = obj.replace_special_character(file_list);

            min_value = -1e-32;

            file_tcl_name_with_path = fullfile(obj.files_folder,file_tcl_name);
            file_path_in = fullfile(pwd,'FEM','PostProcess','Make_Video_characteristic3D.tcl');
            filepath = obj.replace_special_character(file_path_in);

            fid = fopen(file_tcl_name_with_path,'w+');

            fprintf(fid,'GiD_Process PostProcess \n');
            fprintf(fid,'%s\n',['set arg1 "',file_list,'"']);
            fprintf(fid,'%s\n',['set arg2 "',output_video_name,'"']);
            fprintf(fid,'%s\n',['set arg3 "',field2print,'"']);
            fprintf(fid,'%s\n',['set arg4 "',componentfield,'"']);
            fprintf(fid,'%s\n',['set arg5 "',num2str(min_value),'"']);
            fprintf(fid,'%s\n',['set arg6 "',output_photo,'"']);

            fprintf(fid,'%s\n',['source "',filepath,'"']);
            fprintf(fid,'%s\n',['Make_Video_characteristic $arg1 $arg2 $arg3 $arg4 $arg5 $arg6']);
            fprintf(fid,'%s\n',['GiD_Process Mescape Quit']);
            fclose(fid);
            obj.execute_tcl_files(obj.gidPath,file_tcl_name_with_path)
        end

  end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function Snirf2Nirs(rootdir, options)
if ~exist('rootdir','var')
    rootdir = pwd;
end
if ~exist('options','var')
    options = 'sort';
end
rootdir = filesepStandard(rootdir);
files = DataFilesClass(rootdir, 'snirf');
files = files.files;
for ii = 1:length(files)
    if files(ii).IsDir()
        continue
    end
    fname = filesepStandard([files(ii).rootdir, files(ii).name]);
    s = SnirfClass(fname);
    n = NirsClass(s);
    if optionExists(options, 'sort')
        n.SortData();
    end
    [pname, fname] = fileparts(fname);
    fnameNew = [filesepStandard(pname), fname, '.nirs'];
    fprintf('Converting %s to %s\n', fname, fnameNew)
    n.Save(fnameNew);
end

author: BUNPC


text:
%  INTERNAL FUNCTION: checks whether the covariance matrix has converged
%
%  ::
%
%    [is_steady,oldK]=check_steady_state_kalman(is_steady,K,oldK,options,t,no_more_missing)
%
%  Args:
%
%     - **is_steady** [true|false]: flag indicating whether or not the steady
%       state is reached
%     - **K** [matrix]: Kalman gain
%     - **oldK** [matrix]: old Kalman gain
%     - **options** [struct]: options for the Kalman filter
%     - **t** [integer]: iteration number
%     - **no_more_missing** [boolean]: true if there are no more missing
%       observations, false otherwise.
%
%  Returns:
%     :
%
%     - **is_steady** [true|false]: flag indicating whether or not the steady
%       state is reached
%     - **oldK** [matrix]: updated old Kalman gain
%
%

author: jmaih


text:
function [net,opts,x] = invert_nn_dw_pre(net, ref, opts)
% INVERT  Invert a CNN representation
% Parse the input arguments to override the above defaults
% Update the number of iterations using the learning rate if required
if isinf(opts.maxNumIterations)
    opts.maxNumIterations = numel(opts.learningRate) ;
end

net.regu = opts.regu;
net.fmask = opts.fmask;
% The size of the image that we are trying to obtain
x0_size = cat(2,net.normalization.imageSize,opts.numRepeats);
%keyboard

% x0_sigma is computed using a separate dataset.
% This is a useful normalization that helps scale the different terms in the
% optimization.
load('x0_sigma.mat', 'x0_sigma');

% Replicate the feature into a block. This is used for multiple inversions in parallel.
y0 = repmat(ref, [1 1 1 opts.numRepeats]);

% initial inversion image of size x0_size
if  ~isempty(opts.init)
    x= opts.init;
else
    x = randn(x0_size, 'single') ;
    %x= imfilter(x,fspecial('gaussian',50,30));
    %x = ones(x0_size, 'single') ;
    %x = zeros(x0_size, 'single') ;
    x = x / norm(x(:)) * x0_sigma  ;
end
x_momentum = zeros(x0_size, 'single') ;

switch opts.task
    case {0,-1}
        % allow reconstructing a subset of the representation by setting
        % a suitable mask on the features y0
        sf = 1:size(y0,3) ;
        if opts.filterGroup == 1
            sf= vl_colsubset(sf, 0.5, 'beginning') ;
        elseif opts.filterGroup == 2 ;
            sf= vl_colsubset(sf, 0.5, 'ending') ;
        end
        nx = min(opts.neigh, size(y0,2)) ;
        ny = min(opts.neigh, size(y0,1)) ;
        sx = (0:nx-1) + ceil((size(y0,2)-nx+1)/2) ;
        sy = (0:ny-1) + ceil((size(y0,1)-ny+1)/2) ;
        mask = zeros(size(y0), 'single') ;
        mask(sy,sx,sf,:) = 1 ;
        y0_sigma = norm(squeeze(y0(find(mask(:))))) ;
    case {1,-2,2}
        %mask = opts.mask;
        mask = ref;
        y0_sigma = 1;
end

%% Tweak the network by adding a reconstruction loss at the end

% This is saved here just for printing our progress as optimization proceeds
if net.cnn_mode==0
    switch opts.objective
        case 'l2'
            % Add the l2 loss over the network
            ly.type = 'custom' ;
            ly.w = y0 ;
            ly.mask = mask ;
            ly.forward = @nndistance_forward ;
            ly.backward = @nndistance_backward ;
            net.layers{end+1} = ly ;
        case 'l1'
            % The L1 loss might want to use a dropout layer.
            % This is just a guess and hasn't been tried.
            ly.type = 'dropout' ;
            ly.rate = opts.dropout ;
            net.layers{end+1} = ly ;
            ly.type = 'custom' ;
            ly.w = y0 ;
            ly.mask = mask ;
            ly.forward = @nndistance1_forward ;
            ly.backward = @nndistance1_backward ;
            net.layers{end+1} = ly ;
        case 'inner'
            % The inner product loss may be suitable for some networks
            ly.type = 'custom' ;
            ly.w = - y0 .* mask ;
            ly.forward = @nninner_forward ;
            ly.backward = @nninner_backward ;
            net.layers{end+1} = ly ;
        case 'oneclass'
            % The inner product loss may be suitable for some networks
            % maxize the probability of one class
            ly.type = 'custom' ;
            ly.mask =  mask ;
            ly.forward = @oneclass_forward ;
            ly.backward = @oneclass_backward ;
            net.layers{end+1} = ly ;
        otherwise
            error('unknown opts.objective') ;
    end
end

% --------------------------------------------------------------------
function res_ = oneclass_forward(ly, res, res_)
% --------------------------------------------------------------------
% the value for that one class
res_.x = -sum(reshape(res.x.*oneclass(ly.mask,res.x),1,[]));
% --------------------------------------------------------------------
function res = oneclass_backward(ly, res, res_)
% --------------------------------------------------------------------
% positive or negative ?
res.dzdx = -oneclass(ly.mask,res.x);

function out=oneclass(mask,res)
out = reshape(mask,size(res));
%out = repmat(reshape(mask,[1,1,numel(mask)]),[size(res,1),size(res,2),1]);

% --------------------------------------------------------------------
function res_ = nndistance_forward(ly, res, res_)
% --------------------------------------------------------------------
res_.x = nndistance(res.x, ly.w, ly.mask) ;

% --------------------------------------------------------------------
function res = nndistance_backward(ly, res, res_)
% --------------------------------------------------------------------
res.dzdx = nndistance(res.x, ly.w, ly.mask, res_.dzdx) ;

% --------------------------------------------------------------------
function y = nndistance(x,w,mask,dzdy)
% --------------------------------------------------------------------
if nargin <= 3
    d = x - w ;
    y = sum(sum(sum(sum(d.*d.*mask)))) ;
else
    y = dzdy * 2 * (x - w) .* mask ;
end

% --------------------------------------------------------------------
function res_ = l1loss_forward(ly, res, res_)
% --------------------------------------------------------------------
res_.x = ly.w*sum(abs(res.x(:)));

% --------------------------------------------------------------------
function res = l1loss_backward(ly, res, res_)
% --------------------------------------------------------------------
res.dzdx = zeros(size(res.x), 'single');
res.dzdx(res.x > 0) = single(ly.w)*res_.dzdx;
res.dzdx(res.x < 0) = -single(ly.w)*res_.dzdx;
res.dzdx(res.x == 0) = single(0);

% --------------------------------------------------------------------
function res_ = nndistance1_forward(ly, res, res_)
% --------------------------------------------------------------------
res_.x = nndistance1(res.x, ly.w, ly.mask) ;

% --------------------------------------------------------------------
function res = nndistance1_backward(ly, res, res_)
% --------------------------------------------------------------------
res.dzdx = nndistance1(res.x, ly.w, ly.mask, res_.dzdx) ;

% --------------------------------------------------------------------
function y = nndistance1(x,w,mask,dzdy)
% --------------------------------------------------------------------
if nargin <= 3
    d = x - w ;
    y = sum(sum(sum(sum(abs(d).*mask)))) ;
else
    y = dzdy * sign(x - w) .* mask ;
end

% --------------------------------------------------------------------
function res_ = nninner_forward(ly, res, res_)
% --------------------------------------------------------------------
res_.x = nninner(res.x, ly.w) ;

% --------------------------------------------------------------------
function res = nninner_backward(ly, res, res_)
% --------------------------------------------------------------------
res.dzdx = nninner(res.x, ly.w, res_.dzdx) ;

% --------------------------------------------------------------------
function y = nninner(x,w,dzdy)
% --------------------------------------------------------------------
if nargin <= 2
    y = sum(sum(sum(sum(w.*x)))) ;
else
    y = dzdy * w ;
end

author: donglaiw

repo: mNeuron

save_path: github-repos/MATLAB/donglaiw-mNeuron


function plot_frcst_(frcsts,y,time,options)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filippo Ferroni, 6/1/2015
% Revised, 2/15/2017
% Revised, 3/21/2018
% Revised, 8/08/2019

% input : frcsts
% 1st dimension: horizon
% 2nd dimension: variable
% 3rdth dimension: draws
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if length(time) ~= length(y)
    error('Mismatch between time and in sample data (y)');
end
if size(time,2)>size(time,1)
    time=time';
end

hor     = size(frcsts,1);
nvar    = size(frcsts,2);
ndraws  = size(frcsts,3);
nplots  = [ceil(sqrt(nvar)) ceil(sqrt(nvar))];
savefig_yes = 0;
conf_sig    = 0.68;
add_frcst_yes = 0;
add_multiple_bands_yes = 0;
fnam_dir    = '.';
fnam_suffix = 'frcsts';
trasf_yes   =  0;

% retreive the frequency
integerTest = ~mod(time,1);
indx        = find(integerTest==1);
frq         = indx(2)-indx(1);

if frq == 12 % monthly
       timefor = time(end) + 1/12 : 1/12 : time(end) + hor/12;
       time = [time; timefor'];
elseif frq == 4 % quarterly
       timefor = time(end) + 1/4 : 1/4 : time(end) + hor/4;
       time = [time; timefor'];
elseif frq == 1 % annual
       timefor = time(end) + 1 : 1 : time(end) + hor;
       time = [time; timefor'];
elseif frq == 48 % weekly
       timefor = time(end) + 1/48 : 1/48 : time(end) + hor/48;
       time = [time; timefor'];
else
    error('Frequency not defined.')
end
%
% if strmatch(freq,'m') == 1 %#ok<*MATCH2>
%        timefor = time(end) + 1/12 : 1/12 : time(end) + hor/12;
%        time = [time; timefor'];
% elseif strmatch(freq,'q') == 1
%        timefor = time(end) + 1/4 : 1/4 : time(end) + hor/4;
%        time = [time; timefor'];
% elseif strmatch(freq,'a') == 1
%        timefor = time(end) + 1 : 1 : time(end) + hor;
%        time = [time; timefor'];
% else
%     error('You need to provide a frequency: ''m'', ''q'' or ''a''.')
% end
time_start = 1;
time_end = length(time);

if nargin < 4
    disp('You did not provided names for variables.')
    disp('I call them Var 1, Var 2, ... ')
    for v = 1 : nvar
        eval(['varnames{'   num2str(v) '} =  ''Var  ' num2str(v) ''';'])
    end
else
    if isfield(options,'time_start') ==1
        time_start = find(options.time_start==time);
        if isempty(time_start) ==1
            error('''time_start'' is not included in ''T''.')
        end
    end
    if isfield(options,'order_transform') ==1
        trasf_yes       = 1;
        order_trasform  = options.order_transform;
        if length(order_trasform) ~= nvar
            error('Mismatch between the ''order_transform'' size and # of variables.')
        end
    end
    if isfield(options,'varnames') ==1
        varnames = options.varnames;
        if length(varnames) ~= nvar
            error('Mismatch between the # varnames and # of variables to plot')
        end
    else
        disp('You did not provided names for the endogenous variables.')
        disp('I call them Var 1, Var 2, ...')
        for v = 1 : nvar
            eval(['varnames{'   num2str(v) '} =  ''Var  ' num2str(v) ''';'])
        end
    end
    if isfield(options,'nplots') ==1
        nplots = options.nplots;
    end
    if isfield(options,'saveas_strng') ==1
        savefig_yes = 1;
        % setting the names of the figure to save
        fnam_suffix = [ fnam_suffix options.saveas_strng ];
    end
    if isfield(options,'saveas_dir') == 1
        savefig_yes = 1;
        % setting the folder where to save the figure
        fnam_dir = options.saveas_dir;
        if exist(fnam_dir,'dir') == 0
            mkdir(fnam_dir)
        end
    end
    if isfield(options,'add_frcst') ==1
        add_frcst = options.add_frcst;
        if size(add_frcst) ~= [length(time), nvar]
            error('The ''add_frcst'' dimensions must be in-sample + output-of-sample length and the # of variables to plot')
        end
        add_frcst_yes = 1;
    end
    if isfield(options,'conf_sig') ==1
        conf_sig = options.conf_sig;
    end
    if isfield(options,'conf_sig_2') ==1
        if options.conf_sig_2 < conf_sig
            error('Additional confidence bands should be larger than ''options.conf_sig''.')
        end
        add_multiple_bands_yes = 1;
        sort_idx_2   = round((0.5 + [-options.conf_sig_2, options.conf_sig_2, 0]/2) * ndraws);
    end
end

nfigs  = ceil(length(varnames)/( nplots(1)*nplots(2)) );
nplots = repmat(nplots,nfigs,1);
for j=1:size(nplots,1),
    nbofplots(j)=nplots(j,1)*nplots(j,2);
end

ntotplots = sum(nbofplots);
if ntotplots<length(varnames),
    nfigplus = ceil((length(pplotvar)-ntotplots)/nbofplots(end));
    lastrow=nplots(end,:);
    lastrow=repmat(lastrow,nfigplus,1);
    nplots = [nplots;lastrow];
    nbofplots = [nbofplots repmat(nbofplots(end),1,nfigplus)];
end

if ndraws > 1
    frcsts_ = nan(length(time),nvar,ndraws);
    for kk = 1 : ndraws
        frcsts_(:,:,kk) = [y; frcsts(:,:,kk)];
    end
    frcsts = frcsts_;

    sort_idx   = round((0.5 + [-conf_sig, conf_sig, 0]/2) * ndraws);

    if trasf_yes ==1
        frcsts_ = nan(size(frcsts));
        for var = 1 : nvar
            if order_trasform(var) == 1 % period over period
                frcsts_(2:end,var,:) = (frcsts(2:end,var,:) - frcsts(1:end-1,var,:)) ;

            elseif order_trasform(var) == 100 % percentage period over period
                frcsts_(2:end,var,:) = 100*(frcsts(2:end,var,:) - frcsts(1:end-1,var,:));

            elseif order_trasform(var) == 12 % percentage 12 periord over 12 period (year over year % change f  or monthly data)
                frcsts_(13:end,var,:) = 100*(frcsts(13:end,var,:) - frcsts(1:end-12,var,:));

            elseif order_trasform(var) == 4 % percentage 4 periord over 4 period (year over year % change for quarterly data)
                frcsts_(5:end,var,:) = 100*(frcsts(5:end,var,:) - frcsts(1:end-4,var,:));

            elseif order_trasform(var) == 400
                frcsts_(2:end,var,:) = 400*(frcsts(2:end,var,:) - frcsts(1:end-1,var,:));

            elseif order_trasform(var) == 1200
                frcsts_(2:end,var,:) = 1200*(frcsts(2:end,var,:) - frcsts(1:end-1,var,:));

            else
                frcsts_(:,var,:) = frcsts(:,var,:);
            end
        end
        frcsts = frcsts_;
    end

    frcsts_sort   = sort(frcsts,3);
    if sort_idx(1) == 0
        sort_idx(1) = 1;
        warning('Bands not reliable. You have too few draws.')
    end
    irf_Median = squeeze(frcsts_sort(:, :, sort_idx(3) ));
    irf_low    = squeeze(frcsts_sort(:, :, sort_idx(1) ));
    irf_up     = squeeze(frcsts_sort(:, :, sort_idx(2) ));
    if  add_multiple_bands_yes == 1
        irf_low_low  = squeeze(frcsts_sort(:, :, sort_idx_2(1) ));
        irf_up_up    = squeeze(frcsts_sort(:, :, sort_idx_2(2) ));
    end

else
    irf_Median = [y; frcsts];

    if trasf_yes ==1
        irf_Median_ = nan(size(irf_Median));
        for var = 1 : nvar
            if order_trasform(var) == 1
                irf_Median_(2:end,var) = diff(irf_Median(:,var)) ;

            elseif order_trasform(var) == 12 % percentage 12 periord over 12 period (year over year % change f  or monthly data)
                irf_Median_(13:end,var) = 100*(irf_Median(13:end,var) - irf_Median(1:end-12,var));

            elseif order_trasform(var) == 4 % percentage 4 periord over 4 period (year over year % change for quarterly data)
                irf_Median_(5:end,var) = 100*(irf_Median(5:end,var) - irf_Median(1:end-4,var));

            elseif order_trasform(var) == 100
                irf_Median_(2:end,var) = 100*diff(irf_Median(:,var));

            elseif order_trasform(var) == 400
                irf_Median_(2:end,var) = 400*diff(irf_Median(:,var));

            elseif order_trasform(var) == 1200
                irf_Median_(2:end,var) = 1200*diff(irf_Median(:,var));

            else
                irf_Median_(:,var) = irf_Median(:,var);
            end
        end
        irf_Median = irf_Median_;
    end
    irf_low    = irf_Median;
    irf_up     = irf_Median;

end

if trasf_yes ==1 && add_frcst_yes ==1
    add_frcst_ = nan(size(add_frcst));
    for var = 1 : nvar
        if order_trasform(var) == 1
            add_frcst_(2:end,var) = diff(add_frcst(:,var)) ;

        elseif order_trasform(var) == 12 % percentage 12 periord over 12 period (year over year % change for monthly data)
            add_frcst_(13:end,var) = 100*(add_frcst(13:end,var) - add_frcst(1:end-12,var));

%         elseif order_trasform(var) == 4 % percentage 4 periord over 4 period (year over year % change for quarterly data)
%             add_frcst_(4:end,var) = 100*(add_frcst(4:end,var) - add_frcst(1:end-3,var));

        elseif order_trasform(var) == 4 % percentage 4 periord over 4 period (m over m % change for monthly data or YoY for Q data)
            add_frcst_(5:end,var) = 100*(add_frcst(5:end,var) - add_frcst(1:end-4,var));

        elseif order_trasform(var) == 100
            add_frcst_(2:end,var) = 100*diff(add_frcst(:,var));

        elseif order_trasform(var) == 400
            add_frcst_(2:end,var) = 400*diff(add_frcst(:,var));

        elseif order_trasform(var) == 1200
            add_frcst_(2:end,var) = 1200*diff(add_frcst(:,var));

        else
            add_frcst_(:,var) = add_frcst(:,var);
        end
    end
    add_frcst = add_frcst_;
end

jplot = 0;
jfig  = 0;

index_time = time_start:time_end;

for var= 1: nvar

    if jplot==0,
        figure('name',['Forecasts'] );
        jfig=jfig+1;
    end

    jplot=jplot+1;
    subplot(nplots(1),nplots(2),jplot)

    if add_multiple_bands_yes == 1

        h = area([time(index_time)],[irf_low_low(index_time,var),...
            irf_low(index_time,var) - irf_low_low(index_time,var),...
            irf_up(index_time,var) - irf_low(index_time,var),...
            irf_up_up(index_time,var) - irf_up(index_time,var)]);%,'FaceColor',[.85 .85 .85]);
        set(h(4),'FaceColor',[.95 .95 .95])
        set(h(3),'FaceColor',[.85 .85 .85])
        set(h(2),'FaceColor',[.95 .95 .95])
        set(h(1),'FaceColor',[1 1 1])
        set(h,'linestyle','none')
        hold on
        min_=min(irf_low_low(index_time,var));
        max_=max(irf_up_up(index_time,var));
    else
        h = area([time(index_time)],[irf_low(index_time,var),...
            irf_up(index_time,var) - irf_low(index_time,var)]);%,'FaceColor',[.85 .85 .85]);
        set(h(2),'FaceColor',[.85 .85 .85])
        set(h(1),'FaceColor',[1 1 1])
        set(h,'linestyle','none')
        min_=min(irf_low(index_time,var));
        max_=max(irf_up(index_time,var));
    end
    hold on;
    plot(time(index_time),irf_Median(index_time,var),'k');
    hold on;
    plot(time(index_time),zeros(length(index_time),1),'k');
    hold on;
    axis tight
    ylim([min_,max_]);
    if add_frcst_yes == 1
        plot(time(index_time),add_frcst(index_time,var),'b','LineWidth',2);
        ylim([min([min_,min(add_frcst(index_time,var))]) , max([max_,max(add_frcst(index_time,var))])  ] );
    end
    title(varnames{var})
    set(gcf,'position' ,[50 50 800 650])
    if jplot==nbofplots(jfig) || var==length(varnames)
        %legend(legenda);
        if savefig_yes == 1
            STR_RECAP = [ fnam_dir '/' fnam_suffix  '_' int2str(jfig)];
            saveas(gcf,STR_RECAP,'fig');
            if strcmp(version('-release'),'2022b') == 0
                saveas(gcf,STR_RECAP,'eps');
                % savefigure_pdf(STR_RECAP);
                savefigure_pdf([STR_RECAP '.pdf']);
            end
        end
        jplot=0;
    end
end


text:
function r8mat_transpose_print_test ( )

%*****************************************************************************80
%
%% R8MAT_TRANSPOSE_PRINT_TEST tests R8MAT_TRANSPOSE_PRINT;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  m = 7;
  n = 12;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_TRANSPOSE_PRINT_TEST\n' );
  fprintf ( 1, '  R8MAT_TRANSPOSE_PRINT prints a R8MAT,\n' );
  fprintf ( 1, '  transposed.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix row order M =    %d\n', m );
  fprintf ( 1, '  Matrix column order N = %d\n', n );
%
%  Set the matrix.
%
  for i = 1 : m
    for j = 1 : n
      a(i,j) = i * 100 + j;
    end
  end

  r8mat_transpose_print ( m, n, a, '  The transposed matrix A:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [] = aps_weather_model_SAR(model_type)
% [] = aps_weather_model_SAR()
% Script to load weather model data and compute SAR delays
% The DEM file can be .grd file or other. If the latter, the DEM should have an asociated
% ".rsc" file, with the same filename as the DEM. The ".rsc" files should
% contain a WIDTH, LENGTH, X_FIRST, Y_FIRST, X_STEP, Y_STEP and optional a
% FORMAT string. The weather model data is assumed to be structured in date_before(d,:) folders.
%
%
% INPUTS:
% demfile               Full path to the DEM file. The DEM needs to me in meters.
% xlims                 Limits in the x-direction, either in degrees
% ylims                 Limits in the y-direction, either in degrees
% demnull               The value for no DEM data, default is -32768.
% smpres                The output resolution, either in degrees
%                       Units needs to be consistend with xlims and ylims.
%
% OUTPUTS:
% It will give the computed ZENITH hydrostatic and wet delay map in cm for the selected region.
%
%     Copyright (C) 2015  Bekaert David - University of Leeds
%     Email: eedpsb@leeds.ac.uk or davidbekaert.com
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
%
% Modified from Richard Walters - Oxford / Leeds University 2012-2013
% Modifcations:
% DB 	10/2013		Convert script to a function, integrate in aps_toolbox,
%                   add syntax. Retrieve parameters automatic.
% DB    10/2013     Output the combined delay map.
% DB    10/2013     Increased processing efficiency
% DB    10/2013     add compatibility with Doris construct_dem.sh script
% DB    10/2013     Bug fix incase the the netcdf is split over two days
% DB    11/2013     Change the computation f saturated pressure to be
%                   identical witt how it was computed in ERA-I model.
% DB    11/2013     Update filenaming to Hydrostatic as its not dry delay.
%                   Compute the hydrostatic delay directly from int(P/T)
% DB    01/2014     Fix bug for the ERA-I longitude boundary, being set to small.
% DB    01/2014     Allow for tiled ERA-I data
% HB-DB 02/2014     Fix to include ECMWF with scaling and offset
% DB    02/2014     Include a parameter for ECMWF data website
% DB    03/2014     Include an extra check for the DEM grd file and the
%                   selected crop
% DB 	03/2014     Save individual time stamp SAR date_before(d,:) processing
% DB    07/2014     Redefine meris_lat(lon)_range to region_lat(lon)_range
% DB    08/2014     Check if netcdf exist (ERA-I is not an operation service)
% DB    11/2014     Include option to output the 3D delays for each date_before(d,:)
% HW    02/2015     check dem file format automatically
% DB 	06/2015     Bug fix for plotting the support information
% DB    06/2015     Fix typo in error message
% DB    11/2015     Branch of the DEM in to seperate function same for all codes
% DB    11/2015     Add multicore option from matlab
% DB 	02/2016     Close netcdf files
% DB    04/2016     Branch into weather model script and include merra too
% SSS   04/2016     Clear variables such memory need is reduced
% DB    07/2016     redefine hydrostatic delay to be based on surface pressure.
% DB    08/2016     Uncomment a keyboard
% DB    07/2016     expand to include ERA5 test model data
% KM    02/2018     expand to include NARR model data

fig_test = 0;           % when 1 show the dem as debug figure
save_3D_delays = 0;     % When 1 saves the tropopsheric delays for each x,y and with height

if nargin<1
    error('Give at least the model_type: era, era5, narr, merra, or merra2')
end
% change to lower caps for saving and filename generation consistency
model_type = lower(model_type);

%% Constants
% Parameters from Jolivet et al 2011, GRL
Rd = 287.05;            % [j/kg/K] specific gas constants dry air
Rv = 461.495;           % [j/kg/K] water vapour
k1 = 0.776;             % [K/Pa]
k2 = 0.716;             % [K/Pa]
k3 = 3.75e3;            % [K^2/Pa]
coeff.Rd = Rd;
coeff.Rv = Rv;

%% Defaults
zref = 15000;       % zref for integration- tropopause
zincr = 15;         % z increment spacing for vertical interpolation before integration
vertres = 5;        % vertical resolution for comparing dem to vert profiles of delay

%% output file names
% output file for the DEM and look angles
smpdem = 'dem_smp.xyz';

% getting the variables from the parms_aps file
stamps_processed = getparm_aps('stamps_processed',1);

%%% defaults for the weather models. If not applicable it will be changed below for the specific model.
if strcmp(model_type,'narr')
    timelist_model = ['0000' ;'0300'; '0600' ; '0900'; '1200' ;'1500'; '1800' ;'2100'; '0000'];
    model_lag = 8*7;    % days
else
    timelist_model= ['0000' ; '0600' ; '1200' ; '1800' ; '0000'];       % the time interval the model is outputed
    model_lag = 0; % ? check lags
end
era_data_type = [];                                                 % the weather model data type for ERA only.

%%% Updating specific weather model information
if strcmpi(model_type,'era')
    weather_model_datapath = getparm_aps('era_datapath',1);
    era_data_type = getparm_aps('era_data_type');         % the datatype of the model either BADC or ERA
elseif strcmpi(model_type,'merra') || strcmpi(model_type,'merra2')
    weather_model_datapath = getparm_aps('merra_datapath',1);
elseif strcmpi(model_type,'narr')
    weather_model_datapath = getparm_aps('narr_datapath',1);

else
    error(['weather model type not supported, either: wrf, era, narr, merra for now'])
end

lambda = getparm_aps('lambda',1)*100;                       % radar wavelength in cm
demfile = getparm_aps('demfile',1);
demnull = getparm_aps('dem_null',1);
UTC_sat =  getparm_aps('UTC_sat',1);
ifgday_matfile = getparm_aps('ifgday_matfile',1);
ifgs_dates = load(ifgday_matfile);

% loading the data
if strcmp(stamps_processed,'y')
   dates = ifgs_dates.day;
   load psver
   fprintf('Stamps processed structure \n')
else
    psver = 2;
    ifgs_dates = ifgs_dates.ifgday;
    dates = reshape(ifgs_dates,[],1);
    dates = unique(dates);
    dates = datenum(num2str(dates),'yyyymmdd');
end

% getting the dates
n_dates = length(dates);

%% Compute and resample DEM
[dem,xmin,xmax,ymin,ymax,smpres,nncols,nnrows] = get_DEM;

% the region which is cropped from the ERA data and used to make the interpolation.
% Should be  larger than the region to which the delay is computed
lonmin = floor(xmin)-1;
lonmax= ceil(xmax)+1;
latmin = floor(ymin)-1;
latmax = ceil(ymax)+1;

% setting the maximum height of the DEM to limit the range at which ERA-I
% needs to be interpolated to
maxdem = ceil(max(max(dem))/100)*100+50; % max height of dem in metres

fprintf(['Interpolate to a maximum dem height of ', num2str(maxdem) ,' m\n'])

%% Compute based on Satellite pass which weather model outputs that will be used
[time_before,time_after, date_before, date_after,f_before,f_after] = aps_weather_model_times(timelist_model,dates,UTC_sat,model_lag);

%% generating a file
[modelfile_before,modelfile_after] = aps_weather_model_filenames(model_type,time_before,time_after,date_before, date_after,weather_model_datapath);

%% performing the calucluation for each date
for d = 1:n_dates

    % the save filenames
    outfile = [weather_model_datapath filesep date_before(d,:) filesep date_before(d,:) '_ZWD.xyz'];
    hydroutfile = [weather_model_datapath filesep date_before(d,:) filesep date_before(d,:) '_ZHD.xyz'];

    no_data = 0;
    for kk = 1:2
        if kk == 1
            file = modelfile_before(d,:);
        end
        if kk == 2
            file = modelfile_after(d,:);
        end

        % check
        if exist(file,'file')~=2
            no_data = no_data+1;
        else

            % loading the weather model data
            if  strcmpi(model_type,'era')
                 [ Temp,e,Geopot,P,longrid,latgrid,xx,yy,lon0360_flag] =  aps_load_era(file,era_data_type) ;
            elseif  strcmpi(model_type,'era5')
                 [ Temp,e,Geopot,P,longrid,latgrid,xx,yy,lon0360_flag] =  aps_load_era(file,era_data_type) ;
            elseif strcmpi(model_type,'merra') || strcmpi(model_type,'merra2')
                 [ Temp,e,Geopot,P,longrid,latgrid,xx,yy,lon0360_flag] =  aps_load_merra(file,model_type,coeff) ;
            elseif strcmpi(model_type,'narr')
                 [ Temp,e,Geopot,P,longrid,latgrid,xx,yy,lon0360_flag] =  aps_load_narr(file,model_type,coeff) ;
            end
%             longrid=longrid';latgrid=latgrid';
            % verify and cope with NAN's
            [Temp,e,Geopot,P,longrid,latgrid] =  aps_weather_model_nan_check(Temp,e,Geopot,P,longrid,latgrid) ;

            % define weather model grid nodes
            latlist = reshape(latgrid(:,:,1),[],1);
            lonlist = reshape(longrid(:,:,1),[],1);
            xlist = reshape(xx,[],1);
            ylist = reshape(yy,[],1);

            % Limit weather model to those grid points around the user-defined InSAR box
            % Getting the weather model resolution
            lat_res = abs(diff(unique(latgrid)))*1.5;
            lat_res = lat_res(1);
            lon_res = abs(diff(unique(longrid)))*1.5;
            lon_res = lon_res(1);

            % make sure the weather model and the lonlat specified by user
            % is defined in the same way. Weatehr models can be [0 360]
            % degrees. User can be [-180 180] unlikely [0 360]
            if strcmpi(lon0360_flag,'y')
                % fprintf('This needs to be defined better and some checks are needed \n')
                % not all pixels should be shifted.

                if xmin<0
                   xmin= xmin+360;
                end
                if xmax<0
                   xmax = xmax +360;
                end
            end

            % generation a plot
            if fig_test ==1 & d==1 & kk==1
                fontsize = 15;
                hfig = figure('name',[model_type ' weather model nodes and InSAR region'],'position', [ 200   243   610   603]);
                % for the legend generation
                plot(mean([xmin xmax]),mean([ymax ymin]),'wo','markeredgecolor','k','markerfacecolor','w','markersize',15)
                hold on
                plot(mean([xmin xmax]),mean([ymax ymin]),'r.')
                hold on
                plot(mean([xmin xmax]),mean([ymax ymin]),'r-','linewidth',2)

                imagesc([xmin xmax],[ymax ymin],dem)
                cc=colorbar;
                view(0,90)

                axis xy
                ylabel(cc,'Topography [m]','fontsize',fontsize)
                hold on
                plot(lonlist,latlist,'wo','markeredgecolor','k','markerfacecolor','w','markersize',15)
                hold on
                plot([xmin xmin xmax xmax xmin],[ymin ymax ymax ymin ymin],'r-','linewidth',2)

                title([model_type ' weather model nodes'],'fontsize',fontsize)
                set(gca,'fontsize',fontsize)
                axis equal
                axis tight

                legend('location','northoutside',[model_type ' weather model nodes'],['Used ' model_type ' weather model nodes'],'InSAR box')
            end

            % making the weather model grid slightly larger than the InSAR
            % bounding box. Will use the weather model resolution for this
            % to make sure an extra grid node is included.
            ix = find(ymin-lat_res<= latlist & latlist <= ymax+lat_res  & xmin-lon_res<= lonlist & lonlist<= xmax+lon_res) ;
            xlist = xlist(ix);
            ylist = ylist(ix);
            latlist = latlist(ix);
            lonlist = lonlist(ix);
            numy = length(unique(latlist));
            numx = length(unique(lonlist));
            ulatlist = unique(latlist);
            ulonlist = unique(lonlist);
            uxlist = unique(xlist);
            uylist = unique(ylist);

            if fig_test ==1 & d==1 & kk==1
                figure(hfig)
                hold on
                plot(lonlist,latlist,'r.')
                xlim([xmin-4*lon_res  xmax+4*lon_res])
                ylim([ymin-4*lat_res  ymax+4*lat_res])
                str='';
                while ~strcmpi(str,'y') && ~strcmpi(str,'n')
                    fprintf(['Do the red nodes extend (just) outside the red InSAR rectangle? \n'])
                    str = input('Continue? [y: for yes, n: no] \n','s');
                end
                if strcmpi(str,'n')
                    error('Check your lon lat crop, otherwize extend area of downlaoded weather model data.')
                end
            end

            clear lon_res lat_res ylist xlist %%%SSS 4/16

            % saving the information for support plotting
            eval([model_type '.' model_type '_lonlat =[lonlist latlist];']);                            % era.era_lonlat = [lonlist latlist];
            eval([model_type '.region =[[xmin xmin xmax xmax xmin]'' [ymin ymax ymax ymin ymin]''];']); % era.region = [[xmin xmin xmax xmax xmin]' [ymin ymax ymax ymin ymin]'];
            deminfo.xmin = xmin;
            deminfo.xmax = xmax;
            deminfo.ymax = ymax;
            deminfo.ymin = ymin;
            deminfo.dem = dem;
            eval([model_type '.deminfo =deminfo;']);                                                    % era.deminfo = deminfo;
            clear deminfo
            % checking if the file already exist. Yes append otherwiuze create it
            if exist('tca_support.mat','file')==2
                eval(['save(''tca_support.mat'',''-append'',''' model_type ''');']);                    % save('tca_support.mat','-append','era')
            else
                eval(['save(''tca_support.mat'',''' model_type ''');']);                                % save('tca_support.mat','era')
            end
%             if fig_test ==1 & d==1 & kk==1
%                     if exist(['aps_' model_type],'dir')~=7
%                         mkdir(['aps_' model_type]);
%                     end
%                     print(hfig,'-dpng',['aps_' model_type  filesep model_type '_datapoints.png'])
%                     print(hfig,'-depsc',['aps_' model_type filesep model_type '_datapoints.eps'])
%             end
            eval(['clear ' model_type]);

            % Convert Geopotential to Geopotential Height and then to Geometric Height
            g0 = 9.80665;
            % Calculate Geopotential Height, H
            H = Geopot./g0;

            % map of g with latitude
            g = 9.80616.*(1 - 0.002637.*cosd(2.*latgrid) + 0.0000059.*(cosd(2.*latgrid)).^2);
            % map of Re with latitude
            Rmax = 6378137;
            Rmin = 6356752;
            Re = sqrt(1./(((cosd(latgrid).^2)./Rmax^2) + ((sind(latgrid).^2)./Rmin^2)));

            % Calculate Geometric Height, Z
            Z = (H.*Re)./(g/g0.*Re - H);

            % Find middle of scene to work out glocal and Rlocal for use later
            midx = round(mean(uxlist));
            midy = round(mean(uylist));
            glocal = g(midy,midx,1);
            Rlocal = Re(midy,midx,1);

            cdslices = maxdem/vertres +1;
            cdstack = zeros(numy,numx,cdslices);
            cdstack_dry = zeros(numy,numx,cdslices);
            cdstack_wet = zeros(numy,numx,cdslices);

            XI=(0:zincr:zref)';
            gh = glocal.*(Rlocal./(Rlocal+XI)).^2; %gravity with height for XI height range

            % Interpolate Temp P and e from 0:20:15000 m
            % then integrate using trapz to estimate delay as function of height
            for i = 1:numx;
                for j = 1:numy;
                    xn = uxlist(i);
                    yn = uylist(j);

                    %interpolate at res zincr before doing integration
                    X=double(squeeze(Z(yn,xn,:)));
                    Ye=double(squeeze(e(yn,xn,:)));
                    YeI = interp1(X,Ye,XI,'spline')*100; %change from hPa to Pa

                    YP=double(squeeze(P(yn,xn,:)));
                    YPI = interp1(X,YP,XI,'spline')*100; %change from hPa to Pa

                    YT=double(squeeze(Temp(yn,xn,:)));
                    YTI = interp1(X,YT,XI,'spline');

                    tmp1 = ((k2-(Rd*k1/Rv)).*YeI./YTI + k3.*YeI./(YTI.^2));
                    Lw = (10^-6).*-1*flipud(cumtrapz(flipud(XI),flipud(tmp1)));
                    % L is zenith delay one way in metres
                    % tmp2 = k1.*YPI./YTI;
                    %Ld = (10^-6).*-1*flipud(cumtrapz(flipud(XI),flipud(tmp2)));             % This is using P/T expression (Hanssen, 2001)
                    gm = glocal;
                    Ld = (10^-6).*((k1*Rd/gm).*(YPI - YPI(zref/zincr +1)));                 % This is P0 expression (Hanssen, 2001)

                    % Interpolate important part (i.e. total delay at elevations
                    % less than maxdem) at high res i.e. vertres, and put in cdstack.
                    cdI=(0:vertres:maxdem)';
                    LdI=interp1(XI,Ld,cdI,'spline');
                    LwI=interp1(XI,Lw,cdI,'spline');

                    %cdstack(j,i,:) = LsI;
                    cdstack_dry(j,i,:) = LdI;
                    cdstack_wet(j,i,:) = LwI;

                    if save_3D_delays==1
                        cdI3D=(0:100:max(XI))';
                        LdI3D=interp1(XI,Ld,cdI3D,'spline');
                        LwI3D=interp1(XI,Lw,cdI3D,'spline');

                        %cdstack(j,i,:) = LsI;
                        cdstack_dry3D(j,i,:) = LdI3D;
                        cdstack_wet3D(j,i,:) = LwI3D;
                        clear LdI3D LwI3D
                    end
                end
            end
            clear uxlist uylist ulonlist ulatlist %%%SSS 4/16

            % Interpolate each cdstack layer onto a grid given by the DEM extents
            % in UTM m.
            xsmpres = (xmax-xmin)/nncols;
            ysmpres = (ymax-ymin)/nnrows;
            [xi,yi] = meshgrid(xmin+0.5*xsmpres:xsmpres:xmax-0.5*xsmpres,ymin+0.5*ysmpres:ysmpres:ymax-0.5*ysmpres);

            ix_temp = diff(lonlist);
            ix_temp = find(ix_temp~=0);
            ix_temp = ix_temp(1);
            lonlist_matrix = reshape(lonlist,ix_temp,[]) ;
            latlist_matrix = reshape(latlist,ix_temp,[]) ;
            clear lonlist latlist %%%SSS 4/16

            % saving the outputs
            if save_3D_delays==1
                if kk==1
                    clear hgt lon lat dry1 dry2 wet1 dem_temp hgt_topo %%%SSS 4/16
                    hgt = cdI3D;
                    lon = [lonlist_matrix];
                    lat = latlist_matrix;
                    dry1 = cdstack_dry3D;
                    wet1 = cdstack_wet3D;

                    % also give the station topography
                    dem_temp = dem;
                    dem_temp(isnan(dem_temp))=0;
                    hgt_topo = griddata(xi,yi,dem_temp,lon,lat,'linear');
                    clear dem_temp

                elseif kk==2
                    clear dry2 wet2
                    dry2 = cdstack_dry3D;
                    wet2 = cdstack_wet3D;
                end
            end

            clear cdstack_interp_dry %%%SSS 4/16
            cdstack_interp_dry = zeros(nnrows,nncols,cdslices);
            parfor n = 1:cdslices
                newslice = interp2(lonlist_matrix,latlist_matrix,cdstack_dry(:,:,n),xi,yi,'linear');
                cdstack_interp_dry(:,:,n)= flipud(newslice); % flipud due to ypositive downpage for matlab figs, and ypositive uppage for utmy
            end

            clear cdstack_interp_wet %%%SSS 4/16
            cdstack_interp_wet = zeros(nnrows,nncols,cdslices);
            parfor n = 1:cdslices
               newslice = interp2(lonlist_matrix,latlist_matrix,cdstack_wet(:,:,n),xi,yi,'linear');
               cdstack_interp_wet(:,:,n)= flipud(newslice); % flipud due to ypositive downpage for matlab figs, and ypositive uppage for utmy
            end
            clear lonlist_matrix latlist_matrix %%%SSS 4/16
            % keeping the coordinates in the same grid as the data
            xi = flipud(xi);
            yi = flipud(yi);

            % Pull out delays from cdstack layers that match dem heights
            clear wetcorrection hydrcorrection rounddem %%%SSS 4/16
            wetcorrection = ones(nnrows,nncols);
            hydrcorrection = ones(nnrows,nncols);
            rounddem = round(dem/vertres);
            rounddem(dem < 0)=0;

            % cope with the case that NaN are present. This can be sea-level
            rounddem(isnan(dem))=0;

            for i=1:nnrows
                for j=1:nncols
                    wetcorrection(i,j) = cdstack_interp_wet(i,j,rounddem(i,j)+1);
                end
            end

            for i=1:nnrows
                for j=1:nncols
                    hydrcorrection(i,j) = cdstack_interp_dry(i,j,rounddem(i,j)+1);
                end
            end

            if kk==1
                wetcorrection1 = wetcorrection;
                drycorrection1 = hydrcorrection;
            end
            if kk==2
                wetcorrection2 = wetcorrection;
                drycorrection2 = hydrcorrection;
            end
            clear wetcorrection hydrcorrection

        end
    end

    if sum(no_data)==0
        % note that this is a one way Zenith delay and not a slant delay.
        % Units are in cm

        % saving individual estimates based on the time-stamp
        outfile_wet_before = [weather_model_datapath filesep date_before(d,:) filesep date_before(d,:) '_ZWD_before.xyz'];
        outfile_wet_after = [weather_model_datapath filesep date_before(d,:) filesep date_before(d,:) '_ZWD_after.xyz'];
        outfile_dry_before = [weather_model_datapath filesep date_before(d,:) filesep date_before(d,:) '_ZHD_before.xyz'];
        outfile_dry_after = [weather_model_datapath filesep date_before(d,:) filesep date_before(d,:) '_ZHD_after.xyz'];

        fidwet_before = fopen(outfile_wet_before,'w');
        data_write =  [reshape(xi,[],1) reshape(yi,[],1) reshape(wetcorrection1,[],1)]';
        tally = fwrite(fidwet_before,data_write,'double');
        fclose(fidwet_before);
        clear data_write tally %%%SSS 4/16
        fidwet_after = fopen(outfile_wet_after,'w');
        data_write =  [reshape(xi,[],1) reshape(yi,[],1) reshape(wetcorrection2,[],1)]';
        tally = fwrite(fidwet_after,data_write,'double');
        fclose(fidwet_after);
        clear data_write tally %%%SSS 4/16
        fiddry_before = fopen(outfile_dry_before,'w');
        data_write =  [reshape(xi,[],1) reshape(yi,[],1) reshape(drycorrection1,[],1)]';
        tally = fwrite(fiddry_before,data_write,'double');
        fclose(fiddry_before);
        clear data_write tally %%%SSS 4/16
        fiddry_after = fopen(outfile_dry_after,'w');
        data_write =  [reshape(xi,[],1) reshape(yi,[],1) reshape(drycorrection2,[],1)]';
        tally = fwrite(fiddry_after,data_write,'double');
        fclose(fiddry_after);
        clear data_write tally %%%SSS 4/16

        % saving the outputs
        if save_3D_delays==1
            wet = wet1*f_before(d) +  wet2*f_after;
            dry = dry1*f_before(d) +  dry2*f_after;
            save([weather_model_datapath filesep date_before(d,:) filesep date_before(d,:) '_3D.mat'],'dry','wet','hgt','lon','lat','hgt_topo')
            clear wet dry hgt dry1 dry2 wet1 wet2 %%%SSS 4/16
        end

        % Output wet correction
        wetcorrection = wetcorrection1*f_before(d) +  wetcorrection2*f_after(d);
        clear wetcorrection1 wetcorrection2
        wetcorrection = wetcorrection*100;                 % delay in [cm]
        fid = fopen(outfile,'w');
        data_write =  [reshape(xi,[],1) reshape(yi,[],1) reshape(wetcorrection,[],1)]';
        tally = fwrite(fid,data_write,'double');
        fclose(fid);
        clear data_write tally wetcorrection %%%SSS 4/16

        % Output hydrostatic correction
        hydrcorrection = drycorrection1*f_before(d) +  drycorrection2*f_after(d);
        clear drycorrection1 drycorrection2
        hydrcorrection = hydrcorrection*100;                 % delay in [cm]
        fid = fopen(hydroutfile,'w');
        data_write =  [reshape(xi,[],1) reshape(yi,[],1) reshape(hydrcorrection,[],1)]';
        tally = fwrite(fid,data_write,'double');
        fclose(fid);
        clear data_write tally hydrcorrection %%%SSS 4/16

        fprintf([num2str(d) ' completed out of ' num2str(n_dates) '\n' ])

    else
        fprintf([num2str(d) ' completed out of ' num2str(n_dates) ' (NO DATA)\n' ])
    end
end

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
clear all, close all, clc
figpath = '../FIGURES/'; mkdir(figpath)
datapath = '../DATA/EX_LOTKA_Dependencies/';
addpath('../utils');

%% Paramaters
NARXtraining = 'trainlm'; %'trainlm',trainbr
InputSignalType = 'sphs'; %sphs,sine2

% Noise
% Ntrain_vec = 1000; %1000;
% eta_vec = [0.01 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5];
% Nr = 50;

% Training length
Ntrain_vec = [5:15,20:5:95,100:100:1000];%,1500:500:3000];
eta_vec = 0;
Nr = 1;

N_ETA = length(eta_vec);
N_LENGTHS = length(Ntrain_vec);
Nmodels = N_ETA*N_LENGTHS;

Nt = 1000; % length(tv);

LOG_SCALE = 1;
%% Load all models
ResultsALL(3) = struct('err', [], 'RelErr', [], 'Ntrain_vec', [], 'DataTrain', [], 'DataValid', [], 'Ttraining', []);

ModelName = 'DMDc';
datapath1 = [datapath,'DMDc/'];
filename = fullfile(datapath1,['EX_LOTKA_SI_',ModelName,'_',InputSignalType,'_Nevol','_STATS.mat']);
load(filename);
ResultsALL(1) = Results;

ModelName = 'SINDYc';
datapath1 = [datapath,'SINDYc/'];
filename = fullfile(datapath1,['EX_LOTKA_SI_',ModelName,'_',InputSignalType,'_Nevol','_STATS.mat']);
load(filename);
ResultsALL(2) = Results;

ModelName = 'NARX';
datapath1 = [datapath,'NARX/',NARXtraining,'/'];
filename = fullfile(datapath1,['EX_LOTKA_SI_',ModelName,'_',InputSignalType,'_Nevol','_STATS.mat']);
load(filename);
ResultsALL(3) = Results;

%% Time series
clear ph
figure, hold on, box on
ccolors = get(gca,'colororder');
iM = 1;
ph(1) = plot(ResultsALL(iM).DataTrain.t, ResultsALL(iM).DataTrain.x(:,1),'-','Color',ccolors(1,:),'LineWidth',1);
ph(2) = plot(ResultsALL(iM).DataTrain.t, ResultsALL(iM).DataTrain.x(:,2),'-','Color',ccolors(2,:),'LineWidth',1);

xlim([0 max(ResultsALL(iM).DataTrain.t)]), ylim([0 100])
xlabel('Time');
ylabel('Population size')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 150])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_TS','_noleg.eps']);
l1 = legend(ph,'Predator','Prey');
set(l1,'Location','NorthEast')
l1.Position = [l1.Position(1)-0.01 l1.Position(2)+0.02 l1.Position(3)-0.01 l1.Position(4)];
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_TS','.eps']);

clear ph
figure, hold on, box on
ccolors = get(gca,'colororder');
iM = 1;
plot(ResultsALL(iM).DataTrain.x(:,1), ResultsALL(iM).DataTrain.x(:,2),'-k','LineWidth',1);
axis equal, axis tight, %axis off
set(gca,'xtick', [], 'ytick', [])
% xlim([0 max(ResultsALL(iM).DataTrain.t)]), ylim([0 100])
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 150])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_PhasePlot','.eps']);

clear ph
figure, hold on, box on
ccolors = get(gca,'colororder');
iM = 1;
plot(ResultsALL(iM).DataTrain.t, ResultsALL(iM).DataTrain.u,'-k','LineWidth',1);
axis equal, axis tight, %axis off
set(gca,'xtick', [], 'ytick', [])
% xlim([0 max(ResultsALL(iM).DataTrain.t)]), ylim([0 100])
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 50])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_Actuation','.eps']);

%% Mean squared error
clear ph
% Prediction over training/validation
figure, hold on, box on
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
for iM = 1:3
    ph(iM) = plot(ResultsALL(iM).DataTrain.t(ResultsALL(iM).Ntrain_vec), ResultsALL(iM).err(:,1),'-','Color',ccolors(iM,:),'LineWidth',1);
    plot(ResultsALL(iM).DataTrain.t(ResultsALL(iM).Ntrain_vec), ResultsALL(iM).err(:,2),'--','Color',ccolors(iM,:),'LineWidth',1);
end

xlim([0 max(ResultsALL(iM).DataTrain.t)]), ylim([0 12000])
xlabel('Training Time');
ylabel('MSE')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 150])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_MSE','_noleg.eps']);
l1 = legend(ph,'DMDc','SINDYc','NARX');
set(l1,'Location','NorthEast')
l1.Position = [l1.Position(1)-0.01 l1.Position(2)+0.02 l1.Position(3)-0.01 l1.Position(4)];
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_MSE','.eps']);

delete(l1)
xlim([0 30]), ylim([0 2500])
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_MSE','_ZOOM.eps']);

% Prediction over training
figure, hold on, box on
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
for iM = 1:3
    ph(iM) = plot(ResultsALL(iM).DataTrain.t(ResultsALL(iM).Ntrain_vec), ResultsALL(iM).err(:,1),'-','Color',ccolors(iM,:),'LineWidth',1);
end

xlim([0 max(ResultsALL(iM).DataTrain.t)]), ylim([0 12000])
xlabel('Training Time');
ylabel('MSE')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 150])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_MSE_train','_noleg.eps']);
l1 = legend(ph,'DMDc','SINDYc','NARX');
set(l1,'Location','NorthEast')
l1.Position = [l1.Position(1)-0.01 l1.Position(2)+0.02 l1.Position(3)-0.01 l1.Position(4)];
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_MSE_train','.eps']);

delete(l1)
xlim([0 30]), ylim([0 2500])
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_MSE_train','_ZOOM.eps']);

% Prediction over validation
figure, hold on, box on
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];

for iM = 1:3
    data2plot = ResultsALL(iM).err(:,2);
    if any(isnan(data2plot)==1)
        data2plot(find(isnan(ResultsALL(iM).err(:,2)))) = 10^6;
    end
    ph(iM) = plot(ResultsALL(iM).DataTrain.t(ResultsALL(iM).Ntrain_vec), data2plot,'-','Color',ccolors(iM,:),'LineWidth',1);
end

xlim([0 max(ResultsALL(iM).DataTrain.t)]),

% if log scale
if LOG_SCALE
    ylim([10^-1 5*10^4])
    set(gca,'yscale','log')
    set(gca,'ytick',[10.^[-1:2:4]])
end

xlabel('Training Time');
ylabel('MSE')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 150])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_MSE_valid','_noleg.eps']);
l1 = legend(ph,'DMDc','SINDYc','NARX');
set(l1,'Location','NorthEast')
l1.Position = [l1.Position(1)-0.01 l1.Position(2)+0.02 l1.Position(3)-0.01 l1.Position(4)];
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_MSE_valid','.eps']);

delete(l1)
xlim([0 30]), ylim([0 2500])
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_MSE_valid','_ZOOM.eps']);

%% Relative error
clear ph
% Prediction over training/validation
figure, hold on, box on
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
for iM = 1:3
    ph(iM) = plot(ResultsALL(iM).DataTrain.t(ResultsALL(iM).Ntrain_vec), ResultsALL(iM).RelErr(:,1),'-','Color',ccolors(iM,:),'LineWidth',1);
    plot(ResultsALL(iM).DataTrain.t(ResultsALL(iM).Ntrain_vec), ResultsALL(iM).RelErr(:,1),'--','Color',ccolors(iM,:),'LineWidth',1);
end

xlim([0 max(ResultsALL(iM).DataTrain.t)]), ylim([0 6])
xlabel('Training Time');
ylabel('Avg. Rel. Error')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 150])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_RelErr','_noleg.eps']);
l1 = legend(ph,'DMDc','SINDYc','NARX');
set(l1,'Location','NorthEast')
l1.Position = [l1.Position(1)-0.01 l1.Position(2)+0.02 l1.Position(3)-0.01 l1.Position(4)];
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_RelErr','.eps']);

delete(l1)
xlim([0 30]), ylim([0 2])
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_RelErr','_ZOOM.eps']);

% Prediction over training
figure, hold on, box on
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
for iM = 1:3
    ph(iM) = plot(ResultsALL(iM).DataTrain.t(ResultsALL(iM).Ntrain_vec), ResultsALL(iM).RelErr(:,1),'-','Color',ccolors(iM,:),'LineWidth',1);
end

xlim([0 max(ResultsALL(iM).DataTrain.t)]), ylim([0 6])
xlabel('Training Time');
ylabel('Avg. Rel. Error')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 150])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_RelErr_train','_noleg.eps']);
l1 = legend(ph,'DMDc','SINDYc','NARX');
set(l1,'Location','NorthEast')
l1.Position = [l1.Position(1)-0.01 l1.Position(2)+0.02 l1.Position(3)-0.01 l1.Position(4)];
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_RelErr_train','.eps']);

delete(l1)
xlim([0 30]), ylim([0 2])
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_RelErr_train','_ZOOM.eps']);

% Prediction over validation
figure, hold on, box on
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
for iM = 1:3
    data2plot = ResultsALL(iM).RelErr(:,2);
    if any(isnan(ResultsALL(iM).err(:,2))==1)
        data2plot(find(isnan(ResultsALL(iM).err(:,2)))) = 10^6;
    end
    ph(iM) = plot(ResultsALL(iM).DataTrain.t(ResultsALL(iM).Ntrain_vec), data2plot,'-','Color',ccolors(iM,:),'LineWidth',1);
end

xlim([0 max(ResultsALL(iM).DataTrain.t)]), ylim([0 6])

% if log scale
if LOG_SCALE
    ylim([5*10^-3 10^1])
    set(gca,'yscale','log')
    set(gca,'ytick',[10.^[-3:1:1]])
end

xlabel('Training Time');
ylabel('Avg. Rel. Error')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 150])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_RelErr_valid','_noleg.eps']);
l1 = legend(ph,'DMDc','SINDYc','NARX');
set(l1,'Location','NorthEast')
l1.Position = [l1.Position(1)-0.01 l1.Position(2)+0.02 l1.Position(3)-0.01 l1.Position(4)];
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_RelErr_valid','.eps']);

delete(l1)
xlim([0 30]), ylim([0 2])
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_RelErr_valid','_ZOOM.eps']);

%% Training time
clear ph
figure, hold on, box on
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
for iM = 1:3
    ph(iM) = plot(ResultsALL(iM).DataTrain.t(ResultsALL(iM).Ntrain_vec), ResultsALL(iM).Ttraining,'-','Color',ccolors(iM,:),'LineWidth',1);
end

xlim([0 max(ResultsALL(iM).DataTrain.t)]), %ylim([0 6])
xlabel('Training Length');
ylabel('Training Time [s]')

% if log scale
set(gca,'yscale','log')
ylim([10^-3 3*10^0])
set(gca,'ytick',[10^-3 10^-2 10^-1 10^0])

set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 150])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_TrainTime','_noleg.eps']);
l1 = legend(ph,'DMDc','SINDYc','NARX');
set(l1,'Location','NorthWest')
l1.Position = [l1.Position(1)-0.01 l1.Position(2)+0.02 l1.Position(3)-0.01 l1.Position(4)];
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_TrainingLength_TrainTime','.eps']);

author: eurika-kaiser

repo: SINDY-MPC


text:
classdef math_model_opf_dc < mp.math_model_opf
%MP.MATH_MODEL_OPF_DC  MATPOWER mathematical model for DC optimal power flow (OPF) problem.
%   ?
%
%   MP.MATH_MODEL_OPF_DC ... power flow ...
%
%   Properties
%       ? - ?
%
%   Methods
%       ?

%   MATPOWER
%   Copyright (c) 2021-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

    methods
        %% constructor
        function obj = math_model_opf_dc()
            obj@mp.math_model_opf();
            obj.element_classes = { @mp.mme_bus_opf_dc, @mp.mme_gen_opf_dc, ...
                @mp.mme_load_pf_dc, @mp.mme_branch_opf_dc, @mp.mme_shunt_pf_dc };
        end

        function tag = form_tag(obj)
            tag = 'dc';
        end

        function name = form_name(obj)
            name = 'DC';
        end

        function [vx, z, x] = convert_x_m2n(obj, mmx, nm)
            nm_vars = obj.update_nm_vars(mmx, nm);

            %% convert (real) math model x to network model x
            vx = nm_vars.va;
            z  = nm_vars.z;
            if nargout < 2
                vx = [vx; z];
            elseif nargout > 2
                x = [vx; z];
            end
        end

        function obj = add_node_balance_constraints(obj, nm, dm, mpopt)
            [B, K, p] = nm.get_params();

            %% power balance constraints
            C = nm.C;
            Amis = C * [B*C' K*nm.D'];
            bmis = -C * p;
            obj.add_lin_constraint('Pmis', Amis, bmis, bmis, ...
                                [nm.va.order; nm.z.order]);
        end

        function opt = solve_opts(obj, nm, dm, mpopt)
            opt = mpopt2qpopt(mpopt, obj.problem_type());

            switch opt.alg
                case {'MIPS', 'IPOPT'}
                    if mpopt.opf.start < 2      %% initialize interior point
                        opt.x0 = obj.interior_x0(obj, nm, dm);
                    end
                case 'OSQP'
                    opt.x0 = [];        %% disable provided starting point
            end
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function plotLowLevelPerformanceMrf

load('./results/classifierComparisonResult_stage3.mat');
load('./results/classifierComparisonResult_stage3mrf.mat');
load('./results/classifierComparisonResult_stage3mrfg.mat');

figure(1), hold off
plot(pr_pb.r, pr_pb.p, 'b:', 'LineWidth', 4);

hold on
plot(pr_im.r, pr_im.p, 'g--', 'LineWidth', 3);
plot(pr_all.r, pr_all.p, 'r', 'LineWidth', 3);
plot(pr_mrf.r, pr_mrf.p, '--c', 'LineWidth', 4);
plot(pr_mrfg.r, pr_mrfg.p, '--m', 'LineWidth', 4);

legend({'Pb Only', 'Pb+Edge/Region Cues', 'Pb+Edge/Region+3D Cues', 'All Cues + MRF'}, ...
    'Location', 'SouthWest', 'FontSize', 12);
title('Boundary Classification', 'FontSize', 20);
xlabel('Recall', 'FontSize', 18);
ylabel('Precision', 'FontSize', 18);
axis([0 1 0 1])
set(gca, 'FontSize', 14);

author: Cloud-CV

repo: object-proposals


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the RUBBERBAND-EXAMPLE geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Pulsing_Heart()

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx =  64;        % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Ny =  64;        % # of Eulerian Grid Pts. in y-Direction (MUST BE EVEN!!!)
Lx = 1.0;        % Length of Eulerian Grid in x-Direction
Ly = 1.0;        % Length of Eulerian Grid in y-Direction

% Immersed Structure Geometric / Dynamic Parameters %
N = 2.5*Nx;              % Number of Lagrangian Pts. (2x resolution of Eulerian grid)
struct_name = 'heart'; % Name for .vertex, .spring, etc files.

% Call function to construct geometry
frac1 = 0.0135; %0.022
[xA,yA] = give_Me_Immsersed_Boundary_Geometry_1(Lx,Nx,frac1);
Nb1 = length(xA);

% Call function to construct geometry
frac2 = 0.0085; %0.015
[xB,yB] = give_Me_Immsersed_Boundary_Geometry_2(Lx,Nx,frac2,Nb1);

% Plot Geometry to test BEFORE taking out pts.
figure(1)
plot(xA,yA,'r-'); hold on;
plot(xA,yA,'*'); hold on;
plot(xB,yB,'m-'); hold on;
plot(xB,yB,'g*'); hold on;
xlabel('x'); ylabel('y');
axis([0 Lx 0 Ly]);

% Cuts hole into heart geometry
n1 = 145; n2 = 180;
[xA,yA,xB,yB] = please_Take_Out_Points(xA,yA,xB,yB,n1,n2);

% Plot Geometry to test AFTER taking out pts.
figure(2)
plot(xA,yA,'r-'); hold on;
plot(xA,yA,'*'); hold on;
plot(xB,yB,'m-'); hold on;
plot(xB,yB,'g*'); hold on;
xlabel('x'); ylabel('y');
axis([0 Lx 0 Ly]);

% Print STATES A,B, and C:
print_States(xA,yA,'State_A');
print_States(xB,yB,'State_B');

% Print files to .txt files
%please_Print_Vertices_To_File(xA,yA,xB,yB)

% Prints .vertex file!
print_Lagrangian_Vertices(xA,yA,struct_name);

% Prints .spring file!
k_Spring = 1e3; ds_Rest = 0;
print_Lagrangian_Springs(xA,k_Spring,ds_Rest,struct_name,n1)

% Prints .target file!
k_Target = 2e10;
print_Lagrangian_Target_Pts(xA,k_Target,struct_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called State_<j>.pts
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_States(xA,yA,struct_name)

    N = length(xA);

    vertex_fid = fopen([struct_name '.pts'], 'w');

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xA(s);
        Y_v = yA(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called heart.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag);

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints Vertex points to a file called heart.target
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = length(xLag);

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called rubberband.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,k_Spring,ds_Rest,struct_name,n1)

    N = length(xLag);

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %SPRINGS BETWEEN VERTICES
    for s = 1:N
            if s < n1
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            elseif ( (s>=n1+1) && (s<N) )
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            elseif s == N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, 1,   k_Spring, ds_Rest);
            end
    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry for PHASE 1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry_1(Lx,Nx,frac)

% The immsersed structure is a heart %
tVec = 0:Lx/(2*Nx):2*pi;
for i=1:length(tVec)
    t = tVec(i);
    xLag(i)	=	16*sin(t)^3;
    yLag(i)	=	13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t);
end

xLag = frac*xLag + 0.5;
yLag = frac*yLag + 0.5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry for PHASE 1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry_2(Lx,Nx,frac,Nb)

% The immsersed structure is a heart %
ds = 2*pi/(Nb-1);
tVec = 0:ds:2*pi;

for i=1:length(tVec)
    t = tVec(i);
    xLag(i)	=	16*sin(t)^3;
    yLag(i)	=	13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t);
end

xLag = frac*xLag + 0.5;
yLag = frac*yLag + 0.5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry for PHASE 2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry_3(Lx,Nx,frac,Nb)

% The immsersed structure is a heart %
ds = 2*pi/(Nb-1);
tVec = 0:ds:2*pi;
for i=1:length(tVec)
    t = tVec(i);
    r(i) = 1 - sin(t);
    xLag(i)	=	r(i)*cos(t);
    yLag(i)	=	r(i)*sin(t);
end

xLag = frac*xLag;
yLag = frac*yLag;

minY = min(yLag);
yLag = yLag - 3*minY/8;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: take out Lagrangian Pts.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x1,y1,x2,y2] = please_Take_Out_Points(x1,y1,x2,y2,n1,n2)

x1 = [x1(1:n1) x1(n2:end)];
y1 = [y1(1:n1) y1(n2:end)];

x2 = [x2(1:n1) x2(n2:end)];
y2 = [y2(1:n1) y2(n2:end)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints all Vertices to File
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function please_Print_Vertices_To_File(X1,Y1,X2,Y2)

fileID = fopen('All_Positions.txt','w');
for j=1:length(X1)
    fprintf(fileID,'%1.16e %1.16e %1.16e %1.16e\n', X1(j),Y1(j),X2(j),Y2(j) );
end
fclose(fileID);

author: nickabattista


text:
% Display Model headers (MR protocol, model detail, fitting parameters, options)
% Example:
%   qMRinfo('mwf'); % display help for Myelin Water Fraction
%
% This function is to provide Octave compatibility for help text.
% Octave accepts a different convention for help documentation.
% Until we find a through solution for Octave help, we will skip this for
% Octave.

% This is the reason why this function will be used instead of 'help'

function varargout = qMRinfo(varargin)
if nargin==0, varargin{1}='Display Model headers (MR protocol, model detail, fitting parameters, options)'; end
if moxunit_util_platform_is_octave || isdeployed
    hh = load('iqmr_gethelp.mat');
    if ischar(varargin{1})
        varargout{1} = hh.(varargin{1});
    else
        varargout{1} = hh.(class(varargin{1}));
    end
else
    [varargout{1:nargout}] = help(varargin{:});
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [H] = traditional(f);

% TRADITIONAL creates the homogenous spatial transformation matrix
% for a 9 parameter traditional "Talairach-model" transformation
%
% Use as
%   [H] = traditional(f)
%
% The transformation vector f should contain the
%   x-shift
%   y-shift
%   z-shift
% followed by the
%   pitch (rotation around x-axis)
%   roll  (rotation around y-axis)
%   yaw   (rotation around z-axis)
% followed by the
%   x-rescaling factor
%   y-rescaling factor
%   z-rescaling factor
%
% The order in which the transformations are done is exactly opposite as
% the list above, i.e. first z-rescale, ... and finally x-shift.

% Copyright (C) 2000-2005, Robert Oostenveld
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% $Log: traditional.m,v $
% Revision 1.1  2009/01/30 04:02:12  arno
% *** empty log message ***
%
% Revision 1.5  2005/08/15 08:15:33  roboos
% reimplemented the rotate function, which contained an error (the error is in the AIR technical reference)
% changed all functions to be dependent on the rotate, translate and scale function
% all functions now behave consistenly, which also means that they are not compleetly backward compatible w.r.t. the order of the rotations
%
% Revision 1.4  2005/08/11 07:57:14  roboos
% fixed bug in y-rotation
%
% Revision 1.3  2005/04/21 08:28:45  roboos
% fixed bug in rotation matrix (thanks to Arno)
%
% Revision 1.2  2004/05/19 09:57:07  roberto
% added GPL copyright statement, added CVS log item
%

% compute the homogenous transformation matrix for the translation
T = translate(f([1 2 3]));

% compute the homogenous transformation matrix for the rotation
R = rotate(f([4 5 6]));

% compute the homogenous transformation matrix for the scaling
S = scale(f([7 8 9]));

% compute the homogenous transformation matrix for the combination
H = T*R*S;

author: PatternRecognition

repo: OpenBMI


text:
function mr = mrReadOldParamMap(pth);
%
%  mr = mrReadOldParamMap(pth);
%
% Reads in a mrVista 1.0 Parameter map file
% as a mrVista 2.0 mr object. Tries to extract as
% much information as possible from the mrSESSION
% file, returning the rest in the mr.hdr field.
%
% For now, this only works on Original, Inplane TSeries (or
% TSeries files whose params aren't significantly
% changed from the Original tSeries). The tSeries files
% must be in a directory named 'Scan[#]'.
%
% ras, 07/08/05.
% ras, 03/01/07: reads in data units for newer param maps.

[p f ext] = fileparts(pth);

% init mr struct
mr = mrCreateEmpty;
mr.path = pth;
mr.format = '1.0map';
mr.dimUnits = {'mm' 'mm' 'mm' 'scan'};

% read in header info from mrSESSION:
% if we can't do this, we can't read the data
mrSessPath = fullfile(p,'..','..','mrSESSION.mat');
if exist(mrSessPath,'file')
    load(mrSessPath,'mrSESSION','dataTYPES');
else
    % dialog to get file
    error('No mrSESSION file found!')
end

% load the param map file
M = load(pth);
if ~isfield(M, 'map') | ~isfield(M, 'mapName')
	error( sprintf('%s needs ''map'' and ''mapName'' variables.', pth) );
end

mr.name = M.mapName;
map = M.map;  M.map =[]; % quick hack for back-compatibility

% check for an empty map
if isempty(cellfind(map))
    myErrorDlg('Map File is empty!');
end

% since most contrast maps are named w/ the format
% '[active]V[control]', and no other common param
% map has a V in it, make an educated guess about
% units:
if isfield(M, 'units'),				mr.dataUnits = M.units;
elseif findstr(M.mapName,'V'),		mr.dataUnits = '-log_{10}(p)';
else,								mr.dataUnits = M.mapName;
end

% assign each volume that is nonempty to a new time point:
for i = 1:length(map), ind(i) = ~isempty(map{i}); end
ind = find(ind);
for t = 1:length(ind)
    mr.data(:,:,:,t) = map{ind(t)};
	mr.info.scans(t) = ind(t);
end

% get other header info, using first assigned scan as index
% (first, we see if this map reflects data from the Original data type:
%  if so, we can use the scan index as an index into the functional
%  headers. Otherwise, we guess the first scan header is appropriate.)
[p2 f2] = fileparts(p);
if isequal(f2, 'Original')
	nScan = ind(1);
else
	nScan = 1; % guess the 1st scan
end
mr.hdr = mrSESSION.functionals(nScan);
mr.hdr.comments = sprintf('Computed for scans: %s',num2str(ind));
if (~isfield(mr.hdr,'voxelSize')) && (isfield(mr.hdr,'effectiveResolution'))
    mr.hdr.voxelSize=mr.hdr.effectiveResolution;
end

mr.voxelSize = [mr.hdr.voxelSize];
mr.dims = [mr.hdr.cropSize length(mr.hdr.slices)];
mr.extent = mr.dims .* mr.voxelSize;

if isfield(mrSESSION,'alignment')
    mr.spaces(1).name = 'I|P|R';
    mr.spaces(1).xform = mrSESSION.alignment;
    mr.spaces(1).dirLabels =  {'S <--> I' 'A <--> P' 'L <--> R'};
    mr.spaces(1).sliceLabels =  {'Axial' 'Coronal' 'Sagittal'};
    mr.spaces(1).units = 'mm';
    mr.spaces(1).coords = [];
    mr.spaces(1).indices = [];

    % account for differences between inplane / functional resolutions
    % for the alignment:
    % the mrVista 1.0 alignments mapped from inplane coords to volume
    % coords, rather than the map / tSeries:
    szRatio = mr.voxelSize ./ mrSESSION.inplanes.voxelSize;
    [trans rot scale] = affineDecompose(mr.spaces(1).xform);
    mr.spaces(1).xform = affineBuild(trans,rot,scale.*szRatio);

    mr.spaces(2).name = 'R|A|S';
    mr.spaces(2).xform = ipr2ras(mr.spaces(1).xform,mr.dims);
    mr.spaces(2).dirLabels =  {'L <--> R' 'P <--> A' 'I <--> S'};
    mr.spaces(2).sliceLabels =  {'Sagittal' 'Coronal' 'Axial'};
    mr.spaces(2).units = 'mm';
    mr.spaces(2).coords = [];
    mr.spaces(2).indices = [];
end

% pre-pend standard space definitions
mr.spaces = [mrStandardSpaces(mr) mr.spaces];

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function f = p31_f ( n, x )

%*****************************************************************************80
%
%% P31_F evaluates the objective function for problem 31.
%
%  Discussion:
%
%    The minimal function value is -10.15320.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Zbigniew Michalewicz,
%    Genetic Algorithms + Data Structures = Evolution Programs,
%    Third Edition,
%    Springer Verlag, 1996,
%    ISBN: 3-540-60676-9,
%    LC: QA76.618.M53.
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  m = 5;
  a = [ 4.0, 4.0, 4.0, 4.0; ...
        1.0, 1.0, 1.0, 1.0; ...
        8.0, 8.0, 8.0, 8.0; ...
        6.0, 6.0, 6.0, 6.0; ...
        3.0, 7.0, 3.0, 7.0 ]';

  c = [ 0.1, 0.2, 0.2, 0.4, 0.6 ]';

  f = 0.0;
  for j = 1 : m
    f = f - 1.0 / ( c(j) + sum ( ( x(1:n) - a(1:n,j) ).^2 ) );
  end

  return
end


text:
function res = innerprod(X,Y)
%INNERPROD Efficient inner product with a ktensor.
%
%   R = INNERPROD(X,Y) efficiently computes the inner product between
%   two tensors X and Y.  If Y is a ktensor, the inner product is
%   computed using inner products of the factor matrices, X{i}'*Y{i}.
%   Otherwise, the inner product is computed using ttv with all of
%   the columns of X's factor matrices, X{i}.
%
%   See also KTENSOR, KTENSOR/TTV
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if ~isequal(size(X),size(Y))
    error('X and Y must be the same size.');
end

% X is a ktensor
switch class(Y)

  case {'ktensor'}
    M = X.lambda * Y.lambda';
    for n = 1:ndims(X)
        M = M .* (X.u{n}' * Y.u{n});
    end
    res = sum(M(:));

  case {'tensor','sptensor','ttensor'}
    R = length(X.lambda);
    vecs = cell(1,ndims(X));
    res = 0;
    for r = 1:R
      for n = 1:ndims(X)
        vecs{n} = X.u{n}(:,r);
      end
      res = res + X.lambda(r) * ttv(Y,vecs);
    end

  otherwise
    disp(['Inner product not available for class ' class(Y)]);
end

return;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function Z = projectData(X, U, K)
%PROJECTDATA Computes the reduced data representation when projecting only
%on to the top k eigenvectors
%   Z = projectData(X, U, K) computes the projection of
%   the normalized inputs X into the reduced dimensional space spanned by
%   the first K columns of U. It returns the projected examples in Z.
%

% You need to return the following variables correctly.
Z = zeros(size(X, 1), K);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the projection of the data using only the top K
%               eigenvectors in U (first K columns).
%               For the i-th example X(i,:), the projection on to the k-th
%               eigenvector is given as follows:
%                    x = X(i, :)';
%                    projection_k = x' * U(:, k);
%

U_reduced = U(:, 1:K);
Z = X * U_reduced;

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

end

repo: MachineLearning


text:
%
% Nonuniform FFT in R^2.
%   nufft2d1 - Type 1 - Evaluate linearly spaced Fourier modes (Type 1 nuFFT).
%   nufft2d2 - Type 2 - Evaluate function values (inverse of Type 1 nuFFT).
%   nufft2d3 - Type 3 - Evaluate arbitrary located Fourier modes (Type 3 nuFFT).
%
% Direct (slow) evaluation of nonuniform FFT in R^2.
%   dirft2d1 - Type 1 - Evaluate linearly spaced Fourier modes (Type 1 nuFFT).
%   dirft2d2 - Type 2 - Evaluate function values (inverse of Type 1 nuFFT).
%   dirft2d3 - Type 3 - Evaluate arbitrary located Fourier modes (Type 3 nuFFT).
%
% Test for nuFFT in R^2.
%   test_2d1 - Test Type 1 nuFFT.
%   test_2d2 - Test Type 2 nuFFT.
%   test_2d3 - Test Type 3 nuFFT.
%


text:
function n = r8vec_read_size ( input_file )

%*****************************************************************************80
%
%% R8VEC_READ_SIZE reads the size of an R8VEC from a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE, the name of the file to be read.
%
%    Output, integer N, the size of the R8VEC.
%
  n = 0

  input_unit = fopen ( input_file );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_READ_SIZE - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'R8VEC_READ_SIZE - Error!' );
    return;
  end

  while ( 1 )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    n = n + 1;

  end

  fclose ( input_unit );

  return
end


text:
%BUILDPYRAMID  Constructs the Gaussian pyramid for an image
%
%     dst = cv.buildPyramid(src)
%     dst = cv.buildPyramid(src, 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ Source image. Check cv.pyrDown for the list of supported types.
%
% ## Output
% * __dst__ Destination vector of `Maxlevel+1` images of the same type as
%   `src`. A cell array of images. `dst{1}` will be the same as `src`.
%   `dst{2}` is the next pyramid layer, a smoothed and down-sized `src`, and
%   so on.
%
% ## Options
% * __MaxLevel__ 0-based index of the last (the smallest) pyramid layer. It
%   must be non-negative. default 5
% * __BorderType__ Pixel extrapolation method, ('Constant' isn't supported).
%   See cv.copyMakeBorder for details. Default 'Default'
%
% The function constructs a vector of images and builds the Gaussian pyramid
% by recursively applying cv.pyrDown to the previously built pyramid layers,
% starting from `dst{1}==src`.
%
% See also: cv.pyrDown, cv.pyrUp, vision.Pyramid, impyramid
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function g = diagKernGradient(kern, x, varargin)

% DIAGKERNGRADIENT Gradient of DIAG kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% diagonal noise covariance function
% kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO diagKernParamInit, kernGradient, diagKernDiagGradient, kernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2011

% KERN

  if nargin < 4
    trans = str2func([kern.trans, 'Transform']);
    g(1, 1) = sum(diag(varargin{end}).*trans(x, 'atox'));
  else
    g = 0;
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function clf = defineClf(obj, params, symbolic_state)
    v = symbolic_state(2);
    vd = params.vd; % desired velocity.

    clf = (v - vd)^2;
end

author: HybridRobotics

repo: CBF-CLF-Helper


text:
function [g,m] = do_pedro_transfer(models)
error('deprecated');

m = models;
for i = 1:length(models)
  m{i}.models_name = 'pedroxfer';
end

addpath(genpath('/nfs/hn22/tmalisie/ddip/util/voc-release3.1/'))
load bus_final.mat
fg = get_pascal_bg('test','bus');

allbg = cat(1,get_pascal_bg('trainval'),get_pascal_bg('test'));
curids = cell(length(allbg),1);
for i = 1:length(allbg)
  [tmp,curids{i},tmp] = fileparts(allbg{i});
end

%fg = fg(1:5);
g = cell(length(fg),1);

xs = cellfun2(@(x)x.model.x,models);
xs = cat(2,xs{:});

xs2 = xs;
for i = 1:size(xs2,2)
  xs2(:,i) = xs2(:,i) - mean(xs2(:,i));
end

for i = 1:length(fg)
  I = convert_to_I(fg{i});
  boxes = detect(I,model,model.thresh);
  boxes = nms(boxes,.5);

  %if size(boxes,1)==0
  %  continue
  %end

  %newboxes are in my format, boxes is in pedros format
  newboxes = zeros(size(boxes,1),8);
  if size(boxes,1)>0
    [aa,bb] = sort(boxes(:,end),'descend');
    boxes = boxes(bb,:);
    newboxes(:,1:4) = boxes(:,1:4);
    newboxes(:,end) = boxes(:,end);
  end

  for j = 1:size(boxes,1)
    curm = initialize_model_dt(I,boxes(j,:),...
                               models{1}.model.params.sbin,...
                               models{1}.model.hg_size);

    [tmp,ind] = max(xs2'*curm.x);
    newboxes(j,6) = ind;
  end

  g{i}.bboxes = newboxes;
  g{i}.imbb = [1 1 size(I,2) size(I,1)];
  [tmp,g{i}.curid,tmp] = fileparts(fg{i});

  g{i}.index = find(ismember(g{i}.curid,curids));
  g{i}.extras = [];
  %% no extras, no coarse boxes

  figure(1)
  clf
  imagesc(I)
  if size(newboxes,1)>0
    plot_bbox(newboxes(1,:))
    title('Top Detection');
  end
  % figure(2)
  % [aa,bb] = sort(xs2'*curm.x,'descend');
  % for q = 1:9
  %   subplot(3,3,q)
  %   imagesc(get_exemplar_icon(models,bb(q)));
  % end
  drawnow

end

author: quantombone

repo: exemplarsvm


text:
%filename = 'CantileverSquareSYmmetricMesh';
filename = 'CantileverSquareMedium';
%'CantileverSquareSYmmetricMesh';
%filename = 'CantileverSquare';
%filename = 'CantileverSquareSmall';
%filename = 'CantileverSquareNew';
%'CantileverSquareNewFine';
%filename = 'Cantilever_quad_coarse';
%filename = 'LshapeTriFine';
%filename = 'LshapeTri';
%filename = 'LshapeTriSmall';
%filename = 'Lshape';
%filename = 'LshapeFine';
%filename = 'ArchTriFine';
%'Arch_quad_coarse';
%'Bridge_quad_coarse';
%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';
%filename = 'Bridge';

ptype = 'MACRO';
initial_case = 'given';
m1 = 0.0101;
m2 = 0.0101;
%cost = {'compliance'};
cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%weights = [0.55,0.45];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'PDE';
%filterType = 'P1';
constraint_case = 'EQUALITY';

Vfrac_initial = 0.3;
optimality_initial = 1e-4;
constr_initial = 1e-8;

Vfrac_final = 0.3;
optimality_final = 1e-4;
constr_final = 1e-8;

stressNormExponent_initial = 2;
stressNormExponent_final = 32;
%
optimizer = 'DualNestedInPrimal';
%optimizer = 'AlternatingPrimalDual';
optimizerUnconstrained = 'PROJECTED GRADIENT';

%optimizer = 'MMA';
%optimizer = 'IPOPT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
%vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
vademecumFileName = 'SuperEllipseQOptAnalytic';

% designVariable = 'LevelSet';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';
% initial_case = 'full';% optimizerUnconstrained = 'SLERP';

% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';
% initial_case = 'full';
% rho0 = 0.3;

line_search_initiator = 'INCREASING LAST STEP';
incrementFactor = 1.95;
%

%kfrac = 2;
nsteps = 17;

plotting = true;
printing = true;
monitoring = true;
monitoring_interval = 3;
maxiter = 8000;

%
% % %
isDirichletPartX = @(x) x > -1e-12 & x < 0.1;
isDirichletPart1 = @(y) y > 0.20 & y < 0.30;
isDirichletPart2 = @(y) y > 0.70 & y < 0.80;
isDirichletPartY = @(y) isDirichletPart1(y) | isDirichletPart2(y);
isDirichletPart = @(x,y) isDirichletPartX(x) & isDirichletPartY(y);
isNeumannPartX = @(x) x > (1-0.05) & x < (1+1e-12);
isNeumannPartY = @(y) y > 0.40 & y < 0.60;

isNeumannPart = @(x,y) isNeumannPartX(x) & isNeumannPartY(y);
iNotOptimizable = @(coord) isDirichletPart(coord(:,1),coord(:,2)) & ~isNeumannPart(coord(:,1),coord(:,2));

costDomainNotOptimizable       = {iNotOptimizable};
constraintDomainNotOptimizable = {[]};

isDesignVariableFixed.nodes  = iNotOptimizable;
isDesignVariableFixed.values = @(x) [m1*ones(size(x,1),1);m2*ones(size(x,1),1)];

author: SwanLab


text:
function [in2] = km22in2(km2)
% Convert area from square kilometers to square inches.
% Chad A. Greene 2012
in2 = km2*1550003100.006;

author: Sable

repo: mcbench-benchmarks