Pde proj

.gitignore

@@ -0,0 +1,5 @@
+
+stenglib/Fast/fsetop.mexw64
+stenglib/Fast/fsparse.mexw64
+stenglib/Tensor/tprod.mexw64
+stenglib/Tensor/tsum.mexw64

startup.m

@@ -8,4 +8,9 @@
 
 % path to urdme/ folders
 addpath(genpath([link 'urdme/']));
-addpath(genpath([link 'workflows/']));
+addpath(genpath([link 'workflows/']));
+
+% Johannes Dufva 2020-10-30
+% path to other functions used
+run('stenglib/make_all.m')
+run('stenglib/startup.m')

stenglib/Fast/clenshaw.m

@@ -0,0 +1,43 @@
+function Y = clenshaw(A,B,C)
+%CLENSHAW Evaluation of 3-term recurrences.
+%   Y = CLENSHAW(A,B) evaluates the 3-term recurrence as defined by the
+%   matrices A and B.
+%
+%   Let [M,N] be MAX(SIZE(A),SIZE(B)). Then the output Y is M-by-N and, for
+%   1 <= j <= N, the recurrence is defined by Y(1,j) = A(1,j)+B(1,j), Y(2,j)
+%   = A(2,j)*Y(1,j)+B(2,j) and, for i > 2, Y(i,j) =
+%   A(i,j)Y(i-1,j)+B(i,j)Y(i-2,j).
+%
+%   The dimensions of the matrices A and B must either match or be
+%   singletons. Singleton dimensions are as usual silently repeated.
+%
+%   Y = CLENSHAW(A,B,C) evaluates the 3-term recurrence defined by A and B,
+%   multiplies it by the coefficients C and computes the sum of the result
+%   by means of Clenshaws summation formula. C is M-by-K-by-N and the result
+%   is K-by-N. The first dimension of C must be M while the second and third
+%   dimensions of C can be singletons.
+%
+%   Examples:
+%     % Fibonacci numbers
+%     n = 10; clenshaw(1,[0; 0; ones(n-2,1)])'
+%
+%     % Legendre polynomials
+%     x = linspace(-1,1); n = 4; j = [0.5 1:n]';
+%     A = (2*j-1)./j*x;
+%     B = (1-j)./j;
+%     y = clenshaw(A,B);
+%     figure, plot(x,y);
+%
+%     % two functions defined by Chebyshev coefficients
+%     x = linspace(-1,1,50);
+%     n = 4; j = [0; 1; 2*ones(n-1,1)];
+%     A = j*x; B = [1; 0; -ones(n-1,1)];
+%     C = [[1 0.75 0.5 0.25 0.125]' [1 -1.75 1.5 -1.25 2.125]'];
+%     figure, plot(x,clenshaw(A,B,C));
+%
+%   See also FILTER.
+
+% S. Engblom 2007-04-19 (Major revision)
+% S. Engblom 2005-07-28
+
+error('.MEX-file not found on path.');

stenglib/Fast/frepmat.m

@@ -0,0 +1,61 @@
+function B = frepmat(A,rep)
+%FREPMAT Fast replication of array.
+%   FREPMAT does the same job as the MATLAB-function REPMAT but runs
+%   faster since the main part of the code is written in C.
+%
+%   B = FREPMAT(A,[M N]) creates a matrix B consisting of M-by-N
+%   copies of A.
+%
+%   In general, B = FREPMAT(A,REP) repeats A in the first dimension
+%   REP(1) times, in the second dimension REP(2) times and so on for
+%   all of the sizes in REP.
+%
+%   Unlike REPMAT, missing sizes in REP are consistently defined to be
+%   1. This means that if m is a scalar, then FREPMAT(A,m) does not
+%   produce the same result as REPMAT(A,m) (which is defined to be
+%   REPMAT(A,[m m])). Instead, FREPMAT(A,m) produces the same result
+%   as would REPMAT(A,[m 1]). Also, the traditional syntax
+%   FREPMAT(A,M,N,...) is not supported.
+%
+%   Examples:
+%     a = frepmat(1,1000);            % same as ones(1000,1)
+%     b = frepmat(int8(1),[100 100]); % same as ones(100,100,'int8')
+%
+%     s = sprand(10,10,0.5);
+%     a1 = frepmat(s,[100 10]); [nnz(a1) nzmax(a1)]
+%     a2 = repmat(s,[100 10]);  [nnz(a2) nzmax(a2)] % overallocation
+%
+%   See also REPMAT.
+
+%   - For doubles (real or complex), FREPMAT is about 15-20% faster
+%   than REPMAT.
+%
+%   - For other numerical types (single precision floats, integers,
+%   characters and logicals), FREPMAT is about 50% faster than REPMAT.
+%
+%   - For sparse matrices (real, complex or logical), FREPMAT is about
+%   85-95% faster than REPMAT and allocates much less memory.
+%
+%   - For other data-types, such as cell-, structure- and
+%   function-arrays, FREPMAT is about 10% faster than REPMAT.
+
+% S. Engblom 2004-10-28
+
+if isnumeric(A) || ischar(A) || islogical(A)
+  % runs fast
+  B = mexfrepmat(A,rep);
+else
+  % general code using the traditional 'indexing' style
+  rep = rep(:);
+  lenrep = size(rep,1);
+  len = max(lenrep,tndims(A));
+  sizA = int32(tsize(A,1:len));
+  rep = [rep; ones(len-lenrep,1)];
+
+  % while building the indices the speed of mexfrepmat is exploited
+  ix = cell(1,len);
+  for i = 1:len
+    ix{i} = mexfrepmat([1:sizA(i)]',rep(i));
+  end
+  B = A(ix{:});
+end

stenglib/Fast/fsetop.m

@@ -0,0 +1,102 @@
+function [c,ia,ib] = fsetop(op,a,b)
+%FSETOP Fast set operations based on hashing.
+%   FSETOP(OP,A,B) performs the set operation OP on the columns of the
+%   arrays A and B. The main differences from the set operations
+%   implemented in MATLAB are:
+%     - FSETOP runs faster because it is based on hashing instead of on
+%     sorting.
+%
+%     - In FSETOP, the order of the elements in the input arrays are
+%     retained in the output array, in contrast to the corresponding
+%     MATLAB functions where the result is sorted. This means that all
+%     indices for new elements are strictly increasing (except for IB
+%     in the INTERSECT operation and for JA and JB in the UNION
+%     operation). Also, in case of duplications, the first element
+%     (lowest column number) is selected.
+%
+%     - FSETOP always uses the columns of its inputs as elements, like the
+%     'rows' switch in the MATLAB commands, but for columns instead of
+%     rows. In particular, column vectors are not special cases, they are
+%     treated as single elements.
+%
+%     - FSETOP supports the following data types: double, single, char,
+%     logicals, int8, uint8, int16, uint16, int32, uint32, int64, uint64
+%     and cell-arrays containing any of the above data types.
+%
+%   Cell-arrays are treated as vectors and all output indices refer to a
+%   linear ordering. The output C is always a row cell-vector and cells
+%   containing the same data but of different types or shapes are
+%   considered unequal.
+%
+%   C = FSETOP('check',A) computes hash-values for each column in A. The
+%   output C is a row-vector containing uint32-type integers. This is not
+%   a set-operation but is provided as a useful tool for computing
+%   check-sums for large data-sets in a uniform way.
+%
+%   [C,IA,IB] = FSETOP('intersect',A,B) returns the columns common to
+%   both A and B. C = A(:,IA) = B(:,IB).
+%
+%   [C,IA] = FSETOP('setdiff',A,B) returns the columns in A that are
+%   not in B. C = A(:,IA).
+%
+%   [C,IA,IB] = FSETOP('setxor',A,B) returns the columns that are not
+%   in the intersection of A and B. C = [A(:,IA) B(:,IB)].
+%
+%   [C,IA,IB,JA,JB] = FSETOP('union',A,B) returns the combined columns
+%   from A and B but with no repetitions. C = [A(:,IA) B(:,IB)] and A
+%   = C(:,JA), B = C(:,JB). Note that the outputs JA and JB are not
+%   produced by the MATLAB-function UNION.
+%
+%   [B,IA,IB] = FSETOP('unique',A) returns the same columns as in A
+%   but with no repetitions. B = A(:,IA) and A = B(:,IB).
+%
+%   [IA,IB] = FSETOP('ismember',A,B) returns a logical vector IA
+%   containing 1 where the columns of A are also columns of B and 0
+%   otherwise. IB contains the index in B of each column in A and zero if
+%   no such index exists. A(:,IA) = B(:,IB(IA)).
+%
+%   Note: two elements are considered equal if and only if their
+%   bitwise representations are identical. This can sometimes give
+%   unexpected results. For example, with doubles, -0.0 ~= 0.0 and NaN
+%   == NaN.
+%
+%   Examples:
+%     % intersection of integers
+%     a = int8(ceil(3*rand(2,10))), b = int8(ceil(3*rand(2,10)))
+%     aandb = fsetop('intersect',a,b)
+%
+%     % unique strings in cell-array
+%     strs = {'foobar' 'foo' 'foobar' 'bar' 'barfoo'}
+%     strsunq = fsetop('unique',strs)
+%
+%     % find indices ix in b = a(ix)
+%     a = randperm(6), b = ceil(6*rand(1,6))
+%     [foo,ix] = fsetop('ismember',b,a)
+%
+%     % remove indices from set of indices
+%     a = {[1 2 4] [2 3 5 6] [4 2 1] [2 3 5 6]' int32([1 2 4])};
+%     b = {[4 2 1]};
+%     c = fsetop('setdiff',a,b)
+%
+%     % union of structs
+%     a = struct('foo',1,'bar',NaN,'foobar','hello')
+%     b = struct('faa',2,'bar',Inf,'foobar','goodbye')
+%     [cf,ia,ib] = fsetop('union',fieldnames(a),fieldnames(b));
+%     ac = struct2cell(a); bc = struct2cell(b);
+%     c = cell2struct([ac(ia); bc(ib)],cf)
+%
+%     % a single 32-bit checksum from strings
+%     s = {'check','intersect','setdiff','setxor', ...
+%          'union','unique','ismember'};
+%     c = fsetop('check',fsetop('check',s)')
+%
+%   See also INTERSECT, SETDIFF, SETXOR, UNION, UNIQUE, ISMEMBER.
+
+% S. Engblom 2010-02-09 (Revision)
+% S. Engblom 2006-06-13 (Major revision)
+% S. Engblom 2005-06-21
+% Based on a concept by P-O Persson, COMSOL AB. The internal hash
+% function used is based on a code by P. Hsieh, see
+% http://www.azillionmonkeys.com/qed/hash.html.
+
+error('.MEX-file not found on path.');

stenglib/Fast/fsparse.m

@@ -0,0 +1,145 @@
+function S = fsparse(ii,jj,ss,siz,flag,nthreads)
+%FSPARSE Fast assembly of sparse matrix.
+%   FSPARSE does the same job as the MATLAB-function SPARSE but runs faster
+%   since the assembly is based on indirect addressing and on hashing rather
+%   than on sorting. FSPARSE also in general allocates the sparse matrix
+%   exactly, that is, the resulting sparse matrix S satisfies NNZ(S) =
+%   NZMAX(S). In addition, FSPARSE supports an efficient and powerful
+%   'assembly' syntax which makes it easy to create sparse banded matrices
+%   or sparse matrices arising in finite difference computations.
+%
+%   S = FSPARSE(X), where X is a full or sparse matrix, constructs a sparse
+%   matrix S by squeezing out any zero elements.
+%
+%   S = FSPARSE(II,JJ,SS) creates a sparse matrix S from triplet data
+%   (II,JJ,SS). The inputs are all matrices with either matching or
+%   singleton dimensions according to certain rules (see below). For
+%   instance, II may be 4-by-10, JJ 1-by-10 and SS 4-by-1.
+%
+%   The result is a sparse matrix formally satisfying the relation
+%   S(II(i,j),JJ(i,j)) = SS(i,j), except of course for singleton
+%   dimensions. In the example above for instance, we have the relation
+%   S(II(i,j),JJ(1,j)) = SS(i,1). Repeated indices are as usual summed
+%   together.
+%
+%   If II is IM-by-IN, JJ JM-by-JN and SS SM-by-SN, then it is required that
+%   (1) IN = JN or 1, (2) JM = IM or 1, (3) SM = IM or 1 and (4) SN = JN or
+%   1.
+%
+%   S = FSPARSE(II,JJ,SS,[M N Nzmax]) also specifies the dimensions of the
+%   output. All arguments are optional and when they are omitted M =
+%   MAX(II), N = MAX(JJ) and Nzmax = NNZ(S) are used as defaults.
+%
+%   S = FSPARSE(II,JJ,SS,SIZ,'nosort') produces a sparse matrix S
+%   where the columns are not sorted with respect to row-indices (if
+%   the input is sorted with respect to II, then the output is in
+%   order anyway). This runs faster and the resulting unordered matrix
+%   can be used in many, but not all, situations in
+%   MATLAB. Cautionary: on some platforms an unordered sparse matrix
+%   causes MATLAB to crash when performing certain operations.
+%
+%   One or both of the index-matrices II and JJ may be integers instead of
+%   doubles. This is faster still since no typecast or deep copy is
+%   needed. An 'integer' is the type which corresponds to the C declaration
+%   'int' on the platform; -- on most modern platforms this is int32.
+%
+%   Differences between FSPARSE and SPARSE:
+%     - FSPARSE generally allocates the sparse matrix so that it fits
+%     exactly. However, zero elements resulting from cancellation or from
+%     zero data causes a slight over-allocation. Thus, both FSPARSE([1 1],[1
+%     1],[1 -1]) and FSPARSE(1,1,0) allocates space for one entry.
+%
+%     - The call S = FSPARSE(X), where X is sparse, always makes a deep copy
+%     of the sparse matrix so that NNZ(S) = NZMAX(S).
+%
+%     - Integer indices are not supported by SPARSE.
+%
+%     - According to the specification above, the call FSPARSE(1:n,1,1) is
+%     not allowed. The corrected version of this example is simply
+%     FSPARSE([1:n]',1,1) where the row-index is associated with the first
+%     dimension of the input matrix.
+%
+%     - Logical sparse matrices cannot be created by FSPARSE. Use logical
+%     operations for this purpose.
+%
+%     - SPARSE allows the input (II,JJ,SS) itself to be sparse or even
+%     logical sparse. This syntax is not supported by FSPARSE.
+%
+%   Examples:
+%     N = 10;
+%     S1 = fsparse(1:N,1:N,1);    % same as speye(N)
+%     S2 = fsparse([1:N]',1:N,1); % same as sparse(ones(N));
+%
+%     % sparse matrix from triplet
+%     ii = ceil(4*rand(1,N));
+%     jj = ceil(4*rand(1,N));
+%     ss = rand(1,N);
+%     S3 = fsparse(ii,jj,ss); [nnz(S3) nzmax(S3)]
+%     s3 = sparse(ii,jj,ss);  [nnz(s3) nzmax(s3)] % over-allocation
+%
+%     % 3-point stencil
+%     S4 = fsparse([2:N-1]',[1:N-2; 2:N-1; 3:N]',[1 -2 1],[N N]);
+%
+%     % 5-point stencil on an N-by-N grid
+%     N2 = N*N; ix = reshape(1:N2,N,N);
+%     ix([1 N],:) = []; ix(:,[1 N]) = []; ix = ix(:);
+%     S5 = fsparse(ix,[ix-N ix-1 ix ix+1 ix+N],[1 1 -4 1 1],[N2 N2]);
+%
+%     % band-matrix
+%     B = rand(3,N); B(end,1) = 0; B(1,end) = 0;
+%     % same as spdiags(B',[-1 0 1],N,N):
+%     S6 = fsparse([[2:N 1]; [1:N]; [N 1:N-1]],1:N,B);
+%
+%     % circulant matrix
+%     jj = 1+mod(cumsum([0:3; ones(N-1,4)]),N);
+%     S7 = fsparse([1:N]',jj,1:4);
+%
+%   See also SPARSE, NNZ, NZMAX, SPDIAGS, FIND.
+
+
+%   Hidden argument nthreads and OpenMP-support
+%
+%   The full call currently supported is actually
+%     S = FSPARSE(II,JJ,SS,SIZ,[] | 'sort' | 'nosort',nthreads)
+%   with nthreads an integer >= 1. The default value is
+%   omp_get_max_threads(). The value of nthreads is remembered until
+%   the mex-file is cleared from memory:
+%     S1 = fsparse(...,4); % use 4 threads
+%     S2 = fsparse(...);   % continues to use 4 threads
+%     clear all;
+%     S3 = fsparse(...);   % uses omp_get_max_threads()
+%   The threaded version is a beta-release and is only supported for
+%   a subset of the syntaxes covered by FSPARSE. Also, compilation
+%   must be done under OpenMP, see MAKE.
+
+%   Hidden output timing
+%
+%   With #define FSPARSE_TIME an extra output time vector is supported
+%   as follows:
+%     [S,t] = FSPARSE(...)
+%   returns an 1-by-6 vector containing timings of the internal
+%   phases of FSPARSE.
+%     t(1) getix()-function
+%     t(2) Part 1
+%     t(3) Part 2
+%     t(4) Part 3
+%     t(5) Part 4
+%     t(6) sparse_insert()
+%   The suppport for timing is a beta-release and is only supported
+%   for a subset of the syntaxes covered by FSPARSE. See MAKE.
+
+%   Early performance observations for the serial version
+%
+%   In general, when the input indices and values are of the same
+%   sizes, FSPARSE is about 30-70% faster than SPARSE depending on the
+%   number of colliding indices, the usage of integer indices and the
+%   'nosort' option. When converting from full to sparse storage
+%   format, FSPARSE is about 50% faster than MATLAB. Finally, for the
+%   'assembly' syntax, FSPARSE is about 50-60% faster than any
+%   equivalent construction in MATLAB.
+
+% S. Engblom 2013-12-03 (Revision, OpenMP added, hidden argument nthreads).
+% S. Engblom 2010-01-13 (Minor revision, caution with 'nosort' added)
+% S. Engblom 2004-11-12
+
+error('.MEX-file not found on path.');

stenglib/Fast/powerseries.m

@@ -0,0 +1,25 @@
+y = powerseries(c,x,tol)
+%POWERSERIES Sum power series.
+%   Y = POWERSERIES(C,X,tol) evaluates the power series with coefficients C
+%   at points X where both C and X are vectors. The output Y is the same
+%   size as X and contains Y(i) = sum_(j >= 1) C(j)*X(i)^(j-1). The sum is
+%   truncated when the last contribution is relatively less than tol.
+%
+%   Note: use
+%     warning('off','powerseries:w1');
+%   to turn the warning for inaccurate evaluation off. As an
+%   alternative, you may explicitly set C(end) = 0.
+%
+%   Example:
+%     % the complex exponential function
+%     c = 1./cumprod([1 1:10]);
+%     x = complex(linspace(-2,3),linspace(-4,6));
+%     y = powerseries(c,x,0.5e-7);
+%     figure, plot(real(x),real(y),'b',real(x),real(exp(x)),'r');
+%     legend('Series','Exp()');
+
+% S. Engblom 2007-08-17 (Minor revision)
+% S. Engblom 2007-06-13 (Minor revision)
+% S. Engblom 2007-01-22
+
+error('.MEX-file not found on path.');