| ... | ... |
@@ -17,7 +17,9 @@ void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* max |
| 17 | 17 |
|
| 18 | 18 |
//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__; |
| 19 | 19 |
// Parallelization settings |
| 20 |
+ #ifdef _OPENMP |
|
| 20 | 21 |
omp_set_num_threads(*num_threads); |
| 22 |
+ #endif |
|
| 21 | 23 |
|
| 22 | 24 |
// Print some information |
| 23 | 25 |
//FILE_LOG(logINFO) << "number of states = " << *N; |
| ... | ... |
@@ -267,7 +269,9 @@ void multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, doubl |
| 267 | 269 |
|
| 268 | 270 |
//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__; |
| 269 | 271 |
// Parallelization settings |
| 272 |
+ #ifdef _OPENMP |
|
| 270 | 273 |
omp_set_num_threads(*num_threads); |
| 274 |
+ #endif |
|
| 271 | 275 |
|
| 272 | 276 |
// Print some information |
| 273 | 277 |
//FILE_LOG(logINFO) << "number of states = " << *N; |
| ... | ... |
@@ -675,119 +675,6 @@ void NegativeBinomial::update(double* weights) |
| 675 | 675 |
|
| 676 | 676 |
} |
| 677 | 677 |
|
| 678 |
-// void NegativeBinomial::update(double* weight) |
|
| 679 |
-// {
|
|
| 680 |
-// //FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__; |
|
| 681 |
-// double eps = 1e-4, kmax; |
|
| 682 |
-// double numerator, denominator, rhere, dr, Fr, dFrdr, DigammaR, DigammaRplusDR; |
|
| 683 |
-// // Update p |
|
| 684 |
-// numerator=denominator=0.0; |
|
| 685 |
-// // clock_t time, dtime; |
|
| 686 |
-// // time = clock(); |
|
| 687 |
-// for (int t=0; t<this->T; t++) |
|
| 688 |
-// {
|
|
| 689 |
-// numerator+=weight[t]*this->size; |
|
| 690 |
-// denominator+=weight[t]*(this->size+this->obs[t]); |
|
| 691 |
-// } |
|
| 692 |
-// this->prob = numerator/denominator; // Update of r is now done with updated p |
|
| 693 |
-// double logp = log(this->prob); |
|
| 694 |
-// // dtime = clock() - time; |
|
| 695 |
-// // //FILE_LOG(logDEBUG1) << "updateP(): "<<dtime<< " clicks"; |
|
| 696 |
-// // Update of r with Newton Method |
|
| 697 |
-// rhere = this->size; |
|
| 698 |
-// dr = 0.00001; |
|
| 699 |
-// kmax = 20; |
|
| 700 |
-// // time = clock(); |
|
| 701 |
-// // Select strategy for computing digammas |
|
| 702 |
-// if (this->max_obs <= this->T) |
|
| 703 |
-// {
|
|
| 704 |
-// //FILE_LOG(logDEBUG3) << "Precomputing digammas in " << __func__ << " for every obs[t], because max(O)<=T"; |
|
| 705 |
-// std::vector<double> DigammaRplusX(this->max_obs+1); |
|
| 706 |
-// std::vector<double> DigammaRplusDRplusX(this->max_obs+1); |
|
| 707 |
-// for (int k=1; k<kmax; k++) |
|
| 708 |
-// {
|
|
| 709 |
-// Fr=dFrdr=0.0; |
|
| 710 |
-// DigammaR = digamma(rhere); // boost::math::digamma<>(rhere); |
|
| 711 |
-// DigammaRplusDR = digamma(rhere + dr); // boost::math::digamma<>(rhere+dr); |
|
| 712 |
-// // Precompute the digammas by iterating over all possible values of the observation vector |
|
| 713 |
-// for (int j=0; j<=this->max_obs; j++) |
|
| 714 |
-// {
|
|
| 715 |
-// DigammaRplusX[j] = digamma(rhere+j); |
|
| 716 |
-// DigammaRplusDRplusX[j] = digamma(rhere+dr+j); |
|
| 717 |
-// } |
|
| 718 |
-// for(int t=0; t<this->T; t++) |
|
| 719 |
-// {
|
|
| 720 |
-// if(this->obs[t]==0) |
|
| 721 |
-// {
|
|
| 722 |
-// Fr+=weight[t]*logp; |
|
| 723 |
-// //dFrdr+=0; |
|
| 724 |
-// } |
|
| 725 |
-// if(this->obs[t]!=0) |
|
| 726 |
-// {
|
|
| 727 |
-// Fr+=weight[t]*(logp-DigammaR+DigammaRplusX[(int)obs[t]]); |
|
| 728 |
-// dFrdr+=weight[t]/dr*(DigammaR-DigammaRplusDR+DigammaRplusDRplusX[(int)obs[t]]-DigammaRplusX[(int)obs[t]]); |
|
| 729 |
-// } |
|
| 730 |
-// } |
|
| 731 |
-// if(fabs(Fr)<eps) |
|
| 732 |
-// {
|
|
| 733 |
-// break; |
|
| 734 |
-// } |
|
| 735 |
-// if(Fr/dFrdr<rhere) |
|
| 736 |
-// {
|
|
| 737 |
-// rhere=rhere-Fr/dFrdr; |
|
| 738 |
-// } |
|
| 739 |
-// else if (Fr/dFrdr>=rhere) |
|
| 740 |
-// {
|
|
| 741 |
-// rhere=rhere/2.0; |
|
| 742 |
-// } |
|
| 743 |
-// } |
|
| 744 |
-// } |
|
| 745 |
-// else |
|
| 746 |
-// {
|
|
| 747 |
-// //FILE_LOG(logDEBUG2) << "Computing digammas in " << __func__ << " for every t, because max(O)>T"; |
|
| 748 |
-// double DigammaRplusX, DigammaRplusDRplusX; |
|
| 749 |
-// for (int k=1; k<kmax; k++) |
|
| 750 |
-// {
|
|
| 751 |
-// Fr = dFrdr = 0.0; |
|
| 752 |
-// DigammaR = digamma(rhere); // boost::math::digamma<>(rhere); |
|
| 753 |
-// DigammaRplusDR = digamma(rhere + dr); // boost::math::digamma<>(rhere+dr); |
|
| 754 |
-// for(int t=0; t<this->T; t++) |
|
| 755 |
-// {
|
|
| 756 |
-// DigammaRplusX = digamma(rhere+this->obs[t]); //boost::math::digamma<>(rhere+this->obs[ti]); |
|
| 757 |
-// DigammaRplusDRplusX = digamma(rhere+dr+this->obs[t]); // boost::math::digamma<>(rhere+dr+this->obs[ti]); |
|
| 758 |
-// if(this->obs[t]==0) |
|
| 759 |
-// {
|
|
| 760 |
-// Fr+=weight[t]*logp; |
|
| 761 |
-// //dFrdr+=0; |
|
| 762 |
-// } |
|
| 763 |
-// if(this->obs[t]!=0) |
|
| 764 |
-// {
|
|
| 765 |
-// Fr+=weight[t]*(logp-DigammaR+DigammaRplusX); |
|
| 766 |
-// dFrdr+=weight[t]/dr*(DigammaR-DigammaRplusDR+DigammaRplusDRplusX-DigammaRplusX); |
|
| 767 |
-// } |
|
| 768 |
-// } |
|
| 769 |
-// if(fabs(Fr)<eps) |
|
| 770 |
-// {
|
|
| 771 |
-// break; |
|
| 772 |
-// } |
|
| 773 |
-// if(Fr/dFrdr<rhere) |
|
| 774 |
-// {
|
|
| 775 |
-// rhere=rhere-Fr/dFrdr; |
|
| 776 |
-// } |
|
| 777 |
-// else if (Fr/dFrdr>=rhere) |
|
| 778 |
-// {
|
|
| 779 |
-// rhere=rhere/2.0; |
|
| 780 |
-// } |
|
| 781 |
-// } |
|
| 782 |
-// } |
|
| 783 |
-// this->size = rhere; |
|
| 784 |
-// //FILE_LOG(logDEBUG1) << "r = "<<this->size << ", p = "<<this->prob; |
|
| 785 |
-// |
|
| 786 |
-// // dtime = clock() - time; |
|
| 787 |
-// // //FILE_LOG(logDEBUG1) << "updateR(): "<<dtime<< " clicks"; |
|
| 788 |
-// |
|
| 789 |
-// } |
|
| 790 |
- |
|
| 791 | 678 |
void NegativeBinomial::copy(Density* other) |
| 792 | 679 |
{
|
| 793 | 680 |
//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__; |
| ... | ... |
@@ -163,7 +163,9 @@ void ScaleHMM::baumWelch(int* maxiter, int* maxtime, double* eps) |
| 163 | 163 |
double logPnew; |
| 164 | 164 |
|
| 165 | 165 |
// Parallelization settings |
| 166 |
+// #ifdef _OPENMP |
|
| 166 | 167 |
// omp_set_nested(1); |
| 168 |
+// #endif |
|
| 167 | 169 |
|
| 168 | 170 |
// measuring the time |
| 169 | 171 |
this->baumWelchStartTime_sec = time(NULL); |
| ... | ... |
@@ -984,7 +986,7 @@ void ScaleHMM::calc_densities() |
| 984 | 986 |
} |
| 985 | 987 |
|
| 986 | 988 |
// Check if the density for all states is close to zero and correct to prevent NaNs |
| 987 |
- double zero_cutoff = 1e-30; // 32-bit precision is 1.18e-38 |
|
| 989 |
+ double zero_cutoff = 1e-20; // 32-bit precision is 1.18e-38 |
|
| 988 | 990 |
std::vector<double> temp(this->N); |
| 989 | 991 |
// t=0 |
| 990 | 992 |
for (int iN=0; iN<this->N; iN++) |
