Ticket #347: roc.tol

NameBlock ROC =
[[
  Set _.interp.Set = SetOfText("linear", "akima", "cspline");

  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 ySet = EvalSet(interpSet, Real(Text code)
        {
          Code gslInterp = gsl_interp(code, xFunSet, yFunSet);
          gslInterp(0, x)
        });
        SetAvr(ySet)
      });
      SetOfReal(x)<<ySet
    })
  };

  Set Curve(Matrix y, Matrix py, Set sliceSet)
  {
    Set graf = EvalSet(sliceSet, Set(Real slice)
    {
      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 claTFP = Classify(graf, Real(Set reg1, Set reg2)
    {
      Real tfp1 = reg1[1];
      Real tfp2 = reg2[1];
      Compare(tfp1, tfp2)
    });
    EvalSet(claTFP, Set(Set class)
    {
      Real maxTvp = SetMax(Traspose(class)[2]);
      Select(class, Real(Set reg){EQ(reg[2], maxTvp)})[1]
    })
  };

  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)
  }
]];