/* 
vl = pointer from a node to an edge incident with it 
     this edge serves as the first element of the list 
     of edges incident with this node; edges are directed: ab = zero - ba;
ll = pointer to the next element of the list of edges incident with a node 
lv = node at edge's endpoint; 
q =  node degree

*/

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <string.h>

#define V 100000
#define L (V-1)
#define NOBS 6
#define MAXSEED 2000000000

/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/

int main(){

  void init_tree(int, int, int* , int* , int* , int* );
  void reshuffle_leaves(int, int, int*, int*, int*, int* );
  void reshuffle_branches(int, int, int*, int* ,int*, int* );
  void show_graph(int , int* , int* , int* , int*);
  void measure(int, int , int* , int* , int* , int*, double*,
               double*, double* ,double* ,double*);
  void read_in_configuration(int *, int *, long* , int* , int* , int* , int* );
  void save_configuration(int , int , int* , int* , int* , int* );
  void save_results(int, int, double*);
  void reset_cumulative_histograms(int,double*,double*,double*,double*);
  void save_cumulative_histograms(int,double*,double*,double*,double*);
  
  int vl[V];
  int ll[L+L];
  int lv[L+L];
  int q[V];
  double res[NOBS];
  double mh[V],mdd[V],bsd[V],lsd[V];

  int i;
  int nv,nl; 
  int init_option, v_option, mc_option;
  int irep,nrep,n0;
  int nobs=NOBS;
  long seed;
  int flag, sfreq;
  
  scanf("%d %d %d %d %ld %d %d",&init_option,&v_option,&mc_option,&nv,&seed,&nrep,&sfreq);
  srand48(seed);

  if(init_option==0) 
    read_in_configuration(&n0,&nv,&seed,vl,ll,lv,q);
  else {
    init_tree(init_option,nv,vl,ll,lv,q);
    n0=0;
    measure(n0,nv,vl,ll,lv,q,res,mh,mdd,bsd,lsd); 
    save_results(n0,nobs,res);                                       
  }
  
   reset_cumulative_histograms(nv,mh,mdd,bsd,lsd);
  
  //  show_graph(nv,vl,ll,lv,q); 
  
  for(irep=n0+1;irep<=nrep;irep++){  
     switch(v_option){
       case 0: 
         if(init_option<1 || init_option>3){
            printf("v_option=0 applies only if init_option=1/2/3\n");
            exit(5);
          }  
         init_tree(init_option,nv,vl,ll,lv,q);
       break;
       case 1: 
         reshuffle_leaves(mc_option,nv,vl,ll,lv,q);
       break;
       case 2: 
         reshuffle_branches(mc_option,nv,vl,ll,lv,q);
       break;
       default: printf("choose v_option 0-2 \n"); exit(4); 
     }

     if(irep % sfreq == 0) flag = irep; 
     else flag = -1;
     measure(flag,nv,vl,ll,lv,q,res,mh,mdd,bsd,lsd);
     save_results(irep,nobs,res);                                      
  }
  
  if(n0<nrep){                                   
    save_configuration(nrep,nv,vl,ll,lv,q);
    save_cumulative_histograms(nv,mh,mdd,bsd,lsd);
  }
//  show_graph(nv,vl,ll,lv,q); 

}

/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/
void init_tree(int init_option, int nv, int* vl, int* ll, int* lv, int* q){

// init_option 1=random tree (random Pruefer code)
//             2=exponential growing tree
//             3=Barabasi-Albert tree
//             4=line graph
//             5=star graph

  int nl=nv-1;
  int zero=nl+nl-1;
  int edge, next, va,vb,vc, link;
  int i,j,jp,k,flag;
  int aa,a[nv-2],b[nv+1],bb[nv+1];

  for(va=0;va<nv;va++) q[va]=0;
  
  switch(init_option){
  
  case 1: // Pruefer code
  
  for(i=0;i<nv-2;i++){
    aa=nv; while(aa==nv) aa=nv*drand48(); 
    a[i]=aa;
  }
  for(j=0;j<nv;j++){
    bb[j]=j+1; b[j+1]=j;
  }
  for(i=0;i<nv-2;i++){
    edge=i;
    j=0;
    do{
      jp=j;j=bb[jp]; 
      flag=0;
      for(k=i;k<nv-2;k++) if(a[k]==b[j]) flag=1;
    } while(flag==1); 
      
    lv[edge]=b[j];
    lv[zero-edge]=a[i];
    bb[jp]=bb[j];
  }
  edge++;
  lv[edge]=b[bb[0]];
  lv[zero-edge]=b[bb[bb[0]]];
  
  for(i=0;i<nv;i++) vl[i]=-1;
  for(i=0;i<nv;i++) q[i]=0;
  for(edge=0;edge<nl+nl;edge++){
    va=lv[zero-edge];
    if(vl[va]>=0) ll[edge]=vl[va];
    vl[va]=edge;
    q[va]++;
  }

  break;
  
  case 2: // exponential growing tree 
  
    for(edge=0;edge<nl;edge++){
      vb=edge+1;
      va=vb; while(va==vb) va=(int) vb*drand48();
       
      next=vl[va];
      vl[va]=edge;
      ll[edge]=next;  
      lv[edge]=vb;
       
      next=vl[vb];
      vl[vb]=zero-edge;
      ll[zero-edge]=next;
      lv[zero-edge]=va;
       
      q[va]++;
      q[vb]++;
   }
   
   break;

   case 3: // BA tree
  
   lv[0]=1; lv[zero]=0;
   q[0]=q[1]=1;
   vl[0]=0; vl[1]=zero;
  
     for(edge=1;edge<nl;edge++){
       vb=edge+1;
       link=edge; while(link==edge) link=(int) edge*drand48(); 
       if(drand48()<0.5) va=lv[link]; else va=lv[zero-link];
       
       next=vl[va];
       vl[va]=edge;
       ll[edge]=next;  
       lv[edge]=vb;
       
       next=vl[vb];
       vl[vb]=zero-edge;
       ll[zero-edge]=next;
       lv[zero-edge]=va;
       
       q[va]++;
       q[vb]++;
    }
    break;
  
    case 4: //line graph 
      for(edge=0;edge<nl;edge++){
        va=edge; vb=edge+1;
        lv[edge]=vb; lv[zero-edge]=va;
        vl[va]=edge;
        ll[edge]=zero-(edge-1);
        q[vb]=2;
      }  
      q[0]=1;
      q[nv-1]=1;
      vl[nv-1]=zero-(nl-1);
    break;
    
    case 5: //star graph
      for(edge=0;edge<nl;edge++){
        vb=edge+1;
        lv[edge]=vb; lv[zero-edge]=0;
        vl[vb]=zero-edge;
        ll[edge]=(edge+1);
        q[vb]=1;
      }  
      q[0]=nl;
      vl[0]=0;
    break;
  
    default: printf("choose init-option 1-5 \n"); exit(3);  
  }  
}

/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/

void reshuffle_leaves(int mc_option, int nv, int* vl, int* ll, int* lv, int* q){

/* nv = number of vertices 
   nl = number of edges
   zero = reflection point for directed edges 
*/ 

  int nl=nv-1;
  int zero=nl+nl-1;
  int edge, next, va,vb,vc;
  int list_element,leaves;
  double weight;


/* cut and paste nv leaves */

  leaves=0;
  while(leaves<nv){
    va=nv;
    while(va==nv) va = drand48()*nv;
    if(q[va]==1){
      leaves++;
      edge=vl[va];
      vb=lv[edge];
      vc=va; while(vc==va||vc==nv) vc=nv*drand48();
      
      switch(mc_option){
      case 1: weight=1.0; 
                                    break; // Cayley trees 1/e 1/(q-1)!
      case 2: weight=q[vc]/(q[vb]-1.0); 
                                    break; // exponential trees 2^{-q} 
      case 3: weight=
      (q[vb]+2.0)*q[vc]*q[vc]/(q[vc]+3.0)/(q[vb]-1.0)/(q[vb]-1.0);
                                    break; // BA-trees 4/(q(q+1)(q+2)
      default: printf("mc_option must be 1,2, or 3\n"); exit(1);
      }
      
      if(drand48()< weight){
        if(vl[vb]==(zero-edge)) vl[vb]=ll[zero-edge];
        else
          {
           list_element=vl[vb];
           while(ll[list_element]!=(zero-edge)) list_element=ll[list_element]; 
           ll[list_element]=ll[zero-edge];
          }
        lv[edge]=vc;
        next=vl[vc];
        vl[vc]=zero-edge;
        ll[zero-edge]=next;
        q[vc]++;
        q[vb]--;
       }
    }
  }
}

/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/

void reshuffle_branches(int mc_option, int nv, int* vl, int* ll, int* lv, int* q){

/* nv = number of vertices 
   nl = number of edges
   zero = reflection point for directed edges 
*/ 

  int nl=nv-1;
  int zero=nl+nl-1;
  int edge, next, va,vb,vc;
  int list_element;
  int i,ih;
  int black, white;
  int color[nv];
  int white_heap[nv], black_heap[nv], nwh, nbh;
  int v1,v2,v3,i1,i2;
  double weight;


/* cut and paste nv branches */

  black=0;
  for(i=0;i<nv;i++) color[i]=black;

  for(i=0;i<nv;i++){
    white=black+1;
    black=white+1;
    edge=nl+nl;
    while(edge==(nl+nl)) edge = drand48()*(nl+nl);
    vb=lv[edge];
    va=lv[zero-edge];
    if(q[vb]>1){
      vc=va; while(vc==va||vc==nv) vc=nv*drand48();
      
      switch(mc_option){
        case 1: weight=1.0; 
                                    break; // Cayley trees 1/e 1/(q-1)!
        case 2: weight=q[vc]/(q[vb]-1.0); 
                                    break; // exponential trees 2^{-q} 
        case 3: weight=
              (q[vb]+2.0)*q[vc]*q[vc]/(q[vc]+3.0)/(q[vb]-1.0)/(q[vb]-1.0);
                                    break; // BA-trees 4/(q(q+1)(q+2)
        default: printf("mc_option must be 1,2, or 3\n"); exit(1);
      }
      
      if(drand48()< weight){
        color[va]=white;
        color[vb]=black;
        
        nwh=1; white_heap[0]=va;
        nbh=1; black_heap[0]=vb;
         
        ih=0;  
        while(ih<nwh && ih<nbh){
          v1=white_heap[ih]; 
          list_element=vl[v1];
          for(i1=0;i1<q[v1];i1++){
            v3=lv[list_element];
            if(color[v3]<white){color[v3]=white; white_heap[nwh]=v3; nwh++;}
            list_element=ll[list_element];
          }
          v2=black_heap[ih];
          list_element=vl[v2];
          for(i2=0;i2<q[v2];i2++){
            v3=lv[list_element];
            if(color[v3]<white){color[v3]=black;  black_heap[nbh]=v3; nbh++;}
            list_element=ll[list_element];
          }
          ih++;
        }          
        if((ih==nwh && color[vc]!=white) || color[vc]==black){
/* cut + paste */     
          if(vl[vb]==(zero-edge)) vl[vb]=ll[zero-edge];
          else{
            list_element=vl[vb]; 
            while(ll[list_element]!=(zero-edge)) list_element=ll[list_element]; 
            ll[list_element]=ll[zero-edge];
          }
          lv[edge]=vc;   
          next=vl[vc];
          vl[vc]=zero-edge;
          ll[zero-edge]=next;
          q[vc]++;
          q[vb]--;
        }
      }
    }
  }
}

/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/

void show_graph(int nv, int* vl, int* ll, int* lv, int* q){

  int i,j;
  int list_element;

  for(i=0;i<nv;i++){
    printf("%d : %d : ",i, q[i]);
    list_element=vl[i];
    for(j=0;j<q[i];j++){ 
      printf("%d ", lv[list_element]);
      list_element=ll[list_element];
    }
    printf("\n");   
  }
}

/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/
void measure(int flag, int nv, int* vl, int* ll, int* lv, int* q, double *res,
      double *mh, double *mdd, double *bsd, double *lsd){

  void distance_distribution(int, int , int* , int* , int* , int*, 
                             double*, double*);
  void tomography(int, int, int* , int* , int* , int* , 
                  double*, double*);
  void branch_size_distr(int, int, int* , int* , int* , int* , 
                         double*, double*);

  void int2str(int , char *);
  
  int i,qmax,next_to_leaves,va;
  int marked[nv];
  int h[nv];
  double q2,sigma;
  char suffix[80],filename[80];
  FILE *file;

  for(i=0;i<nv;i++) h[i]=0;
  for(i=0;i<nv;i++) h[q[i]]++;
  for(i=0;i<nv;i++) mh[i]+=h[i];
  
  for(i=0;i<80;i++) suffix[i]=0;
  
  if(flag>=0){
    int2str(flag,suffix);
    strcpy(filename,"./NDD/hist");
    strcat(filename,suffix);
    strcat(filename,".dat");
    file=fopen(filename,"w+");
    for(i=0;i<nv;i++) if(h[i]>0) fprintf(file,"%d, %lf\n", i, h[i]/((double) nv));   
    fclose(file);
  }

// width of the degree distribution
  q2=0.0;
  for(i=0;i<nv;i++) q2+=q[i]*q[i];
  q2/=nv;
  
  sigma = sqrt(q2 - (2.-2./nv)*(2.-2./nv));

// range of the degree distribution
  qmax=0;
  for(i=0;i<nv;i++) if(q[i]>qmax) qmax=q[i];
  
// ratio of the number of neighbours of leaves and of leaves 

  for(i=0;i<nv;i++) marked[i]=0;
  next_to_leaves=0;
  
  for(i=0;i<nv;i++){
    if(q[i]==1){
      va=lv[vl[i]];
      if(marked[va]==0){
        marked[va]=1; next_to_leaves++;
      }
    }
  }   

// collect results 
 
  res[0]=sigma;
  res[1]=next_to_leaves/(double) h[1];

  distance_distribution(flag,nv,vl,ll,lv,q,res,mdd);
  tomography(flag,nv,vl,ll,lv,q,res,lsd);
  branch_size_distr(flag,nv,vl,ll,lv,q,res,bsd);
  
// res[2]= mean distance;  res[3]= std(distance); 
// res[4]=number of layers; res[5]= stem length;
 
}
/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/

void 
distance_distribution(int flag, int nv, int* vl, int* ll, int* lv, int* q, 
                      double *res, double *mdd){
  
  void int2str(int , char *);
  char suffix[80],filename[80];
  FILE *file;
  
  int i,j;
  int list_element;
  int va,vb;
  int generation[nv];
  int dd[nv];
  int heap[nv];
  int ih,nh;
  double m1,m2,sigma;
  int max;


  for(i=0;i<nv;i++) dd[i]=0;
  
  for(i=0;i<nv;i++){
    for(j=0;j<nv;j++) generation[j]=-1;
    generation[i]=0; nh=0; heap[nh]=i; 
    for(ih=0;ih<=nh;ih++){
      va=heap[ih];
      list_element=vl[va];
      for(j=0;j<q[va];j++){
        vb=lv[list_element];
        if(generation[vb]==-1){
          generation[vb]=generation[va]+1; nh++; heap[nh]=vb;
        }
        list_element=ll[list_element];
      } 
    }
  for(j=0;j<nv;j++) dd[generation[j]]++;
  }
  
  for(j=0;j<nv;j++) mdd[j]+=dd[j];

  m1=m2=0.0;
  for(j=1;j<nv;j++){ m1 += j*(double) dd[j]; m2 += j*j*(double) dd[j];}
    m1/=(nv*(nv-1.0));
    m2/=(nv*(nv-1.0));
    sigma = sqrt(m2-m1*m1);
    max=nv-1; while(dd[max]==0) max--;

  res[2]=m1;
  res[3]=max;

  for(i=0;i<80;i++) suffix[i]=0;  
  if(flag>=0){
    int2str(flag,suffix);
    strcpy(filename,"./DD/hist");
    strcat(filename,suffix);
    strcat(filename,".dat");
    file=fopen(filename,"w+");
    for(i=0;i<nv;i++) if(dd[i]>0) fprintf(file,"%d, %lf\n", i, dd[i]/(nv-1.0)/nv);   
    fclose(file);
  }
  
}
  
/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/

// nvc = number of vertices in the core (main stem)
// n = number of branches in the pruned tree 


void 
tomography(int flag, int nv, int* vl, int* ll, int* lv, int* q, 
           double *res, double *lsd){  

  void int2str(int , char *);
  char suffix[80],filename[80];
  FILE *file;

  int nl=nv-1;
  int zero=nl+nl-1;
  int vlc[nv];      /* c stands for cloned versions */
  int llc[nl+nl];
  int lvc[nl+nl];
  int qc[nv];
  int leaf[nv], next_to_leaf[nv];
  int i,j,k,n,nvc,ngen;
  int va,vb,edge,list_element;
  int layer_size[nv];

  for(i=0;i<nv;i++) vlc[i]=vl[i];
  for(i=0;i<nl+nl;i++) lvc[i]=lv[i];
  for(i=0;i<nl+nl;i++) llc[i]=ll[i];
  for(i=0;i<nv;i++) qc[i]=q[i];

  j=0;
  for(i=0;i<nv;i++){
    if(qc[i]==1){
      leaf[j]=i;
      qc[i]=0;
      j++;
    }
  }
  n=j;
  
  nvc=nv;
  ngen=0;
  while(n>2){
    for(i=0;i<n;i++){
      va=leaf[i];
      edge=vlc[va];
      vb=lvc[edge];
      next_to_leaf[i]=vb;
      if(vlc[vb]==(zero-edge)) vlc[vb]=llc[zero-edge];
      else{
        list_element=vlc[vb];
        while(llc[list_element]!=(zero-edge)) list_element=llc[list_element]; 
        llc[list_element]=llc[zero-edge];
      }
      qc[vb]--;
    }
    nvc-=n;
    layer_size[ngen]=n;
    ngen++;
    
    j=0;
    for(i=0;i<n;i++){
      k=next_to_leaf[i];
      if(qc[k]==1){
        leaf[j]=k; 
        qc[k]=0;
        j++;
      }
    }
    n=j;
  }
  
  for(i=0;i<ngen;i++) lsd[i]+=layer_size[i];
//  lsd[ngen]+=nvc;
  
  res[4]=ngen;
  res[5]=nvc;
 
  for(i=0;i<80;i++) suffix[i]=0;  
  if(flag>=0){
    int2str(flag,suffix);
    strcpy(filename,"./TOM/hist");
    strcat(filename,suffix);
    strcat(filename,".dat");
    file=fopen(filename,"w+");
    for(i=0;i<ngen;i++) fprintf(file,"%d %d\n", i, layer_size[i]);
    fprintf(file,"%d %d\n", ngen,nvc);
    fclose(file);
  }
}
  
/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/

// nvc = number of vertices in the core (main stem)
// n = number of branches in the pruned tree 


void 
branch_size_distr(int flag, int nv, int* vl, int* ll, int* lv, int* q, 
                  double *res, double *bsd){  

  int nl=nv-1;
  int zero=nl+nl-1;
  int vlc[nv];      /* c stands for cloned versions */
  int llc[nl+nl];
  int lvc[nl+nl];
  int qc[nv];
  int leaf[nv], next_to_leaf[nv];
  int i,j,k,n,nvc,ngen;
  int va,vb,edge,list_element;
  int layer_size[nv];
  int bsl[nl+nl]; /* branch size for links */
  int bsv[nv]; /* auxiliary variable branch size for vertices */
  int x,hbs[nv/2+1];

  for(i=0;i<nv;i++) vlc[i]=vl[i];
  for(i=0;i<nl+nl;i++) lvc[i]=lv[i];
  for(i=0;i<nl+nl;i++) llc[i]=ll[i];
  for(i=0;i<nv;i++) qc[i]=q[i];
  for(i=0;i<nv;i++) bsv[i]=0;
  for(i=0;i<nl+nl;i++) bsl[i]=0;

  j=0;
  for(i=0;i<nv;i++){
    if(qc[i]==1){
      leaf[j]=i;
      qc[i]=0;
      j++;
    }
  }
  
  n=j;
  
  nvc=nv;
  ngen=0;
  while(n>2){
    for(i=0;i<n;i++){
      va=leaf[i];
      edge=vlc[va];
      vb=lvc[edge];
      next_to_leaf[i]=vb;
      bsl[edge]=bsv[va]+1;
      bsv[vb]+=bsl[edge];
      if(vlc[vb]==(zero-edge)) vlc[vb]=llc[zero-edge];
      else{
        list_element=vlc[vb];
        while(llc[list_element]!=(zero-edge)) list_element=llc[list_element]; 
        llc[list_element]=llc[zero-edge];
      }
      qc[vb]--;
    }
    nvc-=n;
    layer_size[ngen]=n;
    ngen++;
    
    j=0;
    for(i=0;i<n;i++){
      k=next_to_leaf[i];
      if(qc[k]==1){
        leaf[j]=k; 
        qc[k]=0;
        j++;
      }
    }
    n=j;
  }
  
  va=leaf[0];
  for(i=0;i<nvc-1;i++){  
    edge=vlc[va];
    vb=lvc[edge];
    bsl[edge]=bsv[va]+1;
    bsv[vb]+=bsl[edge];
    if(vlc[vb]==(zero-edge)) vlc[vb]=llc[zero-edge];
    else{
      list_element=vlc[vb];
      while(llc[list_element]!=(zero-edge)) list_element=llc[list_element]; 
      llc[list_element]=llc[zero-edge];
    }
    va=vb;
  }
    
  for(i=0;i<nl+nl;i++)
    if(bsl[i]==0) bsl[i]=nv-bsl[zero-i];

  for(i=0;i<nv/2+1;i++){
    hbs[i]=0;
  }
  
  for(i=0;i<nl;i++){
    x=bsl[i]; if(x>nv/2) x=nv-x;
    hbs[x]++;
  }
  
  for(i=0;i<nv/2+1;i++){
    bsd[i]+=hbs[i];
//  if(dd(hbs[i]>0) printf("%d %d\n",i,hbs[i]);
  }  

/*
  for(i=0;i<nl;i++)
   if(bsl[i]<bsl[zero-i])
     printf("%d %d \n",i, bsl[i]);
   else
     printf("%d %d \n",i, bsl[zero-i]);
  
  for(i=0;i<nv;i++){
    printf("%d :  ",i);
    edge=vl[i];
    printf("%d ",lv[edge]);
    for(j=1;j<q[i];j++){
      edge=ll[edge];
       printf("%d ",lv[edge]);
    }
    printf("\n");
  }
*/  
}
  
/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/  

void int2str(int i, char *s){

  int dec,n,j,k;
  dec=1;n=1;
  while(10*dec<=i){dec*=10;n++;}
  for(j=0;j<n;j++){
    k=i/dec;
    s[j]=48 + (char) k;
    i=i-dec*k;
    dec=dec/10;
  }
}

/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/

void save_results(int n, int nobs, double *res){                                      
     int i;
     printf("%d ",n);
     for(i=0;i<nobs;i++) printf("%lf ",res[i]);
     printf("\n");
}

/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/

void save_configuration(int n, int nv, int* vl, int* ll, int* lv, int* q){

  FILE *file;
  int i;
  int nl;
  long seed;
  
  nl=nv-1;
  seed=MAXSEED*drand48();
  file = fopen("config.dat","w");
  fprintf(file,"%d ", n);
  fprintf(file,"%d ",nv);
  fprintf(file,"%ld ",seed);
  for(i=0;i<nv;i++) fprintf(file,"%d ", vl[i]);
  for(i=0;i<nl+nl;i++) fprintf(file,"%d ", ll[i]);
  for(i=0;i<nl+nl;i++) fprintf(file,"%d ", lv[i]); 
  for(i=0;i<nv;i++) fprintf(file,"%d ", q[i]); 
  fclose(file);
  
}

/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/

void 
read_in_configuration(int *n, int *nv, long *seed, int* vl, int* ll, int* lv, int* q){

  FILE *file;
  int i;
  int nl;
  int nv_old;
  
  system("cp config.dat old_config.dat");  
  file = fopen("config.dat","r");
  fscanf(file,"%d", n);
  fscanf(file,"%d", &nv_old);
  if(nv_old!=(*nv)){printf("wrong size in config\n"); fclose(file); exit(2);} 
    fscanf(file,"%ld",seed); 
  nl=*nv-1;
  for(i=0;i<*nv;i++) fscanf(file,"%d", &vl[i]); 
  for(i=0;i<nl+nl;i++) fscanf(file,"%d", &ll[i]); 
  for(i=0;i<nl+nl;i++) fscanf(file,"%d", &lv[i]); 
  for(i=0;i<*nv;i++) fscanf(file,"%d", &q[i]); 
  fclose(file);
  
}

/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/

void reset_cumulative_histograms(int nv, 
     double *mh, double *mdd,double *bsd, double *lsd){
 
  for(int i=0;i<nv;i++) mh[i]=mdd[i]=bsd[i]=lsd[i]=0.0;

}

/*------------------------------------------------------------------------*/
/*------------------------------------------------------------------------*/

void save_cumulative_histograms(int nv,  
        double *mh, double *mdd,double *bsd, double *lsd){

  FILE *file;
  int i;
  double norm_mh,norm_mdd,norm_bsd,norm_lsd;
 
  norm_mh=norm_mdd=norm_bsd=norm_lsd=0.0;
  
  for(i=0;i<nv;i++){
    norm_mh+=mh[i];
    norm_mdd+=mdd[i];
    norm_bsd+=bsd[i];
    norm_lsd+=lsd[i];
  }
  
  for(i=0;i<nv;i++){
    mh[i]=mh[i]/norm_mh;
    mdd[i]=mdd[i]/norm_mdd;
    bsd[i]=bsd[i]/norm_bsd;
    lsd[i]=lsd[i]/norm_lsd;
  }
  
  file=fopen("DD/cum_hist.dat","w");
  for(i=0;i<nv;i++){
    if(mdd[i]> 0.5/norm_mdd) fprintf(file,"%d %lf\n",i,mdd[i]);
  }
  fclose(file);
  
  file=fopen("NDD/cum_hist.dat","w");
  for(i=0;i<nv;i++){
    if(mh[i]> 0.5/norm_mh) fprintf(file,"%d %lf\n",i,mh[i]);
  }
  fclose(file);
  
  file=fopen("TOM/cum_hist.dat","w");
  for(i=0;i<nv;i++){
    if(lsd[i]> 0.5/norm_lsd) fprintf(file,"%d %lf\n",i,lsd[i]);
  }
  fclose(file);
  
  file=fopen("BSD/cum_hist.dat","w");
  for(i=0;i<nv;i++){
    if(bsd[i]> 0.5/norm_bsd) fprintf(file,"%d %lf\n",i,bsd[i]);
  }
  fclose(file);
  
}