Ticket #1350: roc_v2.tol

NameBlock ROC =
[[
  Set _.interp.Set = SetOfText("linear", "cspline", "akima");
  Set _.slice.Set  = Range(0, 1, 0.02);

  Set _Interp
  (
    Set xSet,      // Set xSet with we want to know y for all funSet
    Set setFunSet, // Set of tables with same structure SetOfReal(x, y)
    Set interpSet  // Set with gsl valid interpolate method (see gsl_interp)
  )
  {
    Set setxyFunSet = EvalSet(setFunSet, Set(Set funSet)
    {
      SetOfMatrix(SetCol(Traspose(funSet)[1]), SetCol(Traspose(funSet)[2])) 
    });    
  
    EvalSet(xSet, Set(Real x)
    {
      Set ySet = EvalSet(setxyFunSet, Real(Set xyFunSet)
      {
        Matrix xFunSet = xyFunSet[1];
        Matrix yFunSet = xyFunSet[2];

        Set ySetAux = EvalSet(interpSet, Real(Text code)
        {
          Code gslInterp = gsl_interp(code, xFunSet, yFunSet);
          gslInterp(0, x)
        });
        SetAvr(ySetAux)
      });
      SetOfReal(x)<<ySet
    })
  };

  Set Curve(Matrix y, Matrix py_, Real numQ)
  {
    Matrix py = Min(Max(py_, 0), 1);
    Set index = Range(1, 0, -1/Min(Max(numQ, 2), Rows(py)-1));
    Set graf = [[ [[ 0, 0, 0 ]] ]] << EvalSet(index, Set(Real q)
    {
      Real slice = MatQuantile(py, q);
      Matrix yEst = GE(py, Rand(Rows(py), 1 , slice, slice));
      Real VP     = MatSum(And(yEst, y));
      Real FN     = MatSum(And(Not(yEst), y));     
      Real VN     = MatSum(And(Not(yEst), Not(y)));
      Real FP     = MatSum(And(yEst, Not(y)));
  
      Real TVP    = VP/(VP+FN);
      Real TFP    = 1-VN/(VN+FP);//=FP/(VN+FP)  
      SetOfReal(TFP, TVP, slice)
    })
  };

  Set Curve.Abs
  (
    Set ROCSet, // Table with reg = SetOfReal(TFP, TVP, slice, ...)
    Real M, // Total
    Real m  // Infected
  )
  {
    Real p = m/(M-m);
    EvalSet(ROCSet, Set(Set reg)
    {
      Real TFP        = reg[1];
      Real TVP        = reg[2];
      Real slice      = reg[3]; 
      Real infected   = TVP*m;
      Real population = (p*TVP+TFP)*M/(1+p);
      SetOfReal(TFP, TVP, population, infected, slice)
    })
  };

  Set Eval.TFPSet(Set tfpSet, Set ROCSet)
  {
    Set tROCSet = Traspose(ROCSet); 
    Set txFP    = tROCSet[1];
    Set txVP    = tROCSet[2]; 
    Set txSL    = tROCSet[3]; 
    
    Set setFunSet = SetOfSet
    (
      Traspose(SetOfSet(txFP, txVP)), 
      Traspose(SetOfSet(txFP, txSL))
    );
    _Interp(tfpSet, setFunSet, _.interp.Set)
  };

  Set EvalAbs.PopSet(Set popSet, Set ROCAbsSet)
  {
    Set tROCSet = Traspose(ROCAbsSet); 
    Set txPop    = tROCSet[3];
    Set txInf    = tROCSet[4]; 
    Set txSL     = tROCSet[5]; 
    
    Set setFunSet = SetOfSet
    (
      Traspose(SetOfSet(txPop, txInf)), 
      Traspose(SetOfSet(txPop, txSL))
    );
    _Interp(popSet, setFunSet, _.interp.Set)
  };

  Real Get.Area(Set ROCSet)
  {
    Set txFP     = Traspose(ROCSet)[1];
    Set txVP     = Traspose(ROCSet)[2]; 
    Set setFunSet = SetOfSet( Traspose(SetOfSet(txFP, txVP)) );
    Real get.tvp(Real tfp)
    {
      Set tfpSet = SetOfReal(tfp);
      Real tvp = _Interp(tfpSet, setFunSet, _.interp.Set)[1][2];
      tvp
    };
    IntegrateQAG(get.tvp, 0, 1)
  };

  Real Get.AreaTrap(Set xySet)
  {
    Real auc = 0;
    Real k = 2;
    Real While(LE(k,Card(xySet)),
    {
      Real xk   = xySet[k][1];
      Real xk_1 = xySet[k-1][1]; 

      Real yk = xySet[k][2];
      Real yk_1 = xySet[k-1][2];

      Real auc:= (xk-xk_1)*(yk+yk_1)/2 + auc;
      Real k := k+1
    });
    auc
  };

  Real Get.AreaFull(Set ROCSet)
  {
    Set txFP     = Traspose(ROCSet)[1];
    Set txVP     = Traspose(ROCSet)[2]; 
    Set setFunSet = SetOfSet( Traspose(SetOfSet(txFP, txVP)) );
    Real get.tvp(Real tfp)
    {
      Set tfpSet = SetOfReal(tfp);
      Real tvp = _Interp(tfpSet, setFunSet, _.interp.Set)[1][2];
      tvp
    };
    IntegrateQAG(get.tvp, 0, 1)
  }
]];