Mercurial > hg > ishara
changeset 4:e44f49929e56
Adding reorganised general toolbox, now in several subdirectories.
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/algo/gatherer.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,47 @@ +% gatherer - API for autoexpanding arrays for accumulating data +function api=gatherer(size_in,varargin) + opts=prefs('chunk',256,'init',512,'grow',2,'max',1e9,varargin{:}); + n=uint32(0); + cap=opts.init; % initial capacity of buffer + buf=zeros(size_in(1),cap); % buffer + + api.remove=@remove; + api.collect=@collect; + api.length=@curlen; + + if isnan(size_in(2)), + api.append=@append; + elseif size_in(2)>1, + m=uint32(size_in(2)); M=1:m; + api.append=@append_m; + else + api.append=@append_1; + end + + function l=curlen, l=n; end + function remove(rem), n=n-rem; end + function x=collect, x=buf(:,1:n); buf=[]; end + function append_1(x) + if n+1>cap, grow; end % need more room + n=n+1; buf(:,n)=x; + end + + function grow + if n>opts.max, error('maximum capacity exceeded'); end + cap=opts.grow*cap; + buf=repmat(buf,1,opts.grow); + end + + function append_var(x) + chunk=size(x,2); + if n+chunk>cap, grow; end % need more room + buf(:,n+(1:chunk))=x; + n=n+chunk; + end + + function append_m(x) + if n+m>cap, grow; end % need more room + buf(:,n+M)=x; + n=n+m; + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/algo/history_extend.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,51 @@ +function [g,X]=history(X,fields,dim) +% history_extend - return state transformer that logs old state +% +% history_extend :: +% struct A ~'any structure', +% F:{[N]->string} ~'list of field names', +% natural ~'dimension to cat along (defaults to 1)' +% -> (struct hist(A,F) ->struct hist(A,F)), +% struct hist(A,F). +% +% struct hist(A,F) ~ 'structure A with history of fields F' ::= +% struct { +% A ~'original fields from A'; +% history :: struct F ~ 'sub-struct containing fields F' +% }. +% +% This function returns a state transformer function. The state MUST be +% a structure but it can have any fields. The history mechanism adds the +% field 'history' which contains the histories of certain fields from +% the main structure. The returned function concatenates the current +% values of the named fields onto the end of the corresponding fields +% in the history, if X contains a field called 'cost', then +% +% [logcost,X]=history(X,{'cost'},1); +% Y=iterate(compose(logcost,update),X,'its',100); +% +% iterates the update function while maintaining a history of cost values. + + g=@logfields; + if nargin<3, dim=1; end + + % prepare X to receive history + if ~isfield(X,'history'), X.history=struct; end + for i=1:length(fields) + % initialise history fields if not present already + if ~isfield(X.history,fields{i}), + X.history.(fields{i}) = []; + end + end + + function X=logfields(X) + for i=1:length(fields) + field=fields{i}; + X.history.(field) = cat(dim,X.history.(field),X.(field)); + end + end +end + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/algo/history_nest.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,61 @@ +function [f,Y]=history_nest(fields,dim,H) +% history_nest - return state transformer that logs old states +% +% history :: +% F:{[N]->string} ~'list of field names', +% natural ~'dimension to cat along (defaults to 1)', +% struct F ~'old history to build on' +% -> (history(A,F) ->history(A,F)), +% history(A,F). +% +% history(A,F) ~ 'structure A with history of fields F' ::= +% struct { +% current :: struct A ~'original fields from A'; +% history :: struct F ~ 'sub-struct containing fields F' +% }. +% +% This function returns a state transformer function. The state MUST be +% a structure but it can have any fields. The history mechanism adds the +% field 'history' which contains the histories of certain fields from +% the main structure. The returned function concatenates the current +% values of the named fields onto the end of the corresponding fields +% in the history, if X is a sequence whose elements contain a field +% called 'cost', then +% +% [logcost,H0]=history({'cost'},1); +% Y=scandata(logcost,H0,X); +% +% Elements of Y will have fields 'current' and 'history' which +% contain values from X and a history of 'cost' respectively. + + g=@logfields_scanner; + if nargin<2, dim=1; end + if nargin<3, H=struct; end + + for i=1:length(fields) + % initialise history fields if not present already + if ~isfield(H,fields{i}), + H.(fields{i}) = []; + end + end + + f=@logfields; + Y.current=[]; + Y.history=H; + + function Y=logfields(Y,X) + Y.current=X; + Y.history=logfields_scanner(Y.history,X); + end + + function H=logfields_scanner(H,X) + for i=1:length(fields) + field=fields{i}; + H.(field) = cat(dim,H.(field),X.(field)); + end + end +end + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/algo/iterate.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,97 @@ +function x0=iterate(varargin) +% iterate - Simulate recursive iteration of a function +% +% iterate :: +% (A=>(A | empty)) ~'function to iterate, can return [] to terminate', +% A ~'initial state', +% options iterate_timed_options +% -> A ~'final state'. +% +% iterate :: +% cell { +% (A->(A | empty)) ~'function to iterate, can return [] to terminate', +% A ~'initial state' +% } ~'iterated system as a cell array', +% options iterate_timed_options +% -> A ~'final state'. +% +% iterate_timed_options = { +% its :: natural /inf ~'Iteration limit'; +% quiet :: bool /1 ~'controls console display'; +% drawnow :: bool /1 ~'controls graphics updating'; +% pre :: A=>void ~'something to do before each iteration'; +% post :: A=>void ~'something to do after each iteration'; +% save :: natural /0 ~'iterations between state saving (0=off)'; +% id :: string /'iterstate@' ~'name of mat file to save to'; +% recover :: bool /0 ~'recover from saved state' +% }. +% +% Plus the usual OPTPAUSE options. +% +% NOTE: old termination by throwing an exception is no longer supported. +% The function must return an empty array to force termination. + +if iscell(varargin{1}) + [f,x]=varargin{1}{:}; + options=varargin(2:end); +else + f=varargin{1}; + x=varargin{2}; + options=varargin(3:end); +end + +opts=prefs('its',inf,'quiet',1,'save',0,'recover',0,'id','iterstate@',options{:}); +pauserfn=pauser(opts); +post=isfield(opts,'post'); +pre=isfield(opts,'pre'); + +i=0; +if opts.recover + if exist([opts.id,'.mat'],'file') + fprintf('Attempting to recover from saved state %s.mat\n',opts.id); + load(opts.id); + fprintf('Continuing from iteration %d\n',i); + else + fprintf('Recovery file %s.mat does not exist.\n',opts.id); + end +end + +if ~post && ~pre && opts.quiet && opts.save==0 + while ~isempty(x) && i<opts.its, + i=i+1; + x0=x; x=f(x0); + pauserfn(); + end +else + if opts.save>0 + fprintf('Will save every %d iterations to %s.mat\n',opts.save,opts.id); + end + + % how cunning is this... + if isfinite(opts.its), itformat=sprintf('%%%dd. ',ceil(log10(opts.its))); + else itformat='%d. '; end + + + while ~isempty(x) && i<opts.its, + i=i+1; + if ~opts.quiet, fprintf(itformat,i); end + if pre, opts.pre(x); end + x0=x; x=feval(f,x0); + if post, opts.post(x); end + + % save state every opts.save iterations + if mod(i,opts.save)==0, + fprintf('\nsaving at %d.\n',i); + save(opts.id,'x','i'); + end + pauserfn(); + end + if opts.save>0, + fprintf('Deleting recovery file %s.mat\n',opts.id); + delete([opts.id '.mat']); + end +end +if ~isempty(x), x0=x; end + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/algo/iterate_gui.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,47 @@ +function iterate_gui(stfn,S0,varargin) +% iterate_gui - iterate state transformer under GUI control +% +% iterate_gui :: +% (A->A,handle) ~'state transformer function', +% A ~'initial state' +% -> gui. +% +% The idea is that the state transformer returns a handle to a +% GUI element to be used for advancing the sequence. You click +% on that element to move to the next element of the sequence. +% Button 2 rewinds to the beginning of the sequence. + +% to do: user configurable keys to allow multiple plotseqs +% in same figure, each watching different keys. + + opts=prefs('keyctl','set',varargin{:}); + + [S,h]=stfn(S0); + set(h,'ButtonDownFcn',@btndn); + switch opts.keyctl + case 'set', addkbcallback(gcf), addkbcallback(gcf,@keypress); + case 'add', addkbcallback(gcf,@keypress); + end + + function next() + if isempty(S), beep; + else + [S,h]=stfn(S); + set(h,'ButtonDownFcn',@btndn); + end + end + + function rewind(), [S,h]=stfn(S0); set(h,'ButtonDownFcn',@btndn); end + function btndn(a,b) + if strcmp(get(gcf,'SelectionType'),'alt'), rewind(); + else next(); end + end + + function keypress(b) + switch b.Character + case {'b','r'}, rewind(); + case {' ','n'}, next(); + end + end +end % of main function +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/algo/iterate_sched.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,24 @@ +function tm=iterate_sched(nextfn,X,T,varargin) +% iterate_sched - Iterate function under control of timer +% +% iterate_sched :: +% (A=>A) ~'state transformer action', +% A ~'initial state', +% real ~'time between updates in seconds' +% options { +% drawnow :: {0,1} /0 ~'call drawnow after each iteration'; +% busy_mode :: {'queue','drop'} /'queue' ~'See TIMER'; +% its :: natural / inf ~'iteration limit' +% -> timer. + + opts=prefs('its',inf,'drawnow',0,'busy_mode','queue','props',{},varargin{:}); + if opts.drawnow, post=@drawnow; else post=@nop; end + + tm=rsched(@action,nows,T,opts.its,iterdata(nextfn,X)); + function [err,s]=action(t_sched,dt,s) + err=nows-t_sched; + s=nextfn(s); + post(); + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/algo/iterate_timed.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,126 @@ +function api=iterate_timed(varargin) +% iterate_timed - Iterate function under control of timer +% +% iterate_timed :: +% (A=>A) ~'state transformer action', +% A ~'initial state', +% real ~'timer between updates in seconds', +% options iterate_timed_options ~'options (see below)' +% -> sched(A) ~'scheduler API'. +% +% iterate_timed :: +% cell { +% (A=>A) ~'state transformer action', +% A ~'initial state' +% } ~'iterated system as a cell array', +% real ~'timer between updates in seconds', +% options iterate_timed_options ~'options (see below)' +% -> sched(A) ~'scheduler API'. +% +% iterate_timed_options = { +% drawnow :: {0,1} /0 ~'call drawnow after each iteration'; +% busy_mode :: {'queue','drop'} /'queue' ~'See TIMER'; +% exec_mode :: {'fixedRate','fixedDelay','fixedSpacing'} /'fixedRate' ~'See TIMER'; +% its :: natural / inf ~'iteration limit'; +% onstart :: A->void / @nop ~'do this when timer starts'; +% onstop :: A->void / @nop ~'do this when timer stops'; +% onfinish :: A->void / @nop ~'do this when iterator terminates'; +% pre :: A->void / [] ~'something to do before each iteration'; +% post :: A->void / [] ~'something to do after each iteration'; +% drawnow :: bool / 0 ~'whether to update graphics every it'; +% defer :: bool / 0 ~'if true then don't start the timer'; +% props :: {[[_]]} / {} ~'extra name/value pairs for timer' +% } +% +% sched(A) ::= struct { +% dispose :: void => void; +% isrunning :: void => bool; +% startat :: real => void; +% start :: void => void; +% stop :: void => void; +% rewind :: void => void; +% getstate :: void => A; +% setstate :: A => void +% }. +% +% NB: Unlike ITERATE, this does NOT respect the id, save, and recover properties. +% Neither does it respect the OPTPAUSE property 'pause'. + + warning('off','MATLAB:TIMER:STARTDELAYPRECISION'); + warning('off','MATLAB:TIMER:RATEPRECISION'); + + if iscell(varargin{1}) + [nextfn,X0]=varargin{1}{:}; + T=varargin{2}; + options=varargin(3:end); + else + nextfn=varargin{1}; + X0=varargin{2}; + T=varargin{3}; + options=varargin(4:end); + end + + opts=prefs('its',inf,'drawnow',0, 'quiet', 0, 'defer', 0, ... + 'busy_mode','queue','exec_mode','fixedRate', ... + 'onfinish',@nop,'onstart',@nop,'onstop',@nop,options{:}); + + X=X0; %% !!! NB: this is a mutable shared local variable used by timercb + + if ~isfield(opts,'pre') && ~isfield(opts,'post') && ~opts.drawnow + stfn=nextfn; % streamlined version with no checks + else % version with calls to pre,post, and drawnow + prefn=getparam(opts,'pre',@nop); + postfn=getparam(opts,'post',@nop); + stfn=@bracket; + end + + tm=timer('ExecutionMode',opts.exec_mode,'BusyMode',opts.busy_mode, ... + 'TimerFcn',@timercb,'StartFcn',@startfn,'StopFcn',@stopfn, ... + 'TasksToExecute',opts.its,'Period',T); + + api=struct( ... + 'dispose', @()delete(tm), ... + 'isrunning',@()isrunning(tm), ... + 'startat', @(t)startat(tm,t/86400), ... + 'start', @()start(tm), ... + 'stop', @()stop(tm), ... + 'rewind', @()setstate(X0), ... + 'getstate', @getstate, ... + 'setstate', @setstate ... + ); + + if ~opts.defer, start(tm); end + + function setstate(x), X=x; end + function x=getstate, x=X; end + + function X=bracket(X), + prefn(X); X=nextfn(X); postfn(X); + if opts.drawnow, drawnow; end + end + + function startfn(tm,a) + if ~opts.quiet, status(tm,'starting',X,a.Data.time); end + opts.onstart(X); + end + + function stopfn(tm,a) + if ~opts.quiet, status(tm,'stopping',X,a.Data.time); end + opts.onstop(X); + if isempty(X), + if ~opts.quiet, status(tm,'finished',X,a.Data.time); end + opts.onfinish(X); + end + end + + % well, we have mutable variables captured by closure so we may + % as well use them... + function timercb(tm,ev) + if isempty(X), stop(tm); else, X=stfn(X); end + end +end + +function status(tm,msg,x,t), + fprintf('| %s %s at %s with %s\n',get(tm,'name'),msg,mat2str(t),tostring(x)); +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/algo/iterate_timed2.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,65 @@ +function Sched=iterate_timed(nextfn,X0,T,varargin) +% iterate_timed - Iterate function under control of timer +% +% iterate_timed :: +% (A=>A) ~'state transformer action', +% A ~'initial state', +% real ~'timer between updates in seconds' +% options { +% drawnow :: {0,1} /0 ~'call drawnow after each iteration'; +% busy_mode :: {'queue','drop'} /'queue' ~'See TIMER'; +% exec_mode :: {'fixedRate','fixedDelay','fixedSpacing'} /'fixedRate' ~'See TIMER'; +% its :: natural / inf ~'iteration limit'; +% onstart :: (A=>void)/ @nop ~'do this when timer starts'; +% onstop :: (A=>void)/ @nop ~'do this when timer stops'; +% onfinish :: (A=>void)/ @nop ~'do this when iterator terminates'; +% pre :: (A=>void)/ @nop ~'do before each iteration'; +% post :: (A=>void)/ @nop ~'do after each iteration'; +% defer :: bool / 0 ~'if true then don't start the timer'; +% -> timer, (A=>void) ~'function to seek to given state'. +% +% NB: Unlike ITERATE, this does NOT respect the id, save, and recover properties. +% Neither does it respect the OPTPAUSE property 'pause'. + + opts=prefs('its',inf,'drawnow',0, 'quiet', 0, 'defer', 0, ... + 'busy_mode','queue','exec_mode','fixedRate', ... + 'onfinish',@nop,'onstart',@nop,'onstop',@nop,varargin{:}); + + if ~isfield(opts,'pre') && ~isfield(opts,'post') && ~opts.drawnow + stfn=nextfn; % streamlined version with no checks + else % version with calls to pre,post, and drawnow + prefn=getparam(opts,'pre',@nop); + postfn=getparam(opts,'post',@nop); + stfn=@bracket; + end + + Sched=rsched(@schedfn,X0,T,0,opts,'defer',1, ... + 'onstart',@onstart,'onstop',@onstop,'onfinish',@onfinish); + + if ~opts.defer, Sched.start(); end + + function [s,t_sched]=schedfn(s,dt,t_sched,t_actual), s=stfn(s); end + + function X=bracket(X), + prefn(X); X=nextfn(X); postfn(X); + if opts.drawnow, drawnow; end + end + + function onfinish(X,t_actual), opts.onfinish(X); end + + function onstart(X,t_sched,t_actual) + if ~opts.quiet, status('starting',X,t_actual); end + opts.onstart(X); + end + + function onstop(X,t_actual) + if ~opts.quiet, status('stopping',X,t_actual); end + opts.onstop(X); + end +end + +function status(msg,x,t), + fprintf('| %s at %s with %s\n',msg,mat2str(t),tostring(x)); +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/algo/optpause.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,14 @@ +function optpause(opt) +% optpause - Optionally pause and/or update graphics +% +% optpause :: struct { +% pause :: natural/0 ~'0-nopause, 1-pause, N-timed pause in ms'; +% drawnow :: bool/0 ~'flush graphics before pausing' +% } => void. + +ps=getparam(opt,'pause',0); +dn=getparam(opt,'drawnow',1); +if ps, + if ps>=10, pause(ps/1000); else pause; end +elseif dn, drawnow; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/algo/pauser.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,21 @@ +function f=pauser(varargin) +% pauser - Return function to pause or not depending on options +% +% pauser :: options { +% pause :: natural/0 ~'0 for no pause, 1 pause, >10 for ms pause'; +% drawnow :: {0,1}/1 ~'call drawnow to flush graphics' +% } -> (void => void) ~'function to pause or not'. +% +% pauser(opts) returns a function that does exactly what optpause(opts) +% does. + + opt=prefs(varargin{:}); + ps=getparam(opt,'pause',0); + dn=getparam(opt,'drawnow',1); + if ps, + if ps<10, f=@pause; + else wait=ps/1000; f=@()pause(wait); + end + elseif dn, f=@drawnow; + else f=@nop; + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/arrset.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,19 @@ +function x=arrset(x,i,z) +% arrset - set elements of array and return result +% +% arrset :: +% [(S:[[1,E]])->R] ~'E-dim array of size S', +% { K:[E] -> [D->[S(K)]] } ~'cell array of E index arrays', +% [D->R] ~'data to write' +% -> [S->R]. +% +% Linear indexing version: +% +% arrset :: +% [(S:[[1,E]])->R] ~'E-dim array of size S', +% [D->[prod(S)]] ~'linear indices between 1 and prod(S)', +% [D->R] ~'data to write' +% -> [S->R]. +if iscell(i), x(i{:})=z; +else x(i)=z; end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/cat_sep.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,5 @@ +% cat_sep - Format nonempty list of string with comma separators +% cat_sep :: {[N]->string} -> string. +function s=cat_sep(sep,list) + s=list{1}; for i=2:length(list), s=[s,sep,list{i}]; end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/chunk.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,11 @@ +function B=chunk(A,n,step) +% CHUNK: return chunk from inside matrix +% B=chunk(A,n,hop): returns nth block of rows of size hop from A +% A: source matrix +% n: index of chunk +% hop: size of chunk + +i=n*step+1; +B=A(i:i+step-1,:); + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/col.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,5 @@ +function Y=col(X,I) +% col - returns selected columns +% col :: [[N,M]->A], [[L]->[M]] -> [[N,L]->A]. + +Y=X(:,I);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/defmat.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,14 @@ +function X=defmat(I,w,D) +% defmat - specify matrix as a list of non-zero elements +% +% defmat :: +% [[N,E]->natural] ~'N E-dim array subscripts, 1 per row', +% [[N]] ~'values to place at positions specified by subs', +% D:[[1,E]] ~'size of target array' +% -> [[D]] ~'array of size D'. + +X=accumarray(I,w); +if any(size(X)<D) + J=num2cell(D); + X(J{:})=0; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/extract.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,25 @@ +% EXTRACT - Extract a sub-array +% +% extract :: +% [[size:[E]]], +% n:1..E, +% [[2]->natural]~'start and end indices' +% -> [[size1:[E]]. +% +% Examples (assuming A is 3D): +% extract(A,2,[4 20]) = A(:,4:20,:) +% extract(A,3,[2 10]) = A(:,:,2:10) +function y=extract(x,dim,range) + + persistent colons + + n=ndims(x); + if length(colons)<n, + colons=repmat({':'},1,n); + end + + s=colons(1:n); %cell(1,n); s(:)={':'}; + s{dim}=range(1):range(2); + S.type='()'; + S.subs=s; + y=subsref(x,S);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/findmd.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,11 @@ +% findmd - Multidimensional find +% +% findmd :: [D:[[1,E]]] -> [[N,E]->natural]. +% +% Take an E dimensional array and return subscripts of nonzero +% elements. +function S=findmd(A) + +E=ndims(A); +[SX{1:E}]=ind2sub(size(A),find(A)); +S=cellcat(2,SX);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/flatten.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,6 @@ +function X=flatten(X); +% flatten - convert array into column vector +% +% flatten :: [[size]]->[[prod(size)]] + +X=reshape(X,numel(X),1);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/interleave.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,19 @@ +function z=interleave(x,y) +% interleave - interleave elements of an array + +x=x(:); +y=y(:); + +nx=length(x); +ny=length(y); + +if nx>ny + z(1:2:2*ny-1)=x(1:ny); + z(2:2:2*ny)=y; + z(2*ny+1:nx+ny)=x(ny+1:end); +else + z(1:2:2*nx-1)=x; + z(2:2:2*nx)=y(1:nx); + z(2*nx+1:nx+ny)=y(nx+1:end); +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/join.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,12 @@ +function Z=join(X,Y) +% join - Cartesian product of two arrays +% +% join :: [[N,M]], [[L,P]] -> [[N*L,M+P]]. +% +% The result contains every pairing of rows from the two arguments + +%Z=[repmat(X,size(Y,1),1) kron(Y,ones(size(X,1),1))]; + +% this is probably faster (no multiplications!) +m=size(X,1); n=size(Y,1); +Z=[repmat(X,n,1) Y(repmat(1:n,m,1),:)];
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/join_tr.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,12 @@ +function Z=join_tr(X,Y) +% join_tr - Cartesian product of two arrays, transpose of join +% +% join_tr :: [[M,N]], [[P,L]] -> [[M+P,N*L]]. +% +% The result contains every pairing of rows from the two arguments + +%Z=[repmat(X,size(Y,1),1) kron(Y,ones(size(X,1),1))]; + +% this is probably faster (no multiplications!) +m=size(X,2); n=size(Y,2); +Z=[repmat(X,1,n); Y(:,repmat((1:n),m,1))];
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/lags.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,20 @@ +function Y=lags(X,tau) +% lags - returns lagged versions of sequence in X +% +% lags :: [[N,T]], [[K]->integer]~'K lags' -> [[N,T2,K]]. +% +% T2 is the maxmimum length such that all lagged sequences contain valid data. + +% doesn't work for data sequences yet +if isdata(X), X=head(X); end + +tau=tau-min(tau); % zero base lags +[N T0]=size(X); +M=length(tau); +T=T0-max(tau); +Y=zeros(N,T,M); + +tt=1:T; +for i=1:length(tau) + Y(:,:,i)=X(:,tau(i)+tt); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/mapwindows.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,51 @@ +function m=mapwindows(fn,x,L,M) +% mapwindows - Apply function to windowed signal with causality control +% +% mapwindows :: +% F:([[N]]->[[M]]) ~'function to apply to windowed segments', +% [[T]] ~'signal of length T', +% L:natural ~'post-window (ie causal part)', +% M:natural ~'pre-window (ie anti-causal part)' +% -> [[M,T]] ~'sequence of results'. +% +% This function applies a function to a sequence of segments extracted from +% the given signal by sliding a window along the signal. The +% post- and pre- window lengths determine the positioning of the window: +% +% .... x1 x2 x3 x4 x5 x6 x7 x8 x9 ... +% |<--- post --> | now |<----- pre ------> | +% \________________________________________/ +% | +% function +% | +% ... y3 y4 y5 ... +% +% Note: the function F must make sense when applied to windows of +% different lengths since the first and last few windows will +% be shorter than L+M+1. The windows will not be zero padded! +% +% EXAMPLES +% +% To compute a running mean using only past values (ie causally) +% >> Y=mapwindows(@mean,X,6,0); +% +% To compute a running local median +% >> Y=mapwindows(@median,X,4,3); + +N=length(x); + +for i=1:min(L,N) + k=min(i+M,N); + MM=feval(fn,x(1:k)); + m(:,i)=feval(fn,x(1:k)); +end + +if N>(L+M) + m=[m feval(fn,buffer(x,L+M+1,L+M,'nodelay'))]; +end + +for i=N+(1-M:0) + j=max(i-L,1); + m(:,i)=feval(fn,x(j:end)); +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/mzip.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,13 @@ +% mzip - Apply array to multiple index arrays +% +% mzip :: [[N,M]->D], [[L]->[N]], [[L]->[M]] -> [[L]->D]. +% +% The name mzip is meant to suggest Matrix or Multiple Zip, +% with zip referring to the functional programming sort of +% zip where you can apply a function of N arguments to N +% lists of arguments and get a list of return values. +% eg, mzip(x,[2 4 1],[3 2 2]) = [x(2,3) x(4,2) x(1,2)]. + +function B=mzip(A,varargin), B=A(sub2ind(size(A),varargin{:})); + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/numdims.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,14 @@ +function n=numdims(x) +% numdims - Dimension of multidimensional array +% +% numdims :: [Size] -> natural. +% +% n=numdims(x), where n=length(size1(x)), ie not counting any +% trailing singleton dimensions. Eg +% numdims(1)=0 +% numdims([1;2;3])=1 +% numdims([1 2 3])=2 + +sz=size(x); +n=length(sz); +while n>0 && sz(n)==1, n=n-1; end;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/order.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,6 @@ +function I=order(X) +% order - return sort order of given +% +% order :: [[N]] -> [[N]->[N]]. + +[Y,I]=sort(X);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/orderby.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,6 @@ +function b=orderby(Q) +% orderby - sort vector and return permuted index order rather than sorted values +% +% orderby :: [[N]] -> [[N]->[N]]. + +[a,b]=sort(Q);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/pad.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,8 @@ +function b=pad(x,sz,a), +% pad - pad out an array to given size with given elements +% +% pad :: A, D:[[2]], [[S]->A] -> [D->A]. + +[N,M]=size(a); +b=repmat_to(x,max(sz,[N,M])); +b(1:N,1:M)=a;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/pad1s.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,6 @@ +function b=pad1s(n,a), b=[a ones(1,n-length(a))]; +% pad1s - pad a size vector with 1s till it is a given length +% +% pad1s :: E:natural, [[1,D]] -> [[1,E]]. +% +% Note, pad1s has no effect if vector is already longer than n.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/paren.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,12 @@ +function Y=paren(X,varargin) +% paren - array dereferencing as a function +% +% See also SUBS: paren also accepts ':' and 'end' + +n=length(varargin); +for i=1:n + if strcmp(varargin{i},'end') + varargin{i}=size(X,i); + end +end +Y=X(varargin{:});
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/permutegroups.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,69 @@ +function A=permutegroups(A,doms,perm) +% permutegroups - Reorder dimensions of array in groups +% +% Like PERMUTE, PERMUTEGROUPS rearranges the dimensions of a +% multidimensional array in the same way that transpose swaps +% the 2 dimensions of a matrix. However, rather than specifying +% the permutation in terms of individual dimensions, the +% dimensions are matched to a groups template and the then +% the groups are permuted. Eg, suppose A is an array such that +% the index domain size(A)=[6 3 4 7]. Now suppose the these +% are grouped as [ [6 3] [4 7] ] and we wish to swap the groups. +% Rather than saying +% +% permute(A,[3 4 1 2]) +% +% we can say +% +% permutegroups(A,{[6 3] [4 7]},[2 1]) +% +% where the 2nd parameter specifies the two groups, and the third +% indicates that we want to swap them. Actually, the numbers in the +% group specification are not important - only the lengths of the +% groups matters. +% +% Usage: B=permutegroups(A,groups,perm) +% +% A: orginal array +% groups: list of array subdomains (ie sub-sequences of size(A)) +% One of these may be nan, in which case, all remaining dims +% are assigned to that group. +% perm: desired permutation of groups + +a=0; +dims=cell(size(doms)); +nanat=0; +for i=1:length(doms) + if isempty(doms{i}), dims{i}=[]; + else + if ~isnan(doms{i}) + l=length(doms{i}); + dims{i}=a+(1:l); + a=a+l; + else + nanat=i; + dims{i}=a; + end + end +end + +if nanat>0 + gap=length(size1(A))-a; + if gap>0, + dims{nanat}=dims{nanat}+(1:gap); + for j=nanat+1:length(doms) + if ~isempty(dims{j}), dims{j}=dims{j}+gap; end + end + else + dims{nanat}=[]; + end +else + gap=0; +end + +perm=cell2mat(dims(perm)); +maxdim=a+gap; +if maxdim<1, perm(end+1)=1; end +if maxdim<2, perm(end+1)=2; end + +A=permute(A,perm);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/promote.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,14 @@ +function varargout=promote(varargin) +% promote - repmats all args to be the same size +% +% promote :: var {I:[N]->[[Size(I)]->A]} ~'variable argument list of arrays' +% -> var {I:[[N]]->[[MSize]->A]} ~'all arrays of same size'. +% +% Target size is the maximum along each dimensions of all the input sizes. + +maxdims=max(cellcat(1,cellmap(@numdims,varargin))); +target=max(cellcat(1,cellmap(@(x)pad1s(maxdims,size1(x)),varargin)),[],1); +varargout=cellmap(@(x)repmat_to(x,tosize(target)),varargin); + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/replace.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +function y=replace(y,a,b), y(y==a)=b; +% replace - Search and replace in array +% +% replace :: +% [Size->A] ~'array of any size and type', +% A ~'value to search for', +% A ~'value to replace it with' +% -> [Size->A]. +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/repmat_to.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,14 @@ +function b=repmat_to(a,sz) +% repmat_to(a,sz): Like repmat, but sz is TARGET size, not repliction count +% +% repmat_to :: [[Size1]->A], Size2:dim -> [[Size2]->A]. + +% this will barf if any dimension of a is zero. + if any(sz==0) + b=zeros(sz); + else + sa=size1(a); + l=length(sz)-length(sa); + b=repmat(a,tosize(sz./[sa ones(1,l)])); + end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/revdims.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,18 @@ +function B=revdims(A,N) +% revdims - reverse N dimensions of a multidimensional array +% +% revdims :: [[D:[[1,M]]]], N:natural -> [[D2:[[1,M]]]]. +% revdims :: [[D:[[1,M]]]] -> [[D2:[[1,M]]]]. +% +% This takes an array of size [D(1:N) DX] +% and returns an array of size [D(N:-1:1) DX] +% where DX are the rest of the dimensions after the Nth and are +% unaffected. N defaults to the number of dimensions disregarding +% trailing singletons, ie numdims(A)==length(size1(A)) + +if nargin<2, N=numdims(A); end + +ND=max(N,2); +DD=1:ND; +DD(1:N)=fliplr(DD(1:N)); +B=permute(A,DD);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/reverse.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,3 @@ +function y=reverse(x) +% reverse a vector +y = x(length(x):-1:1);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/rotl.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,8 @@ +function y=rotl(x,n) +% rotl - rotate columns of array to left +% +% rotl :: [[N,M]] -> [[N,M]]. +% rotl :: [[N,M]], natural -> [[N,M]]. + +if nargin<2, n=1; end +y = x(:,[n+1:size(x,2),1:n]);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/rotr.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +function y=rotr(x,n) +% rotl - rotate columns of array to left +% +% rotl :: [[N,M]] -> [[N,M]]. +% rotl :: [[N,M]], natural -> [[N,M]]. + +m=size(x,2); +if nargin<2, n=1; end +y = x(:,mod((0:m-1)-n,m)+1);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/row.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,3 @@ +% row - returns selected rows of a 2D array +% row :: [[N,M]->A], [[L]->[N]] -> [[L,M]->A]. +function Y=row(X,I), Y=X(I,:);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/setdiag.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +function A=setdiag(A,x) +% setdiag - set elements on diagonal to something +% +% setdiag :: [[N,M]], real -> [[N,M]]. +% setdiag :: [[N,N]], [[N]] -> [[N,N]]. + +dom=size(A); +i=1:min(dom); +A(sub2ind(dom,i,i))=x;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/shiftl.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +function y=shiftl(x) +% shiftl - shift a column vector up one +% maintaining original length +% eg [ 1 2 3 4 ] --> [ 2 3 4 0] +if isvector(x) + y = [ x(2:length(x)); 0]; +else + y = [ x(:,2:size(x,2)), zeros(size(x,1),1) ]; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/shiftr.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,5 @@ +function y=shiftr(x) +% shiftr - shift a column vector down one +% maintaining original length +% eg [ 1 2 3 4 ] --> [ 0 1 2 3 ] +y = [0; x(1:length(x)-1) ];
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/shuffle.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,7 @@ +function a=shuffle(a) +% shuffle - randomly shuffle the columns of A +% +% shuffle :: [[N,M]] -> [[N,M]]. +[n,i]=sort(rand(1,size(a,2))); +a=a(:,i); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/size1.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,7 @@ +function s=size1(A) +% size1 - size of array with trailing singleton dimensions stripped +% +% size1 :: [[Size]] -> Size:[[1,E]]. +% +% size is defined as size1(A)=strip1(size(A)). See STRIP1 +s=strip1(size(A));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/sizedims.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,8 @@ +function s=sizedims(A,i) +% sizedims - size of array in specified dimensions only +% +% sizedims :: [Size:[[1,E]]], [[M]->[E]] -> [[M]->natural]. +% +% Defined as sizedims(A,I)=paren(size(A),I) +s=size(A); +s=s(i);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/sizes.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,13 @@ +function S=sizes(A) +% sizes - return sizes of elements in cell array +% +% sizes :: {I:[N]->[Size(I):[[1,E]]] -> [[N,E]]. +% +% The sizes of the N E-dimensional arrays in the input are collected +% into an N-by-E array. All arrays must have the same number of dimensions +% as returned by size. + +S=zeros([length(A) 2]); +for k=1:length(A) + S(k,:)=size(A{k}); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/sli.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,15 @@ +function y=sli(x,i,dim) +% sli - Take a slice of a multidimensional array +% +% sli :: +% [Size:[[E]]->A] ~'E dimensional array of size Size, element type A', +% [[M]->1..Size(Dim)] ~'vector of M indices in Dim'th dimension', +% Dim:1..E ~'dimension to act along' +% -> [Size2:[[E]]->A] ~'new array of size Size2 :- +% Size2=set(Size,Dim,M) ~'new size same as old but with Dim'th entry replaced with M'. + + n=ndims(x); + idx=repmat({':'},1,n); + idx{dim}=i; + y=x(idx{:}); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/strip1.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,7 @@ +function sz=strip1(sz) +% strip1 - strip trailing ones from an array size specification +% +% strip1 :: [[1,E]] -> [[1,F]]. +e=length(sz); +while e>0 && sz(e)==1, e=e-1; end +if e>0, sz=sz(1:e); else sz=[]; end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/submat.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,6 @@ +% submat - extract square submatrix +% +% submat :: [[N,M]], [[L]->[min(N,M)]] -> [[L,L]]. +% +% submat(X,I)=X(I,I) +function B=submat(A,I), B=A(I,I); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/subs.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,11 @@ +function Y=subs(X,varargin) +% subs - Array reference as a function +% +% subs(A,...) = A(...) +% +% Colon can be be used if quoted: +% subs(A,':') = A(:) +% End indexing cannot be used, eg subs(A,'end') does not +% work as A(end) + +Y=X(varargin{:});
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/to1s.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,6 @@ +% to1s - return array of ones same size as given +% +% to1s :: [Size->A] -> [Size->real]. +function a=to1s(a) + a=ones(size(a)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/tosize.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,6 @@ +function sz=tosize(varargin) +% tosize - Works out size parameter for calls to rand, randn etc. +% +% The main idea is to override the stupid default: rand(N) returns +% an NxN matrix, not an N-element vector. +sz=pad1s(2,[varargin{:}]);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/tween.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,11 @@ +function I=tween(a,b) +% tween - all indices in a given range (like interval to set) +% +% tween :: interger, integer -> [[N]->integer]. +% tween :: [[2]->interger] -> [[N]->integer]. +% +% tween([A,B])=tween(A,B). + +if nargin<2, b=a(2); a=a(1); end +if b>=a, I=a:b; else I=a:-1:b; end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/vecop.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,16 @@ +function Z=vecop(F,X,Y) +% vecop - apply binary function to different sized arrays +% +% vecop :: +% ([[D]],[[D]]->[[D]]) ~'some function requiring equal size args', +% [[DX]] ~'first arg of size DX', +% [[DY]] ~'second arg of size DY' +% -> [[DZ]] ~'result of size DZ' :- DZ=max(DX,DY). +% + + +DX=size(X); DY=size(Y); +E=max(length(DX),length(DY)); +DZ=max(pad1s(E,DX),pad1s(E,DY)); +Z=feval(F,repmat_to(X,DZ),repmat_to(Y,DZ)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/arrutils/vecop1.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,13 @@ +function Z=vecop1(F,X,Y) +% vecop1 - apply binary function to different sized arrays +% +% vecop1 :: +% ([[D]],[[D]]->[[D]]) ~'some function requiring equal size args', +% [[DX]] ~'first arg of size DX', +% [[DY]] ~'second arg of size DY' +% -> [[DX]] ~'result of size DX' :- all(DX>=DY). +% +% this is a slightly simpler version of vecop that assumes that the +% first argument is at least as big as the second in every dimension. +Z=feval(F,X,repmat_to(Y,size(X))); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/cellutils/@cell/gather.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +function Y=gather(dim,X,varargin) +% gather - Gather for cell arrays +% +% gather :: +% D:natural ~'dimenion along which to gather', +% cells(A) ~'1 D cell array of things' +% -> [Size->A]. + +if isempty(X), Y=[]; else Y=cellcat(dim,X); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/cellutils/@cell/head.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,4 @@ +% head - Head for cell arrays (first element) +% +% head :: cells(A) -> A. +function y=head(x), y=x{1}; % if iscell(x), y=x{1}; else y=x; end;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/cellutils/@cell/next.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,4 @@ +% next - next for cell arrays (remove first element) +% +% next :: cells(A) -> A. +function y=next(x), y=x(2:end);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/cellutils/cdeal.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,4 @@ +function varargout=cdeal(y), varargout=y; +% cdeal - deal cell array input to multiple return values +% +% cdeal :: cell { A, B, ... } -> A, B, ... .
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/cellutils/cellcat.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +function Y=cellcat(dim,X), Y=cat(dim,X{:}); +% cellcat - join contents of cell array into ordinary array +% +% cellcat :: +% E:natural ~'dimension to concatenate along', +% {[N]->A} ~'cell array of N values' +% -> [Size->A] :- Size=[ones(1,1-E),E]. +% +% It's faster than cell2mat
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/cellutils/cellfilt.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,16 @@ +function Y=cellfilt(fn,X) +% cellfilt - Filter cell array with boolean test +% +% cellfilt :: (A->bool), {[N]->A} -> {[M]->A}. +% cellfilt :: (A->bool), {[1,N]->A} -> {[1,M]->A}. +% +% eg, cellfilt(@iseven,{1,2,3,4}) = {2,4} + +% SA 2008-06-27 - Now returns cell array with same orientation as input + +Y={}; +[X,E]=shiftdim(X); +for i=1:numel(X) + if fn(X{i}), Y=vertcat(Y,X(i)); end +end +Y=shiftdim(Y,E);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/cellutils/cellget.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,10 @@ +function B=cellget(A,varargin) +% cellget - Cell array reference as ordinary function +% +% cellget :: {[N]->A}, 1..N -> A. +% cellget :: {[N,M]->A}, 1..N, 1..M -> A. +% +% etc. + +B=A{varargin{:}}; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/cellutils/cellmap.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,10 @@ +function Y=cellmap(fn,X) +% cellmap - Map a function over a cell array +% +% cellmap :: (A->B, {[Size]->A}) -> {[Size]->B} + +% preallocate to fix size +Y=cell(size(X)); +for i=1:numel(X) + Y{i}=feval(fn,X{i}); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/cellutils/cellrep.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,7 @@ +function z=cellrep(n,a) +% cellrep - Create cell array consisting of copies of give value +% +% cellrep :: N:natural, A -> {[1,N]->A}. + +% z=repmat({a},1,n); some Matlabs don't like this +c={a}; z=c(ones(1,n));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/cellutils/cellset.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,13 @@ +function A=cellset(n,B,A) +% cellset - set element of cell array to new value and return new cell array +% +% cellset :: [N], A, {[N]} -> {[N]}. +% cellset :: N:natural, A -> {[N]}. +% +% The order of arguments is a bit weird here. Also, if the third argument +% is not present, this function returns a cell array which is empty apart +% from the last (Nth) element. + + if nargin<3, A=cell(n,1); end; + A{n}=B; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/cellutils/cellzip.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,23 @@ +function Y=cellzip(fn,varargin) +% cellzip - Zip cell arrays with a function +% +% cellzip :: ( +% (D1,...,Dn)->R ~ function of n arguments, +% {[size]->D1}, ..., {[size]->Dn} ~ n cell arrays of appropriate types +% ) -> {[size]->R} ~ cell array of results + +% preallocate to fix size +sizes=cellmap(@size,varargin); +sz=sizes{1}; +for i=1:nargin-1 + if ~all(sizes{i}==sz), error('Sizes do not match'); end +end + +Y=cell(sz); + +for i=1:prod(sz) + args=cellmap(@(c)c{i},varargin); + Y{i}=feval(fn,args{:}); +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/cpfields.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,11 @@ +% cpfields - Copy selected fields from one structure to another +% +% cpfields :: +% C:{[N]->string} ~'fields names to copy', +% struct A ~'struct with field names in A', +% struct B ~'struct with field names in B' +% -> struct union(B,C) ~'struct with fields B union C' :- +% subset(C,A) ~'target fields must be in A'. +function B=cpfields(fields,A,B) + B=foldl(@(b,fn)setfield(b,fn,getfield(A,fn)),B,fields); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/discretise/@dmap/dmap.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,136 @@ +classdef dmap + properties (GetAccess=private, SetAccess=immutable) + cardr_ + map_rn + map_ni + domain_ + end + + methods + % dmap - Create discretisation map + % + % dmap :: + % N:natural~'range of discretisation function will be 1..N', + % (real->natural) ~'discretisation function itself', + % (natural->[[2]]) ~'reverse map, returns upper and lower bin edges' + % -> dmap. + % + % The map consists of a mapping between the integers 1..N and a set + % of consectutive half-open intervals. The discretisation is performed + % by mapping a real number to the index of the half-open interval in + % which it lies. + % + % The map be applied using parentheses (function application), eg + % + % m = linmap(-2,10,24); + % i = m(2.45); + % + % If the input is outside, the map returns -Inf or Inf. + % + % METHODS + % bins - return some or all of the intervals in the map + % cardr - number of intervals in a map + % centres - return the centres of some or all of the intervals + % domain - return the interval which is union of all the intervals + % edges - return the set of all the interval endpoints. + % range - the set 1..N + % map_clamped - apply map with output clamped to valid range + % feval - apply map to real numbers + % + % dmap also implements subsref, tostring and display + % + % linmap, binmap, edgemap - create convenient discretisation maps + + function M=dmap(card,mapfn,revmap) + if nargin==0, M=dmap(1,@floor,@(i)[i-1;i]); % dummy map + elseif isa(card,'dmap'), M=card + else + M.cardr_=card; + M.map_rn=mapfn; + M.map_ni=revmap; + + intervals=revmap([1 M.cardr_]); + M.domain_=[intervals(1,1) intervals(2,2)]; + end + end + + % bins - return bin edges of discretisation map + % + % bins :: dmap(N) -> [[2,N]]~'all the bins'. + % bins :: dmap(N), [[M]->[N]]~'M bin indices'-> [[2,M]]~'selected bins'. + function X=bins(M,I), + if nargin<2,I=range(M); end + X=M.map_ni(I); + end + + % cardr - Cardinality of range (ie size) + % + % cardr :: dmap(N) -> N:natural. + function N=cardr(M), N=M.cardr_; end + + % centres - Return centres of discretisation bins + % + % centres :: dmap(N) -> [[N]] ~'centres of all bins'. + % centres :: dmap(N), [[M]->[N]] -> [[M]] ~'centres of selected bins'. + function X=centres(M,I) + if nargin<2,I=range(M); end + X=mean(M.map_ni(I),1); + end + + % tostring - string describing dmap + % + % tostring :: dmap(N) -> string. + function s=tostring(M) + s=sprintf('dmap::%g--%g->[%d]',M.domain_(1),M.domain_(2),M.cardr_); + end + + % DISPLAY - Display dmap object as string + function display(c), disp(tostring(c)); end + + % domain - return domain of a dmap + % + % domain :: dmap(N) -> [[2]] + function X=domain(M), X=M.domain_; end + + % edges - return edges of discretisation map + % + % edges :: dmap(N) -> [[N+1]]. + % edges :: dmap(N) [[M]->[N]-> [[M+1]]. + function X=edges(M,I), + if nargin<2,I=range(M); end + Y=M.map_ni(I); + X=[Y(1,1) Y(2,:)]; + end + + function I=feval(M,X), I=M.map_rn(X); end +% if isa(M,'dmap'), I=M.map_rn(X); +% else I=builtin('feval',M,X); end +% end + + % map_clamped - Apply clamped discretisation map + % + % map_clamped :: dmap(N), real -> [N]. + function I=map_clamped(M,x) + I=M.map_rn(x); + I(I==-inf)=1; + I(I==inf) =M.cardr_; + end + + % range - range of dmap (set of natural numbers) + % + % range :: dmap(N) -> 1..N:[[N]->natural]. + function R=range(M), R=1:M.cardr_; end + + % SUBSREF - Subscript referencing for DMAP objects + % + % I=MAP(X) - Apply real-to-natural map to all elements of X + function I=subsref(M,S) + s=S(1); + switch s.type + case '()' + I=M.map_rn(s.subs{1}); + end + end + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/discretise/binmap.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,22 @@ +% binmap - linear discretisation map from bin centres +% +% binmap :: real~'first bin centre', real~'last bin centre', N:natural~'num bins'->dmap(N). +function M=binmap(min,max,N) + + dx=(max-min)/(N-1); + rdx=1/dx; + M=dmap(N,@map,@revmap); + + function I=map(X) + I=1+round(rdx*(X-min)); + I(I<1)=-inf; + I(I>N)=inf; + end + + function X=revmap(I) + I=shiftdim(shiftdim(I),-1); + X1=min+(I-1.5)*dx; + X2=min+(I-0.5)*dx; + X=cat(1,X1,X2); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/discretise/data2maps.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,19 @@ +function M=data2maps(RFn,X,N) +% data2maps - build linear maps for covering data ranges +% +% data2maps :: +% ([[N,M]] -> [[2,M]]) ~ 'function to compute min and max of range', +% [[N,M]] ~ 'the data' +% L:natural ~ 'number of bins in output maps' +% -> { [M] -> dmap(L) } ~ 'one dmap per column in a cell array'. +% +% Two argument version computes L=sqrt(N) +% +% data2maps :: +% ([[N,M]] -> [[2,M]]) ~ 'function to compute min and max of range', +% [[N,M]] ~ 'the data' +% -> { [M] -> dmap(L) } ~ 'one dmap per column in a cell array'. + +if nargin<3, N=sqrt(size(X,1)); end +M=maprows(@(r){linmap(r(1),r(2),N)}, closed2hopen(feval(RFn,X)')); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/discretise/edgemap.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,53 @@ +function M=edgemap(edges) +% edgemap - Create arbitrary discretisation map from bin edges +% +% edgemap :: [[N+1]] -> dmap(N). +% +% See also: dmap + +N=length(edges)-1; +M=dmap(N,@(t)map(edges,t),@(t)revmap(edges,t)); + +function I=map(edges,X) + I=X; % make sure I is same shape as X + min=edges(1); + max=edges(end); + for k=1:numel(X) + x=X(k); + if x<min, I(k)=-inf; + elseif x>=max, I(k)=inf; + else I(k)=bsearch_it(edges,x); end + end + +function X=revmap(edges,I) + X=[edges(I);edges(I+1)]; + + +% bsearch :: real, [[N]] -> 1..N. + +% recursive version. +function i=bsearch(y,x) + l=length(y); + if l==2, i=1; + else + half=ceil(l/2); + if x<y(half), i=bsearch(x,y(1:half)); + else i=bsearch(x,y(half:end))+half-1; + end + end + +% iterative version +function i=bsearch_it(y,x) + off=0; + l=length(y); + while l>2 + half=ceil(l/2); + if x<y(half) + y=y(1:half); + else + y=y(half:end); + off=off+half-1; + end + l=length(y); + end + i=off+1;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/discretise/expmap.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,31 @@ +function M=expmap(p1,p2,p3) +% expmap - Create exponentially spaced discretisation map +% +% expmap :: Min:real, Max:real, N:natural -> dmap(N). +% expmap :: [Min,Max]:[[1,2]], N:natural -> dmap(N). +% +% Parameters work almost but not quite like logspace: the +% third parameter is the number of GAPS, not the number of +% points, defining the map, So for a map with edges at 0:10, +% use linmap(0,10,10). To get the same edges from linspace, +% you would need linspace(0,10,11). + +% unpack parameters +if nargin==2, min=log(p1(1)); max=log(p1(2)); bins=p2; +else min=log(p1); max=log(p2); bins=p3; end + +dx=(max-min)/bins; +M=dmap(bins,@(t)map(bins,min,1/dx,t),@(t)revmap(min,dx,t)); + +function I=map(N,min,k,X) + U=log(X); + I=1+floor(k*(U-min)); + I(I<1)=-inf; + I(I>N)=inf; + +function X=revmap(min,ik,I) + I=shiftdim(shiftdim(I),-1); + U1=min+(I-1)*ik; + U2=min+(I)*ik; + X=exp(cat(1,U1,U2)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/discretise/intmap.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,26 @@ +function M=intmap(min,max) +% intmap - Discretisation map for integers (ie alread discrete) +% +% intmap :: integer~'min value', integer~'max value' -> dmap(N). +% +% Maps integers to natural numbers, saves on multiplications and +% roundings as performed by the real->natural maps. + + N=max-min+1; + off=min-1; + M=dmap(N,@map,@revmap); + + function I=map(X) + I=X-off; + I(I<1)=-inf; + I(I>N)=inf; + end + + function X=revmap(I) + I=shiftdim(shiftdim(I),-1); + X1=off+I-0.5; + X2=off+I+0.5; + X=cat(1,X1,X2); + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/discretise/linmap.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,30 @@ +function M=linmap(p1,p2,p3) +% linmap - Create linearly spaced discretisation map +% +% linmap :: Min:real, Max:real, N:natural -> dmap(N). +% linmap :: [Min,Max]:[[1,2]], N:natural -> dmap(N). +% +% Parameters work almost but not quite like linspace: the +% third parameter is the number of GAPS, not the number of +% points, defining the map, So for a map with edges at 0:10, +% use linmap(0,10,10). To get the same edges from linspace, +% you would need linspace(0,10,11). + +% unpack parameters +if nargin==2, min=p1(1); max=p1(2); bins=p2; +else min=p1; max=p2; bins=p3; end + +dx=(max-min)/bins; +M=dmap(bins,@(t)map(bins,min,1/dx,t),@(t)revmap(min,dx,t)); + +function I=map(N,min,k,X) + I=1+floor(k*(X-min)); + I(I<1)=-inf; + I(I>N)=inf; + +function X=revmap(min,ik,I) + I=shiftdim(shiftdim(I),-1); + X1=min+(I-1)*ik; + X2=min+(I)*ik; + X=cat(1,X1,X2); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/discretise/maps2grid.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,6 @@ +function varargout=maps2grid(M) +% maps2grid - Compute multidimensional grid domain from multiple maps + +E=cellmap(@edges,M); +varargout=cell(nargout,1); +[varargout{:}]=meshgrid(E{:});
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/discretise/quant.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,8 @@ +function y=quant(q,x) +% quant - Quantize values +% +% quant :: real ~'quantum', real~'values to quantize' -> real. +% +% Note: this function uses ceil, not floor or round to quantize. + +y=q.*ceil(x./q);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/disposables.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,37 @@ +% disposables - Manage resources that need to be disposed of cleanly +% +% This is implemented as a message handler to be used as follows: +% +% >> ref=disposables('reg',x) % register x as a disposable returning reference ref +% >> disposables('dereg',ref) % deregister (but do not dispose) of registered object +% >> dx=disposables('get') % get cell array of all registered objects +% >> disposables('dispose') % dispose of and deregister all disposables +% +% In most cases, you will not need to use this function as resources should +% be disposed of correctly by the functions that acquire them. However, if +% things go wrong or the user presses Ctrl-C, some objects might not be +% disposed of properly and cause problems later. In this case, a call to +% disposables('dispose') should take care of it. + +function varargout=disposables(msg,varargin) + persistent registry + + switch(msg) + case 'reg', + obj=varargin{1}; + registry = [registry, obj]; + varargout{1}=obj; + case 'dereg', + obj=varargin{1}; + idx=find(registry==obj); + registry=registry([1:idx-1,idx+1:length(registry)]); + case 'get' + varargout{1}=registry; + case 'dispose' + for i=1:length(registry) + registry(i).dispose(); + end + registry=[]; + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/endswith.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,12 @@ +function b=endswith(str,suffix) +% endswith - test to see if a string ends with a certain suffix +% +% endswith :: string ~'string to check', string ~'suffix' -> bool. +% +% returns 1 iff end of string matches suffix (case insensitive) + +if length(str)<length(suffix) b=0; else +b=strcmpi(str(end-length(suffix)+1:end),suffix); +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/fileutils/frecurse.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,32 @@ +function frecurse(fn,dirname,level) +% recurse: recursively descend a directory tree +% +% recurse :: +% (string ~'directory name',natural~'level'=>bool~'descend into?') ~'fn to process directory, +% string ~'directory to start from, +% natural ~'initial level assigned to start directory' +% => void. +% +% usage: recurse(<function>,<dirname>[,<level0>]) +% +% traverses directory tree, starting at <dirname>, +% and calling the named <function> for each directory. +% <function> must be of the form: +% flag=<function>(<thisdirname>,<level>) +% <level> indicates how many directory levels down from +% the original directory <thisdirname> is. +% if <flag> is 1, recurse will descend into that directory + +if nargin<2, dirname=''; end +if nargin<3, level=0; end + +if fn(dirname,level) + list=dir(dirname); + for k=1:length(list) + name=[dirname '/' list(k).name]; + if list(k).isdir & list(k).name(1)~='.' + disp(['--- descending into ' name]); + recurse(fn,name,level+1); + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/fileutils/listfiles.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,47 @@ +function L=listfiles(dirname,varargin) +% listfiles - list files in a directory +% +% listfiles :: text ~'directory path', +% options { +% 'pattern':text ~'pattern to match files against', +% 'regext':text ~'regular expression to match files against', +% 'ext':text ~'file names must in end in this' +% 'rec':bool/0 ~'recursive search?' +% } +% -> {[N]->text} ~'cell array of N paths'. + +opts=prefs('rec',0,varargin{:}); + +if isfield(opts,'pattern'), + dirpat = fullfile(dirname,opts.pattern); +else + dirpat = dirname; +end + +if isfield(opts,'regexp'), + filter=@(f)(~isempty(regexp(f,opts.regexp))); +elseif isfield(opts,'ext'), + filter=@(f)endswith(f,opts.ext); +else + filter=@(f)true; +end + +D=dir(dirpat); +L={}; + +for i=1:length(D) + f=D(i); + if ~f.isdir && feval(filter,f.name), + L=[L {fullfile(dirname,f.name)}]; + end +end + +if opts.rec + if strcmp(dirname,dirpat), DX=D; else DX=dir(dirname); end + for i=1:length(DX) + f=DX(i); + if f.isdir && ~strcmp(f.name,'.') && ~strcmp(f.name,'..') + L=[L,listfiles(fullfile(dirname,f.name),opts)]; + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/fileutils/loadmat.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,31 @@ +function varargout=loadmat(file,varargin) +% loadmat - hygienic loading from mat files +% +% The point of this function is to load certain named variables +% from a Mat file without polluting the current environment with +% other objects from the Mat file or having name clashes. +% +% loadmat :: text ~'filename', text~'variable name' -> object~'something'. +% loadmat :: +% text ~'filename', +% {[N]->text} ~'list of variable names' +% -> {[N]->object} ~'list of things'. +% +% Vararg form (not well typed) +% +% [X,Y,Z,...]=loadmat(Filename,var1,var2,var,...); + +load(file); +for j=1:length(varargin) + var=varargin{j}; + if iscell(var), + X__=cell(1,length(var)); + for k=1:length(var) + eval(sprintf('X__{%d}=%s;',k,var{k})); + end + varargout{j}=X__; + else + eval(['varargout{j}=' var ';']); + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/fileutils/md5file.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +function hash=md5file(fn) + [rc,result]=system(sprintf('md5sum -b "%s"',fn)); + if rc==0 + sp=find(result==' '); + hash=result(1:sp); + else + error('Error running md5sum'); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/fileutils/music_files.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,15 @@ +% music_files - search for audio files under music directory +% +% music_files :: +% string ~'Path based at ~/Music' +% options { +% pattern :: string/'*' ~'filename matching patter'; +% rec :: bool/0 ~'recurse into subdirs?' +% } +% -> cells {string}. +% +% This function assumes that the HOME environment variable points to the +% user's home directory and that there exists there a directory called 'Music'. + +function L=music_files(d,varargin) + L=listfiles([getenv('HOME'),filesep,'Music',filesep,d],'regexp','.*\.(mp3|ogg|m4a|wav|au|flac)',varargin{:});
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/fileutils/read.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,32 @@ +function M=read(name) +% M=read(name): read and return ascii matrix file +% name may be with or without .txt extension + +if name(end-3)~='.' + name=[name '.txt']; +end + +load('-ascii',name); +last=max(findstr(name,'\')); +if ~isempty(last) + name=name(last+1:length(name)); +end +last=max(findstr(name,'/')); +if ~isempty(last) + name=name(last+1:length(name)); +end +if name(1)=='.' + name=name(2:length(name)); +end +if name(1)>='0' & name(1)<='9' + name=['X' name]; +end + +dot=max(findstr(name,'.')); +if ~isempty(dot) + name=name(1:dot-1); +end +rr=findstr(name,'.'); name(rr)=repmat('_',1,length(rr)); +rr=findstr(name,'&'); name(rr)=repmat('_',1,length(rr)); +eval(['M=' name ';']); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/fileutils/removeext.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,8 @@ +% removeext - remove extension from a file name +% +% removeext :: string -> string. +function nm=removeext(path) +dot=max(strfind(path,'.')); +if isempty(dot), nm=path; +else nm=path(1:(dot-1)); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/fileutils/tofile.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,24 @@ +function tofile(outfile,varargin) +% tofile - Run functions that need a file handle to write output +% +% tofile :: +% path ~'file name to write to', +% (fid -> action unit) ~'action taking file handle', +% (fid -> action unit) ~'action taking file handle', +% ... +% -> action unit. +% +% Opens the given file in APPEND mode, then calls each function +% passing the file handle to each. + + +fh=fopen(outfile,'a'); + +fprintf(fh,'\n%% --- New session ---\n\n'); + +for i=1:length(varargin) + fprintf(fh,'\n%% --- Running %s\n\n', func2str(varargin{i})); + feval(varargin{i},fh); +end + +fclose(fh);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/fileutils/write.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,8 @@ +function write(X,path) +% write - Write matrix as ascii text file +% +% write :: matrix, path -> action unit. + +save('-ascii','-double',path,'X'); + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/@func/func.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,208 @@ +classdef func + properties (GetAccess=private, SetAccess=private) + fn + args + slots + end + methods + % A func is a partially applied function: Any given function can + % have some arguments bound to some given parameters, producing a + % new function which needs only the remaining unbound parameters to + % be supplied. In addition, the order in which parameters supplied + % in a subsequent invocation are bound to unbound arguments can be + % controlled using a slot permutation mechanism. Subscripting operation + % () is overloaded so that function application works in the normal + % way except when there are no arguments, in which case you must use FEVAL. + % + % For example, if linspace is a function such that + % + % linspace(a,b,n) + % + % returns an array of n values linearly spaced between a and b, + % then after + % + % fromzero=func(@linspace,0) + % + % fromzero(b,n) returns n linearly spaced values from 0 to b. + % The arguments can be permuted, eg, + % + % perm(func(@linspace),[3 1 2]) + % + % returns a function that takes the arguments in the order (n,a,b) + % instead of (a,b,n). Permutation and binding can be combined, eg + % + % countdown=bind(perm(func(@linspace),[3 2 1]),10) + % countdown(5,50) + % ans= + % 50 45 40 35 30 25 20 15 10 5 + % + % See Also: + % PERM, BIND, CURRY, FEVAL, BINDAT + % + % Internal structure: + % cfn.fn : an ordinary matlab function handle + % cfn.args : a cell array with bound arguments filled in in the order + % that the ordinary function fn expects them. There may be + % empty slots if the arguments have been permuted (see CFPERM) + % cfn.slots : List of empty slots in the order that they are intended to be + % bound when this curried function is evaluated or has more + % parameters bound. Eg, the first parameter will be bound to + % to the slots(1)'th argument of the eventual call to cfn.fn + + % default function is null function + function c=func(fn,varargin) + + if nargin<1, fn=@nop; end + if isa(fn,'func'), c=fn; + elseif isa(fn,'struct'), + c.fn=fn.fn; + c.args=fn.args; + c.slots=fn.slots; + else + c.fn=fn; + c.args=varargin; + c.slots=[]; + end + end + + % bind - Bind arguments to a func to create partially applied function + % + % bind(fn,a1,a2,..) + % fn is curried function returned by FUNC, BIND or PERM, + % the given parameters are bound to the free slots of the function + % and a new curried function is returned. + + % Matlab 'feature': if ANY of the paremeters + % to bind are function objects, then this function gets called, even + % if the first parameter is just an ordinary function pointer. + function c=bind(c,varargin) + if isa(c,'function_handle'), c=func(c); end + + % first, fill available slots in arg list using supplied args + m=min(length(c.slots),length(varargin)); + c.args(c.slots(1:m))=varargin(1:m); + c.slots=c.slots(m+1:end); % remove filled slots + c.args=horzcat(c.args,varargin(m+1:end)); % append any left over args + end + + % BINDAT: Bind arguments at specified positions to create partially applied + % function + % + % BIND(cfn,pos,a1,a2,..) + % If cfn is an ordinary function, the parameters a1, a2 etc + % are bound at the positions specified in pos(1),pos(2), etc. + % Otherwise, the same as BIND. + % + % See BIND + + % need to make sure pos is a complete permutation before + % calling perm + function c=bindat(c,pos,varargin) + pos=[pos setdiff(1:max(pos),pos)]; + c=bind(perm(c,pos),varargin{:}); + end + + % DISPLAY - Display FUNC object as string + function display(c), disp([tostring(c),sprintf(' :: func(%d->%d)',nargin(c),nargout(c))]); end + + % FEVAL: Evaluate partially applied function + % + % FEVAL(cfn,a1,a2,...) + % The given parameters a1, a2, etc. are bound to the arguments of cfn, + % and eventually to the arguments of a matlab function, possibly + % permuted according to cfn.slots (see CFPERM) + function varargout=feval(c,varargin) + % first, fill available slots in arg list using supplied args + if ~isa(c,'func'), + [varargout{1:max(1,nargout)}] = builtin('feval',c,varargin{:}); + return; + % error('This shouldn''t be happening!'); + end + m=min(length(c.slots),length(varargin)); + args=c.args; + args(c.slots(1:m))=varargin(1:m); + + % append any left over args + if nargout==0, c.fn(args{:},varargin{m+1:end}); + else [varargout{1:max(1,nargout)}] = c.fn(args{:},varargin{m+1:end}); + end + end + + + % GETARG - Extract a bound argument from a FUNC object + % getarg :: (func (D1,...,Dn)->R, k:natural) -> Dk + function a=getarg(c,n), a=c.args{n}; end + + % GETFN - Get ordinary function handle from FUNC object + function fn=getfn(c), fn=c.fn; end + + % PERM: permute arguments to a curried function + % PERM(cfn,perm) + % cfn can be an ordinary function handle or a curried function as returned + % by BIND or PERM. + % + % The 1st argument of PERM(F,P) will bind to the P(1)'th argument of F. + % Sequential permutations compose properly. + % P must be a complete permutation: if length(P)=L, then P must + % contain each of the intergers 1 to P exactly once. + % The binding order of arguments after the Lth is not affected. + function c=perm(c,slots) + m=length(c.slots); + n=length(slots); + l=length(c.args); + if m<n, c.slots=[c.slots l+1:l+n-m]; + elseif m>n, slots=[slots n+1:m]; end + c.slots=c.slots(slots); + end + + function varargout=subsref(c,S) + % SUBSREF - Subscript referencing for FUNC objects + % + % f(a,b,c)=feval(f,a,b,c) - parenthesis subsref invokes the function + % f{a,b,c}=bind(f,a,b,c) - braces subsref binds more arguments + + + s=S(1); + switch s.type + case '()' + % first, fill available slots in arg list using supplied args + m=min(length(c.slots),length(s.subs)); + args=c.args; + args(c.slots(1:m))=s.subs(1:m); + [varargout{1:nargout}] = c.fn(args{:},s.subs{m+1:end}); + case '{}' + varargout{1}=bind(c,s.subs{:}); + if length(S)>1, + [varargout{1:nargout}] = subsref(varargout{1},S(2:end)); + end + end + end + + function str=tostring(c) + % TOSTRING - Implementation of string conversion for FUNC objects + % + % tostring :: func -> string + + str=[]; + slots=c.slots; m=length(slots); + args=c.args; l=length(args); + n=max(slots); + if l<n, args=[args cell(1,n-l)]; l=n; end; + for i=1:m + args{slots(i)}=['#' num2str(i)]; + end + + for i=1:l + str=[str tostring(args{i}) ',']; + end + if n>0, str=[str ' ']; end + if isa(c.fn,'function_handle') fstr=['@' func2str(c.fn)]; + elseif isa(c.fn,'inline') fstr=['\' char(c.fn) '\']; + end + str=[ fstr '(' str '...)']; + end + function n=nargout(f), n=nargout(f.fn); end + function n=nargin(f), n=nargin(f.fn)-length(f.args); end + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/@thunk/thunk.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,33 @@ +% thunk(B:arglist(M)) ~ a thunk is something that eventually yields M outputs +classdef thunk + properties (GetAccess=private, SetAccess=immutable) + head % :: (A:arglist(N) -> B | thunk(B)) + args % :: A:arglist(N) + end + methods + % thunk - create delayed evaluation + function f=thunk(h,varargin), f.head=h; f.args=varargin; end + function s=tostring(f), s=sprintf('%s(%s)',tostring(f.head),commalist(map(@tostring,f.args))); end + function display(f), display([tostring(f),' :: thunk']); end + + % fforce - Completely force a thunk, if keep evaluating till not a thunk + function f=fforce(f), while isa(f,'thunk'), f=force(f); end; end + + function varargout=force(f), [varargout{1:nargout}] = feval(f.head,f.args{:}); end + + % SUBSREF - Subscript referencing for FUNC objects + function varargout=subsref(f,S) + s=S(1); + switch s.type + case '()', [varargout{1:nargout}] = feval(f.head,f.args{:}); + end + end + end +end + +function s=commalist(a) + if isempty(a), s=''; + elseif length(a)==1, s=a{1}; + else s=[a{1},',',commalist(a(2:end))]; + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/apply.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +function varargout=apply(fn,args) +% apply - Apply function to list (cell array) of arguments +% +% Not typable. +% +% if [y1,y2,...]=fn(x1,x2,...), then +% [y1,y2,...]=apply(fn,{x1,x2,...}) + +[varargout{1:max(1,nargout)}] = feval(fn,args{:});
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/bind.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,15 @@ +function cfn=bind(fn,varargin) +% bind - Create a partially applied function +% +% BIND(fn,a1,a2,...) +% If fn is an ordinary function, the parameters a1, a2 etc +% are bound as the first few arguments +% +% If fn is a function object as returned by FUNC, BIND, or PERM, +% the given parameters are bound to the free slots of the function +% and a new partially applied function is returned. +% + +if ischar(fn), fn=str2func(fn); end +if ~isa(fn,'func'), fn=func(fn); end +cfn=bind(fn,varargin{:});
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/bind1.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +% bind1 - Bind arguments to a function using Matlab closure +% +% bind1 :: (A1, A2, ... -> B1, B2 ..), A1 -> (A2, ...) -> B1, B2 ... + +function g=bind1(f,varargin) + args=varargin; + g=@(varargin)f(args{:},varargin{:}); +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/bindat.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,19 @@ +function cfn=bindat(fn,varargin) +% bindat - Create function with bound arguments at given positions +% +% BINDAT(fn,positions,arguments) +% If fn is an ordinary function handle, the parameters a1, a2 etc +% are bound at the positions given +% +% If fn is a closure as returned by FUNC, BIND, or PERM, +% the given parameters are bound to the specified slots of the function. +% +% Eg, a function which divides a value by 4 can be given as +% +% quarter = bindat(@rdivide,2,2); +% quarter(8) +% ans = 2 + +if ischar(fn), fn=str2func(fn); end +if ~isa(fn,'func'), fn=func(fn); end +cfn=bindat(fn,varargin{:});
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/compose.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,10 @@ +function h=compose(f,g) +% compose - Constructs composition of two functions +% +% compose :: (B->C), (A->B) -> (A->C). +% +% returns the function h such that h(...) = f(g(...)) + + h=@fg; + function x=fg(varargin), x=f(g(varargin{:})); end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/constfn.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,8 @@ +function f=constfn(a) +% constfn - returns a function which always returns the same thing. +% +% constfn :: A -> (_ -> A). + f=@const; + function x=const(varargin), x=a; end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/doer.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,22 @@ +function h=funseq(varargin) +% funseq - Constructs sequential application of several functions +% +% h=funseq(f,...,g) +% +% returns the function h such that h(...) = g(...) but +% all functions are evaluated in order. + + fns=varargin; + h=@seq; + + function varargout=seq(varargin) + for i=1:length(fns)-1 + feval(fns{i},varargin{:}); + end + if nargout==0, + feval(fns{end},varargin{:}); + else + [varargout{1:nargout}]=feval(fns{end},varargin{:}); + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/doreturn.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +% doreturn - Do an action and then return some other value +% +% doreturn :: (void=>A), B => B. +% doreturn :: (void=>A), B1, ..., BN => B1, ..., BN. +% +% This is useful when you need to construct an action expression +% using sub actions that return no value. +function varargout=doreturn(a,varargin), a(); varargout=varargin; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/feval2cell.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,13 @@ +function Y=feval2cell(I,F,varargin) +% feval2cell - Evaluate function and return multiple values in cell array +% +% Usage: +% +% if F is a function that returns M >= max(I) values, then +% Y = feval2cell(I,F,...) +% evaluates [Z{1:max(I)}]=F(...) +% and returns Z(I). + +Z=cell(max(I),1); +[Z{:}]=F(varargin{:}); +Y=Z(I);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/flip.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,13 @@ +function g=flip(f) +% flip - Flip order of first two arguments to a function +% +% flip :: (A,B->C) -> (B,A->C). + +if isa(f,'func'), + g=perm(f,[2,1]); +else + if nargin(f)==2, f=@(a,b)f(b,a); + else f=@(a,b,varargin)f(b,a,varargin{:}); + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/fold.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,13 @@ +function X=fold(fn,args) +% fold - Fold combinator from functional programming +% +% This function applies an associative operator to a list of arguments, +% as if the list was written out with the operator between consecutive +% elements, eg fold(@plus,{1,2,3,4}) = 1+2+3+4. +% +% fold :: +% (X,X->X) ~'associative binary operator', +% {[N]->X} ~'list (cell array) of arguments (at least 2 elements)' +% -> X. + +if isempty(args), X=[]; else X=foldl(fn,args{1},args(2:end)); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/foldcols.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,24 @@ +% foldcols - fold from left combinator for columns of an array +% +% foldcols :; +% (A,[[N]]->A) ~'folding function', +% A ~'initial value', +% [[N,L]] ~'data to scan, sequence of length L' +% -> A. + +function y=foldcols(f,y,X,varargin) + if nargin>3 + opts=prefs('draw',0,varargin{:}); + for i=1:size(X,2) + y1=f(y,X(:,i)); + if opts.draw + opts.plotfn(i,y,X(:,i),y1); + end + optpause(opts); + y=y1; + end + else % streamlined version + for i=1:size(X,2), y=f(y,X(:,i)); end + end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/foldl.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,13 @@ +function X=foldl(fn,e,args) +% foldl - Fold fom the left combinator +% +% This function applies an associative operator to cell array, +% starting from the left using the given starting element. +% eg fold(@plus,0,{1,2,3,4}) = ((((0+1)+2)+3)+4). +% +% foldl :: +% (X,Y->X) ~'associative binary operator', +% X ~'initial element', +% {[N]->Y} ~'list (cell array) of arguments (at least 2 elements)' +% -> X. +X=e; for i=1:length(args), X=fn(X,args{i}); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/foldr.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,13 @@ +function X=foldr(fn,e,args) +% foldr - Fold combinator from functional programming +% +% This function applies an associative operator to a list of arguments, +% starting from the right, eg +% elements, eg foldr(@plus,0,{1,2,3,4}) = (1+(2+(3+(4+0)))). +% +% foldr :: +% (X,X->X) ~'associative binary operator', +% X ~'initial element', +% {[N]->X} ~'list (cell array) of arguments (at least 2 elements)' +% -> X. +X=e; for i=length(args):-1:1, X=fn(args{i},X); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/foldrcols.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,22 @@ +% foldcols :; +% ([[M]], [[N]] -> [[M]]) ~'folding function', +% [[M]] ~'initial value', +% [[N,L]] ~'data to scan, sequence of length L' +% -> [[M]]. + +function y=foldrcols(f,y,X,varargin) + if nargin>3 + opts=prefs('draw',0,varargin{:}); + for i=size(X,2):-1:1 + y1=f(y,X(:,i)); + if opts.draw + opts.plotfn(i,y,X(:,i),y1); + end + optpause(opts); + y=y1; + end + else % streamlined version + for i=size(X,2):-1:1, y=f(y,X(:,i)); end + end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/foreach.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,15 @@ +function foreach(f,X,varargin) +% foreach - do an action for each element in a cell array in order +% +% foreach :: (A->action), {[N]->A} -> action. +% foreach :: (A->action), {[N]->A}, options {} -> action. +% +% Takes the same options as iterate. + + iterate(@g,X,varargin{:}); + + function x=g(x) + f(x{1}); + x=x(2:end); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/forels.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,23 @@ +function forels(f,X,varargin) +% forels - do an action for each element of an array in order +% +% forels :: +% (A->action) ~'action to apply to each element', +% [[Size]->A] ~'array of elements of type', +% options { +% pause :: bool/0 +% drawnow :: bool/1 +% } ~'see ITERATE for more options' +% -> action. +% +% Note, the array can be multidimensional - the standard order +% cycles through the earlier indices before the later ones, eg +% rows, then columns, then slices etc. + + N=numel(X); + iterate(@g,1,varargin{:}); + + function i=g(i) + if i>N, i=[]; else f(X(i)); i=i+1; end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/groupby.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,21 @@ +function y=groupby_or(f,x) +% groupby_or - collect rows of x into equivalence classes based f +% +% groupby_or :: ([[1,M]]->A), [[N,M]] -> {[K]->[[L,M]]} :- equality(A). +% +% IMPORTANT: this version assumes that data is ordered by f, that is, +% only consecutive rows will be grouped together. + + N=size(x,1); + y={}; + + i=1; j=1; + while j<=N + if i==j, z=f(x(i,:)); end + if j==N || f(x(j+1,:)~=z) + y=vertcat(y,{x(i:j,:)}); + i=j+1; j=i; + else j=j+1; + end + end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/id.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,5 @@ +function y=id(x) +% id - identity function +% +% id :: A->A. +y=x;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/ifx.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,4 @@ +% ifx - 'if' expression +% +% ifx :: bool, A, A -> A. +function v=ifx(f,a,b), if f, v=a; else v=b; end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/lazy.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,27 @@ +% lazy - Lazy evaluator, or memoiser +% +% lazy :: (unit -> A) -> (unit -> A). +% lazy :: (unit -> action A) -> (unit -> action A). +% +% If x=lazy(f), then f is not called until x is +% called. Then x() returns x1=f(). Thereafter, +% x() returns x1 and f is never called again. +% +% Eg, let x=lazy(@()rand). +% Then x() returns the same random number each +% time it is called. +function f=lazy(thunk) + value=thunk; + evalled=false; + f=@get; + + function v=get, + if evalled, + v=value; + else + value=thunk(); + evalled=true; + v=value; + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/map.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,23 @@ +function y=map(f,x,dim) +% map - Map function over collection +% +% map :: (A->B), collection(A) -> collection(B). +% +% collection(A) ::= seq A | {D->A} | [D->A]. +% +% This implementation handles cell arrays and ordinary arrays. +% Sequences should be taken care of by data/map + +if isempty(x) + y=x; +elseif iscell(x), + y=cellmap(f,x); +elseif isnumeric(x) + if nargin<3, y=reshape(f(x),size(x)); + elseif dim==1, y=maprows(f,x); + elseif dim==2, y=mapcols(f,x); + else error(sprintf('Cannot map over dimension %d',dim)); + end +else + error(sprintf('Cannot map over objects of class %s',class(x))); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/mapcols.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,16 @@ +function Y=mapcols(f,X) +% mapcols - Map a function of a vector over the columns of an array +% +% mapcols :: +% ([[N]->A] -> [[M]->B]) ~'function maps a column of A to one of B', +% [[N,L]->A] +% -> [[M,L]->B]. + +n=size(X,2); +if n==0, Y=[]; +else + for i=n:-1:1 + Y(:,i)=flatten(feval(f,flatten(X(:,i)))); + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/maprows.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,16 @@ +function Y=maprows(f,X) +% maprows - Map a function of a vector over the rows of an array +% +% maprows :: +% ([[1,N]->A] -> [[1,M]->B]) ~'function maps a row of A to a row of B', +% [[L,N]->A] +% -> [[L,M]->B]. + +n=size(X,1); +if n==0, Y=[]; +else + for i=n:-1:1 + Y(i,:)=f(X(i,:)); + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/multicall.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,17 @@ +function varargout=multicall(varargin) +% multicall - sequential call to several closures, return values from last +% +% y=do(f,g,...,h,x) +% +% equivalent to +% f(x); g(x); ...; y=h(x); + + for i=1:length(varargin)-1 + feval(varargin{i}); + end + if nargout==0, + feval(varargin{end}); + else + [varargout{1:nargout}]=feval(varargin{end}); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/nop.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,3 @@ +function x=nop(varargin) +% nop - Function that takes any number of args and does nothing +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/nthret.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,13 @@ +function X=nthret(I,F,varargin) +% nthret - Evaluate function and return Nth return value. +% +% Usage: +% +% if F is a function that returns M >= I values, then +% X = nthret(N,F,...) +% evaluates [Z{1:N}]=F(...) +% and returns Z{N}. + +Z=cell(I,1); +[Z{:}]=F(varargin{:}); +X=Z{I};
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/scancols.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,23 @@ +% scancols :; +% ([[M]], [[N]] -> [[M]]) ~'scannning function', +% [[M]] ~'initial value', +% [[N,L]] ~'data to scan, sequence of length L' +% -> [[M,L]]. + +function Y=scancols(f,y,X,varargin) + Y=zeros(size(y,1),size(X,2)); + + if nargin>3 + opts=prefs('draw',0,varargin{:}); + for i=1:size(X,2) + y1=f(y,X(:,i)); + Y(:,i)=y1; + if opts.draw, opts.plotfn(i,y,X(:,i),y1); end + optpause(opts); + y=y1; + end + else + for i=1:size(X,2), y=f(y,X(:,i)); Y(:,i)=y; end + end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/scanl.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,24 @@ +% scanl - scanl for cell arrays +% +% scanl :: +% (S,X->S) ~'scannning function', +% S ~'initial value', +% {[N]->X} ~'data to scan, sequence of length L' +% -> {[N]->S}. + +function Y=scanl(f,y,X,varargin) + Y=cell(size(X)); + if nargin>3 + opts=prefs('draw',0,varargin{:}); + for i=1:size(X,2) + y1=f(y,X{i}); + Y{i}=y1; + if opts.draw, opts.plotfn(i,y,X{i},y1); end + optpause(opts); + y=y1; + end + else + for i=1:size(X,2), y=f(y,X{i}); Y{i}=y; end + end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/scanrcols.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,25 @@ +% scanrcols - scan array columns from the right +% +% scanrcols :; +% ([[M]], [[N]] -> [[M]]) ~'scannning function', +% [[M]] ~'initial value', +% [[N,L]] ~'data to scan, sequence of length L' +% -> [[M,L]]. + +function Y=scanrcols(f,y,X,varargin) + Y=zeros(size(y,1),size(X,2)); + + if nargin>3 + opts=prefs('draw',0,varargin{:}); + for i=size(X,2):-1:1 + y1=f(y,X(:,i)); + Y(:,i)=y1; + if opts.draw, opts.plotfn(i,y,X(:,i),y1); end + optpause(opts); + y=y1; + end + else + for i=size(X,2):-1:1, y=f(y,X(:,i)); Y(:,i)=y; end + end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/tail.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,20 @@ +% tail1 - Call function with tail recursion +% +% tail1 :: (A{1:N}->thunk(B{1:M}), A{1:N} -> B{1:M}. +% +% thunk(B{1:M}) ::= +% exists C{1:L} . cell { C{1:L}->B{1:M}, cell { C{1:L} } } +% | cell { {0}, cell { B{1:L} } } +% | unit -> B{1:M}. + +function varargout=feval_tail1(fn,varargin) + ret=fn(varargin{:}); + keyboard + if ~isnumeric(ret{1}) + [varargout{1:nargout}]=feval(ret{1},ret{2}{:}); + else + varargout=ret{2}(1:nargout); + end +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/tail_call.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +% tail_call - Simulate effect of tail call +% +% tail_call :: (A{1:N}->B{1:M}), A{1:N} -> thunk(B{1:M}). +% +% thunk(B{1:M}) ::= +% exists C{1:L} . cell { C{1:L}->B{1:M}, cell { C{1:L} } } +% | cell { {0}, cell { B{1:L} } } +% | unit -> B{1:M}. +function ret=tail_call(fn,varargin), ret={fn,varargin}; end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/tail_eval.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,11 @@ +% tail_eval - Evaluate thunk +% +% tail_eval :: thunk(B{1:M}) -> B{1:M}. +function varargout=tail_eval(thunk) + if isnumeric(thunk{1}), varargout=thunk{2}; + else + [varargout{1:nargout}]=feval(thunk{1},thunk{2}{:}); + end +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/tail_return.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +% tail_return - Package return values for tail_feval +% +% tail_return :: A{1:N} -> thunk(A{1:N}). +% +% thunk(B{1:M}) ::= +% exists C{1:L} . cell { C{1:L}->B{1:M}, cell { C{1:L} } } +% | cell { {0}, cell { B{1:L} } } +% | unit -> B{1:M}. +function ret=tail_return(varargin), ret={0,varargin}; end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/tailrec.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,16 @@ +% tail_feval - Call function with tail recursion +% +% tail_feval :: (A{1:N}->rthunk(B{1:M}), A{1:N} -> B{1:M}. +% +% rthunk(B{1:M}) ::= +% exists C{1:L} . cell { C{1:L} -> rthunk(B{1::}), C{1:L} }. +% | cell { {0}, cell { B{1:L} } } +function varargout=tail_feval(fn,varargin) + ret=fn(varargin{:}); + while ~isnumeric(ret{1}) + ret=feval(ret{1},ret{2}{:}); + end + varargout=ret{2}(1:nargout); +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/tbind.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,8 @@ + +% tbind - bind args to tail recursive call +% +function g=tbind(varargin), + args=varargin; + g=@(varargin)tailrec(args{:},varargin{:}); +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/zip.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,4 @@ +% zip - Polymorphic zip using zipwith(@tuple,...) +function y=zip(varargin) + y=zipwith(@tuple,varargin{:}); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/zipmap.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,139 @@ +function varargout=zipmap(spec,f,varargin) +% zipmap - Map and zip function over multidimensional arrays +% +% ZIPMAP allows a certain simple functions to be applied to multidimensional +% array arguments, manipulating the dimensions to combine elementwise and +% outer-product-like behaviour. (The name comes from functional programming +% terminology: functions/operators like .* and + which apply elementwise to +% vector arguments are said to zip their arguments, while outer products are +% like mapping a function over a list, or a list of functions over another list.) +% +% Let f be a vectorised function which knows how to zip array arguments. Eg, +% f(x,y)=x+y. If x and y are vectors, they must be the same size, and +% the function is applied to corresponding pairs of elements, resulting in +% a vector of the same size. With ZIPMAP, you can build functions, that, eg, +% add EVERY pair of elements (not just corresponding pairs) resulting in a +% 2D array: zipmap(0,inline('x+y','x','y'),[1:4]',[1:6]') will do this +% +% The domain of a multidimensional array can be described a list containing the +% size of each dimension, as returned by the SIZE function (or SIZE1 which strips +% trailing 1s). The arguments to f must have a certain domain structure: +% +% dom arg1 = [p1 zipdom] +% dom arg2 = [p2 zipdom] +% : +% +% where p1, p2, etc are the domains of the elementary function of which f +% is a vectorisation, eg the elementary domain of .* can be scalars and +% hence p1=p2=[], but the two elementary domains of CROSS are both [3], because it +% it operates on pairs of 3-D vectors. The zipdom part must be the same for +% all arguments, otherwise, we cannot extract matching argument tuples. +% If this is the case, the f must return a result: +% +% dom f(...) : [rdom zipdom] +% +% where rdom is the domain of the elementary function, eg [3] for the vector +% valued CROSS product, but [] for the scalar valued DOT product. Note that +% these domains can be vectors of any length from zero up, eg the determinant +% of a square matrix has a domain of [n n]. +% +% ZIPDIM allows the domains of the arguments and result to be generalised as follows: +% +% dom arg1 : [p1 m1 zipdom] +% dom arg2 : [p2 m2 zipdom] +% : +% +% dom zipdom(...) : [rdom m1 m2 ... zipdom ] +% +% All the mapped dimensions m1, m2 etc go into a sort of outer product where +% every combination of values is evaluated. In addition, the order of the +% m1, m2 etc can be permuted in the result at essentially no extra computational cost. +% +% Usage: r=zipmap(spec,f,<varargs>) +% +% spec defines how the domain of each of the arguments is to be interpreted, ie +% how is is split into elementary, mappable, and zippable domains. +% +% If spec is a number z, then the last z dimensions of each argument are marked for +% zipping, while the rest are for mapping. The elementary domains are assumed to +% be scalar and hence p1, p2 etc =[]. So, zipmap(0,...) does the complete outer +% product type of evaluation. +% +% If spec is structure, the field 'zipdims' determines how many dimensions are zipped. +% The field 'pdims' is a the vector [length(p1) length(p2) ..], ie the dimensionality +% of each elementary domain: 0 for scalars, 1 for vectors, 2 for matrices etc. +% The optional field 'mapord' determines the reordering of the mapped domains, eg +% if mapord=[2 1 3], then the result domain is [m2 m1 m3] instead of [m1 m2 m3]. +% +% If spec is a cell array, it is interpreted as {zipdims pdims mapord}. +% +% Extra options +% +% If spec is a structure, it can contain further options: if spec.shift=1, +% then the domains of all arguments are pre-shifted to remove singleton +% dimensions before anything else is done. Amongst other effects, this means +% that row vectors can be used in place of column vectors. + + +% unpack args from spec +if iscell(spec) + [zipdims,pdims,mapord]=getargs(spec,{0 [] []}); +elseif isstruct(spec) + zipdims=getparam(spec,'zipdims',0); + pdims=getparam(spec,'pdims',[]); + mapord=getparam(spec,'mapord',[]); + if getparam(spec,'shift',0), + varargin=cellmap(@shiftdim,varargin); + end +elseif isscalar(spec) + zipdims=spec; pdims=[]; mapord=[]; +end + + +argdoms=cellmap(@size1,varargin); +if isempty(pdims), pdims=zeros(size(varargin)); end +for i=1:length(varargin) + [pdoms{i},mapdoms{i},zipdoms{i}]=split(argdoms{i},pdims(i),zipdims); +end + + +% check all zipdoms match and use first +zipdom=zipdoms{1}; +for i=2:length(zipdoms) + if ~all(zipdoms{i}==zipdom), error('Zip domains do not match'); end +end + + +if isempty(mapord), mapord=1:length(varargin); end +outerdom=cat(2,mapdoms{mapord}); +outer1=cellmap(@to1s,mapdoms); +for i=1:length(varargin) + % need to reshape each arg and then repmat up to correct size + newsizes=outer1; newsizes{i}=mapdoms{i}; + newsize=[pdoms{i} cat(2,newsizes{mapord}) zipdom]; + varargin{i}=reshape(varargin{i},tosize([pdoms{i} cat(2,newsizes{mapord}) zipdom])); + newsizes=mapdoms; newsizes{i}=outer1{i}; + newsize=[to1s(pdoms{i}) cat(2,newsizes{mapord}) to1s(zipdom)]; + varargin{i}=repmat(varargin{i},tosize(newsize)); +% size(varargin{i}) +end + +[varargout{1:max(1,nargout)}]=feval(f,varargin{:}); + +function [a,b,c]=split(x,n,m) + a=x(1:n); + b=x(n+1:end-m); + c=x(end-m+1:end); + +% getargs - get values of optional parameters with defaults +% +% getargs :: {I:[N]->(A(I) | nil)}, {I:[N]->A(I)} -> A(1), ..., A(N). + +function varargout=getargs(args,defs) + +nout=max(length(args),length(defs)); +varargout=defs; +for i=1:length(args) + if ~isempty(args{i}), varargout(i)=args(i); end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/funutils/zipwith.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,2 @@ +% zipwith - zipwith for cell arrays only +function y=zipwith(f,varargin), y=cellzip(f,varargin{:}); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/general.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,34 @@ +% general - Add subdirectories of general to path +% +% general :: cells { string } => void. +% +% Use general('all') or +% >> general all +% To add all subdirs to path. +% +% Available subdirectories are: +% algo - Algorithms and APIs +% arrutils - General array manipulations +% cellutils - Cell manipulations +% fileutils - File handling +% funutils - Functional programming +% numerical - Numerical analysis (inc scalar, vector, matrix) +% discretise - Discretisation and quantisation tools + +function general(dirs) + if nargin<1, + disp('Type >> general all to add all subdirectories to matlab path'); + return; + end + + if strcmp(dirs,'all') + dirs={'algo','arrutils','cellutils','fileutils','funutils','numerical','discretise'}; + end + thisfile=which('general'); + seps=strfind(thisfile,filesep); + thisdir=thisfile(1:seps(end)); + for i=1:length(dirs) + addpath([thisdir,dirs{i}]); + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/getparam.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +function a=getparam(s,field,def) +% getparam - get value of named field or default if not present +% +% getparam :: struct Fields, string, A -> A. +% +% Typical usage is with structures returned by prefs. + +if isfield(s,field), a=s.(field); else a=def; end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/addnorm.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,23 @@ +function [Y,M]=addnorm(F,X) +% addnorm - additive normalisation with respect to arbitrary function +% +% addnorm :: +% ([[N,M]] -> [[N,1]]) | ([[N,M]] -> [[1,M]]) +% ~'function to compute offsets', +% [[N,M]] ~'data to normalise' +% -> [[N,M]] ~'offset vectors' +% [[N,1]] | [[1,M]] ~'the vector that was subtracted'. +% +% This function can work row-wise OR column-wise depending on what +% the offset computing function returns. +% +% EXAMPLES +% +% To zero-mean each ROW of an N-by-M array X, use +% addnorm(@(t)mean(t,2),X); +% +% To zero the MEDIAN of each COLUMN, use +% addnorm(@median,X); + +Y=vecop1(@minus,X,F(X)); % vecshift(F(X),X); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/alpha_norm.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,12 @@ +function l=alpha_norm(alpha,x) +% alpha_norm - Minkowski norm (actually more like mean) using given exponent +% +% alpha_norm :: nonneg ~'exponent', [[N,M]] -> [[1,M]->nonneg]. +% alpha_norm :: nonneg ~'exponent', [[1,N]] -> nonneg. +% +% This version divides by the dimemnsionality of the space, +% eg +% alpha_norm(2,x) = sqrt(mean(x.^2,1)) +% alpha_norm(1,x) = mean(abs(x),1)) +l=real(mean(abs(x).^alpha).^(1/alpha)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/argmax.m Sat Jan 12 19:21:22 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(:));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/argmax1.m Sat Jan 12 19:21:22 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); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/argmaxv.m Sat Jan 12 19:21:22 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{:}];
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/argmin1.m Sat Jan 12 19:21:22 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); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/arrshift.m Sat Jan 12 19:21:22 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)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/asym.m Sat Jan 12 19:21:22 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]].
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/binary.m Sat Jan 12 19:21:22 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 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/bitmap.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/cauchy_norm.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,21 @@ +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); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/circulant.m Sat Jan 12 19:21:22 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 + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/clamp.m Sat Jan 12 19:21:22 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(x<a)=a; +x(x>b)=b; +%y=min(b,max(a,x));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/closed2hopen.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,12 @@ +function x=closed2hopen(x) +% closed2hopen - convert closed interval to half-open interval +% +% closed2hopen :: [[N,2]] -> [[N,2]] +% +% The second column is assumed to be the upper limit of a closed +% interval. These upper limits are incremented to the next floating +% point number so that the result can be interpreted as a half-open +% interval equivalent to the original closed interval. + +mx=x(:,2); +x(:,2)=mx+eps(mx);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/complexdom.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,4 @@ +function Z=complexdom(reals,imags) + +[X,Y]=meshgrid(reals,imags); +Z=X+i*Y;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/deta.m Sat Jan 12 19:21:22 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 +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/diffdims.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/divnorm.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,15 @@ +function [Y,Z]=divnorm(F,X) +% divnorm - Do divisive normalisation of a sequence of vectors. +% +% divnorm :: +% ([[N,L]] -> [[1,L]]) ~'function to compute norm', +% [[N,L]] ~'array of L N-dim vectors' +% -> [[N,L]] ~'normalised vectors (ie unit 2-norm)', +% [[1,L]] ~'array of normalisation factors'. +% +% ie, if [Y,Z]=divnorm(X), +% then X = repmat(Z,size(X,1),1).*Y +% +% TODO: might like to smooth the envelope before dividing. +Y = vecop1(@rdivide,X,F(X)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/eigsel.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/eigvecs.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,5 @@ +function V=eigvecs(A,I) + [V,D]=eig(A); + V=fliplr(V); + V=V(:,I); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/erfp.m Sat Jan 12 19:21:22 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;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/erfratio.m Sat Jan 12 19:21:22 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)));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/expdecay.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/fixpoint.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,26 @@ +function y=fixpoint(f,x0,varargin) +% fixpoint - Find the fixed point of a function by repeated evaluation +% +% fixpoint :: +% (A->A) ~'function to iterate', +% A ~'initial value', +% options { +% its :: natural/inf ~'maximimum numbert of iterations'; +% testfn :: (A,A->boolean) ~'convergence test' +% } +% -> maybe A ~'final value or nan if none'. + + +opts=prefs('its',inf,'testfn',@epseq,varargin{:}); + +conv=opts.testfn; +its =opts.its; +i=0; while i<its + x=f(x0); + if conv(x0,x), break; end + x0=x; i=i+1; +end + +if i<its, y=x; else y=nan; end + +function b=epseq(x,y), b=abs(y-x)<=eps(x);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/frint.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,22 @@ +function y=frint(r,x,dc) +% frint - Fractional integration +% +% y=frint(r,x[,dc]) +% Returns the rth order integral of x, where r can be fractional +% if dc=0 (default is 1) dc term is not added back in + +x=x(:); +n=length(x); +z=fft(x); z0=z(1); z(1)=0; +K=(0:n-1)'; +D=exp(2*pi*i*K/n)-1; +D(1)=1; % to avoid divide by zero +z2=z./D.^r; +y=real(ifft(z2)); % integral minus DC term + +% put DC term back as a ramp +if nargin<3 || dc==1, + y=y+sign(z0)*(K*abs(z0)/n).^r; +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/fromdbs.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,7 @@ +function X=fromdbs(Y) +% fromdbs - Convert from decibels +% +% fromdbs :: real -> real. + +X=10.^(Y/10); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/fuzzeye.m Sat Jan 12 19:21:22 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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/gauss_fn.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,4 @@ +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'.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/gauss_logpdf.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +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));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/gauss_win.m Sat Jan 12 19:21:22 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=gauss_fn(2*s/(sigma*n));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/gausspdf.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +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));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/hgauss_win.m Sat Jan 12 19:21:22 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=gauss_fn((0:n)'/(sigma*n));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/hshrink.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,7 @@ +% shrink - Shrinkage towards zero +% +% shrink :: real, [Size] -> [Size]. +function y=shrink(th,x) + +y=x; +y(y<th)=0;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/hsphere_surf.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,5 @@ +function S=hsphere_surf(N), S=(2*pi.^(N/2))./gamma(N/2); end +% hsphere_surf - Surface area of unit hypersphere +% +% hsphere_surf :: natural -> nonneg. +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/hsphere_surfln.m Sat Jan 12 19:21:22 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)) +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/hypdecay.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/hzeta.c Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,96 @@ +/*================================================================= + * hzeta.c + * + * Wrapper function for gsl_sf_hzeta_e + * + * This MEX function computes the Hurwitz-zeta-function by calling + * the function gsl_sf_hzeta_e from the GNU Scientific Library. For + * x>1 and a>=0 the Hurwitz-zeta-function is defined as + * + * oo + * ----- + * \ 1 + * hzeta(x, a) = ) -------- . + * / x + * ----- (i + a) + * i = 0 + * + * The MEX function takes two arguments, x and a. If both arguments + * are matrices of the same size, then the MEX function will compute + * the Hurwitz-zeta-function elementwise and the result of the MEX + * function will be a matrix of the same size as the inputs. If one + * input argument is a scalar s, then the function behaves as if this + * argument would be a matrix filled with the value s of the same + * size as the other argument. + * + * Compile this MEX function with "mex hzeta.c -lgsl". See Matlab + * documentation for additional information on compiling MEX functions. + * + * (C) by Heiko Bauke, 2007, heiko.bauke@physics.ox.ac.uk + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License in + * version 2 as published by the Free Software Foundation. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA + * 02111-1307, USA. + * + *=================================================================*/ + +#include "mex.h" +#include <gsl/gsl_errno.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_sf_zeta.h> + +void mexFunction(int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) { + double *x, *y, *a; + int nrows_x, ncols_x, i_x, + nrows_a, ncols_a, i_a, + nrows_y, ncols_y, i_y, + status; + gsl_sf_result result; + + /* check input arguments */ + if (nrhs!=2) + mexErrMsgTxt("Two input arguments required.\n"); + if (!mxIsDouble(prhs[0])) + mexErrMsgTxt("1st argument not a matrix of double."); + nrows_x=mxGetM(prhs[0]); + ncols_x=mxGetN(prhs[0]); + if (!mxIsDouble(prhs[1])) + mexErrMsgTxt("2nd argument not a matrix of double."); + nrows_a=mxGetM(prhs[1]); + ncols_a=mxGetN(prhs[1]); + if (!( (nrows_a==nrows_x && ncols_a==ncols_x) || + (nrows_x==1 && ncols_x==1) || + (nrows_a==1 && ncols_a==1) ) ) + mexErrMsgTxt("arguments must be matrices of the same size or scalars."); + /* create matrix for the return argument */ + nrows_y=nrows_x>nrows_a ? nrows_x : nrows_a; + ncols_y=ncols_x>ncols_a ? ncols_x : ncols_a; + plhs[0]=mxCreateDoubleMatrix(nrows_y, ncols_y, mxREAL); + /* assign pointers to each input and output */ + x=mxGetPr(prhs[0]); + a=mxGetPr(prhs[1]); + y=mxGetPr(plhs[0]); + /* compute zeta function */ + gsl_set_error_handler_off(); + for (i_x=i_a=i_y=0; i_y<nrows_y*ncols_y; ++i_y) { + if (gsl_sf_hzeta_e(x[i_x], a[i_a], &result)==GSL_SUCCESS) + y[i_y]=result.val; + else + y[i_y]=GSL_NAN; + if (!(nrows_x==1 && ncols_x==1)) + ++i_x; + if (!(nrows_a==1 && ncols_a==1)) + ++i_a; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/incr.m Sat Jan 12 19:21:22 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;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/inf2nan.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,6 @@ +function x=inf2nan(x) +% inf2nan - convert infs to nans +% +% inf2nan :: [Size] -> [Size]. + +x(isinf(x))=nan;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/inv_triu.m Sat Jan 12 19:21:22 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 + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/iqform.m Sat Jan 12 19:21:22 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]].
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/iqforma.m Sat Jan 12 19:21:22 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]].
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/iseven.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,2 @@ +% iseven :: integer -> bool. +function b=iseven(n), b=~mod(n,2);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/isodd.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,2 @@ +% isodd :: integer -> bool. +function b=isodd(n), b=mod(n,2);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/issym.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/lambertw.m Sat Jan 12 19:21:22 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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/linterp.m Sat Jan 12 19:21:22 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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/linterp_d.m Sat Jan 12 19:21:22 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;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/log2pie.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,5 @@ +function k=log2pie + persistent K; + if isempty(K) K=1+log(2*pi); end + k=K; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/logabsdet.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/logdet.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/logeps.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,4 @@ +function y=logeps(x), y=logth(eps,x); +% logeps - Natural log with lower threshold of eps +% +% logeps :: real -> real.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/logth.m Sat Jan 12 19:21:22 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<eps)=eps; y=log(x);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/map01.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,13 @@ +function x=map01(x,dim) +% map01 - Map values into the range [0,1]. +% +% map01 :: [Size], D:natural -> [Size]. + +if nargin<2, + sz=size(x); + dims=find(sz>1); + if isempty(dims), dim=1; else dim=dims(1); end +end + +x=divnorm(@(z)max(z,[],dim),addnorm(@(y)min(y,[],dim),x)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/mapinner.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,16 @@ +function C=mapinner(A,B) +% mapinner - inner product for multidim arrays +% +% mapinner :: [[N,M]], [[M,L]] -> [[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)]); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/maxnorm.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +% maxnorm - additive normalisation of vectors +% +% maxnorm :: [[N]] -> [[N]], real. +function [y,m]=maxnorm(x) + k=max(x); + if isfinite(k), y=x-k; + else y=zeros(size(x)); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/maxnormalise.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,6 @@ +function B=maxnormalise(A) +% maxnormalise - columnwise multiplicative normalisation maximum value +% +% maxnormalise :: [[N,M]] -> [[N,M]]. + +B=A./repmat(max(A,[],1),size(A,1),1);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/meandims.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/meanlastdims.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/measure.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,8 @@ +function b=measure(x,l) +% measure - returns vector in given direction but with a different length +% +% measure :: [[N]] -> [[N]] ~'unit vector of length 1'. +% measure :: [[N]], L:nonneg -> [[N]] ~'vector of length L'. + +if nargin<2, l=1; end +b=(l/norm(x))*x;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/minmax.m Sat Jan 12 19:21:22 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));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/mmax.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/mmin.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/mod1.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,3 @@ +% mod1 - Like mod but for ring 1..N instead of 0..N-1 +function z=mod1(x,y), z=1+mod(x-1,y); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/mouter.m Sat Jan 12 19:21:22 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))); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/mse.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +function E=mse(x1,x2) +% mse - Mean-square diff between arrays. +% +% mse :: [[Size]], [[Size]] -> nonneg. + +E=mean(flatten((x1-x2).^2)); + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/mulrows.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,3 @@ +function y=mulrows(x,k) + y=repmat(k,1,size(x,2)).*x; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/nan2x.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,6 @@ +function y=nan2x(x,y) +% nan2x - convert nans to something else +% +% nan2x :: real, [Size] -> [Size]. + +y(isnan(y))=x;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/nary.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,12 @@ +% 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/normalise.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,10 @@ +function x=normalise(x,dim) +% normalise - make x zero mean and unit variance in given dimension +% +% normalise :: [Size], D:natural -> [Size]. +% +% The return value has zero mean and unit variance in the +% Dth dimension. + +x=divnorm(@(z)std(z,1,dim),addnorm(@(y)mean(y,dim),x)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/orthogonalise.m Sat Jan 12 19:21:22 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)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/outer.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,4 @@ +function Y=outer(X), Y=X*X'; +% outer - outer product for 2D matrix +% +% outer :: [[N,M]] -> [[N,N]].
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/packvec.m Sat Jan 12 19:21:22 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{:}); + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/pdev.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,22 @@ +function N=pdev(p,A,D) +% pdev - p-deviation, generalisation of stddev for exponents other than 2 +% +% pdev :: nonneg, [[D1]], I:natural -> [[D2]]. +% pdev :: nonneg, [[N D]] -> [[ 1 D]]. +% +% where D2 = set(D1,I,1), ie the Ith dimension is +% collapsed. I defaults to 1. The p-deviation is defined as +% pdev(p,x,i) = mean(abs(x).^p,i).^(1/p) +% +% See also pnorm. + +if nargin<2, D=1; end; + +% Should these be sum or mean? +if p==2 + N=sqrt(mean(abs(A).^2,D)); +else + N=mean(abs(A).^p,D).^(1/p); +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/pnorm.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,22 @@ +function N=pnorm(p,A,D) +% pnorm - general vector norm wrt exponent p +% +% pnorm :: nonneg, [[D1]], I:natural -> [[D2]]. +% pnorm :: nonneg, [[N D]] -> [[ 1 D]]. +% +% where D2 = set(D1,I,1), ie the Ith dimension is +% collapsed. I defaults to 1. The p-norm is defined as +% pnorm(p,x,i) = sum(abs(x).^p,i).^(1/p) +% +% See also pdev for alternative which uses mean instead of sum + +if nargin<2, D=1; end; + +% Should these be sum or mean? +if p==2 + N=sqrt(sum(abs(A).^2,D)); +else + N=sum(abs(A).^p,D).^(1/p); +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/pnormalise.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,8 @@ +function A=pnormalise(p,A,i) +% pnormalise - normalise A relative to p-norm in given array dimension +% +% pnormalise :: nonneg, [D:[[1 E]]], [E] -> [D]. +if nargin<3, i=1; end +A=divnorm(@(x)pnorm(p,x,i),A); + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/prod1.m Sat Jan 12 19:21:22 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]);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/proddims.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/prodn.m Sat Jan 12 19:21:22 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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/product_iterator.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,26 @@ +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 + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/qform.m Sat Jan 12 19:21:22 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]].
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/qforma.m Sat Jan 12 19:21:22 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]].
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/qlog.m Sat Jan 12 19:21:22 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]]. +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/quadf.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,1 @@ +function y=quadf(x,s), y= -(x.^2)/(2*s); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/quantile.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,13 @@ +function Y=quantile(Q,X) +% quantile - compute quantiles of X +% +% quantile :: +% [[L]->0--1] ~'the L quantiles to compute', +% [[N,M]] ~'M columns of data' +% -> [[L,M]] ~'the L quantiles for each of M columns'. + +j=max(1,round(Q*size(X,1))); +S=sort(X); +Y=S(j,:); + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/range2set.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,1 @@ +function Y=range2set(I), Y=I(1):I(2);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/ratesched.m Sat Jan 12 19:21:22 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)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/rectify.m Sat Jan 12 19:21:22 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;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/rescale.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,8 @@ +function B=rescale(A,mx), B=A/max(A(isfinite(A))); +% rescale - rescale array values so maximum value is MX +% +% rescale :: [DIM->real], real -> [DIM->real] +% +% eg, B=rescale(A,10) +% returns matrix K*A such that maximum finite value is 10 +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/rot45.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/rpsi.m Sat Jan 12 19:21:22 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 +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/safe_exp.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,11 @@ +function [y,k]=safe_exp(x,l) +% safe_exp - high dynamic range exponential +% +% safe_exp :: X:[[N,M]] -> [[N,M]], [[1,M]]. +% +% returns y and k such that exp(X) = y * exp(k) and +% maximum value in y is 1. + +k=max(x,[],1); +y=exp(x-repmat(k,size(x,1),1)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/safe_exp2.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,11 @@ +function [y,k]=safe_exp2(x,l) +% safe_exp2 - high dynamic range exponential +% +% safe_exp2 :: X:[[N,M]], L:real -> [[N,M]], [[1,M]]. +% +% returns y and k such that exp(X) = y * exp(k) and +% maximum value in y is exp(L). +% Maximum sensible value of L is log(realmax)=709.78. + +k=max(x,[],1)-l; +y=exp(x-repmat(k,size(x,1),1));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/sbitmap.m Sat Jan 12 19:21:22 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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/shrink.m Sat Jan 12 19:21:22 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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/slog.m Sat Jan 12 19:21:22 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); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/smooth_with.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,24 @@ +function Z=smooth_with(h,i,Y) +% smooth_with - Smooth a signal by convolution without changing length +% +% smooth_with :: [[M]], 1..M, [[N]] -> [[N]]. + +m=length(h); +h=stoch(h(:))'; +csh=cumsum(h); +rsh=fliplr(cumsum(fliplr(h))); + +if isvector(Y), Z=sm(Y); +else Z=maprows(@sm,Y); +end + +function z=sm(y) + z=conv(h,y); % complete convolution + z=z(i:end-(m-i)); % trim the ends + + % rescale end bits + z(1:m-i)=z(1:m-i)./csh(i:end-1); + z(end-i+2:end)=z(end-i+2:end)./rsh(2:i); +end + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/spdiag.m Sat Jan 12 19:21:22 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));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/squeeze_prod.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/squeeze_sum.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/stddevs.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,6 @@ +function Y=stddevs(Q,X) +% stddevs - returns values a certain number of standard deviations from mean +% +% stddevs :: [[L]], [[N,M]] -> [[L,M]]. + +Y=Q(:)*std(X) + repmat(mean(X),length(Q),1);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/stoch.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,20 @@ +function [Y,Z]=stoch(X), +% stoch- make columns sum to one (like probabilities) +% +% stoch:: [[N,M]] -> [[N,M]], [[M]]. +% +% optionally returns the column sums as well. +% NOTE: this function assumes that X contains no NaNs and +% at most one Inf per column. Eg, the function will map +% +% [2;3;6;Inf;1;8] -> [0;0;0;1;0;0] +% +% See also: the function stochastic works on ROWS not columns. + + +Z=sum(X,1); +Y=X./repmat(Z,size(X,1),1); + +% any nans are assumed to be inf/inf so are replaced with 1 +Y(isnan(Y))=1; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/stochastic.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +function [Y,Z]=stochastic(X), +% stochastic - make rows sum to one (like probabilities) +% optionally returns the row sums as well +% +% stochastic :: [[N,M]] -> [[N,M]], [[N]]. + +Z=sum(X,2); +Y=X./repmat(max(Z,realmin),1,size(X,2)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/stochastic1.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,21 @@ +function [Y,Z]=stochastic1(X), +% stochastic1 - make columns sum to one (like probabilities) +% +% stochastic1 :: [[N,M]] -> [[N,M]], [[M]]. +% +% optionally returns the column sums as well. +% NOTE: this function assumes that X contains no NaNs and +% at most one Inf per column. Eg, the function will map +% +% [2;3;6;Inf;1;8] -> [0;0;0;1;0;0] +% +% See also: the function stochastic works on ROWS not columns. + + +Z=sum(X,1); + +Y=X./repmat(max(Z,realmin),size(X,1),1); + +% any nans are assumed to be inf/inf so are replaced with 1 +Y(isnan(Y))=1; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/sum1.m Sat Jan 12 19:21:22 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]);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/sumdims.m Sat Jan 12 19:21:22 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/sumn.m Sat Jan 12 19:21:22 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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/sym.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,4 @@ +function B=sym(A), B=A+A'; +% sym - Symmetric part of matrix (/2) +% +% sym :: [[N,N]] -> [[N,N]].
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/todbs.m Sat Jan 12 19:21:22 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); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/tozeros.m Sat Jan 12 19:21:22 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)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/ttimes.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,11 @@ +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)])); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/vecshift.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,17 @@ +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)); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/vlinterp.m Sat Jan 12 19:21:22 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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/vlinterp_d.m Sat Jan 12 19:21:22 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;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/zerobase.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,9 @@ +function x=zerobase(x,dim) +% zerobase - shift values so minimum is zero +% +% zerobase :: +% [D:[[E]]] ~'E dimensional array of size D', +% 1..D ~'dimension to normalise along' +% -> [D] ~'shifted data'. +x=addnorm(@(y)min(y,[],dim),x); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/zeromean.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,5 @@ +function x=zeromean(x,dim) +% zeromean(x,dim): make x zero mean +% if x is a matrix, operate on each ROW independently +x=addnorm(@(y)mean(y,dim),x); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/numerical/zeta.m Sat Jan 12 19:21:22 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 +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/popd.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,8 @@ +function D=popd +% popd: Pop directory off stack +% dirname=popd: pop and return directory name +% popd: pop directory and cd into it. +global DirStack +d=DirStack{1}; +DirStack=DirStack{2}; +if nargout<1, cd(d); else D=d; end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/prefs.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,41 @@ +function P=prefs(varargin) +% prefs - manage name-value map as a structure +% Values can be retreived with defaults using GETPARAM +% +% Arguments can be name value pairs, eg +% opts=prefs('foo',5,'bar',7) +% opts.foo=4 +% opts.bar=7 +% +% Or structures +% opts=prefs(defopts,'foo',5) +% +% Later arguments override earlier ones. A typical usage is +% with variable argument lists and GETPARAM, eg +% +% function somefun(a,b,varargin) +% ... +% opts=prefs('foo',5,varargin{:}); +% fooness=opts.foo; +% barness=getparam(opts,'bar',12); +% +% See also: GETPARAM, GETFIELD, ISFIELD + +P=struct; +n=1; +while n<=length(varargin) + arg=varargin{n}; + if isstruct(arg) + P=cpfields(fieldnames(arg),arg,P); + elseif iscell(arg) + pairs=arg; + for k=1:length(pairs) + P=setfield(P,pairs{k}{1},pairs{k}{2}); + end + else + P=setfield(P,varargin{n},varargin{n+1}); + n=n+1; + end + n=n+1; +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/pushd.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,7 @@ +function pushd(d) +% pushd: Push directory on directory stack +% pushd(dirname): push given directory +% pushd push current directory +global DirStack +if nargin<1, d=pwd; end +DirStack={d, DirStack};
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/repl.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,34 @@ +% repl - Read-eval-print loop. +% +% repl :: +% string ~'command line prompt', +% E:struct ~'structure fields to put in current environment', +% -> {...}. +% +% Runs an embedded command interpreter much like the Matlab top-level. +% The second argument E is available as a value with the name E in the environment. +% The environment can be examined using WHO or WHOS as usual, but there will be +% some other variables visible that are internal to repl - they all end with two +% underscores so they can be easily recognised. +% The repl can return any number of values, but the values must be placed +% by calling the ret(...) function within the repl. + +function varargout=repl(pr__,E) + fprintf('\n\nEntering REPL. Type "return","exit", or "quit" to finish.\n'); + fprintf('The function ret(..) sets return values.\n\n'); + cont__=true; + returns__={}; + while cont__, + str__=input([pr__,'>> '],'s'); + if strcmp(str__,'quit') break; end + if strcmp(str__,'exit') break; end + if strcmp(str__,'return') break; end + try, eval(str__) + catch ex__ + fprintf('\n%s\n',getReport(ex__)); + end + end + varargout=returns__; + + function ret(varargin), returns__=varargin; cont__=false; end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/sayf.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,6 @@ +function sayf(S,varargin) +% sayf - Use Mac OS X speech synthesiser like printf +% +% sayf(FormatString,...) + +system(sprintf(['say ',S],varargin{:}));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/sh.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,3 @@ +% sh - drop out to Unix shell +function sh +unix('bash');
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/tostring.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,22 @@ +function s=tostring(varargin) +% tostring - Tries to represent values as a sensible string +% +% tostring :: A -> string. +% tostring :: A, B -> string. +% etc. +% +% Multiple inputs are converted to comma-separated string. + s=cat_sep(',',map(@tostr,varargin)); + +function s=tostr(x) + if ischar(x), s=x; + elseif isnumeric(x), + sz=size(x); + if prod(sz)<16, s=mat2str(x); + else s=sprintf('%s[%s]',class(x),mat2str(sz)); end + elseif isa(x,'function_handle'), s=func2str(x); + if s(1)~='@', s=['@' s]; end + elseif isa(x,'func'), s=tostring(x); + else s=['obj:' class(x)]; + end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/general/tuple.m Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,4 @@ +function x=tuple(varargin), x=varargin; +% tuple - collect arguments into a cell array +% +% ie tuple(a,b,c,...) = {a,b,c,....}
--- a/sched/README Wed Jan 09 22:22:21 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,66 +0,0 @@ -Functions for timed execution - -Some of these functions rely on a high-resolution timer implemented -in Java - you will need to compile HRTimer.java and put it somewhere in -your classpath for these to work. - - -TYPES - - time ::= double ~ a type synonym standing for absolute time in seconds. - A1, A2 ... -> B1, ... ~ a proper function with no side effects. - A1, A2 ... => B1, ... ~ a 'function' which may be stateful and have side effects. - void ~ not really a type, but shorthand for an empty list of types. - - - timed_action(Ins:arglist(N),Outs:arglist(M)) ::= - double ~'scheduled time', Ins{:} - => double ~'timing error', Outs{:}. - - schedule ::= timed_action({},{schedule}). - - struct { F1; F2; .. } ~ a structure with *at-least* the fields specified by F1, F2 etc. - A field specification is like Name :: Type. - struct {} denotes the supertype of all structures. - - options { O1; O2; .. } ~ represents one or more arguments specifying optional named parameters. - - Options are a list of name/value pairs or structures containing name value pairs. I will - attempt to formalise this as follows: - - Option specification (the O1, O2 etc above) looks like - - option_spec --> Name :: Type - ; Name :: Type/Default where - - Option arguments can be passed as arbitrary sequences of name/value pairs or structs. If the - option specification is options Spec (where Spec is the list O1,O2 etc), then - - options_args(Spec) --> void - ; Name::string, Value::Type, options_args(Spec) - ; struct {}, options_args(Spec). - - -now_ns_hr - current time in nanoseconds using high resolution timer in JVM -nows - Current time in seconds -nows_hr - Current time using high resolution timer in JVM -sleeptil_hr - Sleep until absolute time - -rsched - regular state threading scheduler -schedrec_ - Schedule events where each event returns the next event. -schedrec_hr - execute action recursively at given time using high-resolution timer -sched_after - one shot, relative scheduler using Matlab timer -sched_at - Use one-shot timer to schedule action - -list_schedule - Make event action from list of times and parameters -null_schedule - Make empty schedule -st_schedule - Make state-threading event action from state transformer and state - -timed_action - return time-aware action function given ordinary action function - -timer_gc - Deletes all timers that have been released using releasetimer -timer_release - Mark timer for garbage collection -timer_wait - Wait for timer to stop running -isrunning - true if Matlab timer object is running -statuscb - Returns status message timer callback function handle -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sched/README.txt Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,66 @@ +Functions for timed execution + +Some of these functions rely on a high-resolution timer implemented +in Java - you will need to compile HRTimer.java and put it somewhere in +your classpath for these to work. + + +TYPES + + time ::= double ~ a type synonym standing for absolute time in seconds. + A1, A2 ... -> B1, ... ~ a proper function with no side effects. + A1, A2 ... => B1, ... ~ a 'function' which may be stateful and have side effects. + void ~ not really a type, but shorthand for an empty list of types. + + + timed_action(Ins:arglist(N),Outs:arglist(M)) ::= + double ~'scheduled time', Ins{:} + => double ~'timing error', Outs{:}. + + schedule ::= timed_action({},{schedule}). + + struct { F1; F2; .. } ~ a structure with *at-least* the fields specified by F1, F2 etc. + A field specification is like Name :: Type. + struct {} denotes the supertype of all structures. + + options { O1; O2; .. } ~ represents one or more arguments specifying optional named parameters. + + Options are a list of name/value pairs or structures containing name value pairs. I will + attempt to formalise this as follows: + + Option specification (the O1, O2 etc above) looks like + + option_spec --> Name :: Type + ; Name :: Type/Default where + + Option arguments can be passed as arbitrary sequences of name/value pairs or structs. If the + option specification is options Spec (where Spec is the list O1,O2 etc), then + + options_args(Spec) --> void + ; Name::string, Value::Type, options_args(Spec) + ; struct {}, options_args(Spec). + + +now_ns_hr - current time in nanoseconds using high resolution timer in JVM +nows - Current time in seconds +nows_hr - Current time using high resolution timer in JVM +sleeptil_hr - Sleep until absolute time + +rsched - regular state threading scheduler +schedrec_ - Schedule events where each event returns the next event. +schedrec_hr - execute action recursively at given time using high-resolution timer +sched_after - one shot, relative scheduler using Matlab timer +sched_at - Use one-shot timer to schedule action + +list_schedule - Make event action from list of times and parameters +null_schedule - Make empty schedule +st_schedule - Make state-threading event action from state transformer and state + +timed_action - return time-aware action function given ordinary action function + +timer_gc - Deletes all timers that have been released using releasetimer +timer_release - Mark timer for garbage collection +timer_wait - Wait for timer to stop running +isrunning - true if Matlab timer object is running +statuscb - Returns status message timer callback function handle +
--- a/sequences/README Wed Jan 09 22:22:21 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,174 +0,0 @@ -SEQUENCE CLASSES - -These class manage a variety of sequence construction methods. A sequence -object here is thought of as a lazy list or sequence of values, -where subsequent elements of the sequence are created only when -necessary. - -These are mostly inspired by the Haskell list processing library - - - -HASKELL LIST FUNCTIONS - -head :: [a] -> a -tail :: [a] -> [a] -null :: [a] -> bool -length :: [a] -> natural -last :: [a] -> a - -cons :: a -> [a] -> [a] -map :: (a->b) -> [a] -> [b] -append :: [a] -> [a] -> [a] -concat :: [[a]] -> [a] -concatMap :: (a->[b]) -> [a] -> [b] -interleave:: [a] -> [a] -> [a] -scanl :: (a->b->a) -> a -> [b] -> [a] -scanl1 :: (a->a->a) -> [a] -> [a] -mapAccumL :: (a->b->(a,c)) -> a -> [b] -> (a,[c]) -unfoldseq :: (b->maybe(a,b)) -> b -> [a] -zipWith :: (a->b->c) -> [a] -> [b] -> [c] -iterate :: (a->a) -> a -> [a] -repeat :: a -> [a] -replicate :: natural -> a -> [a] -cycle :: [a] -> [a] -take :: natural -> [a] -> [a] -splitAt :: natural -> [a] -> ([a],[a]) -filter :: (a->bool) -> [a] -takeWhile :: (a->bool) -> [a] -> [a] -span :: (a->bool) -> [a] -> ([a],[a]) -break :: (a->bool) -> [a] -> ([a],[a]) -partition :: (a->bool) -> [a] -> ([a],[a]) -init :: [a] -> [a] - -foldl :: (a->b->a) -> a -> [b] -> a -foldl1 :: (a->a->a) -> [a] -> a -drop :: natural -> [a] -> [a] -dropWhile :: (a->bool) -> [a] -> [a] - - -concat = foldl append [] -append x y = concat [x,y] -replicate n = take n . repeat -repeat x = cycle [x] -cycle = concat . repeat -concatMap = concat . map -once = take 1 - -naturals :: [natural] -fromTo :: natural -> natural -> [natural] -from :: natural -> [natural] -linseq :: real -> real -> natural -> [real]. -expseq :: real -> real -> natural -> [real]. -buflinear :: real -> real -> natural -> natural -> [arr(real)]. - - -naturals = iterate (+1) 0 -from = iterate (+1) -fromTo x y = take (y-x) (from x) -fromTo x y = takeWhile (<=y) (from x) - - -MY CLASSES - - The base class defining the interface is seq - seq abstract methods - next head tostring elsize - - The implementations are in the +seq package directory: - - bindcat cons map scanl subsample unfold_inf - buffer cycle mapaccum select take windowdata - cache iterate merge skin takewhile zipaccum - concat lcons repeat slices unfold zipwith - - - -METHODS THAT RETURN SEQUENCES - - Basic list combinators - - take drop dropwhile - map zip zipwith - mapaccum zipaccum scanl - and concat append - select merge cycle - span spanc split - cache once split_gather - diffwith bindcat - buffer unbuffer - - - Stateful processors (scanl, mapaccum, zipaccum) - - cumsum minmax integrate filter dynfilter - - Lifting ordinary functions and operators with maps and zips: - - le minus rdivide times mtimes - ctranspose power uminus paren gt - ldivide transpose mrdivide eq plus - lt mldivide ge subsref - vertcat horzcat cat - - binfun vecop binop unzip - isinf isnan isfinite - tanh cos mod sin - log10 log log2 exp - reshape sqrt abs - min max mean sum - powspec magspec phasespec fft ifft - - -METHODS THAT DO NOT RETURN SEQUENCES - - decons gather gathern length - limit display headfn size - nth1 numel imagesc - plot foldl last foreach - seq2cell meandata - - - -SEQUENCE CREATING FUNCTIONS - - cellseq - Convert cell array directly to sequence - framedata - return abstract data object for buffering signal frames - lazy_cons - Make sequence given head and tail - singleton - infinite sequence repeating one value - naturals - Sequence of natural numbers starting at 0 - integers - Sequence of integers - - rndscanl - Random scanl - rndseq - Sequence of values sampled from a random variable model - rndwindow - get random windows of a signal - rndzip - Random sequence by zipping argument sequences - rndzipaccum - Random sequence by zipping argument sequences with state - - scanseqcols - scandatacols - Scan over columns of array sequence - subseqdata - Sequence of subranges of a large array - - window - Window a sequnce of arrays in a given dimension - window_ns - Window a sequnce of arrays in a given dimension (no strict size) - windowparams - compute windowing parameters to get buffers of a certain size - - expdata - exponential data sequence - lindata - linear data sequence - - -SPECIFICATIONS - - (Some Haskellish stuff) - - 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} - buffer x n dim = cons (gather dim (take n x)) (buffer (drop n x) n dim ); - windowdata = concatScanl window1 (sort of) ?? - -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sequences/README.txt Sat Jan 12 19:21:22 2013 +0000 @@ -0,0 +1,174 @@ +SEQUENCE CLASSES + +These class manage a variety of sequence construction methods. A sequence +object here is thought of as a lazy list or sequence of values, +where subsequent elements of the sequence are created only when +necessary. + +These are mostly inspired by the Haskell list processing library + + + +HASKELL LIST FUNCTIONS + +head :: [a] -> a +tail :: [a] -> [a] +null :: [a] -> bool +length :: [a] -> natural +last :: [a] -> a + +cons :: a -> [a] -> [a] +map :: (a->b) -> [a] -> [b] +append :: [a] -> [a] -> [a] +concat :: [[a]] -> [a] +concatMap :: (a->[b]) -> [a] -> [b] +interleave:: [a] -> [a] -> [a] +scanl :: (a->b->a) -> a -> [b] -> [a] +scanl1 :: (a->a->a) -> [a] -> [a] +mapAccumL :: (a->b->(a,c)) -> a -> [b] -> (a,[c]) +unfoldseq :: (b->maybe(a,b)) -> b -> [a] +zipWith :: (a->b->c) -> [a] -> [b] -> [c] +iterate :: (a->a) -> a -> [a] +repeat :: a -> [a] +replicate :: natural -> a -> [a] +cycle :: [a] -> [a] +take :: natural -> [a] -> [a] +splitAt :: natural -> [a] -> ([a],[a]) +filter :: (a->bool) -> [a] +takeWhile :: (a->bool) -> [a] -> [a] +span :: (a->bool) -> [a] -> ([a],[a]) +break :: (a->bool) -> [a] -> ([a],[a]) +partition :: (a->bool) -> [a] -> ([a],[a]) +init :: [a] -> [a] + +foldl :: (a->b->a) -> a -> [b] -> a +foldl1 :: (a->a->a) -> [a] -> a +drop :: natural -> [a] -> [a] +dropWhile :: (a->bool) -> [a] -> [a] + + +concat = foldl append [] +append x y = concat [x,y] +replicate n = take n . repeat +repeat x = cycle [x] +cycle = concat . repeat +concatMap = concat . map +once = take 1 + +naturals :: [natural] +fromTo :: natural -> natural -> [natural] +from :: natural -> [natural] +linseq :: real -> real -> natural -> [real]. +expseq :: real -> real -> natural -> [real]. +buflinear :: real -> real -> natural -> natural -> [arr(real)]. + + +naturals = iterate (+1) 0 +from = iterate (+1) +fromTo x y = take (y-x) (from x) +fromTo x y = takeWhile (<=y) (from x) + + +MY CLASSES + + The base class defining the interface is seq + seq abstract methods + next head tostring elsize + + The implementations are in the +seq package directory: + + bindcat cons map scanl subsample unfold_inf + buffer cycle mapaccum select take windowdata + cache iterate merge skin takewhile zipaccum + concat lcons repeat slices unfold zipwith + + + +METHODS THAT RETURN SEQUENCES + + Basic list combinators + + take drop dropwhile + map zip zipwith + mapaccum zipaccum scanl + and concat append + select merge cycle + span spanc split + cache once split_gather + diffwith bindcat + buffer unbuffer + + + Stateful processors (scanl, mapaccum, zipaccum) + + cumsum minmax integrate filter dynfilter + + Lifting ordinary functions and operators with maps and zips: + + le minus rdivide times mtimes + ctranspose power uminus paren gt + ldivide transpose mrdivide eq plus + lt mldivide ge subsref + vertcat horzcat cat + + binfun vecop binop unzip + isinf isnan isfinite + tanh cos mod sin + log10 log log2 exp + reshape sqrt abs + min max mean sum + powspec magspec phasespec fft ifft + + +METHODS THAT DO NOT RETURN SEQUENCES + + decons gather gathern length + limit display headfn size + nth1 numel imagesc + plot foldl last foreach + seq2cell meandata + + + +SEQUENCE CREATING FUNCTIONS + + cellseq - Convert cell array directly to sequence + framedata - return abstract data object for buffering signal frames + lazy_cons - Make sequence given head and tail + singleton - infinite sequence repeating one value + naturals - Sequence of natural numbers starting at 0 + integers - Sequence of integers + + rndscanl - Random scanl + rndseq - Sequence of values sampled from a random variable model + rndwindow - get random windows of a signal + rndzip - Random sequence by zipping argument sequences + rndzipaccum - Random sequence by zipping argument sequences with state + + scanseqcols - scandatacols - Scan over columns of array sequence + subseqdata - Sequence of subranges of a large array + + window - Window a sequnce of arrays in a given dimension + window_ns - Window a sequnce of arrays in a given dimension (no strict size) + windowparams - compute windowing parameters to get buffers of a certain size + + expdata - exponential data sequence + lindata - linear data sequence + + +SPECIFICATIONS + + (Some Haskellish stuff) + + 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} + buffer x n dim = cons (gather dim (take n x)) (buffer (drop n x) n dim ); + windowdata = concatScanl window1 (sort of) ?? + +