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

Revision 598, 9.2 kB (checked in by smidl, 15 years ago)

new buffered datalink ( #32 ) and new datasources - all with minor trivial tests

  • 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;
8RVmap RV_MAP;
9Array<string> RV_NAMES ( RV_BUFFER_STEP );
10ivec RV_SIZES ( RV_BUFFER_STEP );
11
12RV RV0 = RV();
13
14void RV::clear_all() {
15        RV_MAP.clear();
16        RV_SIZES.clear();
17        RV_NAMES = Array<string> ( RV_BUFFER_STEP );
18}
19
20int RV::init ( const string &name, int size ) {
21        //Refer
22        int id;
23        RVmap::const_iterator iter = RV_MAP.find ( name );
24        if ( iter == RV_MAP.end() ) {
25                id = RV_MAP.size() + 1;
26                //debug
27                /*              {
28                                        cout << endl;
29                                        RVmap::const_iterator iter = RV_MAP.begin();
30                                        for(RVmap::const_iterator iter=RV_MAP.begin(); iter!=RV_MAP.end(); iter++){
31                                                cout << "key: " << iter->first << " val: " << iter->second <<endl;
32                                        }
33                                }*/
34
35                RV_MAP.insert ( make_pair ( name, id ) ); //add new rv
36                if ( id >= RV_NAMES.length() ) {
37                        RV_NAMES.set_length ( id + RV_BUFFER_STEP, true );
38                        RV_SIZES.set_length ( id + RV_BUFFER_STEP, true );
39                }
40                RV_NAMES ( id ) = name;
41                RV_SIZES ( id ) = size;
42        } else {
43                id = iter->second;
44                bdm_assert_debug ( RV_SIZES ( id ) == size, "RV " + name + " of different size already exists" );
45        }
46        return id;
47};
48
49int RV::countsize() const {
50        int tmp = 0;
51        for ( int i = 0; i < len; i++ ) {
52                tmp += RV_SIZES ( ids ( i ) );
53        }
54        return tmp;
55}
56
57ivec RV::cumsizes() const {
58        ivec szs ( len );
59        int tmp = 0;
60        for ( int i = 0; i < len; i++ ) {
61                tmp += RV_SIZES ( ids ( i ) );
62                szs ( i ) = tmp;
63        }
64        return szs;
65}
66
67void RV::init ( const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times ) {
68        len = in_names.length();
69        bdm_assert_debug ( in_names.length() == in_times.length(), "check \"times\" " );
70        bdm_assert_debug ( in_names.length() == in_sizes.length(), "check \"sizes\" " );
71
72        times.set_length ( len );
73        ids.set_length ( len );
74        int id;
75        for ( int i = 0; i < len; i++ ) {
76                id = init ( in_names ( i ), in_sizes ( i ) );
77                ids ( i ) = id;
78        }
79        times = in_times;
80        dsize = countsize();
81}
82
83RV::RV ( string name, int sz, int tm ) {
84        Array<string> A ( 1 );
85        A ( 0 ) = name;
86        init ( A, vec_1 ( sz ), vec_1 ( tm ) );
87}
88
89bool RV::add ( const RV &rv2 ) {
90        // TODO
91        if ( rv2.len > 0 ) { //rv2 is nonempty
92                ivec ind = rv2.findself ( *this ); //should be -1 all the time
93                ivec index = itpp::find ( ind == -1 );
94
95                if ( index.length() < rv2.len ) { //conflict
96                        ids = concat ( ids, rv2.ids ( index ) );
97                        times = concat ( times, rv2.times ( index ) );
98                } else {
99                        ids = concat ( ids, rv2.ids );
100                        times = concat ( times, rv2.times );
101                }
102                len = ids.length();
103                dsize = countsize();
104                return ( index.length() == rv2.len ); //conflict or not
105        } else { //rv2 is empty
106                return true; // no conflict
107        }
108};
109
110RV RV::subselect ( const ivec &ind ) const {
111        RV ret;
112        ret.ids = ids ( ind );
113        ret.times = times ( ind );
114        ret.len = ind.length();
115        ret.dsize = ret.countsize();
116        return ret;
117}
118
119RV RV::operator() ( int di1, int di2 ) const {
120        ivec sz = cumsizes();
121        int i1 = 0;
122        while ( sz ( i1 ) < di1 ) i1++;
123        int i2 = i1;
124        while ( sz ( i2 ) < di2 ) i2++;
125        return subselect ( linspace ( i1, i2 ) );
126}
127
128void RV::t ( int delta ) {
129        times += delta;
130}
131
132bool RV::equal ( const RV &rv2 ) const {
133        return ( ids == rv2.ids ) && ( times == rv2.times );
134}
135
136shared_ptr<mpdf> epdf::condition (const RV &rv) const {
137        bdm_warning ("Not implemented");
138        return shared_ptr<mpdf>();
139}
140
141shared_ptr<epdf> epdf::marginal (const RV &rv) const {
142        bdm_warning ("Not implemented");
143        return shared_ptr<epdf>();
144}
145
146mat epdf::sample_m ( int N ) const {
147        mat X = zeros ( dim, N );
148        for ( int i = 0; i < N; i++ ) X.set_col ( i, this->sample() );
149        return X;
150}
151
152vec epdf::evallog_m (const mat &Val) const {
153        vec x (Val.cols());
154        for (int i = 0; i < Val.cols(); i++) {
155                x (i) = evallog ( Val.get_col (i) );
156        }
157
158        return x;
159}
160
161vec epdf::evallog_m (const Array<vec> &Avec) const {
162        vec x (Avec.size());
163        for (int i = 0; i < Avec.size(); i++) {
164                x (i) = evallog ( Avec (i) );
165        }
166
167        return x;
168}
169
170mat mpdf::samplecond_m (const vec &cond, int N) {
171        mat M( dimension(), N );
172        for ( int i = 0; i < N; i++ ) {
173                M.set_col ( i, samplecond( cond ) );
174        }
175
176        return M;
177}
178
179void mpdf::from_setting ( const Setting &set ) {
180        shared_ptr<RV> r = UI::build<RV> ( set, "rv", UI::optional );
181        if ( r ) {
182                set_rv ( *r );
183        }
184
185        r = UI::build<RV> ( set, "rvc", UI::optional );
186        if ( r ) {
187                set_rvc ( *r );
188        }
189}
190
191void datalink::set_connection (const RV &rv, const RV &rv_up) {
192        downsize = rv._dsize();
193          upsize = rv_up._dsize();
194          v2v_up = rv.dataind (rv_up);
195          bdm_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up");
196}
197
198void datalink::set_connection (int ds, int us, const ivec &upind) {
199        downsize = ds;
200        upsize = us;
201        v2v_up = upind;
202    bdm_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up");
203}
204
205void datalink_part::set_connection (const RV &rv, const RV &rv_up) {
206         rv.dataind (rv_up, v2v_down, v2v_up);
207         downsize = v2v_down.length();
208         upsize = v2v_up.length();
209}
210
211void datalink_m2e::set_connection (const RV &rv, const RV &rvc, const RV &rv_up) {
212        datalink::set_connection (rv, rv_up);
213        condsize = rvc._dsize();
214        //establish v2c connection
215        rvc.dataind (rv_up, v2c_lo, v2c_up);
216}
217
218vec datalink_m2e::get_cond (const vec &val_up) {
219        vec tmp (condsize);
220        set_subvector (tmp, v2c_lo, val_up (v2c_up));
221        return tmp;
222}
223
224void datalink_m2e::pushup_cond (vec &val_up, const vec &val, const vec &cond) {
225        bdm_assert_debug (downsize == val.length(), "Wrong val");
226        bdm_assert_debug (upsize == val_up.length(), "Wrong val_up");
227        set_subvector (val_up, v2v_up, val);
228        set_subvector (val_up, v2c_up, cond);
229}
230
231std::ostream &operator<< ( std::ostream &os, const RV &rv ) {
232        int id;
233        for ( int i = 0; i < rv.len ; i++ ) {
234                id = rv.ids ( i );
235                os << id << "(" << RV_SIZES ( id ) << ")" <<  // id(size)=
236                "=" << RV_NAMES ( id )  << "_{"  << rv.times ( i ) << "}; "; //name_{time}
237        }
238        return os;
239}
240
241str RV::tostr() const {
242        ivec idlist ( dsize );
243        ivec tmlist ( dsize );
244        int i;
245        int pos = 0;
246        for ( i = 0; i < len; i++ ) {
247                idlist.set_subvector ( pos, pos + size ( i ) - 1, ids ( i ) );
248                tmlist.set_subvector ( pos, pos + size ( i ) - 1, times ( i ) );
249                pos += size ( i );
250        }
251        return str ( idlist, tmlist );
252}
253
254ivec RV::dataind ( const RV &rv2 ) const {
255        ivec res ( 0 );
256        if ( rv2._dsize() > 0 ) {
257                str str2 = rv2.tostr();
258                ivec part;
259                int i;
260                for ( i = 0; i < len; i++ ) {
261                        part = itpp::find ( ( str2.ids == ids ( i ) ) & ( str2.times == times ( i ) ) );
262                        res = concat ( res, part );
263                }
264        }
265
266        bdm_assert_debug ( res.length() == dsize, "this rv is not fully present in crv!" );
267        return res;
268
269}
270
271void RV::dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const {
272        //clean results
273        selfi.set_size ( 0 );
274        rv2i.set_size ( 0 );
275
276        // just in case any rv is empty
277        if ( ( len == 0 ) || ( rv2.length() == 0 ) ) {
278                return;
279        }
280
281        //find comon rv
282        ivec cids = itpp::find ( this->findself ( rv2 ) >= 0 );
283
284        // index of
285        if ( cids.length() > 0 ) {
286                str str1 = tostr();
287                str str2 = rv2.tostr();
288
289                ivec part1;
290                ivec part2;
291                int i, j;
292                // find common rv in strs
293                for ( j = 0; j < cids.length(); j++ ) {
294                        i = cids ( j );
295                        part1 = itpp::find ( ( str1.ids == ids ( i ) ) & ( str1.times == times ( i ) ) );
296                        part2 = itpp::find ( ( str2.ids == ids ( i ) ) & ( str2.times == times ( i ) ) );
297                        selfi = concat ( selfi, part1 );
298                        rv2i = concat ( rv2i, part2 );
299                }
300        }
301        bdm_assert_debug ( selfi.length() == rv2i.length(), "this should not happen!" );
302}
303
304RV RV::subt ( const RV &rv2 ) const {
305        ivec res = this->findself ( rv2 ); // nonzeros
306        ivec valid;
307        if ( dsize > 0 ) {
308                valid = itpp::find ( res == -1 );    //-1 => value not found => it remains
309        }
310        return ( *this ) ( valid ); //keep those that were not found in rv2
311}
312
313ivec RV::findself ( const RV &rv2 ) const {
314        int i, j;
315        ivec tmp = -ones_i ( len );
316        for ( i = 0; i < len; i++ ) {
317                for ( j = 0; j < rv2.length(); j++ ) {
318                        if ( ( ids ( i ) == rv2.ids ( j ) ) & ( times ( i ) == rv2.times ( j ) ) ) {
319                                tmp ( i ) = j;
320                                break;
321                        }
322                }
323        }
324        return tmp;
325}
326
327ivec RV::findself_ids ( const RV &rv2 ) const {
328        int i, j;
329        ivec tmp = -ones_i ( len );
330        for ( i = 0; i < len; i++ ) {
331                for ( j = 0; j < rv2.length(); j++ ) {
332                        if ( ( ids ( i ) == rv2.ids ( j ) ) ) {
333                                tmp ( i ) = j;
334                                break;
335                        }
336                }
337        }
338        return tmp;
339}
340
341void RV::from_setting ( const Setting &set ) {
342        Array<string> A;
343        UI::get ( A, set, "names" );
344
345        ivec szs;
346        if ( !UI::get ( szs, set, "sizes" ) )
347                szs = ones_i ( A.length() );
348
349        ivec tms;
350        if ( !UI::get ( tms, set, "times" ) )
351                tms = zeros_i ( A.length() );
352
353        // TODO tady se bude plnit primo do jeho promennych, a pak se zavola validacnni metoda, takze cele prepsat, ano?
354        init ( A, szs, tms );
355}
356
357RV concat ( const RV &rv1, const RV &rv2 ) {
358        RV pom = rv1;
359        pom.add ( rv2 );
360        return pom;
361}
362
363void mepdf::from_setting ( const Setting &set ) {
364        shared_ptr<epdf> e ( UI::build<epdf> ( set, "epdf", UI::compulsory ) );
365        iepdf = e;
366        set_ep(iepdf.get());
367}
368
369RV get_composite_rv ( const Array<shared_ptr<mpdf> > &mpdfs,
370                      bool checkoverlap ) {
371        RV rv; //empty rv
372        bool rvaddok;
373        for ( int i = 0; i < mpdfs.length(); i++ ) {
374                rvaddok = rv.add ( mpdfs ( i )->_rv() ); //add rv to common rvs.
375                // If rvaddok==false, mpdfs overlap => assert error.
376                bdm_assert_debug ( rvaddok || !checkoverlap, "mprod::mprod() input mpdfs overlap in rv!" );
377        }
378
379        return rv;
380}
381
382void BM::bayesB ( const mat &Data ) {
383        for ( int t = 0; t < Data.cols(); t++ ) {
384                bayes ( Data.get_col ( t ) );
385        }
386}
387}
Note: See TracBrowser for help on using the browser.