Commit 793934aa authored by Mathieu Giraud's avatar Mathieu Giraud

core/kmerstore.h, core/segment.cpp: move getProbabilityAtLeastOrAbove to IKmerStore

parent 100c27ad
......@@ -6,6 +6,7 @@
#include <list>
#include <stdexcept>
#include <stdint.h>
#include <math.h>
#include "fasta.h"
#include "tools.h"
......@@ -84,6 +85,11 @@ public:
*/
float getIndexLoad() const;
/**
* @return probability that the number of kmers is 'at_least' or more in a sequence of length 'length'
*/
double getProbabilityAtLeastOrAbove(int at_least, int length) const;
/**
* @return the value of k
*/
......@@ -235,6 +241,28 @@ float IKmerStore<T>::getIndexLoad() const {
return nb_kmers_inserted*1. / (1 << (2 * k));
}
template<class T>
double IKmerStore<T>::getProbabilityAtLeastOrAbove(int at_least, int length) const {
// n: number of kmers in the sequence
int n = length - getS() + 1;
float index_load = getIndexLoad() ;
double proba = 0;
double probability_having_system = pow(index_load, at_least);
double probability_not_having_system = pow(1 - index_load, n - at_least);
for (int i=at_least; i<=n; i++) {
proba += nChoosek(n, i) * probability_having_system * probability_not_having_system;
probability_having_system *= index_load;
probability_not_having_system /= (1 - index_load);
}
return proba;
}
template<class T>
int IKmerStore<T>::getK() const {
return k;
......
......@@ -27,7 +27,6 @@
#include "affectanalyser.h"
#include <sstream>
#include <string>
#include <math.h>
Segmenter::~Segmenter() {}
......@@ -366,20 +365,8 @@ KmerMultiSegmenter::KmerMultiSegmenter(Sequence seq, MultiGermline *multigermlin
double KmerSegmenter::getProbabilityAtLeastOrAbove(int at_least) const {
// n: number of kmers in the sequence
int n = getSequence().sequence.size() - getKmerAffectAnalyser()->getIndex().getS() + 1;
float index_load = getKmerAffectAnalyser()->getIndex().getIndexLoad() ;
double proba = 0;
double probability_having_system = pow(index_load, at_least);
double probability_not_having_system = pow(1 - index_load, n - at_least);
for (int i=at_least; i<=n; i++) {
proba += nChoosek(n, i) * probability_having_system * probability_not_having_system;
probability_having_system *= index_load;
probability_not_having_system /= (1 - index_load);
}
return proba;
return getKmerAffectAnalyser()->getIndex().getProbabilityAtLeastOrAbove(at_least,
getSequence().sequence.size());
}
double KmerMultiSegmenter::getNbExpected() const {
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment