samer@4: function varargout=zipmap(spec,f,varargin)
samer@4: % zipmap - Map and zip function over multidimensional arrays
samer@4: %
samer@4: % ZIPMAP allows a certain simple functions to be applied to multidimensional
samer@4: % array arguments, manipulating the dimensions to combine elementwise and
samer@4: % outer-product-like behaviour. (The name comes from functional programming
samer@4: % terminology: functions/operators like .* and + which apply elementwise to
samer@4: % vector arguments are said to zip their arguments, while outer products are
samer@4: % like mapping a function over a list, or a list of functions over another list.)
samer@4: %
samer@4: % Let f be a vectorised function which knows how to zip array arguments. Eg,
samer@4: % f(x,y)=x+y. If x and y are vectors, they must be the same size, and
samer@4: % the function is applied to corresponding pairs of elements, resulting in
samer@4: % a vector of the same size. With ZIPMAP, you can build functions, that, eg,
samer@4: % add EVERY pair of elements (not just corresponding pairs) resulting in a
samer@12: % 2D array: zipmap(0,@(x,y)x+y,[1:4]',[1:6]') will do this
samer@4: %
samer@4: % The domain of a multidimensional array can be described a list containing the
samer@4: % size of each dimension, as returned by the SIZE function (or SIZE1 which strips
samer@4: % trailing 1s). The arguments to f must have a certain domain structure:
samer@4: %
samer@4: % 	dom arg1 = [p1 zipdom]
samer@4: % 	dom arg2 = [p2 zipdom]
samer@4: %	         :
samer@4: %
samer@4: % where p1, p2, etc are the domains of the elementary function of which f
samer@4: % is a vectorisation, eg the elementary domain of .* can be scalars and
samer@4: % hence p1=p2=[], but the two elementary domains of CROSS are both [3], because it 
samer@4: % it operates on pairs of 3-D vectors. The zipdom part must be the same for
samer@4: % all arguments, otherwise, we cannot extract matching argument tuples.
samer@4: % If this is the case, the f must return a result:
samer@4: %
samer@4: % 	dom f(...) : [rdom zipdom]
samer@4: %
samer@4: % where rdom is the domain of the elementary function, eg [3] for the vector
samer@4: % valued CROSS product, but [] for the scalar valued DOT product. Note that
samer@4: % these domains can be vectors of any length from zero up, eg the determinant
samer@4: % of a square matrix has a domain of [n n].
samer@4: % 
samer@4: % ZIPDIM allows the domains of the arguments and result to be generalised as follows:
samer@4: %
samer@4: % 	dom arg1 : [p1 m1 zipdom]
samer@4: % 	dom arg2 : [p2 m2 zipdom]
samer@4: %	         :
samer@4: %
samer@4: %	dom zipdom(...) : [rdom m1 m2 ... zipdom ]
samer@4: %
samer@4: % All the mapped dimensions m1, m2 etc go into a sort of outer product where 
samer@4: % every combination of values is evaluated. In addition, the order of the
samer@4: % m1, m2 etc can be permuted in the result at essentially no extra computational cost.
samer@4: %
samer@4: % Usage: r=zipmap(spec,f,<varargs>)
samer@4: % 
samer@4: % spec defines how the domain of each of the arguments is to be interpreted, ie
samer@4: % how is is split into elementary, mappable, and zippable domains.
samer@4: %
samer@4: % If spec is a number z, then the last z dimensions of each argument are marked for
samer@4: % zipping, while the rest are for mapping. The elementary domains are assumed to
samer@4: % be scalar and hence p1, p2 etc =[]. So, zipmap(0,...) does the complete outer
samer@4: % product type of evaluation.
samer@4: %
samer@4: % If spec is structure, the field 'zipdims' determines how many dimensions are zipped.
samer@4: % The field 'pdims' is a the vector [length(p1) length(p2) ..], ie the dimensionality
samer@4: % of each elementary domain: 0 for scalars, 1 for vectors, 2 for matrices etc.
samer@4: % The optional field 'mapord' determines the reordering of the mapped domains, eg
samer@4: % if mapord=[2 1 3], then the result domain is [m2 m1 m3] instead of [m1 m2 m3].
samer@4: %
samer@4: % If spec is a cell array, it is interpreted as {zipdims pdims mapord}.
samer@4: %
samer@4: % Extra options
samer@4: %
samer@4: % If spec is a structure, it can contain further options: if spec.shift=1,
samer@4: % then the domains of all arguments are pre-shifted to remove singleton
samer@4: % dimensions before anything else is done. Amongst other effects, this means
samer@4: % that row vectors can be used in place of column vectors.
samer@4: 
samer@12: 	% unpack args from spec
samer@12: 	if iscell(spec)
samer@12: 		[zipdims,pdims,mapord]=getargs(spec,{0 [] []});
samer@12: 	elseif isstruct(spec)
samer@12: 		zipdims=getparam(spec,'zipdims',0);
samer@12: 		pdims=getparam(spec,'pdims',[]);
samer@12: 		mapord=getparam(spec,'mapord',[]);
samer@12: 		if getparam(spec,'shift',0),
samer@12: 			varargin=cellmap(@shiftdim,varargin);
samer@12: 		end
samer@12: 	elseif isscalar(spec)
samer@12: 		zipdims=spec; pdims=[]; mapord=[];
samer@12: 	end
samer@4: 
samer@12: 
samer@12: 	argdoms=cellmap(@size1,varargin);
samer@12: 	if isempty(pdims), pdims=zeros(size(varargin)); end
samer@12: 	for i=1:length(varargin)
samer@12: 		[pdoms{i},mapdoms{i},zipdoms{i}]=split(argdoms{i},pdims(i),zipdims);
samer@4: 	end
samer@12: 
samer@12: 
samer@12: 	% check all zipdoms match and use first
samer@12: 	zipdom=zipdoms{1};
samer@12: 	for i=2:length(zipdoms)
samer@12: 		if ~all(zipdoms{i}==zipdom), error('Zip domains do not match'); end
samer@12: 	end
samer@12: 
samer@12: 
samer@12: 	if isempty(mapord), mapord=1:length(varargin); end
samer@12: 	outerdom=cat(2,mapdoms{mapord});
samer@12: 	outer1=cellmap(@to1s,mapdoms);
samer@12: 	for i=1:length(varargin)
samer@12: 		% need to reshape each arg and then repmat up to correct size
samer@12: 		newsizes=outer1; newsizes{i}=mapdoms{i};
samer@12: 		newsize=[pdoms{i} cat(2,newsizes{mapord}) zipdom];
samer@12: 		varargin{i}=reshape(varargin{i},tosize([pdoms{i} cat(2,newsizes{mapord}) zipdom]));
samer@12: 		newsizes=mapdoms; newsizes{i}=outer1{i};
samer@12: 		newsize=[to1s(pdoms{i}) cat(2,newsizes{mapord}) to1s(zipdom)];
samer@12: 		varargin{i}=repmat(varargin{i},tosize(newsize));
samer@12: 	%	size(varargin{i})
samer@12: 	end
samer@12: 
samer@12: 	[varargout{1:max(1,nargout)}]=feval(f,varargin{:});
samer@4: end
samer@4: 
samer@4: function [a,b,c]=split(x,n,m)
samer@4: 	a=x(1:n);
samer@4: 	b=x(n+1:end-m);
samer@4: 	c=x(end-m+1:end);
samer@12: end
samer@4: 
samer@4: % getargs - get values of optional parameters with defaults
samer@4: %
samer@4: % getargs :: {I:[N]->(A(I) | nil)}, {I:[N]->A(I)} -> A(1), ..., A(N).
samer@4: function varargout=getargs(args,defs)
samer@12: 	nout=max(length(args),length(defs));
samer@12: 	varargout=defs;
samer@12: 	for i=1:length(args)
samer@12: 		if ~isempty(args{i}), varargout(i)=args(i); end
samer@12: 	end
samer@4: end
samer@4: