projects
/
samtools.git
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
Imported Upstream version 0.1.11
[samtools.git]
/
bcftools
/
prob1.c
diff --git
a/bcftools/prob1.c
b/bcftools/prob1.c
index 176a0fc96d4e0b12f14172bc7c0d68297e01a50a..8bf968f6adffd6a7e22058d6c2c5886f33f7a5c1 100644
(file)
--- a/
bcftools/prob1.c
+++ b/
bcftools/prob1.c
@@
-32,7
+32,7
@@
unsigned char seq_nt4_table[256] = {
};
struct __bcf_p1aux_t {
};
struct __bcf_p1aux_t {
- int n, M, n1, is_indel;
+ int n, M, n1, is_indel
, is_folded
;
double *q2p, *pdg; // pdg -> P(D|g)
double *phi, *phi_indel;
double *z, *zswap; // aux for afs
double *q2p, *pdg; // pdg -> P(D|g)
double *phi, *phi_indel;
double *z, *zswap; // aux for afs
@@
-43,6
+43,13
@@
struct __bcf_p1aux_t {
int PL_len;
};
int PL_len;
};
+static void fold_array(int M, double *x)
+{
+ int k;
+ for (k = 0; k < M/2; ++k)
+ x[k] = x[M-k] = (x[k] + x[M-k]) / 2.;
+}
+
void bcf_p1_indel_prior(bcf_p1aux_t *ma, double x)
{
int i;
void bcf_p1_indel_prior(bcf_p1aux_t *ma, double x)
{
int i;
@@
-155,6
+162,15
@@
int bcf_p1_set_n1(bcf_p1aux_t *b, int n1)
return 0;
}
return 0;
}
+void bcf_p1_set_folded(bcf_p1aux_t *p1a)
+{
+ if (p1a->n1 < 0) {
+ p1a->is_folded = 1;
+ fold_array(p1a->M, p1a->phi);
+ fold_array(p1a->M, p1a->phi_indel);
+ }
+}
+
void bcf_p1_destroy(bcf_p1aux_t *ma)
{
if (ma) {
void bcf_p1_destroy(bcf_p1aux_t *ma)
{
if (ma) {
@@
-373,7
+389,7
@@
int bcf_p1_cal(bcf1_t *b, bcf_p1aux_t *ma, bcf_p1rst_t *rst)
//
rst->rank0 = cal_pdg(b, ma);
rst->f_exp = mc_cal_afs(ma);
//
rst->rank0 = cal_pdg(b, ma);
rst->f_exp = mc_cal_afs(ma);
- rst->p_ref = ma->afs1[ma->M];
+ rst->p_ref = ma->
is_folded? ma->afs1[ma->M] + ma->afs1[0] : ma->
afs1[ma->M];
// calculate f_flat and f_em
for (k = 0, sum = 0.; k <= ma->M; ++k)
sum += (long double)ma->z[k];
// calculate f_flat and f_em
for (k = 0, sum = 0.; k <= ma->M; ++k)
sum += (long double)ma->z[k];
@@
-412,6
+428,7
@@
int bcf_p1_cal(bcf1_t *b, bcf_p1aux_t *ma, bcf_p1rst_t *rst)
void bcf_p1_dump_afs(bcf_p1aux_t *ma)
{
int k;
void bcf_p1_dump_afs(bcf_p1aux_t *ma)
{
int k;
+ if (ma->is_folded) fold_array(ma->M, ma->afs);
fprintf(stderr, "[afs]");
for (k = 0; k <= ma->M; ++k)
fprintf(stderr, " %d:%.3lf", k, ma->afs[ma->M - k]);
fprintf(stderr, "[afs]");
for (k = 0; k <= ma->M; ++k)
fprintf(stderr, " %d:%.3lf", k, ma->afs[ma->M - k]);