ViennaRNA-bindings-0.233.1.1: C/ViennaRNA/data_structures.c
/** \file data_structures.c **/
/*
Data structure creation/destruction
This file contains everything which is necessary to
obtain and destroy datastructures used in the folding
recurrences throughout the VienneRNA paclage
c Ronny Lorenx
Vienna RNA package
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#include <limits.h>
#include "ViennaRNA/utils.h"
#include "ViennaRNA/structure_utils.h"
#include "ViennaRNA/energy_par.h"
#include "ViennaRNA/data_structures.h"
#include "ViennaRNA/fold_vars.h"
#include "ViennaRNA/params.h"
#include "ViennaRNA/gquad.h"
#include "ViennaRNA/aln_util.h"
#include "ViennaRNA/ribo.h"
#include "ViennaRNA/constraints.h"
#include "ViennaRNA/part_func.h"
#include "ViennaRNA/cofold.h"
#include "ViennaRNA/mm.h"
#ifdef _OPENMP
#include <omp.h>
#endif
/*
#################################
# PRIVATE MACROS #
#################################
*/
#define WITH_PTYPE 1L /* passed to set_fold_compound() to indicate that we need to set vc->ptype */
#define WITH_PTYPE_COMPAT 2L /* passed to set_fold_compound() to indicate that we need to set vc->ptype_compat */
/*
#################################
# GLOBAL VARIABLES #
#################################
*/
/*
#################################
# PRIVATE VARIABLES #
#################################
*/
/*
#################################
# PRIVATE FUNCTION DECLARATIONS #
#################################
*/
PRIVATE void set_fold_compound(vrna_fold_compound_t *vc, vrna_md_t *md_p, unsigned int options, unsigned int aux);
PRIVATE void make_pscores(vrna_fold_compound_t *vc);
PRIVATE void add_params(vrna_fold_compound_t *vc, vrna_md_t *md_p, unsigned int options);
/*
#################################
# BEGIN OF FUNCTION DEFINITIONS #
#################################
*/
PUBLIC void
vrna_fold_compound_free(vrna_fold_compound_t *vc){
int s;
if(vc){
/* first destroy common attributes */
vrna_mx_mfe_free(vc);
vrna_mx_pf_free(vc);
free(vc->iindx);
free(vc->jindx);
free(vc->params);
free(vc->exp_params);
free(vc->strand_number);
vrna_hc_free(vc->hc);
vrna_ud_remove(vc);
/* now distinguish the vc type */
switch(vc->type){
case VRNA_FC_TYPE_SINGLE: free(vc->sequence);
free(vc->sequence_encoding);
free(vc->sequence_encoding2);
free(vc->ptype);
free(vc->ptype_pf_compat);
vrna_sc_free(vc->sc);
break;
case VRNA_FC_TYPE_COMPARATIVE: for(s=0;s<vc->n_seq;s++){
free(vc->sequences[s]);
free(vc->S[s]);
free(vc->S5[s]);
free(vc->S3[s]);
free(vc->Ss[s]);
free(vc->a2s[s]);
}
free(vc->sequences);
free(vc->cons_seq);
free(vc->S_cons);
free(vc->S);
free(vc->S5);
free(vc->S3);
free(vc->Ss);
free(vc->a2s);
free(vc->pscore);
free(vc->pscore_pf_compat);
if(vc->scs){
for(s=0;s<vc->n_seq;s++)
vrna_sc_free(vc->scs[s]);
free(vc->scs);
}
break;
default: /* do nothing */
break;
}
/* free Distance Class Partitioning stuff (should be NULL if not used) */
free(vc->reference_pt1);
free(vc->reference_pt2);
free(vc->referenceBPs1);
free(vc->referenceBPs2);
free(vc->bpdist);
free(vc->mm1);
free(vc->mm2);
/* free local folding related stuff (should be NULL if not used) */
if(vc->ptype_local){
for (s=0; (s < vc->window_size + 5) && (s <= vc->length); s++){
free(vc->ptype_local[s]);
}
free(vc->ptype_local);
}
if(vc->free_auxdata)
vc->free_auxdata(vc->auxdata);
free(vc);
}
}
PUBLIC vrna_fold_compound_t *
vrna_fold_compound( const char *sequence,
vrna_md_t *md_p,
unsigned int options){
int i;
unsigned int length, aux_options;
vrna_fold_compound_t *vc;
vrna_md_t md;
if(sequence == NULL) return NULL;
/* sanity check */
length = strlen(sequence);
if(length == 0)
vrna_message_error("vrna_fold_compound@data_structures.c: sequence length must be greater 0");
if(length > vrna_sequence_length_max(options))
vrna_message_error("vrna_fold_compound@data_structures.c: sequence length of %d exceeds addressable range", length);
vc = vrna_alloc(sizeof(vrna_fold_compound_t));
vc->type = VRNA_FC_TYPE_SINGLE;
vc->length = length;
vc->sequence = strdup(sequence);
aux_options = 0L;
/* get a copy of the model details */
if(md_p)
md = *md_p;
else
vrna_md_set_default(&md);
if(options & VRNA_OPTION_WINDOW){ /* sliding window structure prediction */
if(md.window_size <= 0)
md.window_size = (int)vc->length;
else if(md.window_size > (int)vc->length)
md.window_size = (int)vc->length;
vc->window_size = md.window_size;
if((md.max_bp_span <= 0) || (md.max_bp_span > md.window_size))
md.max_bp_span = md.window_size;
set_fold_compound(vc, &md, options, aux_options);
vc->ptype_local = vrna_alloc(sizeof(char *)*(vc->length+1));
for (i = (int)vc->length; ( i > (int)vc->length - vc->window_size - 5) && (i >= 0); i--){
vc->ptype_local[i] = vrna_alloc(sizeof(char)*(vc->window_size+5));
}
if(!(options & VRNA_OPTION_EVAL_ONLY)){
/* add default hard constraints */
/* vrna_hc_init(vc); */ /* no hard constraints in Lfold, yet! */
/* add DP matrices */
vrna_mx_add(vc, VRNA_MX_WINDOW, options);
}
} else { /* regular global structure prediction */
/* set window size to entire sequence */
md.window_size = (int)vc->length;
aux_options |= WITH_PTYPE;
if(options & VRNA_OPTION_PF)
aux_options |= WITH_PTYPE_COMPAT;
set_fold_compound(vc, &md, options, aux_options);
if(!(options & VRNA_OPTION_EVAL_ONLY)){
/* add default hard constraints */
vrna_hc_init(vc);
/* add DP matrices (if required) */
vrna_mx_add(vc, VRNA_MX_DEFAULT, options);
}
}
return vc;
}
PUBLIC vrna_fold_compound_t *
vrna_fold_compound_comparative( const char **sequences,
vrna_md_t *md_p,
unsigned int options){
int s, n_seq, length;
vrna_fold_compound_t *vc;
vrna_md_t md;
unsigned int aux_options;
aux_options = 0U;
if(sequences == NULL) return NULL;
for(s=0;sequences[s];s++); /* count the sequences */
n_seq = s;
length = strlen(sequences[0]);
/* sanity check */
if(length == 0)
vrna_message_error("vrna_fold_compound_comparative@data_structures.c: sequence length must be greater 0");
else if(length > vrna_sequence_length_max(options))
vrna_message_error("vrna_fold_compound_comparative@data_structures.c: sequence length of %d exceeds addressable range", length);
for(s = 0; s < n_seq; s++)
if(strlen(sequences[s]) != length)
vrna_message_error("vrna_fold_compound_comparative@data_structures.c: uneqal sequence lengths in alignment");
vc = vrna_alloc(sizeof(vrna_fold_compound_t));
vc->type = VRNA_FC_TYPE_COMPARATIVE;
vc->n_seq = n_seq;
vc->length = length;
vc->sequences = vrna_alloc(sizeof(char *) * (vc->n_seq + 1));
for(s = 0; sequences[s]; s++)
vc->sequences[s] = strdup(sequences[s]);
/* get a copy of the model details */
if(md_p)
md = *md_p;
else /* this fallback relies on global parameters and thus is not threadsafe */
vrna_md_set_default(&md);
aux_options |= WITH_PTYPE;
if(options & VRNA_OPTION_PF)
aux_options |= WITH_PTYPE_COMPAT;
set_fold_compound(vc, &md, options, aux_options);
make_pscores(vc);
if(!(options & VRNA_OPTION_EVAL_ONLY)){
/* add default hard constraints */
vrna_hc_init(vc);
/* add DP matrices (if required) */
vrna_mx_add(vc, VRNA_MX_DEFAULT, options);
}
return vc;
}
PUBLIC vrna_fold_compound_t *
vrna_fold_compound_TwoD(const char *sequence,
const char *s1,
const char *s2,
vrna_md_t *md_p,
unsigned int options){
int length, l, turn;
vrna_fold_compound_t *vc;
vrna_md_t md;
if(sequence == NULL) return NULL;
/* sanity check */
length = strlen(sequence);
if(length == 0)
vrna_message_error("vrna_fold_compound_TwoD: sequence length must be greater 0");
else if(length > vrna_sequence_length_max(options))
vrna_message_error("vrna_fold_compound_TwoD@data_structures.c: sequence length of %d exceeds addressable range", length);
l = strlen(s1);
if(l != length)
vrna_message_error("vrna_fold_compound_TwoD: sequence and s1 differ in length");
l = strlen(s2);
if(l != length)
vrna_message_error("vrna_fold_compound_TwoD: sequence and s2 differ in length");
vc = vrna_alloc(sizeof(vrna_fold_compound_t));
vc->type = VRNA_FC_TYPE_SINGLE;
vc->length = length;
vc->sequence = strdup(sequence);
/* get a copy of the model details */
if(md_p)
md = *md_p;
else /* this fallback relies on global parameters and thus is not threadsafe */
vrna_md_set_default(&md);
/* always make uniq ML decomposition ! */
md.uniq_ML = 1;
md.compute_bpp = 0;
set_fold_compound(vc, &md, options, WITH_PTYPE | WITH_PTYPE_COMPAT);
if(!(options & VRNA_OPTION_EVAL_ONLY)){
vrna_hc_init(vc); /* add default hard constraints */
/* add DP matrices */
vrna_mx_add(vc, VRNA_MX_2DFOLD, options);
}
/* set all fields that are unique to Distance class partitioning... */
turn = vc->params->model_details.min_loop_size;
vc->reference_pt1 = vrna_ptable(s1);
vc->reference_pt2 = vrna_ptable(s2);
vc->referenceBPs1 = vrna_refBPcnt_matrix(vc->reference_pt1, turn);
vc->referenceBPs2 = vrna_refBPcnt_matrix(vc->reference_pt2, turn);
vc->bpdist = vrna_refBPdist_matrix(vc->reference_pt1, vc->reference_pt2, turn);
/* compute maximum matching with reference structure 1 disallowed */
vc->mm1 = maximumMatchingConstraint(vc->sequence, vc->reference_pt1);
/* compute maximum matching with reference structure 2 disallowed */
vc->mm2 = maximumMatchingConstraint(vc->sequence, vc->reference_pt2);
vc->maxD1 = vc->mm1[vc->iindx[1]-length] + vc->referenceBPs1[vc->iindx[1]-length];
vc->maxD2 = vc->mm2[vc->iindx[1]-length] + vc->referenceBPs2[vc->iindx[1]-length];
return vc;
}
PUBLIC void
vrna_fold_compound_add_auxdata( vrna_fold_compound_t *vc,
void *data,
vrna_callback_free_auxdata *f){
if(vc && data){
if(vc->free_auxdata) /* free pre-existing auxdata */
vc->free_auxdata(vc->auxdata);
vc->auxdata = data;
vc->free_auxdata = f;
}
}
PUBLIC void
vrna_fold_compound_add_callback(vrna_fold_compound_t *vc,
vrna_callback_recursion_status *f){
if(vc && f){
vc->stat_cb = f;
}
}
PUBLIC int
vrna_fold_compound_prepare( vrna_fold_compound_t *vc,
unsigned int options){
int ret = 1; /* success */
/* check maximum sequence length restrictions */
if(vc->length > vrna_sequence_length_max(options)){
vrna_message_warning("vrna_fold_compound_prepare@data_structures.c: sequence length of %d exceeds addressable range", vc->length);
return 0;
}
if(options & VRNA_OPTION_MFE){ /* prepare for MFE computation */
switch(vc->type){
case VRNA_FC_TYPE_SINGLE: if(!vc->ptype)
if(!(options & VRNA_OPTION_WINDOW))
vc->ptype = vrna_ptypes(vc->sequence_encoding2,
&(vc->params->model_details));
break;
case VRNA_FC_TYPE_COMPARATIVE: break;
default: break;
}
}
if(options & VRNA_OPTION_PF){ /* prepare for partition function computation */
switch(vc->type){
case VRNA_FC_TYPE_SINGLE: /* get pre-computed Boltzmann factors if not present*/
if(!vc->exp_params)
vc->exp_params = vrna_exp_params(&(vc->params->model_details));
if(!vc->ptype)
vc->ptype = vrna_ptypes(vc->sequence_encoding2, &(vc->exp_params->model_details));
#ifdef VRNA_BACKWARD_COMPAT
/* backward compatibility ptypes */
if(!vc->ptype_pf_compat)
vc->ptype_pf_compat = get_ptypes(vc->sequence_encoding2, &(vc->exp_params->model_details), 1);
#endif
/* get precomputed Boltzmann factors for soft-constraints (if any) */
if(vc->sc){
if(!vc->sc->exp_energy_up)
vrna_sc_set_up(vc, NULL, VRNA_OPTION_PF);
if(!vc->sc->exp_energy_bp)
vrna_sc_set_bp(vc, NULL, VRNA_OPTION_PF);
if(!vc->sc->exp_energy_stack)
vrna_sc_add_SHAPE_deigan(vc, NULL, 0, 0, VRNA_OPTION_PF);
}
if(vc->domains_up) /* turn on unique ML decomposition with qm1 array */
vc->exp_params->model_details.uniq_ML = 1;
break;
case VRNA_FC_TYPE_COMPARATIVE: /* get pre-computed Boltzmann factors if not present*/
if(!vc->exp_params)
vc->exp_params = vrna_exp_params_comparative(vc->n_seq, &(vc->params->model_details));
break;
default: break;
}
}
/* Add DP matrices, if not they are not present or do not fit current settings */
vrna_mx_prepare(vc, options);
return ret;
}
#ifndef VRNA_DISABLE_C11_FEATURES
PUBLIC void
vrna_C11_features(void){
__asm("nop");
}
#endif
/*
#####################################
# BEGIN OF STATIC HELPER FUNCTIONS #
#####################################
*/
PRIVATE void
add_params( vrna_fold_compound_t *vc,
vrna_md_t *md_p,
unsigned int options){
/* ALWAYS add regular energy parameters */
vc->params = vrna_params(md_p);
if(options & VRNA_OPTION_PF){
vc->exp_params = (vc->type == VRNA_FC_TYPE_SINGLE) ? \
vrna_exp_params(md_p) : \
vrna_exp_params_comparative(vc->n_seq, md_p);
}
}
PRIVATE void
set_fold_compound(vrna_fold_compound_t *vc,
vrna_md_t *md_p,
unsigned int options,
unsigned int aux){
char *sequence, **sequences;
unsigned int length, s, i;
int cp; /* cut point for cofold */
char *seq, *seq2;
sequence = NULL;
sequences = NULL;
cp = -1;
/* some default init values */
vc->params = NULL;
vc->exp_params = NULL;
vc->matrices = NULL;
vc->exp_matrices = NULL;
vc->hc = NULL;
vc->auxdata = NULL;
vc->free_auxdata = NULL;
vc->strand_number = NULL;
vc->domains_struc = NULL;
vc->domains_up = NULL;
vc->aux_grammar = NULL;
switch(vc->type){
case VRNA_FC_TYPE_SINGLE: sequence = vc->sequence;
seq2 = strdup(sequence);
seq = vrna_cut_point_remove(seq2, &cp); /* splice out the '&' if concatenated sequences and
reset cp... this should also be safe for
single sequences */
vc->cutpoint = cp;
if((cp > 0) && (md_p->min_loop_size == TURN))
md_p->min_loop_size = 0; /* is it safe to set this here? */
free(vc->sequence);
vc->sequence = seq;
vc->length = length = strlen(seq);
vc->sequence_encoding = vrna_seq_encode(seq, md_p);
vc->sequence_encoding2 = vrna_seq_encode_simple(seq, md_p);
vc->strand_number = (unsigned int *)vrna_alloc(sizeof(unsigned int) * (vc->length + 1));
if (cp > 0)
for (s = i = 0; i <= vc->length; i++) {
if (i == vc->cutpoint)
s++;
vc->strand_number[i] = s;
}
if(!(options & VRNA_OPTION_EVAL_ONLY)){
vc->ptype = (aux & WITH_PTYPE) ? vrna_ptypes(vc->sequence_encoding2, md_p) : NULL;
/* backward compatibility ptypes */
vc->ptype_pf_compat = (aux & WITH_PTYPE_COMPAT) ? get_ptypes(vc->sequence_encoding2, md_p, 1) : NULL;
} else {
vc->ptype = NULL;
vc->ptype_pf_compat = NULL;
}
vc->sc = NULL;
free(seq2);
break;
case VRNA_FC_TYPE_COMPARATIVE: sequences = vc->sequences;
vc->length = length = vc->length;
vc->strand_number = (unsigned int *)vrna_alloc(sizeof(unsigned int) * (vc->length + 1));
vc->cons_seq = consensus((const char **)sequences);
vc->S_cons = vrna_seq_encode_simple(vc->cons_seq, md_p);
vc->pscore = vrna_alloc(sizeof(int)*((length*(length+1))/2+2));
/* backward compatibility ptypes */
vc->pscore_pf_compat = (aux & WITH_PTYPE_COMPAT) ? vrna_alloc(sizeof(int)*((length*(length+1))/2+2)) : NULL;
oldAliEn = vc->oldAliEn = md_p->oldAliEn;
vc->S = vrna_alloc((vc->n_seq+1) * sizeof(short *));
vc->S5 = vrna_alloc((vc->n_seq+1) * sizeof(short *));
vc->S3 = vrna_alloc((vc->n_seq+1) * sizeof(short *));
vc->a2s = vrna_alloc((vc->n_seq+1) * sizeof(unsigned short *));
vc->Ss = vrna_alloc((vc->n_seq+1) * sizeof(char *));
for (s = 0; s < vc->n_seq; s++) {
vrna_aln_encode(vc->sequences[s],
&(vc->S[s]),
&(vc->S5[s]),
&(vc->S3[s]),
&(vc->Ss[s]),
&(vc->a2s[s]),
md_p);
}
vc->S5[vc->n_seq] = NULL;
vc->S3[vc->n_seq] = NULL;
vc->a2s[vc->n_seq] = NULL;
vc->Ss[vc->n_seq] = NULL;
vc->S[vc->n_seq] = NULL;
vc->scs = NULL;
break;
default: /* do nothing ? */
break;
}
if(vc->length <= vrna_sequence_length_max(options)){
vc->iindx = vrna_idx_row_wise(vc->length);
vc->jindx = vrna_idx_col_wise(vc->length);
} else {
vc->iindx = NULL;
vc->jindx = NULL;
}
/* now come the energy parameters */
add_params(vc, md_p, options);
}
PRIVATE void
make_pscores(vrna_fold_compound_t *vc){
/* calculate co-variance bonus for each pair depending on */
/* compensatory/consistent mutations and incompatible seqs */
/* should be 0 for conserved pairs, >0 for good pairs */
#define NONE -10000 /* score for forbidden pairs */
char *structure = NULL;
int i,j,k,l,s, max_span, turn;
float **dm;
int olddm[7][7]={{0,0,0,0,0,0,0}, /* hamming distance between pairs */
{0,0,2,2,1,2,2} /* CG */,
{0,2,0,1,2,2,2} /* GC */,
{0,2,1,0,2,1,2} /* GU */,
{0,1,2,2,0,2,1} /* UG */,
{0,2,2,1,2,0,2} /* AU */,
{0,2,2,2,1,2,0} /* UA */};
short **S = vc->S;
char **AS = vc->sequences;
int n_seq = vc->n_seq;
vrna_md_t *md = (vc->params) ? &(vc->params->model_details) : &(vc->exp_params->model_details);
int *pscore = vc->pscore; /* precomputed array of pair types */
int *indx = vc->jindx;
int *my_iindx = vc->iindx;
int n = vc->length;
turn = md->min_loop_size;
if (md->ribo) {
if (RibosumFile !=NULL) dm=readribosum(RibosumFile);
else dm=get_ribosum((const char **)AS, n_seq, n);
}
else { /*use usual matrix*/
dm = vrna_alloc(7*sizeof(float*));
for (i=0; i<7;i++) {
dm[i] = vrna_alloc(7*sizeof(float));
for (j=0; j<7; j++)
dm[i][j] = (float) olddm[i][j];
}
}
max_span = md->max_bp_span;
if((max_span < turn+2) || (max_span > n))
max_span = n;
for (i=1; i<n; i++) {
for (j=i+1; (j<i+turn+1) && (j<=n); j++)
pscore[indx[j]+i] = NONE;
for (j=i+turn+1; j<=n; j++) {
int pfreq[8]={0,0,0,0,0,0,0,0};
double score;
for (s=0; s<n_seq; s++) {
int type;
if (S[s][i]==0 && S[s][j]==0) type = 7; /* gap-gap */
else {
if ((AS[s][i] == '~')||(AS[s][j] == '~')) type = 7;
else type = md->pair[S[s][i]][S[s][j]];
}
pfreq[type]++;
}
if (pfreq[0]*2+pfreq[7]>n_seq) { pscore[indx[j]+i] = NONE; continue;}
for (k=1,score=0; k<=6; k++) /* ignore pairtype 7 (gap-gap) */
for (l=k; l<=6; l++)
score += pfreq[k]*pfreq[l]*dm[k][l];
/* counter examples score -1, gap-gap scores -0.25 */
pscore[indx[j]+i] = md->cv_fact *
((UNIT*score)/n_seq - md->nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));
if((j - i + 1) > max_span){
pscore[indx[j]+i] = NONE;
}
}
}
if (md->noLP) /* remove unwanted pairs */
for (k=1; k<n-turn-1; k++)
for (l=1; l<=2; l++) {
int type,ntype=0,otype=0;
i=k; j = i+turn+l;
type = pscore[indx[j]+i];
while ((i>=1)&&(j<=n)) {
if ((i>1)&&(j<n)) ntype = pscore[indx[j+1]+i-1];
if ((otype<md->cv_fact*MINPSCORE)&&(ntype<md->cv_fact*MINPSCORE)) /* too many counterexamples */
pscore[indx[j]+i] = NONE; /* i.j can only form isolated pairs */
otype = type;
type = ntype;
i--; j++;
}
}
if (fold_constrained&&(structure!=NULL)) {
int psij, hx, hx2, *stack, *stack2;
stack = vrna_alloc(sizeof(int)*(n+1));
stack2 = vrna_alloc(sizeof(int)*(n+1));
for(hx=hx2=0, j=1; j<=n; j++) {
switch (structure[j-1]) {
case 'x': /* can't pair */
for (l=1; l<j-turn; l++) pscore[indx[j]+l] = NONE;
for (l=j+turn+1; l<=n; l++) pscore[indx[l]+j] = NONE;
break;
case '(':
stack[hx++]=j;
/* fallthrough */
case '[':
stack2[hx2++]=j;
/* fallthrough */
case '<': /* pairs upstream */
for (l=1; l<j-turn; l++) pscore[indx[j]+l] = NONE;
break;
case ']':
if (hx2<=0) {
vrna_message_error("unbalanced brackets in constraints\n%s", structure);
}
i = stack2[--hx2];
pscore[indx[j]+i]=NONE;
break;
case ')':
if (hx<=0) {
vrna_message_error("unbalanced brackets in constraints\n%s", structure);
}
i = stack[--hx];
psij = pscore[indx[j]+i]; /* store for later */
for (k=j; k<=n; k++)
for (l=i; l<=j; l++)
pscore[indx[k]+l] = NONE;
for (l=i; l<=j; l++)
for (k=1; k<=i; k++)
pscore[indx[l]+k] = NONE;
for (k=i+1; k<j; k++)
pscore[indx[k]+i] = pscore[indx[j]+k] = NONE;
pscore[indx[j]+i] = (psij>0) ? psij : 0;
/* fallthrough */
case '>': /* pairs downstream */
for (l=j+turn+1; l<=n; l++) pscore[indx[l]+j] = NONE;
break;
}
}
if (hx!=0) {
vrna_message_error("unbalanced brackets in constraint string\n%s", structure);
}
free(stack); free(stack2);
}
/*free dm */
for (i=0; i<7;i++) {
free(dm[i]);
}
free(dm);
/* copy over pscores for backward compatibility */
if(vc->pscore_pf_compat){
for(i = 1; i < n; i++)
for(j = i; j <= n; j++){
vc->pscore_pf_compat[my_iindx[i] - j] = (short)pscore[indx[j] + i];
}
}
}