/* -*- mode: C++ -*- */
//////////////////////////////////////////////////////////////////////////////
// FILE   : tools.tol
// PURPOSE: Defines Class BysVecLinReg::@Multivariate
//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
//Converts M in a column vector of m*n rows with stacked columns of M having
//m rows and n columns
VMatrix vmatrix2vec(VMatrix M)
//////////////////////////////////////////////////////////////////////////////
{
  Matrix Mt = VMat2Mat(M,True);
  Mat2VMat(Tra(GetNumeric([[ Mt ]])))
};

//////////////////////////////////////////////////////////////////////////////
//Converts a column vector into a tye 
VMatrix vec2vmatrix(VMatrix vec, VMatrix M)
//////////////////////////////////////////////////////////////////////////////
{
  Matrix Mt = VMat2Mat(M,True);
  Mat2VMat(DeepCopy([[Mt]],VMat2Mat(vec))[1],True)
};

//////////////////////////////////////////////////////////////////////////////
//Returns the position of cell (i,j) of in stacked columns form of a matrix 
//with r rows and c columns
Real vmatrix2vec.getPosition(Real i, Real j, Real r, Real c)
//////////////////////////////////////////////////////////////////////////////
{
  MatDat(VMat2Triplet(vmatrix2vec(Triplet([[ [[i,j,1]] ]],r,c))),1,1)
};

//////////////////////////////////////////////////////////////////////////////
//Draws an standard multinormal x with full rank constraining linear 
//equalities
// 
//  e ~ N(0,I) 
//  C*e = c 
//
// where
//  e <: R^n
//  C <: R^rxn, r<n, rank(C)=r
//  c <: R^r
//
//ALGORITHM:
// 
// 1. Get pivoting order p of columns of diagonal of C*C'
// 2. Select p.1 as gaussian components of x as (n-r) last pivoted columns
// 3. Select p.2 as dependent components of x as first r pivoted columns
// 4. Draw e.1 as standard independent normal 
// 5. Rewrite constraining linear equalities conditioned to e.1
//     
//      C = (C.1 C.2)
//      e' = (e.1' e.2')
//      C.1 * e.1 + C.2 * e.2 = c
//      c.2 = c - C.1*e.1
//      C.2 * e.2 = c.2
// 
// 6. Solve linear system to get e.2
// 7. Checks if solution has the expected rounding error
// 8. Merge e.1 and e.2 following the selected pivoting order to get the
//    wanted draw of e
// 
VMatrix eqConstrainedStdGaussian(VMatrix C, VMatrix c)
//////////////////////////////////////////////////////////////////////////////
{
  Real r = VRows(C);
  Real n = VColumns(C);
  //Calculates diagonal of C*C' in a faster way
  VMatrix C2 = Constant(1,r,1) * (C $* C);
  //Builds pivoting order of columns of diagonal of C*C'
  Matrix p = ColumnPivots(VMat2Mat(C2));
  //Select p.1 as gaussian components of x as (n-r) last pivoted columns
  Set p.1 = MatSet(Sub(p,1,r+1,1,n-r))[1];
  //Select p.2 as dependent components of x as first r pivoted columns
  Set p.2 = MatSet(Sub(p,1,1,  1,r  ))[1];
  //Draw e.1 as standard independent normal 
  VMatrix e.1 = Gaussian(n-r,1,0,1);
  //Rewrite constraining linear equalities conditioned to e.1
  VMatrix C.1 = SubCol(C,p.1);
  VMatrix C.2 = SubCol(C,p.2);
  VMatrix c.2 = c-C.1 * e.1;
  //Solve linear system to get e.2
  VMatrix e.2 = MinimumResidualsSolve(C.2,c.2);
  //Checking method success
  VMatrix err = C.2*e.2-c.2;
  Real errMax = VMatMax(Abs(err));
  If((n*r^2)+errMax^2 != n*r^2,
    WriteLn("[BysVecLinReg::eqConstrainedStdGaussian] errMax="<<errMax,"W"));
  //Merging gaussian and dependent components following the selected 
  //pivoting order
  Triplet(SetCol(p.1<<p.2) | Constant(n,1,1) | VMat2Mat(e.1 << e.2), n, 1)
};