root/library/tests/epdf_harness.cpp @ 725

Revision 725, 6.7 kB (checked in by smidl, 15 years ago)

Sampling from egiw

  • Property svn:eol-style set to native
Line 
1#include "epdf_harness.h"
2#include "base/bdmbase.h"
3#include "base/user_info.h"
4#include "stat/exp_family.h"
5#include "mat_checks.h"
6#include "test_util.h"
7#include "UnitTest++.h"
8#include <memory>
9
10namespace bdm {
11
12void epdf_harness::test_config ( const char *config_file_name ) {
13        RV::clear_all();
14
15        UIFile in ( config_file_name );
16        Array<shared_ptr<epdf_harness> > input;
17        UI::get ( input, in, "data", UI::compulsory );
18        int sz = input.size();
19        CHECK ( sz > 0 );
20        for ( int i = 0; i < sz; ++i ) {
21                input ( i )->test ( config_file_name, i );
22        }
23}
24
25void epdf_harness::from_setting ( const Setting &set ) {
26        hepdf = UI::build<epdf> ( set, "epdf", UI::compulsory );
27        UI::get ( mean, set, "mean", UI::compulsory );
28        UI::get ( variance, set, "variance", UI::compulsory );
29
30        UI::get ( support, set, "support", UI::optional );
31        UI::get ( nbins, set, "nbins", UI::optional );
32        UI::get ( nsamples, set, "nsamples", UI::optional );
33        UI::get ( R, set, "R", UI::optional );
34
35        mrv = UI::build<RV> ( set, "marginal_rv", UI::optional );
36
37        UI::get ( tolerance, set, "tolerance", UI::optional );
38}
39
40void epdf_harness::test ( const char *config_name, int idx ) {
41        CurrentContext cc ( config_name, idx );
42
43        CHECK_CLOSE_EX ( mean, hepdf->mean(), tolerance );
44        CHECK_CLOSE_EX ( variance, hepdf->variance(), tolerance );
45
46        if ( support.rows() == 2 ) {
47                int old_size = nbins.size();
48                if ( old_size < 2 ) {
49                        ivec new_nbins ( "100 100" );
50                        for ( int i = 0; i < old_size; ++i ) {
51                                new_nbins ( i ) = nbins ( i );
52                        }
53
54                        nbins = new_nbins;
55                }
56
57                check_support_mean();
58                check_support_integral();
59        }
60
61        if ( variance.length() > 0 ) {
62                check_sample_mean();
63        }
64        if (R.rows()>0){
65                check_covariance();
66        }
67
68        if ( mrv ) {
69                RV crv = hepdf->_rv().subt ( *mrv );
70                epdf_ptr m = hepdf->marginal ( *mrv );
71                pdf_ptr c = hepdf->condition ( crv );
72
73                pdf_array aa ( 2 );
74                aa ( 0 ) = c;
75                aa ( 1 ) = m;
76                mprod mEp ( aa );
77
78                check_cond_mean ( mEp );
79                if ( R.rows() > 0 ) {
80                        check_cond_covariance ( mEp );
81                }
82
83                // test of pdflog at zero
84                vec zero ( 0 );
85                vec zeron = zeros ( hepdf->dimension() );
86
87                double lpz = hepdf->evallog ( zeron );
88                double lpzc = mEp.evallogcond ( zeron, zero );
89                CHECK_CLOSE_EX ( lpz, lpzc, tolerance );
90
91                vec lpzv ( 1 );
92                lpzv ( 0 ) = lpz;
93
94                mat zero1n ( hepdf->dimension(), 1 );
95                for ( int i = 0; i < zero1n.rows(); ++i ) {
96                        zero1n ( i, 0 ) = 0;
97                }
98
99                vec lpzv_act = hepdf->evallog_mat ( zero1n );
100                CHECK_CLOSE_EX ( lpzv, lpzv_act, tolerance );
101
102                Array<vec> zeroa ( 3 );
103                lpzv = vec ( zeroa.size() );
104                for ( int i = 0; i < zeroa.size(); ++i ) {
105                        zeroa ( i ) = zeron;
106                        lpzv ( i ) = lpz;
107                }
108
109                lpzv_act = hepdf->evallog_mat ( zeroa );
110                CHECK_CLOSE_EX ( lpzv, lpzv_act, tolerance );
111        }
112}
113
114void epdf_harness::check_support_mean() {
115        vec xb = support.get_row ( 0 );
116        vec yb = support.get_row ( 1 );
117
118        vec actual;
119        actual = num_mean2 ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
120
121        CHECK_CLOSE(mean, actual, tolerance);
122}
123
124void epdf_harness::check_support_integral() {
125        vec xb = support.get_row ( 0 );
126        vec yb = support.get_row ( 1 );
127
128        double nc = normcoef ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
129        CHECK_CLOSE(1.0,nc,0.01);
130}
131
132void epdf_harness::check_sample_mean() {
133        vec delta = make_close_tolerance ( variance, nsamples );
134
135        int tc = 0;
136        Array<vec> actual ( CurrentContext::max_trial_count );
137        do {
138                mat smp = hepdf->sample_mat ( nsamples );
139                vec emu = smp * ones ( nsamples ) / nsamples;
140                actual ( tc ) = emu;
141                ++tc;
142        } while ( ( tc < CurrentContext::max_trial_count ) &&
143                  !UnitTest::AreClose ( mean, actual ( tc - 1 ), delta ) );
144                         
145        if ( ( tc == CurrentContext::max_trial_count ) &&
146                ( !UnitTest::AreClose ( mean, actual ( CurrentContext::max_trial_count - 1 ), delta ) ) ) {
147                UnitTest::MemoryOutStream stream;
148                stream << CurrentContext::format_context ( __LINE__ ) << "expected " << mean << " +/- " << delta << " but was " << actual;
149
150                UnitTest::TestDetails details ( *UnitTest::CurrentTest::Details(), 0, false );
151
152                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
153        }
154}
155
156void epdf_harness::check_covariance() {
157        int tc = 0;
158        Array<mat> actual ( CurrentContext::max_trial_count );
159        do {
160                mat smp = hepdf->sample_mat ( nsamples );
161                vec emu = smp * ones ( nsamples ) / nsamples;
162                mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
163                actual ( tc ) = er;
164                ++tc;
165        } while ( ( tc < CurrentContext::max_trial_count ) &&
166                  !UnitTest::AreClose ( R, actual ( tc - 1 ), tolerance ) );
167
168        if ( ( tc == CurrentContext::max_trial_count ) &&
169                ( !UnitTest::AreClose ( R, actual ( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
170                UnitTest::MemoryOutStream stream;
171                stream << CurrentContext::format_context ( __LINE__ ) << "expected " << R << " +/- " << tolerance << " but was " << actual;
172
173                UnitTest::TestDetails details ( *UnitTest::CurrentTest::Details(), 0, false );
174
175                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
176        }
177}
178
179void epdf_harness::check_cond_mean ( mprod &mep ) {
180        vec delta = make_close_tolerance ( variance, nsamples );
181
182        int tc = 0;
183        Array<vec> actual ( CurrentContext::max_trial_count );
184        do {
185                mat smp = mep.samplecond_mat ( vec ( 0 ), nsamples );
186                vec emu = sum ( smp, 2 ) / nsamples;
187                actual ( tc ) = emu;
188                ++tc;
189        } while ( ( tc < CurrentContext::max_trial_count ) &&
190                  !UnitTest::AreClose ( mean, actual ( tc - 1 ), delta ) );
191        if ( ( tc == CurrentContext::max_trial_count ) &&
192                ( !UnitTest::AreClose ( mean, actual ( CurrentContext::max_trial_count - 1 ), delta ) ) ) {
193                UnitTest::MemoryOutStream stream;
194                stream << CurrentContext::format_context ( __LINE__ ) << "expected " << mean << " +/- " << delta << " but was " << actual;
195
196                UnitTest::TestDetails details ( *UnitTest::CurrentTest::Details(), 0, false );
197
198                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
199        }
200}
201
202void epdf_harness::check_cond_covariance ( mprod &mep ) {
203        int tc = 0;
204        Array<mat> actual ( CurrentContext::max_trial_count );
205        do {
206                mat smp = mep.samplecond_mat ( vec ( 0 ), nsamples );
207                vec emu = sum ( smp, 2 ) / nsamples;
208                mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
209                actual ( tc ) = er;
210                ++tc;
211        } while ( ( tc < CurrentContext::max_trial_count ) &&
212                  !UnitTest::AreClose ( R, actual ( tc - 1 ), tolerance ) );
213        if ( ( tc == CurrentContext::max_trial_count ) &&
214                ( !UnitTest::AreClose ( R, actual ( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
215                UnitTest::MemoryOutStream stream;
216                stream << CurrentContext::format_context ( __LINE__ ) << "expected " << mean << " +/- " << tolerance << " but was " << actual;
217
218                UnitTest::TestDetails details ( *UnitTest::CurrentTest::Details(), 0, false );
219
220                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
221        }
222}
223
224}
Note: See TracBrowser for help on using the browser.