#[macro_use]
extern crate lazy_static;
extern crate fnv;
pub mod cpr_cnsts;
use cpr_cnsts::*;
use std::f64::consts::LOG2_E;
use std::collections::HashMap;
use std::hash::BuildHasherDefault;
use fnv::FnvHasher;
use std::mem;
type BElm = u8;
type BSqSlc<'a> = &'a[BElm];
type Prb = Enrgy;
type CntxtDst = [Prb; PRB_DST_DM]; pub type CntxtDstSq = Vec<CntxtDst>;
type HshBSq = Vec<BElmHsh>;
type HshBSqSlc<'a> = &'a[BElmHsh];
type Hshr = BuildHasherDefault<FnvHasher>;
type Exps = HashMap<EnrgyInt, Prb, Hshr>;
const A: BElm = 'A' as BElm;
const U: BElm = 'U' as BElm;
const T: BElm = 'T' as BElm;
const G: BElm = 'G' as BElm;
const C: BElm = 'C' as BElm;
const LWR_A: BElm = 'a' as BElm;
const LWR_U: BElm = 'u' as BElm;
const LWR_T: BElm = 't' as BElm;
const LWR_G: BElm = 'g' as BElm;
const LWR_C: BElm = 'c' as BElm;
const PRB_DST_DM: usize = 6;
const INVRS_LG2_E: Prb = 1. / LOG2_E;
lazy_static! {
static ref HRPN: [Enrgy; 31] = {
let mut hrpn = [0.; 31];
for i in 0 .. 31 {
hrpn[i] = -HRPN_37[i] as Enrgy * 10. / K_T;
}
hrpn
};
static ref BLG: [Enrgy; 31] = {
let mut blg = [0.; 31];
for i in 0 .. 31 {
blg[i] = -BLG_37[i] as Enrgy * 10. / K_T;
}
blg
};
static ref INTRNL: [Enrgy; 31] = {
let mut intrnl = [0.; 31];
for i in 0 .. 31 {
intrnl[i] = -INT_LP_37[i] as Enrgy * 10. / K_T;
}
intrnl
};
static ref MSMTCH_I: [[[Enrgy; 5]; 5]; 7] = {
let mut msmtch_i = [[[0.; 5]; 5]; 7];
for i in 0 .. 7 {
for j in 0 .. 5 {
for k in 0 .. 5 {
msmtch_i[i][j][k] = -MSMTCH_I_37[i][j][k] as Enrgy * 10. / K_T;
}
}
}
msmtch_i
};
static ref MSMTCH_H: [[[Enrgy; 5]; 5]; 7] = {
let mut msmtch_h = [[[0.; 5]; 5]; 7];
for i in 0 .. 7 {
for j in 0 .. 5 {
for k in 0 .. 5 {
msmtch_h[i][j][k] = - MSMTCH_H_37[i][j][k] as Enrgy * 10. / K_T;
}
}
}
msmtch_h
};
static ref STCK: [[Enrgy; 7]; 7] = {
let mut stck = [[0.; 7]; 7];
for i in 0 .. 7 {
for j in 0 .. 7 {
stck[i][j] = -STCK_37[i][j] as Enrgy * 10. / K_T;
}
}
stck
};
static ref DNGL_5: [[Enrgy; 5]; 8] = {
let mut dngl_5 = [[0.; 5]; 8];
for i in 0 .. 8 {
for j in 0 .. 5 {
dngl_5[i][j] = -DNGL_5_37[i][j] as Enrgy * 10. / K_T;
}
}
dngl_5
};
static ref DNGL_3: [[Enrgy; 5]; 8] = {
let mut dngl_3 = [[0.; 5]; 8];
for i in 0 .. 8 {
for j in 0 .. 5 {
dngl_3[i][j] = -DNGL_3_37[i][j] as Enrgy * 10. / K_T;
if i > 2 {
dngl_3[i][j] += TRM_A_U;
}
}
}
dngl_3
};
static ref INT_11: [[[[Enrgy; 5]; 5]; 8]; 8] = {
let mut int_11 = [[[[0.; 5]; 5]; 8]; 8];
for i in 0 .. 8 {
for j in 0 .. 8 {
for k in 0 .. 5 {
for l in 0 .. 5 {
int_11[i][j][k][l] = -INT_11_37[i][j][k][l] as Enrgy * 10. / K_T;
}
}
}
}
int_11
};
static ref INT_21: [[[[[Enrgy; 5]; 5]; 5]; 8]; 8] = {
let mut int_21 = [[[[[0.; 5]; 5]; 5]; 8]; 8];
for i in 0 .. 8 {
for j in 0 .. 8 {
for k in 0 .. 5 {
for l in 0 .. 5 {
for m in 0 .. 5 {
int_21[i][j][k][l][m] = -INT_21_37[i][j][k][l][m] as Enrgy * 10. / K_T;
}
}
}
}
}
int_21
};
static ref INT_22: [[[[[[Enrgy; 5]; 5]; 5]; 5]; 8]; 8] = {
let mut int_22 = [[[[[[0.; 5]; 5]; 5]; 5]; 8]; 8];
for i in 0 .. 8 {
for j in 0 .. 8 {
for k in 0 .. 5 {
for l in 0 .. 5 {
for m in 0 .. 5 {
for n in 0 .. 5 {
int_22[i][j][k][l][m][n] = -INT_22_37[i][j][k][l][m][n] as Enrgy * 10. / K_T;
}
}
}
}
}
}
int_22
};
static ref NINIO: [Enrgy; MX_LP + 1] = {
let mut ninio = [0.; MX_LP + 1];
for i in 0 .. MX_LP + 1 {
ninio[i] = -(MX_NINIO as Enrgy).min((i as EnrgyInt * F_NINIO_37) as Enrgy) * 10. / K_T;
}
ninio
};
static ref EXPS: Exps = {
let mut exps = Exps::default();
for i in (1.1 * NGTV_INF) as EnrgyInt .. INF {
exps.insert(i, (i as Prb).exp());
}
exps
};
static ref LG2_CNST: i64 = (2 as i64).pow(EXP as u32 - 1) - 1;
}
const SGN: usize = 1;
const EXP: usize = 11;
const SGNFCNT: usize = 52;
#[inline]
pub fn cpr(sq: BSqSlc, mx_spn: usize, is_aprxmt: bool) -> CntxtDstSq {
let sq_ln = sq.len();
let mut hsh_sq = sq.iter().map(|&bs| {gt_bs_hsh(bs)}).collect::<HshBSq>();
hsh_sq.insert(0, 0);
let mut alph_enrgy_vrbls = EnrgyVrbls::new(sq_ln, mx_spn);
let mut bt_enrgy_vrbls = EnrgyVrbls::new(sq_ln, mx_spn);
gt_insd_vrbls(&hsh_sq, sq_ln, mx_spn, is_aprxmt, &mut alph_enrgy_vrbls);
gt_otsd_vrbls(&hsh_sq, sq_ln, mx_spn, is_aprxmt, &alph_enrgy_vrbls, &mut bt_enrgy_vrbls);
gt_strct_prf(&hsh_sq, sq_ln, mx_spn, is_aprxmt, &alph_enrgy_vrbls, &bt_enrgy_vrbls)
}
#[inline]
fn gt_bs_hsh(bs: BElm) -> BElmHsh {
match bs {
A => 1,
LWR_A => 1,
C => 2,
LWR_C => 2,
G => 3,
LWR_G => 3,
U => 4,
LWR_U => 4,
T => 4,
LWR_T => 4,
_ => 0,
}
}
#[inline]
fn gt_insd_vrbls(hsh_sq: HshBSqSlc, sq_ln: usize, mx_spn: usize, is_aprxmt: bool, alph_enrgy_vrbls: &mut EnrgyVrbls) {
for j in TRN + 1 .. sq_ln + 1 {
let mn_i = if j < mx_spn + 1 {0} else {j - (mx_spn + 1)};
for i in (mn_i .. j - TRN + 1).rev() {
let mut typ = BP_PR[hsh_sq[i + 1] as usize][hsh_sq[j] as usize];
let mut typ_2 = BP_PR[hsh_sq[i + 2] as usize][hsh_sq[j - 1] as usize];
let mut tmp = 0.;
let mut flg = false;
if typ != 0 {
typ_2 = RTYP[typ_2 as usize];
if alph_enrgy_vrbls.stm[i + 1][j - i - 2] != NGTV_INF {
if typ_2 != 0 {
tmp = alph_enrgy_vrbls.stm[i + 1][j - i - 2] + gt_lp_enrgy(typ, typ_2, i + 1, j, i + 2, j - 1, hsh_sq, is_aprxmt);
}
flg = true;
}
if alph_enrgy_vrbls.stm_end[i + 1][j - i - 2] != NGTV_INF {
tmp = if flg {lg_sm_exp(tmp, alph_enrgy_vrbls.stm_end[i + 1][j - i - 2], is_aprxmt)} else {alph_enrgy_vrbls.stm_end[i + 1][j - i - 2]};
flg = true;
}
alph_enrgy_vrbls.stm[i][j - i] = if !flg {NGTV_INF} else {tmp};
} else {
alph_enrgy_vrbls.stm[i][j - i] = NGTV_INF;
}
tmp = 0.;
flg = false;
for k in i + 1 .. j {
if alph_enrgy_vrbls.mlt_1[i][k - i] != NGTV_INF && alph_enrgy_vrbls.mlt_2[k][j - k] != NGTV_INF {
tmp = if !flg {alph_enrgy_vrbls.mlt_1[i][k - i] + alph_enrgy_vrbls.mlt_2[k][j - k]} else {lg_sm_exp(tmp, alph_enrgy_vrbls.mlt_1[i][k - i] + alph_enrgy_vrbls.mlt_2[k][j - k], is_aprxmt)};
flg = true;
}
}
alph_enrgy_vrbls.mlt_bf[i][j - i] = if !flg {NGTV_INF} else {tmp};
tmp = 0.;
flg = false;
if typ != 0 {
if alph_enrgy_vrbls.stm[i][j - i] != NGTV_INF {
tmp = alph_enrgy_vrbls.stm[i][j - i] + ML_INTRN + gt_dngl_enrgy(typ, i, j, sq_ln, hsh_sq);
flg = true;
}
}
if alph_enrgy_vrbls.mlt_2[i][j - i - 1] != NGTV_INF {
alph_enrgy_vrbls.mlt_2[i][j - i] = alph_enrgy_vrbls.mlt_2[i][j - i - 1] + ML_BS;
if flg {
alph_enrgy_vrbls.mlt_2[i][j - i] = lg_sm_exp(tmp, alph_enrgy_vrbls.mlt_2[i][j - i], is_aprxmt);
}
} else {
alph_enrgy_vrbls.mlt_2[i][j - i] = if !flg {NGTV_INF} else {tmp};
}
if alph_enrgy_vrbls.mlt_2[i][j - i] != NGTV_INF && alph_enrgy_vrbls.mlt_bf[i][j - i] != NGTV_INF {
alph_enrgy_vrbls.mlt_1[i][j - i] = lg_sm_exp(alph_enrgy_vrbls.mlt_2[i][j - i], alph_enrgy_vrbls.mlt_bf[i][j - i], is_aprxmt);
} else if alph_enrgy_vrbls.mlt_2[i][j - i] == NGTV_INF {
alph_enrgy_vrbls.mlt_1[i][j - i] = alph_enrgy_vrbls.mlt_bf[i][j - i];
} else if alph_enrgy_vrbls.mlt_bf[i][j - i] == NGTV_INF {
alph_enrgy_vrbls.mlt_1[i][j - i] = alph_enrgy_vrbls.mlt_2[i][j - i];
} else {
alph_enrgy_vrbls.mlt_1[i][j - i] = NGTV_INF;
}
flg = false;
if alph_enrgy_vrbls.mlt[i + 1][j - i - 1] != NGTV_INF {
alph_enrgy_vrbls.mlt[i][j - i] = alph_enrgy_vrbls.mlt[i + 1][j - i - 1] + ML_BS;
flg = true;
}
if flg {
if alph_enrgy_vrbls.mlt_bf[i][j - i] != NGTV_INF {
alph_enrgy_vrbls.mlt[i][j - i] = lg_sm_exp(alph_enrgy_vrbls.mlt[i][j - i], alph_enrgy_vrbls.mlt_bf[i][j - i], is_aprxmt);
}
} else {
alph_enrgy_vrbls.mlt[i][j - i] = alph_enrgy_vrbls.mlt_bf[i][j - i];
}
if j != sq_ln {
typ = BP_PR[hsh_sq[i] as usize][hsh_sq[j + 1] as usize];
if typ != 0 {
tmp = gt_hrpn_enrgy(typ, i, j + 1, hsh_sq, is_aprxmt);
let mx_p = if j < TRN + 2 || i + MX_LP < j - (TRN + 2) {i + MX_LP} else {j - (TRN + 2)};
for p in i .. mx_p + 1 {
let u_1 = p - i;
let mn_q = if j < MX_LP - u_1 || p + TRN + 2 > j - (MX_LP - u_1) {p + TRN + 2} else {j - (MX_LP - u_1)};
for q in mn_q .. j + 1 {
typ_2 = BP_PR[hsh_sq[p + 1] as usize][hsh_sq[q] as usize];
if alph_enrgy_vrbls.stm[p][q - p] != NGTV_INF {
if typ_2 != 0 && !(p == i && q == j) {
typ_2 = RTYP[typ_2 as usize];
tmp = lg_sm_exp(tmp, alph_enrgy_vrbls.stm[p][q - p] + gt_lp_enrgy(typ, typ_2, i, j + 1, p + 1, q, hsh_sq, is_aprxmt), is_aprxmt);
}
}
}
}
let tt = RTYP[typ as usize];
tmp = lg_sm_exp(tmp, alph_enrgy_vrbls.mlt[i][j - i] + ML_CLS + ML_INTRN + DNGL_3[tt as usize][hsh_sq[i + 1] as usize] + DNGL_5[tt as usize][hsh_sq[j] as usize], is_aprxmt);
alph_enrgy_vrbls.stm_end[i][j - i] = tmp;
} else {
alph_enrgy_vrbls.stm_end[i][j - i] = NGTV_INF;
}
}
}
}
for i in 1 .. sq_ln + 1 {
let mut tmp = alph_enrgy_vrbls.otr[i - 1];
let mn_p = if i < mx_spn + 1 {0} else {i - (mx_spn + 1)};
for p in mn_p .. i {
if alph_enrgy_vrbls.stm[p][i - p] != NGTV_INF {
let typ = BP_PR[hsh_sq[p + 1] as usize][hsh_sq[i] as usize];
let ao = alph_enrgy_vrbls.stm[p][i - p] + gt_dngl_enrgy(typ, p, i, sq_ln, hsh_sq);
tmp = lg_sm_exp(tmp, ao + alph_enrgy_vrbls.otr[p], is_aprxmt);
}
}
alph_enrgy_vrbls.otr[i] = tmp;
}
}
#[inline]
fn gt_otsd_vrbls(hsh_sq: HshBSqSlc, sq_ln: usize, mx_spn: usize, is_aprxmt: bool, alph_enrgy_vrbls: &EnrgyVrbls, bt_enrgy_vrbls: &mut EnrgyVrbls) {
for i in (0 .. sq_ln).rev() {
let mut tmp = bt_enrgy_vrbls.otr[i + 1];
let mx_p = if i + mx_spn + 1 > sq_ln {sq_ln} else {i + mx_spn + 1};
for p in i + 1 .. mx_p + 1 {
if alph_enrgy_vrbls.stm[i][p - i] != NGTV_INF {
let typ = BP_PR[hsh_sq[i + 1] as usize][hsh_sq[p] as usize];
let bo = alph_enrgy_vrbls.stm[i][p - i] + gt_dngl_enrgy(typ, i, p, sq_ln, &hsh_sq);
tmp = lg_sm_exp(tmp, bo + bt_enrgy_vrbls.otr[p], is_aprxmt);
}
}
bt_enrgy_vrbls.otr[i] = tmp;
}
for q in (TRN + 1 .. sq_ln + 1).rev() {
let mn_p = if q < mx_spn + 1 {0} else {q - (mx_spn + 1)};
for p in mn_p .. q - TRN + 1 {
let mut typ;
let mut typ_2;
let mut tmp = 0.;
let mut flg;
if p != 0 && q != sq_ln {
bt_enrgy_vrbls.stm_end[p][q - p] = if q - p >= mx_spn {NGTV_INF} else {bt_enrgy_vrbls.stm[p - 1][q - p + 2]};
flg = false;
if q - p + 1 <= mx_spn + 1 {
if bt_enrgy_vrbls.mlt[p - 1][q - p + 1] != NGTV_INF {
tmp = bt_enrgy_vrbls.mlt[p - 1][q - p + 1] + ML_BS;
flg = true;
}
}
typ = BP_PR[hsh_sq[p] as usize][hsh_sq[q + 1] as usize];
let tt = RTYP[typ as usize];
if flg {
if bt_enrgy_vrbls.stm_end[p][q - p] != NGTV_INF {
tmp = lg_sm_exp(tmp, bt_enrgy_vrbls.stm_end[p][q - p] + ML_CLS + ML_INTRN + DNGL_3[tt as usize][hsh_sq[p + 1] as usize] + DNGL_5[tt as usize][hsh_sq[q] as usize], is_aprxmt);
}
} else {
if bt_enrgy_vrbls.stm_end[p][q - p] != NGTV_INF {
tmp = bt_enrgy_vrbls.stm_end[p][q - p] + ML_CLS + ML_INTRN + DNGL_3[tt as usize][hsh_sq[p + 1] as usize] + DNGL_5[tt as usize][hsh_sq[q] as usize];
} else {
tmp = NGTV_INF;
}
}
bt_enrgy_vrbls.mlt[p][q - p] = tmp;
tmp = 0.;
flg = false;
let mx_k = if sq_ln < p + mx_spn {sq_ln} else {p + mx_spn};
for k in q + 1 .. mx_k + 1 {
if bt_enrgy_vrbls.mlt_bf[p][k - p] != NGTV_INF && alph_enrgy_vrbls.mlt_2[q][k - q] != NGTV_INF {
tmp = if !flg {bt_enrgy_vrbls.mlt_bf[p][k - p] + alph_enrgy_vrbls.mlt_2[q][k - q]} else {lg_sm_exp(tmp, bt_enrgy_vrbls.mlt_bf[p][k - p] + alph_enrgy_vrbls.mlt_2[q][k - q], is_aprxmt)};
flg = true;
}
}
bt_enrgy_vrbls.mlt_1[p][q - p] = if flg {tmp} else {NGTV_INF};
tmp = 0.;
flg = false;
if bt_enrgy_vrbls.mlt_1[p][q - p] != NGTV_INF {
tmp = bt_enrgy_vrbls.mlt_1[p][q - p];
flg = true;
}
if q - p <= mx_spn {
if bt_enrgy_vrbls.mlt_2[p][q - p + 1] != NGTV_INF {
tmp = if flg {lg_sm_exp(tmp, bt_enrgy_vrbls.mlt_2[p][q - p + 1] + ML_BS, is_aprxmt)} else {bt_enrgy_vrbls.mlt_2[p][q - p + 1] + ML_BS};
flg = true;
}
}
let mn_k = if q < mx_spn {0} else {q - mx_spn};
for k in mn_k .. p {
if bt_enrgy_vrbls.mlt_bf[k][q - k] != NGTV_INF && alph_enrgy_vrbls.mlt_1[k][p - k] != NGTV_INF {
tmp = if !flg {bt_enrgy_vrbls.mlt_bf[k][q - k] + alph_enrgy_vrbls.mlt_1[k][p - k]} else {lg_sm_exp(tmp, bt_enrgy_vrbls.mlt_bf[k][q - k] + alph_enrgy_vrbls.mlt_1[k][p - k], is_aprxmt)};
flg = true;
}
}
bt_enrgy_vrbls.mlt_2[p][q - p] = if !flg {NGTV_INF} else {tmp};
if bt_enrgy_vrbls.mlt_1[p][q - p] != NGTV_INF && bt_enrgy_vrbls.mlt[p][q - p] != NGTV_INF {
bt_enrgy_vrbls.mlt_bf[p][q - p] = lg_sm_exp(bt_enrgy_vrbls.mlt_1[p][q - p], bt_enrgy_vrbls.mlt[p][q - p], is_aprxmt);
} else if bt_enrgy_vrbls.mlt[p][q - p] == NGTV_INF {
bt_enrgy_vrbls.mlt_bf[p][q - p] = bt_enrgy_vrbls.mlt_1[p][q - p];
} else if bt_enrgy_vrbls.mlt_1[p][q - p] == NGTV_INF {
bt_enrgy_vrbls.mlt_bf[p][q - p] = bt_enrgy_vrbls.mlt[p][q - p];
} else {
bt_enrgy_vrbls.mlt_bf[p][q - p] = NGTV_INF;
}
}
typ_2 = BP_PR[hsh_sq[p + 1] as usize][hsh_sq[q] as usize];
if typ_2 != 0 {
tmp = alph_enrgy_vrbls.otr[p] + bt_enrgy_vrbls.otr[q] + gt_dngl_enrgy(typ_2, p, q, sq_ln, hsh_sq);
typ_2 = RTYP[typ_2 as usize];
let mn_i = if p < MX_LP || 1 > p - MX_LP {1} else {p - MX_LP};
for i in mn_i .. p + 1 {
let mx_j = if q + MX_LP - p + i < sq_ln - 1 {q + MX_LP - p + i} else {sq_ln - 1};
for j in q .. mx_j + 1 {
typ = BP_PR[hsh_sq[i] as usize][hsh_sq[j + 1] as usize];
if typ != 0 && !(i == p && j == q) {
if j - i <= mx_spn + 1 && bt_enrgy_vrbls.stm_end[i][j - i] != NGTV_INF {
tmp = lg_sm_exp(tmp, bt_enrgy_vrbls.stm_end[i][j - i] + gt_lp_enrgy(typ, typ_2, i, j + 1, p + 1, q, &hsh_sq, is_aprxmt), is_aprxmt);
}
}
}
}
if p != 0 && q != sq_ln {
typ = BP_PR[hsh_sq[p] as usize][hsh_sq[q + 1] as usize];
if typ != 0 {
if q - p + 2 <= mx_spn + 1 && bt_enrgy_vrbls.stm[p - 1][q - p + 2] != NGTV_INF {
tmp = lg_sm_exp(tmp, bt_enrgy_vrbls.stm[p - 1][q - p + 2] + gt_lp_enrgy(typ, typ_2, p, q + 1, p + 1, q, &hsh_sq, is_aprxmt), is_aprxmt);
}
}
}
bt_enrgy_vrbls.stm[p][q - p] = tmp;
if bt_enrgy_vrbls.mlt_2[p][q - p] != NGTV_INF {
typ_2 = RTYP[typ_2 as usize];
tmp = bt_enrgy_vrbls.mlt_2[p][q - p] + ML_INTRN + gt_dngl_enrgy(typ_2, p, q, sq_ln, hsh_sq);
bt_enrgy_vrbls.stm[p][q - p] = lg_sm_exp(tmp, bt_enrgy_vrbls.stm[p][q - p], is_aprxmt);
}
} else {
bt_enrgy_vrbls.stm[p][q - p] = NGTV_INF;
}
}
}
}
#[inline]
fn gt_hrpn_enrgy(typ: BElmHsh, x: usize, y: usize, hsh_sq: HshBSqSlc, is_aprxmt: bool) -> Enrgy {
let z = y - x - 1;
let mut q = if z <= 30 {HRPN[z]} else {HRPN[30] - LXC_37 * gt_lg2(z as Enrgy / 30., is_aprxmt) * INVRS_LG2_E * 10. / K_T};
if z != 3 {
q += MSMTCH_H[typ as usize][hsh_sq[x + 1] as usize][hsh_sq[y - 1] as usize];
} else {
if typ > 2 {q += TRM_A_U;}
}
q
}
#[inline]
fn gt_dngl_enrgy(typ: BElmHsh, x: usize, y: usize, sq_ln: usize, hsh_sq: HshBSqSlc) -> Enrgy {
let mut tmp = 0.;
if typ != 0 {
if x > 0 {tmp += DNGL_5[typ as usize][hsh_sq[x] as usize];}
if y < sq_ln {tmp += DNGL_3[typ as usize][hsh_sq[y + 1] as usize];}
if y == sq_ln && typ > 2 {tmp += TRM_A_U;}
}
tmp
}
#[inline]
fn lg_sm_exp(x: Enrgy, y: Enrgy, is_aprxmt: bool) -> Enrgy {
if x > y {x + gt_lg2(gt_exp(y - x, is_aprxmt) + 1., is_aprxmt) * INVRS_LG2_E} else {y + gt_lg2(gt_exp(x - y, is_aprxmt) + 1., is_aprxmt) * INVRS_LG2_E}
}
#[inline]
fn gt_lp_enrgy(typ: BElmHsh, typ_2: BElmHsh, i: usize, j: usize, p: usize, q: usize, hsh_sq: HshBSqSlc, is_aprxmt: bool) -> Enrgy {
let mut z;
let (u_1, u_2) = (p - i - 1, j - q - 1);
if u_1 == 0 && u_2 == 0 {
z = STCK[typ as usize][typ_2 as usize];
} else {
if u_1 == 0 || u_2 == 0 {
let u = if u_1 == 0 {u_2} else {u_1};
z = if u <= 30 {BLG[u]} else {BLG[30] - LXC_37 * gt_lg2(u as Enrgy / 30., is_aprxmt) * INVRS_LG2_E * 10. / K_T};
if u == 1 {
z += STCK[typ as usize][typ_2 as usize];
} else {
if typ > 2 {z += TRM_A_U;}
if typ_2 > 2 {z += TRM_A_U;}
}
} else {
if u_1 + u_2 == 2 {
z = INT_11[typ as usize][typ_2 as usize][hsh_sq[i + 1] as usize][hsh_sq[j - 1] as usize];
} else if u_1 == 1 && u_2 == 2 {
z = INT_21[typ as usize][typ_2 as usize][hsh_sq[i + 1] as usize][hsh_sq[q + 1] as usize][hsh_sq[j - 1] as usize];
} else if u_1 == 2 && u_2 == 1 {
z = INT_21[typ_2 as usize][typ as usize][hsh_sq[q + 1] as usize][hsh_sq[i + 1] as usize][hsh_sq[p - 1] as usize]
} else if u_1 == 2 && u_2 == 2 {
z = INT_22[typ as usize][typ_2 as usize][hsh_sq[i + 1] as usize][hsh_sq[p - 1] as usize][hsh_sq[q + 1] as usize][hsh_sq[j - 1] as usize];
} else {
z = INTRNL[u_1 + u_2] + MSMTCH_I[typ as usize][hsh_sq[i + 1] as usize][hsh_sq[j - 1] as usize] + MSMTCH_I[typ_2 as usize][hsh_sq[q + 1] as usize][hsh_sq[p - 1] as usize];
z += NINIO[if u_1 > u_2 {u_1 - u_2} else {u_2 - u_1}];
}
}
}
z
}
#[inline]
fn gt_strct_prf(hsh_sq: HshBSqSlc, sq_ln: usize, mx_spn: usize, is_aprxmt: bool, alph_enrgy_vrbls: &EnrgyVrbls, bt_enrgy_vrbls: &EnrgyVrbls) -> CntxtDstSq {
let mut strct_prf = vec![[0.; PRB_DST_DM]; sq_ln];
let pf = alph_enrgy_vrbls.otr[sq_ln];
if pf >= -690. && pf <= 690. {
gt_blg_intrnl_prbs(hsh_sq, sq_ln, mx_spn, is_aprxmt, alph_enrgy_vrbls, bt_enrgy_vrbls, &mut strct_prf);
} else {
gt_lg_sm_blg_intrnl_prbs(hsh_sq, sq_ln, mx_spn, is_aprxmt, alph_enrgy_vrbls, bt_enrgy_vrbls, &mut strct_prf);
}
gt_hrpn_prbs(hsh_sq, sq_ln, mx_spn, is_aprxmt, alph_enrgy_vrbls, bt_enrgy_vrbls, &mut strct_prf);
for i in 1 .. sq_ln + 1 {
strct_prf[i - 1][3] = gt_extrr_prb(i, sq_ln, is_aprxmt, alph_enrgy_vrbls, bt_enrgy_vrbls);
strct_prf[i - 1][4] = gt_mlt_prb(i, sq_ln, mx_spn, is_aprxmt, alph_enrgy_vrbls, bt_enrgy_vrbls);
strct_prf[i - 1][5] = strct_prf[i - 1].iter().fold(1., |acmltr, prb| acmltr - prb);
}
strct_prf
}
#[inline]
fn gt_extrr_prb(i: usize, sq_ln: usize, is_aprxmt: bool, alph_enrgy_vrbls: &EnrgyVrbls, bt_enrgy_vrbls: &EnrgyVrbls) -> Prb {
gt_exp(alph_enrgy_vrbls.otr[i - 1] + bt_enrgy_vrbls.otr[i] - alph_enrgy_vrbls.otr[sq_ln], is_aprxmt)
}
#[inline]
fn gt_mlt_prb(x: usize, sq_ln: usize, mx_spn: usize, is_aprxmt: bool, alph_enrgy_vrbls: &EnrgyVrbls, bt_enrgy_vrbls: &EnrgyVrbls) -> Prb {
let mut tmp = 0.;
let mut flg = false;
let mx_i = if x + mx_spn < sq_ln {x + mx_spn} else {sq_ln};
for i in x .. mx_i + 1 {
if bt_enrgy_vrbls.mlt[x - 1][i - x + 1] != NGTV_INF && alph_enrgy_vrbls.mlt[x][i - x] != NGTV_INF {
tmp = if !flg {bt_enrgy_vrbls.mlt[x - 1][i - x + 1] + alph_enrgy_vrbls.mlt[x][i - x]} else {lg_sm_exp(tmp, bt_enrgy_vrbls.mlt[x - 1][i - x + 1] + alph_enrgy_vrbls.mlt[x][i - x], is_aprxmt)};
flg = true;
}
}
let mn_i = if x < mx_spn {0} else {x - mx_spn};
for i in mn_i .. x {
if bt_enrgy_vrbls.mlt_2[i][x - i] != NGTV_INF && alph_enrgy_vrbls.mlt_2[i][x - i - 1] != NGTV_INF {
tmp = if !flg {bt_enrgy_vrbls.mlt_2[i][x - i] + alph_enrgy_vrbls.mlt_2[i][x - i - 1]} else {lg_sm_exp(tmp, bt_enrgy_vrbls.mlt_2[i][x - i] + alph_enrgy_vrbls.mlt_2[i][x - i - 1], is_aprxmt)};
flg = true;
}
}
if flg {gt_exp(tmp - alph_enrgy_vrbls.otr[sq_ln], is_aprxmt)} else {0.}
}
#[inline]
fn gt_hrpn_prbs(hsh_sq: HshBSqSlc, sq_ln: usize, mx_spn: usize, is_aprxmt: bool, alph_enrgy_vrbls: &EnrgyVrbls, bt_enrgy_vrbls: &EnrgyVrbls, strct_prf: &mut CntxtDstSq) {
for x in 1 .. sq_ln + 1 {
let mut tmp = 0.;
let mut typ;
let mut flg = false;
let mut h_enrgy;
let mn_i = if x < mx_spn || 1 > x - mx_spn {1} else {x - mx_spn};
for i in mn_i .. x {
let mx_j = if i + mx_spn < sq_ln {i + mx_spn} else {sq_ln};
for j in x + 1 .. mx_j + 1 {
typ = BP_PR[hsh_sq[i] as usize][hsh_sq[j] as usize];
if bt_enrgy_vrbls.stm_end[i][j - i - 1] != NGTV_INF {
h_enrgy = bt_enrgy_vrbls.stm_end[i][j - i - 1] + gt_hrpn_enrgy(typ, i, j, hsh_sq, is_aprxmt);
tmp = if flg {lg_sm_exp(tmp, h_enrgy, is_aprxmt)} else {h_enrgy};
flg = true;
}
}
}
if flg {
strct_prf[x - 1][2] = gt_exp(tmp - alph_enrgy_vrbls.otr[sq_ln], is_aprxmt);
} else {
strct_prf[x - 1][2] = 0.;
}
}
}
#[inline]
fn gt_blg_intrnl_prbs(hsh_sq: HshBSqSlc, sq_ln: usize, mx_spn: usize, is_aprxmt: bool, alph_enrgy_vrbls: &EnrgyVrbls, bt_enrgy_vrbls: &EnrgyVrbls, strct_prf: &mut CntxtDstSq) {
let mut tmp;
let mut typ;
let mut typ_2;
for i in 1 .. sq_ln - TRN - 2 {
let mx_j = if i + mx_spn < sq_ln {i + mx_spn} else {sq_ln};
for j in i + TRN + 3 .. mx_j + 1 {
typ = BP_PR[hsh_sq[i] as usize][hsh_sq[j] as usize];
if typ != 0 {
let mx_p = if i + MX_LP + 1 < j - TRN - 2 {i + MX_LP + 1} else {j - TRN - 2};
for p in i + 1 .. mx_p + 1 {
let u_1 = p - i - 1;
let mn_q = if j < MX_LP - u_1 + 1 || p + TRN + 1 > j - (MX_LP - u_1 + 1) {p + TRN + 1} else {j - (MX_LP - u_1 + 1)};
for q in mn_q .. j {
typ_2 = BP_PR[hsh_sq[p] as usize][hsh_sq[q] as usize];
if typ_2 != 0 && !(p == i + 1 && q == j - 1) {
typ_2 = RTYP[typ_2 as usize];
if bt_enrgy_vrbls.stm_end[i][j - i - 1] != NGTV_INF && alph_enrgy_vrbls.stm[p - 1][q - p + 1] != NGTV_INF {
tmp = gt_exp(bt_enrgy_vrbls.stm_end[i][j - i - 1] + gt_lp_enrgy(typ, typ_2, i, j, p, q, hsh_sq, is_aprxmt) + alph_enrgy_vrbls.stm[p - 1][q - p + 1], is_aprxmt);
for k in i + 1 .. p {
if j == q + 1 {
strct_prf[k - 1][0] += tmp;
} else {
strct_prf[k - 1][1] += tmp;
}
}
for k in q + 1 .. j {
if i == p - 1 {
strct_prf[k - 1][0] += tmp;
} else {
strct_prf[k - 1][1] += tmp;
}
}
}
}
}
}
}
}
}
for i in 0 .. sq_ln {
if strct_prf[i][0] != 0. {
strct_prf[i][0] /= gt_exp(alph_enrgy_vrbls.otr[sq_ln], is_aprxmt);
}
if strct_prf[i][1] != 0. {
strct_prf[i][1] /= gt_exp(alph_enrgy_vrbls.otr[sq_ln], is_aprxmt);
}
}
}
#[inline]
fn gt_lg_sm_blg_intrnl_prbs(hsh_sq: HshBSqSlc, sq_ln: usize, mx_spn: usize, is_aprxmt: bool, alph_enrgy_vrbls: &EnrgyVrbls, bt_enrgy_vrbls: &EnrgyVrbls, strct_prf: &mut CntxtDstSq) {
let mut tmp;
let mut typ;
let mut typ_2;
let mut blg_flgs = vec![false; sq_ln];
let mut intrnl_flgs = blg_flgs.clone();
for i in 1 .. sq_ln - TRN - 2 {
let mx_j = if i + mx_spn < sq_ln {i + mx_spn} else {sq_ln};
for j in i + TRN + 3 .. mx_j + 1 {
typ = BP_PR[hsh_sq[i] as usize][hsh_sq[j] as usize];
if typ != 0 {
let mx_p = if i + MX_LP + 1 < j - TRN - 2 {i + MX_LP + 1} else {j - TRN - 2};
for p in i + 1 .. mx_p + 1 {
let u_1 = p - i - 1;
let mn_q = if j < MX_LP - u_1 + 1 || p + TRN + 1 > j - (MX_LP - u_1 + 1) {p + TRN + 1} else {j - (MX_LP - u_1 + 1)};
for q in mn_q .. j {
typ_2 = BP_PR[hsh_sq[p] as usize][hsh_sq[q] as usize];
if typ_2 != 0 && !(p == i + 1 && q == j - 1) {
typ_2 = RTYP[typ_2 as usize];
if bt_enrgy_vrbls.stm_end[i][j - i - 1] != NGTV_INF && alph_enrgy_vrbls.stm[p - 1][q - p + 1] != NGTV_INF {
tmp = bt_enrgy_vrbls.stm_end[i][j - i - 1] + gt_lp_enrgy(typ, typ_2, i, j, p, q, hsh_sq, is_aprxmt) + alph_enrgy_vrbls.stm[p - 1][q - p + 1];
for k in i + 1 .. p {
if j == q + 1 {
strct_prf[k - 1][0] = if blg_flgs[k - 1] {lg_sm_exp(strct_prf[k - 1][0], tmp, is_aprxmt)} else {tmp};
blg_flgs[k - 1] = true;
} else {
strct_prf[k - 1][1] = if intrnl_flgs[k - 1] {lg_sm_exp(strct_prf[k - 1][1], tmp, is_aprxmt)} else {tmp};
intrnl_flgs[k - 1] = true;
}
}
for k in q + 1 .. j {
if i == p - 1 {
strct_prf[k - 1][0] = if blg_flgs[k - 1] {lg_sm_exp(strct_prf[k-1][0], tmp, is_aprxmt)} else {tmp};
blg_flgs[k - 1] = true;
} else {
strct_prf[k - 1][1] = if intrnl_flgs[k - 1] {lg_sm_exp(strct_prf[k - 1][1], tmp, is_aprxmt)} else {tmp};
intrnl_flgs[k - 1] = true;
}
}
}
}
}
}
}
}
}
for i in 0 .. sq_ln {
if blg_flgs[i] {
strct_prf[i][0] = gt_exp(strct_prf[i][0] - alph_enrgy_vrbls.otr[sq_ln], is_aprxmt);
}
if intrnl_flgs[i] {
strct_prf[i][1] = gt_exp(strct_prf[i][1] - alph_enrgy_vrbls.otr[sq_ln], is_aprxmt);
}
}
}
#[inline]
fn gt_exp(x: Prb, is_aprxmt: bool) -> Prb {
if is_aprxmt {
let rnd_x = x.round();
let df = x - rnd_x;
let exp_rnd_x = EXPS[&(rnd_x as EnrgyInt)];
exp_rnd_x * (1. + df * (1. + df * (0.5 + df * (1. / 6. + 1. / 24. * df))))
} else {
x.exp()
}
}
#[inline]
fn gt_lg2(x: Prb, is_aprxmt: bool) -> Prb {
if is_aprxmt {
let (sgn, exp, sgnfcnt) = dcmps(x);
debug_assert!(sgn == 0 && 1 <= exp && exp <= (2 as u64).pow(EXP as u32) - 1);
let hgh_bt = (sgnfcnt >> (SGNFCNT - 1)) & 1;
let ad_exp = (exp + hgh_bt) as i64 - *LG2_CNST;
let nrmls = rcmps(0, 0x3FF ^ hgh_bt, sgnfcnt) - 1.;
const A: Prb = -0.6296735;
const B: Prb = 1.466967;
ad_exp as Prb + nrmls * (B + A * nrmls)
} else {
x.log2()
}
}
#[inline]
fn dcmps(x: Prb) -> (u64, u64, u64) {
let bts: u64 = unsafe {mem::transmute(x)};
macro_rules! msk{
($crnt: expr => $($othr: expr),*) => {
(bts >> (0 $(+ $othr)*)) & ((1 << $crnt) - 1)
}
}
(
msk!(SGN => EXP, SGNFCNT),
msk!(EXP => SGNFCNT),
msk!(SGNFCNT =>)
)
}
#[inline]
fn rcmps(sgn: u64, exp: u64, sgnfcnt: u64) -> Prb {
debug_assert!(sgn <= 1);
debug_assert!(exp <= (2 as u64).pow(EXP as u32) - 1);
debug_assert!(sgnfcnt < 1 << (SGNFCNT + 1));
macro_rules! unmsk {
($x: expr => $($othr: expr),*) => {
$x << (0 $(+ $othr)*)
}
}
let bts =
unmsk!(sgn => EXP, SGNFCNT) |
unmsk!(exp => SGNFCNT) |
unmsk!(sgnfcnt => );
unsafe {mem::transmute(bts)}
}
struct EnrgyVrbls {
otr: Vec<Enrgy>,
stm: Vec<Vec<Enrgy>>,
stm_end: Vec<Vec<Enrgy>>,
mlt: Vec<Vec<Enrgy>>,
mlt_bf: Vec<Vec<Enrgy>>,
mlt_1: Vec<Vec<Enrgy>>,
mlt_2: Vec<Vec<Enrgy>>,
}
impl EnrgyVrbls {
fn new(sq_ln: usize, mx_spn: usize) -> EnrgyVrbls {
let inr_vrbl = vec![vec![NGTV_INF; mx_spn + 2]; sq_ln + 1];
EnrgyVrbls {
otr: vec![0.; sq_ln + 1],
stm: inr_vrbl.clone(),
stm_end: inr_vrbl.clone(),
mlt: inr_vrbl.clone(),
mlt_bf: inr_vrbl.clone(),
mlt_1: inr_vrbl.clone(),
mlt_2: inr_vrbl,
}
}
}