1 // 2 // sampling.cpp 3 // cufflinks 4 // 5 // Created by Cole Trapnell on 12/19/11. 6 // Copyright 2011 __MyCompanyName__. All rights reserved. 7 // 8 9 #include "sampling.h" 10 #include <limits> 11 12 using namespace std; 13 generate_importance_samples(multinormal_generator<double> & generator,std::vector<Eigen::VectorXd> & samples,int num_samples,bool no_zeros)14void generate_importance_samples(multinormal_generator<double>& generator, 15 std::vector<Eigen::VectorXd>& samples, 16 int num_samples, 17 bool no_zeros) 18 { 19 for (int i = 0; i < num_samples; ++i) 20 { 21 // TODO: we should switch the multinormal generator over to Eigen for 22 // consistency as part of the push to drop uBLAS. 23 boost::numeric::ublas::vector<double> r = generator.next_rand(); 24 25 Eigen::VectorXd scaled_sample = Eigen::VectorXd::Zero(r.size()); 26 27 for (int j = 0; j < scaled_sample.size(); ++j) 28 { 29 scaled_sample(j) = r(j); 30 31 if (scaled_sample(j) < 0) 32 { 33 scaled_sample(j) = 0.0; 34 //scaled_sample(j) = -scaled_sample(j); 35 } 36 } 37 38 double m = scaled_sample.sum(); 39 if (m && !isnan(m)) 40 { 41 for (int j = 0; j < scaled_sample.size(); ++j) 42 { 43 scaled_sample(j) = scaled_sample(j) / m; 44 } 45 if (no_zeros) 46 { 47 bool has_zero = false; 48 for (int j = 0; j < scaled_sample.size(); ++j) 49 { 50 if (scaled_sample[j] == 0) 51 { 52 has_zero = true; 53 break; 54 } 55 } 56 57 if (has_zero) 58 continue; 59 } 60 samples.push_back(scaled_sample); 61 } 62 else 63 { 64 samples.push_back(Eigen::VectorXd::Zero(scaled_sample.size())); 65 //cerr << r << endl; 66 //cerr << scaled_sample << endl; 67 } 68 } 69 } 70