pub fn prob_related<G, E>(&mut self, gap_params: &G, emission_params: &E) -> LogProb
    where
        G: GapParameters + StartEndGapParameters,
        E: EmissionParameters,
    {
        for k in 0..2 {
            self.fm[k].clear();
            self.fx[k].clear();
            self.fy[k].clear();
            self.prob_cols.clear();

            self.fm[k].resize(emission_params.len_y() + 1, LogProb::ln_zero());
            self.fx[k].resize(emission_params.len_y() + 1, LogProb::ln_zero());
            self.fy[k].resize(emission_params.len_y() + 1, LogProb::ln_zero());

            if gap_params.free_end_gap_x() {
                let c = (emission_params.len_x() * 3).saturating_sub(self.prob_cols.capacity());
                self.prob_cols.reserve_exact(c);
            }
        }

        // cache probs
        let prob_no_gap = gap_params
            .prob_gap_x()
            .ln_add_exp(gap_params.prob_gap_y())
            .ln_one_minus_exp();
        let prob_no_gap_x_extend = gap_params.prob_gap_x_extend().ln_one_minus_exp();
        let prob_no_gap_y_extend = gap_params.prob_gap_y_extend().ln_one_minus_exp();
        let prob_gap_x = gap_params.prob_gap_x();
        let prob_gap_y = gap_params.prob_gap_y();
        let prob_gap_x_extend = gap_params.prob_gap_x_extend();
        let prob_gap_y_extend = gap_params.prob_gap_y_extend();
        // let do_gap_extend = prob_gap_x_extend != LogProb::ln_zero() &&
        //                     prob_gap_y_extend != LogProb::ln_zero();

        let mut prev = 0;
        let mut curr = 1;
        self.fm[prev][0] = LogProb::ln_one();

        // iterate over x
        for i in 0..emission_params.len_x() {
            // allow alignment to start from offset in x (if prob_start_gap_x is set accordingly)
            self.fm[prev][0] = gap_params.prob_start_gap_x(i);

            let prob_emit_x = emission_params.prob_emit_x(i);

            // iterate over y
            for j in 0..emission_params.len_y() {
                let j_ = j + 1;

                // match or mismatch
                self.fm[curr][j_] = emission_params.prob_emit_xy(i, j)
                    + LogProb::ln_sum_exp(&[
                        // coming from state M
                        prob_no_gap + self.fm[prev][j_ - 1],
                        // coming from state X
                        prob_no_gap_x_extend + self.fx[prev][j_ - 1],
                        // coming from state Y
                        prob_no_gap_y_extend + self.fy[prev][j_ - 1],
                    ]);

                // gap in y
                self.fx[curr][j_] = prob_emit_x
                    + (
                    // open gap
                    prob_gap_y + self.fm[prev][j_]
                ).ln_add_exp(
                        // extend gap
                        prob_gap_y_extend + self.fx[prev][j_],
                    );

                // gap in x
                self.fy[curr][j_] = emission_params.prob_emit_y(j)
                    + (
                    // open gap
                    prob_gap_x + self.fm[curr][j_ - 1]
                ).ln_add_exp(
                        // extend gap
                        prob_gap_x_extend + self.fy[curr][j_ - 1],
                    );
            }

            if gap_params.free_end_gap_x() {
                // Cache column probabilities or simply record the last probability.
                // We can put all of them in one array since we simply have to sum in the end.
                // This is also good for numerical stability.
                self.prob_cols.push(self.fm[curr].last().unwrap().clone());
                self.prob_cols.push(self.fx[curr].last().unwrap().clone());
                // TODO check removing this (we don't want open gaps in x):
                self.prob_cols.push(self.fy[curr].last().unwrap().clone());
            }

            // next column
            mem::swap(&mut curr, &mut prev);
            // reset next column to zeros
            for v in &mut self.fm[curr] {
                *v = LogProb::ln_zero();
            }
        }

        let p = if gap_params.free_end_gap_x() {
            LogProb::ln_sum_exp(&self.prob_cols)
        } else {
            LogProb::ln_sum_exp(&[
                *self.fm[prev].last().unwrap(),
                *self.fx[prev].last().unwrap(),
                *self.fy[prev].last().unwrap(),
            ])
        };
        // take the minimum with 1.0, because sum of paths can exceed probability 1.0
        // especially in case of repeats
        assert!(!p.is_nan());
        if p > LogProb::ln_one() {
            LogProb::ln_one()
        } else {
            p
        }
    }
/// Calculate the probability of sequence x being related to y via any alignment.
Calculate the probability of sequence x being related to y via any alignment.
Calculate the probability of sequence x being related to y via any alignment .
pub fn prob_related<G, E>(&mut self, gap_params: &G, emission_params: &E) -> LogProb
    where
        G: GapParameters + StartEndGapParameters,
        E: EmissionParameters,
    {
        for k in 0..2 {
            self.fm[k].clear();
            self.fx[k].clear();
            self.fy[k].clear();
            self.prob_cols.clear();
            self.fm[k].resize(emission_params.len_y() + 1, LogProb::ln_zero());
            self.fx[k].resize(emission_params.len_y() + 1, LogProb::ln_zero());
            self.fy[k].resize(emission_params.len_y() + 1, LogProb::ln_zero());
            if gap_params.free_end_gap_x() {
                let c = (emission_params.len_x() * 3).saturating_sub(self.prob_cols.capacity());
                self.prob_cols.reserve_exact(c);
            }
        }
        // cache probs
        let prob_no_gap = gap_params
            .prob_gap_x()
            .ln_add_exp(gap_params.prob_gap_y())
            .ln_one_minus_exp();
        let prob_no_gap_x_extend = gap_params.prob_gap_x_extend().ln_one_minus_exp();
        let prob_no_gap_y_extend = gap_params.prob_gap_y_extend().ln_one_minus_exp();
        let prob_gap_x = gap_params.prob_gap_x();
        let prob_gap_y = gap_params.prob_gap_y();
        let prob_gap_x_extend = gap_params.prob_gap_x_extend();
        let prob_gap_y_extend = gap_params.prob_gap_y_extend();
        // let do_gap_extend = prob_gap_x_extend != LogProb::ln_zero() &&
        //                     prob_gap_y_extend != LogProb::ln_zero();
        let mut prev = 0;
        let mut curr = 1;
        self.fm[prev][0] = LogProb::ln_one();
        // iterate over x
        for i in 0..emission_params.len_x() {
            // allow alignment to start from offset in x (if prob_start_gap_x is set accordingly)
            self.fm[prev][0] = gap_params.prob_start_gap_x(i);
            let prob_emit_x = emission_params.prob_emit_x(i);
            // iterate over y
            for j in 0..emission_params.len_y() {
                let j_ = j + 1;
                // match or mismatch
                self.fm[curr][j_] = emission_params.prob_emit_xy(i, j)
                    + LogProb::ln_sum_exp(&[
                        // coming from state M
                        prob_no_gap + self.fm[prev][j_ - 1],
                        // coming from state X
                        prob_no_gap_x_extend + self.fx[prev][j_ - 1],
                        // coming from state Y
                        prob_no_gap_y_extend + self.fy[prev][j_ - 1],
                    ]);
                // gap in y
                self.fx[curr][j_] = prob_emit_x
                    + (
                    // open gap
                    prob_gap_y + self.fm[prev][j_]
                ).ln_add_exp(
                        // extend gap
                        prob_gap_y_extend + self.fx[prev][j_],
                    );
                // gap in x
                self.fy[curr][j_] = emission_params.prob_emit_y(j)
                    + (
                    // open gap
                    prob_gap_x + self.fm[curr][j_ - 1]
                ).ln_add_exp(
                        // extend gap
                        prob_gap_x_extend + self.fy[curr][j_ - 1],
                    );
            }
            if gap_params.free_end_gap_x() {
                // Cache column probabilities or simply record the last probability.
                // We can put all of them in one array since we simply have to sum in the end.
                // This is also good for numerical stability.
                self.prob_cols.push(self.fm[curr].last().unwrap().clone());
                self.prob_cols.push(self.fx[curr].last().unwrap().clone());
                // TODO check removing this (we don't want open gaps in x):
                self.prob_cols.push(self.fy[curr].last().unwrap().clone());
            }
            // next column
            mem::swap(&mut curr, &mut prev);
            // reset next column to zeros
            for v in &mut self.fm[curr] {
                *v = LogProb::ln_zero();
            }
        }
        let p = if gap_params.free_end_gap_x() {
            LogProb::ln_sum_exp(&self.prob_cols)
        } else {
            LogProb::ln_sum_exp(&[
                *self.fm[prev].last().unwrap(),
                *self.fx[prev].last().unwrap(),
                *self.fy[prev].last().unwrap(),
            ])
        };
        // take the minimum with 1.0, because sum of paths can exceed probability 1.0
        // especially in case of repeats
        assert!(!p.is_nan());
        if p > LogProb::ln_one() {
            LogProb::ln_one()
        } else {
            p
        }
    }
pub fn prob_related < G , E > ( & mut self , gap_params : & G , emission_params : & E ) -> LogProb where G : GapParameters + StartEndGapParameters , E : EmissionParameters , { for k in 0 .. 2 { self . fm [ k ] . clear ( ) ; self . fx [ k ] . clear ( ) ; self . fy [ k ] . clear ( ) ; self . prob_cols . clear ( ) ; self . fm [ k ] . resize ( emission_params . len_y ( ) + 1 , LogProb :: ln_zero ( ) ) ; self . fx [ k ] . resize ( emission_params . len_y ( ) + 1 , LogProb :: ln_zero ( ) ) ; self . fy [ k ] . resize ( emission_params . len_y ( ) + 1 , LogProb :: ln_zero ( ) ) ; if gap_params . free_end_gap_x ( ) { let c = ( emission_params . len_x ( ) * 3 ) . saturating_sub ( self . prob_cols . capacity ( ) ) ; self . prob_cols . reserve_exact ( c ) ; } } // cache probs let prob_no_gap = gap_params . prob_gap_x ( ) . ln_add_exp ( gap_params . prob_gap_y ( ) ) . ln_one_minus_exp ( ) ; let prob_no_gap_x_extend = gap_params . prob_gap_x_extend ( ) . ln_one_minus_exp ( ) ; let prob_no_gap_y_extend = gap_params . prob_gap_y_extend ( ) . ln_one_minus_exp ( ) ; let prob_gap_x = gap_params . prob_gap_x ( ) ; let prob_gap_y = gap_params . prob_gap_y ( ) ; let prob_gap_x_extend = gap_params . prob_gap_x_extend ( ) ; let prob_gap_y_extend = gap_params . prob_gap_y_extend ( ) ; // let do_gap_extend = prob_gap_x_extend != LogProb::ln_zero() && //                     prob_gap_y_extend != LogProb::ln_zero(); let mut prev = 0 ; let mut curr = 1 ; self . fm [ prev ] [ 0 ] = LogProb :: ln_one ( ) ; // iterate over x for i in 0 .. emission_params . len_x ( ) { // allow alignment to start from offset in x (if prob_start_gap_x is set accordingly) self . fm [ prev ] [ 0 ] = gap_params . prob_start_gap_x ( i ) ; let prob_emit_x = emission_params . prob_emit_x ( i ) ; // iterate over y for j in 0 .. emission_params . len_y ( ) { let j_ = j + 1 ; // match or mismatch self . fm [ curr ] [ j_ ] = emission_params . prob_emit_xy ( i , j ) + LogProb :: ln_sum_exp ( & [ // coming from state M prob_no_gap + self . fm [ prev ] [ j_ - 1 ] , // coming from state X prob_no_gap_x_extend + self . fx [ prev ] [ j_ - 1 ] , // coming from state Y prob_no_gap_y_extend + self . fy [ prev ] [ j_ - 1 ] , ] ) ; // gap in y self . fx [ curr ] [ j_ ] = prob_emit_x + ( // open gap prob_gap_y + self . fm [ prev ] [ j_ ] ) . ln_add_exp ( // extend gap prob_gap_y_extend + self . fx [ prev ] [ j_ ] , ) ; // gap in x self . fy [ curr ] [ j_ ] = emission_params . prob_emit_y ( j ) + ( // open gap prob_gap_x + self . fm [ curr ] [ j_ - 1 ] ) . ln_add_exp ( // extend gap prob_gap_x_extend + self . fy [ curr ] [ j_ - 1 ] , ) ; } if gap_params . free_end_gap_x ( ) { // Cache column probabilities or simply record the last probability. // We can put all of them in one array since we simply have to sum in the end. // This is also good for numerical stability. self . prob_cols . push ( self . fm [ curr ] . last ( ) . unwrap ( ) . clone ( ) ) ; self . prob_cols . push ( self . fx [ curr ] . last ( ) . unwrap ( ) . clone ( ) ) ; // TODO check removing this (we don't want open gaps in x): self . prob_cols . push ( self . fy [ curr ] . last ( ) . unwrap ( ) . clone ( ) ) ; } // next column mem :: swap ( & mut curr , & mut prev ) ; // reset next column to zeros for v in & mut self . fm [ curr ] { * v = LogProb :: ln_zero ( ) ; } } let p = if gap_params . free_end_gap_x ( ) { LogProb :: ln_sum_exp ( & self . prob_cols ) } else { LogProb :: ln_sum_exp ( & [ * self . fm [ prev ] . last ( ) . unwrap ( ) , * self . fx [ prev ] . last ( ) . unwrap ( ) , * self . fy [ prev ] . last ( ) . unwrap ( ) , ] ) } ; // take the minimum with 1.0, because sum of paths can exceed probability 1.0 // especially in case of repeats assert ! ( p is_nan ( ) ) ; if p > LogProb :: ln_one ( ) { LogProb :: ln_one ( ) } else { p } }
Calculate the probability of sequence x being related to y via any alignment.
Calculate the probability of sequence x being related to y via any alignment .
&mut self
gap_params
&G
emission_params
&E
&mut self
gap_params
&G
emission_params
&E