diff --git a/internal/scripts/docker/Dockerfile.sonarqube b/internal/scripts/docker/Dockerfile.sonarqube index e1e935d857..7a6a7114f5 100644 --- a/internal/scripts/docker/Dockerfile.sonarqube +++ b/internal/scripts/docker/Dockerfile.sonarqube @@ -9,7 +9,7 @@ MAINTAINER John Halley Gotway # SonarQube static code analysis on the specified branch or tag. # https://docs.sonarqube.org/latest/analysis/scan/sonarscanner/ # -ARG SONAR_SCANNER_VERSION=5.0.1.3006 +ARG SONAR_SCANNER_VERSION=6.1.0.4477 ARG SONAR_HOST_URL ARG SONAR_TOKEN ARG SOURCE_BRANCH diff --git a/internal/scripts/docker/build_met_sonarqube.sh b/internal/scripts/docker/build_met_sonarqube.sh index 2eea05654e..78163b5379 100755 --- a/internal/scripts/docker/build_met_sonarqube.sh +++ b/internal/scripts/docker/build_met_sonarqube.sh @@ -100,6 +100,7 @@ time_command ./configure \ BUFRLIB_NAME=${BUFRLIB_NAME} \ GRIB2CLIB_NAME=${GRIB2CLIB_NAME} \ --enable-all \ + MET_CXX_STANDARD=11 \ CPPFLAGS="-I/usr/local/include -I/usr/local/include/freetype2 -I/usr/local/include/cairo" \ LIBS="-ltirpc" diff --git a/internal/scripts/sonarqube/run_sonarqube.sh b/internal/scripts/sonarqube/run_sonarqube.sh index ac439b8a11..9b3557581a 100755 --- a/internal/scripts/sonarqube/run_sonarqube.sh +++ b/internal/scripts/sonarqube/run_sonarqube.sh @@ -109,8 +109,9 @@ run_command "git checkout ${1}" # Otherwise, the SonarQube logic does not work. export MET_DEVELOPMENT=true -# Run the configure script -run_command "./configure --prefix=`pwd` --enable-all" +# Run the configure script. +# Specify the C++ standard to limit the scope of the findings. +run_command "./configure --prefix=`pwd` --enable-all MET_CXX_STANDARD=11" # Define the version string SONAR_PROJECT_VERSION=$(grep "^version" docs/conf.py | cut -d'=' -f2 | tr -d "\'\" ") diff --git a/internal/test_unit/python/unit.py b/internal/test_unit/python/unit.py index 831a540b32..3bc2b222a9 100755 --- a/internal/test_unit/python/unit.py +++ b/internal/test_unit/python/unit.py @@ -124,7 +124,7 @@ def unit(test_xml, file_log=None, cmd_only=False, noexit=False, memchk=False, ca if cmd_only: if 'env' in test.keys(): for key, val in sorted(test['env'].items()): - print(f"export {key}={val}") + print(f"export '{key}={val}'") print(f"{cmd}") if 'env' in test.keys(): for key, val in sorted(test['env'].items()): @@ -143,7 +143,7 @@ def unit(test_xml, file_log=None, cmd_only=False, noexit=False, memchk=False, ca logger.debug(f"Return code: {cmd_return.returncode}") # # check the return status and output files - ret_ok = not cmd_return.returncode + ret_ok = (cmd_return.returncode == test['retval']) if ret_ok: out_ok = True @@ -279,9 +279,16 @@ def build_tests(test_root): logger.error("ERROR: name attribute not found for test") raise + test['retval'] = 0 for el in test_el: if (el.tag=='exec' or el.tag=='param'): test[el.tag] = repl_env(el.text) + elif el.tag=='retval': + try: + test['retval'] = int(el.text) + except ValueError: + logger.error("ERROR: retval must be an integer value") + raise elif el.tag=='output': test['out_pnc'] = [] test['out_gnc'] = [] diff --git a/internal/test_unit/xml/unit_ascii2nc.xml b/internal/test_unit/xml/unit_ascii2nc.xml index 4424fd0e33..0b2d242ef7 100644 --- a/internal/test_unit/xml/unit_ascii2nc.xml +++ b/internal/test_unit/xml/unit_ascii2nc.xml @@ -44,6 +44,18 @@ + + &MET_BIN;/ascii2nc + 1 + \ + &DATA_DIR_OBS;/gauge/2012041012.badfile.ascii \ + &OUTPUT_DIR;/ascii2nc/gauge_2012041012_24hr.nc \ + -v 1 + + + + + &MET_BIN;/ascii2nc \ diff --git a/src/basic/vx_util/num_array.h b/src/basic/vx_util/num_array.h index 6103b4b9f4..9bb7b927d0 100644 --- a/src/basic/vx_util/num_array.h +++ b/src/basic/vx_util/num_array.h @@ -57,7 +57,7 @@ class NumArray { double operator[](int) const; - const double * vals() const; + const std::vector & vals() const; double * buf(); int has(int, bool forward=true) const; @@ -127,12 +127,13 @@ class NumArray { //////////////////////////////////////////////////////////////////////// -inline int NumArray::n_elements() const { return ( e.size() ); } -inline int NumArray::n () const { return ( e.size() ); } -inline const double * NumArray::vals() const { return ( e.data() ); } -inline double * NumArray::buf() { return ( e.data() ); } -inline void NumArray::inc(int i, int v) { e[i] += v; return; } -inline void NumArray::inc(int i, double v) { e[i] += v; return; } +inline int NumArray::n_elements() const { return e.size(); } +inline int NumArray::n() const { return e.size(); } +inline const std::vector & + NumArray::vals() const { return e; } +inline double * NumArray::buf() { return e.data(); } +inline void NumArray::inc(int i, int v) { e[i] += v; return; } +inline void NumArray::inc(int i, double v) { e[i] += v; return; } //////////////////////////////////////////////////////////////////////// diff --git a/src/basic/vx_util/thresh_array.cc b/src/basic/vx_util/thresh_array.cc index 993857678f..1be3bb723b 100644 --- a/src/basic/vx_util/thresh_array.cc +++ b/src/basic/vx_util/thresh_array.cc @@ -617,6 +617,13 @@ ThreshArray define_prob_bins(double beg, double end, double inc, int prec) { v += inc; } + // Add final 1.0 threshold, if needed + v = ta[(ta.n() - 1)].get_value(); + if(v < 1.0 && !is_eq(v, 1.0)) { + cs << cs_erase << ">=1.0"; + ta.add(cs.c_str()); + } + return ta; } diff --git a/src/libcode/vx_statistics/compute_ci.cc b/src/libcode/vx_statistics/compute_ci.cc index 3edbc66a8f..8296cc5d76 100644 --- a/src/libcode/vx_statistics/compute_ci.cc +++ b/src/libcode/vx_statistics/compute_ci.cc @@ -171,9 +171,8 @@ void compute_wilson_ci(double p, int n_int, double alpha, double vif, //////////////////////////////////////////////////////////////////////// void compute_woolf_ci(double odds, double alpha, - int fy_oy, int fy_on, int fn_oy, int fn_on, + double fy_oy, double fy_on, double fn_oy, double fn_on, double &odds_cl, double &odds_cu) { - double cv_normal_l, cv_normal_u, a, b; if(is_bad_data(odds) || fy_oy == 0 || fy_on == 0 || fn_oy == 0 || fn_on == 0) { @@ -185,14 +184,14 @@ void compute_woolf_ci(double odds, double alpha, // Compute the upper and lower critical values from the // normal distribution. // - cv_normal_l = normal_cdf_inv(alpha/2.0, 0.0, 1.0); - cv_normal_u = normal_cdf_inv(1.0 - (alpha/2.0), 0.0, 1.0); + double cv_normal_l = normal_cdf_inv(alpha/2.0, 0.0, 1.0); + double cv_normal_u = normal_cdf_inv(1.0 - (alpha/2.0), 0.0, 1.0); // // Compute the upper and lower bounds of the confidence interval // - a = exp(cv_normal_l*sqrt(1.0/fy_oy + 1.0/fy_on + 1.0/fn_oy + 1.0/fn_on)); - b = exp(cv_normal_u*sqrt(1.0/fy_oy + 1.0/fy_on + 1.0/fn_oy + 1.0/fn_on)); + double a = exp(cv_normal_l*sqrt(1.0/fy_oy + 1.0/fy_on + 1.0/fn_oy + 1.0/fn_on)); + double b = exp(cv_normal_u*sqrt(1.0/fy_oy + 1.0/fy_on + 1.0/fn_oy + 1.0/fn_on)); odds_cl = odds * a; odds_cu = odds * b; @@ -210,19 +209,16 @@ void compute_woolf_ci(double odds, double alpha, //////////////////////////////////////////////////////////////////////// void compute_hk_ci(double hk, double alpha, double vif, - int fy_oy, int fy_on, int fn_oy, int fn_on, + double fy_oy, double fy_on, double fn_oy, double fn_on, double &hk_cl, double &hk_cu) { - double cv_normal, stdev; - double h, h_var, f_var; - int h_n, f_n; // // Get the counts // - h_n = fy_oy + fn_oy; - f_n = fn_on + fy_on; + double h_n = fy_oy + fn_oy; + double f_n = fn_on + fy_on; - if(is_bad_data(hk) || h_n == 0 || f_n == 0) { + if(is_bad_data(hk) || is_eq(h_n, 0.0) || is_eq(f_n, 0.0)) { hk_cl = hk_cu = bad_data_double; return; } @@ -231,26 +227,26 @@ void compute_hk_ci(double hk, double alpha, double vif, // Compute the critical value for the normal distribution based // on the sample size // - cv_normal = normal_cdf_inv(alpha/2.0, 0.0, 1.0); + double cv_normal = normal_cdf_inv(alpha/2.0, 0.0, 1.0); // // Compute the hit rate and false alarm rate // - h = (double) fy_oy/h_n; + double h = fy_oy/h_n; // // Compute a variance for H and F // - h_var = sqrt(h*(1.0-h)/h_n + cv_normal*cv_normal/(4.0*h_n*h_n)) - / (1.0 + cv_normal*cv_normal/h_n); + double h_var = sqrt(h*(1.0-h)/h_n + cv_normal*cv_normal/(4.0*h_n*h_n)) + / (1.0 + cv_normal*cv_normal/h_n); - f_var = sqrt(h*(1.0-h)/f_n + cv_normal*cv_normal/(4.0*f_n*f_n)) - / (1.0 + cv_normal*cv_normal/f_n); + double f_var = sqrt(h*(1.0-h)/f_n + cv_normal*cv_normal/(4.0*f_n*f_n)) + / (1.0 + cv_normal*cv_normal/f_n); // // Compute the standard deviation for HK // - stdev = sqrt(vif*(h_var*h_var + f_var*f_var)); + double stdev = sqrt(vif*(h_var*h_var + f_var*f_var)); // // Compute the upper and lower bounds of the confidence interval diff --git a/src/libcode/vx_statistics/compute_ci.h b/src/libcode/vx_statistics/compute_ci.h index 03dc64d8cf..ddaf68d36a 100644 --- a/src/libcode/vx_statistics/compute_ci.h +++ b/src/libcode/vx_statistics/compute_ci.h @@ -38,11 +38,11 @@ extern void compute_wilson_ci(double p, int n, double alpha, double vif, double &p_cl, double &p_cu); extern void compute_woolf_ci(double odds, double alpha, - int fy_oy, int fy_on, int fn_oy, int fn_on, + double fy_oy, double fy_on, double fn_oy, double fn_on, double &odds_cl, double &odds_cu); extern void compute_hk_ci(double hk, double alpha, double vif, - int fy_oy, int fy_on, int fn_oy, int fn_on, + double fy_oy, double fy_on, double fn_oy, double fn_on, double &hk_cl, double &hk_cu); extern void compute_cts_stats_ci_bca(const gsl_rng *, diff --git a/src/libcode/vx_statistics/contable.cc b/src/libcode/vx_statistics/contable.cc index 8c6ba55c03..a5b2ab49a1 100644 --- a/src/libcode/vx_statistics/contable.cc +++ b/src/libcode/vx_statistics/contable.cc @@ -6,11 +6,8 @@ // ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* - - //////////////////////////////////////////////////////////////////////// - #include #include #include @@ -26,1725 +23,659 @@ using namespace std; - -//////////////////////////////////////////////////////////////////////// - - -static int table_rc_to_n(int r_table, int c_table, int w, int h); - - //////////////////////////////////////////////////////////////////////// - - - // - // Code for class ContingencyTable - // - - +// +// Code for class ContingencyTable +// //////////////////////////////////////////////////////////////////////// - -ContingencyTable::ContingencyTable() -{ - -init_from_scratch(); - +ContingencyTable::ContingencyTable() { + init_from_scratch(); } - //////////////////////////////////////////////////////////////////////// - -ContingencyTable::~ContingencyTable() - -{ - -clear(); - +ContingencyTable::~ContingencyTable() { + clear(); } - //////////////////////////////////////////////////////////////////////// - -ContingencyTable::ContingencyTable(const ContingencyTable & t) -{ - -init_from_scratch(); - -assign(t); - +ContingencyTable::ContingencyTable(const ContingencyTable & t) { + init_from_scratch(); + assign(t); } - //////////////////////////////////////////////////////////////////////// - -ContingencyTable & ContingencyTable::operator=(const ContingencyTable & t) - -{ - -if ( this == &t ) return *this; - -assign(t); - -return *this; - +ContingencyTable & ContingencyTable::operator=(const ContingencyTable & t) { + if(this == &t) return *this; + assign(t); + return *this; } - //////////////////////////////////////////////////////////////////////// +ContingencyTable & ContingencyTable::operator+=(const ContingencyTable & t) { -ContingencyTable & ContingencyTable::operator+=(const ContingencyTable & t) - -{ - -if ( Nrows != t.Nrows || Ncols != t.Ncols ) { - - mlog << Error << "\nContingencyTable::operator+=() -> " - << "table dimensions do not match: (" << Nrows << ", " << Ncols - << ") != (" << t.Nrows << ", " << t.Ncols << ")\n\n"; - - exit ( 1 ); - -} - -if ( !is_eq(ECvalue, t.ECvalue) ) { - - mlog << Error << "\nContingencyTable::operator+=() -> " - << "the expected correct values do not match: " - << ECvalue << " != " << t.ECvalue << "\n\n"; - - exit ( 1 ); - -} + // Check consistent dimensions + if(Nrows != t.Nrows || Ncols != t.Ncols) { + mlog << Error << "\nContingencyTable::operator+=() -> " + << "table dimensions do not match: (" << Nrows << ", " << Ncols + << ") != (" << t.Nrows << ", " << t.Ncols << ")\n\n"; + exit(1); + } -if ( E ) { - for ( int i=0; isize(); ++i ) (*E)[i] += (*t.E)[i]; -} + // Check consistent expected correct + if(!is_eq(ECvalue, t.ECvalue)) { + mlog << Error << "\nContingencyTable::operator+=() -> " + << "the expected correct values do not match: " + << ECvalue << " != " << t.ECvalue << "\n\n"; + exit(1); + } -return *this; + // Increment table entries + for(int i=0; i(); - Name.clear(); - Nrows = Ncols = 0; +void ContingencyTable::init_from_scratch() { + // No pointers to initialize } - //////////////////////////////////////////////////////////////////////// +void ContingencyTable::clear() { -void ContingencyTable::clear() -{ - if (E) delete E; - E = new vector(); - + E.clear(); + Nrows = Ncols = 0; ECvalue = bad_data_double; Name.clear(); - Nrows = Ncols = 0; - + return; - } - //////////////////////////////////////////////////////////////////////// +void ContingencyTable::assign(const ContingencyTable & t) { -void ContingencyTable::assign(const ContingencyTable & t) -{ - - clear(); - - if(t.E->size() == 0) return; - - ContingencyTable::set_size(t.Nrows, t.Ncols); - - if (E) delete E; - E = new vector(*(t.E)); - ECvalue = t.ECvalue; - Name = t.Name; - - // - // done - // - - return; - -} - - -//////////////////////////////////////////////////////////////////////// + clear(); + E = t.E; + Nrows = t.Nrows; + Ncols = t.Ncols; + ECvalue = t.ECvalue; + Name = t.Name; -void ContingencyTable::zero_out() -{ - - int n = Nrows*Ncols; - - if ( n == 0 ) return; - - E->assign(n, 0); - - return; - + return; } //////////////////////////////////////////////////////////////////////// +void ContingencyTable::zero_out() { -void ContingencyTable::dump(ostream & out, int depth) const - -{ + fill(E.begin(), E.end(), 0.0); -int r, c; -Indent prefix(depth); -ConcatString junk; - -out << prefix << "Name = "; + return; +} -if ( Name.length() > 0 ) out << '\"' << Name << "\"\n"; -else out << "(nul)\n"; +//////////////////////////////////////////////////////////////////////// -out << prefix << "Nrows = " << Nrows << "\n"; -out << prefix << "Ncols = " << Ncols << "\n"; -out << prefix << "ECvalue = " << ECvalue << "\n"; -out << prefix << "\n"; +void ContingencyTable::dump(ostream & out, int depth) const { + Indent prefix(depth); + ConcatString msg; -if ( E->empty() ) { out.flush(); return; } + out << prefix << "Name = "; -for (r=0; r col_width(Ncols); + ContingencyTable::set_size(N, N); -for (c=0; c " + << "# rows (" << NR << ") and # cols (" << NC + << ") must be at least 2!\n\n"; + exit(1); + } - if ( k > col_width[c] ) col_width[c] = k; + Nrows = NR; + Ncols = NC; - for (r=0; r col_width[c] ) col_width[c] = k; +void ContingencyTable::set_ec_value(double v) { - } + // Do not override the default value with bad data + if(!is_bad_data(v)) ECvalue = v; + return; } -w = 2*hpad*Ncols + Ncols + 1; - -for (c=0; c table(w*h, ' '); + Name = text; - // - // top, bottom - // + return; +} -for (c_table=0; c_table= Nrows || c < 0 || c >= Ncols) { + mlog << Error << "\nContingencyTable::rc_to_n() -> " + << "range check error requesting (" << r << ", " + << c << ") from table with dimension (" << Nrows + << ", " << Ncols << ")!\n\n"; + exit(1); + } - table[n] = '='; + return r*Ncols + c; +} - n = table_rc_to_n(h - 1, c_table, w, h); +//////////////////////////////////////////////////////////////////////// +void ContingencyTable::set_entry(int row, int col, double value) { - table[n] = '='; + E[(rc_to_n(row, col))] = value; + return; } - // - // left, right - // - -for (r_table=1; r_table<(h - 1); ++r_table) { - - n = table_rc_to_n(r_table, 0, w, h); - - table[n] = v_sep; +//////////////////////////////////////////////////////////////////////// - n = table_rc_to_n(r_table, w - 1, w, h); +void ContingencyTable::inc_entry(int row, int col, double weight) { - table[n] = v_sep; + E[(rc_to_n(row, col))] += weight; + return; } - // - // col separators - // - -for (c=1; c " - << "c_table (" << c_table << ") is greater then w (" << w << ")\n\n"; + for(int col=0; col " - << "# rows (" << NR << ") and # cols (" << NC - << ") must be at least 2!\n\n"; - - exit ( 1 ); - +TTContingencyTable::TTContingencyTable(const TTContingencyTable & t) { + assign(t); } -int n; - -n = NR*NC; - -E->resize(n, 0); - -Nrows = NR; -Ncols = NC; - - // - // if square, set default expected correct value - // - -if ( Nrows == Ncols ) { - - ECvalue = 1.0 / Nrows; +//////////////////////////////////////////////////////////////////////// -} +TTContingencyTable & TTContingencyTable::operator=(const TTContingencyTable & t) { - // - // done - // + if(this == &t) return *this; -return; + assign(t); + return *this; } - //////////////////////////////////////////////////////////////////////// +void TTContingencyTable::set_fn_on(double k) { -void ContingencyTable::set_ec_value(double v) - -{ - - // - // do not override the default value with bad data - // - -if ( !is_bad_data(v) ) ECvalue = v; - -return; + set_entry(FN_row, ON_col, k); + return; } - //////////////////////////////////////////////////////////////////////// +void TTContingencyTable::set_fy_on(double k) { -void ContingencyTable::set_name(const char * text) - -{ - -Name = text; - -return; + set_entry(FY_row, ON_col, k); + return; } - //////////////////////////////////////////////////////////////////////// +void TTContingencyTable::set_fn_oy(double k) { -int ContingencyTable::rc_to_n(int r, int c) const - -{ - -if ( (r < 0) || (r >= Nrows) || (c < 0) || (c >= Ncols) ) { - - mlog << Error << "\nContingencyTable::rc_to_n() -> " - << "range check error!\n\n"; - - exit ( 1 ); + set_entry(FN_row, OY_col, k); + return; } -int n; - - -n = r*Ncols + c; - +//////////////////////////////////////////////////////////////////////// +void TTContingencyTable::set_fy_oy(double k) { -return n; + set_entry(FY_row, OY_col, k); + return; } - //////////////////////////////////////////////////////////////////////// +void TTContingencyTable::inc_fn_on(double weight) { -void ContingencyTable::set_entry(int row, int col, int value) - -{ - -int n; - -n = rc_to_n(row, col); - -(*E)[n] = value; - - -return; + inc_entry(FN_row, ON_col, weight); + return; } - //////////////////////////////////////////////////////////////////////// +void TTContingencyTable::inc_fy_on(double weight) { -void ContingencyTable::inc_entry(int row, int col) - -{ - -int n; - -n = rc_to_n(row, col); - -++((*E)[n]); - - - -return; + inc_entry(FY_row, ON_col, weight); + return; } - //////////////////////////////////////////////////////////////////////// +void TTContingencyTable::inc_fn_oy(double weight) { -int ContingencyTable::total() const - -{ - -const int n = Nrows*Ncols; - -if ( n == 0 ) return 0; + inc_entry(FN_row, OY_col, weight); -int j, sum; + return; +} -sum = 0; +//////////////////////////////////////////////////////////////////////// -for (j=0; j= Nrows) ) { +//////////////////////////////////////////////////////////////////////// - mlog << Error << "\nContingencyTable::row_total() -> " - << "range check error!\n\n"; +double TTContingencyTable::fn_on() const { + return entry(FN_row, ON_col); +} - exit ( 1 ); +//////////////////////////////////////////////////////////////////////// +double TTContingencyTable::fy() const { + return row_total(FY_row); } -int n, col, sum; - +//////////////////////////////////////////////////////////////////////// -sum = 0; +double TTContingencyTable::fn() const { + return row_total(FN_row); +} -for (col=0; col= Ncols) ) { +//////////////////////////////////////////////////////////////////////// - mlog << Error << "\nContingencyTable::col_total() -> " - << "range check error!\n\n"; +double TTContingencyTable::fy_oy_tp() const { + return compute_proportion(fy_oy(), n()); +} - exit ( 1 ); +//////////////////////////////////////////////////////////////////////// +double TTContingencyTable::fy_on_tp() const { + return compute_proportion(fy_on(), n()); } -int n, row, sum; - +//////////////////////////////////////////////////////////////////////// -sum = 0; +double TTContingencyTable::fn_oy_tp() const { + return compute_proportion(fn_oy(), n()); +} -for (row=0; row a ) a = (*E)[n]; +//////////////////////////////////////////////////////////////////////// +double TTContingencyTable::fn_on_op() const { + return compute_proportion(fn_on(), on()); } +//////////////////////////////////////////////////////////////////////// -return a; - +void TTContingencyTable::set_size(int N) { + mlog << Error << "\nTTContingencyTable::set_size(int) -> " + << "2 x 2 tables cannot be resized!\n\n"; + exit(1); } - //////////////////////////////////////////////////////////////////////// +void TTContingencyTable::set_size(int NR, int NC) { + mlog << Error << "\nTTContingencyTable::set_size(int, int) -> " + << "2 x 2 tables cannot be resized!\n\n"; + exit(1); +} -int ContingencyTable::smallest_entry() const +//////////////////////////////////////////////////////////////////////// +// +// Code for misc functions +// +//////////////////////////////////////////////////////////////////////// -{ +// +// Reference table 7.1a, page 242 in wilks +// -int n = Nrows*Ncols; +TTContingencyTable finley() { + TTContingencyTable t; -if ( n == 0 ) return 0; + t.set_fy_oy(28); + t.set_fn_oy(23); + t.set_fy_on(72); + t.set_fn_on(2680); -int j, a; + t.set_name("Finley Tornado Forecasts (1884)"); -a = (*E)[0]; + return t; +} -for (j=1; j " - << "table not square!\n\n"; - - exit ( 1 ); - -} - -if ( (k < 0) || (k >= Nrows) ) { - - mlog << Error << "\nContingencyTable::condition_on() -> " - << "range check error\n\n"; - - exit ( 1 ); - -} - -int r, c; -int n, sum; -TTContingencyTable t; - - // - // - // - -t.set_entry(0, 0, entry(k, k)); - - // - // - // - -sum = 0; - -for (c=0; c " - << "2 x 2 tables cannot be resized!\n\n"; - -exit ( 1 ); - -} - - -//////////////////////////////////////////////////////////////////////// - - -void TTContingencyTable::set_size(int NR, int NC) - -{ - -mlog << Error << "\nTTContingencyTable::set_size(int, int) -> " - << "2 x 2 tables cannot be resized!\n\n"; - -exit ( 1 ); - -} - - -//////////////////////////////////////////////////////////////////////// - - - // - // Code for misc functions - // - - -//////////////////////////////////////////////////////////////////////// - - // - // see table 7.1a, page 242 in wilks - // - -TTContingencyTable finley() - -{ - -TTContingencyTable t; - - -t.set_fy_oy(28); -t.set_fn_oy(23); - -t.set_fy_on(72); -t.set_fn_on(2680); - - -t.set_name("Finley Tornado Forecasts (1884)"); - - -return t; - -} - - -//////////////////////////////////////////////////////////////////////// - - // - // see table 7.1b, page 242 in wilks - // - -TTContingencyTable finley_always_no() - -{ - -TTContingencyTable t; - - -t.set_fy_oy(0); -t.set_fn_oy(51); - -t.set_fy_on(0); -t.set_fn_on(2752); - - -t.set_name("Finley Tornado Forecasts (Always No) (1884)"); - - -return t; - -} - - -//////////////////////////////////////////////////////////////////////// -// r_table < h -// c_table < w - -int table_rc_to_n(int r_table, int c_table, int w, int h) - -{ - -int n; - -n = r_table*w + c_table; - -return n; + // Check for bad data and divide by zero + if(is_bad_data(num) || + is_bad_data(den) || + is_eq(den, 0.0)) { + prop = bad_data_double; + } + else { + prop = num/den; + } + return prop; } - //////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_statistics/contable.h b/src/libcode/vx_statistics/contable.h index 8c7823c4d7..572893d095 100644 --- a/src/libcode/vx_statistics/contable.h +++ b/src/libcode/vx_statistics/contable.h @@ -6,37 +6,28 @@ // ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* - - //////////////////////////////////////////////////////////////////////// - #ifndef __CONTINGENCY_TABLE_H__ #define __CONTINGENCY_TABLE_H__ - //////////////////////////////////////////////////////////////////////// - #include #include "vx_util.h" #include "vx_math.h" - //////////////////////////////////////////////////////////////////////// - -class TTContingencyTable; // forward reference - +// Forward reference +class TTContingencyTable; //////////////////////////////////////////////////////////////////////// - - - // - // general contingency table - // - +// +// General contingency table +// +//////////////////////////////////////////////////////////////////////// class ContingencyTable { @@ -48,7 +39,8 @@ class ContingencyTable { int rc_to_n(int r, int c) const; - std::vector *E; // this is really a two-dimensional array + // This is really a two-dimensional array (Nrows, Ncols) + std::vector E; int Nrows; int Ncols; @@ -71,101 +63,73 @@ class ContingencyTable { virtual void dump(std::ostream & out, int depth = 0) const; - // - // condition on an event - // - - TTContingencyTable condition_on(int) const; - - // - // set attributes - // - + // Set attributes virtual void set_size(int); virtual void set_size(int NR, int NC); void set_ec_value(double); void set_name(const char *); - // - // get attributes - // - + // Get attributes int nrows() const; int ncols() const; double ec_value() const; ConcatString name() const; - // - // set counts - // + // Set table entries + void set_entry(int row, int col, double value); - void set_entry(int row, int col, int value); + // Increment table entries + void inc_entry(int row, int col, double weight=1.0); - // - // increment counts - // + // Get values + double total() const; - void inc_entry(int row, int col); + double row_total(int row) const; + double col_total(int col) const; - // - // get counts - // + double entry(int row, int col) const; - int total() const; + double max() const; + double min() const; - int row_total(int row) const; - int col_total(int col) const; - - int entry(int row, int col) const; - - int largest_entry() const; - int smallest_entry() const; - - // - // statistics - // + bool is_integer() const; + // Statistics virtual double gaccuracy () const; virtual double gheidke () const; virtual double gheidke_ec(double) const; virtual double gkuiper () const; virtual double gerrity () const; - }; //////////////////////////////////////////////////////////////////////// - inline int ContingencyTable::nrows() const { return Nrows; } inline int ContingencyTable::ncols() const { return Ncols; } inline double ContingencyTable::ec_value() const { return ECvalue; } inline ConcatString ContingencyTable::name() const { return Name; } - //////////////////////////////////////////////////////////////////////// - static const int nx2_event_column = 0; static const int nx2_nonevent_column = 1; - //////////////////////////////////////////////////////////////////////// - - - // - // N x 2 contingency table - // - +// +// N x 2 contingency table +// +//////////////////////////////////////////////////////////////////////// class Nx2ContingencyTable : public ContingencyTable { private: - double * Thresholds; // N + 1 count, increasing + // N + 1 count, parametrically increasing or decreasing + std::vector Thresholds; int value_to_row(double) const; @@ -182,54 +146,36 @@ class Nx2ContingencyTable : public ContingencyTable { void clear(); - void set_size(int NR); - void set_size(int NR, int NC); // NC had better be 2 - - void set_thresholds(const double *); - - // - // get thresholds - // - - double threshold(int index) const; // 0 <= index <= Nrows - - // - // set counts - // + void set_size(int NR) override; + void set_size(int NR, int NC) override; // NC must be 2 - void set_event(int row, int count); - void set_nonevent(int row, int count); + void set_thresholds(const std::vector &); - // - // increment counts - // + // Get thresholds + double threshold(int index) const; // 0 <= index <= Nrows - void inc_event (double); - void inc_nonevent (double); + // Increment table entries + void inc_event (double value, double weight=1.0); + void inc_nonevent (double value, double weight=1.0); - // - // get counts - // + // Get table entries + double event_count_by_thresh(double) const; + double nonevent_count_by_thresh(double) const; - int event_count_by_thresh(double) const; - int nonevent_count_by_thresh(double) const; + double event_count_by_row(int row) const; + double nonevent_count_by_row(int row) const; - int event_count_by_row(int row) const; - int nonevent_count_by_row(int row) const; + // Set counts + void set_event(int row, double); + void set_nonevent(int row, double); - int n() const; + double n() const; - // - // column totals - // - - int event_col_total() const; - int nonevent_col_total() const; - - // - // statistics - // + // Column totals + double event_col_total() const; + double nonevent_col_total() const; + // Statistics double baser () const; double baser_ci (double alpha, double &cl, double &cu) const; double brier_score () const; @@ -252,27 +198,21 @@ class Nx2ContingencyTable : public ContingencyTable { TTContingencyTable ctc_by_row (int row) const; double roc_auc() const; - }; - //////////////////////////////////////////////////////////////////////// +inline double Nx2ContingencyTable::event_count_by_row (int row) const { return entry(row, nx2_event_column); } +inline double Nx2ContingencyTable::nonevent_count_by_row (int row) const { return entry(row, nx2_nonevent_column); } -inline int Nx2ContingencyTable::event_col_total () const { return col_total(nx2_event_column); } -inline int Nx2ContingencyTable::nonevent_col_total () const { return col_total(nx2_nonevent_column); } - -inline int Nx2ContingencyTable::event_count_by_row (int row) const { return entry(row, nx2_event_column); } -inline int Nx2ContingencyTable::nonevent_count_by_row (int row) const { return entry(row, nx2_nonevent_column); } - +inline double Nx2ContingencyTable::event_col_total () const { return col_total(nx2_event_column); } +inline double Nx2ContingencyTable::nonevent_col_total () const { return col_total(nx2_nonevent_column); } //////////////////////////////////////////////////////////////////////// - - - // - // 2 x 2 contingency table - // - +// +// 2 x 2 contingency table +// +//////////////////////////////////////////////////////////////////////// class TTContingencyTable : public ContingencyTable { @@ -283,63 +223,47 @@ class TTContingencyTable : public ContingencyTable { TTContingencyTable(const TTContingencyTable &); TTContingencyTable & operator=(const TTContingencyTable &); - void set_size(int); - void set_size(int NR, int NC); - - // - // set counts - // - - void set_fn_on(int); - void set_fy_on(int); + void set_size(int) override; + void set_size(int NR, int NC) override; - void set_fn_oy(int); - void set_fy_oy(int); + // Set table entries + void set_fn_on(double); + void set_fy_on(double); - // - // increment counts - // + void set_fn_oy(double); + void set_fy_oy(double); - void inc_fn_on(); - void inc_fy_on(); + // Increment table entries + void inc_fn_on(double weight=1.0); + void inc_fy_on(double weight=1.0); - void inc_fn_oy(); - void inc_fy_oy(); + void inc_fn_oy(double weight=1.0); + void inc_fy_oy(double weight=1.0); - // - // get counts - // + // Get table entries + double fn_on() const; + double fy_on() const; - int fn_on() const; - int fy_on() const; + double fn_oy() const; + double fy_oy() const; - int fn_oy() const; - int fy_oy() const; + double on() const; + double oy() const; - int on() const; - int oy() const; + double fn() const; + double fy() const; - int fn() const; - int fy() const; + double n() const; - int n() const; - - // - // FHO rates where: - // f_rate = FY/N - // h_rate = fy_oy/N - // o_rate = OY/N - // - - double f_rate () const; - double h_rate () const; - double o_rate () const; - - // - // Raw counts as proportions of the - // total count. - // + // FHO rates where: + // f_rate = FY/N + // h_rate = fy_oy/N + // o_rate = OY/N + double f_rate() const; + double h_rate() const; + double o_rate() const; + // Proportions of the total double fy_oy_tp () const; double fy_on_tp () const; double fn_oy_tp () const; @@ -351,31 +275,19 @@ class TTContingencyTable : public ContingencyTable { double oy_tp () const; double on_tp () const; - // - // Raw counts as proportions of the - // total forecast yes count. - // - + // Proportions of forecast double fy_oy_fp () const; double fy_on_fp () const; double fn_oy_fp () const; double fn_on_fp () const; - // - // Raw counts as proportions of the - // total observation yes count. - // - + // Proportions of observation double fy_oy_op () const; double fy_on_op () const; double fn_oy_op () const; double fn_on_op () const; - // - // Contingency Table Statistics with confidence intervals - // when applicable. - // - + // Contingency Table Statistics and confidence intervals double baser () const; double baser_ci (double alpha, double &cl, double &cu) const; double fmean () const; @@ -416,26 +328,19 @@ class TTContingencyTable : public ContingencyTable { double cost_loss (double) const; }; - //////////////////////////////////////////////////////////////////////// - extern TTContingencyTable finley(); extern TTContingencyTable finley_always_no(); - //////////////////////////////////////////////////////////////////////// - - - // - // this is the layout on page 239 of - // - // "Statistical Methods in the Atmospheric Sciences" (1st ed) - // - // by Daniel S. Wilks - // - +// +// Reference page 239 of +// "Statistical Methods in the Atmospheric Sciences" (1st ed) +// by Daniel S. Wilks +// +//////////////////////////////////////////////////////////////////////// static const int OY_col = 0; static const int ON_col = 1; @@ -443,19 +348,15 @@ static const int ON_col = 1; static const int FY_row = 0; static const int FN_row = 1; - //////////////////////////////////////////////////////////////////////// +extern void calc_gerrity_scoring_matrix(int N, const std::vector &p, + std::vector &s); -extern void calc_gerrity_scoring_matrix(int N, const double * p, double * s); - +extern double compute_proportion(double, double); //////////////////////////////////////////////////////////////////////// - #endif // __CONTINGENCY_TABLE_H__ - //////////////////////////////////////////////////////////////////////// - - diff --git a/src/libcode/vx_statistics/contable_nx2.cc b/src/libcode/vx_statistics/contable_nx2.cc index c69ec38545..15d6a18b67 100644 --- a/src/libcode/vx_statistics/contable_nx2.cc +++ b/src/libcode/vx_statistics/contable_nx2.cc @@ -6,11 +6,8 @@ // ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* - - //////////////////////////////////////////////////////////////////////// - #include #include #include @@ -26,1012 +23,582 @@ using namespace std; - //////////////////////////////////////////////////////////////////////// - static const int use_center = 1; - -//////////////////////////////////////////////////////////////////////// - - - // - // Code for class Nx2ContingencyTable - // - - -//////////////////////////////////////////////////////////////////////// - - -Nx2ContingencyTable::Nx2ContingencyTable() - -{ - -init_from_scratch(); - -} - - //////////////////////////////////////////////////////////////////////// - - -Nx2ContingencyTable::~Nx2ContingencyTable() - -{ - -clear(); - -} - - +// +// Code for class Nx2ContingencyTable +// //////////////////////////////////////////////////////////////////////// - -Nx2ContingencyTable::Nx2ContingencyTable(const Nx2ContingencyTable & t) - -{ - -init_from_scratch(); - -assign(t); - +Nx2ContingencyTable::Nx2ContingencyTable() { + init_from_scratch(); } - - //////////////////////////////////////////////////////////////////////// - -void Nx2ContingencyTable::init_from_scratch() - -{ - -ContingencyTable::init_from_scratch(); - -Thresholds = (double *) nullptr; - -clear(); - -return; - +Nx2ContingencyTable::~Nx2ContingencyTable() { + clear(); } - //////////////////////////////////////////////////////////////////////// - -void Nx2ContingencyTable::clear() - -{ - -ContingencyTable::clear(); - -if ( Thresholds ) { delete [] Thresholds; Thresholds = (double *) nullptr; } - -return; - +Nx2ContingencyTable::Nx2ContingencyTable(const Nx2ContingencyTable & t) { + init_from_scratch(); + assign(t); } - //////////////////////////////////////////////////////////////////////// - -Nx2ContingencyTable & Nx2ContingencyTable::operator=(const Nx2ContingencyTable & t) - -{ - -if ( this == &t ) return *this; - -assign(t); - -return *this; - +void Nx2ContingencyTable::init_from_scratch() { + ContingencyTable::init_from_scratch(); } - //////////////////////////////////////////////////////////////////////// +void Nx2ContingencyTable::clear() { -void Nx2ContingencyTable::assign(const Nx2ContingencyTable & t) - -{ - -clear(); - -ContingencyTable::assign(t); - -if(t.Thresholds) set_thresholds(t.Thresholds); - -return; + ContingencyTable::clear(); + Thresholds.clear(); + return; } - //////////////////////////////////////////////////////////////////////// +Nx2ContingencyTable & Nx2ContingencyTable::operator=(const Nx2ContingencyTable &t) { -int Nx2ContingencyTable::n() const - -{ - -int k; - -k = total(); - -return k; + if(this == &t) return *this; + assign(t); + return *this; } - //////////////////////////////////////////////////////////////////////// +void Nx2ContingencyTable::assign(const Nx2ContingencyTable & t) { -void Nx2ContingencyTable::set_size(int N) - -{ - -ContingencyTable::set_size(N, 2); + clear(); -return; + ContingencyTable::assign(t); + Thresholds = t.Thresholds; + return; } - //////////////////////////////////////////////////////////////////////// - -void Nx2ContingencyTable::set_size(int NR, int NC) - -{ - -if ( NC != 2 ) { - - mlog << Error << "\nNx2ContingencyTable::set_size(int, int) -> " - << "must have 2 columns!\n\n"; - - exit ( 1 ); - -} - -set_size(NR); - -return; - +double Nx2ContingencyTable::n() const { + return total(); } - //////////////////////////////////////////////////////////////////////// - -int Nx2ContingencyTable::value_to_row(double t) const - -{ - -if ( !Thresholds ) { - - mlog << Error << "\nNx2ContingencyTable::value_to_row(double) const -> " - << "thresholds array not set!\n\n"; - - exit ( 1 ); - -} - -if ( t < Thresholds[0] && !is_eq(t, Thresholds[0]) ) return -1; - -if ( t > Thresholds[Nrows] && !is_eq(t, Thresholds[Nrows]) ) return -1; - // Thresholds array is of size Nrows + 1, so - // the last element has index Nrows, not Nrows - 1 - -int j; - -for (j=0; j Thresholds[j ] || is_eq(t, Thresholds[j ]) ) && - ( t < Thresholds[j + 1] && !is_eq(t, Thresholds[j + 1]) ) ) return j; - -} - -if ( is_eq(t, Thresholds[Nrows]) ) return ( Nrows - 1 ); - -return -1; - +void Nx2ContingencyTable::set_size(int N) { + ContingencyTable::set_size(N, 2); + return; } - //////////////////////////////////////////////////////////////////////// +void Nx2ContingencyTable::set_size(int NR, int NC) { -void Nx2ContingencyTable::set_thresholds(const double * Values) - -{ - -if ( E->empty() ) { - - mlog << Error << "\nNx2ContingencyTable::set_thresholds(const double *) -> " - << "table empty!\n\n"; - - exit ( 1 ); - -} - -if ( Thresholds ) { delete [] Thresholds; Thresholds = (double *) nullptr; } - -Thresholds = new double [Nrows + 1]; - -memcpy(Thresholds, Values, (Nrows + 1)*sizeof(double)); + if(NC != 2) { + mlog << Error << "\nNx2ContingencyTable::set_size(int, int) -> " + << "must have 2 columns, not " << NC << "!\n\n"; + exit(1); + } -return; + set_size(NR); + return; } - //////////////////////////////////////////////////////////////////////// +int Nx2ContingencyTable::value_to_row(double t) const { -double Nx2ContingencyTable::threshold(int k) const - -{ - -if ( !Thresholds ) { - - mlog << Error << "\nNx2ContingencyTable::threshold(int) const -> " - << "no thresholds set!\n\n"; - - exit ( 1 ); - -} - -if ( (k < 0) || (k > Nrows) ) { // there are Nrows + 1 thresholds - - mlog << Error << "\nNx2ContingencyTable::threshold(int) const -> " - << "range check error\n\n"; + if(Thresholds.empty()) { + mlog << Error << "\nNx2ContingencyTable::value_to_row(double) const -> " + << "thresholds array not set!\n\n"; + exit(1); + } - exit ( 1 ); + // Thresholds array is of size Nrows + 1, so + // the last element has index Nrows, not Nrows - 1 + if(t < Thresholds[0] && !is_eq(t, Thresholds[0])) return -1; + if(t > Thresholds[Nrows] && !is_eq(t, Thresholds[Nrows]) ) return -1; -} + for(int j=0; j Thresholds[j] || is_eq(t, Thresholds[j]) ) && + (t < Thresholds[j + 1] && !is_eq(t, Thresholds[j + 1]))) return j; + } -return Thresholds[k]; + if(is_eq(t, Thresholds[Nrows])) return (Nrows - 1); + return -1; } //////////////////////////////////////////////////////////////////////// +void Nx2ContingencyTable::set_thresholds(const vector &Values) { -void Nx2ContingencyTable::set_event(int row, int value) - -{ - -if ( row < 0 || row >= Nrows ) { - - mlog << Error << "\nNx2ContingencyTable::set_event(double) -> " - << "bad row index ... " << row << "\n\n"; - - exit ( 1 ); - -} - -set_entry(row, nx2_event_column, value); + if(Values.size() != Nrows + 1) { + mlog << Error << "\nNx2ContingencyTable::set_thresholds(const double *) -> " + << "expected " << Nrows + 1 << " thresholds but only received " + << Values.size() << "!\n\n"; + exit(1); + } -return; + Thresholds = Values; + return; } - //////////////////////////////////////////////////////////////////////// +double Nx2ContingencyTable::threshold(int k) const { -void Nx2ContingencyTable::set_nonevent(int row, int value) - -{ - -if ( row < 0 || row >= Nrows ) { - - mlog << Error << "\nNx2ContingencyTable::set_nonevent(double) -> " - << "bad row index ... " << row << "\n\n"; - - exit ( 1 ); - -} - -set_entry(row, nx2_nonevent_column, value); - -return; + // There are Nrows + 1 thresholds + if(k < 0 || k > Thresholds.size()) { + mlog << Error << "\nNx2ContingencyTable::threshold(int) const -> " + << "range check error!\n\n"; + exit(1); + } + return Thresholds[k]; } - //////////////////////////////////////////////////////////////////////// +void Nx2ContingencyTable::inc_event(double t, double weight) { + int r = value_to_row(t); -void Nx2ContingencyTable::inc_event(double t) - -{ - -int r; - -r = value_to_row(t); - -if ( r < 0 ) { - - mlog << Error << "\nNx2ContingencyTable::inc_event(double) -> " - << "bad value ... " << t << "\n\n"; - - exit ( 1 ); + if(r < 0) { + mlog << Error << "\nNx2ContingencyTable::inc_event(double) -> " + << "bad value ... " << t << "\n\n"; + exit(1); + } -} - -inc_entry(r, nx2_event_column); - -return; + inc_entry(r, nx2_event_column, weight); + return; } - //////////////////////////////////////////////////////////////////////// +void Nx2ContingencyTable::inc_nonevent(double t, double weight) { + int r = value_to_row(t); -void Nx2ContingencyTable::inc_nonevent(double t) - -{ - -int r; - -r = value_to_row(t); - -if ( r < 0 ) { - - mlog << Error << "\nNx2ContingencyTable::inc_nonevent(double) -> " - << "bad value ... " << t << "\n\n"; - - exit ( 1 ); - -} - -inc_entry(r, nx2_nonevent_column); + if(r < 0) { + mlog << Error << "\nNx2ContingencyTable::inc_nonevent(double) -> " + << "bad value ... " << t << "\n\n"; + exit(1); + } -return; + inc_entry(r, nx2_nonevent_column, weight); + return; } - //////////////////////////////////////////////////////////////////////// +double Nx2ContingencyTable::event_count_by_thresh(double t) const { + int r = value_to_row(t); -int Nx2ContingencyTable::event_count_by_thresh(double t) const - -{ - -int r; - -r = value_to_row(t); - -if ( r < 0 ) { - - mlog << Error << "\nNx2ContingencyTable::event_count_by_thresh(double) -> " - << "bad value ... " << t << "\n\n"; - - exit ( 1 ); - -} - -int k; - -k = entry(r, nx2_event_column); - -return k; + if(r < 0) { + mlog << Error << "\nNx2ContingencyTable::event_count_by_thresh(double) -> " + << "bad value ... " << t << "\n\n"; + exit(1); + } + return entry(r, nx2_event_column); } - //////////////////////////////////////////////////////////////////////// +double Nx2ContingencyTable::nonevent_count_by_thresh(double t) const { + int r = value_to_row(t); -int Nx2ContingencyTable::nonevent_count_by_thresh(double t) const - -{ - -int r; - -r = value_to_row(t); - -if ( r < 0 ) { - - mlog << Error << "\nNx2ContingencyTable::nonevent_count_by_thresh(double) -> " - << "bad value ... " << t << "\n\n"; - - exit ( 1 ); - -} - -int k; - -k = entry(r, nx2_nonevent_column); - -return k; + if(r < 0) { + mlog << Error << "\nNx2ContingencyTable::nonevent_count_by_thresh(double) -> " + << "bad value ... " << t << "\n\n"; + exit(1); + } + return entry(r, nx2_nonevent_column); } - //////////////////////////////////////////////////////////////////////// +void Nx2ContingencyTable::set_event(int row, double value) { -double Nx2ContingencyTable::row_obar(int row) const - -{ + if(row < 0 || row >= Nrows) { + mlog << Error << "\nNx2ContingencyTable::set_event(int, double) -> " + << "bad row index ... " << row << "\n\n"; + exit(1); + } -const int obs_count = event_count_by_row(row); -const int Ni = row_total(row); -double x; - -if(Ni == 0) x = bad_data_double; -else x = ((double) obs_count)/((double) Ni); - -return x; + set_entry(row, nx2_event_column, value); + return; } - //////////////////////////////////////////////////////////////////////// +void Nx2ContingencyTable::set_nonevent(int row, double value) { -double Nx2ContingencyTable::obar() const + if(row < 0 || row >= Nrows) { + mlog << Error << "\nNx2ContingencyTable::set_nonevent(int, double) -> " + << "bad row index ... " << row << "\n\n"; + exit(1); + } -{ - -const int obs_count = event_col_total(); -const int N = n(); -double x; - -if (N == 0) x = bad_data_double; -else x = ((double) obs_count)/((double) N); - -return x; + set_entry(row, nx2_nonevent_column, value); + return; } - //////////////////////////////////////////////////////////////////////// - -double Nx2ContingencyTable::row_proby(int row) const - -{ - -if ( (row < 0) || (row >= Nrows) ) { - - mlog << Error << "\nNx2ContingencyTable::row_proby(int) const -> " - << "range check error\n\n"; - - exit ( 1 ); - -} - -double x; - -if ( use_center ) x = 0.5*(Thresholds[row] + Thresholds[row + 1]); -else x = Thresholds[row]; - -return x; - -} - - -//////////////////////////////////////////////////////////////////////// - - double Nx2ContingencyTable::baser() const { - - double v; - - if( n() == 0 ) v = bad_data_double; - else v = (double) event_col_total()/n(); - - return v; + return compute_proportion(event_col_total(), n()); } - //////////////////////////////////////////////////////////////////////// - double Nx2ContingencyTable::baser_ci(double alpha, double &cl, double &cu) const { - double v; - - v = baser(); + double v = baser(); compute_proportion_ci(v, n(), alpha, 1.0, cl, cu); return v; } - //////////////////////////////////////////////////////////////////////// +double Nx2ContingencyTable::brier_score() const { - // - // Reference: Equation 7.40, page 286 in Wilks, 2nd Ed. - // + if(E.empty()) return bad_data_double; -double Nx2ContingencyTable::reliability() const + double sum = 0.0; + double count; + double yi; + double t; -{ + // Terms for event + for(int j=0; j 1 so that degf > 0 in the call to gsl_cdf_tdist_Pinv() + if(is_bad_data(bs) || N <= 1.0) return bad_data_double; - obari = row_obar(row); + double degf = N - 1.0; + double t = gsl_cdf_tdist_Pinv(1.0 - 0.5*alpha, degf); + double ob = obar(); - // When obari is not defined, don't include it in the sum - if(is_bad_data(obari)) continue; + double af1 = 0.0; + double sf2 = 0.0; + double sf3 = 0.0; + double af4 = 0.0; - t = obari - Obar; + for(int j=0; j= Nrows) ) { - - mlog << Error << "\nNx2ContingencyTable::row_calibration(int) const -> " - << "range check error\n\n"; - - exit ( 1 ); - -} - -double num, denom; -double x; - -num = (double) event_count_by_row(row); - -denom = num + nonevent_count_by_row(row); - -if(is_eq(denom, 0.0)) x = bad_data_double; -else x = num/denom; +double Nx2ContingencyTable::uncertainty() const { + double a = obar(); + double v; -return x; + if(is_bad_data(a)) v = bad_data_double; + else v = a*(1.0 - a); + return v; } - +//////////////////////////////////////////////////////////////////////// +// +// Reference: Equation 8.43, page 340 in Wilks, 3rd Ed. +// //////////////////////////////////////////////////////////////////////// +double Nx2ContingencyTable::bss_smpl() const { + double res = resolution(); + double rel = reliability(); + double unc = uncertainty(); + double bss; -double Nx2ContingencyTable::row_refinement(int row) const - -{ - -if ( (row < 0) || (row >= Nrows) ) { - - mlog << Error << "\nNx2ContingencyTable::row_refinement(int) const -> " - << "range check error\n\n"; - - exit ( 1 ); - -} - -int py_o1, py_o2; -const int N = n(); -double x; - -py_o1 = event_count_by_row(row); -py_o2 = nonevent_count_by_row(row); - -x = (double) (py_o1 + py_o2); - -if (N == 0) x = bad_data_double; -else x /= N; - -return x; + if(is_bad_data(res) || is_bad_data(rel) || + is_bad_data(unc) || is_eq(unc, 0.0)) { + bss = bad_data_double; + } + else { + bss = (res - rel)/unc; + } + return bss; } - //////////////////////////////////////////////////////////////////////// - -double Nx2ContingencyTable::row_event_likelihood(int row) const - -{ - -if ( (row < 0) || (row >= Nrows) ) { - - mlog << Error << "\nNx2ContingencyTable::row_event_likelihood(int) const -> " - << "range check error\n\n"; - - exit ( 1 ); - +double Nx2ContingencyTable::row_obar(int row) const { + return compute_proportion(event_count_by_row(row), row_total(row)); } -double x, num, denom; - -denom = (double) event_col_total(); - -num = (double) event_count_by_row(row); - -if(is_eq(denom, 0.0)) x = bad_data_double; -else x = num/denom; - -return x; +//////////////////////////////////////////////////////////////////////// +double Nx2ContingencyTable::obar() const { + return compute_proportion(event_col_total(), n()); } - //////////////////////////////////////////////////////////////////////// +double Nx2ContingencyTable::row_proby(int row) const { -double Nx2ContingencyTable::row_nonevent_likelihood(int row) const - -{ - -if ( (row < 0) || (row >= Nrows) ) { + if(row < 0 || row >= Nrows) { + mlog << Error << "\nNx2ContingencyTable::row_proby(int) const -> " + << "range check error\n\n"; + exit(1); + } - mlog << Error << "\nNx2ContingencyTable::row_nonevent_likelihood(int) const -> " - << "range check error\n\n"; + double v; - exit ( 1 ); + if(use_center) v = 0.5*(Thresholds[row] + Thresholds[row + 1]); + else v = Thresholds[row]; + return v; } -double x, num, denom; - -denom = (double) nonevent_col_total(); +//////////////////////////////////////////////////////////////////////// -num = (double) nonevent_count_by_row(row); +double Nx2ContingencyTable::row_calibration(int row) const { -if(is_eq(denom, 0.0)) x = bad_data_double; -else x = num/denom; + if(row < 0 || row >= Nrows) { + mlog << Error << "\nNx2ContingencyTable::row_calibration(int) const -> " + << "range check error\n\n"; + exit(1); + } -return x; + double num = event_count_by_row(row); + double den = num + nonevent_count_by_row(row); + return compute_proportion(num, den); } - //////////////////////////////////////////////////////////////////////// +double Nx2ContingencyTable::row_refinement(int row) const { -TTContingencyTable Nx2ContingencyTable::ctc_by_row(int row) const - -{ - -TTContingencyTable tt; - -if ( (row < 0) || (row >= Nrows) ) { - - mlog << Error << "\nNx2ContingencyTable::ctc_by_row(int) const -> " - << "range check error\n\n"; - - exit ( 1 ); + if(row < 0 || row >= Nrows) { + mlog << Error << "\nNx2ContingencyTable::row_refinement(int) const -> " + << "range check error\n\n"; + exit(1); + } + return compute_proportion( event_count_by_row(row) + + nonevent_count_by_row(row), n()); } -int j; -int sy, sn; - - /////////////////// - -sy = sn = 0; +//////////////////////////////////////////////////////////////////////// -for (j=(row + 1); j= Nrows) { + mlog << Error << "\nNx2ContingencyTable::row_event_likelihood(int) const -> " + << "range check error\n\n"; + exit(1); + } + return compute_proportion(event_count_by_row(row), + event_col_total()); } -tt.set_fy_oy(sy); -tt.set_fy_on(sn); - - /////////////////// - -sy = sn = 0; - -for (j=0; j<=row; ++j) { - - sy += event_count_by_row(j); - sn += nonevent_count_by_row(j); - -} +//////////////////////////////////////////////////////////////////////// -tt.set_fn_oy(sy); -tt.set_fn_on(sn); - /////////////////// +double Nx2ContingencyTable::row_nonevent_likelihood(int row) const { -return tt; + if(row < 0 || row >= Nrows) { + mlog << Error << "\nNx2ContingencyTable::row_nonevent_likelihood(int) const -> " + << "range check error\n\n"; + exit(1); + } + return compute_proportion(nonevent_count_by_row(row), + nonevent_col_total()); } - //////////////////////////////////////////////////////////////////////// +TTContingencyTable Nx2ContingencyTable::ctc_by_row(int row) const { + TTContingencyTable tt; -double Nx2ContingencyTable::roc_auc() const + if(row < 0 || row >= Nrows) { + mlog << Error << "\nNx2ContingencyTable::ctc_by_row(int) const -> " + << "range check error\n\n"; + exit(1); + } -{ + double sy = 0.0; + double sn = 0.0; -int j; -TTContingencyTable ct; -double area, x_prev, y_prev, x, y; + for(int j=row+1; j 1 so that degf > 0 in the call to gsl_cdf_tdist_Pinv() - -if(is_bad_data(bs = brier_score()) || N <= 1) return bad_data_double; - -degf = N - 1.0; - -t = gsl_cdf_tdist_Pinv(1.0 - 0.5*alpha, degf); - -ob = obar(); -af1 = sf2 = sf3 = af4 = 0.0; +double Nx2ContingencyTable::roc_auc() const { + TTContingencyTable ct; + double x_prev = 1.0; + double y_prev = 1.0; + double area = bad_data_double; -for (j=0; j " << "table not square!\n\n"; - - exit ( 1 ); - + exit(1); } - for(i=0, num=0.0; i " - << "table not square!\n\n"; - - exit ( 1 ); - -} - -const int N = total(); - - // - // MET #2542: return bad data for empty tables rather than erroring out - // - -if ( N == 0 ) return bad_data_double; - -const double DN = (double) N; -int j, k, m, n; -double num, denom, sum, ans; - - // - // first term in numerator - // - -sum = 0.0; - -for (j=0; j " + << "table not square!\n\n"; + exit(1); + } - n = rc_to_n(j, j); + // MET #2542: return bad data for empty tables rather than erroring out + if(E.empty()) return bad_data_double; - sum += ((*E)[n])/DN; + // First term in numerator + const double DN = total(); + double sum = 0.0; -} + for(int j=0; j " - << "table not square!\n\n"; - - exit ( 1 ); - -} - -if ( ec_value < 0.0 || ec_value >= 1.0 ) { - - mlog << Error << "\nContingencyTable::gheidke_ec(double) -> " - << "ec_value (" << ec_value << ") must be >=0 and <1.0!\n\n"; - - exit ( 1 ); - -} - -const int N = total(); - - // - // MET #2542: return bad data for empty tables rather than erroring out - // - -if ( N == 0 ) return bad_data_double; - -int j, sum; -double num, denom, ans; - - // - // sum counts on the diagonal - // - -for (j=0, sum=0; j " + << "table not square!\n\n"; + exit(1); + } -num = (double) sum - ec; -denom = (double) N - ec; + if(ec_value < 0.0 || ec_value >= 1.0) { + mlog << Error << "\nContingencyTable::gheidke_ec() -> " + << "ec_value (" << ec_value << ") must be >=0 and <1.0!\n\n"; + exit(1); + } - // - // result - // + // MET #2542: return bad data for empty tables rather than erroring out + if(E.empty()) return bad_data_double; -if (is_eq(denom, 0.0)) ans = bad_data_double; -else ans = num/denom; + // Sum entries on the diagonal + double sum = 0.0; + for(int j=0; j " - << "table not square!\n\n"; - - exit ( 1 ); - -} - -const int N = total(); - - // - // MET #2542: return bad data for empty tables rather than erroring out - // - -if ( N == 0 ) return bad_data_double; - -const double DN = (double) N; -int j, k, m, n; -double num, denom, sum, t, ans; - - // - // first term in numerator - // - -sum = 0.0; - -for (j=0; j " - << "table not square!\n\n"; - - exit ( 1 ); - -} - -const int N = total(); - - // - // MET #2542: return bad data for empty tables rather than erroring out - // - -if ( N == 0 ) return bad_data_double; - -int j, k, m, n; -const double DN = (double) N; -double t, sum; -double * p = (double *) nullptr; -double * s = (double *) nullptr; + if(Nrows != Ncols) { + mlog << Error << "\nContingencyTable::gkuiper() -> " + << "table not square!\n\n"; + exit(1); + } - // - // can't compute gerrity when the first column contains all zeros - // + // MET #2542: return bad data for empty tables rather than erroring out + if(E.empty()) return bad_data_double; -if ( col_total(0) == 0 ) return bad_data_double; + const double DN = total(); -p = new double [Nrows]; + // First term in numerator + double sum = 0.0; + for(int j=0; j " + << "table not square!\n\n"; + exit(1); } -} + // MET #2542: return bad data for empty tables rather than erroring out + if(E.empty()) return bad_data_double; - // - // replace nan with bad data - // + double DN = total(); - if (std::isnan(sum)) sum = bad_data_double; + // Can't compute gerrity when the first column contains all zeros + if(is_eq(col_total(0), 0.0)) return bad_data_double; - // - // done - // + // the p array + vector p(Nrows); + for(int j=0; j s(Nrows*Nrows); + calc_gerrity_scoring_matrix(Nrows, p, s); + + // Calculate score + double sum = 0.0; + for(int j=0; j &p, + vector &s) { + double b = 1.0/(N - 1.0); -void calc_gerrity_scoring_matrix(int N, const double * p, double * s) - -{ - -int j, k, n; -double b, t, sum; -double * a = (double *) nullptr; -double * recip_sum = (double *) nullptr; -double * direct_sum = (double *) nullptr; - - -a = new double [N]; -recip_sum = new double [N]; -direct_sum = new double [N]; - -b = 1.0/(N - 1.0); - - - // - // the a array - // - -sum = 0.0; - -for (j=0; j a(N); + double sum = 0.0; + for(int j=0; j recip_sum(N); + recip_sum[0] = 0.0; + sum = 0.0; + for(int j=1; j=0; --j) { - - sum += a[j]; - - direct_sum[j] = sum; - -} - - // - // entries of the scoring matrix - // - -for (j=0; j direct_sum(N); + direct_sum[N - 1] = 0.0; + sum = 0.0; + for(int j=(N - 2); j>=0; --j) { + sum += a[j]; + direct_sum[j] = sum; } -} + // Entries of the scoring matrix + for(int j=0; j thresh(n); for(i=0; i