Testeo de informacion a priori coherente: Comprobación de combinaciones lineales

[lmperez]

Comprobar si existe alguna combinación lineal entre parámetros de un nodo
observacional, teniendo en cuenta si éste tiene o no una distribución a
priori y/o está dentro de una jerarquía. Y además debe informar cuál es la
combinación lineal.

Os escribo debajo las siete funciones que hacen este proceso en bsrlayer,
una llama a las otras, pero usan sets de información que no existen en
mms, si decidís usarlas como guía preguntadme cualquier duda que tengáis.
Al final, si detecta una combinación lineal, escribe los parámetros que la
forman en un fichero log, en forma de query.

Una cosa a tener en cuenta. El proceso en general funcionaba en torno a la
función SVD de tol, sin embargo, como esta función no es exacta sino
aproximada, se probaron los diferentes métodos: Golub_Reinsch, Jacobi y
Sparse. De todos ellos, el más fiable fué el primero. El Sparse era muy
rápido cuando las matrices a factorizar tenían densidades pequeñas (<
0.01), y el Jacobi mostraba una salto de nivel muy claro entre los
autovalores de las variables que no son combinaciones lineales de los que
sí lo son. Sin embargo éstos dos últimos fallaban a veces, detectando
combinaciones lineales donde no las había o al revés. Por defecto se
decidió usar Golub_Reinsch con un redondeo en 10-10

//////////////////////////////////////////////////////////////////////////////
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);
//////////////////////////////////////////////////////////////////////////////

[lmperez]

Se ha vuelto a solicitar la necesidad de tener una herramienta para
testear los inputs: #122.