From bc0840d286030d05c76b889674a4d99a7494e247 Mon Sep 17 00:00:00 2001 From: Christian Date: Mon, 21 Jan 2019 12:53:49 +0100 Subject: [PATCH] release 2.0 --- CITATION | 17 +++ CODE | 18 +++ LICENSE | 2 +- README.md | 147 +++++++++++++++++------- RUNME.m | 150 +++++++++++++------------ hapod.m | 326 ++++++++++++++++++++++++++++-------------------------- 6 files changed, 390 insertions(+), 270 deletions(-) create mode 100644 CITATION create mode 100644 CODE diff --git a/CITATION b/CITATION new file mode 100644 index 0000000..cdee4b7 --- /dev/null +++ b/CITATION @@ -0,0 +1,17 @@ +Cite As: + +C. Himpe, T. Leibner and S. Rave. Hierarchical Approximate Proper Orthogonal Decomposition. +SIAM Journal on Scientific Computing, 40(5): A3267--A3292, 2018. doi:10.1137/16M1085413 + +BibTeX Entry: + +@ARTICLE{hapod, + author = {C. Himpe and T. Leibner and S. Rave}, + title = {Hierarchical Approximate Proper Orthogonal Decomposition}, + journal = {SIAM Journal on Scientific Computing}, + volume = {40}, + number = {5}, + pages = {A3267--A3292}, + year = {2018}, + doi = {10.1137/16M1085413} +} diff --git a/CODE b/CODE new file mode 100644 index 0000000..7aad6fe --- /dev/null +++ b/CODE @@ -0,0 +1,18 @@ +# code.ini +name: Hierarchical Approximate Proper Orthogonal Decomposition +shortname: hapod +version: 2.0 +release-date: 2019-01-21 +author: Christian Himpe, Stephan Rave +orcid: 0000-0003-2194-6754, 0000-0003-0439-7212 +topic: Science, Mathematics, Dimension Reduction +type: Toolbox +license: 2-Clause BSD +license-type: open-source +repository: github.com/gramian/hapod +repository-type: git +language: Matlab +dependencies: Octave >=4.2, Matlab >=2013b +systems: Linux, Windows +website: git.io/hapod +keywords: Singular Value Decomposition, Principal Component Analysis, Proper Orthogonal Decomposition diff --git a/LICENSE b/LICENSE index d1b4525..c2f697a 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ BSD 2-Clause License -Copyright (c) 2017--2018, Christian Himpe & Stephan Rave +Copyright (c) 2017--2019, Christian Himpe & Stephan Rave All rights reserved. Redistribution and use in source and binary forms, with or without diff --git a/README.md b/README.md index 8eafb37..d95c51a 100644 --- a/README.md +++ b/README.md @@ -1,73 +1,140 @@ HAPOD - Hierarchical Approximate Proper Orthogonal Decomposition ================================================================ -* HAPOD - Hierarchical Approximate POD ( http://git.io/hapod ) -* version: 1.3 ( 2018-02-09 ) +* HAPOD - Hierarchical Approximate POD +* version: 2.0 ( 2019-01-21 ) * by: Christian Himpe (0000-0003-2194-6754), Stephan Rave (0000-0003-0439-7212) * under: BSD 2-Clause License (open-source) * summary: Distributed or incremental POD / SVD computation -# Scope +## Scope -* (Truncated) Singular Value Decomposition (SVD) * Proper Orthogonal Decomposition (POD) -* Principal Compoenent Analysis (PCA) -* Empirical Orthogoanl Functions (EOF) +* (Truncated) Singular Value Decomposition (SVD) +* (Sparse) Principal Compoenent Analysis (PCA) +* Empirical Orthogonal Functions (EOF) * Empirical Eigenfunctions -* Karhunen Loeve Decomposition +* Karhunen-Loeve Decomposition +* Dimension Reduction * Model Reduction | Model Order Redction -# Features +## Features * Standard POD -* Incremental HAPOD -* Distributed HAPOD +* Incremental (HA)POD +* Distributed (HA)POD * Custom SVD backend +* Method of Snapshots + +## Compatibility + +* GNU Octave >= 4.2 +* Mathworks MATLAB >= 2013b -# Basic Usage +## Basic Usage ``` -[svec,sval,snfo] = hapod(data,bound,topo,relax,config,mysvd) +[svec,sval,snfo] = hapod(data,bound,type,relax,config,mysvd) ``` -## Arguments +## Documentation + +### Arguments + +* `data` {cell array} +* `bound` {scalar} L2 mean projection error bound +* `type` {string} HAPOD tree type + * `'incr'` Incremental HAPOD (Complete) + * `'incr_1'` Incremental HAPOD (Non-root node) + * `'incr_r'` Incremental HAPOD (Root node) + * `'dist'` Distributed HAPOD (Complete) + * `'dist_1'` Distributed HAPOD (Non-root node) + * `'dist_r'` Distributed HAPOD (Root node) + * `'none'` Standard POD +* `relax` {scalar} Relaxation parameter in (0,1] by default: 0.5. +* `config` {structure} Configuration structure by default empty. + * `nLevels` - Total number of levels in tree + * `nSnapshots` - Number of data snapshots (columns) seen at each node + * `nModes` - Number of modes resulting at each node + * `tNode` - Time consumed at each node +* `mysvd` {handle} Function handle to a custom POD method. + * by default an economic SVD is used. + * alternatively the method-of-snapshots can be used via `'mos'`. + * otherwise a function handle can be provided. + +### Return Values * `svec` {matrix} POD modes (column vectors) * `sval` {vector} Singular values * `snfo` {structure} Information structure -## Return Values +### Usage -* `data` {cell array} -* `bound` {scalar} L2 mean projection error bound -* `topo` {string} HAPOD graph topology - * `'none'` Standard POD - * `'incr'` Incremental HAPOD - * `'incr_0'` Incremental HAPOD (First node only) - * `'incr_1'` Incremental HAPOD (One intermediary node only) - * `'incr_r'` Incremental HAPOD (Root node only) - * `'dist'` Distributed HAPOD - * `'dist_0'` Distributed HAPOD (First node only) - * `'dist_1'` Distributed HAPOD (One intermediary node only) - * `'dist_r'` Distributed HAPOD (Root node only) -* `relax` {scalar} Relaxation parameter in [0,1] by default 0.5. -* `config` {structure} Configuration structure by default empty -* `mysvd` {handle} Function handle to a custom POD method - -# Documentation - -## Information and Configuration Structure - -* `nSnapshots` - Number of data columns passed to this hapod and its children. +If all data partitions can be passed as the data argument, the types: `none` +(standard POD), `incr`(emental) HAPOD or `dist`(ributed) HAPOD are applicable. +In case only a single partition can be passed, the types: `incr_1` and `dist_1` +should be used for the non-root nodes of the associated HAPOD tree, while the +types: `incr_r` and `dist_r` should be used for the root node. The returned +information structure (or a cell-array thereof) can be passed to the parent +nodes in the associated HAPOD tree. + +### Information and Configuration Structure + +* `nLevels` - Total number of levels in tree +* `nSnapshots` - Number of data columns passed to this hapod and its children. * `nModes` - Number of intermediate modes * `tNode` - computational time at this hapod's branch -## Custom SVD Backend +Only for `incr_1`, the number of levels `nLevels` in the tree needs to be +provided in a structure. The other fields are filled by the `hapod` function. -Signature, arguments and return values as for Matlab's svd function. +### Custom SVD Backend -# Cite As +Signature, arguments and return values as for MATLAB's svd function. + +``` +[U,D,V] = mysvd(X) +``` + +### Getting Started + +Run the sample code: + +``` +RUNME() +``` + +which demonstrates the different implemented HAPOD variants and can be used +as a template. + +## Cite As C. Himpe, T. Leibner and S. Rave. -"[Hierarchical Approximate Proper Orthogonal Decomposition](http://arxiv.org/abs/1607.05210)". -Preprint, arXiv math.NA: 1607.05210, 2017. +"[Hierarchical Approximate Proper Orthogonal Decomposition](https://doi.org/10.1137/16M1085413)". +SIAM Journal on Scientific Computing, 40(5): A3267--A3292, 2018. + +## Used In + +* P. Benner, C. Himpe. +"[Cross-Gramian-Based Dominant Subspaces](https://arxiv.org/abs/1809.08066)". +arXiv, math.OC: 1809.08066, 2018. + +* T. Taddei, A.T. Patera. +"[A Localization Strategy for Data Assimilation; Application to State Estimation and Parameter Estimation](https://doi.org/10.1137/17M1116830)". +SIAM Journal on Scientific Computing 40(2): B611--B636, 2018. + +* C. Himpe, T. Leibner, S. Rave. +"[HAPOD - Fast, Simple and Reliable Distributed POD Computation](https://doi.org/10.11128/arep.55.a55283)". +ARGESIM Report 55 (MATHMOD 2018 Volume): 119--120, 2018. + +* B.J. Beach. +"[An Implementation-Based Exploration of HAPOD: Hierarchical Approximate Proper Orthogonal Decomposition](http://hdl.handle.net/10919/81938)". +Master Thesis, Virginia Tech, 2018. + +* C. Himpe, T. Leibner, S. Rave, J. Saak. +"[Fast Low-Rank Empirical Cross Gramians](https://doi.org/10.1002/pamm.201710388)". +Proceedings in Applied Mathematics and Mechanics 17: 841--842, 2017. + +* C. Himpe, T. Leibner, S. Rave. +"[Comprehensive Memory-Bound Simulations on Single Board Computers](https://doi.org/10.5281/zenodo.814497)". +Extended Abstract, 2nd Conference on Power Aware Computing (PACO), 2017. diff --git a/RUNME.m b/RUNME.m index 33bef6b..e68ef2e 100644 --- a/RUNME.m +++ b/RUNME.m @@ -1,101 +1,109 @@ function RUNME() -%%% project: hapod - Hierarchical Approximate POD ( http://git.io/hapod ) -%%% version: 1.3 ( 2018-02-09 ) +%%% project: hapod - Hierarchical Approximate POD ( https://git.io/hapod ) +%%% version: 2.0 ( 2019-01-21 ) %%% authors: C. Himpe ( 0000-0003-2194-6754 ), S. Rave ( 0000-0003-0439-7212 ) %%% license: BSD 2-Clause License ( opensource.org/licenses/BSD-2-Clause ) %%% summary: basic test for incremental HAPOD and distributed HAPOD -% %% Generate test data - randn('seed',1009); - n = 16; - [a,b,c] = svd(randn(n*n,n*n)); - s = a*diag(logspace(0,-16,n*n))*b'; - S = mat2cell(s,size(s,1),n*ones(n,1)); - - w = 0.1; - E = 1e-8; + randn('seed',1009); % seed random number generator + n = 16; % set number of partitions + N = n*n; % set test problem size + [a,b,c] = svd(randn(N,N)); % compute SVD of normal random matrix + s = a*diag(logspace(0,-16,N))*b'; % reassign singular values + S = mat2cell(s,size(s,1),n*ones(n,1)); % split data matrix into partitions + w = 0.1; % relaxation parameter omega + E = 1e-8; % target mean L2 projection error + disp(' '); % Compute default POD + disp('REFERENCE POD:'); [U,D,C] = hapod(S,E,'none'); - PRESCRIBED_MEAN_L2 = E - POD_MEAN_L2 = norm(s-U*U'*s,'fro')/sqrt(n*n) - disp(''); - -%% Incremental HAPOD - - % Compute global mode bound POD for incremental HAPOD [ug,dg,cg] = hapod(S,E*w,'none'); - iHAPOD_GLOBAL_MODE_BOUND = size(ug,2) + prescribed_mean_L2 = E + mean_L2 = norm(s-U*U'*s,'fro')/sqrt(N) + global_mode_bound = size(ug,2) + disp(' '); - % Compute local mode bound POD for incremental HAPOD - [ul,dl,cl] = hapod(S,E*sqrt(1.0-w^2)/sqrt(n*n-1),'none'); - iHAPOD_LOCAL_MODE_BOUND = size(ul,2) - disp(''); +%% Incremental HAPOD - % Test full incremental HAPOD + % Test bulk incremental HAPOD + disp('INCREMENTAL HAPOD (BULK):'); [U,D,C] = hapod(S,E,'incr',w); - iHAPOD_MEAN_L2 = norm(s-U*U'*s,'fro')/sqrt(n*n) - iHAPOD_MODES = size(U,2) - iHAPOD_MAX_MODES = max(C.nModes) - disp(''); + mean_L2 = norm(s-U*U'*s,'fro')/sqrt(N) + num_global_modes = size(U,2) + max_local_modes = max(cell2mat(C.nModes(1:end-1))) + disp(' '); - % Test partial incremental HAPOD - c1.nSets = n; - [u1,d1,c1] = hapod(S(1),E,'incr_0',w,c1); + % Test chunk incremental HAPOD + disp('INCREMENTAL HAPOD (CHUNK):'); + u1 = []; + c1 = struct('nLevels',n); - for k=2:n-1 + for k=1:n-1 - c1.nodeIndex = k; [u1,d1,c1] = hapod({S{k},u1},E,'incr_1',w,c1); end - c1.nodeIndex = n; - [u1,d1,c1] = hapod({S{n},u1},E,'incr_r',w,c1); + [U,D,C] = hapod({S{n},u1},E,'incr_r',w,c1); - iHAPOD_MEAN_L2 = norm(s-U*U'*s,'fro')/sqrt(n*n) - iHAPOD_MODES = size(U,2) - iHAPOD_MAX_MODES = max(C.nModes) - disp(''); + mean_L2 = norm(s-U*U'*s,'fro')/sqrt(N) + num_global_modes = size(U,2) + max_local_modes = max(cell2mat(C.nModes(1:end-1))) + disp(' '); %% Distributed HAPOD - % Compute global mode bound POD for incremental HAPOD - [ug,dg,cg] = hapod(S,E*w,'none'); - dHAPOD_GLOBAL_MODE_BOUND = size(ug,2) + % Test bulk distributed HAPOD + disp('DISTRIBUTED HAPOD (BULK):'); + [U,D,C] = hapod(S,E,'dist',w); + mean_L2 = norm(s-U*U'*s,'fro')/sqrt(N) + num_global_modes = size(U,2) + max_local_modes = max(cell2mat(C.nModes(1:end-1))) + disp(' '); - % Compute local mode bound POD for incremental HAPOD - [ul,dl,cl] = hapod(S,E*sqrt(1.0-w^2),'none'); - dHAPOD_LOCAL_MODE_BOUND = size(ul,2) - disp(''); + % Test chunk distributed HAPOD + disp('DISTRIBUTED HAPOD (CHUNK):'); + u2 = cell(1,n); - % Test full distributed HAPOD - [U,D,C] = hapod(S,E,'dist',w); - dHAPOD_MEAN_L2 = norm(s-U*U'*s,'fro')/sqrt(n*n) - dHAPOD_MODES = size(U,2) - dHAPOD_MAX_MODES = max(C.nModes) - disp(''); - - % Test partial distributed HAPOD - u2t = cell(1,n); - d2t = cell(1,n); - - nSnapshots = 0; - nModes = 0; - tNode = 0; - - for k=1:n - - [u2t{k},d2t{k},c2t] = hapod(S(k),E,'dist_1',w); - nSnapshots = nSnapshots + c2t.nSnapshots; - nModes = max(nModes,c2t.nModes); - tNode = tNode + c2t.tNode; + for k = 1:n + + [u2{k},d2,c2{k}] = hapod(S(k),E,'dist_1',w); + end + + [U,D,C] = hapod(u2,E,'dist_r',w,c2); + + mean_L2 = norm(s-U*U'*s,'fro')./sqrt(N) + num_global_modes = size(U,2) + max_local_modes = max(cell2mat(C.nModes(1:end-1))) + disp(' '); + +%% Distributed-of-Incremental HAPOD + + % Test dist-of-inc HAPOD + disp('DIST-OF-INC HAPOD (CHUNK):'); + strand = {[1:ceil(n/2)],ceil(n/2)+1:n}; + M = numel(strand); + u3 = repmat({[]},1,M); + c3 = repmat({struct('nLevels',n+1)},1,M); + + for m = 1:M + + for k = strand{m} + + [u3{m},d3,c3{m}] = hapod({S{k},u3{m}},E,'incr_1',w,c3{m}); + end + + c3{m}.nSnapshots = sum([c3{m}.nSnapshots{:}]); + c3{m}.nModes = max([c3{m}.nModes{:}]); + c3{m}.tNode = sum([c3{m}.tNode{:}]); end - [u2,d2,c2] = hapod(u2t,E,'dist_r',w,struct('nSnapshots',nSnapshots,'nModes',nModes,'tNode',tNode)); + [U,D,C] = hapod(u3,E,'dist_r',w,c3); - dHAPOD_MEAN_L2 = norm(s-u2*u2'*s,'fro')./sqrt(n*n) - dHAPOD_MODES = size(u2,2) - dHAPOD_MAX_MODES = max(c2.nModes) + mean_L2 = norm(s-U*U'*s,'fro')./sqrt(N) + num_global_modes = size(U,2) + max_local_modes = max(cell2mat(C.nModes(1:end-1))) + disp(' '); end diff --git a/hapod.m b/hapod.m index 003789d..814fdd2 100644 --- a/hapod.m +++ b/hapod.m @@ -1,214 +1,224 @@ -function [svec,sval,snfo] = hapod(data,bound,topo,relax,config,mysvd) -%%% project: hapod - Hierarchical Approximate POD ( http://git.io/hapod ) -%%% version: 1.3 ( 2018-02-09 ) +function [svec,sval,snfo] = hapod(data,bound,type,relax,config,mysvd) +%%% project: hapod - Hierarchical Approximate POD ( https://git.io/hapod ) +%%% version: 2.0 ( 2019-01-21 ) %%% authors: C. Himpe ( 0000-0003-2194-6754 ), S. Rave ( 0000-0003-0439-7212 ) %%% license: BSD 2-Clause License ( opensource.org/licenses/BSD-2-Clause ) %%% summary: Distributed or incremental POD / SVD computation % -%% SYNTAX: -% [svec,sval,snfo] = hapod(data,bound,topo,relax,config,mysvd) +% SYNTAX: +% [svec,sval,snfo] = hapod(data,bound,type,relax,config,mysvd) % -%% ABOUT: -% Compatible with OCTAVE and MATLAB. +% DESCRIPTION: +% The hierarhical approximate proper orthogonal decomposition is a tree-based +% algorithm to compute low-rank representations of column partitioned data +% matrices, of which the special cases of Incremental HAPOD and Distributed +% HAPOD are implemented. % -%% ARGUMENTS: -% (cell) data - snapshot data set, partitioned by column (blocks) -% (scalar) bound - mean L_2 projection error bound -% (string) topo - tree topology, default: 'none' -% * 'incr' - incremental HAPOD -% * 'incr_0' - incremental HAPOD -% * 'incr_1' - incremental HAPOD -% * 'incr_r' - incremental HAPOD -% * 'dist' - distributed HAPOD -% * 'dist_1' - distributed HAPOD -% * 'dist_r' - distributed HAPOD -% * 'none' - standard POD -% (scalar) relax - relaxation parameter in (0,1], default: 0.5 -% (struct) config - partial HAPOD configuration, default: [] -% * nSets - number of snapshot sets -% * nSnapshots - vector of data column count per node -% * nModes - vector of mode column count per node -% * nodeIndex - current nodes global index -% * tNode - vector of per node timings -% (handle) mysvd - custom SVD backend, default: [] -% * signature: [leftvec,singval] = mysvd(data) +% ABOUT: +% Compatible with OCTAVE and MATLAB. % -%% RETURNS: -% (matrix) svec - POD modes of the data -% (vector) sval - singular values to the data -% (struct) snfo - information structure (nSets, nSnapshots, nModes) +% ARGUMENTS: +% (cell) data - snapshot data set, partitioned by column (blocks) +% (scalar) bound - mean L_2 projection error bound +% (string) type - tree type, default: 'none' +% * 'incr' - incremental HAPOD +% * 'incr_1' - incremental HAPOD (non-root node) +% * 'incr_r' - incremental HAPOD (root node) +% * 'dist' - distributed HAPOD +% * 'dist_1' - distributed HAPOD (non-root node) +% * 'dist_r' - distributed HAPOD (root node) +% * 'none' - standard POD +% (scalar) relax - relaxation parameter in (0,1], default: 0.5 +% (struct) config - partial HAPOD configuration, default: {} +% * nLevels - total number of levels in tree +% * nSnapshots - cell array of data column count per node +% * nModes - cell array of mode column count per node +% * tNode - cell array of per node timings +% (handle) mysvd - custom SVD backend, default: [] +% * signature: [leftvec,singval] = mysvd(data) % -%% CITATION: -% C. Himpe, T. Leibner and S. Rave. -% "Hierarchical Approximate Proper Orthogonal Decomposition". -% Preprint, arXiv math.NA: 1607.05210, 2018. +% RETURNS: +% (matrix) svec - POD modes of the data +% (vector) sval - singular values to the data +% (struct) snfo - information structure % -%% SEE ALSO: -% svd, svds, princomp +% USAGE: +% If all data partitions can be passed as the data argument, the types: none +% (standard POD), incr(emental) HAPOD or dist(ributed) HAPOD are applicable. +% In case only a single partition can be passed, the types: incr_1 and dist_1 +% should be used for the non-root nodes of the associated HAPOD tree, while +% the types: incr_r and dist_r should be used for the root node. The returned +% information structure (or a cell-array thereof) can be passed to the parent +% nodes in the associated HAPOD tree. % -%% KEYWORDS: -% POD, BPOD, SVD, tSVD, PCA, MOR, Model Reduction +% CITATION: +% C. Himpe, T. Leibner and S. Rave. +% "Hierarchical Approximate Proper Orthogonal Decomposition". +% SIAM Journal on Scientific Computing, 40(5): A3267--A3292, 2018. % -% Further information: -%* - if(strcmp(data,'version')), svec = 1.3; return; end; +% SEE ALSO: +% svd, svds, princomp +% +% KEYWORDS: +% POD, BPOD, SVD, PCA, Model Reduction, Dimension Reduction +% +% Further information: https://git.io/hapod - if( (nargin<3) || isempty(topo) ), topo = 'none'; end; - if( (nargin<4) || isempty(relax) ), relax = 0.5; end; - if( nargin<6 ), mysvd = []; end; + if(strcmp(data,'version')), svec = 2.0; return; end%if - scaledBound = bound * sqrt(1.0 - relax^2); + if( (nargin<3) || isempty(type) ), type = 'none'; end%if + if( (nargin<4) || isempty(relax) ), relax = 0.5; end%if + if( nargin<5 ), config = []; end%if + if( nargin<6 ), mysvd = []; end%if - switch(lower(topo)) + % Precompute common quantities + scaledBound = bound * sqrt(1.0 - relax * relax); + relaxBound = bound * relax; - case 'incr' % Incremental HAPOD +%% Compute HAPOD - nSets = size(data,2); - nSnapshots = zeros(nSets,1); - nModes = zeros(nSets,1); - tNode = zeros(nSets,1); + switch(lower(type)) - leafBase = []; + case 'incr_1' % Incremental HAPOD (Non-Root Node) - for nodeIndex = 1:nSets-1 + [svec,sval,snfo] = incr_1(data,scaledBound,config,mysvd); - tId = tic(); - nSnapshots(nodeIndex) = size(data{nodeIndex},2); - leafBound = scaledBound * sqrt(sum(nSnapshots) / (nSets - 1)); - [leafBase,leafValues] = pod({leafBase,data{nodeIndex}},leafBound,mysvd); - leafBase = leafBase.*leafValues; - nModes(nodeIndex) = size(leafBase,2); - tNode(nodeIndex) = toc(tId); - end + case 'incr_r' % Incremental HAPOD (Root Node Only) - tId = tic(); - nSnapshots(end) = size(data{end},2); - rootBound = bound * relax * sqrt(sum(nSnapshots)); - [svec,sval] = pod({leafBase,data{end}},rootBound,mysvd); - tNode(nSets) = toc(tId); - snfo = struct('nSets',nSets,'nSnapshots',nSnapshots,'nModes',nModes,'tNode',tNode); + [svec,sval,snfo] = incr_r(data,relaxBound,config,mysvd); - case 'incr_0' % Incremental HAPOD (First Node Only) + case 'incr' % Incremental HAPOD (Complete) - config.nSnapshots = zeros(1,config.nSets); - config.nSnapshots(1) = size(data{1},2); - config.nModes = zeros(1,config.nSets); - config.tNode = zeros(1,config.nSets); + svec = []; + snfo.nLevels = size(data,2); - tId = tic(); - leafBound = scaledBound * sqrt(config.nSnapshots(1) / (config.nSets - 1)); - [svec,sval] = pod(data,leafBound,mysvd); - svec = svec.*sval; - config.nModes(1) = size(svec,2); - config.tNode(1) = toc(tId); - snfo = config; + for k = 1:size(data,2)-1 - case 'incr_1' % Incremental HAPOD (One Node Only) + [svec,sval,snfo] = incr_1({data{k},svec},scaledBound,snfo,mysvd); + end%for - tId = tic(); - config.nSnapshots(config.nodeIndex) = size(data{1},2) + size(data{2},2) - config.nModes(config.nodeIndex-1); - leafBound = scaledBound * sqrt(sum(config.nSnapshots) / (config.nSets - 1)); - [svec,sval] = pod(data,leafBound,mysvd); - svec = svec.*sval; - config.nModes(config.nodeIndex) = size(svec,2); - config.tNode(config.nodeIndex) = toc(tId); - snfo = config; + [svec,sval,snfo] = incr_r({data{end},svec},relaxBound,snfo,mysvd); - case 'incr_r' % Incremental HAPOD (Root Node Only) + case 'dist_1' % Distributed HAPOD (Non-Root Node) - tId = tic(); - config.nSnapshots(config.nodeIndex) = size(data{1},2) + size(data{2},2) - config.nModes(config.nodeIndex-1); - rootBound = bound * relax * sqrt(sum(config.nSnapshots)); - [svec,sval] = pod(data,rootBound,mysvd); - config.nModes(config.nodeIndex) = size(svec,2); - config.tNode(config.nodeIndex) = toc(tId); - snfo = config; + [svec,sval,snfo] = dist_1(data,scaledBound,config,mysvd); + + case 'dist_r' % Distributed HAPOD (Root Node Only) + + [svec,sval,snfo] = dist_r(data,relaxBound,config,mysvd); - case 'dist' % Distributed HAPOD + case 'dist' % Distributed HAPOD (Complete) - nSets = size(data,2); - nSnapshots = zeros(nSets,1); - nModes = zeros(nSets,1); - tNode = zeros(nSets,1); + for k = 1:size(data,2) - leafBases = cell(1,nSets); + [svec{k},sval,snfo{k}] = dist_1(data(k),scaledBound,config,mysvd); + end%for - for nodeIndex = 1:nSets + [svec,sval,snfo] = dist_r(svec,relaxBound,snfo,mysvd); - tId = tic(); - nSnapshots(nodeIndex) = size(data{nodeIndex},2); - leafBound = sqrt(nSnapshots(nodeIndex)) * scaledBound; - [leafBases{nodeIndex},leafValues] = pod(data(nodeIndex),leafBound,mysvd); - leafBases{nodeIndex} = leafBases{nodeIndex}.*leafValues; - nModes(nodeIndex) = size(leafBases{nodeIndex},2); - tNode(nodeIndex) = toc(tId); - end + otherwise % Standard POD tId = tic(); - rootBound = sqrt(sum(nSnapshots)) * relax * bound; - [svec,sval] = pod(leafBases,rootBound,mysvd); - tNode(nSets+1) = toc(tId); - snfo = struct('nSets',nSets,'nSnapshots',nSnapshots,'nModes',nModes,'tNode',tNode); + snfo.nSnapshots = sum(cellfun(@(M) size(M,2),data)); + snfo.nLevels = 1; + normBound = sqrt(snfo.nSnapshots) * bound; + [svec,sval] = pod(data,normBound,mysvd); + snfo.nModes = size(svec,2); + snfo.tNode = toc(tId); + end%switch +end - case 'dist_1' % Distributed HAPOD (One Node Only) +% LOCAL FUNCTION: incr_1 +function [svec,sval,config] = incr_1(data,scaledBound,config,mysvd) +% summary: Incremental HAPOD: Non-root node only - if(iscell(data)) - config.nSnapshots = size(data{1},2); - else - config.nSnapshots = size(data,2); - end + if(not(isfield(config,'nSnapshots'))) + config.nSnapshots{1} = size(data{1},2); + else + config.nSnapshots{end+1} = size(data{1},2) + size(data{2},2) - config.nModes{end}; + end%if - tId = tic(); - leafBound = sqrt(config.nSnapshots) * scaledBound; - [svec,sval] = pod(data,leafBound,mysvd); - svec = svec.*sval; - config.nModes = size(svec,2); - config.tNode = toc(tId); - snfo = config; + tId = tic(); + nodeBound = scaledBound * sqrt(sum(cell2mat(config.nSnapshots)) / (config.nLevels - 1)); + [svec,sval] = pod(data,nodeBound,mysvd); + svec = bsxfun(@times,svec,sval); - case 'dist_r' % Distributed HAPOD (Root Node Only) + if(not(isfield(config,'nModes'))) + config.nModes{1} = size(svec,2); + config.tNode{1} = toc(tId); + else + config.nModes{end+1} = size(svec,2); + config.tNode{end+1} = toc(tId); + end%if +end - tId = tic(); - rootBound = sqrt(sum(config.nSnapshots)) * relax * bound; - [svec,sval] = pod(data,rootBound,mysvd); - tNode = toc(tId); - snfo = struct('nSnapshots',[config.nSnapshots],'nModes',[config.nModes],'tNode',tNode + [config.tNode]); +% LOCAL FUNCTION: incr_r +function [svec,sval,config] = incr_r(data,relaxBound,config,mysvd) +% summary: Incremental HAPOD: Root node only - case 'none' % Standard POD + tId = tic(); + config.nSnapshots{end+1} = size(data{1},2) + size(data{2},2) - config.nModes{end}; + rootBound = sqrt(sum(cell2mat(config.nSnapshots))) * relaxBound; + [svec,sval] = pod(data,rootBound,mysvd); + config.nModes{end+1} = size(svec,2); + config.tNode{end+1} = toc(tId); +end - nSets = size(data,2); - nSnapshots = zeros(nSets,1); +% LOCAL FUNCTION: dist_1 +function [svec,sval,snfo] = dist_1(data,scaledBound,config,mysvd) +% summary: Distributed HAPOD: Non-root node only - tId = tic(); - for nodeIndex = 1:nSets - nSnapshots(nodeIndex) = size(data{nodeIndex},2); - end - [svec,sval] = pod(data,sqrt(sum(nSnapshots)) * bound,mysvd); - tNode = toc(tId); - nModes = size(svec,2); - snfo = struct('nSets',nSets,'nSnapshots',nSnapshots,'nModes',nModes,'tNode',tNode); - - otherwise - - error('hapod: unknown HAPOD topology!'); - end + if(iscell(data)) + snfo.nSnapshots = size(data{1},2); + else + snfo.nSnapshots = size(data,2); + end%if + + tId = tic(); + nodeBound = sqrt(snfo.nSnapshots) * scaledBound; + [svec,sval] = pod(data,nodeBound,mysvd); + svec = bsxfun(@times,svec,sval); + snfo.nLevels = 2; + snfo.nModes = size(svec,2); + snfo.tNode = toc(tId); end +% LOCAL FUNCTION: dist_r +function [svec,sval,snfo] = dist_r(data,relaxBound,config,mysvd) +% summary: Distributed HAPOD: Root node only + + tId = tic(); + snfo.nSnapshots = cellfun(@(c) {c.nSnapshots},config); + snfo.nSnapshots{end+1} = sum(cell2mat(snfo.nSnapshots)); + rootBound = sqrt(snfo.nSnapshots{end}) * relaxBound; + [svec,sval] = pod(data,rootBound,mysvd); + snfo.nLevels = 2; + snfo.nModes = [cellfun(@(c) {c.nModes},config),size(svec,2)]; + snfo.tNode = [cellfun(@(c) {c.tNode},config),toc(tId)]; +end + +% LOCAL FUNCTION: pod function [svec,sval] = pod(data,bound,mysvd) +% summary: Generic proper orthogonal decomposition X = cell2mat(data); - if(isempty(mysvd)) + if(isempty(mysvd)) % Plain economic (rank-revealing) SVD [U,D,V] = svd(X,'econ'); - else + D = diag(D); + elseif(strcmp(mysvd,'mos')) % Method of snapshots + [E,V] = eig(X'*X,'vector'); + [D,I] = sort(sqrt(abs(V)),'descend'); + U = X * bsxfun(@rdivide,E(:,I),D'); + else % Custom SVD [U,D,V] = mysvd(X); - end + if(size(D,2) > 1), D = diag(D); end%if + end%if - D = diag(D); d = flipud(cumsum(flipud(D.^2))); - K = find(d <= bound^2,1); - if(isempty(K)), K = size(X,2) + 1; end; + K = find(d <= (bound*bound),1); + if(isempty(K)), K = size(X,2) + 1; end%if sval = D(1:K-1)'; svec = U(:,1:K-1); end