// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
// ** Copyright UCAR (c) 1992 - 2007
// ** University Corporation for Atmospheric Research (UCAR)
// ** National Center for Atmospheric Research (NCAR)
// ** Research Applications Lab (RAL)
// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA
// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*

using namespace std;

#include <iostream>
#include <unistd.h>
#include <stdlib.h>
#include <string.h>
#include <cstdio>
#include <cmath>

#include "vx_contable/vx_contable.h"
#include "vx_util/vx_util.h"
#include "vx_met_util/vx_met_util.h"

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::baser(int &status) const {
   double v;

   v = oy_tp(status);

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::baser_ci(int &status, double alpha,
                                    double &cl, double &cu) const {
   double v;

   v = baser(status);

   compute_proportion_ci(v, n(), alpha, cl, cu);

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::fmean(int &status) const {
   double v;

   v = fy_tp(status);

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::fmean_ci(int &status, double alpha,
                                    double &cl, double &cu) const {
   double v;

   v = fmean(status);

   compute_proportion_ci(v, n(), alpha, cl, cu);

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::accuracy(int &status) const {
   double num, den, v;

   status = 0;

   num = (double) fy_oy() + fn_on();
   den = (double) n();

   if(den == 0.0) {
      status = 1;
      v = bad_data_double;
   }
   else {
      v = num/den;
   }

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::accuracy_ci(int &status, double alpha,
                                       double &cl, double &cu) const {
   double v;

   v = accuracy(status);

   compute_proportion_ci(v, n(), alpha, cl, cu);

   return(v);
}

////////////////////////////////////////////////////////////////////////
//
// Taken from eq. 7.9, page 241 in Wilks
//
////////////////////////////////////////////////////////////////////////

double TTContingencyTable::fbias(int &status) const {
   double num, den, v;

   status = 0;

   num = (double) fy();
   den = (double) oy();

   if(den == 0.0) {
      status = 1;
      v = bad_data_double;
   }
   else {
      v = num/den;
   }

   return(v);
}

////////////////////////////////////////////////////////////////////////
//
// Taken from eq. 7.7, page 240 in Wilks
//
////////////////////////////////////////////////////////////////////////

double TTContingencyTable::pod_yes(int &status) const {
   double num, den, v;

   status = 0;

   num = (double) fy_oy();
   den = (double) oy();

   if(den == 0) {
      status = 1;
      v = bad_data_double;
   }
   else {
      v = num/den;
   }

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::pod_yes_ci(int &status, double alpha,
                                      double &cl, double &cu) const {
   double v;

   v = pod_yes(status);

   compute_proportion_ci(v, n(), alpha, cl, cu);

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::pod_no(int &status) const {
   double num, den, v;

   status = 0;

   num = (double) fn_on();
   den = (double) on();

   if(den == 0.0) {
      status = 1;
      v = bad_data_double;
   }
   else {
      v = num/den;
   }

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::pod_no_ci(int &status, double alpha,
                                     double &cl, double &cu) const {
   double v;

   v = pod_no(status);

   compute_proportion_ci(v, n(), alpha, cl, cu);

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::pofd(int &status) const {
   double d, v;

   d = pod_no(status);

   if(status == 1) v = bad_data_double;
   else            v = 1.0 - d;

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::pofd_ci(int &status, double alpha,
                                   double &cl, double &cu) const {
   double v;

   v = pofd(status);

   compute_proportion_ci(v, n(), alpha, cl, cu);

   return(v);
}

////////////////////////////////////////////////////////////////////////
//
// Taken from eq. 7.8, page 241 in Wilks
//
////////////////////////////////////////////////////////////////////////

double TTContingencyTable::far(int &status) const {
   double num, den, v;

   status = 0;

   num = (double) fy_on();
   den = (double) fy();

   if(den == 0.0) {
      status = 1;
      v = bad_data_double;
   }
   else {
      v = num/den;
   }

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::far_ci(int &status, double alpha,
                                  double &cl, double &cu) const {
   double v;

   v = far(status);

   compute_proportion_ci(v, n(), alpha, cl, cu);

   return(v);
}

////////////////////////////////////////////////////////////////////////
//
// Taken from eq. 7.6, page 240 in Wilks
//
////////////////////////////////////////////////////////////////////////

double TTContingencyTable::csi(int &status) const {
   double num, den, v;

   status = 0;

   num = (double) fy_oy();
   den = (double) (n() - fn_on());

   if(den == 0.0) {
      status = 1;
      v = bad_data_double;
   }
   else {
      v = num/den;
   }

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::csi_ci(int &status, double alpha,
                                  double &cl, double &cu) const {
   double v;

   v = csi(status);

   compute_proportion_ci(v, n(), alpha, cl, cu);

   return(v);
}

////////////////////////////////////////////////////////////////////////
//
// Equations given by Barb Brown
//
////////////////////////////////////////////////////////////////////////

double TTContingencyTable::gss(int &status) const {
   long long a, b, c;
   double q, num, den, v;

   status = 0;

   a = (long long) fy_oy();
   b = (long long) fy_on();
   c = (long long) fn_oy();

   if(n() == 0) {
      status = 1;
      return(bad_data_double);
   }

   q   = (double) (a + b)*(a + c)/(n());
   num = (double) (a - q);
   den = (double) (a + b + c - q);

   if(den == 0.0) {
      status = 1;
      v = bad_data_double;
   }
   else {
      v = num/den;
   }

   return(v);
}

////////////////////////////////////////////////////////////////////////
//
// Taken from eq. 7.12, page 249 in Wilks
//
////////////////////////////////////////////////////////////////////////

double TTContingencyTable::hk(int &status) const {
   long long a, b, c, d;
   long long num, den;
   double v;

   status = 0;

   a = (long long) fy_oy();
   b = (long long) fy_on();
   c = (long long) fn_oy();
   d = (long long) fn_on();

   num = (a*d - b*c);
   den = (a + c)*(b + d);

   if(den == 0.0) {

      status = 1;
      return(bad_data_double);
   }

   v = (double) num/den;

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::hk_ci(int &status, double alpha,
                                 double &cl, double &cu) const {
   double v;

   v = hk(status);

   compute_hk_ci(v, alpha, fy_oy(), fy_on(), fn_oy(), fn_on(), cl, cu);

   return(v);
}

////////////////////////////////////////////////////////////////////////
//
// Taken from eq. 7.10, page 249 in Wilks
//
////////////////////////////////////////////////////////////////////////

double TTContingencyTable::hss(int &status) const {
   long long a, b, c, d;
   long long num, den;
   double v;

   status = 0;

   a = (long long) fy_oy();
   b = (long long) fy_on();
   c = (long long) fn_oy();
   d = (long long) fn_on();

   num = 2*(a*d - b*c);
   den = (a + c)*(c + d) + (a + b)*(b + d);

   if(den == 0.0) {
      status = 1;
      v = bad_data_double;
   }
   else {
      v = (double) num/den;
   }

   return(v);
}


////////////////////////////////////////////////////////////////////////
//
// Taken from eq. 7.10, page 249 in Wilks
//
////////////////////////////////////////////////////////////////////////

double TTContingencyTable::odds(int &status) const {
   int s;
   double num, den, v, py, pn;

   status = 0;

   py = pod_yes(s);
   pn = pofd(s);

   if(py == 1.0 || py == bad_data_double ||
      pn == 1.0 || pn == bad_data_double) {

      status = 1;
      return(bad_data_double);
   }

   num = py/(1 - py);
   den = pn/(1 - pn);

   if(den == 0.0) {

      status = 1;
      return(bad_data_double);
   }

   v = num/den;

   return(v);
}

////////////////////////////////////////////////////////////////////////

double TTContingencyTable::odds_ci(int &status, double alpha,
                                   double &cl, double &cu) const {
   double v;

   v = odds(status);

   compute_woolf_ci(v, alpha, fy_oy(), fy_on(), fn_oy(), fn_on(),
                    cl, cu);

   return(v);
}

////////////////////////////////////////////////////////////////////////