## Copyright (C) 2013 - Juan Pablo Carbajal
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## 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, see <http://www.gnu.org/licenses/>.

## Author: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>

## -*- texinfo -*-
## @deftypefn {Function File} {@var{C} =} ndmult (@var{A},@var{B},@var{dim})
## Multidimensional scalar product
##
## Given multidimensional arrays @var{A} and @var{B} with entries
## A(i1,12,@dots{},in) and B(j1,j2,@dots{},jm) and the 1-by-2 dimesion array @var{dim}
## with entries [N,K]. Assume that
##
## @example
## shape(@var{A},N) == shape(@var{B},K)
## @end example
##
## Then the function calculates the product
##
## @example
## @group
##
## C (i1,@dots{},iN-1,iN+1,@dots{},in,j1,@dots{},jK-1,jK+1,@dots{},jm) =
##  = sum_over_s A(i1,@dots{},iN-1,s,iN+1,@dots{},in)*B(j1,@dots{},jK-1,s,jK+1,@dots{},jm)
##
## @end group
## @end example
##
## For example if @command{size(@var{A}) == [2,3,4]} and @command{size(@var{B}) == [5,3]}
## then the @command{@var{C} = ndmult(A,B,[2,2])} produces @command{size(@var{C}) == [2,4,5]}.
##
## This function is useful, for example,  when calculating grammian matrices of a set fo signals
## produced from different experiments.
## @example
##   nT      = 100;
##   t       = 2*pi*linspace (0,1,nT)';
##   signals = zeros(nT,3,2); % 2 experiments measuring 3 signals at nT timestamps
##
##   signals(:,:,1) = [sin(2*t) cos(2*t) sin(4*t).^2];
##   signals(:,:,2) = [sin(2*t+pi/4) cos(2*t+pi/4) sin(4*t+pi/6).^2];
##
##   sT(:,:,1) = signals(:,:,1)';
##   sT(:,:,2) = signals(:,:,2)';
##   G = ndmult (signals,sT,[1 2]);
##
## @end example
## In the exmaple G contains the scalar product of all the singals against each other.
## This can be verified in the following way:
## @example
##   sA = 1 eA = 1; % First signal in first experiment;
##   sB = 1 eA = 2; % First signal in second experiment;
##   [G(s1,e1,s2,e2)  signals(:,s1,e1)'*signals(:,s2,e2)]
## @end example
## You may want to reoeder the scalar products into a 2-by-2 arrangement (representing pairs of experiments)
## of gramian matrices. The following command @command{G = permute(G,[1 3 2 4])} does it.
##
## @end deftypefn

function M = ndmult (A,B,dim)

  dA = dim(1);
  dB = dim(2);

  sA   = size (A);
  nA   = length (sA);
  perA = [1:(dA-1) (dA+1):(nA-1) nA dA](1:nA);
  Ap = permute (A, perA);
  Ap = reshape (Ap, prod (sA(perA(1:end-1))), sA(perA(end)));

  sB   = size (B);
  nB   = length (sB);
  perB = [dB 1:(dB-1) (dB+1):(nB-1) nB](1:nB);
  Bp = permute (B, perB);
  Bp = reshape (Bp, sB(perB(1)), prod (sB(perB(2:end))));

  M  = Ap * Bp;
  s = [sA(perA(1:end-1)) sB(perB(2:end))];
  M  = squeeze (reshape (M, s));

  #TODO: generalize to mmult (basis,p,dim)
  # n = length(size(basis))
  # per = [[1:(dim-1) (dim+1):(n-1) n dim](1:n)
  # permute (basis, per), reshape, multiply, reshape,
  #
#   [nT dimQ nPhi order] = size(basis);
#    q = reshape (reshape (
#                        permute (basis, [1 2 4 3]), ...
#                                 nT*dimQ*order,nPhi) * p, ...
#                        nT, dimQ,order);

endfunction
