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

#include "bdmbase.h"

//! Space of basic BDM structures
namespace bdm {

const int RV_BUFFER_STEP=1;
RVmap RV_MAP;
Array<string> RV_NAMES ( RV_BUFFER_STEP );
ivec RV_SIZES ( RV_BUFFER_STEP );

RV RV0=RV();

void RV::clear_all()
{
    RV_MAP.clear();
    RV_SIZES.clear();
    RV_NAMES = Array<string>(RV_BUFFER_STEP);
}
  
int RV::init ( const string &name, int size ) {
        //Refer
        int id;
        RVmap::const_iterator iter = RV_MAP.find ( name );
        if ( iter == RV_MAP.end() ) {
                id=RV_MAP.size() +1;
                //debug
/*              {
                        cout << endl;
                        RVmap::const_iterator iter = RV_MAP.begin();
                        for(RVmap::const_iterator iter=RV_MAP.begin(); iter!=RV_MAP.end(); iter++){
                                cout << "key: " << iter->first << " val: " << iter->second <<endl;
                        }
                }*/
                
                RV_MAP.insert ( make_pair ( name,id ) ); //add new rv
                if ( id>=RV_NAMES.length() ) {
                        RV_NAMES.set_length ( id+RV_BUFFER_STEP,true );
                        RV_SIZES.set_length ( id+RV_BUFFER_STEP,true );
                }
                RV_NAMES ( id ) =name;
                RV_SIZES ( id ) =size;
        }
        else {
                id = iter->second;
                it_assert(RV_SIZES(id)==size,"RV "+name+" of different size already exists");
        }
        return id;
};

int RV::countsize() const{ int tmp=0; for ( int i=0;i<len;i++ ) {tmp+=RV_SIZES ( ids ( i ) );} return tmp;}

ivec RV::cumsizes() const {
        ivec szs ( len );
        int tmp=0;
        for ( int i=0;i<len;i++ ) {
                tmp+=RV_SIZES ( ids ( i ) );
                szs ( i ) = tmp;
        }
        return szs;
}

void RV::init(const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times)
{
        len = in_names.length();
        it_assert_debug ( in_names.length() ==in_times.length(), "check \"times\" " );
        it_assert_debug ( in_names.length() ==in_sizes.length(), "check \"sizes\" " );

        times.set_length ( len );
        ids.set_length ( len );
        int id;
        for ( int i=0; i<len; i++ ) {
                id = init ( in_names ( i ), in_sizes ( i ) );
                ids ( i ) = id;
        }
        times = in_times;
        dsize = countsize();
}

RV::RV ( string name, int sz, int tm ) {
        Array<string> A ( 1 ); A ( 0 ) =name;
        init ( A,vec_1 ( sz ),vec_1 ( tm ) );
}

bool RV::add ( const RV &rv2 ) {
        // TODO
        if ( rv2.len>0 ) { //rv2 is nonempty
                ivec ind = rv2.findself ( *this ); //should be -1 all the time
                ivec index = itpp::find ( ind==-1 );

                if ( index.length() < rv2.len ) { //conflict
                        ids = concat ( ids, rv2.ids ( index ) );
                        times = concat ( times, rv2.times ( index ) );
                }
                else {
                        ids = concat ( ids, rv2.ids );
                        times = concat ( times, rv2.times );
                }
                len = ids.length();
                dsize = countsize();
                return ( index.length() ==rv2.len ); //conflict or not
        }
        else { //rv2 is empty
                return true; // no conflict
        }
};

RV RV::subselect ( const ivec &ind ) const {
        RV ret;
        ret.ids = ids ( ind );
        ret.times= times ( ind );
        ret.len = ind.length();
        ret.dsize=ret.countsize();
        return ret;
}

RV RV::operator()(int di1, int di2) const
{
    ivec sz = cumsizes();
    int i1 = 0;
    while (sz(i1) < di1) i1++;
    int i2 = i1;
    while (sz(i2) < di2) i2++;
    return subselect(linspace(i1, i2));
}

void RV::t ( int delta ) { times += delta;}

bool RV::equal ( const RV &rv2 ) const {
        return ( ids == rv2.ids ) && ( times == rv2.times );
}

mat epdf::sample_m ( int N ) const {
        mat X = zeros ( dim, N );
        for ( int i = 0;i < N;i++ ) X.set_col ( i, this->sample() );
        return X;
}

vec mpdf::samplecond(const vec &cond) {
    condition(cond);
    vec temp = shep->sample();
    return temp;
}

mat mpdf::samplecond_m(const vec &cond, int N) {
    condition(cond);
    mat temp(shep->dimension(), N);
    vec smp(shep->dimension());
    for (int i = 0; i < N; i++) {
        smp = shep->sample();
        temp.set_col(i, smp);
    }

    return temp;
}

double mpdf::evallogcond(const vec &dt, const vec &cond) {
    double tmp;
    condition(cond);
    tmp = shep->evallog(dt);
    // it_assert_debug(std::isfinite(tmp), "Infinite value");
    return tmp;
}

vec mpdf::evallogcond_m(const mat &Dt, const vec &cond) {
    condition(cond);
    return shep->evallog_m(Dt);
}

vec mpdf::evallogcond_m(const Array<vec> &Dt, const vec &cond) {
    condition(cond);
    return shep->evallog_m(Dt);
}

void mpdf::from_setting(const Setting &set){
    if (set.exists("rv")) {
        RV *r = UI::build<RV>(set, "rv");
        set_rv(*r);
        delete r;
    }

    if (set.exists("rvc")) {
        RV *r = UI::build<RV>(set, "rvc");
        set_rvc(*r); 
        delete r;
    }
}

std::ostream &operator<< ( std::ostream &os, const RV &rv ) {
        int id;
        for ( int i = 0; i < rv.len ;i++ ) {
                id=rv.ids ( i );
                os << id << "(" << RV_SIZES ( id ) << ")" <<  // id(size)=
                "=" << RV_NAMES ( id )  << "_{"  << rv.times ( i ) << "}; "; //name_{time}
        }
        return os;
}

str RV::tostr() const {
        ivec idlist ( dsize );
        ivec tmlist ( dsize );
        int i;
        int pos = 0;
        for ( i = 0;i < len;i++ ) {
                idlist.set_subvector ( pos, pos + size ( i ) - 1, ids ( i ) );
                tmlist.set_subvector ( pos, pos + size ( i ) - 1, times ( i ) );
                pos += size ( i );
        }
        return str ( idlist, tmlist );
}

ivec RV::dataind ( const RV &rv2 ) const {
        ivec res ( 0 );
        if ( rv2._dsize() >0 ) {
                str str2 = rv2.tostr();
                ivec part;
                int i;
                for ( i = 0;i < len;i++ ) {
                        part = itpp::find ( ( str2.ids == ids ( i ) ) & ( str2.times == times ( i ) ) );
                        res = concat ( res, part );
                }
        }
        it_assert_debug ( res.length() ==dsize,"this rv is not fully present in crv!" );
        return res;

}

void RV::dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const {
        //clean results
        selfi.set_size ( 0 );
        rv2i.set_size ( 0 );

        // just in case any rv is empty
        if ( ( len==0 ) || ( rv2.length() ==0 ) ) {return;}

        //find comon rv
        ivec cids=itpp::find ( this->findself ( rv2 ) >=0 );

        // index of
        if ( cids.length() >0 ) {
                str str1 = tostr();
                str str2 = rv2.tostr();

                ivec part1;
                ivec part2;
                int i,j;
                // find common rv in strs
                for ( j=0; j < cids.length();j++ ) {
                        i = cids ( j );
                        part1 = itpp::find ( ( str1.ids == ids ( i ) ) & ( str1.times == times ( i ) ) );
                        part2 = itpp::find ( ( str2.ids == ids ( i ) ) & ( str2.times == times ( i ) ) );
                        selfi = concat ( selfi, part1 );
                        rv2i = concat ( rv2i, part2 );
                }
        }
        it_assert_debug ( selfi.length() == rv2i.length(),"this should not happen!" );
}

RV RV::subt ( const RV &rv2 ) const {
        ivec res = this->findself ( rv2 ); // nonzeros
        ivec valid;
        if ( dsize>0 ) {valid= itpp::find ( res == -1 );} //-1 => value not found => it remains
        return ( *this ) ( valid ); //keep those that were not found in rv2
}

ivec RV::findself ( const RV &rv2 ) const {
        int i, j;
        ivec tmp = -ones_i ( len );
        for ( i = 0;i < len;i++ ) {
                for ( j = 0;j < rv2.length();j++ ) {
                        if ( ( ids ( i ) == rv2.ids ( j ) ) & ( times ( i ) == rv2.times ( j ) ) ) {
                                tmp ( i ) = j;
                                break;
                        }
                }
        }
        return tmp;
}

void RV::from_setting( const Setting &set ) 
{       
        Array<string> A;
        if( set.exists("names"))
                UI::get( A, set, "names" );
        else
                A.set_length(0);
        
        ivec szs;
        if( set.exists("sizes"))
                UI::get(szs,set,"sizes");
        else
                szs = ones_i(A.length());
        
        ivec tms;
        if( set.exists( "times") )
                UI::get(tms,set,"times");
        else
                tms = zeros_i(A.length());
        
        // TODO tady se bude plnit primo do jeho promennych, a pak se zavola validacnni metoda, takze cele prepsat, ano?
        init( A, szs, tms );
}

/*void RV::to_setting( Setting &set ) const
{       
        Transport::to_setting( set );

        Setting &kilometers_setting = set.add("kilometers", Setting::TypeInt );
        kilometers_setting = kilometers;

        UI::save( passengers, set, "passengers" );
}*/

RV concat ( const RV &rv1, const RV &rv2 ) {
        RV pom = rv1;
        pom.add ( rv2 );
        return pom;
}

void mepdf::condition(const vec &cond) {
}

void mepdf::from_setting(const Setting &set) {
    shared_ptr<epdf> e(UI::build<epdf>(set, "epdf"));
    set_ep(e);
}

RV compositepdf::getrv ( bool checkoverlap ) {
        RV rv; //empty rv
        bool rvaddok;
        for ( int i = 0;i < mpdfs.length();i++ ) {
                rvaddok=rv.add ( mpdfs ( i )->_rv() ); //add rv to common rvs.
                // If rvaddok==false, mpdfs overlap => assert error.
                it_assert_debug ( rvaddok|| ( !checkoverlap ),"mprod::mprod() input mpdfs overlap in rv!" );
        };
        return rv;
}

void compositepdf::setrvc ( const RV &rv, RV &rvc ) {
        for ( int i = 0;i < mpdfs.length();i++ ) {
                RV rvx = mpdfs ( i )->_rvc().subt ( rv );
                rvc.add ( rvx ); //add rv to common rvc
        };
}

void BM::bayesB ( const mat &Data ) {
        for ( int t=0;t<Data.cols();t++ ) {bayes ( Data.get_col ( t ) );}
}
}