1function varargout = gismo_map(pts, in, der)
2  if ( iscell(pts) ) % create cartesian product points from cell infomation
3      ndim = length(pts);
4      npts = prod(cellfun(@length,pts));
5      pts = cartesian_product_from_cell(pts);
6  else
7      [ndim, npts] = size(pts);
8  end
9
10  if (der == 0 || nargout > 1)
11    F = in.eval(pts); % dimension rdim x npts
12    varargout{1} = F;
13  end
14  if (der == 1 || nargout > 2)
15    jac = reshape(in.deriv(pts),[],ndim,npts); % dimension rdim x ndim x npts
16    if nargout == 1
17        varargout{1} = jac;
18    else % nargout == 2 or 3
19        varargout{2} = jac;
20    end
21  end
22  if (der == 2)
23      rdim = length(in.eval(zeros(ndim,1)));
24      hess = zeros(rdim,ndim,ndim,npts); % dimension rdim x ndim x ndim x npts
25%       for dir = 1:rdim
26%         hess(dir,:,:,:) = reshape(in.hess(pts,dir),ndim,ndim,npts);
27%       end
28    if nargout == 1
29        varargout{1} = hess;
30    elseif nargout == 3
31        varargout{3} = hess;
32    end
33  end
34end
35
36function pts_aux = cartesian_product_from_cell(pts)
37  % create cartesian product points from cell infomation
38  s = cellfun(@length,pts);
39  s_cell = cell(length(s),1);
40  for ii = 1:length(s)
41      s_cell{ii} = 1:s(ii);
42  end
43  x = cell(1,numel(s_cell));
44  [x{:}] = ndgrid(s_cell{:});
45  pts_aux = [];
46  for ii=1:length(s)
47      pts_aux = [pts_aux; reshape(pts{ii}(x{ii}),1,[])];
48  end
49end