import io;
import sys;
import math;
import random;

(string o) tcl_eval(string s) "turbine" "0.0" [
  "set <<o>> [ uplevel #0 <<s>> ]"
];

(int o) tcl_eval_int(string s) "turbine" "0.0" [
  "set <<o>> [ uplevel #0 <<s>> ]"
];

(float o) tcl_eval_float(string s) "turbine" "0.0" [
  "set <<o>> [ uplevel #0 <<s>> ]"
];

(string o, float POTENTIAL) run_t (string i, int r, int f, int n, float NEWTEMP, float OLDTEMP) "namdswift" "0.1" "run_t";

main
{
  string ifile = tcl_eval("set ifile");
  int num_replicas = tcl_eval_int("set num_replicas");
  int num_runs = tcl_eval_int("set num_runs");
  int steps_per_run = tcl_eval_int("set steps_per_run");
  float min_temp = tcl_eval_float("set min_temp");
  float max_temp = tcl_eval_float("set max_temp");
  printf("Running %s with %d replicas from %f to %f for %d runs", ifile, num_replicas, min_temp, max_temp, num_runs);
  float TEMPERATURE[int];
  string states[int][int];
  int sources[int][int];
  float POTENTIAL[int][int];
  foreach i in [0:num_replicas-1] {
    TEMPERATURE[i] = min_temp * exp( log(max_temp/min_temp)*(itof(i)/itof(num_replicas-1)) );
    states[1][i], POTENTIAL[1][i] = run_t(ifile, i, 1, steps_per_run, TEMPERATURE[i], 300);
  }
  foreach f in [2:num_runs] {
    if ( f%%2 == 1 ) {
      sources[f][0] = 0;
    }
    if ( (num_replicas+f)%%2 == 1 ) {
      sources[f][num_replicas-1] = num_replicas-1;
    }
    foreach i in [f%%2+1:num_replicas-1:2] {
      BOLTZMAN = 0.001987191;
      dbeta = ((1.0/TEMPERATURE[i-1]) - (1.0/TEMPERATURE[i])) / BOLTZMAN;
      float delta = dbeta * (POTENTIAL[f-1][i] - POTENTIAL[f-1][i-1]);
      boolean doswap = (delta < 0.0) || (exp(-delta) > random());
      printf("frame %d replicas %d %d swap %s\n", f, i-1, i, doswap);
      if ( doswap ) {
        sources[f][i] = i-1;
        sources[f][i-1] = i;
      } else {
        sources[f][i] = i;
        sources[f][i-1] = i-1;
      }
    }
    foreach i in [0:num_replicas-1] {
      int isrc = sources[f][i];
      states[f][i], POTENTIAL[f][i] = run_t(states[f-1][isrc], i, f, steps_per_run, TEMPERATURE[i], TEMPERATURE[isrc]);
    }
  }
}