root/library/tests/epdf_harness.cpp @ 690

Revision 675, 7.9 kB (checked in by mido, 15 years ago)

experiment: epdf as a descendat of mpdf

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 ( R.rows() > 0 ) {
62                check_sample_mean();
63                check_covariance();
64        }
65
66        if ( mrv ) {
67                RV crv = hepdf->_rv().subt ( *mrv );
68                epdf_ptr m = hepdf->marginal ( *mrv );
69                mpdf_ptr c = hepdf->condition ( crv );
70
71                mpdf_array aa ( 2 );
72                aa ( 0 ) = c;
73                aa ( 1 ) = m;
74                mprod mEp ( aa );
75
76                check_cond_mean(mEp);
77                if ( R.rows() > 0 ) {
78                        check_cond_covariance( mEp );
79                }
80
81                // test of pdflog at zero
82                vec zero ( 0 );
83                vec zeron=zeros ( hepdf->dimension() );
84
85                double lpz = hepdf->evallog ( zeron );
86                double lpzc=mEp.evallogcond ( zeron, zero );
87                CHECK_CLOSE_EX ( lpz, lpzc, tolerance );
88
89                vec lpzv(1);
90                lpzv(0) = lpz;
91
92                mat zero1n ( hepdf->dimension(), 1 );
93                for ( int i = 0; i < zero1n.rows(); ++i ) {
94                        zero1n ( i, 0 ) = 0;
95                }
96
97                vec lpzv_act = hepdf->evallog_m ( zero1n );
98                CHECK_CLOSE_EX ( lpzv, lpzv_act, tolerance );
99
100                Array<vec> zeroa(3);
101                lpzv = vec( zeroa.size() );
102                for ( int i = 0; i < zeroa.size(); ++i ) {
103                        zeroa(i) = zeron;
104                        lpzv(i) = lpz;
105                }
106
107                lpzv_act = hepdf->evallog_m ( zeroa );
108                CHECK_CLOSE_EX ( lpzv, lpzv_act, tolerance );
109        }
110}
111
112void epdf_harness::check_support_mean() {
113        vec xb = support.get_row ( 0 );
114        vec yb = support.get_row ( 1 );
115
116        int tc = 0;
117        Array<vec> actual(CurrentContext::max_trial_count);
118        do {
119                vec emu = num_mean2 ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
120                actual( tc ) = emu;
121                ++tc;
122        } while ( ( tc < CurrentContext::max_trial_count ) &&
123                  !UnitTest::AreClose ( mean, actual( tc - 1 ), tolerance ) );
124        if ( ( tc == CurrentContext::max_trial_count ) &&
125             ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
126                UnitTest::MemoryOutStream stream;
127                stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << tolerance << " but was " << actual;
128
129                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);
130
131                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
132       }
133}
134
135void epdf_harness::check_support_integral() {
136        vec xb = support.get_row ( 0 );
137        vec yb = support.get_row ( 1 );
138
139        int tc = 0;
140        Array<double> actual(CurrentContext::max_trial_count);
141        do {
142                double nc = normcoef ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
143                actual( tc ) = nc;
144                ++tc;
145        } while ( ( tc < CurrentContext::max_trial_count ) &&
146                  !UnitTest::AreClose ( 1.0, actual( tc - 1 ), tolerance ) );
147        if ( ( tc == CurrentContext::max_trial_count ) &&
148             ( !UnitTest::AreClose ( 1.0, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
149                UnitTest::MemoryOutStream stream;
150                stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << tolerance << " but was " << actual;
151
152                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);
153
154                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
155       }
156}
157
158void epdf_harness::check_sample_mean() {
159        vec delta = make_close_tolerance ( variance, nsamples );
160
161        int tc = 0;
162        Array<vec> actual(CurrentContext::max_trial_count);
163        do {
164                mat smp = hepdf->sample_m ( nsamples );
165                vec emu = smp * ones ( nsamples ) / nsamples;
166                actual( tc ) = emu;
167                ++tc;
168        } while ( ( tc < CurrentContext::max_trial_count ) &&
169                  !UnitTest::AreClose ( mean, actual( tc - 1 ), delta ) );
170        if ( ( tc == CurrentContext::max_trial_count ) &&
171             ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), delta ) ) ) {
172                UnitTest::MemoryOutStream stream;
173                stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << delta << " but was " << actual;
174
175                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);
176
177                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
178       }
179}
180
181void epdf_harness::check_covariance() {
182        int tc = 0;
183        Array<mat> actual(CurrentContext::max_trial_count);
184        do {
185                mat smp = hepdf->sample_m ( nsamples );
186                vec emu = smp * ones ( nsamples ) / nsamples;
187                mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
188                actual( tc ) = er;
189                ++tc;
190        } while ( ( tc < CurrentContext::max_trial_count ) &&
191                  !UnitTest::AreClose ( R, actual( tc - 1 ), tolerance ) );
192       
193        if ( ( tc == CurrentContext::max_trial_count ) &&
194             ( !UnitTest::AreClose ( R, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
195                UnitTest::MemoryOutStream stream;
196                stream << CurrentContext::format_context(__LINE__) << "expected " << R << " +/- " << tolerance << " but was " << actual;
197
198                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);
199
200                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
201       }
202}
203
204void epdf_harness::check_cond_mean( mprod &mep ) {
205        vec delta = make_close_tolerance ( variance, nsamples );
206
207        int tc = 0;
208        Array<vec> actual(CurrentContext::max_trial_count);
209        do {
210                mat smp = mep.samplecond_m ( vec ( 0 ), nsamples );
211                vec emu = sum ( smp, 2 ) / nsamples;
212                actual( tc ) = emu;
213                ++tc;
214        } while ( ( tc < CurrentContext::max_trial_count ) &&
215                  !UnitTest::AreClose ( mean, actual( tc - 1 ), delta ) );
216        if ( ( tc == CurrentContext::max_trial_count ) &&
217             ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), delta ) ) ) {
218                UnitTest::MemoryOutStream stream;
219                stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << delta << " but was " << actual;
220
221                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);
222
223                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
224       }
225}
226
227void epdf_harness::check_cond_covariance( mprod &mep ) {
228        int tc = 0;
229        Array<mat> actual(CurrentContext::max_trial_count);
230        do {
231                mat smp = mep.samplecond_m ( vec ( 0 ), nsamples );
232                vec emu = sum ( smp, 2 ) / nsamples;
233                mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
234                actual( tc ) = er;
235                ++tc;
236        } while ( ( tc < CurrentContext::max_trial_count ) &&
237                  !UnitTest::AreClose ( R, actual( tc - 1 ), tolerance ) );
238        if ( ( tc == CurrentContext::max_trial_count ) &&
239             ( !UnitTest::AreClose ( R, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
240                UnitTest::MemoryOutStream stream;
241                stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << tolerance << " but was " << actual;
242
243                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);
244
245                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
246       }
247}
248
249}
Note: See TracBrowser for help on using the browser.