binaries/src/ViennaRNA/Kinfold/nachbar.c

   1 /*
   2   Last changed Time-stamp: <2011-03-10 18:02:10 ivo>
   3   c  Christoph Flamm and Ivo L Hofacker
   4   {xtof,ivo}@tbi.univie.ac.at
   5   Kinfold: $Name:  $
   6   $Id: nachbar.c,v 1.8 2008/06/03 21:55:11 ivo Exp $
   7 */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <string.h>
  12 #include <math.h>
  13 #include "globals.h"
  14 #include "assert.h"
  15 #include "fold_vars.h"
  16 #include "energy_const.h"
  17 #include "utils.h"
  18 #include "cache_util.h"
  19 #include "baum.h"
  20
  21 static char UNUSED rcsid[]="$Id: nachbar.c,v 1.8 2008/06/03 21:55:11 ivo Exp $";
  22
  23 /* arrays */
  24 static short *neighbor_list=NULL;
  25 static float *bmf=NULL; /* boltzmann weight of structure */
  26 static const char *costring(const char *str);
  27
  28 /* globals for laplace stuff */
  29 static double L = 0.0;
  30 static double D = 0.0;
  31 static double sumT = 0.0;
  32 static double sumK = 0.0;
  33 static double sumKK = 0.0;
  34 static double sumD = 0.0;
  35 static double *energies=NULL; /* energies of neighbors */
  36
  37 /* variables */
  38 /*  static double highestE = -1000.0; */
  39 /*  static double OhighestE = -1000.0; */
  40 /*  static char *highestS, *OhighestS; */
  41 static int lmin = 1;
  42 static int top = 0;
  43 static int is_from_cache = 0;
  44 /*  static double meanE = 0.0; */
  45 static double totalflux = 0.0;
  46 static double Zeit = 0.0;
  47 static double zeitInc = 0.0;
  48 static double _RT = 0.6;
  49
  50 /* public functiones */
  51 void ini_nbList(int chords);
  52 void update_nbList(int i, int j, int iE);
  53 int sel_nb(void);
  54 void clean_up_nbList(void);
  55 extern void update_tree(int i, int j);
  56
  57 /* privat functiones */
  58 static void reset_nbList(void);
  59 static void grow_chain(void);
  60 static FILE *logFP=NULL;
  61
  62 /**/
  63 void ini_nbList(int chords) {
  64   char logFN[256];
  65
  66   _RT = (((temperature + K0) * GASCONST) / 1000.0);
  67   if (neighbor_list!=NULL) return;
  68   /*
  69     list for move coding
  70     make room for 2*chords neighbors (safe bet)
  71   */
  72   if (chords == 0) chords = 1;
  73   neighbor_list = (short *)calloc(4*chords, sizeof(short));
  74   assert(neighbor_list != NULL);
  75   /*
  76     list for Boltzmann-factors
  77   */
  78   bmf = (float *)calloc(2*chords, sizeof(double));
  79   assert(bmf != NULL);
  80
  81   /* list of neighbor energies */
  82   energies = (double*)calloc(2*chords, sizeof(double));
  83   assert(energies != NULL);
  84
  85   /* open log-file */
  86   logFP = fopen(strcat(strcpy(logFN, GAV.BaseName), ".log"), "a+");
  87   assert(logFP != NULL);
  88
  89   /* log initial condition */
  90   log_prog_params(logFP);
  91   log_start_stop(logFP);
  92 }
  93
  94 /**/
  95 void update_nbList(int i, int j, int iE) {
  96   double E, dE, p;
  97
  98   E = (double)iE/100.;
  99   neighbor_list[2*top] = (short )i;
 100   neighbor_list[2*top+1] = (short )j;
 101
 102   /* compute rates and some statistics */
 103   /*    meanE += E; */
 104   dE = E-GSV.currE;
 105
 106   /* laplace stuff */
 107   energies[top] = E;
 108   L += GSV.currE-E;
 109   D++;
 110   /* fprintf(stderr, ">>%g %g<<\n", L, D); */
 111
 112   if( GTV.mc ) {
 113     /* metropolis rule */
 114     if (dE < 0) p = 1;
 115     else p = exp(-(dE / _RT*GSV.phi));
 116   }
 117   else  /* kawasaki rule */
 118     p = exp(-0.5 * (dE / _RT*GSV.phi));
 119
 120   totalflux += p;
 121   bmf[top++] = (float )p;
 122   if (dE < 0) lmin = 0;
 123   if ((dE == 0) && (lmin==1)) lmin = 2;
 124 }
 125
 126 /**/
 127 void get_from_cache(cache_entry *c) {
 128   top = c->top;
 129   totalflux = c->flux;
 130   GSV.currE = c->energy;
 131   lmin = c->lmin;
 132   memcpy(neighbor_list, c->neighbors, 2*top*sizeof(short));
 133   memcpy(bmf, c->rates, top*sizeof(float));
 134   memcpy(energies, c->energies, top*sizeof(double));
 135   is_from_cache = 1;
 136 }
 137
 138 /**/
 139 void put_in_cache(void) {
 140   cache_entry *c;
 141
 142   if ((c = (cache_entry *) malloc(sizeof(cache_entry)))==NULL) {
 143     fprintf(stderr, "out of memory\n"); exit(255);
 144   }
 145 /*    c->structure = strdup(GAV.currform); */
 146   c->structure = (char *) calloc(GSV.len+1, sizeof(char));
 147   strcpy(c->structure, GAV.currform);
 148   c->neighbors = (short *) malloc(top*2*sizeof(short));
 149   memcpy(c->neighbors,neighbor_list,top*2*sizeof(short));
 150   c->rates = (float *) malloc(top*sizeof(float));
 151   memcpy(c->rates, bmf, top*sizeof(float));
 152   c->energies = (double*)malloc(top*sizeof(double));
 153   memcpy(c->energies, energies, top*sizeof(double));
 154   c->top = top;
 155   c->lmin = lmin;
 156   c->flux = totalflux;
 157   c->energy = GSV.currE;
 158   write_cache(c);
 159 }
 160
 161 /*============*/
 162
 163 int sel_nb(void) {
 164
 165   char trans, **s;
 166   int next, i;
 167   double pegel = 0.0, schwelle = 0.0, zufall = 0.0;
 168   int found_stop=0;
 169
 170   /* before we select a move, store current conformation in cache */
 171   /* ... unless it just came from there */
 172   if ( !is_from_cache ) put_in_cache();
 173   else
 174     /* laplace stuff */
 175     for (i=0; i<top; i++) {
 176       L += (GSV.currE - energies[i]);
 177       D++;
 178     }
 179   is_from_cache = 0;
 180
 181   /* draw 2 different a random number */
 182   schwelle = urn();
 183   while ( zufall==0 ) zufall = urn();
 184
 185   /* advance internal clock */
 186   if (totalflux>0)
 187     zeitInc = (log(1. / zufall) / totalflux);
 188   else {
 189     if (GSV.grow>0) zeitInc=GSV.grow;
 190     else zeitInc = GSV.time;
 191   }
 192
 193   Zeit += zeitInc;
 194
 195   /* laplace stuff */
 196   sumK  += L*zeitInc;
 197   sumKK += L*L*zeitInc;
 198   sumD  += D*zeitInc;
 199
 200   if (GSV.grow>0 && GSV.len < strlen(GAV.farbe_full)) grow_chain();
 201
 202   /* meanE /= (double)top; */
 203
 204   /* normalize boltzmann weights */
 205   schwelle *=totalflux;
 206
 207   /* and choose a neighbour structure next */
 208   for (next = 0; next < top; next++) {
 209     pegel += bmf[next];
 210     if (pegel > schwelle) break;
 211   }
 212
 213   /* in case of rounding errors */
 214   if (next==top) next=top-1;
 215
 216   /*
 217     process termination contitiones
 218   */
 219   /* is current structure identical to a stop structure ?*/
 220   for (found_stop = 0, s = GAV.stopform; *s; s++) {
 221     if (strcmp(*s, GAV.currform) == 0) {
 222       found_stop = (s - GAV.stopform) + 1;
 223       break;
 224     }
 225   }
 226
 227   if ( ((found_stop > 0) && (GTV.fpt == 1)) || (Zeit > GSV.time) ) {
 228     /* met condition to stop simulation */
 229
 230     /* laplace stuff */
 231     double K, KK, N, sigma;
 232     K = sumK/Zeit;
 233     KK = sumKK/Zeit;
 234     N = sumD/Zeit;
 235     /* graph Laplacian is - Laplace-Beltrami operator */
 236     sigma = -1.0*sqrt((KK-K*K)/N)/(K/N);
 237
 238     /* this goes to stdout */
 239     if ( !GTV.silent ) {
 240       printf("%s %6.2f %10.3f", costring(GAV.currform), GSV.currE, Zeit);
 241
 242       /* laplace stuff*/
 243       if (GTV.phi) printf(" %8.3f %8.3f %3g", zeitInc, L, D);
 244
 245       if (GTV.verbose) printf(" %4d _ %d", top, lmin);
 246       if (found_stop) printf(" X%d\n", found_stop);/* found a stop structure */
 247       else printf(" O\n"); /* time for simulation is exceeded */
 248
 249       /* laplace stuff */
 250       if (GTV.phi) printf("Curvature fluctuation sigma = %7.5f\n", sigma);
 251
 252       fflush(stdout);
 253     }
 254
 255     /* this goes to log */
 256     /* comment log steps of simulation as well !!! %6.2f  round */
 257     if ( found_stop ) {
 258       fprintf(logFP," X%02d %12.3f", found_stop, Zeit);
 259
 260       /* laplace stuff */
 261       if (GTV.phi) fprintf(logFP, " %3g %7.5f", GSV.phi, sigma);
 262
 263       fprintf(logFP,"\n");
 264     }
 265     else {
 266       fprintf(logFP," O   %12.3f", Zeit);
 267
 268       /* laplace stuff */
 269       if (GTV.phi) fprintf(logFP, " %3g %7.5f", GSV.phi, sigma);
 270
 271       fprintf(logFP," %d %s\n", lmin, costring(GAV.currform));
 272     }
 273     GAV.subi[0] = xsubi[0];
 274     GAV.subi[1] = xsubi[1];
 275     GAV.subi[2] = xsubi[2];
 276     fprintf(logFP, "(%5hu %5hu %5hu)", GAV.subi[0], GAV.subi[1], GAV.subi[2]);
 277     fflush(logFP);
 278
 279     Zeit = 0.0;
 280
 281     /* reset laplace stuff for next trajectory */
 282     sumT = 0.0;
 283     sumK = 0.0;
 284     sumKK = 0.0;
 285     sumD = 0.0;
 286     L = 0.0;
 287     D = 0.0;
 288
 289     /*  highestE = OhighestE = -1000.0; */
 290     reset_nbList();
 291     costring(NULL);
 292     return(1);
 293   }
 294   else {
 295     /* continue simulation */
 296     int flag = 0;
 297     if( (!GTV.silent) && (GSV.currE <= GSV.stopE+GSV.cut) ) {
 298
 299       if (!GTV.lmin || (lmin==1 && strcmp(GAV.prevform, GAV.currform) != 0)) {
 300         char format[64];
 301         flag = 1;
 302         sprintf(format, "%%-%ds %%6.2f %%10.3f", strlen(GAV.farbe_full)+1);
 303         printf(format, costring(GAV.currform), GSV.currE, Zeit);
 304       }
 305
 306       /* laplace stuff */
 307       if (GTV.phi) {
 308         printf(" %8.3f %8.3f %3g", zeitInc, L, D);
 309         L = D = 0.0; /* reset L and D for next structure */
 310       }
 311
 312       if ( flag && GTV.verbose ) {
 313         int ii, jj;
 314         if (next<0) trans='g'; /* growth */
 315         else {
 316           ii = neighbor_list[2*next];
 317           jj = neighbor_list[2*next+1];
 318           if (abs(ii) < GSV.len) {
 319             if ((ii > 0) && (jj > 0)) trans = 'i';
 320             else if ((ii < 0) && (jj < 0)) trans = 'd';
 321             else if ((ii > 0) && (jj < 0)) trans = 's';
 322             else trans = 'S';
 323           } else {
 324             if ((ii > 0) && (jj > 0)) trans = 'I';
 325             else trans = 'D';
 326           }
 327         }
 328         printf(" %4d %c %d", top, trans, lmin);
 329       }
 330       if (flag) printf("\n");
 331     }
 332   }
 333
 334
 335   /* store last lmin seen, so we can avoid printing the same lmin twice */
 336   if (lmin==1)
 337     strcpy(GAV.prevform, GAV.currform);
 338
 339 #if 0
 340   if (lmin==1) {
 341     /* went back to previous lmin */
 342     if (strcmp(GAV.prevform, GAV.currform) == 0) {
 343       if (OhighestE < highestE) {
 344         highestE = OhighestE;  /* delete loop */
 345         strcpy(highestS, OhighestS);
 346       }
 347     } else {
 348       strcpy(GAV.prevform, GAV.currform);
 349       OhighestE = 10000.;
 350     }
 351   }
 352
 353   if ( strcmp(GAV.currform, GAV.startform)==0 ) {
 354     OhighestE = highestE = -1000.;
 355     highestS[0] = 0;
 356   }
 357
 358   /* log highes energy */
 359   if (GSV.currE > highestE) {
 360     OhighestE = highestE;
 361     highestE = GSV.currE;
 362     strcpy(OhighestS, highestS);
 363     strcpy(highestS, GAV.currform);
 364   }
 365 #endif
 366
 367   if (next>=0) update_tree(neighbor_list[2*next], neighbor_list[2*next+1]);
 368   else {
 369     clean_up_rl(); ini_or_reset_rl();
 370   }
 371
 372   reset_nbList();
 373   return(0);
 374 }
 375
 376 /*==========================*/
 377 static void reset_nbList(void) {
 378
 379   top = 0;
 380   totalflux = 0.0;
 381   /*    meanE = 0.0; */
 382   lmin = 1;
 383 }
 384
 385 /*======================*/
 386 void clean_up_nbList(void){
 387
 388   free(neighbor_list);
 389   free(bmf);
 390   free(energies);
 391   fprintf(logFP,"\n");
 392   fclose(logFP);
 393 }
 394
 395 /*======================*/
 396 static void grow_chain(void){
 397   int newl;
 398   /* note Zeit=0 corresponds to chain length GSV.glen */
 399   if (Zeit<(GSV.len+1-GSV.glen) * GSV.grow) return;
 400   newl = GSV.len+1;
 401   Zeit = (newl-GSV.glen) * GSV.grow;
 402   top=0; /* prevent structure move in sel_nb */
 403
 404   if (GSV.len<newl) {
 405     strncpy(GAV.farbe, GAV.farbe_full, newl);
 406     GAV.farbe[newl] = '\0';
 407     strcpy(GAV.startform, GAV.currform);
 408     strcat(GAV.startform, ".");
 409
 410     GSV.len = newl;
 411   }
 412 }
 413
 414 static const char *costring(const char *str) {
 415   static char* buffer=NULL;
 416   static int size=0;
 417   int n;
 418   if (str==NULL) {
 419     if (buffer) {
 420       /* make it possible to free buffer */
 421       free(buffer);
 422       size = 0; buffer = NULL;
 423     }
 424     return NULL;
 425   }
 426   n=strlen(str);
 427   if (n>=size) {
 428     size = n+2;
 429     buffer = realloc(buffer, size);
 430   }
 431   if ((cut_point>0)&&(cut_point<=n)) {
 432     strncpy(buffer, str, cut_point-1);
 433     buffer[cut_point-1] = '&';
 434     strncpy(buffer+cut_point, str+cut_point-1, n-cut_point+1);
 435     buffer[n+1] = '\0';
 436   } else {
 437     strncpy(buffer, str, n+1);
 438   }
 439   return buffer;
 440 }