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

% Discretization of a 2d Variance Gamma process

% Monte Carlo parameters
NBatches = 1;           % Number of batches
NSim = 10000;%100000;           % Number of paths per batch
Nt = 250;                 % Number of time steps until T!

% Asset parameters
S0 = [1.0 1.0];             % spot price of stock index
r =  0.042;              % risk free rate
d =  [0.00 0.00];              % dividend yield
T = 0.5;                % Time horizon (in years)
rho = 0.5;
lnS1 = zeros(NSim, Nt+1); % log spot prices asset 1
lnS2 = zeros(NSim, Nt+1); % log spot prices asset 2

% Model parameters (theta, sigma, nu representation)
theta = [-0.6094 -0.8301];      % parameter of VG
sigma = [0.0325 0.9406];        % parameter of VG
nu = 0.2570;                    % parameter of VG
omegaT = -1/nu * [log(1-theta(1)*nu - nu*sigma(1)^2/2) log(1-theta(2)*nu - nu*sigma(2)^2/2)];
drift = [r-d(1) r-d(2)];       % parameter of VG

%Corr = [1 .6495; .6495 1];
Corr = [1 rho; rho 1];
%Corr = [1 .9; .9 1];
%Corr = [1 -.9; -.9 1];
R = chol(Corr);

% Option parameters
value = ones(NBatches,1);   % Stores the option value per batch

% precomputed constants
deltaT = T / Nt;                            % delta for time discretization
lnS1(:,1) = log(S0(1));                         % Set the starting spot price
lnS2(:,1) = log(S0(2));

%oNt = ones(Nt,1);                       % used during simulation
oNs = ones(NSim,1);
% Start Monte Carlo here
tic;
for number = 1 : NBatches
    % Time discretization

    for m=1:Nt
        %G = nu*gamrnd(deltaT/nu, oNt);      % Gamma Subordinator
        G = nu * gamrnd(deltaT/nu,oNs);
        W = randn(NSim,2);
        W = W*R;
        lnS1(:,m+1) = lnS1(:,m) + (drift(1)-omegaT(1)) * deltaT ...
                       + theta(1) * G + sqrt(G) * sigma(1) .* W(:,1);
        lnS2(:,m+1) = lnS2(:,m) + (drift(2)-omegaT(2))* deltaT ...
                      + theta(2) * G + sqrt(G) * sigma(2) .* W(:,2);
    end

    S1 = exp(lnS1);         % Simulated prices for asset 1
    S2 = exp(lnS2);         % Simulated prices for asset 2

    if(NSim ==1)
        figure1 = figure;
        axes1 = axes('Parent',figure1);
        hold on;
        plot1 = plot(S1,'Parent',axes1,'MarkerSize',4,'Marker','o','Color',[0 0 1],'DisplayName','Asset 1');
        plot1 = plot(S2, 'Marker','.','Color',[1 0 0],'DisplayName','Asset 2');
        %plot(S1,'b');
        %plot(S2,'r');

        % Create xlabel
        xlabel('step');

        % Create ylabel
        ylabel('S(t)');

        % Create title
        title({'2d Variance Gamma Process','- zero correlation -'},...
        'FontWeight','bold','FontSize',12,'FontName','Arial');

        % Create legend
        legend1 = legend(axes1,'show');
        set(legend1,'Position',[0.8206 0.831 0.08047 0.06204]);

    end

    % Option Pricing
    value(number) = exp(-r*T)* WorstOfCall(S1,S2);
    %value(number) = exp(-r*T)* BestOfCall(S1,S2);
    %value(number) = exp(-r*T)* Spread(S1,S2);

end

mean(value)      % Output of Option price
%Elapsed_Time = toc              % Time spend on MC simulation


text:
function [varargout]=grid2patch(varargin)

% function [F,V,C,CV]=grid2patch(X,Y,Z,C,perdiocOpt);
% ------------------------------------------------------------------------
%
%
% 2019/12/16: Created as alternative to surf2patch (missing in Octave)
% ------------------------------------------------------------------------

%% parse input

switch nargin
    case 2
        X=varargin{1};
        Y=varargin{2};
        Z=zeros(size(X));
        C=[];
        perdiocOpt=false(1,2);
    case 3
        X=varargin{1};
        Y=varargin{2};
        Z=varargin{3};
        C=[];
        perdiocOpt=false(1,2);
    case 4
        X=varargin{1};
        Y=varargin{2};
        Z=varargin{3};
        C=varargin{4};
        perdiocOpt=false(1,2);
    case 5
        X=varargin{1};
        Y=varargin{2};
        Z=varargin{3};
        C=varargin{4};
        perdiocOpt=varargin{5};
end

if isempty(C)
    C=Z(:);
end
siz=size(X);

%%

V=[X(:) Y(:) Z(:)]; %Vertex set

%Create row of faces
f=[1 1+siz(1) 2+siz(1) 2]; %First element
q=(0:1:siz(1)-2)';
qq=q(:,ones(4,1));
ff=f(ones(siz(1)-1,1),:)+qq; %Row of faces by copying first

if perdiocOpt(1)==1
    ff(end+1,:)=[siz(1) siz(1)+siz(1) 1+siz(1) 1];
end

%Create grid of faces by copying column
q=0:siz(1):(siz(1)*(siz(2)-2));
if perdiocOpt(1)==1
    Q=q(ones(siz(1),1),:);
else
    Q=q(ones(siz(1)-1,1),:);
end
Q=Q(:);
Q=Q(:,ones(1,4));
F=repmat(ff,[siz(2)-1 1])+Q; %Grid of faces

if perdiocOpt(1)==1 && perdiocOpt(2)==1
    t=siz(1)*(siz(2)-1);
    f=[2 2+t 1+t 1]; %First element
    q=(0:1:siz(1)-2)';
    qq=q(:,ones(4,1));
    ff=f(ones(siz(1)-1,1),:)+qq; %Row of faces by copying first
    ff(end+1,:)=[1 1+t siz(1)+t siz(1)];
    F=[F;ff];
elseif perdiocOpt(2)==1
    t=(siz(1)*(siz(2)-1));
    f=[2 2+t 1+t 1]; %First element
    q=(0:1:siz(1)-2)';
    qq=q(:,ones(4,1));
    ff=f(ones(siz(1)-1,1),:)+qq; %Row of faces by copying first
    F=[F;ff];
end

%% Collect output

varargout{1}=F;
varargout{2}=V;
if nargout==3
    varargout{3}=vertexToFaceMeasure(F,C(:));
end
if nargout==4
    A=C;
    varargout{4}=A(:);
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%% This file is the main file of using the SINDy-PI method to
% infer the Yeast Glycolysis Model. We will figure out what is the minimum
% data length needed for SINDy-PI to accurately discover the six-th state
% of the Yeast Glycolysis Model. This file will be used to swipe through
% different sparse parameter.
%
% Date: 2019/06/12
% Coded By: K

%% Close all, clear all, clc
close all;clear all; clc;
set(0,'defaulttextInterpreter','latex')
addpath('Functions')
addpath('Datas')
%% Define some parameters
% Define whehter you have control, if you have it, please define it
n_control=0;u=0;

% Run the ODE files and gather the simulation data.
% (We use the same data for both the iSINDy method and SINDy-PI method for better comparision)
% (Baseline data for the data length comparison of state 1,2,3,4,5,6,7)
% The data is noise clean. It is generated using 900 different initial
% conditions with 51 points of each initial condition.
load('TrainingData.mat')

% Choose whether you want to display actual ODE or not
disp_actual_ode=1;

% If the ODEs you want to display is the actual underlyting dynamics of the
% system, please set actual as 1
actual=1;

% Define how many states we have in our example
n_state=7;

% Print the actual ODE we try to discover
digits(4)
Print_ODEs(@(t,y)YeastGlycolysis_ODE(t,y),n_state,n_control,disp_actual_ode,actual);

% Create symbolic states
dz=sym('dz',[n_state,1]);

% Now we first create the parameters of the function right hand side
Highest_Poly_Order_Guess=0;
Highest_Trig_Order_Guess=0;
Highest_U_Order_Guess=0;

% Then create the right hand side library parameters
Highest_Trig_Order=0;
Highest_U_Order=0;
Highest_dPoly_Order=1;

% Set the parameter normLib=1 to normalize the librry
normLib=1;

% Set how amny iterations you want for the sparse regression
N_iter=10;

% Set whether you want to display the ODE or not
disp=0;

% Set the library size
Highest_Poly_Order_Lib=[6;6;3;3;3;6;3];

% Determine how many percent of data you want.
percent_start=0.03;d_percent=0.005;percent_end=0.05;

% Determine the lambda
lam_start=1;lam_end=20;d_lambda=0.02;

% Define a cell matrix to store the result
Xi=cell(round(percent_end-percent_start)/d_percent+1,2,lam_end-lam_start+1);

% Determine which states you want to discover
Which_State=6;

% Determine the left hand side guess number
if Which_State==1 || Which_State==2 || Which_State==6
    LHS_Num=2;
else
    LHS_Num=1;
end

% Get the poly-order
Highest_Poly_Order=Highest_Poly_Order_Lib(Which_State);

% Determine which LHS you want to use, the first one or the second one
LHS_Pin=2;

% Create the new directory to save the result
FolderName=strcat('Result_DL_SINDy_Data_Length_Compare_State_',num2str(Which_State),'_LHS_Guess_',num2str(LHS_Pin),'_V2');
[fld_status, fld_msg, fld_msgID]=mkdir(FolderName);

tic
% Create the library data using all the data points
[SINDy_Data_Full,SINDy_Struct]=SINDyLib(xt,dxt(:,Which_State),Which_State,u,Highest_Poly_Order,Highest_Trig_Order,Highest_U_Order,Highest_dPoly_Order);
toc

% Create left hand side guess
[LHS_Data_Full,LHS_Sym]=GuessLib(xt,dxt(:,Which_State),Which_State,u,Highest_Poly_Order_Guess,Highest_Trig_Order_Guess,Highest_U_Order_Guess);

tic
% Create the right hand side, exclude the guess from SINDy library
[RHS_Data_Full,RHS_Struct]=ExcludeGuess(SINDy_Data_Full,SINDy_Struct,LHS_Sym{LHS_Pin});
toc

parlooooop=parpool(2)
%% Start!
for Total_Run=1:20
    fprintf('\n \n Get the result for the %i time...\n',Total_Run)

    % Set a pin to count which iteration it is now
    pinpin=0;

    % Define a cell matrix to store the value of discovery result
    Xi=cell(round(percent_end-percent_start)/d_percent+1,lam_end-lam_start+1);
    ODE_Guess=cell(round(percent_end-percent_start)/d_percent+1,lam_end-lam_start+1);
    ODE=cell(round(percent_end-percent_start)/d_percent+1,lam_end-lam_start+1);

    for percent=percent_start:d_percent:percent_end

        fprintf('\n Testing the percentage as %i ...\n',percent*100)

        pinpin=pinpin+1;

        % Define the new data length
        new_length=round(percent*length(xt));

        % Shuffel the original data
        Sequence=randperm(size(SINDy_Data_Full,1));
        Sequence_Trimed=Sequence(1:new_length);
        RHS_Data=RHS_Data_Full(Sequence_Trimed,:);
        LHS_Data=LHS_Data_Full(Sequence_Trimed,:);

        fprintf('\n \t Calculating the %i expression...\n',Which_State)

        % Print the left hand side that we are testing
        fprintf('\t Testing the left hand side as %s:\n',char(LHS_Sym{LHS_Pin}))

        % Set up dummy variables for parfor
        LHS_Data_Dum=LHS_Data(:,LHS_Pin);
        LHS_Sym_Dum=LHS_Sym{LHS_Pin};
        dz_dum=dz(Which_State);

        % Start parfor and sweep the different values of lambda
        parfor pin=lam_start:lam_end
            % Perform the sparse regression problem
            tic
            [Xi{pinpin,pin},ODE{pinpin,pin}]=sparsifyDynamics(RHS_Data,LHS_Data_Dum,LHS_Sym_Dum,d_lambda*pin,N_iter,RHS_Struct,disp,normLib);
            toc
            % Perform sybolic calculation and solve for dX
            digits(6)
            ODE_Guess{pinpin,pin}=vpa(solve(LHS_Sym_Dum==ODE{pinpin,pin},dz_dum));
        end
    end

    % Save the calculation result of current iteration
    cc=clock;
    ResultName=strcat(FolderName,'/DL_SINDY_Data_Length_Result_',num2str(Total_Run),'__','LHS',num2str(LHS_Pin),'_',num2str(cc(3)),'_',num2str(cc(4)),'_',num2str(cc(5)),'_',num2str(round(cc(6))),'_P2.mat');
    save(ResultName,'Xi','percent_start','d_percent','percent_end','lam_start','lam_end','d_lambda')

end

author: dynamicslab

repo: SINDy-PI


text:
function [parsed_update,annotedText] = writeCD(parsed,fname_out)
% Write the parsed CD model structure to a CD-compatible XML file.
%
% USAGE:
%
%    [parsed_update, annotedText] = writeCD(fname, parsed, fname_out)
%
% INPUTS:
%    fname:            The original XML file
%    parsed:           A parsed model structure generated by `parseCD` function.
%    fanme_out:        The name of the output XML file.
%
% OUTPUTS:
%    parsed_update:    An updated version of the parsed
%    annotedText:      A matlab variable storing the all the XML lines
%
% EXAMPLE:
%
%    aaa = writeXML('PD_140620_1.xml', parsePD,'text.html')
%
% .. Author: - Longfei Mao Oct/2014

if nargin<2 || isempty(fname_out)

    [fname_out, fpath]=uiputfile('*.xml','CellDesigner SBML Source File');
    if(fname_out==0)
        return;
    end
    f_out=fopen([fpath,fname_out],'w');
else
    f_out=fopen(fname_out,'w');
end

if isfield(parsed.r_info,'XMLtext');
    text=parsed.r_info.XMLtext;
else
    errordlg('The XMLtext doesn''t exit in the parsed model structure');
end

%%%%%%%%%%%%%%%%%%%%%%%%  modify

ref=readCD(parsed); %% read the the list of the reactions and associated information.
% ref=parsed.r_info;

r_c=size(ref.number)

%while ~feof(f_id);

%numOfLine=0;
%rem=fgets(f_id);numOfLine=numOfLine+1;

items={'width','color'};

for r=1:r_c(1,1); % the row number; the reaction number.

    for c=1:r_c(1,2); % the column number
        % numOfLine=numOfLine+1;
        %disp(numOfLine);

        if isnumeric(ref.number(r,c))&&ref.number(r,c)>0;
            num=ref.number(r,c) % the pre-identified row number
            %%% <read>

            %%
            for i=1:length(items); % the first round sets 'width'; the second round sets 'color'.

                subText=ref.(items{i})(r,c);
                %%% <write>
                reText=cellstr(text(num).str)
                %%% find location and replace it

             ToBeRplaced=position(reText,items{i})
                %%%
                %try
                if ~isstr(subText{1})  % test if subText{1} contains a string
                   warning('OK');
                   disp(subText{1});
                   subText{1}=num2str(subText{1}) % convert the double type into string type
                end

                text(num).str=strrep(text(num).str,ToBeRplaced,subText{1}) %% or subText{1}, exchangeable here.
                %catch

                disp(subText{1});

              % end

                % disp([numOfLine,subText,ToBeRplaced,text(numOfLine).str]);

                %disp(text(num).str);
                %disp(subText);
            end

        end
    end
end

annotedText=text;
parsed_update=parsed;
parsed_update.r_info=ref; % update r_info structure
parsed_update.r_info.XMLtext=text; % update XMLtext

%%%%%%%%%%%%%%%%%%%%%%%%  write to the file

for d=1:length(text);

    fprintf(f_out,'%s\n',char(text(d).str));

end

fclose(f_out);

end

author: opencobra

repo: cobratoolbox


text:
function out = read_tobii_tsv(filename, varargin)

% READ_TOBII_TSV
%
% Use as
%   hdr = read_tobii_tsv(filename)
% or
%   dat = read_tobii_tsv(filename, tsv, begsample, endsample)

needhdr = (nargin==1);
needdat = (nargin>1);

if needhdr
  % get the original header details

  fid = fopen_or_error(filename, 'rt');
  tsv = struct();
  ln  = 0;
  while ln<30 && ~feof(fid)
    ln   = ln+1;
    line = fgetl(fid);
    if any(line==':') && line(end)~=':'
      % this looks like a header line
      [key, val] = strtok(line, ':');
      key = fixname(key);
      val = val(2:end); % skip the ':'
      val = strtrim(val);
      tsv.(key) = val;
    end

    if ~isempty(strfind(line, 'Timestamp')) && ~isempty(strfind(line, 'Gaze'))
      % this looks like the line with the headings of all the columns
      tsv.header = tokenize(line, '\t');
      % remember where the data starts
      tsv.headerline  = ln;
      tsv.datapointer = ftell(fid);
    end

  end

  fclose(fid);

  hdr = [];
  hdr.orig = tsv;

  % return the FieldTrip header information
  out = hdr;

elseif needdat

  hdr       = varargin{1};
  begsample = varargin{2};
  endsample = varargin{3};

  % get the original header details
  tsv = hdr.orig;
  fid = fopen_or_error(filename, 'rt');
  fseek(fid, tsv.datapointer, 'bof');

  cursample = 1;
  while cursample<begsample
    line = fgetl(fid);
    curline = curline + 1;
  end
  dat = zeros(length(tsv.header), endsample-begsample+1);
  while cursample<endsample
    line = fgetl(fid);
    val = tokenize(line, 9); % horizontal tab is 9 in ascii table
    num = cellfun(@str2num, val, 'UniformOutput', false);
    valid = cellfun(@isempty, num, 'UniformOutput', true);
    val(~str) = num;
    val( str) = nan;
    dat(:,cursample-begsample+1) = val;
    curline = curline + 1;
  end
  fclose(fid);

  % return the data
  out = hdr;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [yes] = isfilextension(fname, extension)
% See also CHECK_FILE_EXTENSION, CLEAR_FILE_EXTENSION.
%
% File:      isfilextension.m
% Author:    Ioannis Filippidis, jfilippidis@gmail.com
% Date:      2012.06.22
% Language:  MATLAB R2012a
% Purpose:   check if file has the desired extension
% Copyright: Ioannis Filippidis, 2012-

% dot in extension ?
if ~strcmp(extension(1), '.')
    extension = ['.', extension];
end

[~, ~, ext] = fileparts(fname);

% extension = unwanted ?
if strcmp(ext, extension)
    yes = 1;
else
    yes = 0;
end

author: Sable

repo: mcbench-benchmarks


text:
function uknt = kntunclamp (knt, deg, k, dim)
% KNTUNCLAMP: Compute the unclamped knot vector starting from an open one.
%
% Calling Sequence:
%
%   uknt = kntunclamp (knt, deg, k)
%   uknt = kntunclamp (knt, deg, k, dim)
%
% INPUT:
%
%   knt	: open knot vector: see kntrefine
%   deg : polynomial degree of the spline space
%   k   : continuity for the unclamping (from 0 up to p-1)
%   dim : dimensions in which to unclamp (all by default).
%
% OUTPUT:
%
%   uknt: unclamped knot vector, see nrbmak
%
% Description:
%
%     Unclamps directly the open knot vector. See nrbunclamp
%    for further information.
%
%    Copyright (C) 2013, 2014 Rafael Vazquez
%    Copyright (C) 2020, Bernard Kapidani
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

  knt_is_cell = true;
  if (~iscell (knt))
    knt = {knt};
    knt_is_cell = false;
  end
  uknt = knt;

  ndim = numel (knt);
  if (nargin < 4)
    dim = 1:ndim;
  end

% if (iscell (knt))
  if (numel(k) < ndim)
    k = [k(:).', k(end) * ones(1, ndim-numel(k))];
  end

  assert (numel(deg) == ndim, 'degrees and knots must have the same size');
  for idim = dim

    U  = knt{idim};

    p  = deg(idim);
    n  = numel(U) - p - 1;
    m  = n + p + 1;
    kk = k(idim);

    if (kk >= p)
      warning ('Taking the maximum k allowed, degree - 1')
      kk = p - 1;
    end

  % Unclamp at left end
    for ii=0:kk
      U(kk-ii+1) = U(kk-ii+2) - (U(n+1-ii) - U(n-ii));
    end

  % Unclamp at right end
    for ii=0:kk
      U(m-kk+ii) = U(m-kk+ii-1) + U(p+ii+1+1) - U(p+ii+1);
    end

    uknt{idim} = U;

  end

  if (~knt_is_cell)
    uknt = uknt{1};
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [result_edgesize,result_pertg,nvell_num] = Riverflux_distribution(obj)
% [result_edgesize,result_pertg，nvell_num] = Riverflux_distribution(obj)
%
% Input a msh class object to get the representative edge width and flux
% percentage for each node on the riverine flow boundary.
%
% This function is used to distribute a total flux of a cross-section where
% the riverine boundary located to each node of the riverine boundary.
%
% result_edgesize and result_pertgis are the result of the representative
% edge width and flux percentage for each node, respectively. nvell_num
% represents the number of nodes for each riverine boundary.
%
% Each column of result_edgesize and result_pertg represents the result of
% each riverine boundary.
%
% Nodal representative edge width equals to the sum of half the width of
% each of the two edges it is connected to.
%
% Nodeal flow percentage on the riverine boundary equals to the flow area
% of this node divided by the total flow area of the cross-section.
%
% Nodal flow area is calculated via a trapezoidal rule using its representative
% edge width and the bathymetric depths of this node and its neighboring nodes.
%
% Note that the calculation of representative edge width and flow percentage
% for the nodes at either end of the boundary is a little bit of different.
%
% User need to specify river flow boundaries (ibtype=22) before using it.
%
% Author:      Jiangchao Qiu
% Created:     January 7 2021
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bd_dat = obj.bd;
p_dat = obj.p;
b_dat = obj.b;
river_num = find(bd_dat.ibtype == 22);
N = length(river_num);% total number of river boundarys
nvell_num = bd_dat.nvell(river_num);% number of nodes for each river boundary

% Consider the riverine boundaries may have different number of nodes, the
% number of rows for result_edgesize and result_pertg is set to max(nvell_num)
result_edgesize = NaN(max(nvell_num),N); % the final result of edgesize
result_pertg = NaN(max(nvell_num),N);% the final result of percentage

if isempty(river_num)
    error('No riverine boundary information to distribute total flow')
end

for i=1:N
    node_num = bd_dat.nbvv(:,river_num(i));
    node_num(node_num==0)=[];
    bathy = b_dat(node_num);
    location = p_dat(node_num,:);
    J = length(node_num);% total number of nodes for the current river boundary
    flowarea = zeros(J,1);
    edgesize = zeros(J,1);
    %% calculate the projected distance for each edge on the boundary
    proj = 'Mercator';
    m_proj(proj,'lon',[ min(obj.p(:,1))-0.25 max(obj.p(:,1))+0.25 ],...
                'lat',[ min(obj.p(:,2))-0.25 max(obj.p(:,2))+0.25])
    proj_location = zeros(J,2);
    for j=1:J
        [proj_location(j,1),proj_location(j,2)] = m_ll2xy(location(j,1),location(j,2));
    end
    node_distance = 1000*m_xydist(proj_location(:,1),proj_location(:,2));
    %% calculate the respresentive width for each node
    % note: if the current boundary has J nodes, there are J-1 edges
    for j=1:J
		if j==1      % the first edge (the number is J)
            edgesize(j) = 0.5*node_distance(j);
        elseif j==J  % the last edge (the number is J-1)
            edgesize(j) = 0.5*node_distance(j-1);
        else         % the other edges
            edgesize(j) = 0.5*node_distance(j-1)+0.5*node_distance(j);
        end
        result_edgesize(1:J,i) = edgesize;
    end
    %% calculate flow area for each node on the boundary
    for j=1:J
        if j==1      % the first node
            local_z1 = (bathy(j)+bathy(j+1))*0.5;
            flowarea(j) = 0.5*(bathy(j)+local_z1)*(0.5*node_distance(j));
        elseif j==J  % the last node
            local_z1 = (bathy(j-1)+bathy(j))*0.5;
            flowarea(j) = 0.5*(bathy(j)+local_z1)*(0.5*node_distance(j-1));
        else         % the other nodes
            local_z1 = (bathy(j-1)+bathy(j))*0.5;
            local_z2 = (bathy(j)+bathy(j+1))*0.5;
            flowarea(j) = 0.5*(bathy(j)+local_z1)*(0.5*node_distance(j-1))+...
                          0.5*(bathy(j)+local_z2)*(0.5*node_distance(j));
        end
        percent = flowarea/sum(flowarea);  % the flux distribution percentage for each node
        result_pertg(1:J,i) =  percent;
    end
end

author: CHLNDDEV


text:
function[varargout]=transmaxdist(varargin)
%TRANSMAXDIST  Distributions of wavelet transform maxima in noise.
%
%   This function is part of 'element analysis' described in Lilly (2017),
%   "Element analysis: a wavelet-based method for analyzing time-localized
%   events in noisy time series", available at www.jmlilly.net.
%
%   [COUNT,BINS]=TRANSMAXDIST(GAMMA,BETA,ALPHA,FS,R,N,M) returns the
%   histogram of wavelet transform maxima magnitudes, for a length N time
%   series having spectral slope -2*ALPHA transformed at frequencies FS
%   using a (GAMMA,BETA) wavelet, based on a simulation having N*M points.
%
%   Here GAMMA, BETA, ALPHA, and R are all scalars, or are all arrays of
%   the same length as FS.  FS is a frequency array computed by MORSESPACE.
%
%   R is the ratio between each frequency FS and the next.  This will be
%   constant and greater than one when FS is computed by MORSESPACE.  As
%   as described in Appendix C of Lilly (2017), R=FS(n)./FS(n+1) for all n.
%
%   COUNT is the number of transform maxima observed at each frequency in
%   the magnitude bins BINS.  COUNT is a LENGTH(BINS) x LENGTH(FS) matrix.
%
%   Transform maxima values, as output in BINS, are normalized such that
%   the expected squared magnitude of the wavelet transform of noise occurs
%   at unity.  BINS thus corresponds to the normalized event magnitude.
%
%   TRANSMAXDIST works by simulating a vector whose statistical properties
%   mimic those of the wavelet transform and the four adjacted points, thus
%   avoiding the need to explicitly compute the transform.  The choice of
%   e.g. M=1000 simulates a transform 1000 times as long as time series of
%   interest, which itself is of length N.
%
%   [COUNT,BINS,RATE]=TRANSMAXDIST(...) also returns the RATE, the
%   normalized reversed cumulative density function.  RATE gives the
%   expected number of transform maxima occuring in a time series of length
%   N having a magnitude greater than the corresponding bin value.
%
%   [COUNT,BINS,RATE,SIGMA]=TRANSMAXDIST(...) also returns the theoretical
%   covariance matrix SIGMA from which the Monte Carlo simulations are
%   constructed.  SIGMA is an array of length 5 x 5 x LENGTH(FS).
%
%   Note that if the covariance matrix is not positive definite, as can
%   happen due to numerical complications for extreme BETA and GAMMA
%   choices, then COUNT and RATE will both consist entirely of NaNs.
%   _______________________________________________________________________
%
%   Additional options
%
%   TRANSMAXDIST(...,BINS) alternately uses BINS for the bin centers
%   instead of the default choice, which is set to LINSPACE(0,6,200)'.
%
%   By default, TRANSMAXDIST performs a simulation for each of the scale
%   frequencies in FS.  TRANSMAXDIST(...,'extrapolate') instead computes
%   the distribution only for the highest scale frequency, then
%   extrapolates these values to all other scale frequencies with a scaling
%   law.  GAMMA, BETA, ALPHA, and R must all be scalars in this case.
%
%   For details, see Lilly (2017).
%
%   See also MAXPROPS, TRANSMAX, ISOMAX, MAX2EDDY.
%
%   'transmaxdist --t' runs some tests.
%
%   Usage: [count,bins]=transmaxdist(ga,be,al,fs,r,N,M);
%          [count,bins,rate,sigma]=transmaxdist(ga,be,al,fs,r,N,M);
%          [count,bins,rate,sigma]=transmaxdist(ga,be,al,fs,r,N,M,'extrap');
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2017 J.M. Lilly --- type 'help jlab_license' for details

if strcmp(varargin{1}, '--t')
    transmaxdist_test,return
end

%parstr='serial';
%
%   TRANSMAXDIST(...,'parallel') uses a PARFOR loop in the computation of
%   the covariance matrix SIGMA.  The default behavior is 'serial'.

str='all';
normstr='band';
for i=1:3
    if ischar(varargin{end})
        if strcmpi(varargin{end}(1:3),'all')||strcmpi(varargin{end}(1:3),'ext')
            str=varargin{end};
        elseif strcmpi(varargin{end}(1:3),'ser')||strcmpi(varargin{end}(1:3),'par')
            parstr=varargin{end};
        elseif strcmpi(varargin{end}(1:3),'ban')||strcmpi(varargin{end}(1:3),'ene')
            normstr=varargin{end};
        end
        varargin=varargin(1:end-1);
    end
end

gamma=varargin{1};
beta=varargin{2};
alpha=varargin{3};
fs=varargin{4};
r=varargin{5};
N=varargin{6};
M=varargin{7};
if length(varargin)==7
    bins=linspace(0,6,200)';
else
    bins=varargin{8}(:);
end

s=morsefreq(gamma,beta)./fs;

%Some error checking
if ~isscalar(gamma)
    lg=length(gamma);
    lb=length(beta);
    la=length(alpha);
    ls=length(s);
    lr=length(r);
    if ~(lg==lb&&lg==la&&lg==ls&&lg==lr)
        error('TRANSMAXDIST was expecting the first five input arguments to all be the same length.')
    end
    if strcmpi(str(1:3),'ext')&&(length(gamma)>1)
        error('Sorry, ALPHA, BETA, GAMMA, FS, and R must be scalars for ''extrapolate'' option.')
    end
else
    if length(s)>1
        ro=s(2:end)./s(1:end-1);
        if ~allall(abs(ro-r)<1e8)
            disp('Input value of R does not match that computed from FS.');
        end
    end
end

arrayify(s,r,gamma,beta,alpha);
Sigma=zeros(5,5,length(s));

%s,r,gamma,beta,alpha
%if strcmpi(parstr(1:3),'ser')
    for i=1:length(s)
        Sigma(:,:,i)=sigmamat1(s(i),r(i),gamma(i),beta(i),alpha(i),normstr);
    end
%elseif strcmpi(parstr(1:3),'par')
%    parfor i=1:length(s)
%        Sigma(:,:,i)=sigmamat1(s(i),r(i),gamma(i),beta(i),alpha(i),normstr);
%    end
%end

xeps=(randn(5,N*M)+1i*randn(5,N*M))./sqrt(2);%std(xeps(:))
n=zeros(length(bins),length(s));

if strcmpi(str(1:3),'all')
    for k=1:length(s)
        disp(['TRANSMAXDIST performing simulation ' int2str(k) ' of ' int2str(length(s)) '.'])
        try
           L=chol(Sigma(:,:,k),'lower');
        catch
           L=[];
        end

        if ~isempty(L)
            wsim=L*xeps;

            bool=(abs(wsim(1,:))>abs(wsim(2,:)));
            bool=bool&(abs(wsim(1,:))>abs(wsim(3,:)));
            bool=bool&(abs(wsim(1,:))>abs(wsim(4,:)));
            bool=bool&(abs(wsim(1,:))>abs(wsim(5,:)));

            wmax=wsim(1,bool);
            nk=hist(abs(wmax),bins);
            n(:,k)=nk';
        else
            n(:,k)=nan*ones(size(n(:,k)));
        end
    end
elseif strcmpi(str(1:3),'ext')
    L=chol(Sigma(:,:,1),'lower');
    wsim=L*xeps;

    bool=(abs(wsim(1,:))>abs(wsim(2,:)));
    bool=bool&(abs(wsim(1,:))>abs(wsim(3,:)));
    bool=bool&(abs(wsim(1,:))>abs(wsim(4,:)));
    bool=bool&(abs(wsim(1,:))>abs(wsim(5,:)));

    wmax=wsim(1,bool);
    nk=hist(abs(wmax),bins);
    n(:,1)=nk';
    for k=2:length(s)
        n(:,k)=n(:,1)*frac(s(1),s(k));
    end
end

n=n./M;
rate=cumsum(n,1,'reverse');

%rate=n;
%for i=1:size(n,2)
    %L=2*sqrt(2)*sqrt(gamma(i).*beta(i))./fs(i);
 %   rate(:,i)=cumsum(n(:,i),1,'reverse');
%end

%     for i=1:5
%         for j=1:5
%             Sigmahat(i,j)=vmean(wsim(i,:).*conj(wsim(j,:)),2);
%         end
%     end
%     Sigmahat=Sigmahat./vmean(squared(wsim(1,:)),2);

varargout{1}=n;
varargout{2}=bins;
varargout{3}=rate;
varargout{4}=Sigma;
%varargout{4}=Sigmahat;

function[Sigma]=sigmamat1(s,r,gamma,beta,alpha,normstr)

p={gamma,beta,alpha,normstr};
Sigma=zeros(5,5);

Sigma(1,:)=[xi(0,s,1,p) xi(1,s,1,p)  xi(-1,s,1,p)  xi(0,s,r,p)   xi(0,s,1/r,p)     ];
Sigma(2,:)=[nan         xi(0,s,1,p)  xi(-2,s,1,p)  xi(-1,s,r,p)  xi(-1,s,1/r,p)    ];
Sigma(3,:)=[nan         nan          xi(0,s,1,p)   xi(1,s,r,p)   xi(1,s,1/r,p)     ];
Sigma(4,:)=[nan         nan          nan           xi(0,r*s,1,p) xi(0,r*s,1/r.^2,p)];
Sigma(5,:)=[nan         nan          nan           nan           xi(0,s/r,1,p)     ];

for j=1:5
    for i=(j+1):5
          Sigma(i,j)=conj(Sigma(j,i));
    end
end

[m0,ffun]=morsemom(-2*alpha,gamma,beta);
if strcmpi(normstr(1:3),'ene')
    Sigma=frac(Sigma,ffun.*s.^(2*alpha));
elseif strcmpi(normstr(1:3),'ban')
    Sigma=frac(Sigma,ffun*s.^(2*alpha-1));
end

function[x]=xi(tau,s,r,p)

gamma=p{1};
beta=p{2};
alpha=p{3};
normstr=p{4};

rtilde=(1+r.^gamma).^(1./gamma);

%Implicitly set A^2=1 as I will shortly divide by it
%fact1=frac(morseafun(gamma,beta).^2,morseafun(gamma,2*beta-2*alpha));
%fact2=frac((r.^beta).*(s.^(2*alpha-1)),rtilde.^(2*beta-2*alpha+1));

%Note:  use the input normalization in the numerator, and the *amplitude*
%normalization in the denominator.  All the latter does is cancel the
%coefficient coming from MORSEXPAND, which assumes the amplitude normalization
fact1=frac(morseafun(gamma,beta,normstr).^2,morseafun(gamma,2*beta-2*alpha));
if strcmpi(normstr(1:3),'ene')
    fact2=frac((r.^(beta+1/2)).*(s.^(2*alpha)),rtilde.^(2*beta-2*alpha+1));
elseif strcmpi(normstr(1:3),'ban')
    fact2=frac((r.^beta).*(s.^(2*alpha-1)),rtilde.^(2*beta-2*alpha+1));
end

x=fact1.*fact2.*conj(morsexpand(tau./(s.*rtilde),gamma,2*beta-2*alpha,morsefreq(gamma,2*beta-2*alpha)));
%x=fact1.*fact2.*conj(morsexpand(tau./(s.*rtilde),gamma,2*beta-2*alpha,morsefreq(gamma,2*beta-2*alpha),'cumulant'));

%fact1,fact2
%Equivalent
%fact2=frac((r.^beta).*(s.^(2*alpha)),rtilde.^(2*beta-2*alpha));
%fs=morsefreq(gamma,2*beta-2*alpha)./(s.*rtilde);
%x=fact1.*fact2.*conj(morsexpand(tau,gamma,2*beta-2*alpha,fs));

function[]=transmaxdist_test

for k=1:2
    tic
    ga=2;be=2;
    switch k
        case 1
            alpha=0;
        case 2
            alpha=1;
    end
    %fs=morsespace(ga,be,100);
    fs=morsespace(ga,be,{0.01,pi},2*pi/100);
    N=1e6;
    rng(1);
    x=randn(N,1);
    if alpha ==1
        x=cumsum(x);
        x=x./std(x);
    end
    w=wavetrans(x,{ga,be,fs(1:3),'bandpass'});

    clear xvec
    xvec(:,1)=w(:,2);
    xvec(:,2)=w([2:end 1],2);
    xvec(:,3)=w([end 1:end-1],2);
    xvec(:,4)=w(:,3);
    xvec(:,5)=w(:,1);

    for i=1:5
        for j=1:5
            Sigmahat(i,j)=vmean(xvec(:,i).*conj(xvec(:,j)),1);
        end
    end
    Sigmahat=Sigmahat./vmean(squared(xvec(:,1)),1);

    [count,bins,rate,Sigma]=transmaxdist(ga,be,alpha,fs(2),fs(1)./fs(2),N,1);
    switch k
        case 1
            bool=maxmax(abs(Sigma-Sigmahat)./abs(Sigma))<1/100;
            reporttest('TRANSMAXDIST simulated and theoretical covariance matrices agree to within 1%, white noise case',bool)
        case 2
            bool=maxmax(abs(Sigma-Sigmahat)./abs(Sigma))<4/100;
            reporttest('TRANSMAXDIST simulated and theoretical covariance matrices agree to within 4%, red noise case',bool)
    end

    %[index,ww]=transmax(fs(1:3),w);
    %n=hist(abs(ww)./vstd(w(:,2),1),bins)';
    %figure,plot(bins,[count n]) %That looks great!
end

function[]=transmaxdist_other

%[n,sigma]=transmaxdist(ga,be,0,[fs(2) fs(3)],N/10);
N=1e8;
tic;[n,x,Sigma]=transmaxdist(ga,be,0,fs,N);toc
for i=1:size(n,2)
    ntilde(:,i)=n(:,i)./frac(N,2*sqrt(2)*sqrt(ga*be)./fs(i));
    %ntilde(:,i)=n(:,i)./sum(n(:,i));
end

tic
N=1e8;
x=randn(N,1);
w=wavetrans(x,{ga,be,fs(1:3)});
toc

tic; [ii,jj,ww,ff]=transmax(fs(1:3),w);toc
bins=linspace(0,6,200)';
wwtilde=ww./sqrt(vmean(squared(w(:,2)),1));
nn=hist(abs(wwtilde),bins);
nn=nn./frac(N,2*sqrt(2)*sqrt(ga*be)./fs(2));

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


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

%%
clear; close all; clc;

%% Syntax
% |[cMap]=igviridis(varargin);|

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


text:
function a = wathen ( nx, ny, n )

%*****************************************************************************80
%
%% WATHEN returns the WATHEN matrix.
%
%  Discussion:
%
%    The Wathen matrix is a finite element matrix which is sparse.
%
%    The entries of the matrix depend in part on a physical quantity
%    related to density.  That density is here assigned random values between
%    0 and 100.
%
%    The matrix order N is determined by the input quantities NX and NY,
%    which would usually be the number of elements in the X and Y directions.
%    The value of N is
%
%      N = 3*NX*NY + 2*NX + 2*NY + 1,
%
%    and sufficient storage in A must have been set aside to hold
%    the matrix.
%
%    A is the consistent mass matrix for a regular NX by NY grid
%    of 8 node serendipity elements.
%
%    The local element numbering is
%
%      3--2--1
%      |     |
%      4     8
%      |     |
%      5--6--7
%
%    Here is an illustration for NX = 3, NY = 2:
%
%     23-24-25-26-27-28-29
%      |     |     |     |
%     19    20    21    22
%      |     |     |     |
%     12-13-14-15-16-17-18
%      |     |     |     |
%      8     9    10    11
%      |     |     |     |
%      1--2--3--4--5--6--7
%
%    For this example, the total number of nodes is, as expected,
%
%      N = 3 * 3 * 2 + 2 * 2 + 2 * 3 + 1 = 29
%
%  Properties:
%
%    A is symmetric positive definite for any positive values of the
%    density RHO(X,Y), which is here given the value 1.
%
%    The problem could be reprogrammed so that RHO is nonconstant,
%    but positive.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Nicholas Higham,
%    Algorithm 694: A Collection of Test Matrices in MATLAB,
%    ACM Transactions on Mathematical Software,
%    Volume 17, Number 3, September 1991, pages 289-305.
%
%    Andrew Wathen,
%    Realistic eigenvalue bounds for the Galerkin mass matrix,
%    IMA Journal of Numerical Analysis,
%    Volume 7, Number 4, October 1987, pages 449-457.
%
%  Parameters:
%
%    Input, integer NX, NY, values which determine the size of A.
%
%    Input, integer N, the order of N, as determined by NX and NY.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  em =  [ ...
     6.0, -6.0,  2.0, -8.0,  3.0, -8.0,  2.0, -6.0; ...
    -6.0, 32.0, -6.0, 20.0, -8.0, 16.0, -8.0, 20.0; ...
     2.0, -6.0,  6.0, -6.0,  2.0, -8.0,  3.0, -8.0; ...
    -8.0, 20.0, -6.0, 32.0, -6.0, 20.0, -8.0, 16.0; ...
     3.0, -8.0,  2.0, -6.0,  6.0, -6.0,  2.0, -8.0; ...
    -8.0, 16.0, -8.0, 20.0, -6.0, 32.0, -6.0, 20.0; ...
     2.0, -8.0,  3.0, -8.0,  2.0, -6.0,  6.0, -6.0; ...
    -6.0, 20.0, -8.0, 16.0, -8.0, 20.0, -6.0, 32.0 ];

  for j = 1 : ny

    for i = 1 : nx
%
%  For the element (I,J), determine the indices of the 8 nodes.
%
      node(1) = 3 * j * nx + 2 * j + 2 * i + 1;
      node(2) = node(1) - 1;
      node(3) = node(1) - 2;

      node(4) = ( 3 * j - 1 ) * nx + 2 * j + i - 1;
      node(8) = node(4) + 1;

      node(5) = ( 3 * j - 3 ) * nx + 2 * j + 2 * i - 3;
      node(6) = node(5) + 1;
      node(7) = node(5) + 2;

      rho = 0.5;
%     [ rho, seed ] = r8_uniform_01 ( seed );
      rho = 100.0 * rho;

      for krow = 1 : 8
        for kcol = 1 : 8

          a(node(krow),node(kcol)) = a(node(krow),node(kcol)) ...
            + rho * em(krow,kcol);

        end
      end

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef hdf5prop < handle
% HDF5PROP Class for transparent file data access.
% Matlab class to create and access HDF5 datasets transparently as Matlab
% variable. Data can be accessed and written with subscript referencing and
% assignment methods, just like a matlab variable, only size must be
% explicitly set or changed.
%
%   prop = hdf5prop(file, dataset)
%     Creates a 'hdf5prop' object for the given dataset in the given HDF5
%     file.
%
%   prop = hdf5prop(file, dataset, mode)
%     Specify the access mode,
%       mode = 'r' : readonly access (default)
%       mode = 'rw': readwrite access
%
%   prop = hdf5prop(file, dataset, parameter, value, ...)
%     Creates a 'hdf5prop' object for a new dataset with given parameters.
%
%     parameter
%       'size'        - size of the dataset. Necessary.
%       'chunk_size'  - size of the chunks of the dataset.
%       'compression' - 0-9  compression level.
%       'type'        - type of the dataset to create.  This can be a string
%                       such as 'double', in which case the datatype maps to
%                       'H5T_NATIVE_DOUBLE', or it can be a derived HDF5
%                       datatype.
%       'max_size'    - the maximum size of the dataset to create.
%       'fill'        - the fill value.

% copyleft 2011, Piers Titus van der Torren

  % TODO: replace all H5 calls to direct mex calls for speed
  properties( SetAccess = protected )
    file
    dataset
    mode = 'H5F_ACC_RDONLY'
    propsize
  end
  properties( SetAccess = protected, Transient = true )
    dataset_id
  end

  methods
    function self = hdf5prop(file,dataset,varargin)

      self.file = file;
      self.dataset = dataset;
      if nargin>2
        switch varargin{1}
          case {'H5F_ACC_RDWR','rw'}
            self.mode = 'H5F_ACC_RDWR';
            varargin = varargin(2:end);
          case {'H5F_ACC_RDONLY','r'}
            varargin = varargin(2:end);
          otherwise
        end
      end
      if ~exist(file, 'file')
        % create file
        file_id = H5F.create(file, 'H5F_ACC_EXCL', 'H5P_DEFAULT', 'H5P_DEFAULT');
        self.mode = 'H5F_ACC_RDWR';
      else
        file_id = H5F.open(self.file, self.mode, 'H5P_DEFAULT');
      end
      try
        self.dataset_id = H5D.open(file_id,self.dataset);
        if ~isempty(varargin)
          error('hdf5prop:exists', 'dataset already exists, even though creation parameters are given')
        end
      catch ex
        if strcmp(ex.identifier, 'hdf5prop:exists')
          rethrow(ex);
        end
        % create dataset

        if strcmp(self.mode, 'H5F_ACC_RDONLY')
          self.mode = 'H5F_ACC_RDWR';
          H5F.close(file_id);
          file_id = H5F.open(self.file, self.mode, 'H5P_DEFAULT');
        end

        self.dataset_id = h5datacreate(file_id, self.dataset, varargin{:});

      end

      space = H5D.get_space(self.dataset_id);
      [ans, dims] = H5S.get_simple_extent_dims(space);
      self.propsize = fliplr(dims);
    end

    function open( self )
      if ~isa(self.dataset_id, 'H5ML.id') || self.dataset_id.identifier < 0

        file_id=H5F.open(self.file,self.mode,'H5P_DEFAULT');
        self.dataset_id=H5D.open(file_id,self.dataset);

        % get size of dataset
        % TODO: size is not valid when file is changed elsewhere while open
        space = H5D.get_space(self.dataset_id);
        [ndims, dims] = H5S.get_simple_extent_dims(space);
        self.propsize = fliplr(dims);
      end
    end

    function close( self )
      % CLOSE Close file if open.
      if isa(self.dataset_id, 'H5ML.id')
        H5D.close(self.dataset_id);
      end
    end

    function set_extent( self, sz)
      % SET_EXTENT Allocate space
      %
      %   hdf5prop.set_extent(sz)
      %     Extend allocated space to size sz.
      self.open();
      H5D.extend(self.dataset_id, fliplr(sz));
      %H5D.set_extent(self.dataset_id, fliplr(sz));
      self.propsize = sz;
    end

    function set_mode( self, mode )
      self.close()
      switch mode
        case {'H5F_ACC_RDWR','rw'}
          self.mode = 'H5F_ACC_RDWR';
        case {'H5F_ACC_RDONLY','r'}
          self.mode = 'H5F_ACC_RDONLY';
        otherwise
          error('unknown mode')
      end
      self.open()
    end

    function lock( self )
      self.close()
      self.mode = 'H5F_ACC_RDONLY';
      self.open()
    end
  end

  methods( Hidden=true )

    function varargout = subsref(self, S)
      switch S(1).type
        case {'.', '{}'}
          [varargout{1:nargout}] = builtin('subsref',self,S);
        case '()'
          [varargout{1:nargout}] = self.get_data(S(1).subs{:});
%         case '{}'
%           error('{} not allowed')
      end
    end

    function self = subsasgn(self, S, B)
      switch S(1).type
        case {'.', '{}'}
          builtin('subsasgn',self,S,B);
        case '()'
          self.set_data(S(1).subs, B);
%         case '{}'
%           error('{} not allowed')
      end
    end

    function out = get_data(self, varargin)

      self.open()

      % return all data
      if nargin == 1 || all(strcmp(varargin,':'))
        out = H5D.read(self.dataset_id,'H5ML_DEFAULT','H5S_ALL','H5S_ALL','H5P_DEFAULT');
        return;
      end

      outsize = [];

      propsize = size(self);
      dim = numel(propsize);
      % fill singular dimensions if there's only 1 nonsingular dimension
      if dim > 1 && numel(varargin)==1 && sum(propsize ~= 1) <= 1
        a = varargin;
        varargin = repmat({':'},1,dim);
        if any(propsize ~= 1)
          varargin{propsize ~= 1} = a{1};
        end
        outsize = size(a{1});
        if sum(outsize~=1)==1
          outsize = [];
        end
      end
      assert(dim == numel(varargin),'linear or multidimensional boolean indexing is not yet supported, now use all %d dimensions',dim);

      % divide selection into hyperslabs
      start = nan(1,dim);
      count = nan(1,dim);
      stride = nan(1,dim);
      hyperslab = true;
      for n = 1:dim
        sel = varargin{dim+1-n};  % traverse reversed instead of fliplr
        if isa(sel,'char') && sel == ':'
          start(n) = 1;
          count(n) = propsize(dim+1-n);
          stride(n) = 1;
          idxc{n} = sel;
          selc{n} = sel;
        else
          %if isempty(sel)
          %  % if one dimension is empty the result is empty too
          %  out = [];
          %  return
          %end

          sel = sel(:);
          if isa(sel,'logical')
            sel = find(sel);
            idxc{n} = ':';
          elseif numel(sel) <= 1 || ( issorted(sel) && min(diff(sel))>0 )
            idxc{n} = ':';
          else
            [sel, ans, idxc{n}] = unique(sel);
          end

          selc{n} = sel;

          count(n) = numel(sel);

          dsel = diff(sel);
          if isempty(sel)
            start(n) = 1;
            stride(n) = 1;
          elseif numel(sel)==1
            start(n) = sel(1);
            stride(n) = 1;
          elseif all(dsel == dsel(1))
            start(n) = sel(1);
            stride(n) = dsel(1);
          else
            hyperslab = false;
          end
        end
      end

      % TODO: fix datatype
      if any(count==0)
        out = zeros(fliplr(count));
        return
      end

      space=H5D.get_space(self.dataset_id);
      if hyperslab
        %H5ML.hdf5lib2('H5Sselect_hyperslab',space_double,'H5S_SELECT_SET',start,stride,count,[]);
        H5S.select_hyperslab(space,'H5S_SELECT_SET',start-1,stride,count,[]);
      else
        s = find(isnan(start));
        selend = cellfun(@(x) x(end),selc(s));
        sel1 = cellfun(@(x) x(1),selc(s));
        if prod(count(s)) > .1 * prod(selend - sel1)
          % since hyperslab is way faster than select_elements, it is
          % faster to read a full hypeslab and do the selection in
          % matlab. The factor above should match the speed ratio.
          start(s) =  sel1;
          stride(s) = 1;
          count(s) = selend-sel1+1;
          for n=s
            idxc{n} = selc{n}(idxc{n}) - selc{n}(1) + 1;
          end
          H5S.select_hyperslab(space,'H5S_SELECT_SET',start-1,stride,count,[]);
        else
          sel = zeros(dim,prod(count));
          % ndgrid of selc{:}
          for n=1:dim,
            if ischar(selc{n})
              x = (1:count(n))';
            else
              x = selc{n}; % Extract and reshape as a vector.
            end
            s = count; s(n) = []; % Remove i-th dimension
            x = reshape(x(:,ones(1,prod(s))),[length(x) s]); % Expand x
            x = permute(x,[2:dim-n+1 1 dim-n+2:dim]); % Permute to i'th dimension
            sel(n,:) = x(:);
          end
          %warning('using H5S.select_elements, might be slow')
          H5S.select_elements(space,'H5S_SELECT_SET',sel-1);
        end
      end

      % reserve space for reading
      mem_space = H5S.create_simple(dim, count, []);

      % read the actual data
      out = H5D.read(self.dataset_id,'H5ML_DEFAULT',mem_space,space,'H5P_DEFAULT');

      % shuffle output to match selection
      out=out(idxc{end:-1:1});

      if ~isempty(outsize)
        out = reshape(out,outsize);
      end

      H5S.close(space);
      H5S.close(mem_space);
    end

    function set_data(self, argin, B)

      if ~strcmp(self.mode, 'H5F_ACC_RDWR')
        error('this hdf5prop is read only')
      end

      self.open()

      propsize = size(self);
      dim = numel(propsize);
      % fill singular dimensions if there's only 1 nonsingular dimension
      if dim > 1 && numel(argin)==1 && sum(propsize ~= 1) <= 1
        a = argin;
        argin = repmat({':'},1,dim);
        if any(propsize ~= 1)
          argin{propsize ~= 1} = a{1};
        end
      end
      assert(dim == numel(argin),'linear or multidimensional boolean indexing is not yet supported, now use all %d dimensions',dim);

      colondim = strcmp(argin,':');

      % get input indexing dimensions
      indim = 1:dim;
      for n=1:dim
        if colondim(n)
          indim(n) = propsize(n);
        elseif isa(argin{n},'logical')
          indim(n) = sum(argin{n}(:));
        else
          indim(n) = numel(argin{n});
        end
      end

      % check if squeezed dimensions match
      indim = max(cellfun(@numel,argin), colondim.*propsize);
      szB = size(B);
      assert(numel(B)==1 || all(indim(indim~=1)==szB(szB~=1)),'Subscripted assignment dimension mismatch.')

      if any(indim==0)
        % nothing to be done
        return;
      end

      % unsqueeze input
      if numel(B)~=1
        B = reshape(B,indim);
      end

      % write all data
      if all(colondim)
        if numel(B) == 1
          B = repmat(B,propsize);
          % H5D.fill ?
        end

        H5D.write(self.dataset_id,'H5ML_DEFAULT','H5S_ALL','H5S_ALL','H5P_DEFAULT',B);
        return;
      end

      % divide selection into hyperslabs
      start = zeros(1,dim);
      count = zeros(1,dim);
      stride = zeros(1,dim);
      hyperslab = true;
      for n = 1:dim
        sel = argin{dim+1-n};  % traverse reversed instead of fliplr
        if isa(sel,'char') && sel == ':'
          start(n) = 1;
          count(n) = propsize(dim+1-n);
          stride(n) = 1;
          idxc{n} = sel;
          selc{n} = sel;
        else
          sel = sel(:);
          if isa(sel,'logical')
            sel = find(sel);
            idxc{n} = ':';
          elseif numel(sel) == 1 || ( issorted(sel) && min(diff(sel))>0 )
            idxc{n} = ':';
          else
            % this line is different from get_data
            [sel, idxc{n}] = unique(sel);
          end

          selc{n} = sel;

          count(n) = numel(sel);
          if hyperslab
            dsel = diff(sel);
            if numel(sel)==1
              start(n) = sel(1);
              stride(n) = 1;
            elseif all(dsel == dsel(1))
              start(n) = sel(1);
              stride(n) = dsel(1);
            else
              hyperslab = false;
            end
          end
        end
      end

      space=H5D.get_space(self.dataset_id);
      if hyperslab
        %H5ML.HDF5lib2('H5Sselect_hyperslab',space_double,'H5S_SELECT_SET',start,stride,count,[]);
        H5S.select_hyperslab(space,'H5S_SELECT_SET',start-1,stride,count,[]);
      else
        sel = zeros(dim,prod(count));
        % ndgrid of selc{:}
        for n=1:dim,
          if ischar(selc{n})
            x = (1:count(n))';
          else
            x = selc{n}; % Extract and reshape as a vector.
          end
          s = count; s(n) = []; % Remove i-th dimension
          x = reshape(x(:,ones(1,prod(s))),[length(x) s]); % Expand x
          x = permute(x,[2:dim-n+1 1 dim-n+2:dim]); % Permute to i'th dimension
          sel(n,:) = x(:);
        end
        H5S.select_elements(space,'H5S_SELECT_SET',sel-1);
      end

      % reserve space for reading
      mem_space = H5S.create_simple(dim, count, []);

      % prepare input
      if prod(count) > 1
        if numel(B) == 1
          B = repmat(B,count);
        else
          % shuffle input to match selection
          B = B(idxc{end:-1:1});
        end
      end

      % write the actual data
      H5D.write(self.dataset_id,'H5ML_DEFAULT',mem_space,space,'H5P_DEFAULT',B);

      H5S.close(space);
      H5S.close(mem_space);
    end

    function varargout = size(self, dim)
      self.open()
      sz = self.propsize;
      if nargout > 1
        varargout = num2cell(sz);
      elseif nargin > 1
        varargout = {sz(dim)};
      else
        varargout = {sz};
      end
    end

    % overloading numel causes problems with builtin subsref
%     function out = numel(self)
%       out = 1;%prod(size(self));
%     end

    function out = length(self)
      out = max(size(self));
    end

    function out = end(self, dim, n_dim)
      out = size(self,dim);
    end

    function out = vertcat(varargin)
      error('array of hdf5prop objects is not allowed')
    end

    function out = horzcat(varargin)
      error('array of hdf5prop objects is not allowed')
    end

    function out = cat(varargin)
      error('array of hdf5prop objects is not allowed')
    end

    function out = double( self )
      out = double(self.get_data());
    end

    function out = eq( self, other )
      out = eq(self.get_data(), other());
    end

    function out = ne( self, other )
      out = ne(self.get_data(), other());
    end

    function out = le( self, other )
      out = le(self.get_data(), other());
    end

    function out = ge( self, other )
      out = qe(self.get_data(), other());
    end

    function out = lt( self, other )
      out = lt(self.get_data(), other());
    end

    function out = gt( self, other )
      out = gt(self.get_data(), other());
    end
  end
end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function seed = randomseed(seed)
% RANDOMSEED   Get or set the random seed.
% SEED = RANDOMSEED returns the current random seed for lightspeed random
% number routines.  SEED is a vector of 3 integers.
% RANDOMSEED(NEW_SEED) sets a new random seed.
% The seed determines the sequence of numbers that will be generated.
% Only the routines randgamma, randbinom, and sample_hist are affected.
%
% Example:
% randomseed([4 5 6])
% randgamma(repmat(3,1,3))   % 3 random numbers
% randomseed([4 5 6])
% randgamma(repmat(3,1,3))   % the same 3 numbers

error('You must first run install_lightspeed.m');

author: andrewssobral

repo: lrslibrary