Ticket #79: funLinCom.tol

//////////////////////////////////////////////////////////////////////////////
Real CheckModel(Set algSetNode) 
//////////////////////////////////////////////////////////////////////////////
{
  Text nameOutput = GetNameModel(algSetNode);
  @MMS.Model model = MMS::GetModel([[nameOutput,GetVersion(?)]]);
  Set outputs = model::GetOutputs(?);
  Set EvalSet(outputs, Real (NameBlock output){
    Serie outputSer = output::GetData(?);
    Date iniOutput = Max(First(outputSer),output::GetModelBegin(?));
    Date endOutput = Min(Last(outputSer),output::GetModelEnd(?));
    Set setExpTerm = output::GetExpTerms(?);
    Set EvalSet(setExpTerm, Real (NameBlock ET){
      Serie serExpTerm = ET::GetData(?);
      Serie subSerExpTerm = SubSer(serExpTerm,iniOutput,endOutput);
      If(And(
           AvrS(subSerExpTerm)==0,
           StDsS(subSerExpTerm)==0,
           CheckPriorExpTerm(model,ET)==0,
           CheckHierExpTerm(model,ET)==0
           )
        ,
        ET::SetIsActive(0)
        ,
        ET::SetIsActive(1)
      );
      Real 1
    });
 
    Real 1
  });

  Real 1
        
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("Chequea el modelo. Quita los inputs que son cero", 
CheckModel);
//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
Real CheckPriorExpTerm(NameBlock model,NameBlock ET) 
//////////////////////////////////////////////////////////////////////////////
{
  Set priors = model::GetPriors(?);
  If(IsEmpty(priors),0,
    {
    Set hasprior = EvalSet(priors, Real (NameBlock prior){
      //Set setElements = prior::GetElements(?);
      //Set priorElements = EvalSet(setElements, Real (NameBlock element){
        //Real (prior::GetElement::GetParent(?)::GetName(?)==ET::GetName(?))
  
      //}); 
      Real If(prior::GetElement(?)::GetParent(?)::GetName(?)==ET::GetName(?),
        1,0)
    });
  
    SetSum(hasprior)
    }
  )
        
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("Chequea el modelo. 
Comprueba si un termino explicativo tiene priors", 
CheckPriorExpTerm);
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Real CheckHierExpTerm(NameBlock model,NameBlock ET) 
//////////////////////////////////////////////////////////////////////////////
{
  Set hierarchies = model::GetHierarchies(?);

  If(IsEmpty(hierarchies),0,
    {
    Set hashierarchie = EvalSet(hierarchies, Real (NameBlock hierarchie){
      Set setElements = hierarchie::GetElements(?);
      Set hierarchieElements = EvalSet(setElements, Real (NameBlock element){
        Real (element::GetParent(?)::GetName(?)==ET::GetName(?))
  
      }); 
      SetSum(hierarchieElements)
    });
  
    SetSum(hashierarchie)
    }
  )     
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("Chequea el modelo. 
Comprueba si un termino explicativo tiene jerarquias", 
CheckHierExpTerm);
//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
Real LinComWarning(NameBlock model)
//////////////////////////////////////////////////////////////////////////////
{
  WriteLn("Entrando en lincomwarning");
  Set outputs = model::GetOutputs(?);

  WriteLn("Cogiendo outputs en lincomwarning");

  Real SetProd(EvalSet(outputs,Real (NameBlock output)
  {
    // Para cada output definido en el modelo, sacamos sus fechas de domino, 
    // diferencia y terminos explicativos
    Serie outputSer = output::GetData(?);
    WriteLn("Cogiendo datos de outputs en lincomwarning");

    Date iniCalc = Max(First(outputSer),output::GetModelBegin(?));
    Date endCalc = Min(Last(outputSer),output::GetModelEnd(?));
    Polyn dif    = output::GetARIMA(?)::GetDif(?);
    Set expTerms = model::GetExpTerms(?);

    // Para cada término explicativo, miro que sea no nulo, que no tenga prior
    // y en tal caso devuelvo la serie input, si no lo eliminamos del calculo
    Set series = EvalSet(expTerms,Serie (NameBlock expTerm)
    {
     
      Serie serExpTerm = expTerm::GetData(?);
    WriteLn("Cogiendo datos de terminos explicativos en lincomwarning");
          
      Serie subSerExpTerm = SubSer(serExpTerm,iniCalc,endCalc);
      Serie Case(
        And(Or(
           CheckPriorExpTerm(model,expTerm)>=1,
           CheckHierExpTerm(model,expTerm)>=1
           ),Not(MaxS(subSerExpTerm)==1/0))
        ,
        {
          expTerm::SetIsActive(1);
          UnknownSerie
        },
        Or(And(MaxS(subSerExpTerm)==0,
           MinS(subSerExpTerm)==0),MaxS(subSerExpTerm)==1/0),
        {
          expTerm::SetIsActive(0);
          UnknownSerie
        }
        ,
        1
        ,
        {
          expTerm::SetIsActive(1);
          Eval("Serie "+expTerm::GetName(?)+" = subSerExpTerm")
        }
      )

    }) - [[ Serie UnknownSerie ]];


    // Paso las series input a matriz
    Matrix inputs = SerSetMat(series);
    Matrix matrix = inputs*Tra(inputs);
    Matrix inputs:= UnknownMatrix;//--> Ahorro de memoria en matrices grandes

    // Antes que SVD, llamo a GaussInverse, mucho mas rapida, si devuelve 
    // una matriz de 0x0 entonces es que no tiene inversa (hay combinacion 
    // lineal), y entonces llamo a SVD. Está escrito asi para ahorrar memoria

    Set svd  = {
      Matrix gaussInv = GaussInverse(matrix);
      Set svd = If( And(EQ(Rows(gaussInv),0),EQ(Columns(gaussInv),0)),
        //SVD(matrix),
        RoundSVD(SVD(matrix),False,0.5,10^(-12)),
        Empty
      );
      Matrix matrix:= UnknownMatrix;//--> Ahorro de memoria en matrices grandes
      Set svd
    };

    Real If(Not(Card(svd)), Real False,
    {
      // Saco la última columna de la segunda matriz de autovectores
      Matrix  lastCol = SubCol(svd[3],[[ Columns(svd[3]) ]]);
      Set svd := Empty;//--> Ahorro de memoria en matrices grandes

      // Miro los elementos no nulos, que son las ecuaciones de los 
      // parametros que hacen c.l.
      Matrix  LogicLC = Not(Not(lastCol));

      // Transformamos LogicLC en set para buscar
      Set indices = BinGroup("<<",MatSet(LogicLC));

      // Selecionamos las series con indice no nulo
      Real c = 0;
      Set lcw = EvalSet(series,Text (Serie x)
      {
        Real c := c+1;
        Text If(indices[c], Name(x), "")
      }) - [[ "" ]];

      Text txt2File = TxtListItemQuote(lcw,","+NL); 

      // Sacamos los input que tienen combinacion lineal en traza y en fichero
      // .log, por si lo ejecutas en remoto sin tolbase
      Text WriteLn("[LinComWarning] THERE ARE "+IntText(Card(lcw))+" LINEAR"
        +"COMBINATIONS OF INPUT SERIES CHECK LOG DIRECTORY!!"+NL+NL+NL+
        txt2File+NL+NL+NL,
        Text "E");

      Set EvalSet(lcw,Real (Text nameParameter){
        NameBlock expTermPrior = model::GetExpTerm(nameParameter);
        Set parameterPrior = expTermPrior::GetParameters(?);
        Set EvalSet(parameterPrior, Real (NameBlock parameterPrior)
        {
           Text parameterName = parameterPrior::GetIndex(?);
           Real model::CreatePrior(
            [[
              Text _.element = parameterName;
              Real _.average = 0;
              Real _.sigma = 10
            ]]
            )
        });
        Real 1
      });

      Text logRoute = PathProject+"log/lincom/lincom_"+
        model::GetName(?)+
        FormatDate(Now) +".log";

      Text WriteFile(logRoute,txt2File);

      Real 1//Stop
    })
  }))
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("Combprueba si existen combinaciones lineales en la estrategia
 que recibe como argumento",
LinComWarning);
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Set RoundSVD(Set svd,Real trasposed, Real minSparse, Real tolerance)
//////////////////////////////////////////////////////////////////////////////
{/*
Set svd = SVD(matrix);
Real trasposed = False;
Real minSparse = 0.5;
Real tolerance = 10^(-12);
*/
  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) by default
",
RoundSVD);
//////////////////////////////////////////////////////////////////////////////