Ticket #79: funciones.tol

//////////////////////////////////////////////////////////////////////////////
Real CheckLinCom( Set OutputInfo, Set InputsInfo, Set AllHierarchyInfo,
  Text obs.node, Real vl_inf, Real vl_pri, Text logRoute)
//////////////////////////////////////////////////////////////////////////////
{
  // We extract all parameters involved in any hierarchy for node obs.node
  Set sHieNodeInfo = BinGroup("<<",EvalSet(AllHierarchyInfo,Set (Set s)
  {
    Set prev1 = Extract(s,2)|Extract(s,1);

    Set prev2 = Select(prev1,Real (Set x)
    {
      Real If(x[1] == obs.node, True, False)
    }) - [[ Empty ]];

    Set BinGroup("<<",Extract(prev2,2))
  }));


  Real c = 0;

  Set inputs = Select(OutputInfo::InputDB,Real (Anything x)
  {
    Real c := c+1;

    // Real flag is a label that tell us if an input has not any prior, or
    // has not any unknown value, or if it is involved in any hierarchy
tree

    Real flag =
      If(
        Name(InputsInfo[c]) <: sHieNodeInfo,
        Real False,
        If(
          Not(And( BinEQ( ?, (InputsInfo[c])->prior_mu ),
                   BinEQ( ?, (InputsInfo[c])->prior_sigma ))),
          Real False,
          If(
            Grammar(x) == "Serie",
            Real Not(MatSum(IsUnknown(SerMat(x)))),// HasUnknown is slower
            If(
              Grammar(x) == "Matrix",
              Real Not(MatSum(IsUnknown(x))),// HasUnknown is slower
              Real True
            )
          )
        )
      );

    flag
  });

  Real If( Not(Card(inputs)),
    Real True,
    Real LinComWarning(OutputInfo,inputs,logRoute,id_model,obs.node)
  )
};
//////////////////////////////////////////////////////////////////////////////
PutDescription(
"It checks a set of input series (or matrix) to find linear combinations
between the observational node parameters. It needs an OutputInfo BSR set
of
input series, the observational node, the real variable vl_pri (True or
False)
and a Path to store a log file with the report of the parameters that has
linear combinations between them.
If exist linear combination, it writes in log file and kill TOL",
CheckLinCom);
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Real LinComWarning(Set OutputInfo, Set inputs, Text logRoute, Text
id_model,
  Text node)
//////////////////////////////////////////////////////////////////////////////
{
  Real matrix = If( id_dating == "Matrix",True,False);

  Set lcw = If( Not(matrix), LinComWarningSerie( OutputInfo,inputs ),
    LinComWarningMatrix( inputs ));

  Real If( EQ(Card(lcw),0), True,
  {
    Text WriteLn("[LinComWarning] THERE ARE "+IntText(Card(lcw))+" LINEAR
"+
      "COMBINATIONS OF INPUT SERIES CHECK LOG DIRECTORY!!", Text "E");

    Text txt2File =
    "update "+TBsrParameter+" set
co_active="+SqlFormatText("N",GesAct)+""+NL+
    "
"+NL+
    "where co_node  = "+SqlFormatText(node,GesAct)+"
"+NL+
    "  and co_model = "+SqlFormatText(id_model,GesAct)+"
"+NL+
    "
"+NL+
    "and co_parameter in
"+NL+
    "(
"+NL+
    "
"+NL+
    TxtListItemQuote(lcw,","+NL)+"
"+NL+
    "
"+NL+
    ")";

    Text WriteFile(logRoute,txt2File);

    Real False
  })
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("It takes a set of series and check them to find linear
combinations. If thre are linear combinations, it gives an error message
and
writes a query in a log file in the route logRoute.",
LinComWarning);
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Set LinComWarningSerie(Set OutputInfo, Set inputs)
//////////////////////////////////////////////////////////////////////////////
{
  Set  series   = inputs;
  Date TruncIni = OutputInfo::IniEstim;
  Date TruncEnd = OutputInfo::EndEstim;

  Polyn dif = SetProd(EvalSet(OutputInfo::Arima,Polyn (Set p){ p[4] }));

  Set seriesTrunc = EvalSet(series,Serie (Serie ser)
  {
    Serie ser2 = SubSer(ser,TruncIni,TruncEnd);
    Serie dif:(ser2/MaxS(Abs(ser2)))
  });

  Matrix A = MatSetSeries(seriesTrunc);

  // SVD algorithm builds a set of matrix, the second matrix has the
  // eigenvalues.
  /* Set svd1 = SVDDecomposition(A*Tra(A)); */ // Luis -> SVD no es
estable
  Set svd1 = SVD(A*Tra(A));

  // We transform de second matrix into a rounded matrix because we need
  // to search eigenvalues = 0. The 10-10 edge is 1000 times the cholesky
  // rounding error so i understand that something below this limit is
negligible
  Set svd2 = RoundSVD(svd1,False,0.5,10^(-10));

  // we take de main diagonal of the eigenvalues matrix and we check if
exist
  // eigenvalues = 0
  Matrix diag = SubDiag( svd2[2] ,0 );
  Real   prod = MatSet(diag)[1][Card(MatSet(diag)[1])];

  Set lincom = If(BinEQ(prod,?),
  {
    Set Lin.Com.Procedure(series,svd2)
  },
  {
    Set If(prod,
    {
      Set Empty
    },
    {
      Set Lin.Com.Procedure(series,svd2)
    })
  })
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("It takes a set of series and gives you the number of
linear
combinations that happens",
LinComWarningSerie);
//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
Set LinComWarningMatrix(Set inputs)
//////////////////////////////////////////////////////////////////////////////
{
  Set vectors = inputs;

  Matrix A = Tra(BinGroup("|",vectors));

  // SVD algorithm builds a set of matrix, the second matrix has the
  // eigenvalues.
  /* Set svd1 = SVDDecomposition(A*Tra(A)); */ // Luis -> SVD no es
estable
  Set svd1 = SVD(A*Tra(A));

  // we transform de second matrix into a rounded matrix because we need
  // to search eigenvalues = 0
  Set svd2 = RoundSVD(svd1,False,0.5,10^(-10));

  // we take de main diagonal of the eigenvalues matrix and we check if
exist
  // eigenvalues = 0
  Matrix diag = SubDiag( svd2[2] ,0 );
  Real   prod = MatSet(diag)[1][Card(MatSet(diag)[1])];

  Set lincom = If(BinEQ(prod,?),
  {
    Set Lin.Com.Procedure(vectors,svd2)
  },
  {
    Set If(prod,
    {
      Set Empty
    },
    {
      Set Lin.Com.Procedure(vectors,svd2)
    })
  })
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("It takes a set of matrix and gives you the number of
linear
combinations that happens",
LinComWarningMatrix);
//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
Set Lin.Com.Procedure(Set elements, Set svd)
//////////////////////////////////////////////////////////////////////////////
{
  // If exist eigenvalues = 0, we rounded to zero the smallest values of
the
  // last column of the third SVD matrix, and then, the not zero values
show
  // the linear combinations.

  Matrix  lastCol = SubCol(svd[3],[[ Columns(svd[3]) ]]);
  Matrix  LogicLC = Not(Not(lastCol));
  Real    maxDiag = MaxMatrix(SubDiag(svd[2],0));

  // We take the linear combinatioNs and store them into a set
  Set indices = BinGroup("<<",MatSet(LogicLC));

  Real c = 0;

  Set select  = Select(elements,Real (Anything ser)
  {
    Real c := c+1;
    Real If(maxDiag,indices[c],True)
  });

  Set EvalSet(select,Text (Anything x)
  {
    Text Name(x)
  })
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("It needs a set of elements (series or matrix) and the
result
of a svd algorithm, and gives you tha distinct variables of elements Set
that
does linear combinations between them",
Lin.Com.Procedure);
//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
Set RoundSVD(Set svd,Real trasposed, Real minSparse,Real tolerance)
//////////////////////////////////////////////////////////////////////////////
{
  Set EvalSet(svd,Matrix (Matrix matrix)
  {
    VMatrix vmatrix = Mat2VMat(
      matrix,     // matrix to transform
      trasposed,  // does it trasposed?
      minSparse,  // minSparse
      tolerance); // tolerance

    Matrix VMat2Mat(vmatrix)
  })
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("It rounds an exit of SVD function to avoid negligible
values.
It needs:
  matrix     --> matrix to transform
  trasposed  --> does it trasposed? (True or False)
  minSparse  --> minSparse use 0.5 by default
  tolerance  --> round use 10^(-10) recommended
",
RoundSVD);
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Set SVDDecomposition(Matrix m)
//////////////////////////////////////////////////////////////////////////////
{
  VMatrix vm = Mat2VMat(m);

  Real rows    = VRows(vm);
  Real columns = VColumns(vm);
  Real nonNull = VNonNullCells(vm);

  Real density = nonNull/(rows*columns);

 // Set If(density <= 0.05, SVD(m,"Sparse"), SVD(m,"Jacobi"))
  Set If(density <= 0.01, SVD(m,"Sparse"), SVD(m))
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("It does the SVD decomposition choosing betwwen Sparse and
Jacobi method according to the density of de matrix m.
",
SVDDecomposition);
//////////////////////////////////////////////////////////////////////////////