# HG changeset patch # User samer # Date 1358428844 0 # Node ID db7f4afd27c594dacd49e5874086f6a1f72c4c2e # Parent 19ec801ee487772baf6a7332484b0174332bc820 Rearranging numerical toolbox. diff -r 19ec801ee487 -r db7f4afd27c5 dsp/rot45.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/rot45.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,17 @@ +function B = rot45(A) +% rot45 - rotate square array by 45 degrees +% +% rot45 :: [[N,N]] -> [[N,N]]. +% +% Array is SUBSAMPLED, so, you should pre-filter to avoid aliasing. +% Result is padded with zeros in the corners. + +n = length(A); + +B = diag(A,0)'; +pad = [0]; + +for k=2:2:n-1 + B = [ pad diag(A,k)' pad; B; pad diag(A,-k)' pad]; + pad = [pad 0]; +end diff -r 19ec801ee487 -r db7f4afd27c5 general/general.m --- a/general/general.m Wed Jan 16 12:17:09 2013 +0000 +++ b/general/general.m Thu Jan 17 13:20:44 2013 +0000 @@ -22,7 +22,8 @@ end if strcmp(dirs,'all') - dirs={'algo','arrutils','cellutils','fileutils','funutils','numerical','discretise'}; + dirs={'algo','arrutils','cellutils','fileutils','funutils','discretise', ... + 'numerical','numerical/scalar','numerical/array','numerical/matrix'}; end thisfile=which('general'); seps=strfind(thisfile,filesep); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/argmax.m --- a/general/numerical/argmax.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function i=argmax(v) -% argmax - linear subscript of the largest element array. -% -% argmax :: [[Size]->A] -> 1..prod(Size) :- numeric(A). - -[m i] = max(v(:)); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/argmax1.m --- a/general/numerical/argmax1.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,15 +0,0 @@ -function [I,Y]=argmax1(D,X) -% argmax1 - return indices of maxima along a particular dimension -% -% argmax1 :: -% D:1..E ~'dimension over which to argmax', -% [S:[[E]->natural]->real] ~ -% 'array of size S where S is an array of E naturals. -% The array must be real so linear ordering is defined' -% -> [T:[[E]->natural]->1..S(D)] ~ -% 'array of indices of maxima, of size T, where T is like -% 'S except that T(D)=1. Indices are in the range 1..S(D)', -% [T:[[E]->natural]->real] ~'the actual maximum values'. - -[Y,I]=max(X,[],D); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/argmaxv.m --- a/general/numerical/argmaxv.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -function indices = argmaxv(v) -% argmaxv - multidimensional index of maximal element of array -% -% argmaxv :: [Size:[[1,E]]->A] -> [[1,E]] :- numeric(A). -% -% Returns the first maximum in the case of ties. - -[m i] = max(v(:)); -[ind{1:numdims(v)}]=ind2sub(size(v),i); -indices=[ind{:}]; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/argmin1.m --- a/general/numerical/argmin1.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,15 +0,0 @@ -function [I,Y]=argmin1(D,X) -% argmin1 - return index of minimum along a particular dimension -% -% argmin1 :: -% D:1..E ~'dimension over which to argmax', -% [S:[[E]->natural]->real] ~ -% 'array of size S where S is an array of E naturals. -% The array must be real so linear ordering is defined' -% -> [T:[[E]->natural]->1..S(D)] ~ -% 'array of indices of maxima, of size T, where T is like -% 'S except that T(D)=1. Indices are in the range 1..S(D)', -% [T:[[E]->natural]->real] ~'the actual minimum values'. - -[Y,I]=min(X,[],D); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/argmax.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/argmax.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,6 @@ +function i=argmax(v) +% argmax - linear subscript of the largest element array. +% +% argmax :: [[Size]->A] -> 1..prod(Size) :- numeric(A). + +[m i] = max(v(:)); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/argmax1.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/argmax1.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,15 @@ +function [I,Y]=argmax1(D,X) +% argmax1 - return indices of maxima along a particular dimension +% +% argmax1 :: +% D:1..E ~'dimension over which to argmax', +% [S:[[E]->natural]->real] ~ +% 'array of size S where S is an array of E naturals. +% The array must be real so linear ordering is defined' +% -> [T:[[E]->natural]->1..S(D)] ~ +% 'array of indices of maxima, of size T, where T is like +% 'S except that T(D)=1. Indices are in the range 1..S(D)', +% [T:[[E]->natural]->real] ~'the actual maximum values'. + +[Y,I]=max(X,[],D); + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/argmaxv.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/argmaxv.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,10 @@ +function indices = argmaxv(v) +% argmaxv - multidimensional index of maximal element of array +% +% argmaxv :: [Size:[[1,E]]->A] -> [[1,E]] :- numeric(A). +% +% Returns the first maximum in the case of ties. + +[m i] = max(v(:)); +[ind{1:numdims(v)}]=ind2sub(size(v),i); +indices=[ind{:}]; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/argmin1.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/argmin1.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,15 @@ +function [I,Y]=argmin1(D,X) +% argmin1 - return index of minimum along a particular dimension +% +% argmin1 :: +% D:1..E ~'dimension over which to argmax', +% [S:[[E]->natural]->real] ~ +% 'array of size S where S is an array of E naturals. +% The array must be real so linear ordering is defined' +% -> [T:[[E]->natural]->1..S(D)] ~ +% 'array of indices of maxima, of size T, where T is like +% 'S except that T(D)=1. Indices are in the range 1..S(D)', +% [T:[[E]->natural]->real] ~'the actual minimum values'. + +[Y,I]=min(X,[],D); + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/arrshift.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/arrshift.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,24 @@ +function Y=arrshift(O,X) +% arrshift - columnwise or rowwise subtraction for arrays +% +% This SUBTRACTS the first argument (a vector) from each +% of the vectors (row or column) in the second argument. +% Works in two modes: row vector or column vector mode. +% +% arrshift :: +% [Size] ~'values to subtract', +% [[N,M]] ~'array of vectors, array domain is M' +% -> [[N,M]] ~'vectors relative to new origin'. +% +% The first argument is REPMATed up to the size of +% the second and then subtracted. +% +% NB: this used to be called vecshift. vecshift is now +% specialised to subtract a COLUMN vector from an array +% the same size in the first dimension. This version +% retains the generality of the original vecshift. + +% note: this now works for any pair of arrays where size O +% is an integer fraction of size X in any dimension +Y=X-repmat_to(O,size(X)); + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/binary.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/binary.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,14 @@ +% binary - return array contain binary sequece over N columns +% +% binary :: N:natural -> [[2^N,N]->{0,1}]. +% +% note: returns an array of uint32 + +function B=binary(N) + +M=2^N; +B=zeros(M,N,'uint8'); +b=1:N; +for i=1:M + B(i,:)=bitget(uint32(i-1),b); +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/bitmap.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/bitmap.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,9 @@ +function X=bitmap(varargin), X=full(sbitmap(varargin{:})); +% bitmap - turn sequence of integers into binary bitmap +% +% bitmap :: +% Y:[[1,L]->[K]] ~'sequence of L natural numbers in 1..K', +% K:natural ~'height of array to return' +% -> X:[[K,L]->0|1] ~'X(i,j)=1 if Y(j)=i'. +% +% See also SBITMAP diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/complexdom.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/complexdom.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,4 @@ +function Z=complexdom(reals,imags) + +[X,Y]=meshgrid(reals,imags); +Z=X+i*Y; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/diffdims.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/diffdims.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,9 @@ +function X=diffdims(X,Dims) +% diffdims - diff along multiple dimensions +% +% diffdims :: +% [D:[[1,E]->natural]] ~'E dimensional array of size D', +% [[K]->[E]] ~'K dimensions each from 1 to E' +% -> [DD]. ~'array of diffs, smaller by one in each diffed dim'. + +for d=Dims, X=diff(X,1,d); end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/gauss_win.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/gauss_win.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,7 @@ +function w=gauss_win(n,sigma) +% gauss_win - Gaussian window of given length and std deviation +% +% gauss_win :: N:natural, nonneg~'std dev as multiple of N/2' -> [[N]]. + +s=(1:n)'-(n+1)/2; +w=expquadf(2*s,sigma*n); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/hgauss_win.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/hgauss_win.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,6 @@ +function w=hgauss_win(n,sigma) +% w=hgauss_win(n,sigma): half-Gaussian window +% half a gaussian window of length n +% sigma is std dev rel to window length + +w=expquadf((0:n)',sigma*n); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/incr.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/incr.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,9 @@ +function y=incr(y,varargin) +% incr - Increment element of array +% +% incr :: [[N]], 1..N-> [[N]]. +% incr :: [[M,N]], 1..M, 1..N -> [[M,N]]. +% incr :: [[L,M,N]], 1..L, 1..M, 1..N -> [[L,M,N]]. +% etc.. + + y(varargin{:})=y(varargin{:})+1; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/meandims.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/meandims.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,11 @@ +function A=meandims(A,dims) +% meandims - Sum over multiple dimensions +% +% meandims :: +% [Size:[[1,E]]->A] ~'E dimensional array', +% Dims:[[K]->1..E] ~'dims to sum over' +% -> [Size1->A] :- Size1=arrset(Size,Dims,1). +% +% Result is same size as input except for summed dimensions +% which are collapsed to a size of 1. +for i=dims, A=mean(A,i); end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/meanlastdims.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/meanlastdims.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,13 @@ +function x=meanlastdims(x,n) +% meanlastdims - compute average over last n dimensions of multidim array +% +% Usage: y=meanlastdims(x,n) +% +% If size1(x) = [m(1) m(2) ... m(k)], then size1(y) = [m(1) ... m(k-n)] +% MEANLASTDIMS ignores any trailing singleton dimensions including 2nd +% dimension of a column vector. + +sz=size1(x); +for d=length(sz)-(0:n-1); + x=sum(x,d)/sz(d); +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/minmax.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/minmax.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,15 @@ +function R=minmax(X,I) +% minmax - return minimum and maximum along a particular dimension +% +% minmax :: [[N,M]], 1:natural -> [[2,M]]. +% minmax :: [[N,M]], 2:natural -> [[N,2]]. +% minmax :: [D:[[1,E]]], I:[E] -> [set(D,I,2)]. +% +% The most general type means that the return array is the same size +% as the input array except that the size along the Ith dimension +% becomes 2, first element is min, second is max. +% +% The functions is constructed so that it is idemponent: +% minmax(X,I) == minmax(minmax(X,I),I) + +R= cat(I,min(X,[],I),max(X,[],I)); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/mmax.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/mmax.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,13 @@ +function [t,i,j]=mmax(A) +% mmax - maximum over both dimensions +% +% mmax :: [[N,M]-A] -> A, 1..N, 1..M. + +if nargout==1, t=max(A(:)); +else + [t1,i1]=max(A,[],1); + [t2,i2]=max(t1,[],2); + t=t2; + if nargout==3, i=i1(i2); j=i2; + elseif nargout==2, i=[i1(i2) i2]; end +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/mmin.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/mmin.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,12 @@ +function t=mmin(A) +% mmin - minimum over both dimensions +% +% mmin :: [[N,M]-A] -> A, 1..N, 1..M. +if nargout==1, t=min(A(:)); +else + [t1,i1]=min(A,[],1); + [t2,i2]=min(t1,[],2); + t=t2; + if nargout==3, i=i1(i2); j=i2; + elseif nargout==2, i=[i1(i2) i2]; end +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/mulrows.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/mulrows.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,6 @@ +% mulrows - Multiply rews of array by corresponding values in row vector +% +% mulrows :: [[N,M]], [[N]] -> [[N,M]]. +function y=mulrows(x,k) + y=repmat(k,1,size(x,2)).*x; +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/nary.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/nary.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,12 @@ +% nary - return array contain nary sequece over N columns +% +% nary :: M:natural, N:natural -> [[M^N,N]->1..M]. + +function B=nary(M,N) + +if (N==1), B=(1:M)'; +else + b=nary(M,N-1); + m=size(b,1); + B=[ kron((1:M)',ones(m,1)), repmat(b,M,1)]; +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/packvec.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/packvec.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,22 @@ +function X=packvec(varargin) +% packvec - Pack coordinate values in separate arrays into one big array +% +% packvec :: {[K]->[[E]]} -> [[K,E]]. +% +% There is also a variable argument list form: +% +% packvec :: [[E]], [[E]] -> [[2,E]]. +% packvec :: [[E]], [[E]], [[E]] -> [[3,E]]. +% etc.. + + +if nargin==1 && iscell(varargin{1}) + Y=varargin{1}; +else + Y=varargin; +end +[Y{:}]=promote(Y{:}); +Y=cellmap(@(y)reshape(y,[1 size(y)]),Y); +X=cat(1,Y{:}); + + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/prod1.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/prod1.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,8 @@ +function Y=prod1(X) +% prod1 - prod over dimension 1 and shiftdim 1 +% +% prod1 :: [[N M]] -> [[M]]. + + Z=prod(X,1); + S=size(Z); + Y=reshape(Z,[S(2:end) 1]); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/proddims.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/proddims.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,11 @@ +function A=proddims(A,dims) +% proddims - Product over multiple dimensions +% +% proddims :: +% [Size:[[1,E]]->A] ~'E dimensional array', +% Dims:[[K]->1..E] ~'dims to prod over' +% -> [Size1->A] :- Size1=arrset(Size,Dims,1). +% +% Result is same size as input except for summed dimensions +% which are collapsed to a size of 1. +for i=dims, A=prod(A,i); end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/prodn.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/prodn.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,5 @@ +function Y=prodn(X,n) +% prod1 - prod over dimension 1 and shiftdim 1 +% +% prod1 :: [[N M]] -> [[M]]. +Y=shiftdim(proddims(X,1:n),n); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/range2set.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/range2set.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,4 @@ +% range2set - 2 element array specifying integer range to array of integers +% +% range2set :: integer, integer -> [[1,N]->integer]. +function Y=range2set(I), Y=I(1):I(2); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/ratesched.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/ratesched.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,7 @@ +function rates=ratesched(a,b,c) +% rates=ratesched(a,b,c): learning rate schedule +% a: start +% b: time to half +% c: number of steps +rates=b*a./(b:(b+c-1)); + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/sbitmap.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/sbitmap.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,13 @@ +function X=sbitmap(c,K) +% sbitmap - turn sequence of integers into sparse binary bitmap +% +% sbitmap :: +% Y:[[1,L]->[K]] ~'sequence of L natural numbers in 1..K', +% K:natural ~'height of array to return' +% -> X:[[K,L]->0|1] ~'X(i,j)=1 if Y(j)=i'. +% +% See also BITMAP + +L=length(c); +if nargin<2, K=max(c); end +X=sparse(c,1:L,1,K,L); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/squeeze_prod.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/squeeze_prod.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,10 @@ +function A=squeeze_prod(A,dims) +% squeeze_prod - Squeeze out first N dimensions by multiplying +% +% squeeze_prod :: [Size], K:natural -> [Size1] :- Size1=Size(K+1:end). +% +% Defined as squeeze_prod(A,K)=shiftdim(proddims(A,1:K),K) +if ~isempty(dims) + for i=dims, A=prod(A,i); end + A=squeeze(A); +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/squeeze_sum.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/squeeze_sum.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,10 @@ +function A=squeeze_sum(A,dims) +% squeeze_sum - Squeeze out first N dimensions by summing +% +% squeeze_sum :: [Size], K:natural -> [Size1] :- Size1=Size(K+1:end). +% +% Defined as squeeze_sum(A,K)=shiftdim(sumdims(A,1:K),K) +if ~isempty(dims), + for i=dims, A=sum(A,i); end + A=squeeze(A); +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/sum1.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/sum1.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,8 @@ +function Y=sum1(X) +% sum1 - sum over dimension 1 and shiftdim 1 +% +% sum1 :: [[N M]] -> [[M]]. + + Z=sum(X,1); + S=size(Z); + Y=reshape(Z,[S(2:end) 1]); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/sumdims.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/sumdims.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,12 @@ +function A=sumdims(A,dims) +% sumdims - Sum over multiple dimensions +% +% sumdims :: +% [Size:[[1,E]]->A] ~'E dimensional array', +% Dims:[[K]->1..E] ~'dims to sum over' +% -> [Size1->A] :- Size1=arrset(Size,Dims,1). +% +% Result is same size as input except for summed dimensions +% which are collapsed to a size of 1. + +for i=dims, A=sum(A,i); end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/sumn.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/sumn.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,5 @@ +function Y=sumn(X,n) +% sumn - sum over first n dimension and shiftdim n +% +% sumn :: [[N:[[1,L]] M]], L:natural -> [[M]]. +Y=shiftdim(sumdims(X,1:n),n); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/array/vecshift.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/array/vecshift.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,14 @@ +function Y=vecshift(O,X) +% vecshift - refer an array of vectors to a new origin. +% +% This SUBTRACTS the first argument (a vector) from each +% of the vectors (row or column) in the second argument. +% Works in two modes: row vector or column vector mode. +% +% vecshift :: +% [[N,1]] ~'new origin', +% [[N,M]] ~'array of vectors, array domain is M' +% -> [[N,M]] ~'vectors relative to new origin'. +% +Y=X-repmat(O,1,size(X,2)); + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/arrshift.m --- a/general/numerical/arrshift.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -function Y=arrshift(O,X) -% arrshift - columnwise or rowwise subtraction for arrays -% -% This SUBTRACTS the first argument (a vector) from each -% of the vectors (row or column) in the second argument. -% Works in two modes: row vector or column vector mode. -% -% arrshift :: -% [Size] ~'values to subtract', -% [[N,M]] ~'array of vectors, array domain is M' -% -> [[N,M]] ~'vectors relative to new origin'. -% -% The first argument is REPMATed up to the size of -% the second and then subtracted. -% -% NB: this used to be called vecshift. vecshift is now -% specialised to subtract a COLUMN vector from an array -% the same size in the first dimension. This version -% retains the generality of the original vecshift. - -% note: this now works for any pair of arrays where size O -% is an integer fraction of size X in any dimension -Y=X-repmat_to(O,size(X)); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/asym.m --- a/general/numerical/asym.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -function B=asym(A), B=A-A'; -% asym - Skew-symmetric part of matrix (/2) -% -% asym :: [[N,N]] -> [[N,N]]. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/binary.m --- a/general/numerical/binary.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,14 +0,0 @@ -% binary - return array contain binary sequece over N columns -% -% binary :: N:natural -> [[2^N,N]->{0,1}]. -% -% note: returns an array of uint8 - -function B=binary(N) - -M=2^N; -B=zeros(M,N,'uint8'); -b=1:N; -for i=1:M - B(i,:)=bitget(uint32(i-1),b); -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/bitmap.m --- a/general/numerical/bitmap.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,9 +0,0 @@ -function X=bitmap(varargin), X=full(sbitmap(varargin{:})); -% bitmap - turn sequence of integers into binary bitmap -% -% bitmap :: -% Y:[[1,L]->[K]] ~'sequence of L natural numbers in 1..K', -% K:natural ~'height of array to return' -% -> X:[[K,L]->0|1] ~'X(i,j)=1 if Y(j)=i'. -% -% See also SBITMAP diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/cauchy_norm.m --- a/general/numerical/cauchy_norm.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,21 +0,0 @@ -function l=cauchy_norm(DF,eta,maxit) -% cauchy_norm - Quick and dirty function to fit Cauchy density to data -% -% cauchy_norm :: -% [[N,L]] ~'data', -% real ~'learning rate', -% natural ~'max iterations' -% -> [[1,L]] ~'the norms'. - -w=1./mean(abs(DF)); -for k=1:maxit - y=w*DF; - g=1-mean(y.*score(y)); - w=w.*exp(eta*g); - if all(abs(g))<0.002, break; end; -end -l=1./w; - -function g=score(x) -g=2*x./(1+x.^2); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/circulant.m --- a/general/numerical/circulant.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,15 +0,0 @@ -function t=circulant(c) -% circulant - Construct circulant matrix from first column -% -% circulant :: [[N]] -> [[N,N]]. - -m = length(c); -c=c(:); - -x = [c(2:m) ; c(:)]; % build vector of user data -cidx = (0:m-1)'; -ridx = m:-1:1; -t = cidx(:,ones(m,1)) + ridx(ones(m,1),:); % Toeplitz subscripts -t(:) = x(t); % actual data - - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/clamp.m --- a/general/numerical/clamp.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,8 +0,0 @@ -function x=clamp(a,b,x) -% clamp - Clamp values into given (closed) interval -% -% clamp :: real, real, [Size] -> [Size]. - -x(xb)=b; -%y=min(b,max(a,x)); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/complexdom.m --- a/general/numerical/complexdom.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -function Z=complexdom(reals,imags) - -[X,Y]=meshgrid(reals,imags); -Z=X+i*Y; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/deta.m --- a/general/numerical/deta.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,180 +0,0 @@ -function [f] = deta(z,k) -% deta - Calculates Dirichlet functions of the form -% -% f = sum((-1)^n/(k*n+1)^z) -% -% over the entire complex plane -% Z may be complex and any size -% Best accuracy for Abs(z) < 100 -% -% Usage: f = deta(z) -% or f = deta(z,k) -% -% where k determines which Dirichlet function to sum -% For Eta (Zeta, Lambda): k=1 -% For Betad: k=2 -% -% This function can use a LOT of memory when size(z) -% is large. Consider using the Memory and Pack commands. -% Also, consider breaking z up into smaller chunks. -% -% Requires a complex Gamma routine. -% Tested under version 5.3.1 -% -%see also: Zeta, Eta, Lambda, Betad -%see also: sym/zeta.m -%see also: mhelp zeta - -%Andrew Odlyzko has Riemann Zeta critical line zeros listed on: -%http://www.research.att.com/~amo/zeta_tables/index.html - -%Paul Godfrey -%pgodfrey@conexant.com -%March 24, 2001 - -if nargin==1 - k=1; -end -k=k(1); -if k<1 | k>2 - warning('Unknown function being calculated! Results valid only for Real(z)>0.5') -% k=1 --> Eta --> Zeta or Lambda --> Bern numbers -% k=2 --> Betad --> Euler numbers -end - -[sizz]=size(z); -z=z(:).'; % make z a row vector -rz=real(z); -iz=imag(z); -iszero=find(z==0); -iseven=find(z==(round(z/2)*2) & rz<0 & iz==0); -isodd= find(z==(round((z-1)/2)*2+1) & rz<0 & iz==0); - -r=1/2; % reflection point -L=find(rz< r); -if ~isempty(L) - zL=z(L); - z(L)=1-zL; -end - -%series coefficients were calculated using -% c(m)=sum from n=m to 64 of (binomial(n,m)/2^n) for m=0:64 - -%coefficients are symmetrical about the 0.5 value. Each pair sums to +-1 -%abs(coefficients) look like erfc(k*m)/2 due to binomial terms -%sum(cm) must = 0.5 = eta(0) = betad(0) -%worst case error occurs for z = 0.5 + i*large - -cm= [ .99999999999999999997; - -.99999999999999999821; - .99999999999999994183; - -.99999999999999875788; - .99999999999998040668; - -.99999999999975652196; - .99999999999751767484; - -.99999999997864739190; - .99999999984183784058; - -.99999999897537734890; - .99999999412319859549; - -.99999996986230482845; - .99999986068828287678; - -.99999941559419338151; - .99999776238757525623; - -.99999214148507363026; - .99997457616475604912; - -.99992394671207596228; - .99978893483826239739; - -.99945495809777621055; - .99868681159465798081; - -.99704078337369034566; - .99374872693175507536; - -.98759401271422391785; - .97682326283354439220; - -.95915923302922997013; - .93198380256105393618; - -.89273040299591077603; - .83945793215750220154; - -.77148960729470505477; - .68992761745934847866; - -.59784149990330073143; - .50000000000000000000; - -.40215850009669926857; - .31007238254065152134; - -.22851039270529494523; - .16054206784249779846; - -.10726959700408922397; - .68016197438946063823e-1; - -.40840766970770029873e-1; - .23176737166455607805e-1; - -.12405987285776082154e-1; - .62512730682449246388e-2; - -.29592166263096543401e-2; - .13131884053420191908e-2; - -.54504190222378945440e-3; - .21106516173760261250e-3; - -.76053287924037718971e-4; - .25423835243950883896e-4; - -.78585149263697370338e-5; - .22376124247437700378e-5; - -.58440580661848562719e-6; - .13931171712321674741e-6; - -.30137695171547022183e-7; - .58768014045093054654e-8; - -.10246226511017621219e-8; - .15816215942184366772e-9; - -.21352608103961806529e-10; - .24823251635643084345e-11; - -.24347803504257137241e-12; - .19593322190397666205e-13; - -.12421162189080181548e-14; - .58167446553847312884e-16; - -.17889335846010823161e-17; - .27105054312137610850e-19]; - -cm=flipud(cm).'; % sum from small to big -nmax=size(cm,2); -n=(1:k:k*nmax)'; -n=flipud(n); -% z is a LR vector -% n is an UD vector -[Z,N]=meshgrid(z,n); - -% this can take a LOT of memory -f=cm*(N.^-Z); -% but it's really fast to form the series expansion N.^-Z -% and then sum it by an inner product cm*() :) - -%reflect across 1/2 -if ~isempty(L) - zz=z(L); - if k==1 - % Eta function reflection - % for test: deta(1,1) should = log(2) - t=(2-2.^(zz+1))./(2.^zz-2)./pi.^zz; - f(L)=t.*cos(pi/2*zz).*gamma(zz).*f(L); - if ~isempty(iszero) - f(iszero)=0.5; - end - if ~isempty(iseven) - f(iseven)=0; - end - end - if k==2 - % Betad function reflection - %for test: deta(0,2) should = 0.5 - %for test: deta(1,2) should = pi/4 - f(L)=(2/pi).^zz.*sin(pi/2*zz).*gamma(zz).*f(L); - if ~isempty(isodd) - f(isodd)=0; - end - end - if k>2 - % Insert reflection formula for other Dirichlet functions here - f(L)=(1/pi).^zz.*gamma(zz).*f(L); - f(L)=NaN; - end -end - -f = reshape(f,sizz); -return - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/diffdims.m --- a/general/numerical/diffdims.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,9 +0,0 @@ -function X=diffdims(X,Dims) -% diffdims - diff along multiple dimensions -% -% diffdims :: -% [D:[[1,E]->natural]] ~'E dimensional array of size D', -% [[K]->[E]] ~'K dimensions each from 1 to E' -% -> [DD]. ~'array of diffs, smaller by one in each diffed dim'. - -for d=Dims, X=diff(X,1,d); end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/eigsel.m --- a/general/numerical/eigsel.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,16 +0,0 @@ -function [r1,r2]=eigs(T,I) -% eigs - Return selected eigs relative to order by magnitude -% -% eigs :: [[N,N]], [[M]->[N]] ~'M indices between 1 and N' -> [[M]]. -% eigs :: [[N,N]], [[M]->[N]] ~'M indices between 1 and N' -> [[N,M]], [[M]]. - -if nargout==1 - L0=eig(T); - [dummy,ord]=sort(-abs(L0)); - r1=L0(ord(I)); -else - [V0,D0]=eig(T); L0=diag(D0); - [dummy,ord]=sort(-abs(L0)); - r1=V0(:,ord(I)); - r2=L0(ord(I)); -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/eigvecs.m --- a/general/numerical/eigvecs.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,5 +0,0 @@ -function V=eigvecs(A,I) - [V,D]=eig(A); - V=fliplr(V); - V=V(:,I); -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/erfp.m --- a/general/numerical/erfp.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,9 +0,0 @@ -function y=erfp(x) -% erfp - returns Pr(X>x) for a normalised Gaussian random variable -% -% erfp :: X:real -> Z:real. -% -% Y = integeral from X to infinity of -% (1/sqrt(2*pi)) * exp(-0.5*t.^2) - -y = erfc(x/sqrt(2))/2; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/erfratio.m --- a/general/numerical/erfratio.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,5 +0,0 @@ -function y=erfratio(x) -% erfratio - returns gausspdf(x)./erfp(x) but using a better computation -% -% erfratio :: real -> nonneg. -y = 2./(sqrt(2*pi)*erfcx(x/sqrt(2))); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/expdecay.m --- a/general/numerical/expdecay.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,9 +0,0 @@ -function lambda=expdecay(k,d) -% expdecay - compute interpolation ratio for exponential decay -% -% expdecay :: [Size->nonneg], [Size->nonneg] -> [Size->nonneg]. -% -% this gives exponential convergence - -if nargin==2, lambda=1-k; -else lambda=@(d)1-k; end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/fromdbs.m --- a/general/numerical/fromdbs.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -function X=fromdbs(Y) -% fromdbs - Convert from decibels -% -% fromdbs :: real -> real. - -X=10.^(Y/10); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/fuzzeye.m --- a/general/numerical/fuzzeye.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,30 +0,0 @@ -function I=fuzzeye(n,m,sigma) -% fuzzeye Fuzzy nonsquare unit matrix -% -% I=fuzzeye(n,m) -% I=fuzzeye([n m]) -% I=fuzzeye(n,m,sigma) -% I=fuzzeye([n m],sigma) -% -% return n by m 'fuzzy' unit matrix -% spread factor sigma is roughly in elements -% default sigma is m/n or n/m (the larger) - -if length(n)>1, - if nargin<2, - sigma=2*max([n(2)/n(1) n(1)/n(2)]); - else - sigma=m; - end - m=n(2); - n=n(1); -else - if nargin<3, sigma=max([m/n n/m]); end -end - -sigma=sigma/max([n m]); -[X,Y]=meshgrid(0:m-1,0:n-1); -X = X/(m-1); -Y = flipud(Y)/(n-1); -Z = ((X-Y).^2)/(sigma^2); -I = exp(-Z/2); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/gauss_fn.m --- a/general/numerical/gauss_fn.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -function x=gauss_fn(t), x=exp(-0.5*t.^2); -% gauss_fn - Gaussian kernel function, exp(-t^2/2) -% -% gauss_fn :: [D] ~'any size array' -> [D->noneg]~'same size result'. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/gauss_logpdf.m --- a/general/numerical/gauss_logpdf.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,9 +0,0 @@ -function w=gauss_logpdf(x,sigma) -% gausspdf - Gaussian probability density -% -% gausspdf :: [Size->real] -> [Size->nonneg]. -% gausspdf :: [Size->real], nonneg ~'stddev of Guassian' -> [Size->nonneg]. -% returns normal (gaussian) pdf of 1st argument. - -if nargin<2, sigma=1; end; -w=-0.5*(x./sigma).^2 - (log(sigma)+2*log(pi)); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/gauss_win.m --- a/general/numerical/gauss_win.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -function w=gauss_win(n,sigma) -% gauss_win - Gaussian window of given length and std deviation -% -% gauss_win :: N:natural, nonneg~'std dev as multiple of N/2' -> [[N]]. - -s=(1:n)'-(n+1)/2; -w=gauss_fn(2*s/(sigma*n)); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/gausspdf.m --- a/general/numerical/gausspdf.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,9 +0,0 @@ -function w=gausspdf(x,sigma) -% gausspdf - Gaussian probability density -% -% gausspdf :: [Size->real] -> [Size->nonneg]. -% gausspdf :: [Size->real], nonneg ~'stddev of Guassian' -> [Size->nonneg]. -% returns normal (gaussian) pdf of 1st argument. - -if nargin<2, sigma=1; end; -w=exp(-0.5*(x./sigma).^2)./(sigma*sqrt(2*pi)); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/hgauss_win.m --- a/general/numerical/hgauss_win.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function w=hgauss_win(n,sigma) -% w=hgauss_win(n,sigma): half-Gaussian window -% half a gaussian window of length n -% sigma is std dev rel to window length - -w=gauss_fn((0:n)'/(sigma*n)); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/hshrink.m --- a/general/numerical/hshrink.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -% shrink - Shrinkage towards zero -% -% shrink :: real, [Size] -> [Size]. -function y=shrink(th,x) - -y=x; -y(y nonneg. - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/hsphere_surfln.m --- a/general/numerical/hsphere_surfln.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -function S=hsphere_surfln(N), S=log(2)+(N/2)*log(pi)-gammaln(N/2); end -% hsphere_surfln - log of surface area of unit hypersphere -% -% hsphere_surfln :: natural -> real. -% -% hsphere_surfln(N)=log(hsphere_surf(N)) - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/hypdecay.m --- a/general/numerical/hypdecay.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,17 +0,0 @@ -function ret=hypdecay(k,d) -% hypdecay - compute interpolation ratio for hyperbolic decay -% -% hypdecay :: [Size->nonneg], [Size->nonneg] -> [Size->nonneg]. -% hypdecay :: [Size->nonneg] -> ([Size->nonneg] -> [Size->nonneg]). -% -% this gives 'hyperbolic' convergence which is more like -% a sort of diffusion by Brownian motion. The trick is -% to add a constant to the inverse of each natural parameter, -% which is like a temperature or variance. The constant -% is like a diffusion constant -% -% This function supports partial application: if only one -% argument is supplied, it returns a function handle. - -if nargin==2, ret=1./(1+k*d); -else ret=@(d)1./(1+k*d); end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/incr.m --- a/general/numerical/incr.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,9 +0,0 @@ -function y=incr(y,varargin) -% incr - Increment element of array -% -% incr :: [[N]], 1..N-> [[N]]. -% incr :: [[M,N]], 1..M, 1..N -> [[M,N]]. -% incr :: [[L,M,N]], 1..L, 1..M, 1..N -> [[L,M,N]]. -% etc.. - - y(varargin{:})=y(varargin{:})+1; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/inf2nan.m --- a/general/numerical/inf2nan.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function x=inf2nan(x) -% inf2nan - convert infs to nans -% -% inf2nan :: [Size] -> [Size]. - -x(isinf(x))=nan; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/inv_triu.m --- a/general/numerical/inv_triu.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,19 +0,0 @@ -function t=inv_triu(t) -% inv_triu - Inverse of upper triangular matrix -% -% inv_triu :: [[N,N]] -> [[N,N]]. - -n=size(t,1); -for k=1:n - if t(k,k)~=0 - t(k,k) = 1/t(k,k); - t(1:k-1,k) = -t(1:k-1,k) * t(k,k); - - j = k+1:n; - temp = t(k,j); - t(k,j) = 0; - t(1:k,j) = t(1:k,j) + t(1:k,k)*temp; - end -end - - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/iqform.m --- a/general/numerical/iqform.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -function z=iqform(A,u), z=u'*(A\u); -% iqform - Inverse quadratic form: iqform(A,x)=x'*inv(A)*x -% -% iqform :: [[N,N]], [[N,T]] -> [[T,T]]. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/iqforma.m --- a/general/numerical/iqforma.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -function z=iqforma(A,u), z=sum(u.*(A\u),1); -% iqforma - Quadratic form (array version) iqforma(A,x)=diag(x'*inv(A)*x) -% -% iqforma :: [[N,N]], [[N,T]] -> [[1,T]]. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/iseven.m --- a/general/numerical/iseven.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,2 +0,0 @@ -% iseven :: integer -> bool. -function b=iseven(n), b=~mod(n,2); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/isodd.m --- a/general/numerical/isodd.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,2 +0,0 @@ -% isodd :: integer -> bool. -function b=isodd(n), b=mod(n,2); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/issym.m --- a/general/numerical/issym.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,11 +0,0 @@ -function z=issym(x,th) -% issym - Test if matrix is symmetric with optional tolerance -% -% issym :: [[N,N]] -> bool. -% issym :: [[N,N]], nonneg -> bool. - -if nargin<2, - z=all(all(x==x')); -else - z=all(all(abs(x-x')<=th)); -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/lambertw.m --- a/general/numerical/lambertw.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -function W = lambertw(x) -% lambertw - Lambert's W function -% -% lambertw :: real -> real. -% -% evaluates W(y) for -1 branch of Lambert W function where -% y = -exp(-x) -% -% The W function satisfies W(t)*exp(W(t)) = t - -if x<1, error('lambertw: argument must be >= 1'); end -W=fixpoint(@(w)(-x-log(abs(w))),-x,'its',60); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/linterp.m --- a/general/numerical/linterp.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -function z=linterp(k,z0,z) -% linterp - linear interpolation between values -% -% linterp :: -% [Size] ~'fraction of distance to travel from from initial to final point', -% [Size] ~'initial point', -% [Size] ~'final point' -% -> [Size] ~'interpolated point'. - -z = z0+k.*(z-z0); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/linterp_d.m --- a/general/numerical/linterp_d.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,14 +0,0 @@ -function z=linterp_d(k,z0,d) -% linterp_d - linear interpolation between two points -% -% linterp_d :: -% real ~'fraction of distance to travel from from initial to final point', -% [Size] ~'initial point', -% [Size] ~'vector from initial to final point' -% -> [Size] ~'interpolated point'. -% -% Definition is linterp_d(k,z0,d) = linterp(k,z0,z0+d) -% -% See also: LINTERP, VLINTERP - -z = z0+k.*d; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/log2pie.m --- a/general/numerical/log2pie.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,5 +0,0 @@ -function k=log2pie - persistent K; - if isempty(K) K=1+log(2*pi); end - k=K; -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/logabsdet.m --- a/general/numerical/logabsdet.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,22 +0,0 @@ -function y = logabsdet(A) -% logabsdet - logabsdet(X)=log(abs(det(X))) -% -% logabsdet :: [[N,N]] -> nonneg. -% -% This is faster and more stable than using log(det(A)). -% (Except when A is not positive definite, in which case -% we fall back to computing det(A) the usual way.) - -% From Tom Minka's lightspeed toolbox -% Samer: added clause incase A not pos def. - -% this only works for HERMITIAN POS DEF matrix. -% y = 2*sum(log(diag(chol(A)))); - -[R,p]=chol(A*A'); -if p>0, - fprintf('*** logabsdet warning: matrix is singular'); - y=-inf; -else - y = sum(log(full(diag(R)))); -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/logdet.m --- a/general/numerical/logdet.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,14 +0,0 @@ -function y = logdet(A) -% logdet - logdet(X)=log(det(X)) where A is positive-definite and Hermitian. -% -% logdet :: [[N,N]] -> nonneg. -% -% This is faster and more stable than using log(det(A)). -% Samer: Use LOGABSDET for general matrices. - -% From Tom Minka's lightspeed toolbox - -[U,p] = chol(A); -if p>0, y=-inf; -else y = 2*sum(log(full(diag(U)))); -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/logeps.m --- a/general/numerical/logeps.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -function y=logeps(x), y=logth(eps,x); -% logeps - Natural log with lower threshold of eps -% -% logeps :: real -> real. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/logth.m --- a/general/numerical/logth.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function y=logth(eps,x) -% logth - Natural log with lower threshold -% -% logth :: real~'lower threshold', real -> real. - -x(x [[N,L]]. -% -% This is like matrix multiplication except that the N and L -% can stand for multiple dimensions, ie inputs can be more -% than 2-dimensional arrays. - -Adom=size(A); -Bdom=size(B); - -A=reshape(A,[],Bdom(1)); -B=reshape(B,Bdom(1),[]); -C=reshape(A*B,[Adom(1:end-1) Bdom(2:end)]); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/asym.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/asym.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,4 @@ +function B=asym(A), B=A-A'; +% asym - Skew-symmetric part of matrix (/2) +% +% asym :: [[N,N]] -> [[N,N]]. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/circulant.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/circulant.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,15 @@ +function t=circulant(c) +% circulant - Construct circulant matrix from first column +% +% circulant :: [[N]] -> [[N,N]]. + +m = length(c); +c=c(:); + +x = [c(2:m) ; c(:)]; % build vector of user data +cidx = (0:m-1)'; +ridx = m:-1:1; +t = cidx(:,ones(m,1)) + ridx(ones(m,1),:); % Toeplitz subscripts +t(:) = x(t); % actual data + + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/eigsel.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/eigsel.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,16 @@ +function [r1,r2]=eigs(T,I) +% eigs - Return selected eigs relative to order by magnitude +% +% eigs :: [[N,N]], [[M]->[N]] ~'M indices between 1 and N' -> [[M]]. +% eigs :: [[N,N]], [[M]->[N]] ~'M indices between 1 and N' -> [[N,M]], [[M]]. + +if nargout==1 + L0=eig(T); + [dummy,ord]=sort(-abs(L0)); + r1=L0(ord(I)); +else + [V0,D0]=eig(T); L0=diag(D0); + [dummy,ord]=sort(-abs(L0)); + r1=V0(:,ord(I)); + r2=L0(ord(I)); +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/eigvecs.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/eigvecs.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,8 @@ +% eigvec - Return selected eigenvectors of a matrix +% +% eigvec :: [[N,N]], [[M]->[N]] -> [[N,M]]. +function V=eigvecs(A,I) + [V,D]=eig(A); + V=fliplr(V); + V=V(:,I); +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/fuzzeye.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/fuzzeye.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,30 @@ +function I=fuzzeye(n,m,sigma) +% fuzzeye Fuzzy nonsquare unit matrix +% +% I=fuzzeye(n,m) +% I=fuzzeye([n m]) +% I=fuzzeye(n,m,sigma) +% I=fuzzeye([n m],sigma) +% +% return n by m 'fuzzy' unit matrix +% spread factor sigma is roughly in elements +% default sigma is m/n or n/m (the larger) + +if length(n)>1, + if nargin<2, + sigma=2*max([n(2)/n(1) n(1)/n(2)]); + else + sigma=m; + end + m=n(2); + n=n(1); +else + if nargin<3, sigma=max([m/n n/m]); end +end + +sigma=sigma/max([n m]); +[X,Y]=meshgrid(0:m-1,0:n-1); +X = X/(m-1); +Y = flipud(Y)/(n-1); +Z = ((X-Y).^2)/(sigma^2); +I = exp(-Z/2); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/inv_triu.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/inv_triu.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,19 @@ +function t=inv_triu(t) +% inv_triu - Inverse of upper triangular matrix +% +% inv_triu :: [[N,N]] -> [[N,N]]. + +n=size(t,1); +for k=1:n + if t(k,k)~=0 + t(k,k) = 1/t(k,k); + t(1:k-1,k) = -t(1:k-1,k) * t(k,k); + + j = k+1:n; + temp = t(k,j); + t(k,j) = 0; + t(1:k,j) = t(1:k,j) + t(1:k,k)*temp; + end +end + + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/iqform.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/iqform.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,4 @@ +function z=iqform(A,u), z=u'*(A\u); +% iqform - Inverse quadratic form: iqform(A,x)=x'*inv(A)*x +% +% iqform :: [[N,N]], [[N,T]] -> [[T,T]]. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/iqforma.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/iqforma.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,4 @@ +function z=iqforma(A,u), z=sum(u.*(A\u),1); +% iqforma - Quadratic form (array version) iqforma(A,x)=diag(x'*inv(A)*x) +% +% iqforma :: [[N,N]], [[N,T]] -> [[1,T]]. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/issym.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/issym.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,11 @@ +function z=issym(x,th) +% issym - Test if matrix is symmetric with optional tolerance +% +% issym :: [[N,N]] -> bool. +% issym :: [[N,N]], nonneg -> bool. + +if nargin<2, + z=all(all(x==x')); +else + z=all(all(abs(x-x')<=th)); +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/logabsdet.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/logabsdet.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,22 @@ +function y = logabsdet(A) +% logabsdet - logabsdet(X)=log(abs(det(X))) +% +% logabsdet :: [[N,N]] -> nonneg. +% +% This is faster and more stable than using log(det(A)). +% (Except when A is not positive definite, in which case +% we fall back to computing det(A) the usual way.) + +% From Tom Minka's lightspeed toolbox +% Samer: added clause incase A not pos def. + +% this only works for HERMITIAN POS DEF matrix. +% y = 2*sum(log(diag(chol(A)))); + +[R,p]=chol(A*A'); +if p>0, + fprintf('*** logabsdet warning: matrix is singular'); + y=-inf; +else + y = sum(log(full(diag(R)))); +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/logdet.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/logdet.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,14 @@ +function y = logdet(A) +% logdet - logdet(X)=log(det(X)) where A is positive-definite and Hermitian. +% +% logdet :: [[N,N]] -> nonneg. +% +% This is faster and more stable than using log(det(A)). +% Samer: Use LOGABSDET for general matrices. + +% From Tom Minka's lightspeed toolbox + +[U,p] = chol(A); +if p>0, y=-inf; +else y = 2*sum(log(full(diag(U)))); +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/minner.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/minner.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,13 @@ +function C=minner(A,B) +% minner - inner product for arbitrary dimension arrays +% +% ttimes :: [[I,J]], [[J,K]] -> [[I,K]]. +% +% ie just like mtimes, but I and K can be row vectors + +Adom=size1(A); +Bdom=size1(B); +A=reshape(A,[],Bdom(1)); +B=reshape(B,Bdom(1),[]), +C=reshape(A*B,tosize([Adom(1:end-1) Bdom(2:end)])); + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/mouter.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/mouter.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,14 @@ +function B=mouter(varargin) +% mouter - Multidimensional outer product of multiple arrays +% Index domain of result is the concatenation of the index domains +% of the arguments, (with trailing 1s removed). +% +% mouter :: [Size1->A], [Size2->A] -> [[Size1,Size2]->A]. +% mouter :: [Size1->A], [Size2->A], [Size3->A] -> [[Size1,Size2,Size3]->A]. +% etc. +B=1; +for i=1:length(varargin) + B=kron(varargin{i},B); +end +B=reshape(B,cell2mat(cellmap(@size1,varargin))); + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/orthogonalise.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/orthogonalise.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,13 @@ +function B1=orthogonalise(B), +% orthogonalise - Orthogonalise a basis matrix +% +% orthogonalise :: [[N,M]] -> [[N,M]]. +% +% Works using SVD. + +[U,S,V] = svd(B,0); +B1 = U*V'; + +% alternative method, seems to be slower +% B1 = B*real((B'*B)^(-0.5)); + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/outer.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/outer.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,4 @@ +% outer - outer product for 2D matrix +% +% outer :: [[N,M]] -> [[N,N]]. +function Y=outer(X), Y=X*X'; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/qform.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/qform.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,4 @@ +function z=qform(A,u), z=u'*A*u; +% qform - Quadratic form: qform(A,x)=x'*A*x +% +% qform :: [[N,N]], [[N,T]] -> [[T,T]]. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/qforma.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/qforma.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,4 @@ +function z=qforma(A,u), z=sum(u.*(A*u)); +% qforma - Quadratic form (array version) qforma(A,x)=diag(x'*A*x) +% +% qforma :: [[N,N]], [[N,T]] -> [[1,T]]. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/spdiag.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/spdiag.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,6 @@ +function A=spdiag(H) +% spdiag - Return elementwise multiplier as sparse diagonal matrix +% +% spdiag :: [[N]] -> [[N,N]]. + +A=diag(sparse(H)); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/matrix/sym.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/matrix/sym.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,4 @@ +function B=sym(A), B=A+A'; +% sym - Symmetric part of matrix (/2) +% +% sym :: [[N,N]] -> [[N,N]]. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/meandims.m --- a/general/numerical/meandims.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,11 +0,0 @@ -function A=meandims(A,dims) -% meandims - Sum over multiple dimensions -% -% meandims :: -% [Size:[[1,E]]->A] ~'E dimensional array', -% Dims:[[K]->1..E] ~'dims to sum over' -% -> [Size1->A] :- Size1=arrset(Size,Dims,1). -% -% Result is same size as input except for summed dimensions -% which are collapsed to a size of 1. -for i=dims, A=mean(A,i); end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/meanlastdims.m --- a/general/numerical/meanlastdims.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,13 +0,0 @@ -function x=meanlastdims(x,n) -% meanlastdims - compute average over last n dimensions of multidim array -% -% Usage: y=meanlastdims(x,n) -% -% If size1(x) = [m(1) m(2) ... m(k)], then size1(y) = [m(1) ... m(k-n)] -% MEANLASTDIMS ignores any trailing singleton dimensions including 2nd -% dimension of a column vector. - -sz=size1(x); -for d=length(sz)-(0:n-1); - x=sum(x,d)/sz(d); -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/minmax.m --- a/general/numerical/minmax.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,15 +0,0 @@ -function R=minmax(X,I) -% minmax - return minimum and maximum along a particular dimension -% -% minmax :: [[N,M]], 1:natural -> [[2,M]]. -% minmax :: [[N,M]], 2:natural -> [[N,2]]. -% minmax :: [D:[[1,E]]], I:[E] -> [set(D,I,2)]. -% -% The most general type means that the return array is the same size -% as the input array except that the size along the Ith dimension -% becomes 2, first element is min, second is max. -% -% The functions is constructed so that it is idemponent: -% minmax(X,I) == minmax(minmax(X,I),I) - -R= cat(I,min(X,[],I),max(X,[],I)); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/mmax.m --- a/general/numerical/mmax.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,13 +0,0 @@ -function [t,i,j]=mmax(A) -% mmax - maximum over both dimensions -% -% mmax :: [[N,M]-A] -> A, 1..N, 1..M. - -if nargout==1, t=max(A(:)); -else - [t1,i1]=max(A,[],1); - [t2,i2]=max(t1,[],2); - t=t2; - if nargout==3, i=i1(i2); j=i2; - elseif nargout==2, i=[i1(i2) i2]; end -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/mmin.m --- a/general/numerical/mmin.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -function t=mmin(A) -% mmin - minimum over both dimensions -% -% mmin :: [[N,M]-A] -> A, 1..N, 1..M. -if nargout==1, t=min(A(:)); -else - [t1,i1]=min(A,[],1); - [t2,i2]=min(t1,[],2); - t=t2; - if nargout==3, i=i1(i2); j=i2; - elseif nargout==2, i=[i1(i2) i2]; end -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/mod1.m --- a/general/numerical/mod1.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,3 +0,0 @@ -% mod1 - Like mod but for ring 1..N instead of 0..N-1 -function z=mod1(x,y), z=1+mod(x-1,y); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/mouter.m --- a/general/numerical/mouter.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,14 +0,0 @@ -function B=mouter(varargin) -% mouter - Multidimensional outer product of multiple arrays -% Index domain of result is the concatenation of the index domains -% of the arguments, (with trailing 1s removed). -% -% mouter :: [Size1->A], [Size2->A] -> [[Size1,Size2]->A]. -% mouter :: [Size1->A], [Size2->A], [Size3->A] -> [[Size1,Size2,Size3]->A]. -% etc. -B=1; -for i=1:length(varargin) - B=kron(varargin{i},B); -end -B=reshape(B,cell2mat(cellmap(@size1,varargin))); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/mulrows.m --- a/general/numerical/mulrows.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,3 +0,0 @@ -function y=mulrows(x,k) - y=repmat(k,1,size(x,2)).*x; -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/nan2x.m --- a/general/numerical/nan2x.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function y=nan2x(x,y) -% nan2x - convert nans to something else -% -% nan2x :: real, [Size] -> [Size]. - -y(isnan(y))=x; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/nary.m --- a/general/numerical/nary.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -% binary - return array contain binary sequece over N columns -% -% binary :: M:natural, N:natural -> [[M^N,N]->1..M]. - -function B=nary(M,N) - -if (N==1), B=(1:M)'; -else - b=nary(M,N-1); - m=size(b,1); - B=[ kron((1:M)',ones(m,1)), repmat(b,M,1)]; -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/orthogonalise.m --- a/general/numerical/orthogonalise.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,13 +0,0 @@ -function B1=orthogonalise(B), -% orthogonalise - Orthogonalise a basis matrix -% -% orthogonalise :: [[N,M]] -> [[N,M]]. -% -% Works using SVD. - -[U,S,V] = svd(B,0); -B1 = U*V'; - -% alternative method, seems to be slower -% B1 = B*real((B'*B)^(-0.5)); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/outer.m --- a/general/numerical/outer.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -function Y=outer(X), Y=X*X'; -% outer - outer product for 2D matrix -% -% outer :: [[N,M]] -> [[N,N]]. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/packvec.m --- a/general/numerical/packvec.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,22 +0,0 @@ -function X=packvec(varargin) -% packvec - Pack coordinate values in separate arrays into one big array -% -% packvec :: {[K]->[[E]]} -> [[K,E]]. -% -% There is also a variable argument list form: -% -% packvec :: [[E]], [[E]] -> [[2,E]]. -% packvec :: [[E]], [[E]], [[E]] -> [[3,E]]. -% etc.. - - -if nargin==1 && iscell(varargin{1}) - Y=varargin{1}; -else - Y=varargin; -end -[Y{:}]=promote(Y{:}); -Y=cellmap(@(y)reshape(y,[1 size(y)]),Y); -X=cat(1,Y{:}); - - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/prod1.m --- a/general/numerical/prod1.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,8 +0,0 @@ -function Y=prod1(X) -% prod1 - prod over dimension 1 and shiftdim 1 -% -% prod1 :: [[N M]] -> [[M]]. - - Z=prod(X,1); - S=size(Z); - Y=reshape(Z,[S(2:end) 1]); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/proddims.m --- a/general/numerical/proddims.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,11 +0,0 @@ -function A=proddims(A,dims) -% proddims - Product over multiple dimensions -% -% proddims :: -% [Size:[[1,E]]->A] ~'E dimensional array', -% Dims:[[K]->1..E] ~'dims to prod over' -% -> [Size1->A] :- Size1=arrset(Size,Dims,1). -% -% Result is same size as input except for summed dimensions -% which are collapsed to a size of 1. -for i=dims, A=prod(A,i); end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/prodn.m --- a/general/numerical/prodn.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,5 +0,0 @@ -function Y=prodn(X,n) -% prod1 - prod over dimension 1 and shiftdim 1 -% -% prod1 :: [[N M]] -> [[M]]. -Y=shiftdim(proddims(X,1:n),n); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/product_iterator.m --- a/general/numerical/product_iterator.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,26 +0,0 @@ -function h=product_iterator(f,g) -% product_iterator - compose two iterators to run in parallel -% -% product_iterator :: -% cell { A->A, A } ~'first iterator transformer and initial value', -% cell { B->B, B } ~'second iterator transformer and initial value', -% -> cell { (cell {A,B} -> cell {A,B}), cell {A,B} } ~'combined iterator'. -% -% The state transformer function returns [] to signal termination -% if EITHER of the sub-transformer functions returns []. - - ff=f{1}; gg=g{1}; - h={@prodit,{f{2},g{2}}}; - - function h=prodit(s) - s1=ff(s{1}); - s2=gg(s{2}); - if isempty(s1) || isempty(s2) - h=[]; - else - h={s1,s2}; - end - end -end - - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/qform.m --- a/general/numerical/qform.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -function z=qform(A,u), z=u'*A*u; -% qform - Quadratic form: qform(A,x)=x'*A*x -% -% qform :: [[N,N]], [[N,T]] -> [[T,T]]. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/qforma.m --- a/general/numerical/qforma.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -function z=qforma(A,u), z=sum(u.*(A*u)); -% qforma - Quadratic form (array version) qforma(A,x)=diag(x'*A*x) -% -% qforma :: [[N,N]], [[N,T]] -> [[1,T]]. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/qlog.m --- a/general/numerical/qlog.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,5 +0,0 @@ -function y=qlog(x), I=(x==0); x(I)=1; y=log(x); y(I)=-inf; end -% qlog - quiet log (no warning when x=0) -% -% qlog :: [[D]] -> [[D]]. - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/quadf.m --- a/general/numerical/quadf.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,1 +0,0 @@ -function y=quadf(x,s), y= -(x.^2)/(2*s); end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/range2set.m --- a/general/numerical/range2set.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,1 +0,0 @@ -function Y=range2set(I), Y=I(1):I(2); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/ratesched.m --- a/general/numerical/ratesched.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -function rates=ratesched(a,b,c) -% rates=ratesched(a,b,c): learning rate schedule -% a: start -% b: time to half -% c: number of steps -rates=b*a./(b:(b+c-1)); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/rectify.m --- a/general/numerical/rectify.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,8 +0,0 @@ -function y=rectify(x) -% rectify - 'half wave' rectification, rectify(x) x>0 ? x : 0 -% -% rectify :: [Size->real] -> [Size->nonneg]. - -y = (x>0) .* x; -% alternative implementation? -% i=(x<0); y=x; y(i)=0; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/rot45.m --- a/general/numerical/rot45.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,17 +0,0 @@ -function B = rot45(A) -% rot45 - rotate square array by 45 degrees -% -% rot45 :: [[N,N]] -> [[N,N]]. -% -% Array is SUBSAMPLED, so, you should pre-filter to avoid aliasing. -% Result is padded with zeros in the corners. - -n = length(A); - -B = diag(A,0)'; -pad = [0]; - -for k=2:2:n-1 - B = [ pad diag(A,k)' pad; B; pad diag(A,-k)' pad]; - pad = [pad 0]; -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/rpsi.m --- a/general/numerical/rpsi.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -function Y=rpsi(a,b) -% rpsi - Version of PSI that returns the same shape array as supplied -% -% rpsi :: [[Size]->real] -> [[Size]->real]. -if nargin==1 - Y=reshape(psi(a),size(a)); -else - Y=reshape(psi(a,b),size(b)); -end - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/sbitmap.m --- a/general/numerical/sbitmap.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,13 +0,0 @@ -function X=sbitmap(c,K) -% sbitmap - turn sequence of integers into sparse binary bitmap -% -% sbitmap :: -% Y:[[1,L]->[K]] ~'sequence of L natural numbers in 1..K', -% K:natural ~'height of array to return' -% -> X:[[K,L]->0|1] ~'X(i,j)=1 if Y(j)=i'. -% -% See also BITMAP - -L=length(c); -if nargin<2, K=max(c); end -X=sparse(c,1:L,1,K,L); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/clamp.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/clamp.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,8 @@ +function x=clamp(a,b,x) +% clamp - Clamp values into given (closed) interval +% +% clamp :: real, real, [Size] -> [Size]. + +x(xb)=b; +%y=min(b,max(a,x)); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/deta.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/deta.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,180 @@ +function [f] = deta(z,k) +% deta - Calculates Dirichlet functions of the form +% +% f = sum((-1)^n/(k*n+1)^z) +% +% over the entire complex plane +% Z may be complex and any size +% Best accuracy for Abs(z) < 100 +% +% Usage: f = deta(z) +% or f = deta(z,k) +% +% where k determines which Dirichlet function to sum +% For Eta (Zeta, Lambda): k=1 +% For Betad: k=2 +% +% This function can use a LOT of memory when size(z) +% is large. Consider using the Memory and Pack commands. +% Also, consider breaking z up into smaller chunks. +% +% Requires a complex Gamma routine. +% Tested under version 5.3.1 +% +%see also: Zeta, Eta, Lambda, Betad +%see also: sym/zeta.m +%see also: mhelp zeta + +%Andrew Odlyzko has Riemann Zeta critical line zeros listed on: +%http://www.research.att.com/~amo/zeta_tables/index.html + +%Paul Godfrey +%pgodfrey@conexant.com +%March 24, 2001 + +if nargin==1 + k=1; +end +k=k(1); +if k<1 | k>2 + warning('Unknown function being calculated! Results valid only for Real(z)>0.5') +% k=1 --> Eta --> Zeta or Lambda --> Bern numbers +% k=2 --> Betad --> Euler numbers +end + +[sizz]=size(z); +z=z(:).'; % make z a row vector +rz=real(z); +iz=imag(z); +iszero=find(z==0); +iseven=find(z==(round(z/2)*2) & rz<0 & iz==0); +isodd= find(z==(round((z-1)/2)*2+1) & rz<0 & iz==0); + +r=1/2; % reflection point +L=find(rz< r); +if ~isempty(L) + zL=z(L); + z(L)=1-zL; +end + +%series coefficients were calculated using +% c(m)=sum from n=m to 64 of (binomial(n,m)/2^n) for m=0:64 + +%coefficients are symmetrical about the 0.5 value. Each pair sums to +-1 +%abs(coefficients) look like erfc(k*m)/2 due to binomial terms +%sum(cm) must = 0.5 = eta(0) = betad(0) +%worst case error occurs for z = 0.5 + i*large + +cm= [ .99999999999999999997; + -.99999999999999999821; + .99999999999999994183; + -.99999999999999875788; + .99999999999998040668; + -.99999999999975652196; + .99999999999751767484; + -.99999999997864739190; + .99999999984183784058; + -.99999999897537734890; + .99999999412319859549; + -.99999996986230482845; + .99999986068828287678; + -.99999941559419338151; + .99999776238757525623; + -.99999214148507363026; + .99997457616475604912; + -.99992394671207596228; + .99978893483826239739; + -.99945495809777621055; + .99868681159465798081; + -.99704078337369034566; + .99374872693175507536; + -.98759401271422391785; + .97682326283354439220; + -.95915923302922997013; + .93198380256105393618; + -.89273040299591077603; + .83945793215750220154; + -.77148960729470505477; + .68992761745934847866; + -.59784149990330073143; + .50000000000000000000; + -.40215850009669926857; + .31007238254065152134; + -.22851039270529494523; + .16054206784249779846; + -.10726959700408922397; + .68016197438946063823e-1; + -.40840766970770029873e-1; + .23176737166455607805e-1; + -.12405987285776082154e-1; + .62512730682449246388e-2; + -.29592166263096543401e-2; + .13131884053420191908e-2; + -.54504190222378945440e-3; + .21106516173760261250e-3; + -.76053287924037718971e-4; + .25423835243950883896e-4; + -.78585149263697370338e-5; + .22376124247437700378e-5; + -.58440580661848562719e-6; + .13931171712321674741e-6; + -.30137695171547022183e-7; + .58768014045093054654e-8; + -.10246226511017621219e-8; + .15816215942184366772e-9; + -.21352608103961806529e-10; + .24823251635643084345e-11; + -.24347803504257137241e-12; + .19593322190397666205e-13; + -.12421162189080181548e-14; + .58167446553847312884e-16; + -.17889335846010823161e-17; + .27105054312137610850e-19]; + +cm=flipud(cm).'; % sum from small to big +nmax=size(cm,2); +n=(1:k:k*nmax)'; +n=flipud(n); +% z is a LR vector +% n is an UD vector +[Z,N]=meshgrid(z,n); + +% this can take a LOT of memory +f=cm*(N.^-Z); +% but it's really fast to form the series expansion N.^-Z +% and then sum it by an inner product cm*() :) + +%reflect across 1/2 +if ~isempty(L) + zz=z(L); + if k==1 + % Eta function reflection + % for test: deta(1,1) should = log(2) + t=(2-2.^(zz+1))./(2.^zz-2)./pi.^zz; + f(L)=t.*cos(pi/2*zz).*gamma(zz).*f(L); + if ~isempty(iszero) + f(iszero)=0.5; + end + if ~isempty(iseven) + f(iseven)=0; + end + end + if k==2 + % Betad function reflection + %for test: deta(0,2) should = 0.5 + %for test: deta(1,2) should = pi/4 + f(L)=(2/pi).^zz.*sin(pi/2*zz).*gamma(zz).*f(L); + if ~isempty(isodd) + f(isodd)=0; + end + end + if k>2 + % Insert reflection formula for other Dirichlet functions here + f(L)=(1/pi).^zz.*gamma(zz).*f(L); + f(L)=NaN; + end +end + +f = reshape(f,sizz); +return + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/erfp.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/erfp.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,9 @@ +function y=erfp(x) +% erfp - returns Pr(X>x) for a normalised Gaussian random variable +% +% erfp :: X:real -> Z:real. +% +% Y = integeral from X to infinity of +% (1/sqrt(2*pi)) * exp(-0.5*t.^2) + +y = erfc(x/sqrt(2))/2; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/erfratio.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/erfratio.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,5 @@ +function y=erfratio(x) +% erfratio - returns gausspdf(x)./erfp(x) but using a better computation +% +% erfratio :: real -> nonneg. +y = 2./(sqrt(2*pi)*erfcx(x/sqrt(2))); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/expdecay.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/expdecay.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,9 @@ +function lambda=expdecay(k,d) +% expdecay - compute interpolation ratio for exponential decay +% +% expdecay :: [Size->nonneg], [Size->nonneg] -> [Size->nonneg]. +% +% this gives exponential convergence + +if nargin==2, lambda=1-k; +else lambda=@(d)1-k; end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/expquadf.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/expquadf.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,5 @@ +% gauss_fn - Gaussian kernel function, exp(-t^2/2) +% +% gauss_fn :: [D] ~'any size array', nonneg ~'variance' -> [D->noneg]~'same size result'. +% gauss_fn :: [D] ~'any size array' -> [D->noneg]~'same size result'. +function x=gauss_fn(t,s), if nargin<2, x=exp(quadf(t)); else x=exp(quadf(t,s)); end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/fromdbs.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/fromdbs.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,7 @@ +function X=fromdbs(Y) +% fromdbs - Convert from decibels +% +% fromdbs :: real -> real. + +X=10.^(Y/10); + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/hshrink.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/hshrink.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,7 @@ +% shrink - Shrinkage towards zero +% +% shrink :: real, [Size] -> [Size]. +function y=shrink(th,x) + +y=x; +y(y nonneg. + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/hsphere_surfln.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/hsphere_surfln.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,7 @@ +function S=hsphere_surfln(N), S=log(2)+(N/2)*log(pi)-gammaln(N/2); end +% hsphere_surfln - log of surface area of unit hypersphere +% +% hsphere_surfln :: natural -> real. +% +% hsphere_surfln(N)=log(hsphere_surf(N)) + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/hypdecay.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/hypdecay.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,17 @@ +function ret=hypdecay(k,d) +% hypdecay - compute interpolation ratio for hyperbolic decay +% +% hypdecay :: [Size->nonneg], [Size->nonneg] -> [Size->nonneg]. +% hypdecay :: [Size->nonneg] -> ([Size->nonneg] -> [Size->nonneg]). +% +% this gives 'hyperbolic' convergence which is more like +% a sort of diffusion by Brownian motion. The trick is +% to add a constant to the inverse of each natural parameter, +% which is like a temperature or variance. The constant +% is like a diffusion constant +% +% This function supports partial application: if only one +% argument is supplied, it returns a function handle. + +if nargin==2, ret=1./(1+k*d); +else ret=@(d)1./(1+k*d); end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/inf2nan.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/inf2nan.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,6 @@ +function x=inf2nan(x) +% inf2nan - convert infs to nans +% +% inf2nan :: [Size] -> [Size]. + +x(isinf(x))=nan; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/iseven.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/iseven.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,2 @@ +% iseven :: integer -> bool. +function b=iseven(n), b=~mod(n,2); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/isodd.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/isodd.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,2 @@ +% isodd :: integer -> bool. +function b=isodd(n), b=mod(n,2); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/lambertw.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/lambertw.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,12 @@ +function W = lambertw(x) +% lambertw - Lambert's W function +% +% lambertw :: real -> real. +% +% evaluates W(y) for -1 branch of Lambert W function where +% y = -exp(-x) +% +% The W function satisfies W(t)*exp(W(t)) = t + +if x<1, error('lambertw: argument must be >= 1'); end +W=fixpoint(@(w)(-x-log(abs(w))),-x,'its',60); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/linterp.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/linterp.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,10 @@ +function z=linterp(k,z0,z) +% linterp - linear interpolation between values +% +% linterp :: +% [Size] ~'fraction of distance to travel from from initial to final point', +% [Size] ~'initial point', +% [Size] ~'final point' +% -> [Size] ~'interpolated point'. + +z = z0+k.*(z-z0); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/linterp_d.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/linterp_d.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,14 @@ +function z=linterp_d(k,z0,d) +% linterp_d - linear interpolation between two points +% +% linterp_d :: +% real ~'fraction of distance to travel from from initial to final point', +% [Size] ~'initial point', +% [Size] ~'vector from initial to final point' +% -> [Size] ~'interpolated point'. +% +% Definition is linterp_d(k,z0,d) = linterp(k,z0,z0+d) +% +% See also: LINTERP, VLINTERP + +z = z0+k.*d; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/log2pie.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/log2pie.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,5 @@ +function k=log2pie + persistent K; + if isempty(K) K=1+log(2*pi); end + k=K; +end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/logth.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/logth.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,6 @@ +function y=logth(eps,x) +% logth - Natural log with lower threshold +% +% logth :: real~'lower threshold', real -> real. + +x(x [Size]. + +y(isnan(y))=x; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/qlog.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/qlog.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,5 @@ +function y=qlog(x), I=(x==0); x(I)=1; y=log(x); y(I)=-inf; end +% qlog - quiet log (no warning when x=0) +% +% qlog :: [[D]] -> [[D]]. + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/quadf.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/quadf.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,7 @@ +% quadf - scalar quadratic form +% +% quadf :: [Size], real -> [Size]. +% quadf :: [Size] -> [Size]. +% +% Scalar function: quadf(x,s) = -x^2/(2*s), s defaults to 1. +function y=quadf(x,s), if nargin<2, y=-(x.^2)/2; else y= -(x.^2)/(2*s); end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/rectify.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/rectify.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,8 @@ +function y=rectify(x) +% rectify - 'half wave' rectification, rectify(x) x>0 ? x : 0 +% +% rectify :: [Size->real] -> [Size->nonneg]. + +y = (x>0) .* x; +% alternative implementation? +% i=(x<0); y=x; y(i)=0; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/rpsi.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/rpsi.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,10 @@ +function Y=rpsi(a,b) +% rpsi - Version of PSI that returns the same shape array as supplied +% +% rpsi :: [[Size]->real] -> [[Size]->real]. +if nargin==1 + Y=reshape(psi(a),size(a)); +else + Y=reshape(psi(a,b),size(b)); +end + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/shrink.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/shrink.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,6 @@ +function y=shrink(th,x) +% shrink - Shrink values towards zero +% +% shrink :: nonneg~'threshold', [Size] -> [Size]. + +y=max(0,x-th); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/slog.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/slog.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,8 @@ +function y=slog(x) +% slog - safe log, x<=0 => slog(x)=0 +% +% slog :: real -> real. + + x(x<=0)=1; + y=log(x); + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/todbs.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/todbs.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,14 @@ +function Y=h_todbs(X,dBs) +% todbs - Convert values to decibels +% +% todbs :: +% real ~'dimensionless value', +% real ~'dynamic range in dB' +% -> real ~'value in dBs'. + +if nargin>1, + mmax=max(X(:)); + X=max(X,10^(-dBs/10)*mmax); +end +Y=10*log10(X); + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/tozeros.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/tozeros.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,9 @@ +function a=tozeros(a) +% tozeros - Return array of zeros same size as given +% +% tozeros :: [Size->A] -> [Size->real]. +% +% All entries are zero. + + a=zeros(size(a)); + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/vlinterp.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/vlinterp.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,15 @@ +function z=vlinterp(k,z0,z) +% vlinterp - vectorial linear interpolation between points +% +% vlinterp :: +% real ~'fraction of distance to travel from from initial to final point', +% [Size] ~'initial point', +% [Size] ~'final point' +% -> [Size] ~'interpolated point'. +% +% This returns a point (vector,array etc) on the straight line +% between the given points. +% +% See also: LINTERP which treats each component independently. + +z = z0+k*(z-z0); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/vlinterp_d.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/vlinterp_d.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,15 @@ +function z=vlinterp_d(k,z0,d) +% vlinterp_d - linear interpolation between two points +% +% vlinterp_d :: +% real ~'fraction of distance to travel from from initial to final point', +% [Size] ~'initial point', +% [Size] ~'vector from initial to final point' +% -> [Size] ~'interpolated point'. +% +% Definition is vlinterp_d(k,z0,d) = vlinterp(k,z0,z0+d) +% +% See also: LINTERP, VLINTERP + + +z = z0+k*d; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/scalar/zeta.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/scalar/zeta.m Thu Jan 17 13:20:44 2013 +0000 @@ -0,0 +1,34 @@ +function [f] = zeta(z) +% zeta - Riemann Zeta function +% +% zeta :: complex -> complex. +% +%tested on version 5.3.1 +% +% This program calculates the Riemann Zeta function +% for the elements of Z using the Dirichlet deta function. +% Z may be complex and any size. Best accuracy for abs(z)<80. +% +% Has a pole at z=1, zeros for z=(-even integers), +% infinite number of zeros for z=1/2+i*y +% +% +%see also: Eta, Deta, Lambda, Betad, Bern, Euler +%see also: mhelp zeta + +%Paul Godfrey +%pgodfrey@conexant.com +%3-24-01 + +zz=2.^z; +k = zz./(zz-2); + +f=k.*deta(z,1); + +p=find(z==1); +if ~isempty(p) + f(p)=Inf; +end + +return + diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/shrink.m --- a/general/numerical/shrink.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function y=shrink(th,x) -% shrink - Shrink values towards zero -% -% shrink :: nonneg~'threshold', [Size] -> [Size]. - -y=max(0,x-th); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/slog.m --- a/general/numerical/slog.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,8 +0,0 @@ -function y=slog(x) -% slog - safe log, x<=0 => slog(x)=0 -% -% slog :: real -> real. - - x(x<=0)=1; - y=log(x); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/spdiag.m --- a/general/numerical/spdiag.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function A=spdiag(H) -% spdiag - Return elementwise multiplier as sparse diagonal matrix -% -% spdiag :: [[N]] -> [[N,N]]. - -A=diag(sparse(H)); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/squeeze_prod.m --- a/general/numerical/squeeze_prod.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -function A=squeeze_prod(A,dims) -% squeeze_prod - Squeeze out first N dimensions by multiplying -% -% squeeze_prod :: [Size], K:natural -> [Size1] :- Size1=Size(K+1:end). -% -% Defined as squeeze_prod(A,K)=shiftdim(proddims(A,1:K),K) -if ~isempty(dims) - for i=dims, A=prod(A,i); end - A=squeeze(A); -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/squeeze_sum.m --- a/general/numerical/squeeze_sum.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -function A=squeeze_sum(A,dims) -% squeeze_sum - Squeeze out first N dimensions by summing -% -% squeeze_sum :: [Size], K:natural -> [Size1] :- Size1=Size(K+1:end). -% -% Defined as squeeze_sum(A,K)=shiftdim(sumdims(A,1:K),K) -if ~isempty(dims), - for i=dims, A=sum(A,i); end - A=squeeze(A); -end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/sum1.m --- a/general/numerical/sum1.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,8 +0,0 @@ -function Y=sum1(X) -% sum1 - sum over dimension 1 and shiftdim 1 -% -% sum1 :: [[N M]] -> [[M]]. - - Z=sum(X,1); - S=size(Z); - Y=reshape(Z,[S(2:end) 1]); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/sumdims.m --- a/general/numerical/sumdims.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -function A=sumdims(A,dims) -% sumdims - Sum over multiple dimensions -% -% sumdims :: -% [Size:[[1,E]]->A] ~'E dimensional array', -% Dims:[[K]->1..E] ~'dims to sum over' -% -> [Size1->A] :- Size1=arrset(Size,Dims,1). -% -% Result is same size as input except for summed dimensions -% which are collapsed to a size of 1. - -for i=dims, A=sum(A,i); end diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/sumn.m --- a/general/numerical/sumn.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,5 +0,0 @@ -function Y=sumn(X,n) -% sumn - sum over first n dimension and shiftdim n -% -% sumn :: [[N:[[1,L]] M]], L:natural -> [[M]]. -Y=shiftdim(sumdims(X,1:n),n); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/sym.m --- a/general/numerical/sym.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -function B=sym(A), B=A+A'; -% sym - Symmetric part of matrix (/2) -% -% sym :: [[N,N]] -> [[N,N]]. diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/todbs.m --- a/general/numerical/todbs.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,14 +0,0 @@ -function Y=h_todbs(X,dBs) -% todbs - Convert values to decibels -% -% todbs :: -% real ~'dimensionless value', -% real ~'dynamic range in dB' -% -> real ~'value in dBs'. - -if nargin>1, - mmax=max(X(:)); - X=max(X,10^(-dBs/10)*mmax); -end -Y=10*log10(X); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/tozeros.m --- a/general/numerical/tozeros.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,9 +0,0 @@ -function a=tozeros(a) -% tozeros - Return array of zeros same size as given -% -% tozeros :: [Size->A] -> [Size->real]. -% -% All entries are zero. - - a=zeros(size(a)); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/ttimes.m --- a/general/numerical/ttimes.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,11 +0,0 @@ -function C=ttimes(A,B) -% ttimes - inner product for arbitrary dimension arrays -% -% ttimes :: [[I,J]], [[J,K]] -> [[I,J]]. -% -% ie just like mtimes, but I and J can be row vectors - -SA=size1(A); -SB=size1(B); -C=reshape(reshape(A,[],SB(1))*reshape(B,SB(1),[]),tosize([SA(1:end-1) SB(2:end)])); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/vecshift.m --- a/general/numerical/vecshift.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,17 +0,0 @@ -function Y=vecshift(O,X) -% vecshift - refer an array of vectors to a new origin. -% -% This SUBTRACTS the first argument (a vector) from each -% of the vectors (row or column) in the second argument. -% Works in two modes: row vector or column vector mode. -% -% vecshift :: -% [[N,1]] ~'new origin', -% [[N,M]] ~'array of vectors, array domain is M' -% -> [[N,M]] ~'vectors relative to new origin'. -% - -% note: this now works for any pair of arrays where size O -% is an integer fraction of size X in any dimension -Y=X-repmat(O,1,size(X,2)); - diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/vlinterp.m --- a/general/numerical/vlinterp.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,15 +0,0 @@ -function z=vlinterp(k,z0,z) -% vlinterp - vectorial linear interpolation between points -% -% vlinterp :: -% real ~'fraction of distance to travel from from initial to final point', -% [Size] ~'initial point', -% [Size] ~'final point' -% -> [Size] ~'interpolated point'. -% -% This returns a point (vector,array etc) on the straight line -% between the given points. -% -% See also: LINTERP which treats each component independently. - -z = z0+k*(z-z0); diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/vlinterp_d.m --- a/general/numerical/vlinterp_d.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,15 +0,0 @@ -function z=vlinterp_d(k,z0,d) -% vlinterp_d - linear interpolation between two points -% -% vlinterp_d :: -% real ~'fraction of distance to travel from from initial to final point', -% [Size] ~'initial point', -% [Size] ~'vector from initial to final point' -% -> [Size] ~'interpolated point'. -% -% Definition is vlinterp_d(k,z0,d) = vlinterp(k,z0,z0+d) -% -% See also: LINTERP, VLINTERP - - -z = z0+k*d; diff -r 19ec801ee487 -r db7f4afd27c5 general/numerical/zeta.m --- a/general/numerical/zeta.m Wed Jan 16 12:17:09 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,34 +0,0 @@ -function [f] = zeta(z) -% zeta - Riemann Zeta function -% -% zeta :: complex -> complex. -% -%tested on version 5.3.1 -% -% This program calculates the Riemann Zeta function -% for the elements of Z using the Dirichlet deta function. -% Z may be complex and any size. Best accuracy for abs(z)<80. -% -% Has a pole at z=1, zeros for z=(-even integers), -% infinite number of zeros for z=1/2+i*y -% -% -%see also: Eta, Deta, Lambda, Betad, Bern, Euler -%see also: mhelp zeta - -%Paul Godfrey -%pgodfrey@conexant.com -%3-24-01 - -zz=2.^z; -k = zz./(zz-2); - -f=k.*deta(z,1); - -p=find(z==1); -if ~isempty(p) - f(p)=Inf; -end - -return - diff -r 19ec801ee487 -r db7f4afd27c5 general/pair.m --- a/general/pair.m Wed Jan 16 12:17:09 2013 +0000 +++ b/general/pair.m Thu Jan 17 13:20:44 2013 +0000 @@ -1,1 +1,6 @@ +% pair - Make 2-element cell array +% +% pair :: A,B -> cell {A,B}. +% +% The type pair(A,B) can be used as a synonym for cell {A,B} when convenient. function p=pair(a,b), p={a,b}; diff -r 19ec801ee487 -r db7f4afd27c5 sequences/README.txt --- a/sequences/README.txt Wed Jan 16 12:17:09 2013 +0000 +++ b/sequences/README.txt Thu Jan 17 13:20:44 2013 +0000 @@ -162,9 +162,7 @@ once = take 1 drop = drop - repeat = loop . (take 1) foldseq = foldl - meandata = foldl accum init scandatacols f = scanl (y\scancols f (last y)) gather dim = foldl (cat dim) [] seq2cell = foldl \x\y[x {y}] {} = gather . map box where box x = {x}