root/bdm/stat/libBM.cpp @ 263

#include "libBM.h"
#include "../itpp_ext.h"

//! Space of basic BDM structures
namespace bdm {

using std::cout;

static int RVcounter=0;

RV RV0=RV();

void RV::init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times ) {
        //Refer
        int i;
        len = in_ids.length();
        //PRUDENT_MODE
        // All vetors should be of same length
        if ( ( len != in_names.length() ) || \
                ( len != in_sizes.length() ) || \
                ( len != in_times.length() ) ) {
                it_error ( "RV::RV inconsistent length of input vectors." );
        }

        ids = in_ids;
        names = in_names;
        sizes = in_sizes;
        times = in_times;
        tsize = 0;
        for ( i = 0;i < len;i++ ) {tsize += sizes ( i );}
};

RV::RV ( Array<std::string> in_names, ivec in_sizes, ivec in_times ) {
        int len = in_names.length();
        init ( linspace ( RVcounter+1, RVcounter+len ), in_names, in_sizes, in_times );
        RVcounter+=len;
}

RV::RV ( Array<std::string> in_names, ivec in_sizes ) {
        int len = in_names.length();
        init ( linspace ( RVcounter+1, RVcounter+len ), in_names, in_sizes, zeros_i ( in_names.length() ) );
        RVcounter+=len;
}

RV::RV ( Array<std::string> in_names ) {
        int len = in_names.length();
        init ( linspace ( RVcounter+1, RVcounter+len ), in_names, ones_i ( in_names.length() ), zeros_i ( in_names.length() ) );
        RVcounter+=len;
}

RV::RV() : tsize ( 0 ), len ( 0 ), ids ( 0 ),  sizes ( 0 ), times ( 0 ),names ( 0 ) {};

RV::RV(string name, int id){
        it_assert(id>=RVcounter,"RV::This id is already taken");
        RVcounter=id;
        Array<string> A(1); A(0)=name;
        init(vec_1(id),A,vec_1(1),vec_1(0));
}

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 ) );
                        sizes = concat ( sizes, rv2.sizes ( index ) );
                        times = concat ( times, rv2.times ( index ) );
                        names = concat ( names, rv2.names ( to_Arr ( index ) ) );
                }
                else {
                        ids = concat ( ids, rv2.ids );
                        sizes = concat ( sizes, rv2.sizes );
                        times = concat ( times, rv2.times );
                        names = concat ( names, rv2.names );
                }
                tsize = sum ( sizes );
                len = ids.length();
                return ( index.length() ==rv2.len );
        }
        else { //rv2 is empty
                return true;
        }
//      return *this;
};

// RV::RV ( ivec in_ids ) {
//
//      len = in_ids.length();
//      Array<std::string> A ( len );
//      std::string rvstr = "rv";
//
//      for ( int i = 0; i < len;i++ ) {
//              A ( i ) = rvstr + to_str ( i );
//      }
//
//      init ( in_ids, A, ones_i ( len ), zeros_i ( len ) );
// }

RV RV::subselect (const ivec &ind ) const {
        RV ret;
        ret.init ( ids ( ind ), names ( to_Arr ( ind ) ), sizes ( ind ), times ( ind ) );
        return ret;
}

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

RV RV::operator() (const ivec &ind ) const {
        RV ret;
        if ( ind.length() >0 ) {
                ret.init ( ids ( ind ), names ( to_Arr ( ind ) ), sizes ( ind ), times ( ind ) );
        }
        return ret;
}

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

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


std::ostream &operator<< ( std::ostream &os, const RV &rv ) {

        for ( int i = 0; i < rv.len ;i++ ) {
                os << rv.ids ( i ) << "(" << rv.sizes ( i ) << ")" <<  // id(size)=
                "=" << rv.names ( i )  << "_{"  << rv.times ( i ) << "}; "; //name_{time}
        }
        return os;
}

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

ivec RV::dataind (const RV &rv2 ) const {
        ivec res ( 0 );
        if ( rv2.count()>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()==tsize,"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 (tsize>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;
}

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

void RV::newids(){ids=linspace ( RVcounter+1, RVcounter+len ),RVcounter+=len;}

RV compositepdf::getrv(bool checkoverlap){
        RV rv; //empty rv
        bool rvaddok;
        for (int i = 0;i < n;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 < n;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));}
}
}