Ticket #63: def_strategy_vlogit.tol

//////////////////////////////////////////////////////////////////////////////
// FILE:    def_strategy_logit.tol
// PURPOSE: Define la estrategia de estimación asociada al método 
//          Logit de TOL que está basado en un método de estimación 
//          puntual máximo verosímil.
//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
Class @MMS.StrategyVLogit : @MMS.Strategy
//////////////////////////////////////////////////////////////////////////////
{
  // Especificaciones propias:
  //   (susceptibles de convertirse en clases)
  // * NameBlock de comunicación entre _BuildMatrices y _DoEstimation
  //   [[ Text name, Matrix outData, Matrix inpData ]]

  ////////////////////////////////////////////////////////////////////////////
  // Métodos no triviales:
  // * Execute()
  // ¿? Real SetEstimationIntervalMinMax()
  //
  // * Real _AfterSetModel()
  // * Set _BuildMatrices()
  // * Set _DoEstimation(Set AllModelDef)

  ////////////////////////////////////////////////////////////////////////////  
  Real SetEstimationIntervalMinMax(Real void)
  ////////////////////////////////////////////////////////////////////////////
  {
    Real 0
  };
  
  ////////////////////////////////////////////////////////////////////////////
  //  Prepara los parametros configurables para cada output a partir del
  //  modelo establecido en SetModel. Este metodo es invocado
  //  automaticamente desde el metodo SetModel implementado en la clase
  //  base.
  Real _AfterSetModel(Real void)
  ////////////////////////////////////////////////////////////////////////////
  {
    // (pgea) Considero interesante obviar el uso de esta función...
    Real void
  };

  ////////////////////////////////////////////////////////////////////////////
  // recorre todos los outputs en model y construye su ModelDef
  // correspondiente. Advierte de posible informacion contenida en
  // @MMS.Model que no es usada, por ejemplo informacion de jerarquia.

  ////////////////////////////////////////////////////////////////////////////
  Set _BuildMatrices(Real void)
  ////////////////////////////////////////////////////////////////////////////
  {
    @MMS.ModelAdapter model = @MMS.ModelAdapter::Default(_.model[1]);
    Set submodels_ = model::GetSubmodels(?); // active

    // Comprueba que son VLogit
    Set submodels = Select(submodels_, Real (@MMS.SubmodelAdapter submodel) {
      If(submodel::GetModelType(?)=="Logit", 1, { 
        WriteLn("El submodelo '"<<submodel::GetIndex(?)
          <<"' no es del tipo Logit.","E");
      0})
    });

    Set all_modeldef = EvalSet(submodels, NameBlock (@MMS.SubmodelAdapter submodel) {

      @MMS.OutputAdapter output = submodel::GetOutput(?);

      Text outputName = output::GetIndex(0);
      Matrix outputData = output::GetData(?);

      Set expTerms = submodel::GetExpTermsLinear(?); // active

      Set matrices = EvalSet(expTerms, Matrix (@MMS.ExpTermLinearAdapter expTerm) {
        expTerm::GetData(?)
      });
      Matrix inputData = BinGroup("|", matrices);

      //! Apaño para mantener el nombre de los parámetros...
      Set parameterNames = EvalSet(expTerms, Text(@MMS.ExpTermLinearAdapter expTerm) {
        expTerm::GetParameter(?)::GetIndex(?)
      });
      Set parameterIndices = For(1, Card(parameterNames), Real (Real i) {
        Eval("Real "<<parameterNames[i]<<" = Copy(i)")
      });

      NameBlock [[
        outputName,
        parameterIndices,
        outputData,
        inputData
      ]]    
    })
  };

//_.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-.
NameBlock Logit =
//_.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-.
[[
   Real _.ctInfo           = 0;
   Real _.dEpsilon         = DiffDist;
   Real _.maxIter          = MaxIter; 
   Real _.tolerance        = 10^(-8);
   Real _.tolerance.Rec    = 10^(-4);
   Real _.probability.Init = 1/2;
   Real _.num.Step         = 8;

   Matrix _Error(Matrix y, Matrix p)
   { 
     y-p //p = Probability(X, B) 
   };
   
   Matrix _Gradient(Matrix X, Matrix error)
   {
     Tra(Tra(error)*X) //error = y - Probability(X, B)
   };

   Matrix _Hessian(VMatrix vX, Matrix p)
   {
     VMatrix W = Mat2VMat(p$*RSum(-p,1));
     VMatrix V = Eye(Rows(p), Rows(p), 0, W);
     VMat2Mat(-Tra(V*vX)*vX)     
   };

   Matrix _Probability(Matrix X, Matrix B, Matrix LogitGroupProb)
   {
     Matrix p = RPow(RSum(Exp(-X*B-LogitGroupProb), 1), -1);
     p
   };

   Matrix _MLnLikelyhood(Matrix y, Matrix p)
   {
     y$*Log(p)+RSum(-y,1)$*Log(RSum(-p,1))
   };

   Matrix _MLikelyhood(Matrix mLnLikelyhood)
   {
     Exp(mLnLikelyhood)
   };

   Set Estimate.MaxLikelyhood
   (
     Matrix y,        // Matriz y de variable salida binaria
     Matrix XIni,     // Matriz x de variables entrada
     Matrix B0Ini,    // Matriz de parametros iniciales
     Anything ctInfo, // Un real 0 sin constante, 
                      // Un real 1 estima constante,
                      // Una matriz  
     Real dEpsilon,   // Diferencia de paso
     Real maxIter,    // Numero maximo de iteraciones
     Real tolerance   // Tolerancia al error
   )
   {
     Text grCtInfo = Grammar(ctInfo);
     
     Matrix X = Case
     (
/*
       And(grCtInfo=="Real", EQ(ctInfo, 0)), XIni,
       And(grCtInfo=="Real", EQ(ctInfo, 1)), 
       {
         XIni|Rand(Rows(XIni), 1, 1, 1)
       },
*/
       grCtInfo=="Matrix", XIni,
       1, XIni
     );

     Matrix IdCt = Case
     (
//       And(grCtInfo=="Real", EQ(ctInfo, 0)), Rand(Rows(XIni), 1, 0, 0),
//       And(grCtInfo=="Real", EQ(ctInfo, 1)), Rand(Rows(XIni), 1, 0, 0),
       grCtInfo=="Matrix", ctInfo,
       1, Rand(Rows(XIni), 1, 0, 0)
     );
        
     Matrix B0 = Case
     (
       grCtInfo=="Matrix", B0Ini,
       1, B0Ini
     );
     
     Real n = Columns(X); // Number of variables
     Real N = Rows(X);    // Number of data

     Text iniMsg = 
      "Model Logit Init ("+IntText(N)+"x"+IntText(n)+")"+Time;
     Real WriteLn(iniMsg);  

     VMatrix vX       = Mat2VMat(X);

     Matrix zeroMat   = Rand(N, 1, 0, 0);
     Matrix oneMat    = Rand(N, 1, 1, 1);
     Matrix negInfMat = Rand(N, 1, -1/0, -1/0);  
     Matrix tolMat    = Rand(N, 1, tolerance, tolerance);   
     
     Real completeTime    = 0;
     Real difTimeMin      = 0;
     Real exit            = 0;
     Matrix B             = Copy(B0);
     Real iter            = 0;
     Real oldLnLikelyhood = 0; 
   
     Set cycle = Empty;
     
     Real While(Not(exit), 
     {
       Real time0 = Copy(Time);

//Real CMsg::Trace::show(1, "Printing oldLnLikelyhood="<<oldLnLikelyhood);
       Matrix p     = _Probability(X, B, IdCt); 

//Real CMsg::Trace::show(1, "Printing p="<<p);
       Matrix error = _Error(y, p);

       Matrix G     = _Gradient(X, error);
       Real time0.1 = Copy(Time);
       Matrix H     = _Hessian(vX, p);
       Real time0.2 = Copy(Time);  
       Matrix dif   = MinimumResidualsSolve(H, G);
 
       Real norm      = MatFrobeniusNorm(G);
       Real advance   = MatFrobeniusNorm(dif);
       Real maxAbsDif = MatMax(Abs(dif));
       

       Matrix mLnLikelyhood = _MLnLikelyhood(y, p);
       Matrix mLLCorrect    = IfMat(EQ(Abs(error), oneMat), negInfMat,
       IfMat(LT(Abs(error), tolMat), zeroMat, mLnLikelyhood));

       Real lnLikelyhood     = MatSum(mLLCorrect);
//Real CMsg::Trace::show(1, "Printing dif"<<dif);
//Real CMsg::Trace::show(1, "Printing mLnLikelyhood"<<mLnLikelyhood);
//Real CMsg::Trace::show(1, "Printing mLLCorrect"<<mLLCorrect);
//Real CMsg::Trace::show(1, "Printing error"<<error);

//Real CMsg::Trace::show(1, "Printing lnLikelyhood="<<lnLikelyhood);

       Real difTimeMin := (time0.2-time0.1)+difTimeMin;
       
       Real difTime = Copy(Time)-time0; 
       Real completeTime := completeTime+difTime;
       Real exitTreatment = Case
       (  
         IsUnknown(advance),
         { 
           Real CMsg::Trace::show(1, "Advance is unknown");
           Real exit:=1 
         },
         GE(iter, maxIter), 
         { 
           Real CMsg::Trace::show(1, "MaxIter reached"); 
           Real exit:=1
         },
         LT(norm, dEpsilon),
         {
           Real CMsg::Trace::show(1, "Norm lower than "<<dEpsilon); 
           Real exit:=1
         },
         LT(maxAbsDif, tolerance), 
         { 
           Real CMsg::Trace::show(1, "MaxAbsDif lower than "<<tolerance); 
           Real exit:=1
         },
         LT(Abs(oldLnLikelyhood - lnLikelyhood), dEpsilon), 
         { 
           Real CMsg::Trace::show(1, 
           "LnLikelihood is stable: lnLikelyhood    ="<<lnLikelyhood+
                                  " oldLnLikelyhood ="<<oldLnLikelyhood);
           Real exit:=1
         },         
         1,
         {
           Text exitMsg = " Logit model iteration("+
            FormatReal(iter, "%0"+IntText(Floor(Log10(maxIter)+1))+".lf")+
           ")"+"        "+
           " LogLikelyhood = "+FormatReal(lnLikelyhood, "%.4E")+"       "+
           " MaxAbsDif = "+FormatReal(maxAbsDif, "%.4E")+"      "+
           " Gradient Norm = "+FormatReal(norm, "%.4E")+"       "+
           " Time = "+FormatReal(difTime, "%.4lf");
           Real CMsg::Trace::show(1, exitMsg); 
           Real iter := iter+1;   
           Real oldLnLikelyhood := Copy(lnLikelyhood);
           Matrix B:=(Matrix B-dif);
           Real 0           
         } 
       );       
       Real exitCycle = If(EQ(exitTreatment, 1), 
       {
         Set cycle := SetOfAnything
         (B, p, G, H, dif, norm, advance, 
          maxAbsDif, error, lnLikelyhood, oldLnLikelyhood, ctInfo, X, XIni, y,
          mLLCorrect);
         1
       }, 0)
     }); 
     Real time2     = Copy(Time); 
/*
     Matrix p       = _Probability(X, B); 

     Matrix error   = _Error(y, p);

     Matrix G       = _Gradient(X, error);
     Matrix H       = _Hessian(vX, p);
     Matrix dif     = MinimumResidualsSolve(H, G);

     Real norm            = MatFrobeniusNorm(G);
     Real advance         = MatFrobeniusNorm(dif);
     Real maxAbsDif       = MatMax(Abs(dif));
     Matrix mLnLikelyhood = _MLnLikelyhood(y, p);
     Real lnLikelyhood    = MatSum(mLnLikelyhood);
*/
     Real difTime2 = Copy(Time)-time2;  
     Real totalTime = completeTime+difTime2;
     Text endMsg = 
      "Model Logit Ended. Time:"<<totalTime+" SplitTime:"<<difTimeMin+NL+NL;
     Real WriteLn(endMsg);  

/*
     SetOfAnything
     (B, p, G, H, dif, norm, advance, 
      maxAbsDif, error, lnLikelyhood, oldLnLikelyhood)
*/
     cycle
   };

   Set Estimate.MaxLikelyhood.Default(Matrix y, Matrix X)
   { 
     Matrix B0Ini = Rand(Columns(X), 1, 0, 0);
     Estimate.MaxLikelyhood
     (
       y,      // Matriz y de variable salida binaria
       X,      // Matriz x de variables entrada
       B0Ini,
       _.ctInfo,
       _.dEpsilon, // Diferencia de paso
       _.maxIter,  // Numero maximo de iteraciones
       _.tolerance // Tolerancia al error
     )
   };
   Set Estimate.MaxLikelyhood.Constant(Matrix y, Matrix X, Matrix constant)
   { 
     Matrix B0Ini = Rand(Columns(X), 1, 0, 0); 
     Estimate.MaxLikelyhood
     (
       y,      // Matriz y de variable salida binaria
       X,      // Matriz x de variables entrada
       B0Ini,
       constant,
       _.dEpsilon, // Diferencia de paso
       _.maxIter,  // Numero maximo de iteraciones
       _.tolerance // Tolerancia al error
     )
   };
   Set Estimate.MaxLikelyhood.ProbB0
   (
     Matrix y, 
     Matrix X, 
     Matrix B0Ini, 
     Real prob
   )
   { 
     Matrix probMat  = Rand(Rows(y), 1, prob, prob);
     Matrix constant = Log(probMat$/RSum(-probMat, 1));

     Estimate.MaxLikelyhood
     (
       y,      // Matriz y de variable salida binaria
       X,      // Matriz x de variables entrada
       B0Ini, 
       constant,
       _.dEpsilon, // Diferencia de paso
       _.maxIter,  // Numero maximo de iteraciones
       _.tolerance // Tolerancia al error
     )
   };

   Set Estimate.MaxLikelyhood.ProbRec
   (
     Matrix y, 
     Matrix X, 
     Matrix B0Cur,
     Real step,
     Real tolerance,
     Real probCur,
     Real probEnd
   )
   { 
     Text iniMsg = NL+NL+"Model Logit (p="<<probCur+") "+Time;
     Real WriteLn(iniMsg);  

     Set logitResult = Estimate.MaxLikelyhood.ProbB0(y, X, B0Cur, probCur);
     If(LT(Abs(probCur - probEnd), tolerance), logitResult,
     {
       Real prob = probCur*step;
       Matrix B0 = logitResult::B;
       Estimate.MaxLikelyhood.ProbRec 
       (y, X, B0, step, tolerance, prob, probEnd)
     })
   };

   Set Estimate.MaxLikelyhood.Prob
   (
     Matrix y, 
     Matrix X, 
     Real prob
   )
   {
     Real step    = (prob/_.probability.Init)^(1/_.num.Step);
     Matrix B0Ini = Rand(Columns(X), 1, 0, 0); 
     Estimate.MaxLikelyhood.ProbRec
     (y, X, B0Ini, step, _.tolerance.Rec, _.probability.Init, prob)
   };
   
   Set Diagnosis
   ( 
     Matrix y,
     Matrix X,
     Matrix B,
     Matrix error,
     Real lnLikelyhood,
     Set names,
     Matrix H  //Hessian
   ) 
   {
     Text iniMsg = 
      "Model Logit Diagnosis. Init Time:"+Time;
     Real WriteLn(iniMsg); 

     Real N = Rows(X);    //
     Real n = Columns(X); // Parameter number

     Matrix FIM = -H; // Fisher Information Matrix
     Matrix COV = SVDInverse(FIM); // Varianze Covarianza Parameter Matrix
     
     Set Parameters = For(1, Columns(X), Set(Real k)
     {
       Text name  = names[k];
       Real value = MatDat(B, k, 1);
       Real var   = MatDat(COV, k, k);
       Real std   = SqRt(var);
       Real tStudent = value/std;
       Real refProb  = 2*(1-DistT(Abs(tStudent), N-n-1));
       ParameterInf
       (
         name,     //Name
         0,        //Factor
         0,        //Order
         value,    //Value
         std,      //StDs
         tStudent, //TStudent
         refProb   //RefuseProb
       )
     });
  
     Real VarTot    = MatVar(y);
     Real VarError  = MatVar(error);
     Real R2 = 1-VarError/VarTot;

     Real MaxProb   = MatSum(y)/N;
     Real lnLikelyhoodIntercep = 
                      MatSum(_MLnLikelyhood(y, Rand(N, 1, MaxProb, MaxProb)));

     Real Nagelkerke.R2    = 1-(Exp((2/N)*(lnLikelyhoodIntercep-lnLikelyhood)));
     Real Nagelkerke.R2Max = 1-Exp((2/N)*lnLikelyhoodIntercep); 
     Real Nagelkerke.R2MaxRescaled = Nagelkerke.R2/Nagelkerke.R2Max;    
     Text endMsg = 
      "Model Logit Diagnosis. End Time:"+Time+NL;
     Real WriteLn(endMsg); 
 
     SetOfAnything
     (
       Parameters,
       FIM,
       COV,
       MaxProb,
       R2,
       Nagelkerke.R2,
       Nagelkerke.R2Max,
       Nagelkerke.R2MaxRescaled,
       lnLikelyhood, 
       lnLikelyhoodIntercep
     )
   };

  Set PreTesting(Matrix Y, Matrix X, Set varNames)
  {
    Text iniMsg = 
     "Model Logit PreTesting. Init Time:"+Time;
    Real WriteLn(iniMsg);
 
    Real WriteLn("  Checking column stability..."+Time);
    Matrix unkX    = IsUnknown(X);
    Matrix posInfX = IsPosInf(X);
    Matrix negInfX = IsNegInf(X);
  
    Matrix unkY    = IsUnknown(Y);
    Matrix posInfY = IsPosInf(Y);
    Matrix negInfY = IsNegInf(Y);
  
    Real isUnkX    = MatSum(unkX);
    Real isPosInfX = MatSum(posInfX);
    Real isNegInfX = MatSum(negInfX);
  
    Real isUnkY    = MatSum(unkY);
    Real isPosInfY = MatSum(posInfY);
    Real isNegInfY = MatSum(negInfY);
  
    Real n = Rows(Y);
    Real balanced = MatSum(Y)/n;
  
    Real valid = 
     Not(Or(isUnkX, isPosInfX, isNegInfX, isUnkY, isPosInfY, isNegInfY));

    Set checkValid = If(EQ(valid, 1), Empty,
    {
      Set data = For(1, Card(varNames), Set(Real k)
      {
        Text name = varNames[k];
        Real kUnkX    = MatSum(SubCol(unkX, [[k]]));  
        Real kPosInfX = MatSum(SubCol(posInfX, [[k]]));
        Real kNegInfX = MatSum(SubCol(negInfX, [[k]]));
        SetOfAnything(name, kUnkX, kPosInfX, kNegInfX)
      });
      Set header = SetOfText("VarName", "Unk", "PosInf", "NegInf");
      SetOfSet(header)<<
      SetOfSet(SetOfAnything("Y", isUnkY, isPosInfY, isNegInfY))<<
      data 
    });
    Real WriteLn("  Adjusting X matrix..."+Time);

    Matrix YPre = 
     If(isUnkY, IfMat(unkY, VMat2Mat(Eye(Rows(Y), 1, 0, 0)), Y), Y);
    Matrix XPre = 
     If(isUnkX, IfMat(unkX, VMat2Mat(Eye(Rows(X), Columns(X), 0, 0)), X), X);    
     
    VMatrix vXPre = Mat2VMat(XPre);
    Set index    = Range(1, Columns(X), 1); 
//    Real WriteLn("  Correlation Y|X matrix..."+Time);
//    Matrix corVarX = Cor(Tra(YPre|XPre));

    Real WriteLn("  X information column..."+Time);  
    Set ColInfo = For(1, Card(varNames), Set(Real k)
    {
      Text name  = varNames[k];
      Matrix col = SubCol(XPre, [[k]]);
      Real min   = MatMin(col);
      Real max   = MatMax(col);
      Real stds  = MatStDs(col); 
      Real avr   = MatAvr(col);
      Matrix freq   = Frequency(col, 100, min, max);
      Matrix freq01 = Frequency(col$*Y, 100, min, max); 
      Matrix ratioDisc = freq01$/freq; 
      Real minValue    = MatDat(freq01, 1, 2);
      Real maxValue    = MatDat(freq01, 100, 2); 
      Real maxRatio    = MatDat(ratioDisc, 100, 2)/balanced; 
      Real dicotomicRatio = (minValue+maxValue)/n ; 

      VMatrix y = SubCol(vXPre, [[k]]);
      VMatrix x = SubCol(vXPre, index-[[k]]);

      Set linReg   = LinReg::Get.GeneralInformation(y, x);
      Real R2MultiColinearity = linReg::R2;
  
      SetOfAnything
      (
        name, 
        min, 
        max, 
        stds, 
        avr, 
        freq, 
        freq01, 
        ratioDisc, 
        maxRatio, 
        maxValue,
        dicotomicRatio,
        R2MultiColinearity
      )
    });

    Text endMsg = 
     "Model Logit PreTesting. End Time:"+Time+NL;
    Real WriteLn(endMsg); 

    SetOfAnything
    (unkX, posInfX, negInfX, unkY, posInfY, negInfY, valid, checkValid, 
     balanced, /*corVarX,*/
     ColInfo, YPre, XPre)
  }
]];
//_.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-.


  ////////////////////////////////////////////////////////////////////////////
  //  Aplica el metodo Logit a todos los outputs.
  Set _DoEstimation(Set allOutMat)
  ////////////////////////////////////////////////////////////////////////////
  {
    Set EvalSet(allOutMat, Set(NameBlock args) {
      Set results = Logit::Estimate.MaxLikelyhood.Default(
        args::outputData, 
        args::inputData
      );

      Set resultsMod = [[
        Matrix Parameters = results["B"];
        Matrix Residuals = results["error"];
        Matrix Probabilities = results["p"];
        Matrix Gradient = results["G"];
        Matrix Hessian = results["H"];
        Matrix Likelihoods = results["mLLCorrect"];
        Real Likelihood = Exp(results["lnLikelyhood"]);
        Real LogLikelihood = results["lnLikelyhood"];
        Set ParameterIndices = args::parameterIndices 
      ]];
      Set PutName(args::outputName, resultsMod)
    })
  };

  ////////////////////////////////////////////////////////////////////////////
  @MMS.ResultsAdapter Execute(Real void)
  ////////////////////////////////////////////////////////////////////////////
  {
    Set argsEstim = _BuildMatrices(0);
    Real ApplySettings(0);
    Set _results = _DoEstimation(argsEstim);
    Real RestoreSettings(0);
    @MMS.ResultsAdapterLogit resultsAdapter = 
      @MMS.ResultsAdapterLogit::New(_results)
  };

  ////////////////////////////////////////////////////////////////////////////
  Static @MMS.StrategyVLogit New(NameBlock args)
  ////////////////////////////////////////////////////////////////////////////
  {
    @MMS.StrategyVLogit obj = [[
      // ¿necesita nombre y descripción?
      Text _.name = getOptArg(args, "_.name", "VLogit");
      Text _.description = getOptArg(args, "_.description", "")
    ]];
    // Settings de Logit construidas de acuerdo al ticket de TOL-Project #736
    Real obj::AddSetting.FromObject(MaxIter);
    Real obj::AddSetting.FromObject(Tolerance);
    obj
  };

  ////////////////////////////////////////////////////////////////////////////
  Static @MMS.StrategyVLogit Default(Real void)
  ////////////////////////////////////////////////////////////////////////////
  {
    @MMS.StrategyVLogit::New([[ Text _.name = "VLogit" ]])
  }
};