/* -*- mode: C++ -*- */
//////////////////////////////////////////////////////////////////////////////
// FILE   : BaseC2.LinComb.tol
// PURPOSE: Defines class @BaseC2.LinComb
//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
Text _.autodoc.member.BaseC2.LinComb = "Defines a new base B2 of "
"order K2 as a linear combination B2=B1*V of functions of a previous base "
"B1 of order K1>=K2. The matrix of change V must have rank K2 to ensure "
"linear independence.";
Class @BaseC2.LinComb : @BaseC2.Cache
//////////////////////////////////////////////////////////////////////////////
{
  @BaseC2 _.base; 
  VMatrix _.V;
  VMatrix _.Vi;

  ////////////////////////////////////////////////////////////////////////////
  Static Text _.autodoc.member.New = "Creates an instance B2 of "
  "@BaseC2.LinComb as B2 = B1*V";
  Static @BaseC2.LinComb New(@BaseC2 B1, VMatrix V, VMatrix Vi)
  ////////////////////////////////////////////////////////////////////////////
  {
    Real n = B1::_.n;
    
    @BaseC2.LinComb aux = [[
      Real _.n = B1::_.n;
      Real _.K = VColumns(V);
      @BaseC2 _.base = base; 
      VMatrix _.V = V; 
      VMatrix _.Vi = Vi
  ]]};

  ////////////////////////////////////////////////////////////////////////////
  Static Text _.autodoc.member.Compact.Ratio = "Creates an instance of "
  "@BaseC2.LinComb that is close to be orthonormal and has at least "
  "the specified ratio of information of original base";
  Static @BaseC2.LinComb Compact.Ratio(
    @BaseC2 base, 
    Real informationRatio)
  ////////////////////////////////////////////////////////////////////////////
  {
    Real K1 = base::_.K;

    VMatrix u = Rand(K1*10, K1, -1, 1);

    VMatrix Pu = Group("ConcatColumns", For(1, K1, VMatrix(Real k)
    {
      base::function.v(k,u)
    }));

    Set Pu.svd = SVD(VMat2Mat(Pu));

    Matrix Pu.eigenvalue = Tra(SubDiag(Pu.svd[2],0));
    Matrix Pu.eigenvalueCum = 
    {
       Matrix aux = DifEq(1/(1-B),Pu.eigenvalue);
       aux*(1/MatSum(Pu.eigenvalue))
    };
    Matrix ratio = Rand(K1,1,informationRatio,informationRatio);
    Real K = MatSum(LE(Pu.eigenvalueCum,ratio));
    Set range.K = Range(1,K,1);
    VMatrix V = Mat2VMat(SubCol(Pu.svd[3],range.K));
    VMatrix D = Eye(K,K,0,RPow(Mat2VMat(SubRow(Pu.eigenvalue,range.K)),+1));
    VMatrix Di = Eye(K,K,0,RPow(Mat2VMat(SubRow(Pu.eigenvalue,range.K)),-1));
    VMatrix W = V*Di;
    VMatrix Wi = D*Tra(V);
    New(base, W, Wi)
  };

  ////////////////////////////////////////////////////////////////////////////
  Static Text _.autodoc.member.Compact.K = "Creates an instance of "
  "@BaseC2.LinComb that is close to be orthonormal and has exactly K "
  "elements";
  Static @BaseC2.LinComb Compact.K(
    @BaseC2 base, 
    Real K)
  ////////////////////////////////////////////////////////////////////////////
  {
    Set range.K = Range(1,K,1);
    Real K1 = base::_.K;

    VMatrix u = Rand(K1*10, K1, -1, 1);

    VMatrix Pu = Group("ConcatColumns", For(1, K1, VMatrix(Real k)
    {
      base::function.v(k,u)
    }));

    Set Pu.svd = SVD(VMat2Mat(Pu));

    Matrix Pu.eigenvalue = Tra(SubDiag(Pu.svd[2],0));
    VMatrix V = Mat2VMat(SubCol(Pu.svd[3],range.K));
    VMatrix D = Eye(K,K,0,RPow(Mat2VMat(SubRow(Pu.eigenvalue,range.K)),+1));
    VMatrix Di = Eye(K,K,0,RPow(Mat2VMat(SubRow(Pu.eigenvalue,range.K)),-1));
    VMatrix W = V*Di;
    VMatrix Wi = D*Tra(V);
    New(base, W, Wi)
  };

  ////////////////////////////////////////////////////////////////////////////
  Real _build.cache(VMatrix X, Real order) 
  ////////////////////////////////////////////////////////////////////////////
  { 
    Real _.base::upgrade.cache(X,order);
    VMatrix If(order>=0, { _.function := 
      _.base::function.m(X)*_.V
    });
    Set If(order>=1, { _.partial.derivative.1 := 
      For(1, _.n, VMatrix(Real i)
      { 
        _.base::partial.derivative.1.m(X,i)*_.V 
      })
    });
    Set If(order>=2, { _.partial.derivative.2 := 
      For(1, _.n, Set(Real i)
      { 
        For(1, i, VMatrix(Real j)
        { 
          _.base::partial.derivative.2.m(X,i,j)*_.V 
        })
      })
    });
    True
  }

};