/** \file unstructured_domains.c **/

/*
                  Unstructured domains

                  This file contains everything necessary to
                  deal with the default implementation for unstructured
                  domains in secondary structures. This feature enables,
                  for instance, ligand binding to unpaired stretches of
                  an RNA secondary structure.

                  c 2016 Ronny Lorenz

                  ViennaRNA package
*/

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <float.h>
#include <math.h>

#include "ViennaRNA/utils.h"
#include "ViennaRNA/alphabet.h"
#include "ViennaRNA/unstructured_domains.h"

/*
#################################
# PRIVATE MACROS                #
#################################
*/

/*
#################################
# GLOBAL VARIABLES              #
#################################
*/

/*
#################################
# PRIVATE VARIABLES/STRUCTS     #
#################################
*/

struct default_outside {
  int         motif_num;
  FLT_OR_DBL  exp_energy;
};

/*
 *  Default data structure for ligand binding to unpaired stretches
 */
struct ligands_up_data_default {

  /*
  **********************************
    pre-computed position-wise
    motif list
  **********************************
  */
  int           n;
  int           **motif_list_ext;
  int           **motif_list_hp;
  int           **motif_list_int;
  int           **motif_list_mb;

  int           *dG;
  FLT_OR_DBL    *exp_dG;
  int           *len;

  /*
  **********************************
    below are DP matrices to store
    the production rule results
  **********************************
  */
  int           *energies_ext;
  int           *energies_hp;
  int           *energies_int;
  int           *energies_mb;
  FLT_OR_DBL    *exp_energies_ext;
  FLT_OR_DBL    *exp_energies_hp;
  FLT_OR_DBL    *exp_energies_int;
  FLT_OR_DBL    *exp_energies_mb;

  /*
  **********************************
    below are lists to store the
    outside partition function for
    each motif starting at each
    position
  **********************************
  */
  unsigned int            *outside_ext_count;
  struct default_outside  **outside_ext;
  unsigned int            *outside_hp_count;
  struct default_outside  **outside_hp;
  unsigned int            *outside_int_count;
  struct default_outside  **outside_int;
  unsigned int            *outside_mb_count;
  struct default_outside  **outside_mb;

  FLT_OR_DBL  (*default_cb[32])(int, int, struct ligands_up_data_default *);
  FLT_OR_DBL  *exp_e_mx[32];
};

/*
#################################
# PRIVATE FUNCTION DECLARATIONS #
#################################
*/

PRIVATE void        remove_ligands_up(vrna_fold_compound_t *vc);
PRIVATE void        init_ligands_up(vrna_fold_compound_t *vc);

PRIVATE void        add_ligand_motif(vrna_fold_compound_t *vc, const char *motif, double motif_en, unsigned int loop_type);
PRIVATE void        remove_default_data(void *d);

/* default implementations for unstructured domains feature */
PRIVATE void        default_prod_rule(vrna_fold_compound_t *vc, void *d);
PRIVATE void        default_exp_prod_rule(vrna_fold_compound_t *vc, void *d);
PRIVATE int         default_energy(vrna_fold_compound_t *vc, int i, int j, unsigned int loop_type, void *d);
PRIVATE FLT_OR_DBL  default_exp_energy( vrna_fold_compound_t *vc, int i, int j, unsigned int loop_type, void *d);
PRIVATE void        default_probs_add(vrna_fold_compound_t *vc, int i, int j, unsigned int loop_type, FLT_OR_DBL exp_energy, void *data);
PRIVATE FLT_OR_DBL  default_probs_get(vrna_fold_compound_t *vc, int i, int j, unsigned int loop_type, int motif, void *data);

/* helper functions for default implementatations of unstructured domains feature */
PRIVATE int         default_energy_ext_motif(int i, int j, struct ligands_up_data_default *data);
PRIVATE int         default_energy_hp_motif(int i, int j, struct ligands_up_data_default *data);
PRIVATE int         default_energy_int_motif(int i, int j, struct ligands_up_data_default *data);
PRIVATE int         default_energy_mb_motif(int i, int j, struct ligands_up_data_default *data);
PRIVATE FLT_OR_DBL  default_exp_energy_ext_motif(int i, int j, struct ligands_up_data_default *data);
PRIVATE FLT_OR_DBL  default_exp_energy_hp_motif(int i, int j, struct ligands_up_data_default *data);
PRIVATE FLT_OR_DBL  default_exp_energy_int_motif(int i, int j, struct ligands_up_data_default *data);
PRIVATE FLT_OR_DBL  default_exp_energy_mb_motif(int i, int j, struct ligands_up_data_default *data);
PRIVATE void        free_default_data_matrices(struct ligands_up_data_default *data);
PRIVATE void        free_default_data_exp_matrices(struct ligands_up_data_default *data);
PRIVATE void        prepare_matrices( vrna_fold_compound_t *vc, struct ligands_up_data_default *data);
PRIVATE void        prepare_exp_matrices( vrna_fold_compound_t *vc, struct ligands_up_data_default *data);
PRIVATE struct ligands_up_data_default *get_default_data(void);
PRIVATE void        prepare_default_data(vrna_fold_compound_t *vc, struct ligands_up_data_default *data);
PRIVATE void        free_default_data(struct ligands_up_data_default *data);
PRIVATE int         *get_motifs(vrna_fold_compound_t *vc, int i, unsigned int loop_type);
PRIVATE void        annotate_ud(vrna_fold_compound_t *vc, int start, int end, char l, vrna_ud_motif_t **list, int *list_size, int *list_pos);


/*
#################################
# BEGIN OF FUNCTION DEFINITIONS #
#################################
*/
PUBLIC void
vrna_ud_remove( vrna_fold_compound_t *vc){

  if(vc && vc->domains_up)
    remove_ligands_up(vc);
}

PUBLIC void
vrna_ud_set_data( vrna_fold_compound_t                *vc,
                          void                        *data,
                          vrna_callback_free_auxdata  *free_cb){

  if(vc){
    /* init if not already present */
    if(!vc->domains_up)
      init_ligands_up(vc);

    /* free previous data if 'free_data' function present */
    if(vc->domains_up->free_data)
      vc->domains_up->free_data(vc->domains_up->data);

    /* set new data and free callback */
    vc->domains_up->free_data = free_cb;
    vc->domains_up->data      = data;
  }
}

PUBLIC void
vrna_ud_set_prod_rule_cb( vrna_fold_compound_t        *vc,
                          vrna_callback_ud_production *pre_cb,
                          vrna_callback_ud_energy     *e_cb){

  if(vc){
    /* init if not already present */
    if(!vc->domains_up)
      init_ligands_up(vc);

    /* set new callback */
    vc->domains_up->prod_cb   = pre_cb;
    vc->domains_up->energy_cb = e_cb;
  }
}

PUBLIC void
vrna_ud_set_exp_prod_rule_cb( vrna_fold_compound_t            *vc,
                              vrna_callback_ud_exp_production *pre_cb,
                              vrna_callback_ud_exp_energy     *exp_e_cb){

  if(vc){
    /* init if not already present */
    if(!vc->domains_up)
      init_ligands_up(vc);

    /* set new callback */
    vc->domains_up->exp_prod_cb   = pre_cb;
    vc->domains_up->exp_energy_cb = exp_e_cb;
  }
}


PUBLIC void
vrna_ud_set_prob_cb(vrna_fold_compound_t        *vc,
                    vrna_callback_ud_probs_add  *setter,
                    vrna_callback_ud_probs_get  *getter){

  if(vc){
    /* init if not already present */
    if(!vc->domains_up)
      init_ligands_up(vc);

    /* set new callback */
    vc->domains_up->probs_add = setter;
    vc->domains_up->probs_get = getter;
  }
}

PUBLIC void
vrna_ud_add_motif(vrna_fold_compound_t  *vc,
                  const char            *motif,
                  double                motif_en,
                  unsigned int          loop_type){

  if(vc){
    if(!vc->domains_up){
      /* set all default callbacks */
      vrna_ud_set_prod_rule_cb(vc, &default_prod_rule, &default_energy);
      vrna_ud_set_exp_prod_rule_cb(vc, &default_exp_prod_rule, &default_exp_energy);
      vrna_ud_set_data(vc, get_default_data(), &remove_default_data);
      vrna_ud_set_prob_cb(vc, &default_probs_add, &default_probs_get);
    }
    add_ligand_motif(vc, motif, motif_en, loop_type);
  }
}


PUBLIC int *
vrna_ud_get_motif_size_at(vrna_fold_compound_t  *vc,
                          int                   i,
                          unsigned int          loop_type){

  if(vc && vc->domains_up){
    int k, l, cnt, *ret, *ptr;

    ret = NULL;
    if((i > 0) && (i <= vc->length)){
      ptr = get_motifs(vc, i, loop_type);
      if(ptr){
        for(k = 0; ptr[k] != -1; k++) /* replace motif number with its size */
          ptr[k] = vc->domains_up->motif_size[ptr[k]];
        /* make the list unique */
        ret     = (int *)vrna_alloc(sizeof(int) * (k + 1));
        ret[0]  = -1;
        cnt     = 0;
        for(k = 0; ptr[k] != -1; k++){
          for(l = 0; l < cnt; l++){
            if(ptr[k] == ret[l])
              break; /* we've already seen this size */
          }
          if(l == cnt){ /* we've not seen this size before */
            ret[cnt]      = ptr[k];
            ret[cnt + 1]  = -1;
            cnt++;
          }
        }
        /* resize ret array */
        ret = (int *)vrna_realloc(ret, sizeof(int) * (cnt + 1));
      }
      free(ptr);
    }
    return ret;
  }
  return NULL;
}

PUBLIC int *
vrna_ud_get_motifs_at(vrna_fold_compound_t  *vc,
                      int                   i,
                      unsigned int          loop_type){

  if(vc && vc->domains_up){
    if((i > 0) && (i <= vc->length)){
      return get_motifs(vc, i, loop_type);
    }
  }
  return NULL;
}

vrna_ud_motif_t *
vrna_ud_detect_motifs(vrna_fold_compound_t  *vc,
                      const char            *structure){

  int list_size, list_pos;
  vrna_ud_motif_t *motif_list;

  motif_list = NULL;

  if(structure && vc->domains_up){
    int   l, start, end;
    char  last, *loops;

    l           = 0;
    list_pos    = 0;
    list_size   = 15;
    motif_list  = (vrna_ud_motif_t *)vrna_alloc(sizeof(vrna_ud_motif_t) * list_size);
    loops       = vrna_db_to_element_string(structure);

    while(l < vc->length){
      /* skip uppercase encodings */
      while(l < vc->length){
        if(islower(loops[l]))
          break;
        l++;
      }

      if(l < vc->length){
        start = 1 + l;
        last  = loops[l];
        while(loops[l++] == last){
          if(l == vc->length)
            break;
        }
        end = l - 1;
        annotate_ud(vc, start, end, last, &motif_list, &list_size, &list_pos);
      }
    }

    motif_list = (vrna_ud_motif_t *)vrna_realloc(motif_list, sizeof(vrna_ud_motif_t) * (list_pos + 1));
    motif_list[list_pos].start = 0;
    motif_list[list_pos].number = -1;
    free(loops);
  }

  return motif_list;
}

/*
#####################################
# BEGIN OF STATIC HELPER FUNCTIONS  #
#####################################
*/
PRIVATE struct ligands_up_data_default *
get_default_data(void){

  struct ligands_up_data_default *data = vrna_alloc(sizeof(struct ligands_up_data_default));
  data->n                 = 0;
  data->motif_list_ext    = NULL;
  data->motif_list_hp     = NULL;
  data->motif_list_int    = NULL;
  data->motif_list_mb     = NULL;
  data->dG                = NULL;
  data->exp_dG            = NULL;
  data->energies_ext      = NULL;
  data->energies_hp       = NULL;
  data->energies_int      = NULL;
  data->energies_mb       = NULL;
  data->exp_energies_ext  = NULL;
  data->exp_energies_hp   = NULL;
  data->exp_energies_int  = NULL;
  data->exp_energies_mb   = NULL;
  data->outside_ext       = NULL;
  data->outside_hp        = NULL;
  data->outside_int       = NULL;
  data->outside_mb        = NULL;
  data->outside_ext_count = NULL;
  data->outside_hp_count  = NULL;
  data->outside_int_count = NULL;
  data->outside_mb_count  = NULL;
  return data;
}

PRIVATE void
remove_ligands_up(vrna_fold_compound_t *vc){

  int i;

  /* free auxiliary data */
  if(vc->domains_up->free_data)
    vc->domains_up->free_data(vc->domains_up->data);

  for( i = 0; i < vc->domains_up->motif_count; i++ ){
    free(vc->domains_up->motif[i]);
  }
  free(vc->domains_up->motif);
  free(vc->domains_up->motif_size);
  free(vc->domains_up->motif_en);
  free(vc->domains_up->motif_type);

  free(vc->domains_up->uniq_motif_size);

  free(vc->domains_up);

  vc->domains_up = NULL;
}

PRIVATE void
init_ligands_up(vrna_fold_compound_t *vc){

  vc->domains_up = (vrna_ud_t *)vrna_alloc(sizeof(vrna_ud_t));

  vc->domains_up->uniq_motif_count  = 0;
  vc->domains_up->uniq_motif_size   = NULL;
  vc->domains_up->motif_count       = 0;
  vc->domains_up->motif             = NULL;
  vc->domains_up->motif_size        = NULL;
  vc->domains_up->motif_en          = NULL;
  vc->domains_up->motif_type        = NULL;
  vc->domains_up->prod_cb           = NULL;
  vc->domains_up->exp_prod_cb       = NULL;
  vc->domains_up->energy_cb         = NULL;
  vc->domains_up->exp_energy_cb     = NULL;
  vc->domains_up->data              = NULL;
  vc->domains_up->free_data         = NULL;
  vc->domains_up->probs_add         = NULL;
  vc->domains_up->probs_get         = NULL;
}

/*
**********************************
  Default implementation for
  ligand binding to unpaired
  stretches follows below
**********************************
*/

PRIVATE void
add_ligand_motif( vrna_fold_compound_t *vc,
                  const char *motif,
                  double motif_en,
                  unsigned int loop_type){

  unsigned int  i, n, same_size;
  vrna_ud_t     *ud;

  n   = (unsigned int)strlen(motif);
  ud  = vc->domains_up;

  /* First, we update the list of unique motif lengths */
  for(same_size = i = 0; i < ud->uniq_motif_count; i++){
    if(ud->uniq_motif_size[i] == n){
      same_size = 1;
      break;
    }
  }

  if(!same_size){
    ud->uniq_motif_size = (unsigned int *)vrna_realloc(ud->uniq_motif_size, sizeof(unsigned int *) * (ud->uniq_motif_count + 1));
    ud->uniq_motif_size[ud->uniq_motif_count] = n;
    ud->uniq_motif_count++;
  }

  /* And finally, we store the motif */
  ud->motif                       = (char **)vrna_realloc(ud->motif, sizeof(char *) * (ud->motif_count + 1));
  ud->motif[ud->motif_count]      = strdup(motif);
  ud->motif_size                  = (unsigned int *)vrna_realloc(ud->motif_size, sizeof(unsigned int *) * (ud->motif_count + 1));
  ud->motif_size[ud->motif_count] = n;
  ud->motif_en                    = (double *)vrna_realloc(ud->motif_en, sizeof(double) * (ud->motif_count + 1));
  ud->motif_en[ud->motif_count]   = motif_en;
  ud->motif_type                  = (unsigned int *)vrna_realloc(ud->motif_type, sizeof(double) * (ud->motif_count + 1));
  ud->motif_type[ud->motif_count] = loop_type;
  ud->motif_count++;

}

PRIVATE void
remove_default_data(void *d){

  struct ligands_up_data_default  *data;

  data = (struct ligands_up_data_default *)d;

  free_default_data_matrices(data);
  free_default_data_exp_matrices(data);
  free_default_data(data);

  free(data->dG);
  free(data->exp_dG);
}

PRIVATE void
free_default_data(struct ligands_up_data_default *data){

  int i;
  if(data->motif_list_ext){
    for(i=0; i <= data->n; i++)
      free(data->motif_list_ext[i]);
    free(data->motif_list_ext);
  }
  if(data->motif_list_hp){
    for(i=0; i <= data->n; i++)
      free(data->motif_list_hp[i]);
    free(data->motif_list_hp);
  }
  if(data->motif_list_int){
    for(i=0; i <= data->n; i++)
      free(data->motif_list_int[i]);
    free(data->motif_list_int);
  }
  if(data->motif_list_mb){
    for(i=0; i <= data->n; i++)
      free(data->motif_list_mb[i]);
    free(data->motif_list_mb);
  }

  free(data->len);

}

PRIVATE void
free_default_data_matrices(struct ligands_up_data_default *data){

  /* the following four pointers may point to the same memory */
  if(data->energies_ext){
    /* check whether one of the other b* points to the same memory location */
    if(data->energies_ext == data->energies_hp)
      data->energies_hp = NULL;
    if(data->energies_ext == data->energies_int)
      data->energies_int = NULL;
    if(data->energies_ext == data->energies_mb)
      data->energies_mb = NULL;
    free(data->energies_ext);
    data->energies_ext = NULL;
  }
  if(data->energies_hp){
    /* check whether one of the other b* points to the same memory location */
    if(data->energies_hp == data->energies_int)
      data->energies_int = NULL;
    if(data->energies_hp == data->energies_mb)
      data->energies_mb = NULL;
    free(data->energies_hp);
    data->energies_hp = NULL;
  }
  if(data->energies_int){
    /* check whether one of the other b* points to the same memory location */
    if(data->energies_int == data->energies_mb)
      data->energies_mb = NULL;
    free(data->energies_int);
    data->energies_int = NULL;
  }
  free(data->energies_mb);
  data->energies_mb = NULL;
}

PRIVATE void
free_default_data_exp_matrices(struct ligands_up_data_default *data){

  int i;

  /* the following four pointers may point to the same memory */
  if(data->exp_energies_ext){
    /* check whether one of the other b* points to the same memory location */
    if(data->exp_energies_ext == data->exp_energies_hp)
      data->exp_energies_hp = NULL;
    if(data->exp_energies_ext == data->exp_energies_int)
      data->exp_energies_int = NULL;
    if(data->exp_energies_ext == data->exp_energies_mb)
      data->exp_energies_mb = NULL;
    free(data->exp_energies_ext);
    data->exp_energies_ext = NULL;
  }
  if(data->exp_energies_hp){
    /* check whether one of the other b* points to the same memory location */
    if(data->exp_energies_hp == data->exp_energies_int)
      data->exp_energies_int = NULL;
    if(data->exp_energies_hp == data->exp_energies_mb)
      data->exp_energies_mb = NULL;
    free(data->exp_energies_hp);
    data->exp_energies_hp = NULL;
  }
  if(data->exp_energies_int){
    /* check whether one of the other b* points to the same memory location */
    if(data->exp_energies_int == data->exp_energies_mb)
      data->exp_energies_mb = NULL;
    free(data->exp_energies_int);
    data->exp_energies_int = NULL;
  }
  free(data->exp_energies_mb);
  data->exp_energies_mb = NULL;

  if(data->outside_ext)
    for(i = 0; i <= data->n; i++)
      if(data->outside_ext[i])
        free(data->outside_ext[i]);
  free(data->outside_ext);
  free(data->outside_ext_count);

  if(data->outside_hp)
    for(i = 0; i <= data->n; i++)
      if(data->outside_hp[i])
        free(data->outside_hp[i]);
  free(data->outside_hp);
  free(data->outside_hp_count);

  if(data->outside_int)
    for(i = 0; i <= data->n; i++)
      if(data->outside_int[i])
        free(data->outside_int[i]);
  free(data->outside_int);
  free(data->outside_int_count);

  if(data->outside_mb)
    for(i = 0; i <= data->n; i++)
      if(data->outside_mb[i])
        free(data->outside_mb[i]);
  free(data->outside_mb);
  free(data->outside_mb_count);
}

PRIVATE int *
get_motifs(vrna_fold_compound_t *vc, int i, unsigned int loop_type){

  int               k, j, u, n, *motif_list, cnt, guess;
  char              *sequence;
  vrna_ud_t         *domains_up;

  sequence    = vc->sequence;
  n           = (int)vc->length;
  domains_up  = vc->domains_up;

  cnt         = 0;
  guess       = domains_up->motif_count;
  motif_list  = (int *)vrna_alloc(sizeof(int) * (guess + 1));

  /* collect list of motif numbers we find that start at position i */
  for(k = 0; k < domains_up->motif_count; k++){

    if(!(domains_up->motif_type[k] & loop_type))
      continue;

    j = i + domains_up->motif_size[k] - 1;
    if(j <= n){ /* only consider motif that does not exceed sequence length (does not work for circular RNAs!) */
      for(u = i; u <= j; u++){
        if(!vrna_nucleotide_IUPAC_identity(sequence[u-1], domains_up->motif[k][u-i]))
          break;
      }
      if(u > j) /* got a complete motif match */
        motif_list[cnt++] = k;
    }
  }

  if(cnt == 0){
    free(motif_list);
    return NULL;
  }
  
  motif_list = (int *)vrna_realloc(motif_list, sizeof(int) * (cnt + 1));
  motif_list[cnt] = -1; /* end of list marker */

  return motif_list;
}


static void
annotate_ud(vrna_fold_compound_t *vc,
            int start,
            int end,
            char l,
            vrna_ud_motif_t **list,
            int *list_size,
            int *list_pos){

  int i,j;

  /* get motifs in segment [start,end] */
  for(i = start; i <= end; i++){
    unsigned int type = 0;
    switch(l){
      case 'e': type = VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP;
                break;
      case 'h': type = VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP;
                break;
      case 'i': type = VRNA_UNSTRUCTURED_DOMAIN_INT_LOOP;
                break;
      case 'm': type = VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP;
                break;
    }

    int *m = vrna_ud_get_motifs_at(vc, i, type);
    if(m){
      for(j = 0; m[j] != -1; j++){
        int size = vc->domains_up->motif_size[m[j]];

        if(i + size - 1 <= end){
          if(*list_pos == *list_size){
            *list_size *= 1.2;
            *list = (vrna_ud_motif_t *)vrna_realloc(*list, sizeof(vrna_ud_motif_t) * (*list_size));
          }
          (*list)[*list_pos].start = i;
          (*list)[*list_pos].number = m[j];
          (*list_pos)++;
        }
      }
    }
    free(m);
  }
}


PRIVATE void
prepare_matrices( vrna_fold_compound_t *vc,
                  struct ligands_up_data_default *data){

  int         i,j,k,n,size;
  vrna_ud_t   *domains_up;

  n             = (int)vc->length;
  size          = ((n+1)*(n+2))/2 + 1;
  domains_up    = vc->domains_up;

  free_default_data_matrices(data);

  /* here we save memory by re-using DP matrices */
  unsigned int lt[4] = {  VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP,
                          VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP,
                          VRNA_UNSTRUCTURED_DOMAIN_INT_LOOP,
                          VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP };
  int **m[4], *mx;
  m[0] = &data->energies_ext;
  m[1] = &data->energies_hp;
  m[2] = &data->energies_int;
  m[3] = &data->energies_mb;

  for(i=0; i<4; i++){
    unsigned int *col,*col2;
    if(*(m[i]))
      continue;

    mx      = (int *)vrna_alloc(sizeof(int) * size);
    col     = (unsigned int *)vrna_alloc(sizeof(unsigned int) * domains_up->motif_count);
    col2    = (unsigned int *)vrna_alloc(sizeof(unsigned int) * domains_up->motif_count);
    *(m[i]) = mx;

    for(k = 0; k < domains_up->motif_count; k++)
      col[k] = domains_up->motif_type[k] & lt[i];

    /* check if any of the remaining DP matrices can point to the same location */
    for(j=i+1;j<4;j++){
      for(k = 0; k < domains_up->motif_count; k++){
        col2[k] = domains_up->motif_type[k] & lt[j];
        if(col[k] != col2[k])
          break;
      }
      if(k == domains_up->motif_count){
        *(m[j]) = mx;
      }
    }

    free(col);
    free(col2);
  }
}

PRIVATE void
prepare_exp_matrices( vrna_fold_compound_t *vc,
                      struct ligands_up_data_default *data){

  int         i,j,k,n,size;
  vrna_ud_t   *domains_up;

  n             = (int)vc->length;
  size          = ((n+1)*(n+2))/2 + 1;
  domains_up    = vc->domains_up;

  free_default_data_exp_matrices(data);

  /* here we save memory by re-using DP matrices */
  unsigned int lt[4] = {  VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP,
                          VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP,
                          VRNA_UNSTRUCTURED_DOMAIN_INT_LOOP,
                          VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP };
  FLT_OR_DBL **m[4], *mx;
  m[0] = &data->exp_energies_ext;
  m[1] = &data->exp_energies_hp;
  m[2] = &data->exp_energies_int;
  m[3] = &data->exp_energies_mb;

  for(i=0; i<4; i++){
    unsigned int *col,*col2;
    if(*(m[i]))
      continue;

    mx      = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * size);
    col     = (unsigned int *)vrna_alloc(sizeof(unsigned int) * domains_up->motif_count);
    col2    = (unsigned int *)vrna_alloc(sizeof(unsigned int) * domains_up->motif_count);
    *(m[i]) = mx;

    for(k = 0; k < domains_up->motif_count; k++)
      col[k] = domains_up->motif_type[k] & lt[i];

    /* check if any of the remaining DP matrices can point to the same location */
    for(j=i+1;j<4;j++){
      for(k = 0; k < domains_up->motif_count; k++){
        col2[k] = domains_up->motif_type[k] & lt[j];
        if(col[k] != col2[k])
          break;
      }
      if(k == domains_up->motif_count){
        *(m[j]) = mx;
      }
    }

    free(col);
    free(col2);
  }

  /* now prepate memory for outside partition function */
  data->outside_ext = (struct default_outside **)vrna_alloc(sizeof(struct default_outside *) * (n + 2));
  data->outside_hp  = (struct default_outside **)vrna_alloc(sizeof(struct default_outside *) * (n + 2));
  data->outside_int = (struct default_outside **)vrna_alloc(sizeof(struct default_outside *) * (n + 2));
  data->outside_mb  = (struct default_outside **)vrna_alloc(sizeof(struct default_outside *) * (n + 2));
  data->outside_ext_count = (unsigned int *)vrna_alloc(sizeof(unsigned int) * (n + 2));
  data->outside_hp_count  = (unsigned int *)vrna_alloc(sizeof(unsigned int) * (n + 2));
  data->outside_int_count = (unsigned int *)vrna_alloc(sizeof(unsigned int) * (n + 2));
  data->outside_mb_count  = (unsigned int *)vrna_alloc(sizeof(unsigned int) * (n + 2));
}

PRIVATE void
prepare_default_data( vrna_fold_compound_t *vc,
                      struct ligands_up_data_default *data){

  int         i, n;
  vrna_ud_t   *domains_up;

  n           = (int)vc->length;
  domains_up  = vc->domains_up;

  data->n = n;
  free_default_data(data);

  /*
      create motif_list for associating a nucleotide position with all
      motifs that start there
  */
  data->motif_list_ext  = (int **)vrna_alloc(sizeof(int *) * (n+1));
  data->motif_list_hp   = (int **)vrna_alloc(sizeof(int *) * (n+1));
  data->motif_list_int  = (int **)vrna_alloc(sizeof(int *) * (n+1));
  data->motif_list_mb   = (int **)vrna_alloc(sizeof(int *) * (n+1));
  data->motif_list_ext[0] = NULL;
  data->motif_list_hp[0]  = NULL;
  data->motif_list_int[0] = NULL;
  data->motif_list_mb[0]  = NULL;
  for(i = 1; i <= n; i++){
    data->motif_list_ext[i] = get_motifs(vc, i, VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP);
    data->motif_list_hp[i]  = get_motifs(vc, i, VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP);
    data->motif_list_int[i] = get_motifs(vc, i, VRNA_UNSTRUCTURED_DOMAIN_INT_LOOP);
    data->motif_list_mb[i]  = get_motifs(vc, i, VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP);
  }

  /*  store length of motifs in 'data' */
  data->len = (int *)vrna_alloc(sizeof(int) * domains_up->motif_count);
  for(i = 0; i < domains_up->motif_count; i++)
    data->len[i] = domains_up->motif_size[i];

}

PRIVATE void
default_prod_rule(vrna_fold_compound_t *vc,
                  void *d){

  int                             i,j,k,l,u,n,size,e_ext, e_hp, e_int, e_mb,en,en2,*idx;
  unsigned int                    loop_type;
  vrna_ud_t                       *domains_up;
  struct ligands_up_data_default  *data;

  int           *energies_ext;
  int           *energies_hp;
  int           *energies_int;
  int           *energies_mb;

  n             = (int)vc->length;
  size          = ((n+1)*(n+2))/2 + 1;
  idx           = vc->jindx;
  domains_up    = vc->domains_up;
  data          = (struct ligands_up_data_default *)d;

  prepare_default_data(vc, data);
  prepare_matrices(vc, data);

  energies_ext  = data->energies_ext;
  energies_hp   = data->energies_hp;
  energies_int  = data->energies_int;
  energies_mb   = data->energies_mb;

  data->default_cb[VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP] = default_exp_energy_ext_motif;
  data->default_cb[VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP]  = default_exp_energy_hp_motif;
  data->default_cb[VRNA_UNSTRUCTURED_DOMAIN_INT_LOOP] = default_exp_energy_int_motif;
  data->default_cb[VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP]  = default_exp_energy_mb_motif;

  /*  precompute energy contributions of the motifs */
  data->dG = (int *)vrna_alloc(sizeof(int) * domains_up->motif_count);
  for(i = 0; i < domains_up->motif_count; i++)
    data->dG[i] = (int)roundf(domains_up->motif_en[i] * 100.);

  /* now we can start to fill the DP matrices */
  for(i=n; i>0; i--){
    int *list_ext = data->motif_list_ext[i];
    int *list_hp  = data->motif_list_hp[i];
    int *list_int = data->motif_list_int[i];
    int *list_mb  = data->motif_list_mb[i];
    for(j=i;j<=n;j++){
      if(i < j){
        e_ext = energies_ext[idx[j]+i+1];
        e_hp  = energies_hp[idx[j]+i+1];
        e_int = energies_int[idx[j]+i+1];
        e_mb  = energies_mb[idx[j]+i+1];
      } else {
        e_ext = INF;
        e_hp  = INF;
        e_int = INF;
        e_mb  = INF;
      }
      if(list_ext){
        for(k = 0; -1 != (l = list_ext[k]); k++){
          u   = i + data->len[l] - 1;
          en  = data->dG[l];
          if(u <= j){
            e_ext = MIN2(e_ext, en);
            if(u < j){
              en2 = en + energies_ext[idx[j]+u+1];
              e_ext = MIN2(e_ext, en2);
            }
          }
        }
      }
      if(list_hp){
        for(k = 0; -1 != (l = list_hp[k]); k++){
          u   = i + data->len[l] - 1;
          en  = data->dG[l];
          if(u <= j){
            e_hp = MIN2(e_hp, en);
            if(u < j){
              en2 = en + energies_hp[idx[j]+u+1];
              e_hp = MIN2(e_hp, en2);
            }
          }
        }
      }
      if(list_int){
        for(k = 0; -1 != (l = list_int[k]); k++){
          u   = i + data->len[l] - 1;
          en  = data->dG[l];
          if(u <= j){
            e_int = MIN2(e_int, en);
            if(u < j){
              en2 = en + energies_int[idx[j]+u+1];
              e_int = MIN2(e_int, en2);
            }
          }
        }
      }
      if(list_mb){
        for(k = 0; -1 != (l = list_mb[k]); k++){
          u   = i + data->len[l] - 1;
          en  = data->dG[l];
          if(u <= j){
            e_mb = MIN2(e_mb, en);
            if( u < j){
              en2 = en + energies_mb[idx[j]+u+1];
              e_mb = MIN2(e_mb, en2);
            }
          }
        }
      }

      energies_ext[idx[j]+i]  = e_ext;
      energies_hp[idx[j]+i]   = e_hp;
      energies_int[idx[j]+i]  = e_int;
      energies_mb[idx[j]+i]   = e_mb;
    }
  }
}

PRIVATE void
default_exp_prod_rule(vrna_fold_compound_t *vc,
                      void *d){

  int                             i,j,k,l,u,n,size,*idx;
  unsigned int                    loop_type;
  FLT_OR_DBL                      q_ext, q_hp, q_int, q_mb, q, qq;
  vrna_ud_t                       *domains_up;
  struct ligands_up_data_default  *data;

  FLT_OR_DBL    *exp_energies_ext;
  FLT_OR_DBL    *exp_energies_hp;
  FLT_OR_DBL    *exp_energies_int;
  FLT_OR_DBL    *exp_energies_mb;
  double        kT;

  n             = (int)vc->length;
  size          = ((n+1)*(n+2))/2 + 1;
  idx           = vc->iindx;
  domains_up    = vc->domains_up;
  data          = (struct ligands_up_data_default *)d;
  kT            = vc->exp_params->kT;

  prepare_default_data(vc, data);
  prepare_exp_matrices(vc, data);

  exp_energies_ext  = data->exp_energies_ext;
  exp_energies_hp   = data->exp_energies_hp;
  exp_energies_int  = data->exp_energies_int;
  exp_energies_mb   = data->exp_energies_mb;

  data->default_cb[VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP] = default_exp_energy_ext_motif;
  data->default_cb[VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP]  = default_exp_energy_hp_motif;
  data->default_cb[VRNA_UNSTRUCTURED_DOMAIN_INT_LOOP] = default_exp_energy_int_motif;
  data->default_cb[VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP]  = default_exp_energy_mb_motif;

  data->exp_e_mx[VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP] = data->exp_energies_ext;
  data->exp_e_mx[VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP]  = data->exp_energies_hp;
  data->exp_e_mx[VRNA_UNSTRUCTURED_DOMAIN_INT_LOOP] = data->exp_energies_int;
  data->exp_e_mx[VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP]  = data->exp_energies_mb;

  /*  precompute energy contributions of the motifs */
  data->exp_dG = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * domains_up->motif_count);
  for(i = 0; i < domains_up->motif_count; i++){
    double GT = domains_up->motif_en[i] * 1000.; /* in cal/mol */
    data->exp_dG[i] = (FLT_OR_DBL)exp( -GT / kT);
  }

  /* now we can start to fill the DP matrices */
  for(i=n; i>0; i--){
    int *list_ext = data->motif_list_ext[i];
    int *list_hp  = data->motif_list_hp[i];
    int *list_int = data->motif_list_int[i];
    int *list_mb  = data->motif_list_mb[i];
    for(j = i; j <= n; j++){
      if(i < j){
        q_ext = exp_energies_ext[idx[i + 1] - j];
        q_hp  = exp_energies_hp[idx[i + 1] - j];
        q_int = exp_energies_int[idx[i + 1] - j];
        q_mb  = exp_energies_mb[idx[i + 1] - j];
      } else {
        q_ext = 0;
        q_hp  = 0;
        q_int = 0;
        q_mb  = 0;
      }
      if(list_ext){
        for(k = 0; -1 != (l = list_ext[k]); k++){
          u = i + data->len[l] - 1;
          q = data->exp_dG[l];
          if(u <= j){
            q_ext += q;
            if(u < j)
              q_ext += q * exp_energies_ext[idx[u + 1] - j];
          }
        }
      }
      if(list_hp){
        for(k = 0; -1 != (l = list_hp[k]); k++){
          u   = i + data->len[l] - 1;
          q  = data->exp_dG[l];
          if(u <= j){
            q_hp += q;
            if(u < j)
              q_hp += q * exp_energies_hp[idx[u + 1] - j];
          }
        }
      }
      if(list_int){
        for(k = 0; -1 != (l = list_int[k]); k++){
          u   = i + data->len[l] - 1;
          q  = data->exp_dG[l];
          if(u <= j){
            q_int += q;
            if(u < j)
              q_int += q * exp_energies_int[idx[u + 1] - j];
          }
        }
      }
      if(list_mb){
        for(k = 0; -1 != (l = list_mb[k]); k++){
          u   = i + data->len[l] - 1;
          q  = data->exp_dG[l];
          if(u <= j){
            q_mb += q;
            if(u < j)
              q_mb += q * exp_energies_mb[idx[u + 1] - j];
          }
        }
      }

      exp_energies_ext[idx[i] - j]  = q_ext;
      exp_energies_hp[idx[i] - j]   = q_hp;
      exp_energies_int[idx[i] - j]  = q_int;
      exp_energies_mb[idx[i] - j]   = q_mb;
    }
  }
}

PRIVATE int
default_energy( vrna_fold_compound_t *vc,
                int i,
                int j,
                unsigned int loop_type,
                void *d){

  int en, ij, *idx = vc->jindx;
  struct ligands_up_data_default *data = (struct ligands_up_data_default *)d;

  en = INF;
  ij = idx[j] + i;

  if(j < i)
    return INF;

  if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_MOTIF){
    if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP)
      en = default_energy_ext_motif(i, j, data);
    else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP)
      en = default_energy_hp_motif(i, j, data);
    else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_INT_LOOP)
      en = default_energy_int_motif(i, j, data);
    else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP)
      en = default_energy_mb_motif(i, j, data);
  } else {
    if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP){
      if(data->energies_ext)
        en = data->energies_ext[ij];
    } else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP){
      if(data->energies_hp)
        en = data->energies_hp[ij];
    } else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_INT_LOOP){
      if(data->energies_int)
        en = data->energies_int[ij];
    } else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP){
      if(data->energies_mb)
        en = data->energies_mb[ij];
    }
  }

  return en;
}

PRIVATE FLT_OR_DBL
default_exp_energy( vrna_fold_compound_t *vc,
                    int i,
                    int j,
                    unsigned int loop_type,
                    void *d){

  FLT_OR_DBL                      q;
  int                             ij, *idx;
  struct ligands_up_data_default  *data;

  q     = 0;
  data  = (struct ligands_up_data_default *)d;

  if(j < i)
    return 0.;

  if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_MOTIF){
    q = data->default_cb[loop_type & ~(VRNA_UNSTRUCTURED_DOMAIN_MOTIF)](i, j, data);
  } else {
    idx = vc->iindx;
    ij  = idx[i] - j;
    q   = data->exp_e_mx[loop_type][ij];
  }

  return q;
}

PRIVATE int
default_energy_ext_motif( int i,
                          int j,
                          struct ligands_up_data_default *data){

  int k, m;
  int e = INF;

  if(data->motif_list_ext[i]){
    k = 0;
    while(-1 != (m = data->motif_list_ext[i][k])){
      if((i + data->len[m] - 1) == j)
        e = MIN2(e, data->dG[m]);
      k++;
    }
  }
  return e;
}

PRIVATE int
default_energy_hp_motif(int i,
                        int j,
                        struct ligands_up_data_default *data){

  int k, m;
  int e = INF;

  if(data->motif_list_hp[i]){
    k = 0;
    while(-1 != (m = data->motif_list_hp[i][k])){
      if((i + data->len[m] - 1) == j)
        e = MIN2(e, data->dG[m]);
      k++;
    }
  }
  return e;
}

PRIVATE int
default_energy_int_motif( int i,
                          int j,
                          struct ligands_up_data_default *data){

  int k, m;
  int e = INF;

  if(data->motif_list_int[i]){
    k = 0;
    while(-1 != (m = data->motif_list_int[i][k])){
      if((i + data->len[m] - 1) == j)
        e = MIN2(e, data->dG[m]);
      k++;
    }
  }
  return e;
}

PRIVATE int
default_energy_mb_motif(int i,
                        int j,
                        struct ligands_up_data_default *data){

  int k, m;
  int e = INF;

  if(data->motif_list_mb[i]){
    k = 0;
    while(-1 != (m = data->motif_list_mb[i][k])){
      if((i + data->len[m] - 1) == j)
        e = MIN2(2, data->dG[m]);
      k++;
    }
  }
  return e;
}

PRIVATE FLT_OR_DBL
default_exp_energy_ext_motif( int i,
                              int j,
                              struct ligands_up_data_default *data){

  int         k, m;
  FLT_OR_DBL  q = 0;

  if(data->motif_list_ext[i]){
    k = 0;
    while(-1 != (m = data->motif_list_ext[i][k])){
      if((i + data->len[m] - 1) == j)
        q += data->exp_dG[m];
      k++;
    }
  }

  return q;
}

PRIVATE FLT_OR_DBL
default_exp_energy_hp_motif(int i,
                            int j,
                            struct ligands_up_data_default *data){

  int         k, m;
  FLT_OR_DBL  q = 0;

  if(data->motif_list_hp[i]){
    k = 0;
    while(-1 != (m = data->motif_list_hp[i][k])){
      if((i + data->len[m] - 1) == j)
        q += data->exp_dG[m];
      k++;
    }
  }

  return q;
}

PRIVATE FLT_OR_DBL
default_exp_energy_int_motif( int i,
                              int j,
                              struct ligands_up_data_default *data){

  int         k, m;
  FLT_OR_DBL  q = 0;

  if(data->motif_list_int[i]){
    k = 0;
    while(-1 != (m = data->motif_list_int[i][k])){
      if((i + data->len[m] - 1) == j)
        q += data->exp_dG[m];
      k++;
    }
  }

  return q;
}

PRIVATE FLT_OR_DBL
default_exp_energy_mb_motif(int i,
                            int j,
                            struct ligands_up_data_default *data){

  int         k, m;
  FLT_OR_DBL  q = 0;

  if(data->motif_list_mb[i]){
    k = 0;
    while(-1 != (m = data->motif_list_mb[i][k])){
      if((i + data->len[m] - 1) == j)
        q += data->exp_dG[m];
      k++;
    }
  }

  return q;
}

PRIVATE void
default_probs_add(vrna_fold_compound_t *vc,
                  int i,
                  int j,
                  unsigned int loop_type,
                  FLT_OR_DBL exp_energy,
                  void *data){

  int                             **motif_list, k, l, m;
  unsigned int                    *size, *cnt, o;
  struct ligands_up_data_default  *d;
  struct default_outside          **storage, **st;

  d = (struct ligands_up_data_default  *)data;

  if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_MOTIF){
    if(j < i)
      return;

    if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP){
      motif_list  = d->motif_list_ext;
      storage     = &(d->outside_ext[i]);
      size        = &(d->outside_ext_count[i]);
    } else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP) {
      motif_list = d->motif_list_hp;
      storage     = &(d->outside_hp[i]);
      size        = &(d->outside_hp_count[i]);
    } else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_INT_LOOP){
      motif_list = d->motif_list_int;
      storage     = &(d->outside_int[i]);
      size        = &(d->outside_int_count[i]);
    } else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP){
      motif_list = d->motif_list_mb;
      storage     = &(d->outside_mb[i]);
      size        = &(d->outside_mb_count[i]);
    } else{
      vrna_message_warning("Unknown unstructured domain loop type");
      return;
    }

    k = 0;
    while(-1 != (m = motif_list[i][k])){
      if((i + d->len[m] - 1) == j){

        /* check for addition first */
        for(o = 0; o < *size; o++)
          if((*storage)[o].motif_num == m){ /* found previously added motif constribution */
            (*storage)[o].exp_energy += exp_energy;
            break;
          }

        /* if we haven't added yet, create new list entry */
        if(o == *size){
          *storage = (struct default_outside *)vrna_realloc(*storage, sizeof(struct default_outside) * (*size + 1));
          (*storage)[*size].motif_num   = m;
          (*storage)[*size].exp_energy  = exp_energy;
          (*size)++;
        }
      }
      k++;
    }

  } else {
    if(j < i)
      return;

    FLT_OR_DBL pf, exp_e;
    pf = default_exp_energy(vc, i, j, loop_type, data);

    if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP){
      motif_list  = d->motif_list_ext;
      storage     = d->outside_ext;
      size        = d->outside_ext_count;
    } else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP) {
      motif_list = d->motif_list_hp;
      storage     = d->outside_hp;
      size        = d->outside_hp_count;
    } else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_INT_LOOP){
      motif_list = d->motif_list_int;
      storage     = d->outside_int;
      size        = d->outside_int_count;
    } else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP){
      motif_list = d->motif_list_mb;
      storage     = d->outside_mb;
      size        = d->outside_mb_count;
    } else{
      vrna_message_warning("Unknown unstructured domain loop type");
      return;
    }

    /* check for each motif starting at any k with i <= k <= j */
    for(k = i; k <= j; k++){
      if(motif_list[k]){
        st  = &(storage[k]);
        cnt = &(size[k]);
        for(l = 0; motif_list[k][l] != -1; l++){
          m = motif_list[k][l];

          if(j < k + d->len[m] - 1) /* motifs must be sorted be length */
            continue;

          exp_e         = d->exp_dG[m];
          FLT_OR_DBL p  = exp_e / pf;

          /* add/insert contribution */

          /* check for addition first */
          for(o = 0; o < *cnt; o++)
            if((*st)[o].motif_num == m){ /* found previously added motif constribution */
              (*st)[o].exp_energy += p * exp_energy;
              break;
            }

          /* if we haven't added yet, create new list entry */
          if(o == *cnt){
            *st = (struct default_outside *)vrna_realloc(*st, sizeof(struct default_outside) * (*cnt + 1));
            (*st)[*cnt].motif_num   = m;
            (*st)[*cnt].exp_energy  = p * exp_energy;
            (*cnt)++;
          }
        }
      }
    }
  }

}

PRIVATE FLT_OR_DBL
default_probs_get(vrna_fold_compound_t *vc,
                  int i,
                  int j,
                  unsigned int loop_type,
                  int motif,
                  void *data){

  FLT_OR_DBL                      outside = 0.;
  unsigned int                    *size, k;
  struct ligands_up_data_default  *d;
  struct default_outside          **storage;

  d = (struct ligands_up_data_default  *)data;

  if(j < i)
    return 0.;

  if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP){
    storage     = &(d->outside_ext[i]);
    size        = &(d->outside_ext_count[i]);
  } else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP) {
    storage     = &(d->outside_hp[i]);
    size        = &(d->outside_hp_count[i]);
  } else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_INT_LOOP){
    storage     = &(d->outside_int[i]);
    size        = &(d->outside_int_count[i]);
  } else if(loop_type & VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP){
    storage     = &(d->outside_mb[i]);
    size        = &(d->outside_mb_count[i]);
  } else{
    vrna_message_warning("Unknown unstructured domain loop type");
    return 0.;
  }

  for(k = 0; k < *size; k++){
    /* check for motif number match */
    if((*storage)[k].motif_num == motif)
      /* check for length match */
      if(i + d->len[motif] - 1 == j){
        outside += (*storage)[k].exp_energy;
      }
  }

  return outside;
}