//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [C] = %hm_m_s(A,B)

// A is considered as a set of matrices: 
// A = [A1, A2, ...Ak], k=1:K, with size of A = NxPxK
// B is considered as a set of (column) vectors : 
// B = [B1, B2, ...Bk], k=1:K with size of B = PxK
// Then C is defined by: 
// C = [A1*B1, A2*B2, ... Ak*Bk]; size(C) = (N,K) 
// special cases : 
// - B scalar => C = A * B(1,1)
// - B is a column vector => 
// B considered as [B,B,..., B] (K times)
//
// WARNING: %hm_m_s overloads Scilab's standard definition


// Declarations:

// Code:

if (size(A.dims,'*') ~= 3)
  CL__error("Only hypermatrices with 3 dimensions are handled");
end

[N,P,K] = size(A);
[P2,K2] = size(B);


// --- special case : B = [] 

if (P2 == 0 & K2 == 0)  
  // A = hypermatrix , B = []
  C = [];
  return;  // <== RETURN
end 

// --- special case : B = scalar 

if (P2 == 1 & K2 == 1)  
  // A = hypermatrix , B = scalar (1x1)
  C = zeros(A);
  C(:) = A(:) * B(1,1);
  return;  // <== RETURN
end 


// --- A = hypermatrix , B = matrix (or vector)

if (K2 == 1)
  B = B * ones(1, K); // B now has K columns
  K2 = K;
end

if (P2 <> P | K2 <> K)
  CL__error("Hypermatrix and matrix have incompatible sizes");
end

C = zeros(N,K);

// Note : 
// we note Aik = A(i,:,k)(:), and Bk = B(:,k)(:)
// (Aik and Bjk = column vectors of P elements)
// => A(i,:,:) contains IN THIS ORDER: Ai1; Ai2;... Aik
// which can be reformatted COLUMN BY COLUMN: 
// matrix(A(i,:,:),P,K) = [Ai1, Ai2... ]  -> size = PxK
// => B(:,:) contains IN THIS ORDER: B1; B2;... Bk
// which can be reformatted COLUMN BY COLUMN:
// matrix(B(:,:),P,K) = B = [B1, B2,... ] -> size = PxK
// => easy to compute sum(Ai1 .* B1) = C(i,1), etc... 

for (i = 1:N)
  C(i,:) = sum( matrix(A(i,:,:),P,K) .* B, 'r' );
end


endfunction