bddem.c

/*

EMBLEM and SLIPCASE

Copyright (c) 2013, Fabrizio Riguzzi and Elena Bellodi

This package uses the library cudd, see http://vlsi.colorado.edu/~fabio/CUDD/
for the relative license.

*/
#include <math.h>
#include <float.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <locale.h>
#include "cudd.h"
#include "st.h"
#include <SWI-Prolog.h>
#include <unistd.h>
#include <sys/types.h>

#ifdef _WIN32
#include <Windows.h>
#endif

#define BUFSIZE 200000
#define LOGZERO log(0.01)
#define CACHE_SLOTS 1
#define UNIQUE_SLOTS 1
#define RETURN_IF_FAIL if (ret!=TRUE) return ret;
#define DEBUG 0


typedef struct
{
  int nVal,nRule;
  int firstBoolVar;
  int abducible;
  int query;
  int decision;
} variable;


typedef struct
{
  DdNode *key;
  double value;
} rowel;

typedef struct
{
  int cnt;
  rowel *row;
} tablerow;


typedef struct environment
{
  DdManager * mgr; //Cudd manager
  int * bVar2mVar; //array that maps Boolean vars to multi-valued vars
  variable * vars; // multivalued variables
  int nVars;  // number of multivalued variables
  double * probs; // probabilities of Boolean variables
  int  boolVars;  // number of Boolean variables
  int nRules;  // number of rules
  int * rules; // array with the number of head atoms for each rule
  int n_abd;
  int n_abd_boolVars;
} environment;

typedef struct
{
  environment * env; // one environment for each example
  int ex;  // number of examples
  double * sigma; // sigma array for taking into account deleted paths
  double ***eta;  // eta array: for each rule, each Bool var stores two doubles
  double ***eta_temp; // eta array storing the contribution of the current example
  int * rules; // array with the number of head atoms for each rule
  int * tunable_rules; // array with 1 if the parameters of the rule are tunable, 0 otherwise
  int nRules; // number of rules
  double **arrayprob; //value of paramters. One value ofr each rule and Bool var
  double * nodes_probs;
  tablerow * nodesB; // tables of probabilities for nodes in Backward step
  tablerow * nodesFE; // tables of probabilities for nodes in Forward step
  tablerow * nodesFO; // tables of probabilities for nodes in Forward step
  double * example_prob; // probability (frequency) of examples in the data
  double alpha; // type of parameter initialization in EM: 
                // 0 for truncated Dirichlet process
                // >0 for symmetric Dirichlet distribution with values alpha
} example_data;

typedef struct
{
  int var,val;
} assign;

typedef struct explan
{
  assign a; // for abduction, an int is sufficient, since there are only 1 in a->val
  struct explan * next;
} explan_t;

typedef struct
{
  double prob;
  explan_t * mpa;
} prob_abd_expl;

typedef struct
{
  int n_elements; // quanti elementi ci sono nella lista
  double prob;
  explan_t **mpa; // lista di liste
} prob_abd_expl_list;


typedef struct
{
  DdNode *node;
  int comp;
} explkey;

typedef struct
{
  explkey key;
  prob_abd_expl value;
} explrowel;

typedef struct
{
  explkey key;
  prob_abd_expl_list value;
} explrowel_abd;

typedef struct
{
  int cnt;
  explrowel *row;
} expltablerow;

typedef struct
{
  int cnt;
  explrowel_abd *row;
} expltablerow_abd;

typedef struct {
  double impact;
  DdNode *term;
} term_entry;


static foreign_t ret_prob(term_t,term_t,term_t);
static foreign_t ret_abd_prob(term_t,term_t,term_t,term_t);
static foreign_t ret_map_prob(term_t,term_t,term_t,term_t);

double Prob(DdNode *node,environment *env,tablerow *table);
// prob_abd_expl abd_Prob(DdNode *node,environment *env,expltablerow expltable, tablerow *table, int comp_par);
prob_abd_expl map_Prob(DdNode *node, environment * env,
    expltablerow * maptable, tablerow * table,
    int comp_par);

static foreign_t end_ex(term_t);
static foreign_t init(term_t);
static foreign_t end(term_t arg1);
static foreign_t init_ex(term_t arg1, term_t arg2);
static foreign_t add_var(term_t,term_t,term_t,term_t);
static foreign_t add_query_var(term_t,term_t,term_t,term_t);
static foreign_t add_abd_var(term_t,term_t,term_t,term_t);
static foreign_t init_em(term_t);
static foreign_t end_em(term_t);
static foreign_t EM(term_t,term_t,term_t,term_t,
  term_t,term_t,term_t,term_t,term_t);
static foreign_t reorder(term_t arg1);
static foreign_t make_query_var(term_t arg1, term_t arg2, term_t arg3);


static foreign_t init_par(example_data * ex_d, term_t ruleHeadsTerm);
double ProbPath(example_data * ex_d,DdNode *node, int nex);
//static int rec_deref(void);
void Forward(example_data * ex_d,DdNode *node, int nex);
void UpdateForward(example_data * ex_d,DdNode * node, int nex,
  DdNode *** nodesToVisit,int * NnodesToVisit);
double GetOutsideExpe(example_data *ex_d,DdNode *root,double ex_prob, int nex);
void Maximization(example_data * ex_d);
static double Expectation(example_data *ex_d,DdNode **nodes_ex, int lenNodes);
int reorder_int(environment *env);

FILE *open_file(char *filename, const char *mode);
tablerow* init_table(int varcnt);
double * get_value(tablerow *tab,  DdNode *node);
void add_or_replace_node(tablerow *tab, DdNode *node, double value);
void add_node(tablerow *tab, DdNode *node, double value);
void destroy_table(tablerow *tab,int varcnt);
expltablerow* expl_init_table(int varcnt);
prob_abd_expl * expl_get_value(expltablerow *tab,  DdNode *node, int comp);
void expl_add_node(expltablerow *tab, DdNode *node, int comp, prob_abd_expl value);
void expl_destroy_table(expltablerow *tab,int varcnt);
DdNode* get_node(DdNode *node,tablerow *tab);

install_t install(void);
void write_dot(environment * env, DdNode * bdd, FILE * file);

explan_t * insert(assign assignment,explan_t * head);
explan_t * duplicate(explan_t * head);
void free_list(explan_t * head);

term_t clist_to_pllist(explan_t *mpa, environment * env);
term_t abd_clist_to_pllist(explan_t *mpa);


double uniform_sample();
double gauss_sample(double mean,double var);
double gamma_sample(double shape, double scale);
double gamma_sample_gt1(double shape);
void dirichlet_sample(double * alpha,int K, double * theta);
void symmetric_dirichlet_sample(double alpha,int K, double * theta);

static foreign_t gamma_sample_pl(term_t arg1,term_t arg2,term_t arg3);
static foreign_t gauss_sample_pl(term_t arg1,term_t arg2,term_t arg3);
static foreign_t uniform_sample_pl(term_t arg1);
static foreign_t dirichlet_sample_pl(term_t arg1,term_t arg2);
static foreign_t symmetric_dirichlet_sample_pl(term_t arg1,term_t arg2, term_t arg3);
static foreign_t discrete_sample_pl(term_t arg1,term_t arg2);
static foreign_t initial_values_pl(term_t arg1, term_t arg2);

DdNode* probability_dd_int(environment *env, DdNode *current_node);
static DdNode* probability_dd_int_impl(environment *env, DdNode *current_node);
int find_path(DdNode *node, double value, int **array, int *len);
void debug_cudd_env(environment *env, int i);
void dump_var(variable *var);
void dump_env(environment *env);
DdNode *setLowerBound(DdManager *dd, DdNode *current_node, double lowerBound, double floor_value);
static DdNode *setLowerBound_impl(DdManager *dd, DdNode *current_node, double lowerBound, double floor_value);
static foreign_t add_decision_var(term_t env_ref, term_t rule, term_t nvals, term_t var_out);
static foreign_t probability_dd(term_t env_ref, term_t bdd_ref, term_t add_out);
static foreign_t add_prod(term_t env_ref, term_t add_in, term_t value, term_t add_out);
static foreign_t add_sum(term_t env_ref, term_t add_A, term_t add_B, term_t add_out);
static foreign_t ret_strategy(term_t env_ref, term_t add_A, term_t strategy_list, term_t util);
static foreign_t compute_optimal_strategy_sp(term_t env_ref, term_t b_list, term_t u_list, term_t strategy_list, term_t util);
static term_t debug_cudd_var(term_t env_ref, term_t out_null);
void print_prob_abd_expl(prob_abd_expl *pae);
void print_abd_explan(explan_t *et);
int in_list(explan_t *list, assign element);
explan_t *merge_explain(explan_t *a, explan_t *b);

static DdNode *add_const_ref(DdManager *mgr, double value);
static DdNode *add_apply_ref(DdManager *mgr, DdNode *a, DdNode *b, DD_AOP op);
static DdNode *add_ite_ref(DdManager *mgr, DdNode *f, DdNode *g, DdNode *h);

typedef struct {
  int was_on;
  Cudd_ReorderingType method;
} autodyn_guard_t;

/* Disable Autodyn around code that recursively traverses existing nodes.
   If reordering is active there, child pointers captured earlier can become stale. */
static void autodyn_guard_begin(DdManager *mgr, autodyn_guard_t *guard) {
  guard->method = CUDD_REORDER_NONE;
  guard->was_on = Cudd_ReorderingStatus(mgr, &guard->method);
  if (guard->was_on) {
    Cudd_AutodynDisable(mgr);
  }
}

static void autodyn_guard_end(DdManager *mgr, const autodyn_guard_t *guard) {
  if (guard->was_on) {
    Cudd_AutodynEnable(mgr, guard->method);
  }
}

prob_abd_expl_list abd_Prob(DdNode *node, environment *env, expltablerow_abd *expltable, tablerow *table, int comp_par); 
void print_prob_abd_expl_list(const prob_abd_expl_list *pael);
void expl_destroy_table_abd(expltablerow_abd *tab,int varcnt);
prob_abd_expl_list * expl_get_value_abd(expltablerow_abd *tab,  DdNode *node, int comp);
expltablerow_abd *expl_init_table_abd(int varcnt);
explan_t ** insert_abd(assign assignment,explan_t ** head, int n_elements);
void expl_add_node_abd(expltablerow_abd *tab, DdNode *node, int comp, prob_abd_expl_list value);
explan_t **split_explan(assign assignment, explan_t **head_true, explan_t **head_false, int n_true, int n_false, int *n_new_elements,  int no_false);

explan_t *insord_abd(assign assignment, explan_t *head);
explan_t **insert_abd_ordered(assign assignment, explan_t **head, int n_elements);
explan_t **split_explan_v2(assign assignment, explan_t **head_true, explan_t **head_false, int *n_true, int *n_false, int *n_new_elements);
int is_subset(explan_t *e01, explan_t *e02);

static foreign_t uniform_sample_pl(term_t arg1)
{
  double sample;
  int ret;
  term_t out;

  sample=uniform_sample();
  out=PL_new_term_ref();
  ret=PL_put_float(out,sample);
  RETURN_IF_FAIL
  return PL_unify(out,arg1);
}
double uniform_sample()
{
  return ((double)rand())/RAND_MAX;
}

static foreign_t gauss_sample_pl(term_t arg1,term_t arg2,term_t arg3)
{
  double mean, var, sample;
  int ret;
  term_t out;

  ret=PL_get_float(arg1,&mean);
  RETURN_IF_FAIL
  ret=PL_get_float(arg2,&var);
  RETURN_IF_FAIL
  sample=gauss_sample(mean,var);
  out=PL_new_term_ref();
  ret=PL_put_float(out,sample);
  RETURN_IF_FAIL
  return PL_unify(out,arg3);
}

double gauss_sample(double mean,double var)
{
  double u1,u2,r,theta,s; 

  u1= uniform_sample();
  u2= uniform_sample();
  r= sqrt(-2*log(u1));
  theta=2*M_PI*u2;
  s=r*cos(theta);
  return sqrt(var)*s+mean;
}
static foreign_t gamma_sample_pl(term_t arg1,term_t arg2,term_t arg3)
{
  double shape, scale, sample;
  int ret;
  term_t out;

  ret=PL_get_float(arg1,&shape);
  RETURN_IF_FAIL
  ret=PL_get_float(arg2,&scale);
  RETURN_IF_FAIL
  sample=gamma_sample(shape,scale);
  out=PL_new_term_ref();
  ret=PL_put_float(out,sample);
  RETURN_IF_FAIL
  return PL_unify(out,arg3);
}
double gamma_sample(double shape, double scale)
{
  double u,s;
  if (shape>=1)
    return gamma_sample_gt1(shape)*scale;
  else
  {
    u=uniform_sample();
    s=gamma_sample_gt1(shape+1);
    return pow(s*u,1/shape)*scale;
  }
}
double gamma_sample_gt1(double shape)
{
  double c,d,x,v,u;

  d=shape-1.0/3.0;
  c =1.0/sqrt(9.0*d);

  do
  {
    do
    {
      x=gauss_sample(0.0,1.0);
      v=pow(1+c*x,3);
    } while (v<=0);
    u=uniform_sample();
  } while (u>=1-0.0331*pow(x,4) && log(u)>=0.5*pow(x,2)+d*(1-v+log(v)));
  return d*v;
}
static foreign_t symmetric_dirichlet_sample_pl(term_t arg1,term_t arg2, term_t arg3)
{
  double  alpha, * sample;

  int ret, i, K;
  term_t out, head;

  head=PL_new_term_ref();
  out=PL_new_term_ref();

  ret=PL_get_integer(arg2,&K);
  RETURN_IF_FAIL
  sample=malloc(sizeof(double)*K);

  ret=PL_get_float(arg1,&alpha);
  RETURN_IF_FAIL
  symmetric_dirichlet_sample(alpha,K,sample);

  ret=PL_put_nil(out);
  RETURN_IF_FAIL
  for (i=0;i<K;i++)
  {
    ret=PL_put_float(head,sample[i]);
    RETURN_IF_FAIL
    ret=PL_cons_list(out,head,out);
    RETURN_IF_FAIL
  }
  return PL_unify(out,arg3);

}
static foreign_t dirichlet_sample_pl(term_t arg1,term_t arg2)
{
  double * alpha, * sample;

  int ret, i;
  size_t K;
  term_t out,alphaterm, head;


  head=PL_new_term_ref();
  out=PL_new_term_ref();
  alphaterm=PL_copy_term_ref(arg1);
  ret=PL_skip_list(alphaterm,0,&K);
  if (ret!=PL_LIST) return FALSE;
  alpha=malloc(sizeof(double)*K);
  sample=malloc(sizeof(double)*K);

  for (i=0;i<K;i++)
  {
    ret=PL_get_list(alphaterm,head,alphaterm);
    RETURN_IF_FAIL
    ret=PL_get_float(head,&alpha[i]);
    RETURN_IF_FAIL
  }
  dirichlet_sample(alpha,K,sample);
  ret=PL_put_nil(out);
  RETURN_IF_FAIL
  for (i=0;i<K;i++)
  {
    ret=PL_put_float(head,sample[i]);
    RETURN_IF_FAIL
    ret=PL_cons_list(out,head,out);
    RETURN_IF_FAIL
  }
  return PL_unify(out,arg2);

}
static foreign_t discrete_sample_pl(term_t arg1,term_t arg2)
{
  double * theta;
  double u, p;

  int ret, i;
  size_t K;
  term_t out,thetaterm, head;


  head=PL_new_term_ref();
  out=PL_new_term_ref();
  thetaterm=PL_copy_term_ref(arg1);
  ret=PL_skip_list(thetaterm,0,&K);
  if (ret!=PL_LIST) return FALSE;
  theta=malloc(sizeof(double)*K);

  for (i=0;i<K;i++)
  {
    ret=PL_get_list(thetaterm,head,thetaterm);
    RETURN_IF_FAIL
    ret=PL_get_float(head,&theta[i]);    
    RETURN_IF_FAIL
  }
  u=uniform_sample();
  i=0;
  p=theta[0];
  while (u>p && i<K)
  {
    i++;
    p=p+theta[i];
  }
  ret=PL_put_integer(out,i);
  RETURN_IF_FAIL
  free(theta);
  return PL_unify(out,arg2);
}
void symmetric_dirichlet_sample(double alpha,int K, double * theta)
{
  int i;
  double * alphas;

  alphas=malloc(sizeof(double)*K);

  for (i=0;i<K;i++)
    alphas[i]=alpha;
  dirichlet_sample(alphas,K,theta);
  free(alphas);
}

void dirichlet_sample(double * alpha,int K, double * theta)
{
  int i;
  double sum;
  double * gamma;

  gamma=malloc(sizeof(double)*K);

  sum=0.0;
  for (i=0;i<K;i++)
  {
    gamma[i]=gamma_sample(alpha[i],1.0);
    sum=sum+gamma[i];
  }
  for (i=0;i<K;i++)
    theta[i]=gamma[i]/sum;
  free(gamma);
}

static foreign_t init_em(term_t arg1)
{
  int ret;
  example_data * ex_d;

  term_t ex_d_t=PL_new_term_ref();
  ex_d=(example_data *)malloc(sizeof(example_data));

  ex_d->ex=0;
  ex_d->nRules=0;
  ex_d->env=NULL;
  ex_d->eta=NULL;
  ex_d->eta_temp=NULL;
  ex_d->rules=NULL;
  ex_d->nodes_probs=NULL;
  ex_d->tunable_rules=NULL;
  ex_d->arrayprob=NULL;
  ex_d->alpha=0.0;

  ret=PL_put_pointer(ex_d_t,(void *)ex_d);
  RETURN_IF_FAIL
  return(PL_unify(ex_d_t,arg1));

}

static foreign_t initial_values_pl(term_t arg1, term_t arg2)
{
  example_data * ex_d;

  int ret;

  ret=PL_get_pointer(arg1,(void **)&ex_d);
  RETURN_IF_FAIL
  ret=PL_get_float(arg2,&(ex_d->alpha)); // <- MOD arg1 -> arg2
  RETURN_IF_FAIL
  PL_succeed;
}
static foreign_t init_ex(term_t arg1, term_t arg2)
{
  example_data * ex_d;
  DdManager * mgr;
  term_t env_t;
  int ex,ret;

  env_t=PL_new_term_ref();
  ret=PL_get_pointer(arg1,(void **)&ex_d);
  RETURN_IF_FAIL
  ex=ex_d->ex;
  ex_d->env=(environment *) realloc(ex_d->env, (ex+1)*sizeof(environment));
  ex_d->env[ex].mgr=Cudd_Init(0,0,UNIQUE_SLOTS,CACHE_SLOTS,5120);
  mgr=ex_d->env[ex].mgr;
  // Cudd_AutodynEnable(mgr, CUDD_REORDER_GROUP_SIFT);
  Cudd_SetMaxCacheHard(mgr, 0);
  Cudd_SetLooseUpTo(mgr, 0);
  Cudd_SetMinHit(mgr, 15);

  ex_d->env[ex].bVar2mVar=NULL;

  ex_d->env[ex].vars=NULL;

  ex_d->env[ex].nVars=0;

  ex_d->env[ex].probs=NULL;

  ex_d->env[ex].boolVars=0;

  ex_d->env[ex].nRules=ex_d->nRules;

  ex_d->env[ex].rules=ex_d->rules;

  ret=PL_put_pointer(env_t,(void *) (ex_d->env+ex));
  RETURN_IF_FAIL
  return(PL_unify(env_t,arg2));

}

static foreign_t end_ex(term_t arg1)
{
  int ret;
  example_data *ex_d;

  ret=PL_get_pointer(arg1,(void **)&ex_d);
  RETURN_IF_FAIL
  ex_d->ex=ex_d->ex+1;
  PL_succeed;
}

static foreign_t init(term_t arg1)
{
  term_t env_t;
  environment * env;
  int ret;

  env_t=PL_new_term_ref();
  env=(environment *)malloc(sizeof(environment));
  // env->mgr=Cudd_Init(0,0,UNIQUE_SLOTS,CACHE_SLOTS,0);
  env->mgr=Cudd_Init(0,0,CUDD_UNIQUE_SLOTS,CACHE_SLOTS,0);

  Cudd_AutodynEnable(env->mgr, CUDD_REORDER_GROUP_SIFT);
  Cudd_SetMaxCacheHard(env->mgr, 0);
  Cudd_SetLooseUpTo(env->mgr, 0);
  Cudd_SetMinHit(env->mgr, 15);
  // Cudd_SetMaxMemory(env->mgr,8000000);

  env->bVar2mVar=NULL;
  env->vars=NULL;
  env->nVars=0;
  env->probs=NULL;
  env->boolVars=0;
  env->nRules=0;
  env->rules= NULL;
  env->n_abd=0;
  env->n_abd_boolVars=0;

  ret=PL_put_pointer(env_t,(void *) env);
  RETURN_IF_FAIL
  return(PL_unify(env_t,arg1));
}

static foreign_t end(term_t arg1)
{
  int ret;
  environment *env;

  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL

  // printf("Pre - Cudd_ReadMemoryInUse(env->mgr) = %ld bytes\n",Cudd_ReadMemoryInUse(env->mgr));
  // printf("Cudd_CheckZeroRef(env->mgr) = %d\n",Cudd_CheckZeroRef(env->mgr));
  
  // debug_cudd_env(env,0);
  
  Cudd_Quit(env->mgr);

  free(env->bVar2mVar);
  free(env->vars);
  free(env->probs);
  free(env->rules);
  free(env);

  PL_succeed;
}



static double Expectation(example_data * ex_d,DdNode **nodes_ex,int lenNodes)
{
  int i;
  double rootProb,CLL=0;

  for(i=0;i<lenNodes;i++)
  {
    if (!Cudd_IsConstant(nodes_ex[i]))
    {
      ex_d->nodesB=init_table(ex_d->env[i].boolVars);
      ex_d->nodesFE=init_table(ex_d->env[i].boolVars);
      ex_d->nodesFO=init_table(ex_d->env[i].boolVars);

      Forward(ex_d,nodes_ex[i],i);
      rootProb=GetOutsideExpe(ex_d,nodes_ex[i],ex_d->example_prob[i],i);

      if (rootProb<=0.0)
        CLL = CLL + LOGZERO*ex_d->example_prob[i];
      else
        CLL = CLL + log(rootProb)*ex_d->example_prob[i];

      ex_d->nodes_probs[i]=rootProb;
      destroy_table(ex_d->nodesB,ex_d->env[i].boolVars);
      destroy_table(ex_d->nodesFE,ex_d->env[i].boolVars);
      destroy_table(ex_d->nodesFO,ex_d->env[i].boolVars);
    }
    else
      if (nodes_ex[i]==Cudd_ReadLogicZero(ex_d->env[i].mgr))
      {
        CLL=CLL+LOGZERO*ex_d->example_prob[i];
	ex_d->nodes_probs[i]=0.0;
      }
      else
        ex_d->nodes_probs[i]=1.0;
  }
  return CLL;
}

static foreign_t end_em(term_t arg1)
{
  int r,i,ret;
  example_data * ex_d;
  ret=PL_get_pointer(arg1,(void **)&ex_d);
  RETURN_IF_FAIL

  for (i=0;i<ex_d->ex;i++)
  {
    Cudd_Quit(ex_d->env[i].mgr);
    free(ex_d->env[i].bVar2mVar);
    free(ex_d->env[i].vars);
    free(ex_d->env[i].probs);
  }

  free(ex_d->env);
  for (r=0;r<ex_d->nRules;r++)
  {
    if (ex_d->tunable_rules[r])
    {
      for (i=0;i<ex_d->rules[r]-1;i++)
      {
        free(ex_d->eta[r][i]);
        free(ex_d->eta_temp[r][i]);
      }
      free(ex_d->eta[r]);
      free(ex_d->eta_temp[r]);
    }
  }
  free(ex_d->eta);
  free(ex_d->eta_temp);
  free(ex_d->rules);
  free(ex_d);
  PL_succeed;
}


static foreign_t ret_prob(term_t arg1, term_t arg2, term_t arg3)
{
  term_t out;
  environment * env;
  DdNode * node;
  tablerow * table;
  double prob;
  int ret;

  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL
  ret=PL_get_pointer(arg2,(void **)&node);
  RETURN_IF_FAIL
  out=PL_new_term_ref();

  if (!Cudd_IsConstant(node))
  {
    table=init_table(env->boolVars);
    prob=Prob(node,env,table);
    if (Cudd_IsComplement(node))
      prob=1.0-prob;
    ret=PL_put_float(out,prob);
    RETURN_IF_FAIL
    destroy_table(table,env->boolVars);
  }
  else
  {
    if (node==Cudd_ReadOne(env->mgr))
    {
      ret=PL_put_float(out,1.0);
      RETURN_IF_FAIL
    }
    else
    {
      ret=PL_put_float(out,0.0);
      RETURN_IF_FAIL
    }
  }

  return(PL_unify(out,arg3));
}

int reorder_int(environment *env)
{
  int i,j,var_ind,abd_ind=0,ind=env->n_abd_boolVars;
  variable var,* vars=env->vars;
  DdManager *mgr=env->mgr;
  // int boolVars=env->boolVars;
  int * permutation;
  int * bVar2mVar=env->bVar2mVar;

  // permutation=malloc((boolVars)*sizeof(int));
  permutation = malloc((Cudd_ReadSize(env->mgr))*sizeof(int));

  for (i=0;i<Cudd_ReadSize(env->mgr);i++)
  {
    j=Cudd_ReadInvPerm(mgr,i);
    var_ind=bVar2mVar[j];
    var=vars[var_ind];
    if (var.abducible || var.query)
    {
      permutation[abd_ind]=j;
      abd_ind++;
    }
    else
    {
      permutation[ind]=j;
      ind++;
    }

  }

  j = Cudd_ShuffleHeap(mgr,permutation);

  free(permutation);

  return j;
}

static foreign_t reorder(term_t arg1)
{
  environment * env;
  int ret;

  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL
  ret=reorder_int(env);
  RETURN_IF_FAIL
  return 1;
}

double Prob_given_expl(DdNode *node, environment *env, tablerow *table, explan_t *expl, int comp_par)
{
  int index, comp, pos;
  double p, pf, pt, res;
  double * value_p;
  explan_t *exp = expl; // copia perché scorro 
  comp = Cudd_IsComplement(node);
  comp = (comp && !comp_par) ||(!comp && comp_par);
  // printf("Comp: %d\n",comp);
  if(Cudd_IsConstant(node)) {
    p = 1;
    // if (comp)
      // p = 1.0 - 1;
    // printf("Terminal: %lf\n",p);

    return p;
  }
  index = Cudd_NodeReadIndex(node);
  pos = Cudd_ReadPerm(env->mgr,index);
  if (pos >= env->n_abd_boolVars)
  {
    p = Prob(node,env,table);
    // if (comp)
    //   p = 1.0 - p;

    // printf("Computed prob: %lf\n",p);

    return p;
  }
  
  value_p = get_value(table,Cudd_Regular(node));
  if (value_p != NULL) {
    return *value_p;
  }

  p = env->probs[index];
  // printf("cerco indice: %d che ha prob: %lf\n",index,p);
  while(exp != NULL && exp->a.var != index) {
    exp = exp->next;
  }

  pt = 1;
  pf = 1;
  if(exp != NULL) {
    if(exp->a.val == 1) {
      // printf("selected true: pt %lf\n",pt);
      pt = Prob_given_expl(Cudd_T(node),env,table,expl,comp);
      pf = 0;
    }
    else {
      // printf("selected false pf: %lf\n",pf);
      pf = Prob_given_expl(Cudd_E(node),env,table,expl,comp);
      if(Cudd_IsComplement(Cudd_E(node))) {
        pf = 1 - pf;
      }
      pt = 0;
    }
  }
  else {
    // printf("Ass not found - ERRORE\n");
  }
  
  // printf("Indice: %d, p: %lf, pf: %lf, pt: %lf\n",index,p,pf,pt);
  res = pf*(1-p) + pt*p;
  // printf("res: %lf\n",res);
  add_node(table,Cudd_Regular(node),res);

  return res;
}

static foreign_t ret_abd_prob(term_t arg_env, term_t arg_bdd, term_t arg_prob, term_t arg_expl)
{
  term_t out,outass,out_assign;
  environment * env;
  DdNode * node;
  expltablerow_abd * expltable;
  tablerow * table;
  //abdtablerow * abdtable;
  prob_abd_expl_list delta;
  // double *prob_ics = NULL;
  int ret;
  // double p;
  // explan_t ** mpa;
  int i;
  // assign *array_mpa;

  ret=PL_get_pointer(arg_env,(void **)&env);
  RETURN_IF_FAIL
  ret=PL_get_pointer(arg_bdd,(void **)&node);
  RETURN_IF_FAIL
  // ret=PL_get_pointer(arg_bdd_ic,(void **)&node_ic);
  // RETURN_IF_FAIL

  out=PL_new_term_ref();

  ret=reorder_int(env);
  RETURN_IF_FAIL

  // array_mpa = malloc(sizeof(assign) * env->boolVars);

  if (!Cudd_IsConstant(node))
  {
    expltable = expl_init_table_abd(env->boolVars);
    table=init_table(env->boolVars);

    delta.n_elements = 0;
    delta.mpa = NULL;
    delta.prob = -1;

    delta = abd_Prob(node,env,expltable,table,0);

    ret = PL_put_float(out,delta.prob);
    RETURN_IF_FAIL

    out_assign = PL_new_term_ref();
    PL_put_nil(out_assign);
    RETURN_IF_FAIL
    for(i = 0; i < delta.n_elements; i++) {
      outass=abd_clist_to_pllist(delta.mpa[i]);
      
      ret = PL_cons_list(out_assign,outass,out_assign);
      RETURN_IF_FAIL
    }

    // expl_destroy_table_abd(expltable,env->boolVars);
    destroy_table(table,env->boolVars);
  }
  else
  {
    if (node==Cudd_ReadOne(env->mgr))
    {
      ret=PL_put_float(out,1.0);
      RETURN_IF_FAIL
    }
    else
    {
      ret=PL_put_float(out,0.0);
      RETURN_IF_FAIL
    }
      out_assign=PL_new_term_ref();
      ret=PL_put_nil(out_assign);
      RETURN_IF_FAIL
  }

  return(PL_unify(out,arg_prob)&&PL_unify(out_assign,arg_expl));
}

static foreign_t ret_map_prob(term_t arg1, term_t arg2,
  term_t arg3, term_t arg4)
{
  term_t out,outass;
  environment * env;
  DdNode * node;
  expltablerow * maptable;
  tablerow * table;
  //abdtablerow * abdtable;
  prob_abd_expl delta;
  int ret;
  double p;
  explan_t * mpa;

  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL
  ret=PL_get_pointer(arg2,(void **)&node);
  RETURN_IF_FAIL
  out=PL_new_term_ref();

  ret=reorder_int(env);

  RETURN_IF_FAIL

  if (!Cudd_IsConstant(node))
  {
    maptable=expl_init_table(env->boolVars);
    table=init_table(env->boolVars);
    //abdtable=init_abd_table(env->n_abd);
    delta=map_Prob(node,env,maptable,table,0);
    p=delta.prob;
    mpa=delta.mpa;
    ret=PL_put_float(out,p);
    RETURN_IF_FAIL
    //destroy_table(abdtable,env->n_abd);
    outass=clist_to_pllist(mpa,env);
    RETURN_IF_FAIL
    expl_destroy_table(maptable,env->boolVars);
    destroy_table(table,env->boolVars);
  }
  else
  {
    if (node==Cudd_ReadOne(env->mgr))
    {
      ret=PL_put_float(out,1.0);
      RETURN_IF_FAIL
    }
    else
    {
      ret=PL_put_float(out,0.0);
      RETURN_IF_FAIL
    }
    outass=PL_new_term_ref();
    ret=PL_put_nil(outass);
    RETURN_IF_FAIL
  }

  return(PL_unify(out,arg3)&&PL_unify(outass,arg4));
}

static foreign_t make_query_var(term_t arg1, term_t arg2, term_t arg3)
{
  environment * env;
  int ret,varIndex,i,j;
  DdNode * cons, * tmp0,* tmp1, * tmp2, * vari, * varlast, * varj, * or, * tmpor;
  term_t out;
  variable var;

  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL
  ret=PL_get_integer(arg2,&varIndex);
  /*  printf(" V = varIndex = %d\n",varIndex); */
  RETURN_IF_FAIL

  var=env->vars[varIndex];

  cons=Cudd_ReadOne(env->mgr);
  or=Cudd_ReadLogicZero(env->mgr);

  for (i=var.firstBoolVar; i<var.firstBoolVar+var.nVal-1; i++)
  {    
    vari=Cudd_bddIthVar(env->mgr,i);
    tmpor=Cudd_bddOr(env->mgr,or,vari);
    Cudd_Ref(tmpor);
    Cudd_RecursiveDeref(env->mgr,or);
    or=tmpor;
    for(j=i+1; j<var.firstBoolVar+var.nVal; j++)
    {      
      varj=Cudd_bddIthVar(env->mgr,j);
      tmp0=Cudd_bddAnd(env->mgr,vari,varj);
      Cudd_Ref(tmp0);//added
      tmp1=Cudd_Not(tmp0);
      Cudd_Ref(tmp1);//added
      tmp2=Cudd_bddAnd(env->mgr,cons,tmp1);
      Cudd_Ref(tmp2);
      cons=tmp2;
      Cudd_Ref(cons);//added
    }
  }
  varlast=Cudd_bddIthVar(env->mgr,var.firstBoolVar+var.nVal-1);
  tmpor=Cudd_bddOr(env->mgr,or,varlast);
  Cudd_Ref(tmpor);
  Cudd_RecursiveDeref(env->mgr,or);
  tmp1=Cudd_bddAnd(env->mgr,cons,tmpor);
  Cudd_Ref(tmp1);
  Cudd_RecursiveDeref(env->mgr,cons);
  cons=tmp1;

  out=PL_new_term_ref();
  ret=PL_put_pointer(out,(void *)cons);
  RETURN_IF_FAIL
  return(PL_unify(out,arg3));
}

term_t abd_clist_to_pllist(explan_t *mpa)
{
  term_t out,tail,head,var,val;
  functor_t minus;
  assign a;
  minus=PL_new_functor(PL_new_atom("-"), 2);
  out=PL_new_term_ref();
  head=PL_new_term_ref();
  int ret;
  if (mpa==NULL)
  {
    ret=PL_put_nil(out);
    RETURN_IF_FAIL
  }
  else
  {
    tail=abd_clist_to_pllist(mpa->next);
    a=mpa->a;
    var=PL_new_term_ref();
    val=PL_new_term_ref();
    ret=PL_put_integer(var,a.var);
    RETURN_IF_FAIL
    ret=PL_put_integer(val,a.val);
    RETURN_IF_FAIL
    ret=PL_cons_functor(head, minus,var,val);
    RETURN_IF_FAIL
    ret=PL_cons_list(out,head,tail);
    RETURN_IF_FAIL
  }
  return out;
}
term_t clist_to_pllist(explan_t *mpa, environment * env)
{
  term_t out,tail,head,var,val;
  functor_t minus;
  assign a;
  int value,bvar, ret, mvari, mval;
  variable mvar;

  minus=PL_new_functor(PL_new_atom("-"), 2);
  tail=PL_new_term_ref();
  ret=PL_put_nil(tail);
  RETURN_IF_FAIL
  if (mpa==NULL)
  {
    out=PL_new_term_ref();
    ret=PL_put_nil(out);
    RETURN_IF_FAIL
  }
  else
  {

    for (; mpa; mpa=mpa->next)
    {
      a=mpa->a;
      bvar=a.var;
      value=a.val;
      if (value)
      {
        mvari=env->bVar2mVar[bvar];
        mvar=env->vars[mvari];
        mval=a.var-mvar.firstBoolVar+1;
        var=PL_new_term_ref();
        val=PL_new_term_ref();
        ret=PL_put_integer(var,mvari);
        RETURN_IF_FAIL
        ret=PL_put_integer(val,mval);
        RETURN_IF_FAIL
        head=PL_new_term_ref();
        ret=PL_cons_functor(head, minus,var,val);
        RETURN_IF_FAIL
        out=PL_new_term_ref();
        ret=PL_cons_list(out,head,tail);
        RETURN_IF_FAIL
        tail=out;
      }
    }
    out=tail;
  }
  return out;
}



double Prob(DdNode *node, environment * env, tablerow * table)
/* compute the probability of the expression rooted at node.
table is used to store nodeB for which the probability has already been computed
so that it is not recomputed
 */
{
  int index;
  double res;
  double p,pt,pf,BChild0,BChild1;
  double * value_p;
  DdNode *nodekey,*T,*F;
  //comp=(comp && !comp_par) ||(!comp && comp_par);
  if (Cudd_IsConstant(node))
  {
      return 1.0;
  }
  else
  {
    nodekey=Cudd_Regular(node);
    value_p=get_value(table,nodekey);
    if (value_p!=NULL)
        return *value_p;
    else
    {
      index=Cudd_NodeReadIndex(node);  //Returns the index of the node. The node pointer can be either regular or complemented.
      //The index field holds the name of the variable that labels the node. The index of a variable is a permanent attribute that reflects the order of creation.
      p=env->probs[index];
      T = Cudd_T(node);
      F = Cudd_E(node);
      pf=Prob(F,env,table);
      pt=Prob(T,env,table);
      if (Cudd_IsComplement(F))
        pf=1.0-pf;

      BChild0=pf*(1-p);
      BChild1=pt*p;
      res=BChild0+BChild1;
      add_node(table,nodekey,res);
      return res;
    }
  }
}

prob_abd_expl_list abd_Prob(DdNode *node, environment * env, expltablerow_abd * expltable, tablerow * table, int comp_par)
/* compute the probability of the expression rooted at node.
table is used to store nodeB for which the probability has already been computed
so that it is not recomputed
 */
{
  int index,comp,pos,n_new_elements = 0;
  double p1;
  DdNode *nodekey,*T,*F;
  prob_abd_expl_list deltat,deltaf,delta,*deltaptr;
  assign assignment;
  explan_t ** mpa; //** mpa0,** mpa1;
  // explan_t *mptemp;

  comp = Cudd_IsComplement(node);
  comp=(comp && !comp_par) ||(!comp && comp_par);
  if(Cudd_IsConstant(node)) {
    p1 = 1;
    if (comp)
      p1= 1.0-p1;

    // aggiungo una nuova spiegazione
    delta.prob=p1;
    delta.n_elements = 0;
    delta.mpa=NULL;

    return delta;
  }
  index=Cudd_NodeReadIndex(node);
  pos=Cudd_ReadPerm(env->mgr,index);
  if (pos>=env->n_abd_boolVars)
  {
    p1=Prob(node,env,table);
    if (comp)
      p1= 1.0-p1;

    delta.prob=p1;
    delta.n_elements = 0;
    delta.mpa=NULL;

    return delta;
  }
  else
  {
    nodekey=Cudd_Regular(node);
    deltaptr=expl_get_value_abd(expltable,nodekey,comp);
    if (deltaptr!=NULL)
    {
      // printf("found in table\n");
      return *deltaptr;
    }
    else {

      T = Cudd_T(node);
      F = Cudd_E(node);
      deltaf=abd_Prob(F,env,expltable,table,comp);
      deltat=abd_Prob(T,env,expltable,table,comp);

      assignment.var=env->bVar2mVar[index];
      // printf("variabile: %d\n",assignment.var);
      if (deltat.prob > deltaf.prob) {
        assignment.val=1;
        // printf("T (%lf) > F (%lf) \n",deltat.prob, deltaf.prob);
        // printf("inserisco: %d\n",assignment.var);
        // printf("Prima di insert_abd true\n");
        // print_prob_abd_expl_list(&deltat);
        mpa = insert_abd_ordered(assignment, deltat.mpa, deltat.n_elements);
        delta.prob = deltat.prob;

        // printf("deltat.n_elements: %d\n",deltat.n_elements);

        if(deltat.n_elements == 0) {
          delta.n_elements = 1;
        }
        else {
          delta.n_elements = deltat.n_elements;
        }
      }
      else
      {
        if(deltat.prob == deltaf.prob) {
          // printf("p0 == p1 == %lf\n", deltaf.prob);
          // printf("Prima di insert_abd true and false\ndeltat\n");
          // print_prob_abd_expl_list(&deltat);
          // printf("deltaf\n");
          // print_prob_abd_expl_list(&deltaf);
          
          
          mpa = split_explan_v2(assignment,deltat.mpa,deltaf.mpa,&deltat.n_elements,&deltaf.n_elements, &n_new_elements);
          delta.n_elements = n_new_elements;
        }
        else {
          // printf("Prima di insert_abd false\n");
          // printf("ora non la inserisco\n");
          mpa = deltaf.mpa;
          delta.n_elements = deltaf.n_elements;
        }
        // in entrambi i casi metto p0, visto che p0 e p1 sono uguali
        delta.prob = deltaf.prob;
      }
      
      // printf("delta.n_elements: %d\n",delta.n_elements);
      
      delta.mpa=mpa;

      
      // printf("delta dopo inserimento: \n");
      // print_prob_abd_expl_list(&delta);
      
      expl_add_node_abd(expltable,nodekey,comp,delta);
      return delta;
    }
  }
}

prob_abd_expl map_Prob(DdNode *node, environment * env,
  expltablerow * maptable, tablerow * table,
  int comp_par)
/* compute the probability of the expression rooted at node.
table is used to store nodeB for which the probability has already been computed
so that it is not recomputed
 */
{
  int index,comp,pos;
  double p,p0,p1;
  DdNode *nodekey,*T,*F;
  prob_abd_expl deltat,deltaf,delta,*deltaptr;
  assign assignment;
  explan_t * mpa0,* mpa1,* mpa;


  index=Cudd_NodeReadIndex(node);
  pos=Cudd_ReadPerm(env->mgr,index);
  comp=Cudd_IsComplement(node);
  comp=(comp && !comp_par) ||(!comp && comp_par);

  if (pos>=env->n_abd_boolVars)
  {
    p1=Prob(node,env,table);
    if (comp)
      p1= 1.0-p1;
    delta.prob=p1;
    delta.mpa=NULL;


    return delta;
  }
  else
  {
    nodekey=Cudd_Regular(node);
    deltaptr=expl_get_value(maptable,nodekey,comp);
    if (deltaptr!=NULL)
    {
      return *deltaptr;
    }
    p=env->probs[index];
    T = Cudd_T(node);
    F = Cudd_E(node);
    deltaf=map_Prob(F,env,maptable,table,comp);
    deltat=map_Prob(T,env,maptable,table,comp);

    p0=deltaf.prob;
    mpa0=deltaf.mpa;

    p1=deltat.prob*p;
    mpa1=deltat.mpa;

    if (p1>p0)
    {
      assignment.var=index;
      assignment.val=1;
      mpa=insert(assignment,mpa1);
      delta.prob=p1;
      delta.mpa=mpa;
    }
    else
    {
      assignment.var=index;
      assignment.val=0;
      mpa=insert(assignment,mpa0);
      delta.prob=p0;
      delta.mpa=mpa;
    }
    expl_add_node(maptable,nodekey,comp,delta);
    return delta;
  }

}

int in_list(explan_t *list, assign element) {
  while (list != NULL) {
    if(list->a.var == element.var) {
      return 1;
    }
    list = list->next;
  }
  return 0;
}

explan_t *merge_explain(explan_t *a, explan_t *b) {
  explan_t *init = NULL;
  explan_t *root = NULL;
  
  if(a == NULL && b == NULL) {
    return NULL;
  }
  if(a == NULL) {
    return duplicate(b);
  }
  if(b == NULL) {
    return duplicate(a);
  }

  // printf("merging a\n");
  while(a != NULL) {
    init = malloc(sizeof(explan_t));
    init->a = a->a;
    init->next = root;
    root = init;
    a = a->next;
  }

  // printf("merging b\n");
  // not super fast, but it is simple to implement
  while(b != NULL) {
    if(in_list(root,b->a) <= 0) {
      init = malloc(sizeof(explan_t));
      init->a = b->a;
      init->next = root;
      root = init;
    } 
    b = b->next;
  }

  return root;
}

explan_t **split_explan(assign assignment, explan_t **head_true, explan_t **head_false, int n_true, int n_false, int *n_new_elements, int no_false) {

  explan_t **newhead;

  int i;
  int dim = 0;

  if(n_true == 0) {
    n_true++;
  }
  if(n_false == 0) {
    n_false++;
  }
  dim = n_true + n_false;

  newhead = malloc(sizeof(explan_t*) * (dim + 2));

  // inserisco lato true
  assignment.val = 1 - no_false;

  for(i = 0; i < n_true; i++) {
    newhead[i] = malloc(sizeof(explan_t));
    newhead[i]->a = assignment;
    if(head_true != NULL) {
      newhead[i]->next = duplicate(head_true[i]);
    }
    else {
      newhead[i]->next = NULL;
    }
  }

  assignment.val = 0;
  for(i = 0; i < n_false; i++) {
    newhead[n_true + i] = malloc(sizeof(explan_t));
    newhead[n_true + i]->a = assignment;
    if(head_false != NULL) {
      newhead[n_true + i]->next = duplicate(head_false[i]);
    }
    else {
      newhead[n_true + i]->next = NULL;
    }
  }

  // printf("dim: %d\n",dim);

  *n_new_elements = dim;

  return newhead;

}

// returns 1 if e01 is subset of e02
// e02 ha più elementi
int is_subset(explan_t *e01, explan_t *e02) {
  explan_t *e1, *e2;
  int moved = 0;
  e1 = e01;
  e2 = e02;

  if(e1 == NULL) {
    return 1;
  }
  else if(e2 == NULL) {
    return 0;
  }
  else {
    while(e1 != NULL && e2 != NULL && e2->a.var < e1->a.var) {
      // allineo gli insiemi
      // printf("scorro\n");
      e2 = e2->next;
      moved++;
    }
    while(e1 != NULL && e2 != NULL && (e1->a.var == e2->a.var)) {
      e1 = e1->next;
      e2 = e2->next;
    }
    // if(e1 == NULL) {
    if((e1 == NULL && e2 != NULL) || (e1 == NULL && e2 == NULL && moved != 0)) {
      return 1;
    }
    else {
      return 0;
    }
  }
}

explan_t **split_explan_v2(assign assignment, explan_t **head_true, explan_t **head_false, int *n_true, int *n_false, int *n_new_elements) {
  explan_t **newhead;

  int i,ii,iii;
  int dim = 0;
  int pos = 0;
  int removed = 0;

  // check whether one of the two is empty: if so, empty explanation
  // this does not apply to n_true since i need to insert the element 1
  // if(*n_true == 0) {
    // *n_new_elements = 0;
    // return head_true;
  // }
  if(*n_false == 0) {
    *n_new_elements = 0;
    return head_false;
  }
  else {
    // checks dominated sets in true
    for(i = 0; i < *n_true; i++) {
      for(ii = 0; ii < *n_false; ii++) {
        if(is_subset(head_false[ii],head_true[i])) {
          // remove head_true[i]
          // printf("rimuovo spiegazione %d da true\n",i);
          removed++;
          for(iii = i; iii < *n_true - 1; iii++) {
            head_true[iii] = head_true[iii + 1];
          }
          // free(head_true[iii]);
          *n_true = *n_true - 1;
        }
      }
    }

    // add a node in true, as before, and merge the sets true and false
    assignment.val = 1;
    dim = *n_true + *n_false;
    if(*n_true == 0 && removed > 0) {
      *n_new_elements = *n_false;
      return head_false;
    }
    else if(*n_true == 0) {
      dim++;
    }

    newhead = malloc(sizeof(explan_t*) * (dim + 2)); // perché + 2?

    if(*n_true == 0) {
      // insert empty assignment
      newhead[0] = insord_abd(assignment,NULL);
      pos++;
    }
    else {
      for (i = 0; i < *n_true; i++) {
        // insert empty assignment
        newhead[i] = insord_abd(assignment,head_true[i]);
        pos++;
      }
    }

    // add false
    for(i = 0; i < *n_false; i++) {
      newhead[pos + i] = head_false[i];
    }

    *n_new_elements = dim;

    return newhead;
  }
}

explan_t **insert_abd(assign assignment, explan_t **head, int n_elements) {
  int i;
  explan_t **newhead;

  newhead = malloc(sizeof(explan_t*) * (n_elements == 0 ? 1 : n_elements + 1));
  
  if(n_elements == 0) {
    newhead[0] = malloc(sizeof(explan_t));
    newhead[0]->a = assignment;
    newhead[0]->next = NULL;
  }
  else {
    for(i = 0; i < n_elements; i++) {
      newhead[i] = malloc(sizeof(explan_t));
      newhead[i]->a = assignment;
      if(head != NULL && head[i] != NULL) {
        newhead[i]->next = duplicate(head[i]);
      }
      else {
        newhead[i]->next = NULL;
      }
      // printf("newhead[%d]->a.var = %d\n",i,newhead[i]->a.var);
      // printf("newhead[%d]->a.val = %d\n",i,newhead[i]->a.val);
    }
  }

  return newhead;
}

explan_t **insert_abd_ordered(assign assignment, explan_t **head, int n_elements) {
  int i;
  explan_t **newhead;
  
  // questa si può risparmiare alcune volte
  newhead = malloc(sizeof(explan_t*) * (n_elements == 0 ? 1 : n_elements + 1));
  
  if(n_elements == 0) {
    newhead[0] = malloc(sizeof(explan_t));
    newhead[0]->a = assignment;
    newhead[0]->next = NULL;
  }
  else {
    for(i = 0; i < n_elements; i++) {
      newhead[i] = insord_abd(assignment,head[i]);
    }
    // newhead = head;
  }

  return newhead;
}

explan_t *insord_abd(assign assignment, explan_t *head) {
  if(head == NULL) {
    return insert(assignment, head);
  }
  else if(assignment.var < head->a.var) {
    return insert(assignment,head);
  }
  else{
    head->next = insert(assignment,head->next);
    return head;
  }
} 

explan_t *insord_abd_iterative(assign assignment, explan_t *head) {
  explan_t* prec = NULL, *patt = head, *aux;
  int found = 0;
  while(patt != NULL && (found == 0)) {
    if(assignment.var < patt->a.var) {
      found = 1;
    }
    else {
      prec = patt;
      patt = patt->next;
    }
  }
  
  aux = malloc(sizeof(explan_t));
  aux->a = assignment;
  aux->next = patt;

  if(patt == head) {
    return aux;
  }
  else {
    prec->next = aux;
    return head;
  }
}

explan_t * insert(assign assignment,explan_t * head)
{
  explan_t * newhead;

  newhead=malloc(sizeof(explan_t));
  newhead->a=assignment;
  newhead->next=duplicate(head);
  return newhead;
}

explan_t * duplicate(explan_t * head)
{
  explan_t * newhead;
  if (head)
  {
    newhead=malloc(sizeof(explan_t));
    newhead->a=head->a;
    newhead->next=duplicate(head->next);
  }
  else
    newhead=NULL;
  return newhead;
}

void free_list(explan_t * head)
{
  if (head)
  {
    free_list(head->next);
    free(head);
  }
}
static foreign_t add_var(term_t arg1,term_t arg2,term_t arg3,term_t arg4)
{
  term_t out,head,probTerm;
  variable * v;
  int i,ret,nRules;
  size_t lenProbs;
  double p,p0;
  environment * env;

  head=PL_new_term_ref();
  out=PL_new_term_ref();
  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL
  env->nVars=env->nVars+1;
  env->vars=(variable *) realloc(env->vars,env->nVars * sizeof(variable));

  v=&env->vars[env->nVars-1];
  v->abducible=0;
  v->query=0;
  v->decision=0;

  ret=PL_get_integer(arg3,&v->nRule);
  RETURN_IF_FAIL
  probTerm=PL_copy_term_ref(arg2);
  ret=PL_skip_list(probTerm,0,&lenProbs);
  if (ret!=PL_LIST) return FALSE;

  v->nVal=lenProbs;

  if(v->nVal == 1) {
    v->nVal++;
  }

  nRules=env->nRules;
  if (v->nRule>=nRules)
  {
    env->rules=(int *)  realloc(env->rules,((v->nRule+1)* sizeof(int)));
    for (i=nRules;i<v->nRule;i++)
      env->rules[i]=0;
    env->rules[v->nRule]=lenProbs;
    env->nRules=v->nRule+1;
  }
  v->firstBoolVar=env->boolVars;
  env->probs=(double *) realloc(env->probs,(((env->boolVars+v->nVal-1)* sizeof(double))));
  env->bVar2mVar=(int *) realloc(env->bVar2mVar,((env->boolVars+v->nVal-1)* sizeof(int)));
  
  p0=1;
  for (i=0;i<v->nVal-1;i++)
  {
    ret=PL_get_list(probTerm,head,probTerm);
    RETURN_IF_FAIL
    ret=PL_get_float(head,&p);
    RETURN_IF_FAIL
    env->bVar2mVar[env->boolVars+i]=env->nVars-1;
    env->probs[env->boolVars+i]=p/p0;
    p0=p0*(1-p/p0);
  }
  env->boolVars=env->boolVars+v->nVal-1;
  env->rules[v->nRule]= v->nVal;

  ret=PL_put_integer(out,env->nVars-1);
  RETURN_IF_FAIL

  return(PL_unify(out,arg4));
}

static foreign_t add_query_var(term_t arg1,term_t arg2,term_t arg3,term_t arg4)
{
  term_t out,head,probTerm;
  variable * v;
  int i,ret,nRules;
  size_t lenProbs;
  double p;
  environment * env;

  head=PL_new_term_ref();
  out=PL_new_term_ref();
  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL
  env->nVars=env->nVars+1;
  env->n_abd++;
  env->vars=(variable *) realloc(env->vars,env->nVars * sizeof(variable));

  v=&env->vars[env->nVars-1];
  v->query=1;
  v->abducible=0;
  v->decision=0;

  probTerm=PL_copy_term_ref(arg2);

  ret=PL_skip_list(probTerm,0,&lenProbs);
  if (ret!=PL_LIST) return FALSE;
  v->nVal=lenProbs;

  ret=PL_get_integer(arg3,&v->nRule);
  RETURN_IF_FAIL
  nRules=env->nRules;
  if (v->nRule>=nRules)
  {
    env->rules=(int *)  realloc(env->rules,((v->nRule+1)* sizeof(int)));
    for (i=nRules;i<v->nRule;i++)
      env->rules[i]=0;
    env->rules[v->nRule]=lenProbs;
    env->nRules=v->nRule+1;
  }

  env->n_abd_boolVars=env->n_abd_boolVars+v->nVal;
  v->firstBoolVar=env->boolVars;
  env->probs=(double *) realloc(env->probs,(((env->boolVars+v->nVal)* sizeof(double))));
  env->bVar2mVar=(int *) realloc(env->bVar2mVar,((env->boolVars+v->nVal)* sizeof(int)));

  for (i=0;i<v->nVal;i++)
  {
    ret=PL_get_list(probTerm,head,probTerm);
    RETURN_IF_FAIL
    ret=PL_get_float(head,&p);
    RETURN_IF_FAIL
    env->bVar2mVar[env->boolVars+i]=env->nVars-1;
    env->probs[env->boolVars+i]=p;
  }
  env->boolVars=env->boolVars+v->nVal;
  env->rules[v->nRule]= v->nVal;

  ret=PL_put_integer(out,env->nVars-1);
  RETURN_IF_FAIL

  return(PL_unify(out,arg4));
}

static int add_decision_var_impl(environment *env, int rule, int nvals, term_t vout) {
  term_t out;
  variable *v;
  int i, ret, nRules;

  out = PL_new_term_ref();

  env->nVars = env->nVars+1; // increase the counter of multival var
  env->vars = (variable *) realloc(env->vars,env->nVars * sizeof(variable));
  if (!env->vars) return FALSE;

  v = &env->vars[env->nVars-1];
  v->query = 0;
  v->abducible = 0;
  v->decision = 1; // set decision var

  v->nVal = nvals;
  v->nRule = rule;
  
  nRules=env->nRules;
  if (v->nRule>=nRules) {
    env->rules = (int *)  realloc(env->rules,((v->nRule+1)* sizeof(int)));
    if (!env->rules) return FALSE;
    for (i=nRules;i<v->nRule;i++)
      env->rules[i] = 0;
    env->rules[v->nRule] = v->nVal;
    env->nRules = v->nRule + 1;
  }

  v->firstBoolVar = env->boolVars;
  env->probs = (double *) realloc(env->probs,(((env->boolVars+v->nVal-1)* sizeof(double))));
  env->bVar2mVar = (int *) realloc(env->bVar2mVar,((env->boolVars+v->nVal-1)* sizeof(int)));
  if (!env->probs || !env->bVar2mVar) return FALSE;

  for (i = 0; i < v->nVal - 1; i++) {
    env->bVar2mVar[env->boolVars + i] = env->nVars - 1;
    env->probs[env->boolVars + i] = -1; // set the prob of decision var to -1
  }
  
  env->boolVars = env->boolVars + (v->nVal - 1);
  env->rules[v->nRule] = v->nVal;

  ret=PL_put_integer(out,env->nVars-1);
  RETURN_IF_FAIL

  return(PL_unify(out,vout));
}

static foreign_t add_decision_var(term_t env_ref, term_t current_n_rule, term_t nvals_term, term_t vout) {
  int ret, rule, nvals;
  environment *env;

  ret = PL_get_pointer(env_ref,(void **)&env);
  RETURN_IF_FAIL
  ret = PL_get_integer(current_n_rule, &rule);
  RETURN_IF_FAIL
  ret = PL_get_integer(nvals_term, &nvals);
  RETURN_IF_FAIL

  if (nvals < 2) return FALSE;
  return add_decision_var_impl(env, rule, nvals, vout);
}

static int compare_term_entry(const void *a, const void *b) {
  const term_entry *ta = (const term_entry *)a;
  const term_entry *tb = (const term_entry *)b;
  if (ta->impact > tb->impact) return -1;
  if (ta->impact < tb->impact) return 1;
  return 0;
}

DdNode* setLowerBound(DdManager *dd, DdNode *current_node, double lowerBound, double floor_value) {
  autodyn_guard_t guard;
  autodyn_guard_begin(dd, &guard);
  DdNode *res = setLowerBound_impl(dd, current_node, lowerBound, floor_value);
  autodyn_guard_end(dd, &guard);
  return res;
}

static DdNode* setLowerBound_impl(DdManager *dd, DdNode *current_node, double lowerBound, double floor_value) {
  DdNode *res, *then_branch, *else_branch, *T, *E;
  DdNode *reg = Cudd_Regular(current_node);
  if (Cudd_IsConstant(reg)) {
    double v = Cudd_V(reg);
    if (v < lowerBound) {
      res = Cudd_addConst(dd, floor_value);
      if (!res) return NULL;
      Cudd_Ref(res);
      return res;
    }
    Cudd_Ref(current_node);
    return current_node;
  }

  then_branch = Cudd_T(reg);
  else_branch = Cudd_E(reg);

  T = setLowerBound_impl(dd, then_branch, lowerBound, floor_value);
  if (!T) return NULL;
  E = setLowerBound_impl(dd, else_branch, lowerBound, floor_value);
  if (!E) {
    Cudd_RecursiveDeref(dd, T);
    return NULL;
  }

  if (T == E) {
    Cudd_RecursiveDeref(dd, E); // drop extra ref
    return T; // already protected
  }

  if (T == then_branch && E == else_branch) {
    Cudd_RecursiveDeref(dd, T);
    Cudd_RecursiveDeref(dd, E);
    Cudd_Ref(current_node);
    return current_node;
  }

  DdNode *var = Cudd_addIthVar(dd, Cudd_NodeReadIndex(reg));
  if (!var) {
    Cudd_RecursiveDeref(dd, T);
    Cudd_RecursiveDeref(dd, E);
    return NULL;
  }
  Cudd_Ref(var);
  res = Cudd_addIte(dd, var, T, E);
  Cudd_RecursiveDeref(dd, var);
  Cudd_RecursiveDeref(dd, T);
  Cudd_RecursiveDeref(dd, E);
  if (!res) {
    return NULL;
  }
  Cudd_Ref(res);
  return res;
}


static DdNode *add_const_ref(DdManager *mgr, double value) {
  DdNode *c = Cudd_addConst(mgr, (CUDD_VALUE_TYPE)value);
  if (!c) return NULL;
  Cudd_Ref(c);
  return c;
}

static DdNode *add_apply_ref(DdManager *mgr, DdNode *a, DdNode *b, DD_AOP op) {
  if (!a || !b) return NULL;
  DdNode *res = Cudd_addApply(mgr, op, a, b);
  if (!res) return NULL;
  Cudd_Ref(res);
  return res; // caller manages a/b lifetimes
}

static DdNode *add_ite_ref(DdManager *mgr, DdNode *f, DdNode *g, DdNode *h) {
  if (!f || !g || !h) return NULL;
  DdNode *res = Cudd_addIte(mgr, f, g, h);
  if (!res) return NULL;
  Cudd_Ref(res);
  return res; // caller manages f/g/h lifetimes
}

static int build_decision_choice_list(environment *env, const int *bool_path, int bool_path_len, term_t out_list) {
  int ret;
  term_t list = PL_new_term_ref();
  term_t pair = PL_new_term_ref();
  term_t new_list = PL_new_term_ref();
  term_t mvar_term = PL_new_term_ref();
  term_t val_term = PL_new_term_ref();
  functor_t pair_functor = PL_new_functor(PL_new_atom("-"), 2);

  unsigned char *selected = (unsigned char *)calloc((size_t)env->boolVars, 1);
  if (!selected) return FALSE;
  for (int i = 0; i < bool_path_len; i++) {
    int idx = bool_path[i];
    if (idx >= 0 && idx < env->boolVars) selected[idx] = 1;
  }

  ret = PL_put_nil(list);
  if (ret != TRUE) {
    free(selected);
    return ret;
  }

  for (int m = env->nVars - 1; m >= 0; m--) {
    variable v = env->vars[m];
    if (!v.decision) continue;

    int chosen = v.nVal - 1;
    for (int k = 0; k < v.nVal - 1; k++) {
      int b = v.firstBoolVar + k;
      if (b >= 0 && b < env->boolVars && selected[b]) {
        chosen = k;
        break;
      }
    }

    ret = PL_put_integer(mvar_term, m);
    if (ret != TRUE) {
      free(selected);
      return ret;
    }
    ret = PL_put_integer(val_term, chosen);
    if (ret != TRUE) {
      free(selected);
      return ret;
    }
    ret = PL_cons_functor(pair, pair_functor, mvar_term, val_term);
    if (ret != TRUE) {
      free(selected);
      return ret;
    }
    ret = PL_cons_list(new_list, pair, list);
    if (ret != TRUE) {
      free(selected);
      return ret;
    }
    ret = PL_put_term(list, new_list);
    if (ret != TRUE) {
      free(selected);
      return ret;
    }
  }

  free(selected);
  return PL_unify(out_list, list);
}

static foreign_t compute_optimal_strategy_sp(term_t env_ref, term_t b_list, term_t u_list, term_t strategy_list, term_t util) {
  int ret, i = 0, len_array_of_parents = 0; // i: term count
  int *array_of_parents = NULL; // decision variable indices on the best path
  double current_util, value; // current_util: utility for one term; value: best leaf value
  DdNode *current_root = NULL;
  DdNode *add_sum = NULL;
  long double *suffix_remaining = NULL; // suffix_remaining[k] = sum of impacts from entries[k] to the end
  term_t head_bdd = PL_new_term_ref();
  term_t head_util = PL_new_term_ref();
  term_t bdd_list = PL_copy_term_ref(b_list); // local copy of BDD list
  term_t util_list = PL_copy_term_ref(u_list); // local copy of utility list
  term_t oututil = PL_new_term_ref();
  term_t list;
  size_t nutils; // number of utilities / terms expected
  term_entry *entries = NULL; // scaled ADDs paired with their impact (Im(u_i))
  environment *env;

  ret = PL_get_pointer(env_ref,(void **)&env);
  RETURN_IF_FAIL

  ret = PL_skip_list(bdd_list, 0, &nutils);
  if (ret != PL_LIST) return FALSE;

  list = PL_new_term_ref();

  entries = (term_entry *)calloc(nutils, sizeof(term_entry));
  if (!entries) return FALSE;

  while (PL_get_list(bdd_list, head_bdd, bdd_list) &&
         PL_get_list(util_list, head_util, util_list)) {
    if ((size_t)i >= nutils) goto fail;
    ret = PL_get_pointer(head_bdd, (void **)&current_root);
    RETURN_IF_FAIL
    ret = PL_get_float(head_util, &current_util);
    RETURN_IF_FAIL

    DdNode *root = probability_dd_int(env, current_root);
    if (!root) goto fail;
    DdNode *c = add_const_ref(env->mgr, current_util);
    if (!c) {
      Cudd_RecursiveDeref(env->mgr, root);
      goto fail;
    }
    DdNode *scaled = add_apply_ref(env->mgr, root, c, Cudd_addTimes);
    Cudd_RecursiveDeref(env->mgr, root);
    Cudd_RecursiveDeref(env->mgr, c);
    if (!scaled) goto fail;
    entries[i].term = scaled;
    i++;
  }

  if (!PL_get_nil(bdd_list) || !PL_get_nil(util_list)) goto fail; // safe check
  if (i == 0) goto fail;

  for (int k = 0; k < i; k++) {
    DdNode *maxNode = Cudd_addFindMax(env->mgr, entries[k].term);
    DdNode *minNode = Cudd_addFindMin(env->mgr, entries[k].term);
    if (!maxNode || !minNode) goto fail;
    entries[k].impact = Cudd_V(maxNode) - Cudd_V(minNode);
  }
  qsort(entries, (size_t)i, sizeof(term_entry), compare_term_entry);

  double prune_floor = -DBL_MAX;
  suffix_remaining = (long double *)calloc((size_t)i, sizeof(long double));
  if (!suffix_remaining) goto fail;
  for (int k = i - 1; k >= 0; k--) {
    suffix_remaining[k] = (long double)entries[k].impact;
    if (k + 1 < i) suffix_remaining[k] += suffix_remaining[k + 1];
  }

  add_sum = entries[0].term;
  entries[0].term = NULL;
  for (int k = 1; k < i; k++) {
    DdNode *maxNode = Cudd_addFindMax(env->mgr, add_sum);
    if (!maxNode) goto fail;
    double current_max = Cudd_V(maxNode);
    double remaining_impact = (double)suffix_remaining[k];
    double lb = nextafter(current_max - remaining_impact, -INFINITY); // keep pruning conservative at floating-point boundaries
    DdNode *pruned = setLowerBound(env->mgr, add_sum, lb, prune_floor);
    if (!pruned) goto fail;
    Cudd_RecursiveDeref(env->mgr, add_sum);
    add_sum = pruned;
    DdNode *lhs = add_sum;
    DdNode *rhs = entries[k].term;
    DdNode *next = add_apply_ref(env->mgr, lhs, rhs, Cudd_addPlus);
    if (!next) goto fail;
    if (rhs == lhs) {
      Cudd_RecursiveDeref(env->mgr, lhs);
    } else {
      Cudd_RecursiveDeref(env->mgr, lhs);
      Cudd_RecursiveDeref(env->mgr, rhs);
    }
    add_sum = next;
    entries[k].term = NULL;
  }

  DdNode *bestNode = Cudd_addFindMax(env->mgr, add_sum);
  if (bestNode == NULL) {
    Cudd_RecursiveDeref(env->mgr, add_sum);
    free(entries);
    printf("Cudd_addFindMax returned NULL\n");
    return FALSE;
  }

  value = Cudd_V(bestNode);
  array_of_parents = (int *)malloc(sizeof(int));
  if (!array_of_parents) goto fail;

  ret = find_path(add_sum, value, &array_of_parents, &len_array_of_parents);
  if (ret != 1) goto fail;

  ret = build_decision_choice_list(env, array_of_parents, len_array_of_parents, list);
  RETURN_IF_FAIL
  ret = PL_put_float(oututil, value);
  RETURN_IF_FAIL

  free(array_of_parents);
  Cudd_RecursiveDeref(env->mgr, add_sum);
  free(suffix_remaining);
  free(entries);
  return (PL_unify(list, strategy_list) && PL_unify(util, oututil));

fail:
  if (entries) {
    for (size_t k = 0; k < nutils; k++) {
      if (entries[k].term) Cudd_RecursiveDeref(env->mgr, entries[k].term);
    }
    free(entries);
  }
  if (array_of_parents) free(array_of_parents);
  if (add_sum) Cudd_RecursiveDeref(env->mgr, add_sum);
  free(suffix_remaining);
  return FALSE;
}


static foreign_t probability_dd(term_t env_ref, term_t bdd_ref, term_t add_out) {
  int ret;
  term_t out;
  DdNode *bdd = NULL, *root = NULL;
  environment *env;

  ret = PL_get_pointer(env_ref,(void **)&env);
  RETURN_IF_FAIL

  ret = PL_get_pointer(bdd_ref,(void **)&bdd);
  RETURN_IF_FAIL

  out = PL_new_term_ref();
  root = probability_dd_int(env, bdd);
  if(root == NULL) {
    printf("computed add NULL\n");
    return -1;
  }
  ret = PL_put_pointer(out,(void *)root);
  RETURN_IF_FAIL

  return(PL_unify(out,add_out));
}

DdNode* probability_dd_int(environment *env, DdNode *current_node) {
  autodyn_guard_t guard;
  autodyn_guard_begin(env->mgr, &guard);
  DdNode *res = probability_dd_int_impl(env, current_node);
  autodyn_guard_end(env->mgr, &guard);
  return res;
}

static DdNode* probability_dd_int_impl(environment *env, DdNode *current_node) {
  if(current_node == NULL) {
    printf("NODE NULL\n");
    return NULL;
  }

  DdNode *reg = Cudd_Regular(current_node);
  int is_comp = Cudd_IsComplement(current_node);

  if (Cudd_IsConstant(reg)) {
    double v = Cudd_V(reg);
    if (is_comp) v = 1.0 - v;
    return add_const_ref(env->mgr, v);
  }

  int index = (int)Cudd_NodeReadIndex(reg);
  if(index >= env->boolVars) {
    printf("INDEX BOOL ERROR: index=%d boolVars=%d mgrSize=%d node=%p\n",
           index, env->boolVars, Cudd_ReadSize(env->mgr), (void *)current_node);
    return NULL;
  }

  int indexMultival = env->bVar2mVar[index];
  if(indexMultival >= env->nVars) {
    printf("INDEX MULTIVAR ERROR\n");
    return NULL;
  }

  DdNode *addh = probability_dd_int_impl(env, Cudd_T(reg));
  if (!addh) return NULL;
  DdNode *addl = probability_dd_int_impl(env, Cudd_E(reg));
  if (!addl) {
    Cudd_RecursiveDeref(env->mgr, addh);
    return NULL;
  }

  if (addh == addl) { // addh and addl are the same, therefore merged
    Cudd_RecursiveDeref(env->mgr, addl);
    if (!is_comp) return addh;
    DdNode *one = add_const_ref(env->mgr, 1.0);
    if (!one) {
      Cudd_RecursiveDeref(env->mgr, addh);
      return NULL;
    }
    DdNode *res = add_apply_ref(env->mgr, one, addh, Cudd_addMinus);
    Cudd_RecursiveDeref(env->mgr, one);
    Cudd_RecursiveDeref(env->mgr, addh);
    if (!res) {
      return NULL;
    }
    return res;
  }

  if(env->vars[indexMultival].decision == 1) {
    DdNode *var = Cudd_addIthVar(env->mgr, index);
    if (!var) {
      Cudd_RecursiveDeref(env->mgr, addh);
      Cudd_RecursiveDeref(env->mgr, addl);
      return NULL;
    }
    Cudd_Ref(var);
    DdNode *base = add_ite_ref(env->mgr, var, addh, addl);
    Cudd_RecursiveDeref(env->mgr, var);
    Cudd_RecursiveDeref(env->mgr, addh);
    Cudd_RecursiveDeref(env->mgr, addl);
    if (!base) {
      return NULL;
    }
    if (!is_comp) return base;
    DdNode *one = add_const_ref(env->mgr, 1.0);
    if (!one) {
      Cudd_RecursiveDeref(env->mgr, base);
      return NULL;
    }
    DdNode *res = add_apply_ref(env->mgr, one, base, Cudd_addMinus);
    Cudd_RecursiveDeref(env->mgr, one);
    Cudd_RecursiveDeref(env->mgr, base);
    if (!res) {
      return NULL;
    }
    return res;
  } else if(env->vars[indexMultival].decision == 0) {
    double p = env->probs[index];
    if(p < 0.0 || p > 1.0) {
      printf("ERROR IN PROBABILITY\n");
      Cudd_RecursiveDeref(env->mgr, addh);
      Cudd_RecursiveDeref(env->mgr, addl);
      return NULL;
    }
    if (p == 1.0) {
      Cudd_RecursiveDeref(env->mgr, addl);
      if (!is_comp) return addh;
      DdNode *one = add_const_ref(env->mgr, 1.0);
      if (!one) {
        Cudd_RecursiveDeref(env->mgr, addh);
        return NULL;
      }
      DdNode *res = add_apply_ref(env->mgr, one, addh, Cudd_addMinus);
      Cudd_RecursiveDeref(env->mgr, one);
      Cudd_RecursiveDeref(env->mgr, addh);
      if (!res) {
        return NULL;
      }
      return res;
    }
    if (p == 0.0) {
      Cudd_RecursiveDeref(env->mgr, addh);
      if (!is_comp) return addl;
      DdNode *one = add_const_ref(env->mgr, 1.0);
      if (!one) {
        Cudd_RecursiveDeref(env->mgr, addl);
        return NULL;
      }
      DdNode *res = add_apply_ref(env->mgr, one, addl, Cudd_addMinus);
      Cudd_RecursiveDeref(env->mgr, one);
      Cudd_RecursiveDeref(env->mgr, addl);
      if (!res) {
        return NULL;
      }
      return res;
    }
    DdNode *pnode = add_const_ref(env->mgr, p);
    if (!pnode) {
      Cudd_RecursiveDeref(env->mgr, addh);
      Cudd_RecursiveDeref(env->mgr, addl);
      return NULL;
    }
    DdNode *pa = add_apply_ref(env->mgr, pnode, addh, Cudd_addTimes);
    if (!pa) {
      Cudd_RecursiveDeref(env->mgr, pnode);
      Cudd_RecursiveDeref(env->mgr, addh);
      Cudd_RecursiveDeref(env->mgr, addl);
      return NULL;
    }
    Cudd_RecursiveDeref(env->mgr, pnode);
    Cudd_RecursiveDeref(env->mgr, addh);
    DdNode *p1node = add_const_ref(env->mgr, 1.0 - p);
    if (!p1node) {
      Cudd_RecursiveDeref(env->mgr, pa);
      Cudd_RecursiveDeref(env->mgr, addl);
      return NULL;
    }
    DdNode *pb = add_apply_ref(env->mgr, p1node, addl, Cudd_addTimes);
    if (!pb) {
      Cudd_RecursiveDeref(env->mgr, p1node);
      Cudd_RecursiveDeref(env->mgr, addl);
      Cudd_RecursiveDeref(env->mgr, pa);
      return NULL;
    }
    Cudd_RecursiveDeref(env->mgr, p1node);
    Cudd_RecursiveDeref(env->mgr, addl);
    DdNode *base = add_apply_ref(env->mgr, pa, pb, Cudd_addPlus);
    Cudd_RecursiveDeref(env->mgr, pa);
    Cudd_RecursiveDeref(env->mgr, pb);
    if (!base) {
      return NULL;
    }
    if (!is_comp) return base;
    DdNode *one = add_const_ref(env->mgr, 1.0);
    if (!one) {
      Cudd_RecursiveDeref(env->mgr, base);
      return NULL;
    }
    DdNode *res = add_apply_ref(env->mgr, one, base, Cudd_addMinus);
    Cudd_RecursiveDeref(env->mgr, one);
    Cudd_RecursiveDeref(env->mgr, base);
    if (!res) {
      return NULL;
    }
    return res;
  } else {
    printf("ERROR IN VAR\n");
    Cudd_RecursiveDeref(env->mgr, addh);
    Cudd_RecursiveDeref(env->mgr, addl);
    return NULL;
  }
}

static foreign_t add_prod(term_t env_ref, term_t add_in, term_t value, term_t add_out) {
  int ret;
  double current_value;
  DdNode *add_const, *current_add, *add_ret;
  environment *env;
  term_t out;

  ret = PL_get_pointer(env_ref,(void **)&env);
  RETURN_IF_FAIL

  ret = PL_get_float(value,&current_value);
  RETURN_IF_FAIL

  ret = PL_get_pointer(add_in,(void **)&current_add);
  RETURN_IF_FAIL

  // printf("Current value: %lf\n",current_value);

  add_const = add_const_ref(env->mgr, current_value);
  if (!add_const) return FALSE;
  add_ret = add_apply_ref(env->mgr, add_const, current_add, Cudd_addTimes);
  if (!add_ret) {
    Cudd_RecursiveDeref(env->mgr, add_const);
    return FALSE;
  }
  Cudd_RecursiveDeref(env->mgr, add_const);
  Cudd_RecursiveDeref(env->mgr, current_add);

  out = PL_new_term_ref();
  ret = PL_put_pointer(out,(void *)add_ret);
  RETURN_IF_FAIL

  return(PL_unify(out,add_out)); 
}

static foreign_t add_sum(term_t env_ref, term_t add_A, term_t add_B, term_t add_out) {
  int ret;
  term_t out;
  environment *env;
  DdNode *addA, *addB, *add_ret;

  ret = PL_get_pointer(env_ref,(void **)&env);
  RETURN_IF_FAIL
  ret = PL_get_pointer(add_A,(void **)&addA);
  RETURN_IF_FAIL
  ret = PL_get_pointer(add_B,(void **)&addB);
  RETURN_IF_FAIL

  add_ret = add_apply_ref(env->mgr, addA, addB, Cudd_addPlus);
  if (!add_ret) return FALSE;
  Cudd_RecursiveDeref(env->mgr, addA);
  if (addB != addA) Cudd_RecursiveDeref(env->mgr, addB);

  out = PL_new_term_ref();
  ret = PL_put_pointer(out,(void *)add_ret);
  RETURN_IF_FAIL

  return(PL_unify(out,add_out));  
}


static foreign_t ret_strategy(term_t env_ref, term_t add, term_t strategy_list, term_t util) {
  int ret;
  int *array_of_parents = NULL; // array containing the indexes of the parents
  int len_array_of_parents = 0;
  double value = -1;
  term_t opt_util;
  term_t list;
  DdNode *root = NULL, *bestNode = NULL;
  environment *env;

  ret = PL_get_pointer(env_ref,(void **)&env);
  RETURN_IF_FAIL
  ret = PL_get_pointer(add,(void **)&root);
  RETURN_IF_FAIL
  
  list = PL_new_term_ref();
  opt_util = PL_new_term_ref();

  if(DEBUG) {  
    printf("Cudd_CheckZeroRef: %d\n",Cudd_CheckZeroRef(env->mgr));
    printf("Max memory: %lu\n",Cudd_ReadMaxMemory(env->mgr));
    printf("DdManager nodes: %ld | ", Cudd_ReadNodeCount(env->mgr)); /*Reports the number of live nodes in BDDs and ADDs*/
    printf("DdManager vars: %d | ", Cudd_ReadSize(env->mgr) ); /*Returns the number of BDD variables in existance*/
    printf("DdNode nodes: %d | ", Cudd_DagSize(root)); /*Reports the number of nodes in the BDD*/
    printf("DdNode vars: %d | ", Cudd_SupportSize(env->mgr, root) ); /*Returns the number of variables in the BDD*/
    printf("DdManager reorderings: %d | ", Cudd_ReadReorderings(env->mgr) ); /*Returns the number of times reordering has occurred*/
    printf("DdManager memory: %ld |\n\n", Cudd_ReadMemoryInUse(env->mgr) ); /*Returns the memory in use by the manager measured in bytes*/
    printf("Cudd_CountLeaves: %d |\n\n", Cudd_CountLeaves(root)); /*Returns the memory in use by the manager measured in bytes*/
    debug_cudd_env(env,0);
    Cudd_PrintDebug(env->mgr, root, 2, 4);
  }

  array_of_parents = malloc(sizeof(int));
  bestNode = Cudd_addFindMax(env->mgr,root);

  // check if found
  if(bestNode == NULL) {
    // no solution found -> return empty list and -1 as utility
    ret = PL_put_integer(opt_util,(long)-1);
    RETURN_IF_FAIL
    if(array_of_parents) {
      free(array_of_parents);
    }
    if(root) {
      Cudd_RecursiveDeref(env->mgr,root);
    }
    return PL_unify(opt_util,util);
  } 
  else {
    // find path: root -> terminal
    value = Cudd_V(bestNode);
    ret = find_path(root,value,&array_of_parents,&len_array_of_parents);
    if(ret != 1) {
      goto fail;
    }
    ret = build_decision_choice_list(env, array_of_parents, len_array_of_parents, list);
    RETURN_IF_FAIL
    ret = PL_put_float(opt_util,value);
    RETURN_IF_FAIL
  }

  if(array_of_parents) {
    free(array_of_parents);
  }

  Cudd_RecursiveDeref(env->mgr,root);

  return(PL_unify(list,strategy_list) && (PL_unify(opt_util,util)));

fail:
  if(array_of_parents) {
    free(array_of_parents);
  }
  if(root) {
    Cudd_RecursiveDeref(env->mgr,root);
  }
  return FALSE;
}

int find_path(DdNode *node, double value, int **array, int *len) { //CUDD nodes do not have parent links :(
  int ret;

  if(node == NULL) {
    return 0;
  }
  if(Cudd_IsConstant(node)) {
    if(Cudd_V(node) == value) {
      return 1;
    }
    else {
      return 0;
    }
  }
  // if found int the branch then (1) -> chosen, add one element
  ret = find_path(Cudd_T(node),value, array, len);
  if(ret < 0) {
    return ret;
  }
  if(ret == 1) {
      int *tmp = realloc(*array, ((*len)+1)*sizeof(int));
      if(!tmp) {
        return -1;
      }
      *array = tmp;
      (*array)[*len] = Cudd_NodeReadIndex(node);
      *len = (*len) + 1;
    return 1;
  }
  ret = find_path(Cudd_E(node),value, array, len);
  if(ret != 0) {
    return ret;
  }
  return 0;
}


static foreign_t add_abd_var(term_t arg1,term_t arg2,term_t arg3,term_t arg4)
{
  term_t out,head,probTerm;
  variable * v;
  int i,ret,nRules;
  double p,p0;
  environment * env;
  size_t lenProbs;

  head=PL_new_term_ref();
  out=PL_new_term_ref();
  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL
  env->nVars=env->nVars+1;
  env->vars=(variable *) realloc(env->vars,env->nVars * sizeof(variable));

  env->n_abd=env->n_abd+1;

  v=&env->vars[env->nVars-1];

  v->abducible=1;
  v->query=0;
  v->decision = 0;
  probTerm=PL_copy_term_ref(arg2);
  ret=PL_skip_list(probTerm,0,&lenProbs);
  if (ret!=PL_LIST) return FALSE;
  v->nVal=lenProbs;

  ret=PL_get_integer(arg3,&v->nRule);
  RETURN_IF_FAIL
  nRules=env->nRules;
  if (v->nRule>=nRules)
  {
    env->rules=(int *)  realloc(env->rules,((v->nRule+1)* sizeof(int)));
    for (i=nRules;i<v->nRule;i++)
      env->rules[i]=0;
    env->rules[v->nRule]=lenProbs;
    env->nRules=v->nRule+1;
  }
  env->n_abd_boolVars=env->n_abd_boolVars+v->nVal-1;

  v->firstBoolVar=env->boolVars;
  env->probs=(double *) realloc(env->probs,(((env->boolVars+v->nVal-1)* sizeof(double))));
  env->bVar2mVar=(int *) realloc(env->bVar2mVar,((env->boolVars+v->nVal-1)* sizeof(int)));

  p0=1;
  for (i=0;i<v->nVal-1;i++)
  {
    ret=PL_get_list(probTerm,head,probTerm);
    RETURN_IF_FAIL
    ret=PL_get_float(head,&p);
    RETURN_IF_FAIL
    env->bVar2mVar[env->boolVars+i]=env->nVars-1;
    env->probs[env->boolVars+i]=p/p0;
    p0=p0*(1-p/p0);
  }
  env->boolVars=env->boolVars+v->nVal-1;
  env->rules[v->nRule]= v->nVal;
  ret=PL_put_integer(out,env->nVars-1);
  RETURN_IF_FAIL

  return(PL_unify(out,arg4));
}

static foreign_t equality(term_t arg1,term_t arg2,term_t arg3, term_t arg4)
{
  term_t out;
  int varIndex;
  int value;
  int i,ret;
  variable v;
  DdNode * node, * tmp,*var;
  environment * env;

  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL

  ret=PL_get_integer(arg2,&varIndex);
  RETURN_IF_FAIL
  ret=PL_get_integer(arg3,&value);
  RETURN_IF_FAIL
  v=env->vars[varIndex];
  i=v.firstBoolVar;
  tmp=Cudd_ReadOne(env->mgr);
  Cudd_Ref(tmp);
  node=NULL;
  if (v.query)
  {
    var=Cudd_bddIthVar(env->mgr,v.firstBoolVar+value);
    node=Cudd_bddAnd(env->mgr,tmp,var);
    Cudd_Ref(node);
  }
  else
  {
    for (i=v.firstBoolVar;i<v.firstBoolVar+value;i++)
    {
      var=Cudd_bddIthVar(env->mgr,i);
      node=Cudd_bddAnd(env->mgr,tmp,Cudd_Not(var));
      Cudd_Ref(node);
      Cudd_RecursiveDeref(env->mgr,tmp);
      tmp=node;
    }
    if (!(value==v.nVal-1))
    {
      var=Cudd_bddIthVar(env->mgr,v.firstBoolVar+value);
      node=Cudd_bddAnd(env->mgr,tmp,var);
      Cudd_Ref(node);
      Cudd_RecursiveDeref(env->mgr,tmp);
    }
  }
  out=PL_new_term_ref();
  ret=PL_put_pointer(out,(void *)node);
  RETURN_IF_FAIL
  return(PL_unify(out,arg4));
}

static foreign_t one(term_t arg1, term_t arg2)
{
  term_t out;
  DdNode * node;
  environment *env;
  int res,ret;

  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL

  node =  Cudd_ReadOne(env->mgr);
  Cudd_Ref(node);
  out=PL_new_term_ref();
  ret=PL_put_pointer(out,(void *) node);
  RETURN_IF_FAIL
  res=PL_unify(out,arg2);
  return res;

//  return(PL_unify(out,arg2));
}

static foreign_t zero(term_t arg1, term_t arg2)
{
  term_t out;
  DdNode * node;
  environment *env;
  int ret;

  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL

  node = Cudd_ReadLogicZero(env->mgr);
  Cudd_Ref(node);
  out=PL_new_term_ref();
  ret=PL_put_pointer(out,(void *) node);
  RETURN_IF_FAIL
  return(PL_unify(out,arg2));
}

// arg1 (env ref) unused
static foreign_t bdd_not(term_t arg1,term_t arg2, term_t arg3)
{
  term_t out;
  DdNode * node;
  int ret;

  ret=PL_get_pointer(arg2,(void **)&node);
  RETURN_IF_FAIL
  node=Cudd_Not(node);
  out=PL_new_term_ref();
  ret=PL_put_pointer(out,(void *) node);
  RETURN_IF_FAIL
  return(PL_unify(out,arg3));
}

static foreign_t and(term_t arg1,term_t arg2,term_t arg3, term_t arg4)
{
  term_t out;
  DdNode * node1, *node2,*nodeout;
  environment *env;
  int res,ret;
  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL

  ret=PL_get_pointer(arg2,(void **)&node1);
  RETURN_IF_FAIL
  ret=PL_get_pointer(arg3,(void **)&node2);
  RETURN_IF_FAIL
  nodeout=Cudd_bddAnd(env->mgr,node1,node2);
  Cudd_Ref(nodeout);
  out=PL_new_term_ref();
  ret=PL_put_pointer(out,(void *) nodeout);
  RETURN_IF_FAIL
  res=PL_unify(out,arg4);
  return res;
}

static foreign_t or(term_t arg1,term_t arg2,term_t arg3, term_t arg4)
{
  term_t out;
  DdNode * node1, *node2,*nodeout;
  environment *env;
  int ret;

  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL

  ret=PL_get_pointer(arg2,(void **)&node1);
  RETURN_IF_FAIL
  ret=PL_get_pointer(arg3,(void **)&node2);
  RETURN_IF_FAIL

  nodeout=Cudd_bddOr(env->mgr,node1,node2);
  Cudd_Ref(nodeout);
  out=PL_new_term_ref();
  ret=PL_put_pointer(out,(void *) nodeout);
  RETURN_IF_FAIL
  return(PL_unify(out,arg4));
}

/*
static int garbage_collect(void)
{
  YAP_Term arg1,arg2,out;
  YAP_Int nodes,clearCache;

  arg1=YAP_ARG1;
  arg2=YAP_ARG2;
  clearCache=YAP_IntOfTerm(arg1);
  nodes=(YAP_Int)cuddGarbageCollect(mgr_ex[ex],clearCache);
  out=YAP_MkIntTerm(nodes);
  return(YAP_Unify(out,arg2));
}

static int bdd_to_add(void)
{
  YAP_Term arg1,arg2,out;
  DdNode * node1,*node2;

  arg1=YAP_ARG1;
  arg2=YAP_ARG2;
  node1=(DdNode *)YAP_IntOfTerm(arg1);
  node2= Cudd_BddToAdd(mgr_ex[ex],node1);
  out=YAP_MkIntTerm((YAP_Int) node2);
  return(YAP_Unify(out,arg2));
}
*/
static foreign_t create_dot(term_t arg1, term_t arg2, term_t arg3)
{
  DdNode * node;
  environment *env;
  char *filename;
  FILE * file;
  int ret;
  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL
  ret=PL_get_pointer(arg2,(void **)&node);
  RETURN_IF_FAIL
  ret=PL_get_file_name(arg3,&filename,0);
  RETURN_IF_FAIL
  file = open_file(filename, "w");
  write_dot(env,node,file);
  fclose(file);
  return TRUE;
}

static foreign_t create_dot_string(term_t arg1, term_t arg2, term_t arg3)
{
  term_t out;
  DdNode * node;
  environment *env;
  FILE * file;
  char *buffer=NULL;
  int ret;

  ret=PL_get_pointer(arg1,(void **)&env);
  RETURN_IF_FAIL
  ret=PL_get_pointer(arg2,(void **)&node);
  RETURN_IF_FAIL
  out=PL_new_term_ref();

#ifndef _WIN32
  file=tmpfile();
#else
  char filename[MAX_PATH];
  GetTempFileName(".","temp",0,filename);
  file = fopen(filename,"w+bTD");
#endif
  if (file==NULL) {perror("Error in temporary file opening");}
  write_dot(env,node,file);

  if (fseek(file, 0L, SEEK_END) == 0) {
    /* Get the size of the file. */
    long bufsize = ftell(file);
    if (bufsize == -1) { perror("Error in getting the size of the temporary file");}

      /* Allocate our buffer to that size. */
        buffer = malloc(sizeof(char) * (bufsize + 1));

        /* Go back to the start of the file. */
        if (fseek(file, 0L, SEEK_SET) != 0) { perror("Error going back to the start of the file");}

        /* Read the entire file into memory. */
        size_t newLen = fread(buffer, sizeof(char), bufsize, file);
        if ( ferror( file ) != 0 ) {
            perror("Error reading file");
        } else {
            buffer[newLen++] = '\0'; /* Just to be safe. */
        }
  }
  fclose(file);
  ret=PL_put_string_chars(out,buffer);
  RETURN_IF_FAIL
  return(PL_unify(out,arg3));
}

void write_dot(environment * env, DdNode * bdd, FILE * file)
{
  const char * onames[]={"Out"};
  char ** inames;
  int i,b,index,nv;
  variable v;
  char numberVar[11],numberBit[11];
  inames= (char **) malloc(sizeof(char *)*(env->boolVars));
  index=0;
  for (i=0;i<env->nVars;i++)
  {
    v=env->vars[i];
      if (v.query)
        nv=v.nVal;
      else
        nv=v.nVal-1;
    for (b=0;b<nv;b++)
    {
      inames[b+index]=(char *) malloc(sizeof(char)*20);
      strcpy(inames[b+index],"X");
      sprintf(numberVar,"%d",i);
      strcat(inames[b+index],numberVar);
      strcat(inames[b+index],"_");
      sprintf(numberBit,"%d",b);
      strcat(inames[b+index],numberBit);
    }
    index=index+nv;
  }
  Cudd_DumpDot(env->mgr,1,&bdd,(const char * const *)inames,(const char * const *)onames,file);
  index=0;
  for (i=0;i<env->nVars;i++)
  {
    v=env->vars[i];
    if (v.query)
      nv=v.nVal;
    else
      nv=v.nVal-1;
    for (b=0;b<nv;b++)
    {
      free(inames[b+index]);
    }
    index=index+nv;
  }
  free(inames);
}

/*
static int rec_deref(void)
{
  YAP_Term arg1;
  DdNode * node;

  arg1=YAP_ARG1;
  node=(DdNode *) YAP_IntOfTerm(arg1);
  Cudd_RecursiveDeref(mgr_ex[ex], node);
  return 1;
}

*/

double ProbPath(example_data * ex_d,DdNode *node, int nex)
{
  int index,mVarIndex,pos,position;//,boolVarIndex;
  variable v;
  double res;
  double p,pt,pf,BChild0e,BChild1e,BChild0o,BChild1o,e0,e1;
  double *value_p, * value_p_e, *value_p_o,** eta_rule;
  DdNode *nodekey,*T,*F;

  // printf("node %p comp %d\n",node,comp );
  if (Cudd_IsConstant(node))
  {
    // printf("constant\n");
    return 1.0;
  }
  else
  {
    nodekey=Cudd_Regular(node);
    value_p=get_value(ex_d->nodesB,nodekey);
    if (value_p!=NULL)
    {
      // printf("found %f\n", *value_p);
      return *value_p;
    }
    else
    {
      index=Cudd_NodeReadIndex(node);
      p=ex_d->env[nex].probs[index];
      T = Cudd_T(node);
      F = Cudd_E(node);
      pf=ProbPath(ex_d,F,nex);
      pt=ProbPath(ex_d,T,nex);
      // printf("pt %f pf %f\n",pt,pf );
      if (Cudd_IsComplement(F))
        pf=1.0-pf;
      // printf("pt %f pf %f\n",pt,pf );

      BChild0e=pf*(1-p);
      BChild0o=(1-pf)*(1-p);
      BChild1e=pt*p;
      BChild1o=(1-pt)*p;
      value_p_e=get_value(ex_d->nodesFE,nodekey);
      value_p_o=get_value(ex_d->nodesFO,nodekey);
      e0 = (*value_p_e)*BChild0e+(*value_p_o)*BChild0o;
      e1 = (*value_p_e)*BChild1e+(*value_p_o)*BChild1o;
    //  printf("e node %p %f %f %f %f\n",node,*value_p_e,*value_p_o,e0,e1 );
      mVarIndex=ex_d->env[nex].bVar2mVar[index];
      v=ex_d->env[nex].vars[mVarIndex];
      pos=index-v.firstBoolVar;
      if (ex_d->tunable_rules[v.nRule])
      {
        eta_rule=ex_d->eta_temp[v.nRule];
        eta_rule[pos][0]=eta_rule[pos][0]+e0;
        eta_rule[pos][1]=eta_rule[pos][1]+e1;
      }
      res=BChild0e+BChild1e;
      add_node(ex_d->nodesB,nodekey,res);
      position=Cudd_ReadPerm(ex_d->env[nex].mgr,index);
      position=position+1;
//      boolVarIndex=Cudd_ReadInvPerm(ex_d->env[nex].mgr,position);//Returns the index of the variable currently in the i-th position of the order.
      if (position<ex_d->env[nex].boolVars)
      {
        ex_d->sigma[position]=ex_d->sigma[position]+e0+e1;
      }
      if(!Cudd_IsConstant(T))
      {
        index=Cudd_NodeReadIndex(T);
        position=Cudd_ReadPerm(ex_d->env[nex].mgr,index);
        ex_d->sigma[position]=ex_d->sigma[position]-e1;
      }

      if(!Cudd_IsConstant(F))
      {
        index=Cudd_NodeReadIndex(F);
        position=Cudd_ReadPerm(ex_d->env[nex].mgr,index);
        ex_d->sigma[position]=ex_d->sigma[position]-e0;
      }

      return res;
    }
  }
}




void Forward(example_data * ex_d,DdNode *root, int nex)
{
  DdNode *** nodesToVisit;
  int * NnodesToVisit,comp;
  double FrootE,FrootO;

  environment env;
  int i,j;
  env=ex_d->env[nex];
  comp=Cudd_IsComplement(root);
  if (comp)
  {
    FrootE=0.0;
    FrootO=1.0;
  }
  else
  {
    FrootE=1.0;
    FrootO=0.0;
  }
  add_node(ex_d->nodesFE,Cudd_Regular(root),FrootE);
  add_node(ex_d->nodesFO,Cudd_Regular(root),FrootO);
  if (env.boolVars)
  {
    nodesToVisit= (DdNode ***)malloc(sizeof(DdNode **)* env.boolVars);
    NnodesToVisit= (int *)malloc(sizeof(int)* env.boolVars);
    nodesToVisit[0]=(DdNode **)malloc(sizeof(DdNode *));
    nodesToVisit[0][0]=root;
    NnodesToVisit[0]=1;
    for(i=1;i<env.boolVars;i++)
    {
      nodesToVisit[i]=NULL;
      NnodesToVisit[i]=0;
    }
    for(i=0;i<env.boolVars;i++)
    {
      for(j=0;j<NnodesToVisit[i];j++)
      UpdateForward(ex_d,nodesToVisit[i][j],nex,nodesToVisit,NnodesToVisit);
    }
    for(i=0;i<env.boolVars;i++)
    {
      free(nodesToVisit[i]);
    }
    free(nodesToVisit);
    free(NnodesToVisit);
  }
}

void UpdateForward(example_data *ex_d,DdNode *node, int nex,
  DdNode *** nodesToVisit, int * NnodesToVisit)
{
  int index,position;
  DdNode *T,*E,*nodereg;
  double *value_p_E,*value_p_O,*value_p_T_E,*value_p_T_O,
    *value_p_F_E,*value_p_F_O,p,value_par_E,value_par_O;

// printf("F node %p comp %d\n",node,Cudd_IsComplement(node) );
  if (Cudd_IsConstant(node))
  {
    return;
  }
  else
  {
    index=Cudd_NodeReadIndex(node);
    p=ex_d->env[nex].probs[index];
    nodereg=Cudd_Regular(node);
    value_p_E=get_value(ex_d->nodesFE,nodereg);
    value_p_O=get_value(ex_d->nodesFO,nodereg);
    if (value_p_E== NULL)
    {
      printf("Error\n");
      return;
    }
    else
    {
      T = Cudd_T(node);
      E = Cudd_E(node);
      if (!Cudd_IsConstant(T))
      {
        value_p_T_E=get_value(ex_d->nodesFE,T);
        value_p_T_O=get_value(ex_d->nodesFO,T);
        if (value_p_T_E!= NULL)
        {
           *value_p_T_E= *value_p_T_E+*value_p_E*p;
           *value_p_T_O= *value_p_T_O+*value_p_O*p;
          // printf("update f t %p %f %f %f %f\n",T,*value_p_T_E,
          //  *value_p_E*p,*value_p_T_O,*value_p_O*p);
        }
        else
        {
          // printf("new f t %p %f %f \n",T,*value_p_E*p,*value_p_O*p );

          add_or_replace_node(ex_d->nodesFE,Cudd_Regular(T),*value_p_E*p);
          add_or_replace_node(ex_d->nodesFO,Cudd_Regular(T),*value_p_O*p);
          index=Cudd_NodeReadIndex(T);
          position=Cudd_ReadPerm(ex_d->env[nex].mgr,index);
          nodesToVisit[position]=(DdNode **)realloc(nodesToVisit[position],
	    (NnodesToVisit[position]+1)* sizeof(DdNode *));
          nodesToVisit[position][NnodesToVisit[position]]=T;
          NnodesToVisit[position]=NnodesToVisit[position]+1;
        }
      }
      if (!Cudd_IsConstant(E))
      {
        value_p_F_E=get_value(ex_d->nodesFE,Cudd_Regular(E));
        value_p_F_O=get_value(ex_d->nodesFO,Cudd_Regular(E));
        // if (Cudd_IsComplement(E))
        //   value_par=1 - *value_p;
        // else
        //   value_par= *value_p;
        // if (Cudd_IsComplement(E))
        //   p=1 -p;
        if (Cudd_IsComplement(E))
        {
          value_par_E= *value_p_O;
          value_par_O= *value_p_E;
        }
        else
        {
          value_par_E= *value_p_E;
          value_par_O= *value_p_O;
        }
        // printf("f child %d %f %f\n",Cudd_IsComplement(E),value_par_E,value_par_O );

        if (value_p_F_E!= NULL)
        {

          *value_p_F_E= *value_p_F_E+value_par_E*(1-p);
          *value_p_F_O= *value_p_F_O+value_par_O*(1-p);
          // printf("update f f %p %f %f %f %f\n",E,*value_p_F_E, value_par_E*(1-p),
        // *value_p_F_O, value_par_O*(1-p));
        }
        else
        {
                              // printf("new f f %p %f\n",E,value_par_E*(1-p) );

          add_or_replace_node(ex_d->nodesFE,Cudd_Regular(E),value_par_E*(1-p));
          add_or_replace_node(ex_d->nodesFO,Cudd_Regular(E),value_par_O*(1-p));
          index=Cudd_NodeReadIndex(E);
          position=Cudd_ReadPerm(ex_d->env[nex].mgr,index);
          nodesToVisit[position]=(DdNode **)realloc(nodesToVisit[position],
	    (NnodesToVisit[position]+1)* sizeof(DdNode *));
          nodesToVisit[position][NnodesToVisit[position]]=E;
          NnodesToVisit[position]=NnodesToVisit[position]+1;
        }
      }
      return;
    }
  }
}




double GetOutsideExpe(example_data * ex_d,DdNode *root,double ex_prob, int nex)
{
  int i,j,mVarIndex,bVarIndex,firstBoolVarOfRule,nRule;
  double **eta_rule;
  double theta,rootProb, T=0;


  ex_d->sigma=(double *)malloc(ex_d->env[nex].boolVars * sizeof(double));

  for (j=0; j<ex_d->env[nex].boolVars; j++)
  {
    ex_d->sigma[j]=0;
  }
  for (j=0; j<ex_d->nRules; j++)
  {
    if (ex_d->tunable_rules[j])
      for (i=0; i<ex_d->rules[j]-1; i++)
      {
        ex_d->eta_temp[j][i][0]=0;
        ex_d->eta_temp[j][i][1]=0;
      }
  }
  rootProb=ProbPath(ex_d,root,nex);
  if (Cudd_IsComplement(root))
    rootProb=1.0-rootProb;
  if (rootProb>0.0)
  {
    for (j=0; j<ex_d->env[nex].boolVars; j++)
    {
      T += ex_d->sigma[j];
      bVarIndex=Cudd_ReadInvPerm(ex_d->env[nex].mgr,j);
      if (bVarIndex==-1)
      {
        bVarIndex=j;
      }

      mVarIndex=ex_d->env[nex].bVar2mVar[bVarIndex];

      firstBoolVarOfRule=ex_d->env[nex].vars[mVarIndex].firstBoolVar;
      i=bVarIndex-firstBoolVarOfRule;
      theta=ex_d->env[nex].probs[bVarIndex];
      nRule=ex_d->env[nex].vars[mVarIndex].nRule;
      if (ex_d->tunable_rules[nRule])
      {
        eta_rule=ex_d->eta_temp[nRule];
        eta_rule[i][0]=eta_rule[i][0]+T*(1-theta);
        eta_rule[i][1]=eta_rule[i][1]+T*theta;
      }
    }

    for (j=0; j<ex_d->nRules; j++)
    {
      if (ex_d->tunable_rules[j])
        for (i=0; i<ex_d->rules[j]-1; i++)
        {
          ex_d->eta[j][i][0]=ex_d->eta[j][i][0]+
      ex_d->eta_temp[j][i][0]*ex_prob/rootProb;
          ex_d->eta[j][i][1]=ex_d->eta[j][i][1]+
      ex_d->eta_temp[j][i][1]*ex_prob/rootProb;
        }
    }
  }
  free(ex_d->sigma);
  return rootProb;
}


void Maximization(example_data * ex_d)
{
  int r,i,j,e;
  double sum=0;
  double *probs_rule,**eta_rule;

  for (r=0;r<ex_d->nRules;r++)
  {
    if (ex_d->tunable_rules[r])
    {
      eta_rule=ex_d->eta[r];
      for (i=0;i<ex_d->rules[r]-1;i++)
      {
        sum=(eta_rule[i][0]+eta_rule[i][1]);
        if (sum==0.0)
        {
          ex_d->arrayprob[r][i]=0;
        }
        else
          ex_d->arrayprob[r][i]=eta_rule[i][1]/sum;
      }
    }
  }

  for(e=0;e<ex_d->ex;e++)
  {
    for (j=0;j<ex_d->env[e].nVars;j++)
    {
      r=ex_d->env[e].vars[j].nRule;
      if (ex_d->tunable_rules[r])
      {
        probs_rule=ex_d->arrayprob[r];
        for(i=0;i<ex_d->rules[r]-1;i++)
        {
          ex_d->env[e].probs[ex_d->env[e].vars[j].firstBoolVar+i]=probs_rule[i];
        }
      }
    }
  }
}



static foreign_t init_par(example_data * ex_d, term_t ruleHeadsArg)
{
  double * theta,p0;
  double pmass,par;
  double **Theta_rules;
  double ***eta;
  double ***eta_temp;
  int ret,i,j,e,rule;
  term_t head=PL_new_term_ref();
  term_t p=PL_new_term_ref();
  term_t ruleHeadsTerm;
  size_t nHeads,nRules;
  int *rules, *tun_rules;

  ruleHeadsTerm=PL_copy_term_ref(ruleHeadsArg);
  ret=PL_skip_list(ruleHeadsTerm,0,&nRules);
  if (ret!=PL_LIST) return FALSE;

  ex_d->nRules=nRules;
  ex_d->rules= (int *) malloc(nRules * sizeof(int));
  rules=ex_d->rules;
  ex_d->tunable_rules= (int *) malloc(nRules * sizeof(int));
  tun_rules=ex_d->tunable_rules;
  ex_d->eta= (double ***) malloc(nRules * sizeof(double **));
  eta=ex_d->eta;
  ex_d->eta_temp= (double ***) malloc(nRules * sizeof(double **));
  eta_temp=ex_d->eta_temp;
  ex_d->nodes_probs=NULL;
  ex_d->arrayprob=(double **) malloc(nRules * sizeof(double *));

  Theta_rules=(double **)malloc(nRules *sizeof(double *));

  for (j=0;j<nRules;j++)
  {
    
    ret=PL_get_list(ruleHeadsTerm,head,ruleHeadsTerm);
    RETURN_IF_FAIL
    pmass=0;
    if (PL_is_list(head)) 
    { 
      ret=PL_skip_list(head,0,&nHeads);
      if (ret!=PL_LIST) return FALSE;
      if (nHeads==1) // fixed parameters
      {
        ret=PL_get_head(head,head);
        RETURN_IF_FAIL 
        ret=PL_skip_list(head,0,&nHeads);
        if (ret!=PL_LIST) return FALSE;
        tun_rules[j]=0;
      }
      else // initial parameters
        tun_rules[j]=1;
        
      Theta_rules[j]=(double *)malloc(nHeads*sizeof(double));
      theta=Theta_rules[j];
      rules[j]=nHeads,
      ex_d->arrayprob[j]= (double *) malloc((nHeads-1)*sizeof(double));
      eta[j]= (double **) malloc((nHeads-1)*sizeof(double *));
      eta_temp[j]= (double **) malloc((nHeads-1)*sizeof(double *));

      for (i=0;i<nHeads-1;i++)
      {
        ret=PL_get_list(head,p,head);
        RETURN_IF_FAIL
        ret=PL_get_float(p, &par);
        RETURN_IF_FAIL
        eta[j][i]=(double *) malloc(2*sizeof(double));
        eta_temp[j][i]=(double *) malloc(2*sizeof(double));
        pmass=pmass+par;
        theta[i]=par;
        ex_d->arrayprob[j][i]=par;
      }
      theta[nHeads-1]=1-pmass;
    }
    else
    {
      ret=PL_get_integer(head,&rules[j]);
      nHeads=rules[j];
      
      RETURN_IF_FAIL
      Theta_rules[j]=(double *)malloc(nHeads*sizeof(double));
      theta=Theta_rules[j];
      ex_d->arrayprob[j]= (double *) malloc((nHeads-1)*sizeof(double));
      
      eta[j]= (double **) malloc((nHeads-1)*sizeof(double *));
      eta_temp[j]= (double **) malloc((nHeads-1)*sizeof(double *));
      tun_rules[j]=1;
      for (i=0;i<rules[j]-1;i++)
      {
        eta[j][i]=(double *) malloc(2*sizeof(double));
        eta_temp[j][i]=(double *) malloc(2*sizeof(double));
      }
      if (ex_d->alpha==0.0)
      {
        for (i=0;i<rules[j]-1;i++)
        {
          par=uniform_sample()*(1-pmass);
          pmass=pmass+par;
          theta[i]=par;
          ex_d->arrayprob[j][i]=par;
        }
        theta[nHeads-1]=1-pmass;
      }
      else
      {
        symmetric_dirichlet_sample(ex_d->alpha,rules[j],theta);
        for (i=0;i<rules[j]-1;i++)
          ex_d->arrayprob[j][i]=theta[i];
      }
    }
  }

  for(e=0;e<ex_d->ex;e++)
  {
    for (j=0; j<ex_d->env[e].nVars; j++)
    {
      rule=ex_d->env[e].vars[j].nRule;
      theta=Theta_rules[rule];
      p0=1;
      for (i=0; i<ex_d->env[e].vars[j].nVal-1;i++)
      {
        ex_d->env[e].probs[ex_d->env[e].vars[j].firstBoolVar+i]=theta[i]/p0;
        p0=p0*(1-theta[i]/p0);
      }
    }
  }
  for (j=0;j<ex_d->nRules;j++)
  {
    free(Theta_rules[j]);
  }
  free(Theta_rules);
  return TRUE;
}
static foreign_t rand_seed(term_t arg1)
{
  int seed;
  int ret;

  ret=PL_get_integer(arg1,&seed);
  RETURN_IF_FAIL

  srand((unsigned)seed);

  PL_succeed;
}

static foreign_t EM(term_t arg1,term_t arg2,term_t arg3,term_t arg4,term_t arg5,
term_t arg6,term_t arg7,term_t arg8,term_t arg9)
{
  term_t pterm,nil,out1,out2,out3,nodesTerm,ruleTerm,tail,pair,compoundTerm;
  DdNode * node1,**nodes_ex;
  int r,i,iter,cycle,ret;
  long iter1;
  size_t lenNodes;
  example_data * ex_d;



  double CLL0= -2.2*pow(10,10); //-inf
  double CLL1= -1.7*pow(10,8);  //+inf
  double p,p0,**eta_rule,ea,er;
  double ratio,diff;

  term_t head=PL_new_term_ref();

  ret=PL_get_pointer(arg1,(void **)&ex_d);
  RETURN_IF_FAIL
  ret=init_par(ex_d,arg2);
  RETURN_IF_FAIL
  pair=PL_new_term_ref();
  nodesTerm=PL_copy_term_ref(arg3);

  ret=PL_skip_list(nodesTerm,0,&lenNodes);
  if (ret!=PL_LIST) return FALSE;

  out1=PL_new_term_ref();
  out2=PL_new_term_ref();
  out3=PL_new_term_ref();
  ruleTerm=PL_new_term_ref();
  tail=PL_new_term_ref();
  pterm=PL_new_term_ref();
  nil=PL_new_term_ref();
  compoundTerm=PL_new_term_ref();

  ret=PL_get_float(arg4,&ea);
  RETURN_IF_FAIL

  ret=PL_get_float(arg5,&er);
  RETURN_IF_FAIL
  ret=PL_get_integer(arg6,&iter);
  RETURN_IF_FAIL

  nodes_ex=(DdNode **)malloc(lenNodes*sizeof(DdNode*));
  ex_d->nodes_probs=(double *)malloc(lenNodes*sizeof(double));
  ex_d->example_prob=(double *)malloc(lenNodes*sizeof(double));

  for (i=0;i<lenNodes;i++)
  {
    ret=PL_get_list(nodesTerm,pair,nodesTerm);
    RETURN_IF_FAIL
    ret=PL_get_list(pair,head,pair);
    RETURN_IF_FAIL
    ret=PL_get_pointer(head,(void **)&node1);
    RETURN_IF_FAIL
    nodes_ex[i]=node1;
    ret=PL_get_list(pair,head,pair);
    RETURN_IF_FAIL
    ret=PL_get_float(head,&(ex_d->example_prob[i]));
    RETURN_IF_FAIL
  }
  diff=CLL1-CLL0;
  ratio=diff/fabs(CLL0);
  if (iter==-1)
    iter1= 2147000000;
  else iter1=iter;

  cycle=0;
  while  ( (diff>ea) && (ratio>er) && (cycle<iter1) )
  {
    cycle++;
    for (r=0;r<ex_d->nRules;r++)
    {
      if (ex_d->tunable_rules[r])
        for (i=0;i<ex_d->rules[r]-1;i++)
        {
          eta_rule=ex_d->eta[r];
          eta_rule[i][0]=0;
          eta_rule[i][1]=0;
        }
    }
    CLL0 = CLL1;
    CLL1 = Expectation(ex_d,nodes_ex,lenNodes);
    Maximization(ex_d);
    diff=CLL1-CLL0;
    ratio=diff/fabs(CLL0);
  }
  ret=PL_put_nil(out2);
  RETURN_IF_FAIL
  for (r=0; r<ex_d->nRules; r++)
  {
    ret=PL_put_nil(tail);
    RETURN_IF_FAIL
    p0=1;
    for (i=0;i<ex_d->rules[r]-1;i++)
    {
      p=ex_d->arrayprob[r][i]*p0;
      ret=PL_put_float(pterm,p);
      RETURN_IF_FAIL
      ret=PL_cons_list(tail,pterm,tail);
      RETURN_IF_FAIL
      p0=p0*(1-ex_d->arrayprob[r][i]);
    }
    ret=PL_put_float(pterm,p0);
    RETURN_IF_FAIL
    ret=PL_cons_list(tail,pterm,tail);
    RETURN_IF_FAIL
    ret=PL_put_integer(ruleTerm,r);
    RETURN_IF_FAIL
    ret=PL_put_nil(nil);
    RETURN_IF_FAIL
    ret=PL_cons_list(tail,tail,nil);
    RETURN_IF_FAIL
    ret=PL_cons_list(compoundTerm,ruleTerm,tail);
    RETURN_IF_FAIL
    ret=PL_cons_list(out2,compoundTerm,out2);
    RETURN_IF_FAIL
  }
  ret=PL_put_nil(out3);
  RETURN_IF_FAIL
  for (i=0;i<lenNodes;i++)
  {
    ret=PL_put_float(pterm,ex_d->nodes_probs[i]);
    RETURN_IF_FAIL
    ret=PL_cons_list(out3,pterm,out3);
    RETURN_IF_FAIL
  }
  ret=PL_unify(out3,arg9);
  RETURN_IF_FAIL

  ret=PL_put_float(out1,CLL1);
  RETURN_IF_FAIL
  ret=PL_unify(out1,arg7);
  RETURN_IF_FAIL
  free(nodes_ex);
  free(ex_d->example_prob);
  free(ex_d->nodes_probs);

  return (PL_unify(out2,arg8));
}


/*

static int paths_to_non_zero(void)
{
  double paths;
  YAP_Term arg1,arg2,out;
  DdNode * node;

  arg1=YAP_ARG1;
  arg2=YAP_ARG2;
  node=(DdNode *)YAP_IntOfTerm(arg1);
  paths=Cudd_CountPathsToNonZero(node);
  out=YAP_MkFloatTerm(paths);
  return(YAP_Unify(out,arg2));
}

static int paths(void)
{
  double paths;
  YAP_Term arg1,arg2,out;
  DdNode * node;

  arg1=YAP_ARG1;
  arg2=YAP_ARG2;
  node=(DdNode *)YAP_IntOfTerm(arg1);
  paths=Cudd_CountPath(node);
  out=YAP_MkFloatTerm(paths);
  return(YAP_Unify(out,arg2));
}

static int dag_size(void)
{
  int size;
  YAP_Term arg1,arg2,out;
  DdNode * node;

  arg1=YAP_ARG1;
  arg2=YAP_ARG2;
  node=(DdNode *)YAP_IntOfTerm(arg1);
  size=Cudd_DagSize(node);
  out=YAP_MkIntTerm(size);
  return(YAP_Unify(out,arg2));
}
*/

FILE * open_file(char *filename, const char *mode)
/* opens a file */
{
  FILE *fp;

  if ((fp = fopen(filename, mode)) == NULL)
  {
    perror(filename);
    exit(1);
  }
  return fp;
}


tablerow* init_table(int varcnt) {
  int i;
  tablerow *tab;

  tab = (tablerow *) malloc(sizeof(rowel) * varcnt);
  for (i = 0; i < varcnt; i++)
  {
    tab[i].row = NULL;
    tab[i].cnt = 0;
  }
  return tab;
}


void add_node(tablerow *tab, DdNode *node, double value) {
  int index = Cudd_NodeReadIndex(node);

  tab[index].row = (rowel *) realloc(tab[index].row,
    (tab[index].cnt + 1) * sizeof(rowel));
  tab[index].row[tab[index].cnt].key = node;
  tab[index].row[tab[index].cnt].value = value;
  tab[index].cnt += 1;
}

void add_or_replace_node(tablerow *tab, DdNode *node, double value)
{
  int i;
  int index = Cudd_NodeReadIndex(node);
  for(i = 0; i < tab[index].cnt; i++)
  {
    if (tab[index].row[i].key == node)
    {
      tab[index].row[i].value=value;
      return;
    }
  }
  tab[index].row = (rowel *) realloc(tab[index].row,
    (tab[index].cnt + 1) * sizeof(rowel));
  tab[index].row[tab[index].cnt].key = node;
  tab[index].row[tab[index].cnt].value = value;
  tab[index].cnt += 1;
}

double * get_value(tablerow *tab,  DdNode *node) {
  int i;
  int index = Cudd_NodeReadIndex(node);

  for(i = 0; i < tab[index].cnt; i++)
  {
    if (tab[index].row[i].key == node)
    {
      return &tab[index].row[i].value;
    }
  }
  return NULL;
}

void destroy_table(tablerow *tab,int varcnt)
{
  int i;

  for (i = 0; i < varcnt; i++)
  {
    free(tab[i].row);
  }
  free(tab);
}

expltablerow_abd *expl_init_table_abd(int varcnt) {
  int i;
  expltablerow_abd *tab;

  tab = malloc(sizeof(explrowel_abd) * varcnt);
  for (i = 0; i < varcnt; i++)
  {
    tab[i].row = NULL;
    tab[i].cnt = 0;
  }
  return tab;
}

expltablerow* expl_init_table(int varcnt) {
  int i;
  expltablerow *tab;

  tab = (expltablerow *) malloc(sizeof(explrowel) * varcnt);
  for (i = 0; i < varcnt; i++)
  {
    tab[i].row = NULL;
    tab[i].cnt = 0;
  }
  return tab;
}

void expl_add_node_abd(expltablerow_abd *tab, DdNode *node, int comp, prob_abd_expl_list value) {
  int index = Cudd_NodeReadIndex(node);

  tab[index].row = (explrowel_abd *) realloc(tab[index].row,(tab[index].cnt + 1) * sizeof(explrowel_abd));
  tab[index].row[tab[index].cnt].key.node = node;
  tab[index].row[tab[index].cnt].key.comp = comp;
  tab[index].row[tab[index].cnt].value = value;
  tab[index].cnt += 1;
}

void expl_add_node(expltablerow *tab, DdNode *node, int comp, prob_abd_expl value) {
  int index = Cudd_NodeReadIndex(node);

  tab[index].row = (explrowel *) realloc(tab[index].row,
    (tab[index].cnt + 1) * sizeof(explrowel));
  tab[index].row[tab[index].cnt].key.node = node;
  tab[index].row[tab[index].cnt].key.comp = comp;
  tab[index].row[tab[index].cnt].value = value;
  tab[index].cnt += 1;
}

// try to reuse
DdNode* get_node(DdNode *node, tablerow *tab) {
  int index = Cudd_NodeReadIndex(node);
  int i;

  for(i = 0; i < tab[index].cnt; i++) {
    if(tab[index].row[i].key == node) {
      return tab[index].row[i].key;
    }
  }
  return NULL;
}

// uguale a qeulla senza abd, cambia solo il tipo di dato
prob_abd_expl_list * expl_get_value_abd(expltablerow_abd *tab,  DdNode *node, int comp) {
  int i;
  int index = Cudd_NodeReadIndex(node);

  for(i = 0; i < tab[index].cnt; i++)
  {
    if (tab[index].row[i].key.node == node &&
       tab[index].row[i].key.comp == comp)
    {
      return &tab[index].row[i].value;
    }
  }
  return NULL;
}

prob_abd_expl * expl_get_value(expltablerow *tab,  DdNode *node, int comp) {
  int i;
  int index = Cudd_NodeReadIndex(node);

  for(i = 0; i < tab[index].cnt; i++)
  {
    if (tab[index].row[i].key.node == node &&
       tab[index].row[i].key.comp == comp)
    {
      return &tab[index].row[i].value;
    }
  }
  return NULL;
}

void expl_destroy_table_abd(expltablerow_abd *tab,int varcnt) {
  int i,j,k;

  for (i = 0; i < varcnt; i++) {
    for (j = 0; j< tab[i].cnt; j++) {
      for(k = 0; k < tab[i].row[j].value.n_elements; k++) {
        free_list(tab[i].row[j].value.mpa[k]);
      }
      free(tab[i].row[j].value.mpa);
    }
    free(tab[i].row);
  }
  free(tab);
}

void expl_destroy_table(expltablerow *tab,int varcnt)
{
  int i,j;

  for (i = 0; i < varcnt; i++)
  {
    for (j = 0; j< tab[i].cnt; j++)
    {
      free_list(tab[i].row[j].value.mpa);
    }
    free(tab[i].row);
  }
  free(tab);
}

// debugging functions

void print_abd_explan(explan_t *et) {
  explan_t *current_explan = et;
  while(current_explan != NULL) {
    printf("assign var: %d\nassign val: %d\n",current_explan->a.var,current_explan->a.val);
    current_explan = current_explan->next;
  }
}

void print_prob_abd_expl_list(const prob_abd_expl_list *pael) {
  int i;
  printf("--- print_prob_abd_expl_list --- \nprob: %lf\nn_explan: %d\nExplans:\n",pael->prob, pael->n_elements);
  if(pael->mpa == NULL) {
    printf("mpa NULL\n");
  }
  if(pael->mpa == NULL && pael->n_elements != 0) {
    printf("ERRORE: mpa null ed elementi != 0\n");
  }
  for(i = 0; i < pael->n_elements; i++) {
    printf("--- Explan %d ---\n",i);
    print_abd_explan(pael->mpa[i]); 
  }
  printf("-----\n");
}

void print_prob_abd_expl(prob_abd_expl *pae) {
  printf("prob: %lf\nExplan:\n",pae->prob);
  print_abd_explan(pae->mpa);
}

// prints all the fields of a variable data structure
void dump_var(variable *var) {
  printf("\t\tnVal: %d\n", var->nVal);
  printf("\t\tnRule: %d\n", var->nRule);
  printf("\t\tfirstBoolVar: %d\n", var->firstBoolVar);
  printf("\t\tabducible: %d\n", var->abducible);
  printf("\t\tquery: %d\n", var->query);
  printf("\t\tdecision: %d\n", var->decision);
}

// prints all the variables of the environment
void dump_env(environment *env) {
  int i;

  printf("n_abd_boolVars: %d\n",env->n_abd_boolVars);
  printf("n_abd: %d\n",env->n_abd);

  printf("nVars: %d\n",env->nVars);
  for(i = 0; i < env->nVars; i++) {
    printf("\tvars[%d]: \n",i);
    dump_var(&env->vars[i]);
  }

  printf("boolVars: %d\n",env->boolVars);
  for(i = 0; i < env->boolVars; i++) {
    printf("\tbVar2mVar[%d] = %d\n",i,env->bVar2mVar[i]);
    printf("\tprobs[%d] = %lf\n",i,env->probs[i]);
  }

  printf("nRules: %d\n",env->nRules);
  for(i = 0; i < env->nRules; i++) {
    printf("\trules[%d] = %d\n",i,env->bVar2mVar[i]);
  }
}

// prints CUDD venv variables
void debug_cudd_env(environment *env, int i) {
  printf("----- %d -----\n",i);
  printf("Dead nodes (Cudd_ReadDead): %d\n",Cudd_ReadDead(env->mgr));
  printf("Cudd check zero ref (node with non 0 ref, Cudd_CheckZeroRef): %d\n",Cudd_CheckZeroRef(env->mgr));
  printf("Cudd check keys (Cudd_CheckKeys): %d\n",Cudd_CheckKeys(env->mgr));
  printf("Cudd debug check (Cudd_DebugCheck): %d\n",Cudd_DebugCheck(env->mgr));
  printf("Cudd_ReadMaxMemory: %lu\n",Cudd_ReadMaxMemory(env->mgr));
  printf("DdManager vars: %d | ", Cudd_ReadSize(env->mgr)); /*Returns the number of BDD variables in existence*/
  printf("DdManager nodes: %ld | ", Cudd_ReadNodeCount(env->mgr)); /*Reports the number of live nodes in BDDs and ADDs*/
  printf("DdManager reorderings: %d | ", Cudd_ReadReorderings(env->mgr) ); /*Returns the number of times reordering has occurred*/
  printf("DdManager memory: %ld |\n\n", Cudd_ReadMemoryInUse(env->mgr) ); /*Returns the memory in use by the manager measured in bytes*/
  printf("-------------\n");
}

// predicate equivalent to debug_cudd_env
// out_null is empty, here just in case it's needed
static foreign_t debug_cudd_var(term_t env_ref, term_t out_null) {
  int ret;
  term_t out;
  environment *env;

  ret = PL_get_pointer(env_ref,(void **)&env);
  RETURN_IF_FAIL
  printf("----- DEBUG CUDD -----\n");
  printf("Dead nodes (Cudd_ReadDead): %d\n",Cudd_ReadDead(env->mgr));
  printf("Cudd check zero ref (node with non 0 ref, Cudd_CheckZeroRef): %d\n",Cudd_CheckZeroRef(env->mgr));
  printf("Cudd check keys (Cudd_CheckKeys): %d\n",Cudd_CheckKeys(env->mgr));
  printf("Cudd debug check (Cudd_DebugCheck): %d\n",Cudd_DebugCheck(env->mgr));
  printf("----- END DEBUG CUDD -----\n");
  out = PL_new_term_ref();
  ret = PL_put_nil(out);
  RETURN_IF_FAIL
  return(PL_unify(out,out_null));
}

install_t install()
/* function required by YAP for intitializing the predicates defined by a C function*/
{
  srand(10);

  PL_register_foreign("init_em",1,init_em,0);
  PL_register_foreign("init_ex",2,init_ex,0);
  PL_register_foreign("end_em",1,end_em,0);
  PL_register_foreign("end_ex",1,end_ex,0);
  PL_register_foreign("add_var",4,add_var,0);
  PL_register_foreign("add_query_var",4,add_query_var,0);
  PL_register_foreign("add_abd_var",4,add_abd_var,0);
  PL_register_foreign("equality",4,equality,0);
  PL_register_foreign("and",4,and,0);
  PL_register_foreign("one",2,one,0);
  PL_register_foreign("zero",2,zero,0);
  PL_register_foreign("or",4,or,0);
  PL_register_foreign("bdd_not",3,bdd_not,0);
  PL_register_foreign("create_dot",3,create_dot,0);
  PL_register_foreign("create_dot_string",3,create_dot_string,0);
  PL_register_foreign("init",1,init,0);
  PL_register_foreign("end",1,end,0);
  PL_register_foreign("ret_prob",3,ret_prob,0);
  PL_register_foreign("ret_abd_prob",4,ret_abd_prob,0);
  PL_register_foreign("ret_map_prob",4,ret_map_prob,0);
  PL_register_foreign("reorder",1,reorder,0);
  PL_register_foreign("make_query_var",3,make_query_var,0);
  PL_register_foreign("em",9,EM,0);
  PL_register_foreign("rand_seed",1,rand_seed,0);
  PL_register_foreign("gamma_sample",3,gamma_sample_pl,0);
  PL_register_foreign("gauss_sample",3,gauss_sample_pl,0);
  PL_register_foreign("uniform_sample",1,uniform_sample_pl,0);
  PL_register_foreign("dirichlet_sample",2,dirichlet_sample_pl,0);
  PL_register_foreign("symmetric_dirichlet_sample",3,symmetric_dirichlet_sample_pl,0);
  PL_register_foreign("discrete_sample",2,discrete_sample_pl,0);
  PL_register_foreign("initial_values",2,initial_values_pl,0);
  PL_register_foreign("add_decision_var",4,add_decision_var,0);
  PL_register_foreign("probability_dd",3,probability_dd,0);
  PL_register_foreign("add_prod",4,add_prod,0);
  PL_register_foreign("add_sum",4,add_sum,0);
  PL_register_foreign("ret_strategy",4,ret_strategy,0);
  PL_register_foreign("compute_optimal_strategy_sp",5,compute_optimal_strategy_sp,0);
  PL_register_foreign("debug_cudd_var",2,debug_cudd_var,0);
  // PL_register_foreign("add_const",3,add_const,0);

//  PL_register_foreign("deref",1,rec_deref,0);
//  PL_register_foreign("garbage_collect",2,garbage_collect,0);
//  PL_register_foreign("bdd_to_add",2,bdd_to_add,0);
//  PL_register_foreign("paths_to_non_zero",2,paths_to_non_zero,0);
//  PL_register_foreign("paths",2,paths,0);
//  PL_register_foreign("dag_size",2,dag_size,0);
}