root/library/tests/epdf_harness.cpp @ 493

#include "epdf_harness.h"
#include "base/bdmbase.h"
#include "base/user_info.h"
#include "stat/exp_family.h"
#include "mat_checks.h"
#include "test_util.h"
#include "UnitTest++.h"
#include <memory>

namespace bdm {

void epdf_harness::test_config ( const char *config_file_name ) {
        RV::clear_all();

        UIFile in ( config_file_name );
        Array<epdf_harness *> input;
        UI::get ( input, in, "data", UI::compulsory );
        int sz = input.size();
        CHECK ( sz > 0 );
        for ( int i = 0; i < sz; ++i ) {
                input ( i )->test ( config_file_name, i );
        }
}

void epdf_harness::from_setting ( const Setting &set ) {
        hepdf = UI::build<epdf> ( set, "epdf", UI::compulsory );
        UI::get ( mean, set, "mean", UI::compulsory );
        UI::get ( variance, set, "variance", UI::compulsory );

        UI::get (support, set, "support", UI::optional );
        UI::get (nbins, set, "nbins", UI::optional );
        UI::get (nsamples, set, "nsamples", UI::optional );
        UI::get (R, set, "R", UI::optional );

        mrv = UI::build<RV> (set, "marginal_rv", UI::optional );

        UI::get ( tolerance, set, "tolerance", UI::optional );
}

void epdf_harness::test ( const char *config_name, int idx ) {
        CurrentContext cc ( config_name, idx );

        CHECK_CLOSE_EX ( mean, hepdf->mean(), tolerance );
        CHECK_CLOSE_EX ( variance, hepdf->variance(), tolerance );

        if ( support.rows() == 2 ) {
                int old_size = nbins.size();
                if ( old_size < 2 ) {
                        ivec new_nbins ( "100 100" );
                        for ( int i = 0; i < old_size; ++i ) {
                                new_nbins ( i ) = nbins ( i );
                        }

                        nbins = new_nbins;
                }

                check_support_mean();
                check_support_integral();
        }

        if ( R.rows() > 0 ) {
                check_sample_mean();
                check_covariance();
        }

        if ( mrv.get() ) {
                RV crv = hepdf->_rv().subt ( *mrv );
                shared_ptr<epdf> m = hepdf->marginal ( *mrv );
                shared_ptr<mpdf> c = hepdf->condition ( crv );
                mepdf mm ( m );

                Array<mpdf *> aa ( 2 );
                aa ( 0 ) = c.get();
                aa ( 1 ) = &mm;
                mprod mEp ( aa );

                check_cond_mean(mEp);
                if ( R.rows() > 0 ) {
                        check_cond_covariance( mEp );
                }

                // test of pdflog at zero
                vec zero ( 0 );
                vec zeron ( hepdf->dimension() );
                for ( int i = 0; i < zeron.size(); ++i ) {
                        zeron ( i ) = 0;
                }

                CHECK_CLOSE_EX ( hepdf->evallog ( zeron ), mEp.evallogcond ( zeron, zero ), tolerance );
        }
}

void epdf_harness::check_support_mean() {
        vec xb = support.get_row ( 0 );
        vec yb = support.get_row ( 1 );

        int tc = 0;
        Array<vec> actual(CurrentContext::max_trial_count);
        do {
                vec emu = num_mean2 ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
                actual( tc ) = emu;
                ++tc;
        } while ( ( tc < CurrentContext::max_trial_count ) &&
                  !UnitTest::AreClose ( mean, actual( tc - 1 ), tolerance ) );
        if ( ( tc == CurrentContext::max_trial_count ) &&
             ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
                UnitTest::MemoryOutStream stream;
                stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << tolerance << " but was " << actual;

                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);

                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
       }
}

void epdf_harness::check_support_integral() {
        vec xb = support.get_row ( 0 );
        vec yb = support.get_row ( 1 );

        int tc = 0;
        Array<double> actual(CurrentContext::max_trial_count);
        do {
                double nc = normcoef ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
                actual( tc ) = nc;
                ++tc;
        } while ( ( tc < CurrentContext::max_trial_count ) &&
                  !UnitTest::AreClose ( 1.0, actual( tc - 1 ), tolerance ) );
        if ( ( tc == CurrentContext::max_trial_count ) &&
             ( !UnitTest::AreClose ( 1.0, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
                UnitTest::MemoryOutStream stream;
                stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << tolerance << " but was " << actual;

                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);

                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
       }
}

void epdf_harness::check_sample_mean() {
        vec delta = make_close_tolerance ( variance, nsamples );

        int tc = 0;
        Array<vec> actual(CurrentContext::max_trial_count);
        do {
                mat smp = hepdf->sample_m ( nsamples );
                vec emu = smp * ones ( nsamples ) / nsamples;
                actual( tc ) = emu;
                ++tc;
        } while ( ( tc < CurrentContext::max_trial_count ) &&
                  !UnitTest::AreClose ( mean, actual( tc - 1 ), delta ) );
        if ( ( tc == CurrentContext::max_trial_count ) &&
             ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), delta ) ) ) {
                UnitTest::MemoryOutStream stream;
                stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << delta << " but was " << actual;

                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);

                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
       }
}

void epdf_harness::check_covariance() {
        int tc = 0;
        Array<mat> actual(CurrentContext::max_trial_count);
        do {
                mat smp = hepdf->sample_m ( nsamples );
                vec emu = smp * ones ( nsamples ) / nsamples;
                mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
                actual( tc ) = er;
                ++tc;
        } while ( ( tc < CurrentContext::max_trial_count ) &&
                  !UnitTest::AreClose ( R, actual( tc - 1 ), tolerance ) );
        
        if ( ( tc == CurrentContext::max_trial_count ) &&
             ( !UnitTest::AreClose ( R, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
                UnitTest::MemoryOutStream stream;
                stream << CurrentContext::format_context(__LINE__) << "expected " << R << " +/- " << tolerance << " but was " << actual;

                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);

                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
       }
}

void epdf_harness::check_cond_mean( mprod &mep ) {
        vec delta = make_close_tolerance ( variance, nsamples );

        int tc = 0;
        Array<vec> actual(CurrentContext::max_trial_count);
        do {
                mat smp = mep.samplecond ( vec ( 0 ), nsamples );
                vec emu = sum ( smp, 2 ) / nsamples;
                actual( tc ) = emu;
                ++tc;
        } while ( ( tc < CurrentContext::max_trial_count ) &&
                  !UnitTest::AreClose ( mean, actual( tc - 1 ), delta ) );
        if ( ( tc == CurrentContext::max_trial_count ) &&
             ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), delta ) ) ) {
                UnitTest::MemoryOutStream stream;
                stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << delta << " but was " << actual;

                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);

                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
       }
}

void epdf_harness::check_cond_covariance( mprod &mep ) {
        int tc = 0;
        Array<mat> actual(CurrentContext::max_trial_count);
        do {
                mat smp = mep.samplecond ( vec ( 0 ), nsamples );
                vec emu = sum ( smp, 2 ) / nsamples;
                mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
                actual( tc ) = er;
                ++tc;
        } while ( ( tc < CurrentContext::max_trial_count ) &&
                  !UnitTest::AreClose ( R, actual( tc - 1 ), tolerance ) );
        if ( ( tc == CurrentContext::max_trial_count ) &&
             ( !UnitTest::AreClose ( R, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
                UnitTest::MemoryOutStream stream;
                stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << tolerance << " but was " << actual;

                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);

                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
       }
}

}