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

Revision 1064, 16.4 kB (checked in by mido, 14 years ago)

astyle applied all over the library

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