root/library/bdm/base/bdmbase.cpp @ 949

Revision 949, 14.9 kB (checked in by mido, 14 years ago)

another patch of log_levels in user_info.h where a small piece of code was forgotten in the previous facelift of UI
also, DS adjusted to a new UI::GET policy concerning RVs


  • Property svn:eol-style set to native
Line 
1
2#include "bdmbase.h"
3
4//! Space of basic BDM structures
5namespace bdm {
6
7const int RV::BUFFER_STEP = 1;
8
9Array<string> RV::NAMES ( RV::BUFFER_STEP );
10
11ivec RV::SIZES ( RV::BUFFER_STEP );
12
13RV::str2int_map RV::MAP;
14
15void RV::clear_all() {
16        MAP.clear();
17        SIZES.clear();
18        NAMES = Array<string> ( BUFFER_STEP );
19}
20
21string RV::show_all() {
22        ostringstream os;
23        for ( str2int_map::const_iterator iter = MAP.begin(); iter != MAP.end(); iter++ ) {
24                os << "key: " << iter->first << " val: " << iter->second << endl;
25        }
26        return os.str();
27};
28
29int RV::assign_id( const string &name, int size ) {
30        //Refer
31        int id;
32        str2int_map::const_iterator iter = MAP.find ( name );
33        if ( iter == MAP.end() || name.length() == 0 ) { //add new RV
34                id = MAP.size() + 1;
35                //debug
36                /*              {
37                                        cout << endl;
38                                        str2int_map::const_iterator iter = MAP.begin();
39                                        for(str2int_map::const_iterator iter=MAP.begin(); iter!=MAP.end(); iter++){
40                                                cout << "key: " << iter->first << " val: " << iter->second <<endl;
41                                        }
42                                }*/
43
44                MAP.insert ( make_pair ( name, id ) ); //add new rv
45                if ( id >= NAMES.length() ) {
46                        NAMES.set_length ( id + BUFFER_STEP, true );
47                        SIZES.set_length ( id + BUFFER_STEP, true );
48                }
49                NAMES ( id ) = name;
50                SIZES ( id ) = size;
51                bdm_assert ( size > 0, "RV " + name + " does not exists. Default size (-1) can not be assigned " );
52        } else {
53                id = iter->second;
54                if ( size > 0 && name.length() > 0 ) {
55                        bdm_assert ( SIZES ( id ) == size, "RV " + name + " of size " + num2str ( SIZES ( id ) ) + " exists, requested size " + num2str ( size ) + "can not be assigned" );
56                }
57        }
58        return id;
59};
60
61int RV::countsize() const {
62        int tmp = 0;
63        for ( int i = 0; i < len; i++ ) {
64                tmp += SIZES ( ids ( i ) );
65        }
66        return tmp;
67}
68
69ivec RV::cumsizes() const {
70        ivec szs ( len );
71        int tmp = 0;
72        for ( int i = 0; i < len; i++ ) {
73                tmp += SIZES ( ids ( i ) );
74                szs ( i ) = tmp;
75        }
76        return szs;
77}
78
79void RV::init ( const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times ) {
80        len = in_names.length();
81        bdm_assert ( in_names.length() == in_times.length(), "check \"times\" " );
82        bdm_assert ( in_names.length() == in_sizes.length(), "check \"sizes\" " );
83
84        times.set_length ( len );
85        ids.set_length ( len );
86        int id;
87        for ( int i = 0; i < len; i++ ) {
88                id = assign_id ( in_names ( i ), in_sizes ( i ) );
89                ids ( i ) = id;
90        }
91        times = in_times;
92        dsize = countsize();
93}
94
95RV::RV ( string name, int sz, int tm ) {
96        Array<string> A ( 1 );
97        A ( 0 ) = name;
98        init ( A, vec_1 ( sz ), vec_1 ( tm ) );
99}
100
101RV::RV ( int sz, int tm ) {
102        Array<string> A ( 1 );
103        A ( 0 ) = "";
104        init ( A, vec_1 ( sz ), vec_1 ( tm ) );
105}
106
107bool RV::add ( const RV &rv2 ) {
108        if ( rv2.len > 0 ) { //rv2 is nonempty
109                ivec ind = rv2.findself ( *this ); //should be -1 all the time
110                ivec index = itpp::find ( ind == -1 );
111
112                if ( index.length() < rv2.len ) { //conflict
113                        ids = concat ( ids, rv2.ids ( index ) );
114                        times = concat ( times, rv2.times ( index ) );
115                } else {
116                        ids = concat ( ids, rv2.ids );
117                        times = concat ( times, rv2.times );
118                }
119                len = ids.length();
120                dsize = countsize();
121                return ( index.length() == rv2.len ); //conflict or not
122        } else { //rv2 is empty
123                return true; // no conflict
124        }
125};
126
127RV RV::subselect ( const ivec &ind ) const {
128        RV ret;
129        ret.ids = ids ( ind );
130        ret.times = times ( ind );
131        ret.len = ind.length();
132        ret.dsize = ret.countsize();
133        return ret;
134}
135
136RV RV::operator() ( int di1, int di2 ) const {
137        ivec sz = cumsizes();
138        int i1 = 0;
139        while ( sz ( i1 ) < di1 ) i1++;
140        int i2 = i1;
141        while ( sz ( i2 ) < di2 ) i2++;
142        return subselect ( linspace ( i1, i2 ) );
143}
144
145void RV::t_plus ( int delta ) {
146        times += delta;
147}
148
149bool RV::equal ( const RV &rv2 ) const {
150        return ( ids == rv2.ids ) && ( times == rv2.times );
151}
152
153shared_ptr<pdf> epdf::condition ( const RV &rv ) const NOT_IMPLEMENTED( shared_ptr<pdf>() );
154
155
156shared_ptr<epdf> epdf::marginal ( const RV &rv ) const NOT_IMPLEMENTED( shared_ptr<epdf>() );
157
158mat epdf::sample_mat ( int N ) const {
159        mat X = zeros ( dim, N );
160        for ( int i = 0; i < N; i++ ) X.set_col ( i, this->sample() );
161        return X;
162}
163
164vec epdf::evallog_mat ( const mat &Val ) const {
165        vec x ( Val.cols() );
166        for ( int i = 0; i < Val.cols(); i++ ) {
167                x ( i ) = evallog ( Val.get_col ( i ) );
168        }
169
170        return x;
171}
172
173vec epdf::evallog_mat ( const Array<vec> &Avec ) const {
174        vec x ( Avec.size() );
175        for ( int i = 0; i < Avec.size(); i++ ) {
176                x ( i ) = evallog ( Avec ( i ) );
177        }
178
179        return x;
180}
181
182mat pdf::samplecond_mat ( const vec &cond, int N ) {
183        mat M ( dimension(), N );
184        for ( int i = 0; i < N; i++ ) {
185                M.set_col ( i, samplecond ( cond ) );
186        }
187
188        return M;
189}
190
191void pdf::from_setting ( const Setting &set ) {
192        root::from_setting( set );
193        shared_ptr<RV> r = UI::build<RV> ( set, "rv", UI::optional );
194        if ( r ) {
195                set_rv ( *r );
196        }
197
198        r = UI::build<RV> ( set, "rvc", UI::optional );
199        if ( r ) {
200                set_rvc ( *r );
201        }
202}
203
204void pdf::to_setting ( Setting &set ) const {   
205        root::to_setting( set );
206        UI::save( &rv, set, "rv" );
207        UI::save( &rvc, set, "rvc" );
208}
209
210void datalink::set_connection ( const RV &rv, const RV &rv_up ) {
211        downsize = rv._dsize();
212        upsize = rv_up._dsize();
213        v2v_up = rv.dataind ( rv_up );
214        bdm_assert_debug ( v2v_up.length() == downsize, string("rv (" + rv.to_string() + ") is not fully in rv_up(" + rv_up.to_string()) );
215}
216
217void datalink::set_connection ( int ds, int us, const ivec &upind ) {
218        downsize = ds;
219        upsize = us;
220        v2v_up = upind;
221        bdm_assert_debug ( v2v_up.length() == downsize, "rv is not fully in rv_up" );
222}
223
224void datalink_part::set_connection ( const RV &rv, const RV &rv_up ) {
225        rv.dataind ( rv_up, v2v_down, v2v_up );
226        downsize = v2v_down.length();
227        upsize = v2v_up.length();
228}
229
230void datalink_m2e::set_connection ( const RV &rv, const RV &rvc, const RV &rv_up ) {
231        datalink::set_connection ( rv, rv_up );
232        condsize = rvc._dsize();
233        //establish v2c connection
234        rvc.dataind ( rv_up, v2c_lo, v2c_up );
235}
236
237vec datalink_m2e::get_cond ( const vec &val_up ) {
238        vec tmp ( condsize );
239        set_subvector ( tmp, v2c_lo, val_up ( v2c_up ) );
240        return tmp;
241}
242
243void datalink_m2e::pushup_cond ( vec &val_up, const vec &val, const vec &cond ) {
244        bdm_assert_debug ( downsize == val.length(), "Wrong val" );
245        bdm_assert_debug ( upsize == val_up.length(), "Wrong val_up" );
246        set_subvector ( val_up, v2v_up, val );
247        set_subvector ( val_up, v2c_up, cond );
248}
249
250std::ostream &operator<< ( std::ostream &os, const RV &rv ) {
251        int id;
252        for ( int i = 0; i < rv.len ; i++ ) {
253                id = rv.ids ( i );
254                os << id << "(" << RV::SIZES ( id ) << ")" <<  // id(size)=
255                "=" << RV::NAMES ( id )  << "_{"  << rv.times ( i ) << "}; "; //name_{time}
256        }
257        return os;
258}
259
260RV RV::expand_delayes() const {
261        RV rvt = this->remove_time(); //rv at t=0
262        RV tmp = rvt;
263        int td = mint();
264        for ( int i = -1; i >= td; i-- ) {
265                rvt.t_plus ( -1 );
266                tmp.add ( rvt ); //shift u1
267        }
268        return tmp;
269}
270
271str RV::tostr() const {
272        ivec idlist ( dsize );
273        ivec tmlist ( dsize );
274        int i;
275        int pos = 0;
276        for ( i = 0; i < len; i++ ) {
277                idlist.set_subvector ( pos, pos + size ( i ) - 1, ids ( i ) );
278                tmlist.set_subvector ( pos, pos + size ( i ) - 1, times ( i ) );
279                pos += size ( i );
280        }
281        return str ( idlist, tmlist );
282}
283
284ivec RV::dataind ( const RV &rv2 ) const {
285        ivec res ( 0 );
286        if ( rv2._dsize() > 0 ) {
287                str str2 = rv2.tostr();
288                ivec part;
289                int i;
290                for ( i = 0; i < len; i++ ) {
291                        part = itpp::find ( ( str2.ids == ids ( i ) ) & ( str2.times == times ( i ) ) );
292                        res = concat ( res, part );
293                }
294        }
295
296        //bdm_assert_debug ( res.length() == dsize, "this rv is not fully present in crv!" );
297        return res;
298
299}
300
301void RV::dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const {
302        //clean results
303        selfi.set_size ( 0 );
304        rv2i.set_size ( 0 );
305
306        // just in case any rv is empty
307        if ( ( len == 0 ) || ( rv2.length() == 0 ) ) {
308                return;
309        }
310
311        //find comon rv
312        ivec cids = itpp::find ( this->findself ( rv2 ) >= 0 );
313
314        // index of
315        if ( cids.length() > 0 ) {
316                str str1 = tostr();
317                str str2 = rv2.tostr();
318
319                ivec part1;
320                ivec part2;
321                int i, j;
322                // find common rv in strs
323                for ( j = 0; j < cids.length(); j++ ) {
324                        i = cids ( j );
325                        part1 = itpp::find ( ( str1.ids == ids ( i ) ) & ( str1.times == times ( i ) ) );
326                        part2 = itpp::find ( ( str2.ids == ids ( i ) ) & ( str2.times == times ( i ) ) );
327                        selfi = concat ( selfi, part1 );
328                        rv2i = concat ( rv2i, part2 );
329                }
330        }
331        bdm_assert_debug ( selfi.length() == rv2i.length(), "this should not happen!" );
332}
333
334RV RV::subt ( const RV &rv2 ) const {
335        ivec res = this->findself ( rv2 ); // nonzeros
336        ivec valid;
337        if ( dsize > 0 ) {
338                valid = itpp::find ( res == -1 );    //-1 => value not found => it remains
339        }
340        return ( *this ) ( valid ); //keep those that were not found in rv2
341}
342
343std::string RV::scalarname ( int scalat ) const {
344        bdm_assert ( scalat < dsize, "Wrong input index" );
345        int id = 0;
346        int scalid = 0;
347        while ( scalid + SIZES ( ids ( id ) ) <= scalat )  {
348                scalid += SIZES ( ids ( id ) );
349                id++;
350        };
351        //now id is the id of variable of interest
352        if ( size ( id ) == 1 )
353                return  NAMES ( ids ( id ) );
354        else
355                return  NAMES ( ids ( id ) ) + "_" + num2str ( scalat - scalid );
356
357}
358
359ivec RV::findself ( const RV &rv2 ) const {
360        int i, j;
361        ivec tmp = -ones_i ( len );
362        for ( i = 0; i < len; i++ ) {
363                for ( j = 0; j < rv2.length(); j++ ) {
364                        if ( ( ids ( i ) == rv2.ids ( j ) ) & ( times ( i ) == rv2.times ( j ) ) ) {
365                                tmp ( i ) = j;
366                                break;
367                        }
368                }
369        }
370        return tmp;
371}
372
373ivec RV::findself_ids ( const RV &rv2 ) const {
374        int i, j;
375        ivec tmp = -ones_i ( len );
376        for ( i = 0; i < len; i++ ) {
377                for ( j = 0; j < rv2.length(); j++ ) {
378                        if ( ( ids ( i ) == rv2.ids ( j ) ) ) {
379                                tmp ( i ) = j;
380                                break;
381                        }
382                }
383        }
384        return tmp;
385}
386
387void RV::from_setting ( const Setting &set ) {
388        Array<string> A;
389        UI::get ( A, set, "names" );
390
391        ivec szs;
392        if ( !UI::get ( szs, set, "sizes" ) )
393                szs = ones_i ( A.length() );
394
395        ivec tms;
396        if ( !UI::get ( tms, set, "times" ) )
397                tms = zeros_i ( A.length() );
398
399        init ( A, szs, tms );
400}
401
402void RV::to_setting ( Setting &set ) const {
403        Array<string> names ( len );
404        ivec sizes ( len );
405        for ( int i = 0; i < len; i++ ) {
406                names ( i ) = name ( i );
407                sizes ( i ) = size ( i );
408        }
409        UI::save ( names, set, "names" );
410        UI::save ( sizes, set, "sizes" );
411        UI::save ( times, set, "times" );
412}
413
414RV concat ( const RV &rv1, const RV &rv2 ) {
415        RV pom = rv1;
416        pom.add ( rv2 );
417        return pom;
418}
419
420RV get_composite_rv ( const Array<shared_ptr<pdf> > &pdfs,
421                      bool checkoverlap ) {
422        RV rv; //empty rv
423        bool rvaddok;
424        for ( int i = 0; i < pdfs.length(); i++ ) {
425                bdm_assert( pdfs(i)->isnamed(), "Can not extract RV from pdf no. " + num2str(i));
426                rvaddok = rv.add ( pdfs ( i )->_rv() ); //add rv to common rvs.
427                // If rvaddok==false, pdfs overlap => assert error.
428                bdm_assert_debug ( rvaddok || !checkoverlap, "mprod::mprod() input pdfs overlap in rv!" );
429        }
430
431        return rv;
432}
433
434void epdf::log_register ( logger &L, const string &prefix ) {
435        root::log_register ( L, prefix );
436
437        RV r;
438        if ( isnamed() ) {
439                r = _rv();
440        } else {
441                r = RV (  dimension() );
442        };
443
444        if ( log_level[logfull] ) {
445                // log full data
446                L.add_setting ( log_level, logfull, prefix );
447        } else {
448                // log only
449
450                if ( log_level[logmean] )
451                        L.add_vector ( log_level, logmean, r, prefix );                 
452                if ( log_level[loglbound] )
453                        L.add_vector ( log_level, loglbound, r, prefix );
454                if ( log_level[logubound] )
455                        L.add_vector ( log_level, logubound, r, prefix );
456        }
457}
458
459void epdf::log_write() const {
460        if ( log_level[logfull] ) {
461                log_level.store( logfull, this );
462        } else {
463                if ( log_level[logmean] ) {
464                        log_level.store( logmean, mean() );
465                }
466                if ( log_level[loglbound] || log_level[logubound] ) {
467                                vec lb;
468                                vec ub;
469                                qbounds ( lb, ub );
470                                if (log_level[loglbound])
471                                        log_level.store( loglbound, lb );
472                                if (log_level[logubound])
473                                        log_level.store( logubound, ub );
474                        }
475                }
476        }
477
478
479void datalink_buffered::set_connection ( const RV &rv, const RV &rv_up ) {
480        // create link between up and down
481        datalink_part::set_connection ( rv, rv_up); // only non-delayed version
482
483        RV needed_from_hist = rv.subt(rv_up); //rv_up already copied by v2v
484       
485        // we can store only what we get in rv_up - everything else is removed
486        ivec valid_ids = needed_from_hist.findself_ids ( rv_up ); // return on which position the required id is in rv_up
487        RV rv_hist = needed_from_hist.subselect ( find ( valid_ids >= 0 ) ); // select only rvs that are in rv_up, ie ind>0
488        RV rv_hist0 = rv_hist.remove_time(); // these RVs will form history at time =0
489        // now we need to know what is needed from Up
490        rv_hist = rv_hist.expand_delayes(); // full regressor - including time 0
491        Hrv = rv_hist.subt ( rv_hist0 );   // remove time 0
492        history = zeros ( Hrv._dsize() );
493
494        // decide if we need to copy val to history
495        if ( Hrv._dsize() > 0 ) {
496                v2h_up = rv_hist0.dataind ( rv_up ); // indices of elements of rv_up to be copied
497        } // else v2h_up is empty
498
499        Hrv.dataind ( rv, h2v_hist, h2v_down );
500
501        downsize = v2v_down.length() + h2v_down.length();
502        upsize = v2v_up.length();
503}
504
505void datalink_buffered::set_history ( const RV& rv1, const vec &hist0 ) {
506        bdm_assert ( rv1._dsize() == hist0.length(), "hist is not compatible with given rv1" );
507        ivec ind_H;
508        ivec ind_h0;
509        Hrv.dataind ( rv1, ind_H, ind_h0 ); // find indices of rv in
510        set_subvector ( history, ind_H, hist0 ( ind_h0 ) ); // copy given hist to appropriate places
511}
512
513void DS::log_register ( logger &L,  const string &prefix ) {
514        bdm_assert ( dtsize == Drv._dsize(), "invalid DS: dtsize (" + num2str ( dtsize ) + ") different from Drv " + num2str ( Drv._dsize() ) );
515        //bdm_assert ( utsize == Urv._dsize(), "invalid DS: utsize (" + num2str ( utsize ) + ") different from Urv " + num2str ( Urv._dsize() ) );
516
517        root::log_register ( L, prefix );
518
519
520        if ( log_level[logdt] )
521                L.add_vector ( log_level, logdt, Drv, prefix ); 
522        if ( log_level[logut] )
523                L.add_vector ( log_level, logut, Urv, prefix );
524}
525
526void DS::log_write ( ) const {
527        if( log_level[logdt] ) {
528                vec tmp ( Drv._dsize());
529                getdata ( tmp );
530                // d is first in getdata
531                log_level.store( logdt, tmp );
532        }
533        if( log_level[logut] ) { 
534                // NOT_IMPLEMENTED
535        }
536}
537
538
539void DS::from_setting ( const Setting &set ) {
540        RV rv;
541        if( UI::get( rv, set, "drv", UI::optional ) )
542                set_drv ( rv, RV() );
543}
544
545void DS::validate() {
546        if( !Drv._dsize() && dtsize > 0 )
547        {
548                RV *r = new RV();
549                for ( int i = 0; i < dtsize; i++ ) {
550                        r->add ( RV ( "ch" + num2str ( i ), 1, 0 ) );
551                }
552                set_drv ( *r, RV() );
553                delete r;
554        }
555}
556
557void BM::log_register ( logger &L, const string &prefix ) {
558        root::log_register ( L, prefix );
559
560        if ( log_level[logevidence] )
561                L.add_vector ( log_level, logevidence, RV ( 1 ), prefix );     
562
563        if (log_level[logbounds]){
564                prior().log_level[epdf::loglbound]=true;
565                prior().log_level[epdf::logubound]=true;
566        }
567        if (log_level[logfull]){
568                prior().log_level[epdf::logfull]=true;
569        }
570        const_cast<epdf&> ( posterior() ).log_register ( L, prefix + L.separator + "apost" );
571}
572
573void BM::log_write ( ) const {
574        posterior().log_write();
575
576        if( log_level[logevidence] ) 
577                log_level.store( logevidence, ll );
578}
579
580void BM::bayes_batch ( const mat &Data, const vec &cond ) {
581        for ( int t = 0; t < Data.cols(); t++ ) {
582                bayes ( Data.get_col ( t ), cond );
583        }
584}
585
586void BM::bayes_batch ( const mat &Data, const mat &Cond ) {
587        for ( int t = 0; t < Data.cols(); t++ ) {
588                bayes ( Data.get_col ( t ), Cond.get_col ( t ) );
589        }
590}
591
592}
Note: See TracBrowser for help on using the browser.