[Swift-commit] r5502 - SwiftApps

wilde at ci.uchicago.edu wilde at ci.uchicago.edu
Fri Jan 20 15:26:26 CST 2012


Author: wilde
Date: 2012-01-20 15:26:24 -0600 (Fri, 20 Jan 2012)
New Revision: 5502

Added:
   SwiftApps/Makefile
   SwiftApps/README
   SwiftApps/annealing.swift
   SwiftApps/basiclocal.xml
   SwiftApps/beagle.xml
   SwiftApps/cf
   SwiftApps/evolve.sh
   SwiftApps/local.xml
   SwiftApps/math.swift
   SwiftApps/mathtest.swift
   SwiftApps/movie_graph.txt
   SwiftApps/optimizer.cpp
   SwiftApps/optimizer.orig.cpp
   SwiftApps/optimizer.sh
   SwiftApps/optimizer.snap01.cpp
   SwiftApps/optirun.swift
   SwiftApps/script-smaller-a.py
   SwiftApps/sumloss.sh
   SwiftApps/t1.py
   SwiftApps/t2.cpp
   SwiftApps/t2.py
   SwiftApps/t3.py
   SwiftApps/tc
   SwiftApps/testO.sh
Log:
Initial version.

Added: SwiftApps/Makefile
===================================================================
--- SwiftApps/Makefile	                        (rev 0)
+++ SwiftApps/Makefile	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,10 @@
+all:	toptimizer optimizer Optimizer
+
+optimizer:	optimizer.snap01.cpp
+	g++ -static -fopenmp -I boost_1_47_0 -o optimizer optimizer.snap01.cpp
+
+Optimizer:	optimizer.snap01.cpp
+	g++ -static -O3 -fopenmp -I boost_1_47_0 -o Optimizer optimizer.snap01.cpp
+
+toptimizer:	optimizer.cpp
+	g++ -static -fopenmp -I boost_1_47_0 -o toptimizer optimizer.cpp

Added: SwiftApps/README
===================================================================
--- SwiftApps/README	                        (rev 0)
+++ SwiftApps/README	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,42 @@
+=== Files
+
+optirun.swift: replaces top-level py script for the outermost loops
+
+=== How to build
+
+=== How to Run
+
+
+
+=== C++ app flow logic ===
+
+for target in range(58, 1009 (used 209), 50):  // 20 values
+  for i in range(15):
+    #         P0:                 P1:               D:  P2:                   P3:
+    optimizer |0 0 4 50 -1 target | 40000 20 1000 2 | 1 | 2. 0.01 100 0.3 2.3 | 1 1 0 0 0
+
+    multi_annealing( un[NW], T_start, T_end, Target_rejection, Annealing_repeats, 
+                       starting_jump, Results, Counters, params0, Annealing_repeats);
+
+        Res = multi_loss(un, Results, Counters, x); // Initial
+        curr_x   = Res.first;
+        curr_err = Res.second;
+
+        for i = 0 to annealing_cycles - 1 (100 from .py script)
+          for j = 0 to 4 // 5X: 0..4 is fixed constant
+            setParameters (for a round of parallel Annealing_repeats
+            Res = multi_loss(un, Results, Counters, Annealing_repeats);
+              In parallel: for 1 to AnnealingRepeats (1,000; 10,000 desired)
+                group repeats among NWorkers
+                  evolve()
+            setParameters  again here, conditionally? (check)
+
+20 targets - py
+  15 repeats - py
+    1 initial multi_loss: 1000 to 10000 annealing_repeats
+    100 Annealing_cycles (groups of 10? : cycle=10 ) (fast:50)
+       5 repeats (fast: 1)
+         multi_loss: 1000 to 10000 annealing_repeats
+           evolve()  => 2 mins to 10 mins
+
+=== END

Added: SwiftApps/annealing.swift
===================================================================
--- SwiftApps/annealing.swift	                        (rev 0)
+++ SwiftApps/annealing.swift	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,278 @@
+import "math";
+
+type file;
+
+type Res {
+  float loss;
+  float sdev;
+}
+
+global boolean FIX_VARIABLES = true;
+global int var_fixed[] = [1,1,0,0,0];
+global int Nworkers = 1;
+
+(float nx) newx(float x, float dx) 
+{
+    float r = (random()) / (pow(2.0,31.0)-1.0);
+    if (r > 0.5){
+        nx = x + (random())*dx/(pow(2.0,31.0)-1.0);
+    } else {
+        nx = x - (random())*dx/(pow(2.0,31.0)-1.0);
+    }
+}
+
+app (file outfile, file loss) evolve (string args[], file graph)
+{
+  evolve @loss args stdout=@outfile ;  // graph is passed implicitly
+}
+
+app (file x) sumloss(file loss[])
+{
+  sumloss @filenames(loss) stdout=@x;
+}
+
+/*
+
+Program structure:
+
+  main
+    optimizer_sweep() - formerly, python script
+      multi_annealing()
+        multi_loss()
+          evolve()
+        sumloass()
+*/
+
+(file bestfile, file maxfile) multi_annealing (
+    float T_start,
+    float T_end, 
+    float Target_rejection, 
+    int evolve_reruns, 
+    float starting_jump,
+    float params0[],
+    float target_innov,
+    int annealing_cycles)
+{
+    int cycle=10;     # const
+    int NEVOPARAMS=5; # const - 5 params, alpha 1,m through delta, does not include target_innovation
+
+    float rejection[][];  // [i][j] where i is cycle and j is evolve-parameter (alpha_i, alpha_m, beta, gamma, delta)
+
+    float x[][], dx[][], curr_loss[], curr_sdev[];
+    
+    Res mlres[][];
+    mlres[0][0] = multi_loss( params0, target_innov, evolve_reruns ); // Only done once, not 5x; serves for all evolve-params ???
+
+    foreach j in [0:NEVOPARAMS-1] {
+        x[0][j]=params0[j];
+        dx[0][j] = starting_jump;
+	rejection[0][j] = 0.0;
+        curr_loss[j] = mlres[0][0].loss;
+        curr_sdev[j] = mlres[0][0].sdev;
+    }
+
+    foreach i in [1:annealing_cycles] {
+        // set new temperature, rejection threshold, and dx values for this cycle
+        float temperature = T_start*exp( @tofloat(i)*(jlog(T_end)-jlog(T_start))/@tofloat(annealing_cycles));
+        tracef("....T =%f\n", temperature);
+	// On each new "major" cycle within the annealing_cycles (other than the first) set new rejection and dx values
+        if ( i %% cycle == 1 && i > 1 ){ 
+            foreach  k in [0:NEVOPARAMS-1] {
+                float newrejection = rejection[i-1][k] / @tofloat(cycle);
+                if (newrejection > 0.0){
+                    dx[i][k] = dx[i-1][k] / (newrejection / Target_rejection);
+                    rejection[i][k]=0.0;
+                }
+                else{
+                    dx[i][k] = dx[i-1][k] * 2.0;
+                    rejection[i][k]=rejection[i-1][k];
+                }
+                trace ("Recomputed rejection: i=%d k=%d dx[i][k]=%f\n", i, k, dx[i][k]);
+            }
+        }
+        foreach j in [0:NEVOPARAMS-1] { // Try a new value for each non-fixed param; then write results and accept or reject
+            float try_x[];
+            int curr = (i * NEVOPARAMS) + j;
+            int prev = curr-1;
+            if ( FIX_VARIABLES || var_fixed[j]==0) {
+                foreach k in [0:NEVOPARAMS] { // Select the evolve params to try
+                    if ( k < j ) {
+                        try_x[k] = x[i][k]; // already set x[i][k]
+                    }
+                    else {
+                        if ( k == j ) {
+                            try_x[k] = newx(x[i][j],dx[i][j]); // permute x[i][j]
+                        }
+                        else {
+                            try_x[k] = x[i-1][k]; // use x[i-1][k] (from prior cycle)
+                        }
+                    }
+                }
+
+                mlres[i][j] = multi_loss(try_x, target_innov, evolve_reruns);
+
+                float ratio = min(1.0, exp( -(mlres[i][j].loss-curr_loss[prev]) /temperature));
+                float r = (random()) / (pow(2.0,31.0)-1.0);  // why all the 2^31's ???
+
+                float ALOT=100000000000.0; // 100,000,000,000. = 10^11
+                if (mlres[i][j].loss < ALOT) {  // does this ever occur? if so did we want to still do the ratio computation above???
+                    trace("filestr.open (best_opt_some.txt, ofstream::app);");
+                    # un[0]->get_target() Res.loss un[0]->get_parameter[0:4] Res.error
+                    # filestr.open ("max_dist.txt", ofstream::app);
+                }
+                else { tracef("Loss %f > ALOT at [i][j] = [%d][%d]\n", mlres[i][j].loss, i ,j); }
+                if (r > ratio) {  // Reject
+                    x[i][j] = x[i-1][j];
+                    rejection[i][j] = rejection[i][j] + 1.0;  // Is this correct? incr rejection? 
+                    curr_loss[curr] = curr_loss[prev];
+                    curr_sdev[curr] = curr_sdev[prev];
+                }
+                else {           // Accept
+trace("Accepting try_x[j],j=",j);
+                    x[i][j] = try_x[j];  # FIXME: is this right?  Just assign the j'th try_x, or the whole try_x param set?
+trace("Accepting try_x[j],x[i][j]=",x[i][j]);
+                    curr_loss[curr] = mlres[i][j].loss;
+                    curr_sdev[curr] = mlres[i][j].sdev;
+                }
+            }
+            else {
+                x[i][j] = x[i-1][j];
+                curr_loss[curr] = curr_loss[prev];
+                curr_sdev[curr] = curr_sdev[prev];
+           }
+        }
+    }
+}
+
+(Res r) multi_loss( float x[], float target_innov, int evolve_reruns )
+{
+  file rfile[];
+  foreach i in [1:evolve_reruns] {  // repeats of the evolove() - same as n_reruns
+    file outfile; // FIXME: map and save in future 
+    string args[] = [ // FIXME: move this to a setargs() function
+      //    alpha_i        alpha_m        beta           gamma          delta          target_innov
+            @strcat(x[0]), @strcat(x[1]), @strcat(x[2]), @strcat(x[3]), @strcat(x[4]), @strcat(target_innov),
+      
+      //     n_epochs n_steps evolve_reruns           range
+      //    "40000",  "20",   @strcat(evolve_reruns), "2",
+            "40000",  "20",   "1", "                   2",  # Set reruns to 1 for the cpp code; we do the re-runs in Swift for now.
+      
+      //    verbose_level
+            "1",
+      
+      //    T_start T_end   Annealing_steps Target_rejection Starting_jump
+            "2.",   "0.01", "2",            "0.3",           "2.3",
+      
+      //    FREEZE: alpha_i alpha_m beta gamma delta
+                    "1",    "1",    "0", "0",  "0",
+      
+      //   operation-code:(m,a) Nworkers
+           "m",                 @strcat(Nworkers) ];
+
+    file graph <"movie_graph.txt">;
+    (outfile, rfile[i]) = evolve(args,graph); 
+  }
+  file sumfile = sumloss(rfile);
+  r = readData(sumfile);
+
+}
+
+optimizer_sweep() // Implements logic of python driver script
+{
+  int minrange=58;
+  int maxrange=59;
+  //int maxrange=1009;
+  //int maxrange=209;
+  int rangeinc=50;
+
+  // fixme: add provision for random priming and random param values when x[i] == -100 (see optimizer.cpp main())
+
+  int nreps=1; # 15
+
+//    file bestfile <single_file_mapper; file=@strcat("output/T",target,".R",rep,".best_opt_some")>;
+//    file maxfile <single_file_mapper; file=@strcat("output/T",target,".R",rep,".max_dist")>;
+
+  foreach target_innov in [minrange:maxrange:rangeinc] {
+    foreach rep in [1:nreps] {
+      file outfile;  // <single_file_mapper; file=@strcat("output/T",target_innov,".R",rep,".out")>;
+      file lossfile; // <single_file_mapper; file=@strcat("output/T",target_innov,".R",rep,".loss_data")>;
+/*
+      (outfile,lossfile) = multi_annealing(
+         T_start          = 2.0,
+         T_end            = 0.01,
+         Target_rejection = 0.3,
+         evolve_reruns    = 2,
+         starting_jump    = 2.3,
+         params0[]        = [0.0, 0.0, 4.0, 50.0, -1.0],
+         @tofloat(@strcat(target_innov)), // fixme: would @tofloat(int) work?
+         annealing_cycles = 1);   // fixme: clarify: annealing repeats, steps, and cycles
+*/
+      (outfile,lossfile) = multi_annealing(
+         2.0,
+         0.01,
+         0.3,
+         2,
+         2.3,
+         [0.0, 0.0, 4.0, 50.0, -1.0],
+         @tofloat(@strcat(target_innov)), // fixme: would @tofloat(int) work?
+         1);   // fixme: clarify: annealing repeats, steps, and cycles
+    }
+  }
+}
+
+main()
+{
+  optimizer_sweep();
+}
+
+main();
+
+/*
+
+Program structure:
+
+  main
+    optimizer_sweep()
+      multi_annealing()
+        multi_loss()
+          evolve()
+        sumloass()
+
+Example parameter sets:
+
+for target in range(58,59,50):
+  for i in range(1):
+    args="./toptimizer 0 0 4 50 -1 "+target+" 40000 20 75    2 1 2. 0.01 2 0.3 2.3 1 1 1 0 0 m # > out.T"+str(target)+".i"+str(i)
+    os.system(args);
+
+      string fastargs1[] = [
+        "0", "0", "4", "50", "-1", @strcat(target),
+        "40000", "20", "1000", "2",
+        "1",
+        "2.", "0.01", "100", "0.3", "2.3",
+        "1", "1", "0", "0", "0"];
+      string fastargs2[] = [
+        "0", "0", "4", "50", "-1", @strcat(target),
+        "40000", "20", "1000", "2",
+        "1",
+        "2.", "0.01",  "5", "0.3", "2.3",
+        "1", "1", "0", "0", "0", "m"];
+      string fastargs3[] = [
+        "0", "0", "4", "50", "-1", @strcat(target),
+        "40000", "20", @strcat(repeats), "2",
+        "1",
+        "2.", "0.01",  "2", "0.3", "2.3",
+        "1", "1", "0", "0", "0", "m"];
+*/
+
+(string args[]) setargs()
+{
+  // string longargs[] = @strcat("0 0 4 50 -1 ",target," 40000 20 1000 2 1 2. 0.01 100 0.3 2.3 1 1 0 0 0 m");
+
+  //  [alpha_i alpha_m beta gamma delta target_innov
+  //  [n_epochs n_steps n_reruns] [range]
+  //  [verbose_level]
+  //  [T_start T_end Annealing_steps Target_rejection Starting_jump]
+  //  [FREEZE_alpha_i FREEZE_alpha_m FREEZE_beta FREEZE_gamma FREEZE_delta] [operation-code:(m,a) Nworkers]
+}
+

Added: SwiftApps/basiclocal.xml
===================================================================
--- SwiftApps/basiclocal.xml	                        (rev 0)
+++ SwiftApps/basiclocal.xml	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,10 @@
+<config>
+
+   <pool handle="localhost">
+     <execution provider="local" url="none" />
+     <filesystem provider="local"/>
+     <workdirectory>/home/wilde/swift/lab/swiftwork</workdirectory>
+     <profile namespace="karajan" key="jobThrottle">0</profile>
+   </pool>
+
+</config>

Added: SwiftApps/beagle.xml
===================================================================
--- SwiftApps/beagle.xml	                        (rev 0)
+++ SwiftApps/beagle.xml	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,60 @@
+<config>
+
+<!-- From Glen: used w/ "Swift svn swift-r4813 (swift modified locally) cog-r3175" -->
+
+  <pool handle="beagle">
+
+    <execution provider="coaster" jobmanager="local:pbs" url="none"/>
+    <profile namespace="globus" key="providerAttributes">pbs.aprun;pbs.mpp;depth=24</profile>
+    <profile key="jobsPerNode" namespace="globus">1</profile>
+
+    <profile namespace="env" key="OMP_NUM_THREADS">24</profile>
+    <profile namespace="globus" key="maxwalltime">02:00:00</profile>
+    <profile namespace="globus" key="maxTime">14400</profile>
+    <profile namespace="globus" key="slots">20</profile>
+    <profile namespace="globus" key="nodeGranularity">1</profile>
+    <profile namespace="globus" key="maxNodes">1</profile>
+    <profile namespace="globus" key="lowOverAllocation">100</profile>
+    <profile namespace="globus" key="highOverAllocation">100</profile>
+    <profile namespace="karajan" key="jobThrottle">.15</profile>
+    <profile namespace="karajan" key="initialScore">10000</profile>
+
+    <profile namespace="globus" key="project">CI-MCB000119</profile>
+    <profile namespace="globus" key="queue">route</profile>
+
+    <filesystem provider="local"/>
+    <workdirectory >/lustre/beagle/wilde/swiftwork</workdirectory>
+  </pool>
+
+<!-- From Justin's swift-devel page: 
+
+<import file="sys.xml"/>
+<set name="wdir" value="/lustre/beagle/{user.name}/work"/>
+<echo message="setting workDirectory to: {wdir}"/>
+
+
+<pool handle="beagle-pbs">
+  <execution jobmanager="local:pbs" provider="coaster" url="none"/>
+  <profile namespace="globus" key="maxWallTime">1</profile>
+  <profile namespace="globus" key="maxTime">7200</profile>
+
+  <profile namespace="globus" key="providerAttributes">
+    pbs.aprun;pbs.mpp;depth=24
+  </profile>
+  <profile key="jobsPerNode" namespace="globus">24</profile>
+  <profile key="slots" namespace="globus">1</profile>
+  <profile key="nodeGranularity" namespace="globus">1</profile>
+  <profile key="maxNodes" namespace="globus">1</profile>
+  <profile key="queue" namespace="globus">batch</profile>
+  <profile key="jobThrottle" namespace="karajan">5.99</profile>
+  <profile key="initialScore" namespace="karajan">10000</profile>
+  <profile namespace="globus" key="project">_PROJECT_</profile>
+  <profile namespace="globus" key="project">_QUEUE_</profile>
+
+  <filesystem provider="local" url="none" />
+  <workdirectory>{wdir}</workdirectory>
+</pool>
+
+-->
+
+</config>

Added: SwiftApps/cf
===================================================================
--- SwiftApps/cf	                        (rev 0)
+++ SwiftApps/cf	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,7 @@
+wrapperlog.always.transfer=true
+sitedir.keep=true
+execution.retries=3
+lazy.errors=true
+status.mode=provider
+use.provider.staging=false
+provider.staging.pin.swiftfiles=false

Added: SwiftApps/evolve.sh
===================================================================
--- SwiftApps/evolve.sh	                        (rev 0)
+++ SwiftApps/evolve.sh	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,6 @@
+#1 /bin/sh
+datafile=$1
+touch multi_loss.data
+shift 1
+$(dirname $0)/toptimizer $* 2>&1
+mv multi_loss.data $datafile


Property changes on: SwiftApps/evolve.sh
___________________________________________________________________
Added: svn:executable
   + 

Added: SwiftApps/local.xml
===================================================================
--- SwiftApps/local.xml	                        (rev 0)
+++ SwiftApps/local.xml	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,9 @@
+<config>
+  <pool handle="localhost" >
+    <execution provider="local" url="none" />
+    <profile namespace="karajan" key="jobThrottle">0.23</profile>
+    <profile namespace="karajan" key="initialScore">10000</profile>
+    <filesystem provider="local"/>
+    <workdirectory>/tmp/wilde/swiftwork</workdirectory>
+  </pool>
+</config>

Added: SwiftApps/math.swift
===================================================================
--- SwiftApps/math.swift	                        (rev 0)
+++ SwiftApps/math.swift	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,53 @@
+
+(float result) sin(float x)
+{
+   result = @java("java.lang.Math", "sin", x);
+}
+
+(float result) exp(float x)
+{
+   result = @java("java.lang.Math", "exp", x);
+}
+
+(float result) jlog(float x)
+{
+   result = @java("java.lang.Math", "log", x);
+}
+
+(float result) log10(float x)
+{
+   result = @java("java.lang.Math", "log10", x);
+}
+
+(float result) ceil(float x)
+{
+   result = @java("java.lang.Math", "ceil", x);
+}
+
+(float result) floor (float x)
+{
+   result = @java("java.lang.Math", "floor", x);
+}
+
+(float result) min (float a, float b)
+{
+   //result = @java("java.lang.Math", "min", a, b);
+   if ( a < b ) {
+       result = a;
+   }
+   else {
+       result = b;
+   }
+   tracef("min: result=%f\n",result);
+}
+
+(float result) pow (float x, float y)
+{
+   result = @java("java.lang.Math", "pow", x, y);
+}
+
+(float result) random ()
+{
+  result = @java("java.lang.Math","random");
+}
+

Added: SwiftApps/mathtest.swift
===================================================================
--- SwiftApps/mathtest.swift	                        (rev 0)
+++ SwiftApps/mathtest.swift	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,18 @@
+import "math";
+
+float a = 0.5;
+tracef("random(): %f\n", random());
+tracef("sin(%f): %f\n", a, sin(a));
+tracef("jlog(%f): %f\n", a, jlog(a));
+tracef("log10(%f): %f\n", a, log10(a));
+tracef("exp(%f): %f\n", a, exp(a));
+tracef("ceil(%f): %f\n", 1.23, ceil(1.23));
+tracef("floor(%f): %f\n", 1.23, floor(1.23));
+tracef("pow(%f,%f): %f\n", 2.0, 4.0, pow(2.0,4.0));
+tracef("min(%f,%f): %f\n", 2.0, 4.0, min(2.0,4.0));
+tracef("min(%f,%f): %f\n", 2.0, 1.5, min(2.0,1.5));
+
+
+
+
+

Added: SwiftApps/movie_graph.txt
===================================================================
--- SwiftApps/movie_graph.txt	                        (rev 0)
+++ SwiftApps/movie_graph.txt	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,1010 @@
+500 1008
+0 56
+0 436
+0 446
+1 256
+1 476
+1 260
+1 262
+1 9
+1 266
+1 13
+1 142
+1 273
+1 388
+1 30
+1 36
+1 37
+1 423
+1 168
+1 425
+1 255
+1 175
+1 307
+1 58
+1 316
+1 319
+1 323
+1 68
+1 455
+1 216
+1 82
+1 88
+1 346
+1 348
+1 225
+1 228
+1 236
+1 365
+1 111
+1 424
+1 370
+1 246
+1 376
+1 122
+1 191
+1 252
+1 382
+1 127
+2 3
+3 262
+3 263
+3 392
+3 172
+3 138
+3 395
+3 13
+3 18
+3 150
+3 408
+3 154
+3 156
+3 30
+3 291
+3 37
+3 168
+3 300
+3 52
+3 158
+3 310
+3 265
+3 187
+3 61
+3 319
+3 65
+3 322
+3 54
+3 130
+3 71
+3 332
+3 461
+3 397
+3 419
+3 473
+3 346
+3 95
+3 352
+3 98
+3 103
+3 107
+3 414
+3 368
+3 498
+3 499
+4 224
+4 449
+4 419
+4 356
+4 72
+4 330
+4 173
+4 80
+4 242
+4 180
+4 24
+4 138
+5 33
+5 459
+5 484
+5 106
+5 299
+5 430
+5 16
+5 336
+5 50
+5 147
+5 406
+5 87
+5 91
+5 190
+6 7
+6 264
+6 137
+6 276
+6 149
+6 279
+6 153
+6 29
+6 415
+6 293
+6 166
+6 296
+6 42
+6 44
+6 48
+6 49
+6 438
+6 55
+6 314
+6 69
+6 454
+6 463
+6 84
+6 85
+6 470
+6 349
+6 94
+6 229
+6 233
+6 364
+6 117
+6 248
+6 380
+6 426
+7 129
+7 131
+7 470
+7 136
+7 471
+7 271
+7 402
+7 19
+7 404
+7 282
+7 286
+7 491
+7 35
+7 421
+7 166
+7 295
+7 427
+7 428
+7 301
+7 349
+7 436
+7 439
+7 56
+7 445
+7 446
+7 320
+7 123
+7 451
+7 69
+7 454
+7 78
+7 463
+7 212
+7 213
+7 342
+7 343
+7 474
+7 93
+7 483
+7 357
+7 102
+7 235
+7 237
+7 239
+7 371
+7 379
+8 112
+8 164
+8 10
+8 12
+8 398
+8 176
+8 120
+8 178
+8 373
+8 24
+8 124
+9 262
+10 11
+10 12
+10 47
+10 335
+10 124
+11 96
+11 289
+11 66
+11 283
+11 97
+11 361
+11 47
+11 432
+11 335
+11 215
+11 27
+11 124
+11 189
+12 112
+12 434
+12 47
+13 32
+13 130
+13 260
+13 37
+13 39
+13 202
+13 305
+13 461
+13 145
+13 275
+13 448
+13 54
+13 183
+13 344
+13 36
+14 387
+14 267
+14 155
+14 36
+14 165
+14 51
+14 308
+14 440
+14 441
+14 59
+14 444
+14 199
+14 328
+14 462
+14 81
+14 210
+14 214
+14 91
+14 99
+14 60
+14 106
+14 109
+14 110
+14 244
+15 193
+15 73
+15 400
+15 341
+15 27
+15 31
+16 33
+16 258
+16 353
+16 329
+16 204
+16 207
+16 336
+16 108
+16 50
+16 22
+16 119
+16 169
+16 381
+16 223
+17 321
+17 482
+17 197
+17 276
+17 149
+17 315
+17 94
+18 32
+18 461
+18 305
+18 498
+18 86
+18 25
+19 320
+19 421
+19 334
+19 314
+19 349
+19 446
+20 48
+20 28
+20 83
+20 44
+20 45
+21 288
+21 176
+21 347
+21 325
+21 392
+21 490
+21 366
+21 272
+21 209
+21 115
+21 182
+21 247
+21 184
+21 25
+21 217
+21 27
+21 86
+21 478
+22 353
+23 257
+23 100
+23 25
+24 134
+24 146
+24 302
+24 151
+24 290
+24 291
+24 43
+24 173
+24 46
+24 176
+24 180
+24 245
+24 453
+24 457
+24 208
+24 338
+24 363
+24 116
+24 373
+24 126
+25 32
+25 128
+25 257
+25 86
+26 161
+26 357
+26 193
+26 78
+26 341
+26 69
+27 361
+27 490
+27 333
+27 272
+27 312
+27 468
+27 41
+27 184
+27 429
+27 124
+27 341
+28 64
+28 354
+28 67
+28 135
+28 200
+28 42
+28 75
+28 45
+28 206
+28 48
+28 359
+28 53
+28 238
+29 229
+29 199
+29 296
+29 233
+29 42
+29 274
+29 117
+29 152
+29 155
+29 62
+30 386
+30 225
+30 76
+30 322
+30 367
+30 318
+31 200
+31 73
+31 333
+31 494
+31 45
+31 312
+32 257
+32 305
+32 101
+32 486
+32 39
+32 231
+32 79
+32 145
+32 275
+33 459
+33 336
+33 372
+33 119
+33 409
+34 283
+35 491
+35 243
+35 237
+35 326
+36 260
+36 266
+36 285
+36 292
+36 39
+36 424
+36 171
+36 181
+36 183
+36 440
+36 58
+36 323
+36 328
+36 81
+36 82
+36 87
+36 99
+36 228
+36 101
+36 111
+37 168
+37 202
+38 92
+39 448
+39 101
+39 486
+39 231
+39 145
+39 344
+40 345
+40 355
+40 179
+40 59
+40 125
+41 160
+41 226
+41 230
+41 104
+41 361
+41 492
+41 341
+41 215
+41 121
+41 189
+42 64
+42 67
+42 199
+42 359
+42 109
+42 206
+42 143
+42 48
+42 274
+42 117
+44 193
+44 69
+44 161
+44 73
+44 48
+44 83
+44 279
+45 200
+45 73
+45 83
+46 242
+46 485
+46 126
+46 173
+48 64
+48 206
+50 259
+50 261
+50 329
+50 299
+50 204
+50 430
+50 207
+50 405
+50 59
+50 159
+51 67
+51 101
+51 135
+51 210
+51 244
+51 214
+51 441
+51 249
+52 368
+52 217
+52 347
+52 86
+53 101
+53 327
+53 75
+53 238
+53 79
+53 249
+54 130
+55 84
+57 162
+57 331
+58 266
+58 82
+59 384
+59 259
+59 132
+59 267
+59 496
+59 355
+59 152
+59 91
+59 221
+59 159
+60 91
+60 165
+62 152
+62 155
+62 221
+63 456
+63 334
+63 303
+63 472
+63 314
+63 380
+65 154
+65 499
+65 374
+65 103
+66 96
+66 160
+66 215
+67 359
+67 135
+67 109
+67 143
+67 214
+69 161
+69 166
+69 239
+69 404
+69 85
+69 279
+69 153
+70 410
+70 467
+71 130
+72 138
+73 193
+74 241
+74 82
+74 188
+74 133
+74 469
+75 140
+75 79
+75 90
+75 350
+76 386
+77 120
+77 306
+77 452
+77 254
+77 311
+78 239
+78 357
+79 257
+79 100
+79 101
+80 138
+81 99
+81 171
+82 256
+82 133
+82 348
+82 205
+82 372
+82 181
+82 119
+82 409
+82 255
+82 316
+82 122
+82 469
+84 277
+84 264
+84 149
+84 380
+84 415
+85 153
+85 166
+86 247
+86 217
+87 181
+87 147
+87 285
+89 192
+89 200
+89 333
+89 468
+89 443
+89 222
+90 272
+90 100
+90 350
+91 165
+91 106
+91 430
+91 406
+91 159
+92 258
+93 136
+93 464
+93 493
+93 326
+94 482
+94 233
+94 114
+94 276
+95 107
+95 265
+95 186
+95 332
+96 289
+96 324
+96 102
+96 401
+97 283
+97 220
+97 434
+97 164
+99 101
+99 171
+99 244
+99 251
+100 272
+101 231
+101 466
+101 244
+101 249
+101 251
+102 141
+102 237
+102 369
+102 243
+102 212
+102 121
+102 340
+102 411
+104 121
+104 230
+105 396
+105 309
+106 147
+106 308
+106 406
+106 444
+106 190
+109 199
+109 143
+109 387
+109 214
+110 267
+112 164
+113 195
+114 233
+114 482
+114 315
+114 412
+116 465
+117 296
+118 163
+119 188
+119 409
+120 297
+120 363
+120 398
+120 254
+121 160
+121 230
+121 492
+121 369
+121 212
+122 256
+122 316
+124 176
+125 394
+125 497
+125 179
+125 351
+125 479
+126 465
+126 290
+126 339
+129 491
+130 397
+131 464
+131 136
+132 383
+133 273
+133 232
+133 205
+133 241
+133 316
+134 203
+135 249
+135 354
+136 464
+136 474
+137 276
+138 419
+138 414
+138 158
+139 169
+140 200
+140 222
+140 350
+144 334
+146 363
+147 308
+147 285
+148 175
+149 456
+149 407
+151 422
+152 488
+152 233
+152 221
+152 351
+154 187
+154 358
+154 158
+155 199
+155 267
+155 221
+155 437
+156 332
+156 318
+156 310
+157 436
+157 301
+158 396
+158 187
+158 414
+159 430
+160 215
+161 193
+163 220
+164 297
+164 398
+164 220
+166 427
+167 201
+167 442
+167 170
+167 383
+168 319
+169 280
+169 270
+169 198
+170 353
+170 403
+170 269
+170 207
+170 211
+173 242
+174 194
+174 195
+174 362
+174 177
+174 250
+174 219
+174 412
+175 460
+175 246
+175 447
+176 291
+176 294
+176 366
+176 178
+176 182
+177 219
+179 345
+179 394
+181 372
+182 184
+182 429
+182 294
+184 429
+185 358
+188 280
+188 469
+189 215
+192 272
+192 443
+192 468
+192 287
+193 400
+195 488
+195 233
+196 389
+196 317
+197 321
+198 280
+199 387
+200 333
+200 222
+201 284
+203 311
+207 353
+207 329
+207 269
+207 240
+212 230
+212 341
+212 439
+213 349
+214 462
+215 226
+216 388
+216 236
+217 347
+218 465
+219 450
+219 394
+219 250
+219 298
+219 479
+220 297
+220 283
+221 267
+221 437
+222 443
+222 350
+225 236
+227 472
+227 321
+227 315
+227 412
+227 390
+228 424
+228 260
+230 341
+232 280
+232 241
+233 480
+233 412
+234 243
+234 340
+238 249
+238 354
+240 353
+241 280
+242 356
+243 385
+243 326
+244 251
+246 382
+249 354
+249 327
+249 466
+250 488
+250 479
+250 351
+252 376
+253 384
+253 269
+254 363
+254 302
+254 311
+258 353
+258 317
+258 477
+259 384
+259 405
+259 269
+260 424
+262 352
+262 346
+263 392
+263 300
+263 325
+264 438
+267 437
+268 413
+269 384
+269 383
+270 381
+270 431
+272 443
+272 350
+272 287
+273 316
+276 377
+278 390
+280 417
+281 376
+282 391
+283 289
+283 432
+285 308
+285 292
+288 325
+288 478
+289 401
+289 413
+291 325
+291 478
+292 328
+292 308
+297 398
+297 495
+302 453
+302 311
+304 493
+307 388
+308 328
+309 396
+311 453
+312 333
+312 494
+313 375
+314 487
+314 364
+314 380
+315 482
+315 412
+319 346
+321 456
+326 491
+326 435
+327 466
+328 440
+331 472
+332 408
+332 489
+333 468
+334 446
+337 380
+338 465
+341 400
+343 436
+345 394
+345 399
+347 392
+347 368
+349 364
+349 470
+351 496
+351 479
+353 393
+358 360
+358 396
+365 375
+366 478
+368 392
+372 409
+374 467
+378 411
+379 436
+380 456
+380 415
+382 433
+383 384
+394 481
+394 479
+396 449
+416 420
+416 422
+417 447
+418 446
+450 481
+456 475
+458 471
+463 470
+472 475
+496 497
+

Added: SwiftApps/optimizer.cpp
===================================================================
--- SwiftApps/optimizer.cpp	                        (rev 0)
+++ SwiftApps/optimizer.cpp	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,1814 @@
+//
+//  main.cpp
+//  optimizer
+//
+//  Created by Andrey Rzhetsky on 4/11/11.
+//  Copyright 2011 University of Chicago. All rights reserved.
+//
+
+#define MAXNworkers 24
+int Nworkers=MAXNworkers;
+char operation = 'n'; // n: normal; m: multi_loss; a: analyze and generate next annealing parameter set.
+
+#include <fstream>
+#include <sstream>
+#include <iostream>
+#include <stdio.h>
+#include <time.h>
+#include <ctime>    
+#include <algorithm>
+#include <string>
+
+#include <stdio.h>
+#include <sys/param.h>
+#include <sys/time.h>
+#include <sys/types.h>
+
+// #include <dispatch/dispatch.h>
+#include <fstream>
+
+
+#include <stdlib.h>
+#include <boost/numeric/ublas/io.hpp>
+#include <boost/graph/graph_traits.hpp>
+#include <boost/graph/dijkstra_shortest_paths.hpp>
+#include <boost/graph/loop_erased_random_walk.hpp>
+#include <boost/graph/random.hpp>
+#include <boost/property_map/property_map.hpp>
+#include <boost/graph/graph_concepts.hpp>
+#include <boost/graph/properties.hpp>
+
+#include <boost/graph/graph_traits.hpp>
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/adjacency_matrix.hpp>
+
+#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
+#define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
+#include <boost/graph/random.hpp>
+#include <boost/random/geometric_distribution.hpp>
+#include <boost/random/uniform_01.hpp>
+#include <boost/random.hpp>
+#include <boost/random/linear_congruential.hpp>
+#include <boost/random/uniform_int.hpp>
+#include <boost/random/uniform_real.hpp>
+#include <boost/random/variate_generator.hpp>
+#include <boost/generator_iterator.hpp>
+#include <boost/lexical_cast.hpp>
+
+#define INT_INFINITY 2147483647
+#define NEVOPARAMS 5
+
+#define FIX_VARIABLES 1
+
+using namespace boost;
+using namespace std;
+using namespace boost::numeric::ublas;
+
+static int max_dist=0;
+
+typedef boost::adjacency_matrix<boost::directedS> Graph;
+typedef std::pair<int,int> Edge;
+typedef boost::graph_traits<Graph> GraphTraits;
+typedef boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::strict_upper> prob;
+typedef boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::strict_upper> pathlength;
+typedef boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
+
+namespace std {
+	using ::time;
+}
+
+static int var_fixed[NEVOPARAMS] = {1, 0, 1, 1, 0};
+
+typedef boost::minstd_rand base_generator_type;
+typedef adjacency_list < listS, vecS, directedS,
+no_property, property < edge_weight_t, int > > graph_t;
+typedef graph_traits < graph_t >::vertex_descriptor vertex_descriptor;
+typedef graph_traits < graph_t >::edge_descriptor edge_descriptor;
+
+
+//================================================
+string strDouble(double number)
+{
+    stringstream ss;//create a stringstream
+    ss << number;//add number to the stream
+    return ss.str();//return a string with the contents of the stream
+}
+
+//================================================
+
+double gaussian(double sigma)
+{
+    double GaussNum = 0.0;
+    int NumInSum = 10;
+    for(int i = 0; i < NumInSum; i++)
+    {
+        GaussNum += ((double)rand()/(double)RAND_MAX - 0.5);
+    }
+    GaussNum = GaussNum*sqrt((double)12/(double)NumInSum);
+    
+    
+    return GaussNum;
+    
+}
+
+
+
+//=================================================
+double diffclock(clock_t clock1,clock_t clock2)
+{
+	double diffticks=clock1-clock2;
+	double diffms=(diffticks)/CLOCKS_PER_SEC;
+	return diffms;
+}
+
+//================================================
+//================================================================
+double get_new_x(double x, double dx){
+    
+    double new_x;
+    // boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
+    double r = rand()/(double)(pow(2.,31)-1.);
+    
+    if (r > 0.5){            
+        new_x = x + rand()*dx/(double)(pow(2.,31)-1.);
+    } else {            
+        new_x = x - rand()*dx/(double)(pow(2.,31)-1.);
+    }
+    
+    return new_x;
+    
+}
+
+
+//===============================================   
+string string_wrap(string ins, int mode){
+    
+    std::ostringstream s;
+    
+    switch(mode){
+        case 0:
+            s << "\033[1;29m" << ins << "\033[0m";
+            break;
+        case 1:
+            s << "\033[1;34m" << ins << "\033[0m";
+            break;
+        case 2:
+            s << "\033[1;44m" << ins << "\033[0m";
+            break;
+        case 3:
+            s << "\033[1;35m" << ins << "\033[0m";
+            break;
+        case 4:
+            s << "\033[1;33;44m" << ins << "\033[0m";
+            break;
+        case 5:
+            s << "\033[1;47;34m" << ins << "\033[0m";
+            break;
+        case 6:
+            s << "\033[1;1;31m" << ins << "\033[0m";
+            break;
+        case 7:
+            s << "\033[1;1;33m" << ins << "\033[0m";
+            break;
+        case 8:
+            s << "\033[1;1;43;34m" << ins << "\033[0m";
+            break;
+        case 9:
+            s << "\033[1;1;37m" << ins << "\033[0m";
+            break;
+        case 10:
+            s << "\033[1;30;47m" << ins << "\033[0m";
+            break;
+        default:
+            s << ins;
+    }
+    
+    return s.str();
+}
+
+
+//===============================================
+string wrap_double(double val, int mode){
+    
+    std::ostringstream s;
+    s << string_wrap(strDouble(val),mode);
+    
+    return s.str();
+}
+
+
+
+//===============================================
+const     
+string i2string(int i){
+    
+    std::ostringstream s;
+    s << "worker" 
+    << lexical_cast<std::string>(i);
+    
+    return s.str();
+    
+}
+
+//===============================================
+char* i2char(int i){
+    
+    std::ostringstream s;
+    s << "worker" 
+    << lexical_cast<std::string>(i);
+    
+    char* a=new char[s.str().size()+1];
+    memcpy(a,s.str().c_str(), s.str().size());
+    
+    return a;
+}
+
+
+template <class T>
+bool from_string(T& t, 
+                 const std::string& s, 
+                 std::ios_base& (*f)(std::ios_base&))
+{
+  std::istringstream iss(s);
+  return !(iss >> f >> t).fail();
+}
+
+//================================================
+class Universe {
+	
+private:
+	
+	double alpha_i;
+	double alpha_m;
+	double beta;
+	double gamma;
+	double delta;
+	
+    double TargetNovelty;
+    double CumulativeRelativeLoss;
+    double CRLsquare;
+    string id;
+    
+    
+	int N_nodes;
+	int M_edges;
+	
+	int N_epochs;
+	int N_steps;
+	int N_repeats;
+	
+	int current_epoch;
+	double current_loss;
+	int current_repeat;
+    double current_novelty;
+	
+	int mode_identify_failed;
+    int verbose_level; // 0 is silent, higher is more
+	
+	double k_max;
+	
+	graph_t Full_g;
+	
+	double **Prob;
+	double **Tried;
+	double **Dist;
+	double **Final;
+    double **EdgeIndex;
+	double *Rank;
+	
+    base_generator_type generator;	
+    boost::uniform_real<> uni_dist;
+    boost::geometric_distribution<double> geo;
+    
+public:
+	
+    
+    
+	//======  Constructor ======
+	Universe(const std::string FileToOpen, int Epochs, int Steps, int Repeats, int identify_failed, double target, const std::string idd)
+	{
+		//typedef array_type2::index index2;
+		
+		
+		std::ifstream inFile;
+        //string line;
+        
+        //-------------------------------
+        
+        base_generator_type gene(42u);
+        generator = gene;
+        generator.seed(static_cast<unsigned int>(std::time(0)));
+        boost::uniform_real<> uni_d(0,1);
+        uni_dist = uni_d;
+        
+        //--------------------------------
+        
+		int i, k;
+        int x, y;
+		Edge* edge_array_mine;
+		int num_arcs_mine, num_nodes_mine;
+		int* weights_mine;
+		
+        TargetNovelty = target;
+        CumulativeRelativeLoss = 0.;
+        CRLsquare = 0.;
+        
+		
+		N_epochs  = Epochs;
+		N_steps   = Steps;
+		N_repeats = Repeats;
+		
+		current_epoch = 0;
+		current_loss = 0.;
+		current_repeat = 0;
+        
+        id = idd;
+        
+        verbose_level = 1;
+		
+		mode_identify_failed = identify_failed;
+        
+		
+		//-------------------------------
+		// The first pass though file with the graph
+		inFile.open(FileToOpen.c_str());
+		if (inFile.fail()) {
+			cout << "Unable to open file";
+			exit(1); // terminate with error
+		}else {
+            
+            if (verbose_level > 2){
+                std::cout <<  " Opened <" << FileToOpen << ">"<<std::endl;
+            }
+        }
+		
+		i=0;
+        std::string line;
+		//while (! inFile.eof() && ! inFile.fail()) {
+        while (1==1) {
+            
+            inFile >> x;
+            inFile >> y;
+            
+            if (verbose_level > 2){
+                std::cout << " x: " << x;
+                std::cout << " y: " << y << std::endl;
+            }
+            
+			if (i==0){
+				N_nodes=x;
+				M_edges=y;	
+                break;
+			}
+			i++;
+            
+			
+		}
+		inFile.close();
+        
+        if (verbose_level == 2){
+            std::cout << N_nodes <<  " nodes, " << M_edges << " edges"<<std::endl;
+        }
+		
+		// k_max is the longest distance possible
+		
+        //k_max = M_edges;
+		k_max = 70;
+        
+		//------------------------------------
+		// Get memory allocated for all class members
+		
+		Prob = allocate_2Dmatrix(N_nodes, N_nodes);
+		Tried = allocate_2Dmatrix(N_nodes, N_nodes);
+		Dist = allocate_2Dmatrix(N_nodes, N_nodes);
+		Final = allocate_2Dmatrix(N_nodes, N_nodes);
+        EdgeIndex = allocate_2Dmatrix(N_nodes, N_nodes);
+		Rank = allocate_1Dmatrix(N_nodes);
+		
+        //The second pass through file with the graph
+        
+		for(int i = 0; i < N_nodes; ++i) {
+			Rank[i]=0.;
+			for(int j = 0; j < N_nodes; ++j) {
+				Final[i][j] = 0.;
+				Prob[i][j]=0.;
+				Dist[i][j]=-1.;
+				Tried[i][j]=0.;
+                EdgeIndex[i][j]=-1;
+			}
+		}
+        
+		
+		// Fill in the final graph -- and we are ready to go!
+        
+	    inFile.open(FileToOpen.c_str());
+		if (!inFile) {
+            std::cout << "Unable to open file";
+			exit(1); // terminate with error
+		}
+		else {
+            
+            if (verbose_level > 2){
+                std::cout <<  " Opened <" << FileToOpen << ">"<<std::endl;
+            }
+        }
+        
+		i=0;  
+		while (inFile >> x && inFile >>y) {
+			if (i > 0) {
+				Final[x][y]=1.;
+				Final[y][x]=1.;
+                
+                
+                if (verbose_level == 2){
+                    std::cout << ".";
+                }
+			}
+			i++;
+			
+		}
+        if (verbose_level == 2){
+            std::cout << std::endl;
+        }
+		inFile.close(); 
+		
+        k=0;
+        for (int i=0; i<N_nodes-1; i++){
+            for (int j=i+1;j<N_nodes; j++){
+                if(Final[i][j] > 0.){
+                    EdgeIndex[i][j]=k;
+                    k++;
+                }
+            }
+        }
+        
+        
+		
+		//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+		// create graph -- hopefully, we can keep it, just modifying edge weights
+		
+		
+		edge_array_mine = new Edge[2*M_edges];
+		num_arcs_mine = 2*M_edges;
+		num_nodes_mine = N_nodes;
+		weights_mine = new int[2*M_edges];
+		for (int i=0; i<2*M_edges; i++){ weights_mine[i]=1;}
+		
+		k=0;
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				if (Final[i][j]>0.){
+					edge_array_mine[2*k]  =Edge(i,j);
+					edge_array_mine[2*k+1]=Edge(j,i);
+					k++;
+				}
+			}
+		}
+		graph_t g(edge_array_mine, edge_array_mine + num_arcs_mine, weights_mine, num_nodes_mine);
+		
+		Full_g = g;
+		delete edge_array_mine;
+		delete weights_mine;
+		
+		//===========================================================================
+		std::vector<edge_descriptor> p(num_edges(Full_g));
+		std::vector<int> d(num_edges(Full_g));
+		edge_descriptor s;
+		boost::graph_traits<graph_t>::vertex_descriptor u, v;
+		
+		for (int i=0; i<N_nodes-1; i++){
+			for (int j=i+1; j<N_nodes; j++){
+				if (Final[i][j] > 0.){
+					u = vertex(i, Full_g);
+					v = vertex(j, Full_g);
+					remove_edge(u,v,Full_g);
+					remove_edge(v,u,Full_g);
+					
+				}
+			}
+		}
+		
+        
+    }
+	
+	
+	//=====================================================================
+	int sample_failed_number(double pfail){
+		
+		//boost::geometric_distribution<double> geo(pfail);
+		//boost::variate_generator<base_generator_type&, geometric_distribution<double> > geom(generator, geo);
+		
+		double r, u, g;
+        
+        r=0.;
+		for(int i=0; i<N_steps; i++){
+            
+            u=(double)rand();
+            u = 1.-u /(double)(pow(2.,31)-1.);
+            g=(int)(ceil(log(u) / log(pfail)));
+            
+			//r += geom();
+            
+            r+=g;
+		}
+        
+        if (verbose_level>=3){
+            std::cout << id << " failed " << r << std::endl;
+		}
+		return r;
+		
+	}
+    
+    //=============================================
+    double get_target(void){
+        return TargetNovelty;
+    }
+	
+    //=============================================
+    void set_target(double target){
+        TargetNovelty=target;
+    }
+	
+	//=============================================
+	int sample(){
+		
+        //boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
+        // double r = uni(), Summa = 0.;
+        
+        
+        
+        double r = rand(), Summa = 0.;
+        r /= (double)(pow(2.,31)-1.);
+		int result = 0;
+		int finished = 0;
+        
+        if (verbose_level==4){
+            std::cout << id << " sampled " << r << std::endl;
+        }
+		
+		for(int i=0; i<N_nodes-1 && finished==0; i++){			
+			for( int j=i+1; j<N_nodes && finished==0; j++){
+				
+				Summa += Prob[i][j];
+				
+				if (Summa > r){
+					
+					Tried[i][j]+=1.;
+					
+					if (Final[i][j] > 0.){
+						result = 1;
+					}
+					finished = 1;										
+				}
+			}
+		}
+		
+		return result;
+		
+	}
+	
+	//===============================
+	void update_current_graph(void){
+		
+		std::vector<edge_descriptor> p(num_edges(Full_g));
+		std::vector<int> d(num_edges(Full_g));
+		edge_descriptor s;
+		boost::graph_traits<graph_t>::vertex_descriptor u, v;
+		
+		//property_map<graph_t, edge_weight_t>::type weightmap = get(edge_weight, Full_g);
+		for (int i=0; i<N_nodes-1; i++){
+			for (int j=i+1; j<N_nodes; j++){
+				if (Final[i][j] > 0. && Tried[i][j]>0){
+					//s = edge(i, j, Full_g);	
+					boost::graph_traits<graph_t>::edge_descriptor e1,e2;
+					bool found1, found2;
+					u = vertex(i, Full_g);
+					v = vertex(j, Full_g);
+					tie(e1, found1) = edge(u, v, Full_g);
+					tie(e2, found2) = edge(v, u, Full_g);
+					if (!found1 && !found2){
+						add_edge(u,v,1,Full_g);
+					    add_edge(v,u,1,Full_g);
+					}
+					
+				}
+			}
+			
+		}
+	}
+	
+	//===============================
+	void update_distances(void){
+		// put shortest paths to the *Dist[][]
+		std::vector<vertex_descriptor> p(num_vertices(Full_g));
+		std::vector<int> d(num_vertices(Full_g));
+		vertex_descriptor s;
+		
+		
+		// put shortest paths to the *Dist[][]
+		for (int j=0; j<num_vertices(Full_g); j++){
+			
+			if(Rank[j] > 0.){
+				s = vertex(j, Full_g);	 
+				dijkstra_shortest_paths(Full_g, s, predecessor_map(&p[0]).distance_map(&d[0]));
+				
+				//std::cout <<" Vertex "<< j << std::endl;
+				graph_traits < graph_t >::vertex_iterator vi, vend;
+				
+				for (boost::tie(vi, vend) = vertices(Full_g); vi != vend; ++vi) {
+					
+					if (p[*vi]!=*vi){
+						Dist[*vi][j]=d[*vi];
+						Dist[j][*vi]=d[*vi];
+                        
+                        if (Dist[*vi][j]>max_dist){
+                            max_dist=Dist[*vi][j];
+                        }
+                        
+                        
+					} else {
+						Dist[*vi][j]=-1.;
+						Dist[j][*vi]=-1.;
+					}
+				}
+			}
+			
+		}
+		
+		
+	}
+	
+	//======================================================
+	void update_ranks(void){
+		
+		for(int i=0; i<N_nodes; i++){
+			Rank[i]=0.;
+		}
+		
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				if (Tried[i][j]>0. && Final[i][j] >0.){
+					Rank[i]++;
+					Rank[j]++;
+				}
+			}
+		}
+		
+	}
+	
+	//====================================================================
+	void set_world(double a_i, double a_m, double b, double g, double d){
+		
+		alpha_i=a_i;
+		alpha_m=a_m;
+		gamma=g;
+		beta=b;
+		delta=d;
+		
+	}
+	
+	//====================================================================
+	void reset_world(){
+		
+        //====================================================
+		std::vector<edge_descriptor> p(num_edges(Full_g));
+		std::vector<int> d(num_edges(Full_g));
+		edge_descriptor s;
+		boost::graph_traits<graph_t>::vertex_descriptor u, v;
+        
+		
+		for (int i=0; i<N_nodes-1; i++){
+			for (int j=i+1; j<N_nodes; j++){
+				if (Final[i][j] > 0. && Tried[i][j] > 0){
+					u = vertex(i, Full_g);
+					v = vertex(j, Full_g);
+					remove_edge(u,v,Full_g);
+					remove_edge(v,u,Full_g);
+					
+				}
+			}
+		}
+        
+        //==================================================
+        
+		current_loss=0;
+		current_epoch=0;
+		current_repeat++;
+        current_novelty=0;
+		
+		for(int i = 0; i < N_nodes; ++i) {
+			Rank[i]=0.;
+			for(int j = 0; j < N_nodes; ++j) {
+				Prob[i][j]=0.;
+				Dist[i][j]=-1.;
+				Tried[i][j]=0.;
+			}
+		}
+	}
+	
+	
+    //==============================================
+    void show_parameters(void){
+        
+        std::cout << "Parameters: " 
+        << alpha_i << " "
+        << alpha_m << " | "
+        << beta << " "
+        << gamma << " | "
+        << delta << std::endl;
+        
+    }
+    
+    
+    
+    //===============================================
+    string file_name(){
+        
+        std::ostringstream s;
+        s << "world_" 
+        << lexical_cast<std::string>(alpha_i) << "_" 
+        << lexical_cast<std::string>(alpha_m) << "_"
+        << lexical_cast<std::string>(beta) << "_"
+        << lexical_cast<std::string>(gamma) << "_"
+        << lexical_cast<std::string>(delta) << "_"
+        << lexical_cast<std::string>(N_epochs) << "_"
+        << lexical_cast<std::string>(N_steps) << "_"
+        << lexical_cast<std::string>(N_repeats) << ".txt";
+        
+        return s.str();
+        
+    }
+    
+    
+    
+    
+    //=================================================
+    void set_verbose(int verbose){
+        
+        verbose_level = verbose;
+    }
+    
+    
+    //=============================================================
+    void update_probabilities(void){
+        
+        
+        //=========================
+		// Compute sampling probabilities
+		// first pass: \xi_i,j
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				
+				double bg = 0.;
+				
+				Prob[i][j] = alpha_i*log(min(Rank[i]+1.,Rank[j]+1.)) + 
+                alpha_m*log(max(Rank[i]+1.,Rank[j]+1.));
+				
+                if (Dist[i][j] > 0.){
+                    
+                    double k = Dist[i][j];
+                    if (k >= k_max){
+                        k = k_max-1;
+                    }
+					
+                    bg = beta * log(k/k_max) + gamma * log(1. - k/k_max);
+					
+                } else {
+                    bg = delta;
+                }
+				
+				Prob[i][j] = exp(Prob[i][j] + bg);
+			}
+		}
+        
+		
+		// second pass: sum
+		double Summa = 0.;
+		
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				Summa += Prob[i][j];
+			}
+		}
+		
+		// third pass: normalize
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				Prob[i][j] /= Summa;
+			}
+		}
+        
+    }
+    
+	// Now we are ready for simulations
+	//==============================================
+	void update_world(){
+		
+		int failed = 0;
+        
+		// Given current universe compute shortest paths
+		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+		
+		update_current_graph();
+		update_ranks();				
+		update_distances();
+		update_probabilities();
+		
+		//===============================
+		// sampling
+		int result;
+		double cost=0., novel=0.;
+		int publishable = 0;
+		
+		
+		//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+		if (mode_identify_failed == 1){
+			
+			while(publishable < N_steps){
+				
+		    	result = sample();
+			    publishable += result;
+			    failed += (1-result);
+				
+			}
+			
+			for(int i=0; i<N_nodes-1; i++){
+				for( int j=i+1; j<N_nodes; j++){
+					
+					cost+=Tried[i][j];
+					
+					if (Tried[i][j]>0. && Final[i][j]>0.){
+						novel+=1.;
+					}
+				}
+			}
+			
+		}
+		//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+		else {
+			
+			double pfail=0.;
+			int n_failed;
+			//, n_check = 0;
+			
+			for(int i=0; i<N_nodes-1; i++){
+				for( int j=i+1; j<N_nodes; j++){
+					if (Final[i][j] == 0.){
+						pfail += Prob[i][j];
+						Prob[i][j] = 0.;
+					}
+					
+				}
+			}
+			
+			for(int i=0; i<N_nodes-1; i++){
+				for( int j=i+1; j<N_nodes; j++){
+					Prob[i][j] /= (1.-pfail);
+				}
+				//std::cout << std::endl;
+			}			
+			
+			n_failed = sample_failed_number(pfail);
+			while(publishable < N_steps){
+				
+		    	result = sample();
+			    publishable += result;					
+			}
+            
+            
+			current_loss += (n_failed + N_steps);
+			cost = current_loss;
+			
+			for(int i=0; i<N_nodes-1; i++){
+				for( int j=i+1; j<N_nodes; j++){
+					
+					if (Tried[i][j]>0. && Final[i][j]>0.){
+						novel+=1.;
+					}
+				}
+			}
+		}
+        
+        current_novelty = novel;
+        
+        
+		//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+		if (verbose_level == 2){
+            std::cout << (current_repeat+1) << "  epoch=" << (current_epoch+1) 
+            
+		    << "  cost=" << cost 
+		    << " novel=" << novel 
+		    << " rel_loss=" << cost/novel
+		    << std::endl;
+        }
+        
+		current_epoch++;
+	}
+	
+	
+	//======  Destructor ======
+	~Universe(){
+		
+		delete_2Dmatrix(Final, N_nodes);
+		delete_2Dmatrix(Dist, N_nodes);
+		delete_2Dmatrix(Tried, N_nodes);
+		delete_2Dmatrix(Prob, N_nodes);
+        delete_2Dmatrix(EdgeIndex, N_nodes);
+		delete_1Dmatrix(Rank);		
+	}
+    
+    //================================================
+    // Allocate memory
+    double** allocate_2Dmatrix(int N, int M)
+    {
+        double **pointer;
+        
+        if (verbose_level == 2){
+            std::cout<< "["<<N<<"|"<<M<<"]"<<std::endl;
+        }
+        pointer = new double*[N];
+        for (int i = 0; i < N; ++i)
+            pointer[i] = new double[M];	
+        
+        return pointer;
+    }
+    //===================
+    double* allocate_1Dmatrix(int N)
+    {
+        double *pointer;
+        
+        if(N > 0){
+            
+            pointer = new double[N];
+            
+        }else {
+            
+            pointer = NULL;
+        }
+        
+        return pointer;
+        
+    }
+    
+    //==============================================
+    // De-Allocate memory to prevent memory leak
+    void delete_2Dmatrix(double **pointer, int N){
+        
+        if (pointer != NULL){
+            
+            for (int i = 0; i < N; ++i){
+                delete [] pointer[i];
+            }
+            delete [] pointer;
+        }
+    }
+    //====================
+    void delete_1Dmatrix(double *pointer){
+        
+        delete [] pointer;
+    }
+    
+    //===========================================
+    double get_rel_loss(){
+        
+        return CumulativeRelativeLoss ;
+    }
+    
+    //===========================================
+    double get_rel_loss_err(){
+        
+        return CRLsquare ;
+    }
+    
+    
+    
+    //==================================================================================
+    void evolve_to_target_and_save(int istart, int iend, double* storage, int* counters){
+        
+        double ALOT=100000000000.;
+        
+	// std::cout<<" evolve_to_target_and_save: istart=" << istart << "iend=" << iend << "\n";
+
+        reset_world();
+        
+        for (int k = istart; k < iend; k++){
+
+	  // std::cout<<" evolve: k=" << k << "\n";
+
+            
+            
+            for(int i=0; i< N_epochs &&  current_novelty < TargetNovelty; i++){
+	      // std::cout<<" evolve: k=" << k << " i=" << i << " cur=" << current_novelty << " Target=" << TargetNovelty << "\n";
+                update_world();
+            }
+            
+            storage[k]=current_loss/current_novelty;
+            counters[k]=1;
+            
+            
+            reset_world();
+        }
+        
+    }
+    //==============================================
+    int get_reruns(void){
+        return N_repeats;
+    }
+    
+    //==============================================
+    double get_parameter(int i){
+        
+        switch(i){
+            case 0:
+                return alpha_i;
+            case 1:
+                return alpha_m;
+            case 2:
+                return beta;
+            case 3:
+                return gamma;
+            case 4:
+                return delta;
+            default:
+                
+                std::cout << "Erroneous parameter id!!!!\n\n\n";
+                return 0.;
+        }
+    }
+    
+    
+    //==============================================
+    void evolve_to_target(){
+        
+        reset_world();
+        if (beta < -1. || gamma < -1.){
+            CumulativeRelativeLoss = 100000000000.;
+            CRLsquare = 0.;
+            return;
+        }
+        
+        
+        for (int k=0; k< N_repeats; k++){
+            
+            
+            for(int i=0; i<N_epochs &&  current_novelty < TargetNovelty; i++){
+                update_world();
+            }
+            
+            CumulativeRelativeLoss += current_loss/current_novelty;
+            CRLsquare += (current_loss/current_novelty)*(current_loss/current_novelty);
+            if(verbose_level==3){
+                std::cout <<  CumulativeRelativeLoss << " | " << CRLsquare << std::endl;
+            }
+            
+            if(verbose_level==1){
+                std::cout <<  "." ;
+                
+            }
+            else if(verbose_level==2){
+                std::cout <<  "**" << (k+1) <<  "**  curr loss " << current_loss << "; curr novelty " << current_novelty << std::endl;
+            }
+            
+            
+            reset_world();
+        }
+        
+        CumulativeRelativeLoss /= double(N_repeats);
+        CRLsquare /= double(N_repeats);
+        
+        if(verbose_level==1){
+            std::cout << std::endl;
+        }
+        
+        if(verbose_level==2){
+            std::cout <<  CumulativeRelativeLoss << " || " << CRLsquare << std::endl;
+        }
+        
+        CRLsquare = 2*sqrt((CRLsquare - CumulativeRelativeLoss*CumulativeRelativeLoss)/double(N_repeats));
+        
+    }
+    
+    
+    //================================================================    
+    int set_parameter(double value, int position){
+        
+        if (position < 0 || position > 4) {return 0;}
+        
+        else {
+            
+            switch(position){
+                case 0:
+                    alpha_i=value;
+                    return 1;
+                case 1:
+                    alpha_m=value;
+                    return 1;
+                case 2:
+                    beta=value;
+                    return 1;
+                case 3:
+                    gamma=value;
+                    return 1;
+                case 4:
+                    delta=value;
+                    return 1;
+            }
+            
+        }
+        
+        return 0;
+    }
+
+#ifdef notdef
+    //=================================================================
+    void try_annealing(double starting_jump, int iterations, 
+                       double temp_start, double temp_end, double target_rejection){
+        
+        double dx[5]={0.,0.,0.,0.,0};
+        double x[5]={0.,0.,0.,0.,0};
+        double rejection[5]={0., 0., 0., 0., 0.};
+        double curr_x, curr_err, x_tmp;
+        double temperature;
+        double ratio, r;
+        int cycle=10;
+        boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
+        
+        // set up parameter for annealing
+        
+        x[0]=alpha_i;
+        x[1]=alpha_m;
+        x[2]=beta;
+        x[3]=gamma;
+        x[4]=delta;
+        
+        for(int i=0;i<5;i++){
+            dx[i] = starting_jump;
+        }
+        
+        // establish the current value
+        
+        //..........................................
+        evolve_to_target();        
+        std::cout << CumulativeRelativeLoss << " +- " << CRLsquare << std::endl;
+        
+        curr_x   = CumulativeRelativeLoss;
+        curr_err = CRLsquare;
+        CumulativeRelativeLoss = 0;
+        CRLsquare = 0;
+        //...........................................
+        
+        // optimization cycle
+        for(int i=0; i<iterations; i++){
+            
+            temperature = temp_start*exp( i*(log(temp_end)-log(temp_start))/(double)iterations);
+            std::cout  << std::endl << "....T = " << wrap_double(temperature,3) << std::endl << std::endl;
+            
+            if (i % cycle == 0 && i > 0){
+                
+                for (int k=0; k<5; k++){
+                    
+                    rejection[k]/=(double)cycle;
+                    if (rejection[k] > 0){
+                        dx[k] = dx[k]/(rejection[k]/target_rejection);
+                        rejection[k]=0.;
+                    }
+                    else{
+                        dx[k]*=2.;
+                    }
+                    std::cout  << dx[k] << " ";
+                }
+                std::cout  << std::endl;
+            }
+            
+            
+            for (int j=0; j<5; j++){
+                
+                // get new value of x[j]
+                x_tmp = get_new_x(x[j],dx[j]);
+                
+                
+                
+                //.............................................
+                set_parameter(x_tmp, j);
+                
+                
+                evolve_to_target(); 
+                
+                std::cout  << std::endl << "......... " << std::endl;
+                std::cout << "Trying... " << CumulativeRelativeLoss << " +- " << CRLsquare << std::endl;
+                
+                ratio = min(1.,exp(-(CumulativeRelativeLoss-curr_x)/temperature));
+                r = uni();
+                std::cout << r << " vs " << ratio << std::endl;
+                
+                if (r > ratio){
+                    
+                    std::cout << string_wrap(id, 4) <<" "<< (i+1) << ","<< (j) 
+                    <<" "<< (i+1) << " Did not accept " 
+                    << x_tmp << "(" << j << ")" << std::endl;
+                    std::cout << alpha_i << " "<< alpha_m << " "
+                    << beta << " " << gamma << " " 
+                    << delta << " " << std::endl;
+                    set_parameter(x[j], j);
+                    CumulativeRelativeLoss = 0;
+                    CRLsquare = 0;
+                    
+                    rejection[j]+=1.;
+                }
+                
+                else {
+                    
+                    curr_x   = CumulativeRelativeLoss;
+                    curr_err = CRLsquare;
+                    x[j] = x_tmp;
+                    CumulativeRelativeLoss = 0;
+                    CRLsquare = 0;
+                    std::cout << (i+1) << string_wrap((string) " Rejection counts: ", 8) 
+                    << wrap_double(rejection[0],2) 
+                    << " "<< wrap_double(rejection[1], 7) << " "
+                    << wrap_double(rejection[2],5) << " " << wrap_double(rejection[2],9) << " " 
+                    << wrap_double(rejection[4],6) << " " 
+                    << std::endl << std::endl;
+                    
+                    std::cout << string_wrap(id, 4) <<" "<< (i+1) <<","<< (j) 
+                    <<" "
+                    << string_wrap((string) "***** Did accept! ", 3) 
+                    << wrap_double(alpha_i,2) 
+                    << " "<< wrap_double(alpha_m, 7) << " "
+                    << wrap_double(beta,5) << " " 
+                    << wrap_double(gamma,9) << " " 
+                    << wrap_double(delta,6) << " " 
+                    << std::endl << std::endl;
+                    
+                }
+                //........................................................ 
+                
+            }
+            
+        }
+        
+    }
+    
+#endif // notdef
+
+	
+};
+
+
+//============================================================
+
+std::pair<double,double> multi_loss( // dispatch_group_t group, 
+                                    Universe* un[], 
+                                    // dispatch_queue_t* CustomQueues,
+                                    double* Results,
+                                    int*    Counters,
+                                    double* params){
+    
+    int N = un[0]->get_reruns();
+    int step = (int)(double)N/(double)(Nworkers);
+    int istart=0;
+    int iend = istart+step;
+    
+    double Loss=0., LossSquare=0.;
+    
+    timeval startTime, endTime;
+    double elapsedTime;
+    // start timer
+    gettimeofday(&startTime, NULL);
+
+    //err:    for(int i=0; i<Nworkers; i++){
+
+    std::cout << "Entry to multi_loss, params: ";
+    for(int i=0; i<Nworkers; i++){
+        for(int j=0; j<NEVOPARAMS; j++){
+	  std::cout << "[" << un[i]->get_parameter(j) << "," << params[j] << "] ";
+            un[i]->set_parameter(params[j],j);
+        }
+    }
+    std::cout << "\n";
+    int i;
+    #pragma omp parallel for private (i)
+    for(i=0; i<Nworkers; i++){
+
+      // dispatch_group_async(group, CustomQueues[i], ^{
+      std::cout<<"multi_loss: Calling evolve_to_target_and_save i=" << i << " N=" << N << " step=" << step << " istart=" << i*step << " iend=" << (i+1)*step << "\n";
+      //un[i]->evolve_to_target_and_save(istart, iend, Results, Counters);
+      un[i]->evolve_to_target_and_save(i*step, min((i+1)*step,N), Results, Counters);
+      //});
+
+      // std::cout<<"multi_loss: Returned from evolve_to_target_and_save " << i << "\n";
+      
+
+      //istart += step;
+      //iend = min(istart+step,N);
+            
+    }
+    // err    }
+    // dispatch_group_wait(group, DISPATCH_TIME_FOREVER);
+    //dispatch_release(group);
+    
+
+    for (int i=0; i<N; i++){
+        
+        Loss+=Results[i]/(double)N;
+        LossSquare+=Results[i]*Results[i]/(double)N;
+        
+        std::cout<<i<<":"<< Results[i] << " ";
+    }
+    
+    std::cout<<" \n\n\n";
+    double two_std = ((LossSquare - Loss*Loss)/(double)N);
+    
+    two_std = 2.*sqrt(two_std);
+    std::pair<double,double> Res;
+    Res.first=Loss;
+    Res.second=two_std;
+
+    gettimeofday(&endTime, NULL);
+    elapsedTime = (endTime.tv_sec - startTime.tv_sec) * 1000.0;      // sec to ms
+    elapsedTime += (endTime.tv_usec - startTime.tv_usec) / 1000.0;   // us to ms
+    elapsedTime /= 1000.;
+    cout << "multi_loss(N=" << N << ") elapsed time: " << elapsedTime << " seconds " << elapsedTime/60. << " minutes\n\n";
+    
+    return Res;
+}
+//============================================================
+
+
+//============================================================
+void multi_annealing( // dispatch_group_t group, 
+                     Universe* un[], 
+                     // dispatch_queue_t* CustomQueues, 
+                     double T_start, double T_end, 
+                     double Target_rejection, 
+                     int Annealing_repeats, 
+                     double starting_jump,
+                     double* Results,
+                     int*    Counters,
+                     double* params0,
+                     double annealing_cycles){
+    //.................................
+    // re-implement annealing
+    
+    double dx[NEVOPARAMS]={0.,0.,0.,0.,0};
+    double x[NEVOPARAMS]={0.,0.,0.,0.,0};
+    double rejection[NEVOPARAMS]={0., 0., 0., 0., 0.};
+    double curr_x, curr_err, x_tmp;
+    double temperature;
+    double ratio, r;
+    int cycle=10;
+    //boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
+    
+    // set up parameter for annealing
+
+    for(int i=0;i<NEVOPARAMS;i++){
+        x[i]=params0[i];
+        dx[i] = starting_jump;
+	for(int w=0; w<Nworkers; w++){
+	  un[w]->set_parameter(x[i], i);
+	}
+    }
+    
+    // establish the current value
+    std::pair<double,double>Res;
+    
+    if ( operation == 'm' ) {
+      // Nworkers = 1;
+    }
+    else if (operation == 'g') {
+      // generate params
+    }
+    else if (operation == 'a') {
+      // analyze multi_loss() results
+
+      string line;
+      ifstream mlossdata ("multi_loss.data");
+      double d, Loss, LossSquare, two_std;
+      bool b;
+      int n=0;
+      if (mlossdata.is_open()) {
+	while ( getline (mlossdata,line) ) {
+	  b = from_string<double>(d, std::string(line), std::dec);
+	  cout << line << " d=" << d << endl;
+	  Loss += d;
+          LossSquare += (d*d);
+          n++;
+	}
+        Loss /= double(n);
+	LossSquare /= double(n);
+	two_std = ((LossSquare - Loss*Loss)/(double)n);
+	two_std = 2.*sqrt(two_std);
+	std::cout<<"n="<<n<<" Loss="<<Loss<<" LossSquare="<<LossSquare<<" two_std="<<two_std<<"\n\n\n";
+	mlossdata.close();
+	FILE *f=fopen("multi_loss_stats.txt","w");
+	fprintf(f,"%d n\n",n);
+	fprintf(f,"%.20e Loss\n",Loss);
+	fprintf(f,"%.20e LossSquare\n",LossSquare);
+	fprintf(f,"%.20e two_std\n",two_std);
+        fclose(f);
+	exit(0);
+      }
+      else {
+	cout << "Unable to open file multi_loss.data"; 
+	exit(1);
+      }
+    }
+
+    std::cout << "Calling initial multi_loss:\n";
+    Res = multi_loss( /* group,*/ un, /* CustomQueues,*/ Results, Counters, x);
+    std::cout << "Ret from initial multi_loss:\n";
+    std::cout << Res.first << " +- " << Res.second << std::endl;
+
+    if ( operation == 'm' ) {
+      FILE *f;
+      int N = un[0]->get_reruns();
+
+      f = fopen("multi_loss.data","w");
+      for(int i=0; i<N; i++) {
+	fprintf(f,"%.20e\n",Results[i]);
+      }
+      fclose(f);
+      exit(0);
+    }
+    
+    curr_x   = Res.first;
+    curr_err = Res.second;
+    
+    // optimization cycle
+    
+    for(int i=0; i<annealing_cycles; i++){
+        
+        temperature = T_start*exp( i*(log(T_end)-log(T_start))/(double)annealing_cycles);
+        std::cout  << std::endl << "....T = " << wrap_double(temperature,3) << std::endl << std::endl;
+        
+        if (i % cycle == 0 && i > 0){
+            
+            for (int k=0; k<NEVOPARAMS; k++){
+                rejection[k]/=(double)cycle;
+                
+                if (rejection[k] > 0){
+                    dx[k] = dx[k]/(rejection[k]/Target_rejection);
+                    rejection[k]=0.;
+                }
+                else{
+                    dx[k]*=2.;
+                }
+                std::cout  << dx[k] << " ";
+            }
+            std::cout  << std::endl;
+        }
+
+        for (int j=0; j<NEVOPARAMS; j++){
+            
+            ///////////////////////////////
+            if (FIX_VARIABLES==0 || var_fixed[j]==0){
+                
+                
+                
+                // get new value of x[j]
+                double x_hold=x[j];
+                x_tmp = get_new_x(x[j],dx[j]);
+                x[j]=x_tmp;
+                
+                std::cout << wrap_double(x_tmp,10) << " " << wrap_double(j,9) << "\n\n"; 
+                //=======================================
+                //.............................................
+                for(int w=0; w<Nworkers; w++){
+                    un[w]->set_parameter(x_tmp, j);
+                }
+                
+
+		// WRITE OUT PARAMS HERE; then exit.
+
+		std::cout << "Calling multi_loss: i=" << i << " j=" << j << "\n";
+                Res = multi_loss(/* group, */ un, /* CustomQueues, */ Results, Counters, x);
+		std::cout << "Ret from multi_loss: i=" << i << " j=" << j << "\n";
+                std::cout << Res.first << " +- " << Res.second << std::endl;
+                
+                ratio = min(1.,exp(-(Res.first-curr_x)/temperature));
+                r = rand()/(double)(pow(2.,31)-1.);
+                std::cout << r << " vs " << ratio << std::endl;
+                
+                double ALOT=100000000000.;
+                
+                if (Res.first < ALOT)
+                {
+                    ofstream filestr;
+                    
+                    filestr.open ("best_opt_some.txt", ofstream::app);
+                    
+                    // >> i/o operations here <<
+                    filestr << un[0]->get_target() << "," 
+                    << Res.first 
+                    << "," << un[0]->get_parameter(0) 
+                    << "," << un[0]->get_parameter(1) 
+                    << "," << un[0]->get_parameter(2) 
+                    << "," << un[0]->get_parameter(3) 
+                    << "," << un[0]->get_parameter(4) << "," << Res.second << ",\n";
+                    
+                    filestr.close();
+                    
+                    
+                    filestr.open ("max_dist.txt", ofstream::app);
+                    
+                    // >> i/o operations here <<
+                    filestr << max_dist << ",\n";
+                    
+                    filestr.close();
+                    
+                }
+                
+                
+                if (r > ratio){
+                    
+                    std::cout << " "<< (i+1) << ","<< (j) 
+                    <<" "<< (i+1) << " Did not accept " 
+                    << x_tmp << "(" << j << ")" << std::endl;
+                    std::cout << un[0]->get_parameter(0) 
+                    << " " << un[0]->get_parameter(1) 
+                    << " " << un[0]->get_parameter(2) 
+                    << " " << un[0]->get_parameter(3) 
+                    << " " << un[0]->get_parameter(4) << " " << std::endl;
+                    
+                    x[j]=x_hold;
+                    for(int w=0; w<Nworkers; w++){
+                        un[w]->set_parameter(x[j], j);
+                    }
+                    
+                    
+                    //set_parameter(x[j], j);     
+                    rejection[j]+=1.;
+                }
+                
+                else {
+                    
+                    curr_x   = Res.first;
+                    curr_err = Res.second;
+                    x[j] = x_tmp;
+                    
+                    for(int w=0; w<Nworkers; w++){
+                        un[w]->set_parameter(x[j], j);
+                    }
+                    
+                    std::cout << (i+1) << string_wrap((string) " Rejection counts: ", 8) 
+                    << wrap_double(rejection[0],2) << " " 
+                    << wrap_double(rejection[1],7) << " "
+                    << wrap_double(rejection[2],5) << " " 
+                    << wrap_double(rejection[3],9) << " " 
+                    << wrap_double(rejection[4],6) << " " 
+                    << std::endl << std::endl;
+                    
+                    std::cout << " "<< (i+1) <<","<< (j) 
+                    <<" "
+                    << string_wrap((string) "***** Did accept! ", 3) 
+                    << wrap_double(un[0]->get_parameter(0),2) << " "
+                    << wrap_double(un[0]->get_parameter(1),7) << " "
+                    << wrap_double(un[0]->get_parameter(2),5) << " " 
+                    << wrap_double(un[0]->get_parameter(3),9) << " " 
+                    << wrap_double(un[0]->get_parameter(4),6) << " " 
+                    << std::endl << std::endl;
+                    
+                    
+                    
+                }
+                //........................................................ 
+                
+            }
+        }
+        
+    }
+    
+}
+
+
+//================================================
+int
+main(int argc, char* argv[])
+{
+    
+    double params0[6] = {0., 0., 0., 0., 0., 0.2}, target=50., range;
+    string par_names0[6] = {"alpha_i", "alpha_m", "beta", "gamma", "delta", "target"};
+    string par_names1[4] = {"n_epochs", "n_steps", "n_reruns", "range"};
+    string par_names2[5] = {"T_start", "T_end", "Annealing_steps","Target_rejection","Starting_jump"};
+    string par_names3[5] = {"FREEZE_alpha_i", "FREEZE_alpha_m", "FREEZE_beta", "FREEZE_gamma", "FREEZE_delta"};    
+    string par_names4[2] = {"Operation", "Nworkers"};    
+    int params1[4] = {300, 50, 1000, 10};
+    int params3[5] = { 0, 0, 0, 0, 0};
+    
+    //          temperature_start,  temperature_end,  annealing_steps target_rejection  Starting_jump
+    double params2[5] = {1,             0.001,               100,              0.3,           1.5};
+    
+    int verbose_level = 2;
+    const std::string one="one", two="two";
+    static Universe* un[MAXNworkers];
+    // static dispatch_queue_t CustomQueues[MAXNworkers];
+    
+    static double* Results;
+    static int*    Counters;
+    
+    timeval t1, t2;
+    double elapsedTime;
+    // start timer
+    gettimeofday(&t1, NULL);
+    
+    
+    if (argc < 8) {
+        std::cout << "Usage: super_optimizer alpha_i alpha_m beta gamma delta target_innov [n_epochs n_steps n_reruns] [range] [verbose_level]\n"; 
+        std::cout << "         [T_start T_end Annealing_steps Target_rejection Starting_jump]\n"; 
+        std::cout << "         [FREEZE_alpha_i FREEZE_alpha_m FREEZE_beta FREEZE_gamma FREEZE_delta]\n"; 
+        
+        system("pwd");
+        
+        
+        return(1);
+    }
+    else {
+      std::cout << "argc=" << argc << std::endl;
+
+        for (int nArg=0; nArg < argc; nArg++){
+            //std::cout << nArg << " " << argv[nArg] << std::endl;
+            if (nArg > 0 && nArg < 7){
+                params0[nArg-1]= atof(argv[nArg]);
+                std::cout << par_names0[nArg-1] << ": " << params0[nArg-1] <<  std::endl;
+            }
+            if (nArg > 6 && nArg < 11){
+                params1[nArg-7]= atoi(argv[nArg]);
+                std::cout << par_names1[nArg-7] << ": " << params1[nArg-7] <<  std::endl;
+            }
+            if (nArg == 11){
+                verbose_level = atoi(argv[nArg]);
+                std::cout << "verbose level: " << verbose_level <<  std::endl;
+            }
+            if (nArg > 11 && nArg < 17){
+                params2[nArg-12]= atof(argv[nArg]);
+                std::cout << par_names2[nArg-12] << ": " << params2[nArg-12] <<  std::endl;
+            }
+            if (nArg > 16 && nArg < 22){
+                params3[nArg-17]= atof(argv[nArg]);
+                var_fixed[nArg-17]= atof(argv[nArg]);
+                std::cout << par_names3[nArg-17] << ": " << var_fixed[nArg-17] <<  std::endl;
+            }
+            if (nArg == 22 ){
+                operation = *argv[nArg];
+		std::cout << par_names4[0] << ": " << operation <<  std::endl;
+            }
+            if (nArg == 23 ){
+	        Nworkers = atoi(argv[nArg]);
+		std::cout << par_names4[1] << ": " << Nworkers <<  std::endl;
+            }
+        }
+    }
+    
+    /*
+    for target in range(58,1009,50):
+	s = ("%d" % target)
+	print s
+
+	for i in range(15):
+          # Param groups separated by "|" below for documentation. NOTE that | is not used on command line!
+	  os.system("./supe_duper_optimizer |0 0 4 50 -1 "+s+" | 40000 20 1000 2 | 1 | 2. 0.01 100 0.3 2.3 | 1 1 0 0 0")
+    */
+
+    /*   Parameters, re-iterated:
+
+                        "alpha_i", "alpha_m", "beta", "gamma", "delta", "target"};
+    double params0[6] = {0.,        0.,        0.,     0.,      0.,      0.2}, target=50., range;
+
+                    {"n_epochs", "n_steps", "n_reruns", "range"};
+    int params1[4] = {300,        50,       1000,        10};
+
+                       {"T_start",          "T_end",          "Annealing_steps","Target_rejection","Starting_jump"};
+                         temperature_start,  temperature_end,  annealing_steps target_rejection  Starting_jump
+    double params2[5] = {1,                  0.001,            100,            0.3,              1.5};
+
+                     {"FREEZE_alpha_i", "FREEZE_alpha_m", "FREEZE_beta", "FREEZE_gamma", "FREEZE_delta"};    
+    int params3[5] = { 0,                0,                0,             0,              0};
+
+
+    */
+
+    for (int j=0; j<NEVOPARAMS; j++){
+        
+        cout << j << " | " << var_fixed[j] << " (fixed) \n";
+    }
+	
+    target=params0[NEVOPARAMS];
+    range = (double)params1[3];
+	int identify_failed = 0;
+	char* filename= (char *)"movie_graph.txt";
+ 	int n_ep=params1[0], n_st=params1[1], n_rep=params1[2];
+    
+    //...............................
+    
+    for(int i=0; i<Nworkers; i++){
+        un[i] = new Universe((char *)filename,n_ep,n_st,
+                             (int)n_rep,
+                             identify_failed, target, i2string(i));
+        // CustomQueues[i] = dispatch_queue_create(i2char(i), NULL);
+    }
+    
+    //...............................
+    // mw: n_rep == n_reruns = # times evolve is done within multi_loss, spread across Nworkers.
+    if(n_rep > 0){
+        
+        Results = new double[n_rep];
+        Counters = new int[n_rep];
+        
+    }else {
+        
+        Results =  NULL;
+        Counters = NULL;
+        std::cout << " Number of reruns should be positive! " <<  std::endl;
+        return 0;
+        
+    }
+    //...............................
+
+    //srand(time(0));
+    //srandomdev();
+    {
+      timeval t; 
+      gettimeofday(&t, NULL);
+      srand(t.tv_usec);
+    }
+    
+    {
+        double r=0;
+        for (int j=0; j<100; j++){
+            
+            
+            
+            r = rand()/(double)(pow(2.,31)-1.);
+            std::cout << r << " ";
+        }
+        std::cout << "\n ";
+    }
+  	//random initiation of starting parameters
+    
+    if (range > 0.){
+        
+        for (int i=0; i < 5; i++){
+            
+            if (params0[i]==-100.){
+                
+                double r1 = (rand()/(double)(pow(2.,31)-1.));
+                double r2 = (rand()/(double)(pow(2.,31)-1.));
+                double sign = 1.;
+                
+                if(r1 > 0.5){
+                    sign=-1.;
+                }
+                
+                params0[i] = sign*r2*range;
+                
+                std::cout << par_names0[i] << ": " << params0[i] <<  std::endl;
+            }
+        }
+        
+    }
+    
+    
+    double T_start=params2[0], T_end=params2[1], Target_rejection=params2[3], starting_jump=params2[4];
+    int Annealing_repeats = (int) params2[2];
+    
+    
+    // dispatch_group_t group = dispatch_group_create();
+    
+    //.............................
+    multi_annealing( /* group, */ un, /* CustomQueues, */ T_start, T_end, Target_rejection, Annealing_repeats, 
+		     starting_jump, Results, Counters, params0, Annealing_repeats);
+    
+    //dispatch_group_wait(group, DISPATCH_TIME_FOREVER);
+    // dispatch_release(group);
+    //.............................
+    
+    
+    // stop timer
+    gettimeofday(&t2, NULL);
+    
+    // compute and print the elapsed time in millisec
+    elapsedTime = (t2.tv_sec - t1.tv_sec) * 1000.0;      // sec to ms
+    elapsedTime += (t2.tv_usec - t1.tv_usec) / 1000.0;   // us to ms
+    elapsedTime /= 1000.;
+    cout << elapsedTime << " seconds " << elapsedTime/60. << " minutes\n\n";
+    
+    //.....................
+    
+    for(int i=0; i<Nworkers; i++){
+        delete un[i];
+    }
+    
+    //....................
+    if(n_rep > 0){
+        
+        delete [] Results;
+        delete [] Counters;
+        
+    }
+    
+    return 0;
+	
+	
+	
+}
+


Property changes on: SwiftApps/optimizer.cpp
___________________________________________________________________
Added: svn:executable
   + 

Added: SwiftApps/optimizer.orig.cpp
===================================================================
--- SwiftApps/optimizer.orig.cpp	                        (rev 0)
+++ SwiftApps/optimizer.orig.cpp	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,1710 @@
+//
+//  main.cpp
+//  optimizer
+//
+//  Created by Andrey Rzhetsky on 4/11/11.
+//  Copyright 2011 University of Chicago. All rights reserved.
+//
+
+#define Nworkers 1 // 24
+
+#include <fstream>
+#include <iostream>
+#include <stdio.h>
+#include <time.h>
+#include <ctime>    
+#include <algorithm>
+#include <string>
+
+#include <stdio.h>
+#include <sys/param.h>
+#include <sys/time.h>
+#include <sys/types.h>
+
+// #include <dispatch/dispatch.h>
+#include <fstream>
+
+
+#include <stdlib.h>
+#include <boost/numeric/ublas/io.hpp>
+#include <boost/graph/graph_traits.hpp>
+#include <boost/graph/dijkstra_shortest_paths.hpp>
+#include <boost/graph/loop_erased_random_walk.hpp>
+#include <boost/graph/random.hpp>
+#include <boost/property_map/property_map.hpp>
+#include <boost/graph/graph_concepts.hpp>
+#include <boost/graph/properties.hpp>
+
+#include <boost/graph/graph_traits.hpp>
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/adjacency_matrix.hpp>
+
+#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
+#define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
+#include <boost/graph/random.hpp>
+#include <boost/random/geometric_distribution.hpp>
+#include <boost/random/uniform_01.hpp>
+#include <boost/random.hpp>
+#include <boost/random/linear_congruential.hpp>
+#include <boost/random/uniform_int.hpp>
+#include <boost/random/uniform_real.hpp>
+#include <boost/random/variate_generator.hpp>
+#include <boost/generator_iterator.hpp>
+#include <boost/lexical_cast.hpp>
+
+#define INT_INFINITY 2147483647
+
+#define FIX_VARIABLES 1
+
+using namespace boost;
+using namespace std;
+using namespace boost::numeric::ublas;
+
+static int max_dist=0;
+
+typedef boost::adjacency_matrix<boost::directedS> Graph;
+typedef std::pair<int,int> Edge;
+typedef boost::graph_traits<Graph> GraphTraits;
+typedef boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::strict_upper> prob;
+typedef boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::strict_upper> pathlength;
+typedef boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
+
+namespace std {
+	using ::time;
+}
+
+static int var_fixed[5] = {1, 0, 1, 1, 0};
+
+typedef boost::minstd_rand base_generator_type;
+typedef adjacency_list < listS, vecS, directedS,
+no_property, property < edge_weight_t, int > > graph_t;
+typedef graph_traits < graph_t >::vertex_descriptor vertex_descriptor;
+typedef graph_traits < graph_t >::edge_descriptor edge_descriptor;
+
+
+//================================================
+string strDouble(double number)
+{
+    stringstream ss;//create a stringstream
+    ss << number;//add number to the stream
+    return ss.str();//return a string with the contents of the stream
+}
+
+//================================================
+
+double gaussian(double sigma)
+{
+    double GaussNum = 0.0;
+    int NumInSum = 10;
+    for(int i = 0; i < NumInSum; i++)
+    {
+        GaussNum += ((double)rand()/(double)RAND_MAX - 0.5);
+    }
+    GaussNum = GaussNum*sqrt((double)12/(double)NumInSum);
+    
+    
+    return GaussNum;
+    
+}
+
+
+
+//=================================================
+double diffclock(clock_t clock1,clock_t clock2)
+{
+	double diffticks=clock1-clock2;
+	double diffms=(diffticks)/CLOCKS_PER_SEC;
+	return diffms;
+}
+
+//================================================
+//================================================================
+double get_new_x(double x, double dx){
+    
+    double new_x;
+    // boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
+    double r = rand()/(double)(pow(2.,31)-1.);
+    
+    if (r > 0.5){            
+        new_x = x + rand()*dx/(double)(pow(2.,31)-1.);
+    } else {            
+        new_x = x - rand()*dx/(double)(pow(2.,31)-1.);
+    }
+    
+    return new_x;
+    
+}
+
+
+//===============================================   
+string string_wrap(string ins, int mode){
+    
+    std::ostringstream s;
+    
+    switch(mode){
+        case 0:
+            s << "\033[1;29m" << ins << "\033[0m";
+            break;
+        case 1:
+            s << "\033[1;34m" << ins << "\033[0m";
+            break;
+        case 2:
+            s << "\033[1;44m" << ins << "\033[0m";
+            break;
+        case 3:
+            s << "\033[1;35m" << ins << "\033[0m";
+            break;
+        case 4:
+            s << "\033[1;33;44m" << ins << "\033[0m";
+            break;
+        case 5:
+            s << "\033[1;47;34m" << ins << "\033[0m";
+            break;
+        case 6:
+            s << "\033[1;1;31m" << ins << "\033[0m";
+            break;
+        case 7:
+            s << "\033[1;1;33m" << ins << "\033[0m";
+            break;
+        case 8:
+            s << "\033[1;1;43;34m" << ins << "\033[0m";
+            break;
+        case 9:
+            s << "\033[1;1;37m" << ins << "\033[0m";
+            break;
+        case 10:
+            s << "\033[1;30;47m" << ins << "\033[0m";
+            break;
+        default:
+            s << ins;
+    }
+    
+    return s.str();
+}
+
+
+//===============================================
+string wrap_double(double val, int mode){
+    
+    std::ostringstream s;
+    s << string_wrap(strDouble(val),mode);
+    
+    return s.str();
+}
+
+
+
+//===============================================
+const     
+string i2string(int i){
+    
+    std::ostringstream s;
+    s << "worker" 
+    << lexical_cast<std::string>(i);
+    
+    return s.str();
+    
+}
+
+//===============================================
+char* i2char(int i){
+    
+    std::ostringstream s;
+    s << "worker" 
+    << lexical_cast<std::string>(i);
+    
+    char* a=new char[s.str().size()+1];
+    memcpy(a,s.str().c_str(), s.str().size());
+    
+    return a;
+}
+
+//================================================
+class Universe {
+	
+private:
+	
+	double alpha_i;
+	double alpha_m;
+	double beta;
+	double gamma;
+	double delta;
+	
+    double TargetNovelty;
+    double CumulativeRelativeLoss;
+    double CRLsquare;
+    string id;
+    
+    
+	int N_nodes;
+	int M_edges;
+	
+	int N_epochs;
+	int N_steps;
+	int N_repeats;
+	
+	int current_epoch;
+	double current_loss;
+	int current_repeat;
+    double current_novelty;
+	
+	int mode_identify_failed;
+    int verbose_level; // 0 is silent, higher is more
+	
+	double k_max;
+	
+	graph_t Full_g;
+	
+	double **Prob;
+	double **Tried;
+	double **Dist;
+	double **Final;
+    double **EdgeIndex;
+	double *Rank;
+	
+    base_generator_type generator;	
+    boost::uniform_real<> uni_dist;
+    boost::geometric_distribution<double> geo;
+    
+public:
+	
+    
+    
+	//======  Constructor ======
+	Universe(const std::string FileToOpen, int Epochs, int Steps, int Repeats, int identify_failed, double target, const std::string idd)
+	{
+		//typedef array_type2::index index2;
+		
+		
+		std::ifstream inFile;
+        //string line;
+        
+        //-------------------------------
+        
+        base_generator_type gene(42u);
+        generator = gene;
+        generator.seed(static_cast<unsigned int>(std::time(0)));
+        boost::uniform_real<> uni_d(0,1);
+        uni_dist = uni_d;
+        
+        //--------------------------------
+        
+		int i, k;
+        int x, y;
+		Edge* edge_array_mine;
+		int num_arcs_mine, num_nodes_mine;
+		int* weights_mine;
+		
+        TargetNovelty = target;
+        CumulativeRelativeLoss = 0.;
+        CRLsquare = 0.;
+        
+		
+		N_epochs  = Epochs;
+		N_steps   = Steps;
+		N_repeats = Repeats;
+		
+		current_epoch = 0;
+		current_loss = 0.;
+		current_repeat = 0;
+        
+        id = idd;
+        
+        verbose_level = 1;
+		
+		mode_identify_failed = identify_failed;
+        
+		
+		//-------------------------------
+		// The first pass though file with the graph
+		inFile.open(FileToOpen.c_str());
+		if (inFile.fail()) {
+			cout << "Unable to open file";
+			exit(1); // terminate with error
+		}else {
+            
+            if (verbose_level > 2){
+                std::cout <<  " Opened <" << FileToOpen << ">"<<std::endl;
+            }
+        }
+		
+		i=0;
+        std::string line;
+		//while (! inFile.eof() && ! inFile.fail()) {
+        while (1==1) {
+            
+            inFile >> x;
+            inFile >> y;
+            
+            if (verbose_level > 2){
+                std::cout << " x: " << x;
+                std::cout << " y: " << y << std::endl;
+            }
+            
+			if (i==0){
+				N_nodes=x;
+				M_edges=y;	
+                break;
+			}
+			i++;
+            
+			
+		}
+		inFile.close();
+        
+        if (verbose_level == 2){
+            std::cout << N_nodes <<  " nodes, " << M_edges << " edges"<<std::endl;
+        }
+		
+		// k_max is the longest distance possible
+		
+        //k_max = M_edges;
+		k_max = 70;
+        
+		//------------------------------------
+		// Get memory allocated for all class members
+		
+		Prob = allocate_2Dmatrix(N_nodes, N_nodes);
+		Tried = allocate_2Dmatrix(N_nodes, N_nodes);
+		Dist = allocate_2Dmatrix(N_nodes, N_nodes);
+		Final = allocate_2Dmatrix(N_nodes, N_nodes);
+        EdgeIndex = allocate_2Dmatrix(N_nodes, N_nodes);
+		Rank = allocate_1Dmatrix(N_nodes);
+		
+        //The second pass through file with the graph
+        
+		for(int i = 0; i < N_nodes; ++i) {
+			Rank[i]=0.;
+			for(int j = 0; j < N_nodes; ++j) {
+				Final[i][j] = 0.;
+				Prob[i][j]=0.;
+				Dist[i][j]=-1.;
+				Tried[i][j]=0.;
+                EdgeIndex[i][j]=-1;
+			}
+		}
+        
+		
+		// Fill in the final graph -- and we are ready to go!
+        
+	    inFile.open(FileToOpen.c_str());
+		if (!inFile) {
+            std::cout << "Unable to open file";
+			exit(1); // terminate with error
+		}
+		else {
+            
+            if (verbose_level > 2){
+                std::cout <<  " Opened <" << FileToOpen << ">"<<std::endl;
+            }
+        }
+        
+		i=0;  
+		while (inFile >> x && inFile >>y) {
+			if (i > 0) {
+				Final[x][y]=1.;
+				Final[y][x]=1.;
+                
+                
+                if (verbose_level == 2){
+                    std::cout << ".";
+                }
+			}
+			i++;
+			
+		}
+        if (verbose_level == 2){
+            std::cout << std::endl;
+        }
+		inFile.close(); 
+		
+        k=0;
+        for (int i=0; i<N_nodes-1; i++){
+            for (int j=i+1;j<N_nodes; j++){
+                if(Final[i][j] > 0.){
+                    EdgeIndex[i][j]=k;
+                    k++;
+                }
+            }
+        }
+        
+        
+		
+		//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+		// create graph -- hopefully, we can keep it, just modifying edge weights
+		
+		
+		edge_array_mine = new Edge[2*M_edges];
+		num_arcs_mine = 2*M_edges;
+		num_nodes_mine = N_nodes;
+		weights_mine = new int[2*M_edges];
+		for (int i=0; i<2*M_edges; i++){ weights_mine[i]=1;}
+		
+		k=0;
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				if (Final[i][j]>0.){
+					edge_array_mine[2*k]  =Edge(i,j);
+					edge_array_mine[2*k+1]=Edge(j,i);
+					k++;
+				}
+			}
+		}
+		graph_t g(edge_array_mine, edge_array_mine + num_arcs_mine, weights_mine, num_nodes_mine);
+		
+		Full_g = g;
+		delete edge_array_mine;
+		delete weights_mine;
+		
+		//===========================================================================
+		std::vector<edge_descriptor> p(num_edges(Full_g));
+		std::vector<int> d(num_edges(Full_g));
+		edge_descriptor s;
+		boost::graph_traits<graph_t>::vertex_descriptor u, v;
+		
+		for (int i=0; i<N_nodes-1; i++){
+			for (int j=i+1; j<N_nodes; j++){
+				if (Final[i][j] > 0.){
+					u = vertex(i, Full_g);
+					v = vertex(j, Full_g);
+					remove_edge(u,v,Full_g);
+					remove_edge(v,u,Full_g);
+					
+				}
+			}
+		}
+		
+        
+    }
+	
+	
+	//=====================================================================
+	int sample_failed_number(double pfail){
+		
+		//boost::geometric_distribution<double> geo(pfail);
+		//boost::variate_generator<base_generator_type&, geometric_distribution<double> > geom(generator, geo);
+		
+		double r, u, g;
+        
+        r=0.;
+		for(int i=0; i<N_steps; i++){
+            
+            u=(double)rand();
+            u = 1.-u /(double)(pow(2.,31)-1.);
+            g=(int)(ceil(log(u) / log(pfail)));
+            
+			//r += geom();
+            
+            r+=g;
+		}
+        
+        if (verbose_level>=3){
+            std::cout << id << " failed " << r << std::endl;
+		}
+		return r;
+		
+	}
+    
+    //=============================================
+    double get_target(void){
+        return TargetNovelty;
+    }
+	
+    //=============================================
+    void set_target(double target){
+        TargetNovelty=target;
+    }
+	
+	//=============================================
+	int sample(){
+		
+        //boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
+        // double r = uni(), Summa = 0.;
+        
+        
+        
+        double r = rand(), Summa = 0.;
+        r /= (double)(pow(2.,31)-1.);
+		int result = 0;
+		int finished = 0;
+        
+        if (verbose_level==4){
+            std::cout << id << " sampled " << r << std::endl;
+        }
+		
+		for(int i=0; i<N_nodes-1 && finished==0; i++){			
+			for( int j=i+1; j<N_nodes && finished==0; j++){
+				
+				Summa += Prob[i][j];
+				
+				if (Summa > r){
+					
+					Tried[i][j]+=1.;
+					
+					if (Final[i][j] > 0.){
+						result = 1;
+					}
+					finished = 1;										
+				}
+			}
+		}
+		
+		return result;
+		
+	}
+	
+	//===============================
+	void update_current_graph(void){
+		
+		std::vector<edge_descriptor> p(num_edges(Full_g));
+		std::vector<int> d(num_edges(Full_g));
+		edge_descriptor s;
+		boost::graph_traits<graph_t>::vertex_descriptor u, v;
+		
+		//property_map<graph_t, edge_weight_t>::type weightmap = get(edge_weight, Full_g);
+		for (int i=0; i<N_nodes-1; i++){
+			for (int j=i+1; j<N_nodes; j++){
+				if (Final[i][j] > 0. && Tried[i][j]>0){
+					//s = edge(i, j, Full_g);	
+					boost::graph_traits<graph_t>::edge_descriptor e1,e2;
+					bool found1, found2;
+					u = vertex(i, Full_g);
+					v = vertex(j, Full_g);
+					tie(e1, found1) = edge(u, v, Full_g);
+					tie(e2, found2) = edge(v, u, Full_g);
+					if (!found1 && !found2){
+						add_edge(u,v,1,Full_g);
+					    add_edge(v,u,1,Full_g);
+					}
+					
+				}
+			}
+			
+		}
+	}
+	
+	//===============================
+	void update_distances(void){
+		// put shortest paths to the *Dist[][]
+		std::vector<vertex_descriptor> p(num_vertices(Full_g));
+		std::vector<int> d(num_vertices(Full_g));
+		vertex_descriptor s;
+		
+		
+		// put shortest paths to the *Dist[][]
+		for (int j=0; j<num_vertices(Full_g); j++){
+			
+			if(Rank[j] > 0.){
+				s = vertex(j, Full_g);	 
+				dijkstra_shortest_paths(Full_g, s, predecessor_map(&p[0]).distance_map(&d[0]));
+				
+				//std::cout <<" Vertex "<< j << std::endl;
+				graph_traits < graph_t >::vertex_iterator vi, vend;
+				
+				for (boost::tie(vi, vend) = vertices(Full_g); vi != vend; ++vi) {
+					
+					if (p[*vi]!=*vi){
+						Dist[*vi][j]=d[*vi];
+						Dist[j][*vi]=d[*vi];
+                        
+                        if (Dist[*vi][j]>max_dist){
+                            max_dist=Dist[*vi][j];
+                        }
+                        
+                        
+					} else {
+						Dist[*vi][j]=-1.;
+						Dist[j][*vi]=-1.;
+					}
+				}
+			}
+			
+		}
+		
+		
+	}
+	
+	//======================================================
+	void update_ranks(void){
+		
+		for(int i=0; i<N_nodes; i++){
+			Rank[i]=0.;
+		}
+		
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				if (Tried[i][j]>0. && Final[i][j] >0.){
+					Rank[i]++;
+					Rank[j]++;
+				}
+			}
+		}
+		
+	}
+	
+	//====================================================================
+	void set_world(double a_i, double a_m, double b, double g, double d){
+		
+		alpha_i=a_i;
+		alpha_m=a_m;
+		gamma=g;
+		beta=b;
+		delta=d;
+		
+	}
+	
+	//====================================================================
+	void reset_world(){
+		
+        //====================================================
+		std::vector<edge_descriptor> p(num_edges(Full_g));
+		std::vector<int> d(num_edges(Full_g));
+		edge_descriptor s;
+		boost::graph_traits<graph_t>::vertex_descriptor u, v;
+        
+		
+		for (int i=0; i<N_nodes-1; i++){
+			for (int j=i+1; j<N_nodes; j++){
+				if (Final[i][j] > 0. && Tried[i][j] > 0){
+					u = vertex(i, Full_g);
+					v = vertex(j, Full_g);
+					remove_edge(u,v,Full_g);
+					remove_edge(v,u,Full_g);
+					
+				}
+			}
+		}
+        
+        //==================================================
+        
+		current_loss=0;
+		current_epoch=0;
+		current_repeat++;
+        current_novelty=0;
+		
+		for(int i = 0; i < N_nodes; ++i) {
+			Rank[i]=0.;
+			for(int j = 0; j < N_nodes; ++j) {
+				Prob[i][j]=0.;
+				Dist[i][j]=-1.;
+				Tried[i][j]=0.;
+			}
+		}
+	}
+	
+	
+    //==============================================
+    void show_parameters(void){
+        
+        std::cout << "Parameters: " 
+        << alpha_i << " "
+        << alpha_m << " | "
+        << beta << " "
+        << gamma << " | "
+        << delta << std::endl;
+        
+    }
+    
+    
+    
+    //===============================================
+    string file_name(){
+        
+        std::ostringstream s;
+        s << "world_" 
+        << lexical_cast<std::string>(alpha_i) << "_" 
+        << lexical_cast<std::string>(alpha_m) << "_"
+        << lexical_cast<std::string>(beta) << "_"
+        << lexical_cast<std::string>(gamma) << "_"
+        << lexical_cast<std::string>(delta) << "_"
+        << lexical_cast<std::string>(N_epochs) << "_"
+        << lexical_cast<std::string>(N_steps) << "_"
+        << lexical_cast<std::string>(N_repeats) << ".txt";
+        
+        return s.str();
+        
+    }
+    
+    
+    
+    
+    //=================================================
+    void set_verbose(int verbose){
+        
+        verbose_level = verbose;
+    }
+    
+    
+    //=============================================================
+    void update_probabilities(void){
+        
+        
+        //=========================
+		// Compute sampling probabilities
+		// first pass: \xi_i,j
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				
+				double bg = 0.;
+				
+				Prob[i][j] = alpha_i*log(min(Rank[i]+1.,Rank[j]+1.)) + 
+                alpha_m*log(max(Rank[i]+1.,Rank[j]+1.));
+				
+                if (Dist[i][j] > 0.){
+                    
+                    double k = Dist[i][j];
+                    if (k >= k_max){
+                        k = k_max-1;
+                    }
+					
+                    bg = beta * log(k/k_max) + gamma * log(1. - k/k_max);
+					
+                } else {
+                    bg = delta;
+                }
+				
+				Prob[i][j] = exp(Prob[i][j] + bg);
+			}
+		}
+        
+		
+		// second pass: sum
+		double Summa = 0.;
+		
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				Summa += Prob[i][j];
+			}
+		}
+		
+		// third pass: normalize
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				Prob[i][j] /= Summa;
+			}
+		}
+        
+    }
+    
+	// Now we are ready for simulations
+	//==============================================
+	void update_world(){
+		
+		int failed = 0;
+        
+		// Given current universe compute shortest paths
+		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+		
+		update_current_graph();
+		update_ranks();				
+		update_distances();
+		update_probabilities();
+		
+		//===============================
+		// sampling
+		int result;
+		double cost=0., novel=0.;
+		int publishable = 0;
+		
+		
+		//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+		if (mode_identify_failed == 1){
+			
+			while(publishable < N_steps){
+				
+		    	result = sample();
+			    publishable += result;
+			    failed += (1-result);
+				
+			}
+			
+			for(int i=0; i<N_nodes-1; i++){
+				for( int j=i+1; j<N_nodes; j++){
+					
+					cost+=Tried[i][j];
+					
+					if (Tried[i][j]>0. && Final[i][j]>0.){
+						novel+=1.;
+					}
+				}
+			}
+			
+		}
+		//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+		else {
+			
+			double pfail=0.;
+			int n_failed;
+			//, n_check = 0;
+			
+			for(int i=0; i<N_nodes-1; i++){
+				for( int j=i+1; j<N_nodes; j++){
+					if (Final[i][j] == 0.){
+						pfail += Prob[i][j];
+						Prob[i][j] = 0.;
+					}
+					
+				}
+			}
+			
+			for(int i=0; i<N_nodes-1; i++){
+				for( int j=i+1; j<N_nodes; j++){
+					Prob[i][j] /= (1.-pfail);
+				}
+				//std::cout << std::endl;
+			}			
+			
+			n_failed = sample_failed_number(pfail);
+			while(publishable < N_steps){
+				
+		    	result = sample();
+			    publishable += result;					
+			}
+            
+            
+			current_loss += (n_failed + N_steps);
+			cost = current_loss;
+			
+			for(int i=0; i<N_nodes-1; i++){
+				for( int j=i+1; j<N_nodes; j++){
+					
+					if (Tried[i][j]>0. && Final[i][j]>0.){
+						novel+=1.;
+					}
+				}
+			}
+		}
+        
+        current_novelty = novel;
+        
+        
+		//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+		if (verbose_level == 2){
+            std::cout << (current_repeat+1) << "  epoch=" << (current_epoch+1) 
+            
+		    << "  cost=" << cost 
+		    << " novel=" << novel 
+		    << " rel_loss=" << cost/novel
+		    << std::endl;
+        }
+        
+		current_epoch++;
+	}
+	
+	
+	//======  Destructor ======
+	~Universe(){
+		
+		delete_2Dmatrix(Final, N_nodes);
+		delete_2Dmatrix(Dist, N_nodes);
+		delete_2Dmatrix(Tried, N_nodes);
+		delete_2Dmatrix(Prob, N_nodes);
+        delete_2Dmatrix(EdgeIndex, N_nodes);
+		delete_1Dmatrix(Rank);		
+	}
+    
+    //================================================
+    // Allocate memory
+    double** allocate_2Dmatrix(int N, int M)
+    {
+        double **pointer;
+        
+        if (verbose_level == 2){
+            std::cout<< "["<<N<<"|"<<M<<"]"<<std::endl;
+        }
+        pointer = new double*[N];
+        for (int i = 0; i < N; ++i)
+            pointer[i] = new double[M];	
+        
+        return pointer;
+    }
+    //===================
+    double* allocate_1Dmatrix(int N)
+    {
+        double *pointer;
+        
+        if(N > 0){
+            
+            pointer = new double[N];
+            
+        }else {
+            
+            pointer = NULL;
+        }
+        
+        return pointer;
+        
+    }
+    
+    //==============================================
+    // De-Allocate memory to prevent memory leak
+    void delete_2Dmatrix(double **pointer, int N){
+        
+        if (pointer != NULL){
+            
+            for (int i = 0; i < N; ++i){
+                delete [] pointer[i];
+            }
+            delete [] pointer;
+        }
+    }
+    //====================
+    void delete_1Dmatrix(double *pointer){
+        
+        delete [] pointer;
+    }
+    
+    //===========================================
+    double get_rel_loss(){
+        
+        return CumulativeRelativeLoss ;
+    }
+    
+    //===========================================
+    double get_rel_loss_err(){
+        
+        return CRLsquare ;
+    }
+    
+    
+    
+    //==================================================================================
+    void evolve_to_target_and_save(int istart, int iend, double* storage, int* counters){
+        
+        double ALOT=100000000000.;
+        
+	std::cout<<" evolve_to_target_and_save: istart=" << istart << "iend=" << iend << "\n";
+
+        reset_world();
+        
+        for (int k = istart; k < iend; k++){
+
+	  std::cout<<" evolve: k=" << k << "\n";
+
+            
+            
+            for(int i=0; i< N_epochs &&  current_novelty < TargetNovelty; i++){
+	      std::cout<<" evolve: k=" << k << " i=" << i << " cur=" << current_novelty << " Target=" << TargetNovelty << "\n";
+                update_world();
+            }
+            
+            storage[k]=current_loss/current_novelty;
+            counters[k]=1;
+            
+            
+            reset_world();
+        }
+        
+    }
+    //==============================================
+    int get_reruns(void){
+        return N_repeats;
+    }
+    
+    //==============================================
+    double get_parameter(int i){
+        
+        switch(i){
+            case 0:
+                return alpha_i;
+            case 1:
+                return alpha_m;
+            case 2:
+                return beta;
+            case 3:
+                return gamma;
+            case 4:
+                return delta;
+            default:
+                
+                std::cout << "Erroneous parameter id!!!!\n\n\n";
+                return 0.;
+        }
+    }
+    
+    
+    //==============================================
+    void evolve_to_target(){
+        
+        reset_world();
+        if (beta < -1. || gamma < -1.){
+            CumulativeRelativeLoss = 100000000000.;
+            CRLsquare = 0.;
+            return;
+        }
+        
+        
+        for (int k=0; k< N_repeats; k++){
+            
+            
+            for(int i=0; i<N_epochs &&  current_novelty < TargetNovelty; i++){
+                update_world();
+            }
+            
+            CumulativeRelativeLoss += current_loss/current_novelty;
+            CRLsquare += (current_loss/current_novelty)*(current_loss/current_novelty);
+            if(verbose_level==3){
+                std::cout <<  CumulativeRelativeLoss << " | " << CRLsquare << std::endl;
+            }
+            
+            if(verbose_level==1){
+                std::cout <<  "." ;
+                
+            }
+            else if(verbose_level==2){
+                std::cout <<  "**" << (k+1) <<  "**  curr loss " << current_loss << "; curr novelty " << current_novelty << std::endl;
+            }
+            
+            
+            reset_world();
+        }
+        
+        CumulativeRelativeLoss /= double(N_repeats);
+        CRLsquare /= double(N_repeats);
+        
+        if(verbose_level==1){
+            std::cout << std::endl;
+        }
+        
+        if(verbose_level==2){
+            std::cout <<  CumulativeRelativeLoss << " || " << CRLsquare << std::endl;
+        }
+        
+        CRLsquare = 2*sqrt((CRLsquare - CumulativeRelativeLoss*CumulativeRelativeLoss)/double(N_repeats));
+        
+    }
+    
+    
+    //================================================================    
+    int set_parameter(double value, int position){
+        
+        if (position < 0 || position > 4) {return 0;}
+        
+        else {
+            
+            switch(position){
+                case 0:
+                    alpha_i=value;
+                    return 1;
+                case 1:
+                    alpha_m=value;
+                    return 1;
+                case 2:
+                    beta=value;
+                    return 1;
+                case 3:
+                    gamma=value;
+                    return 1;
+                case 4:
+                    delta=value;
+                    return 1;
+            }
+            
+        }
+        
+        return 0;
+    }
+    
+    
+    //=================================================================
+    void try_annealing(double starting_jump, int iterations, 
+                       double temp_start, double temp_end, double target_rejection){
+        
+        double dx[5]={0.,0.,0.,0.,0};
+        double x[5]={0.,0.,0.,0.,0};
+        double rejection[5]={0., 0., 0., 0., 0.};
+        double curr_x, curr_err, x_tmp;
+        double temperature;
+        double ratio, r;
+        int cycle=10;
+        boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
+        
+        // set up parameter for annealing
+        
+        x[0]=alpha_i;
+        x[1]=alpha_m;
+        x[2]=beta;
+        x[3]=gamma;
+        x[4]=delta;
+        
+        for(int i=0;i<5;i++){
+            dx[i] = starting_jump;
+        }
+        
+        // establish the current value
+        
+        //..........................................
+        evolve_to_target();        
+        std::cout << CumulativeRelativeLoss << " +- " << CRLsquare << std::endl;
+        
+        curr_x   = CumulativeRelativeLoss;
+        curr_err = CRLsquare;
+        CumulativeRelativeLoss = 0;
+        CRLsquare = 0;
+        //...........................................
+        
+        // optimization cycle
+        for(int i=0; i<iterations; i++){
+            
+            temperature = temp_start*exp( i*(log(temp_end)-log(temp_start))/(double)iterations);
+            std::cout  << std::endl << "....T = " << wrap_double(temperature,3) << std::endl << std::endl;
+            
+            if (i % cycle == 0 && i > 0){
+                
+                for (int k=0; k<5; k++){
+                    
+                    rejection[k]/=(double)cycle;
+                    if (rejection[k] > 0){
+                        dx[k] = dx[k]/(rejection[k]/target_rejection);
+                        rejection[k]=0.;
+                    }
+                    else{
+                        dx[k]*=2.;
+                    }
+                    std::cout  << dx[k] << " ";
+                }
+                std::cout  << std::endl;
+            }
+            
+            
+            for (int j=0; j<5; j++){
+                
+                // get new value of x[j]
+                x_tmp = get_new_x(x[j],dx[j]);
+                
+                
+                
+                //.............................................
+                set_parameter(x_tmp, j);
+                
+                
+                evolve_to_target(); 
+                
+                std::cout  << std::endl << "......... " << std::endl;
+                std::cout << "Trying... " << CumulativeRelativeLoss << " +- " << CRLsquare << std::endl;
+                
+                ratio = min(1.,exp(-(CumulativeRelativeLoss-curr_x)/temperature));
+                r = uni();
+                std::cout << r << " vs " << ratio << std::endl;
+                
+                if (r > ratio){
+                    
+                    std::cout << string_wrap(id, 4) <<" "<< (i+1) << ","<< (j) 
+                    <<" "<< (i+1) << " Did not accept " 
+                    << x_tmp << "(" << j << ")" << std::endl;
+                    std::cout << alpha_i << " "<< alpha_m << " "
+                    << beta << " " << gamma << " " 
+                    << delta << " " << std::endl;
+                    set_parameter(x[j], j);
+                    CumulativeRelativeLoss = 0;
+                    CRLsquare = 0;
+                    
+                    rejection[j]+=1.;
+                }
+                
+                else {
+                    
+                    curr_x   = CumulativeRelativeLoss;
+                    curr_err = CRLsquare;
+                    x[j] = x_tmp;
+                    CumulativeRelativeLoss = 0;
+                    CRLsquare = 0;
+                    std::cout << (i+1) << string_wrap((string) " Rejection counts: ", 8) 
+                    << wrap_double(rejection[0],2) 
+                    << " "<< wrap_double(rejection[1], 7) << " "
+                    << wrap_double(rejection[2],5) << " " << wrap_double(rejection[2],9) << " " 
+                    << wrap_double(rejection[4],6) << " " 
+                    << std::endl << std::endl;
+                    
+                    std::cout << string_wrap(id, 4) <<" "<< (i+1) <<","<< (j) 
+                    <<" "
+                    << string_wrap((string) "***** Did accept! ", 3) 
+                    << wrap_double(alpha_i,2) 
+                    << " "<< wrap_double(alpha_m, 7) << " "
+                    << wrap_double(beta,5) << " " 
+                    << wrap_double(gamma,9) << " " 
+                    << wrap_double(delta,6) << " " 
+                    << std::endl << std::endl;
+                    
+                }
+                //........................................................ 
+                
+            }
+            
+        }
+        
+    }
+    
+	
+};
+
+//============================================================
+
+std::pair<double,double> multi_loss( // dispatch_group_t group, 
+                                    Universe* un[], 
+                                    // dispatch_queue_t* CustomQueues,
+                                    double* Results,
+                                    int*    Counters,
+                                    double* params){
+    
+    int N = un[0]->get_reruns();
+    int step = (int)(double)N/(double)(Nworkers);
+    int istart=0;
+    int iend = istart+step;
+    
+    double Loss=0., LossSquare=0.;
+    
+    for(int i=0; i<Nworkers; i++){
+        
+        
+        for(int j=0; j<5; j++){
+            un[i]->set_parameter(params[j],j);
+        }
+        
+        
+        for(int i=0; i<Nworkers; i++){
+            
+	  // dispatch_group_async(group, CustomQueues[i], ^{
+	  std::cout<<" Calling evolve_to_target_and_save\n";
+                
+                un[i]->evolve_to_target_and_save(istart, iend, Results, Counters);
+		//});
+            
+	  std::cout<<" Returned from evolve_to_target_and_save\n";
+            istart += step;
+            iend = min(istart+step,N);
+            
+        }
+    }
+    // dispatch_group_wait(group, DISPATCH_TIME_FOREVER);
+    //dispatch_release(group);
+    
+    for (int i=0; i<N; i++){
+        
+        Loss+=Results[i]/(double)N;
+        LossSquare+=Results[i]*Results[i]/(double)N;
+        
+        std::cout<<" RESULT: " << Results[i];
+    }
+    
+    std::cout<<" \n\n\n";
+    double two_std = ((LossSquare - Loss*Loss)/(double)N);
+    
+    two_std = 2.*sqrt(two_std);
+    std::pair<double,double> Res;
+    Res.first=Loss;
+    Res.second=two_std;
+    
+    return Res;
+    
+    
+}
+//============================================================
+
+
+//============================================================
+void multi_annealing( // dispatch_group_t group, 
+                     Universe* un[], 
+                     // dispatch_queue_t* CustomQueues, 
+                     double T_start, double T_end, 
+                     double Target_rejection, 
+                     int Annealing_repeats, 
+                     double starting_jump,
+                     double* Results,
+                     int*    Counters,
+                     double* params0,
+                     double annealing_cycles){
+    //.................................
+    // re-implement annealing
+    
+    double dx[5]={0.,0.,0.,0.,0};
+    double x[5]={0.,0.,0.,0.,0};
+    double rejection[5]={0., 0., 0., 0., 0.};
+    double curr_x, curr_err, x_tmp;
+    double temperature;
+    double ratio, r;
+    int cycle=10;
+    //boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
+    
+    // set up parameter for annealing
+    
+    x[0]=params0[0];
+    x[1]=params0[1];
+    x[2]=params0[2];
+    x[3]=params0[3];
+    x[4]=params0[4];
+    
+    for(int i=0;i<5;i++){
+        dx[i] = starting_jump;
+    }
+    
+    // establish the current value
+    std::pair<double,double>Res;
+    
+    Res = multi_loss( /* group,*/ un, /* CustomQueues,*/ Results, Counters, x);
+    std::cout << Res.first << " +- " << Res.second << std::endl;
+    
+    curr_x   = Res.first;
+    curr_err = Res.second;
+    
+    // optimization cycle
+    
+    for(int i=0; i<annealing_cycles; i++){
+        
+        temperature = T_start*exp( i*(log(T_end)-log(T_start))/(double)annealing_cycles);
+        std::cout  << std::endl << "....T = " << wrap_double(temperature,3) << std::endl << std::endl;
+        
+        if (i % cycle == 0 && i > 0){
+            
+            for (int k=0; k<5; k++){
+                rejection[k]/=(double)cycle;
+                
+                if (rejection[k] > 0){
+                    dx[k] = dx[k]/(rejection[k]/Target_rejection);
+                    rejection[k]=0.;
+                }
+                else{
+                    dx[k]*=2.;
+                }
+                std::cout  << dx[k] << " ";
+            }
+            std::cout  << std::endl;
+        }
+        
+        
+        for (int j=0; j<5; j++){
+            
+            ///////////////////////////////
+            if (FIX_VARIABLES==0 || var_fixed[j]==0){
+                
+                
+                
+                // get new value of x[j]
+                double x_hold=x[j];
+                x_tmp = get_new_x(x[j],dx[j]);
+                x[j]=x_tmp;
+                
+                std::cout << wrap_double(x_tmp,10) << " " << wrap_double(j,9) << "\n\n"; 
+                //=======================================
+                //.............................................
+                for(int w=0; w<Nworkers; w++){
+                    un[w]->set_parameter(x_tmp, j);
+                }
+                
+                
+                Res = multi_loss(/* group, */ un, /* CustomQueues, */ Results, Counters, x);
+                std::cout << Res.first << " +- " << Res.second << std::endl;
+                
+                ratio = min(1.,exp(-(Res.first-curr_x)/temperature));
+                r = rand()/(double)(pow(2.,31)-1.);
+                std::cout << r << " vs " << ratio << std::endl;
+                
+                double ALOT=100000000000.;
+                
+                if (Res.first < ALOT)
+                {
+                    ofstream filestr;
+                    
+                    filestr.open ("best_opt_some.txt", ofstream::app);
+                    
+                    // >> i/o operations here <<
+                    filestr << un[0]->get_target() << "," 
+                    << Res.first 
+                    << "," << un[0]->get_parameter(0) 
+                    << "," << un[0]->get_parameter(1) 
+                    << "," << un[0]->get_parameter(2) 
+                    << "," << un[0]->get_parameter(3) 
+                    << "," << un[0]->get_parameter(4) << "," << Res.second << ",\n";
+                    
+                    filestr.close();
+                    
+                    
+                    filestr.open ("max_dist.txt", ofstream::app);
+                    
+                    // >> i/o operations here <<
+                    filestr << max_dist << ",\n";
+                    
+                    filestr.close();
+                    
+                }
+                
+                
+                if (r > ratio){
+                    
+                    std::cout << " "<< (i+1) << ","<< (j) 
+                    <<" "<< (i+1) << " Did not accept " 
+                    << x_tmp << "(" << j << ")" << std::endl;
+                    std::cout << un[0]->get_parameter(0) 
+                    << " " << un[0]->get_parameter(1) 
+                    << " " << un[0]->get_parameter(2) 
+                    << " " << un[0]->get_parameter(3) 
+                    << " " << un[0]->get_parameter(4) << " " << std::endl;
+                    
+                    x[j]=x_hold;
+                    for(int w=0; w<Nworkers; w++){
+                        un[w]->set_parameter(x[j], j);
+                    }
+                    
+                    
+                    //set_parameter(x[j], j);     
+                    rejection[j]+=1.;
+                }
+                
+                else {
+                    
+                    curr_x   = Res.first;
+                    curr_err = Res.second;
+                    x[j] = x_tmp;
+                    
+                    for(int w=0; w<Nworkers; w++){
+                        un[w]->set_parameter(x[j], j);
+                    }
+                    
+                    std::cout << (i+1) << string_wrap((string) " Rejection counts: ", 8) 
+                    << wrap_double(rejection[0],2) << " " 
+                    << wrap_double(rejection[1],7) << " "
+                    << wrap_double(rejection[2],5) << " " 
+                    << wrap_double(rejection[3],9) << " " 
+                    << wrap_double(rejection[4],6) << " " 
+                    << std::endl << std::endl;
+                    
+                    std::cout << " "<< (i+1) <<","<< (j) 
+                    <<" "
+                    << string_wrap((string) "***** Did accept! ", 3) 
+                    << wrap_double(un[0]->get_parameter(0),2) << " "
+                    << wrap_double(un[0]->get_parameter(1),7) << " "
+                    << wrap_double(un[0]->get_parameter(2),5) << " " 
+                    << wrap_double(un[0]->get_parameter(3),9) << " " 
+                    << wrap_double(un[0]->get_parameter(4),6) << " " 
+                    << std::endl << std::endl;
+                    
+                    
+                    
+                }
+                //........................................................ 
+                
+            }
+        }
+        
+    }
+    
+}
+
+
+
+//================================================
+int
+main(int argc, char* argv[])
+{
+    
+    double params0[6] = {0., 0., 0., 0., 0., 0.2}, target=50., range;
+    string par_names0[6] = {"alpha_i", "alpha_m", "beta", "gamma", "delta", "target"};
+    string par_names1[4] = {"n_epochs", "n_steps", "n_reruns", "range"};
+    string par_names2[5] = {"T_start", "T_end", "Annealing_steps","Target_rejection","Starting_jump"};
+    string par_names3[5] = {"FREEZE_alpha_i", "FREEZE_alpha_m", "FREEZE_beta", "FREEZE_gamma", "FREEZE_delta"};    
+    int params1[4] = {300, 50, 1000, 10};
+    int params3[5] = { 0, 0, 0, 0, 0};
+    
+    //          temperature_start,  temperature_end,  annealing_steps target_rejection  Starting_jump
+    double params2[5] = {1,             0.001,               100,              0.3,           1.5};
+    
+    int verbose_level = 2;
+    const std::string one="one", two="two";
+    static Universe* un[Nworkers];
+    // static dispatch_queue_t CustomQueues[Nworkers];
+    
+    static double* Results;
+    static int*    Counters;
+    
+    timeval t1, t2;
+    double elapsedTime;
+    // start timer
+    gettimeofday(&t1, NULL);
+    
+    
+    if (argc < 8) {
+        std::cout << "Usage: super_optimizer alpha_i alpha_m beta gamma delta target_innov [n_epochs n_steps n_reruns] [range] [verbose_level]\n"; 
+        std::cout << "         [T_start T_end Annealing_steps Target_rejection Starting_jump]\n"; 
+        std::cout << "         [FREEZE_alpha_i FREEZE_alpha_m FREEZE_beta FREEZE_gamma FREEZE_delta]\n"; 
+        
+        system("pwd");
+        
+        
+        return(1);
+    }
+    else {
+        for (int nArg=0; nArg < argc; nArg++){
+            //std::cout << nArg << " " << argv[nArg] << std::endl;
+            if (nArg > 0 && nArg < 7){
+                params0[nArg-1]= atof(argv[nArg]);
+                std::cout << par_names0[nArg-1] << ": " << params0[nArg-1] <<  std::endl;
+            }
+            if (nArg > 6 && nArg < 11){
+                params1[nArg-7]= atoi(argv[nArg]);
+                std::cout << par_names1[nArg-7] << ": " << params1[nArg-7] <<  std::endl;
+            }
+            if (nArg == 11){
+                verbose_level = atoi(argv[nArg]);
+                std::cout << "verbose level: " << verbose_level <<  std::endl;
+            }
+            if (nArg > 11 && nArg < 17){
+                params2[nArg-12]= atof(argv[nArg]);
+                std::cout << par_names2[nArg-12] << ": " << params2[nArg-12] <<  std::endl;
+            }
+            if (nArg > 16 && nArg < 22){
+                params3[nArg-17]= atof(argv[nArg]);
+                var_fixed[nArg-17]= atof(argv[nArg]);
+                std::cout << par_names3[nArg-17] << ": " << var_fixed[nArg-17] <<  std::endl;
+            }
+            
+            
+        }
+        
+    }
+    
+    /*
+    for target in range(58,1009,50):
+	s = ("%d" % target)
+	print s
+
+	for i in range(15):
+          # Param groups separated by "|" below for documentation. NOTE that | is not used on command line!
+	  os.system("./supe_duper_optimizer |0 0 4 50 -1 "+s+" | 40000 20 1000 2 | 1 | 2. 0.01 100 0.3 2.3 | 1 1 0 0 0")
+    */
+
+    /*   Parameters, re-iterated:
+
+                        "alpha_i", "alpha_m", "beta", "gamma", "delta", "target"};
+    double params0[6] = {0.,        0.,        0.,     0.,      0.,      0.2}, target=50., range;
+
+                    {"n_epochs", "n_steps", "n_reruns", "range"};
+    int params1[4] = {300,        50,       1000,        10};
+
+                       {"T_start",          "T_end",          "Annealing_steps","Target_rejection","Starting_jump"};
+                         temperature_start,  temperature_end,  annealing_steps target_rejection  Starting_jump
+    double params2[5] = {1,                  0.001,            100,            0.3,              1.5};
+
+                     {"FREEZE_alpha_i", "FREEZE_alpha_m", "FREEZE_beta", "FREEZE_gamma", "FREEZE_delta"};    
+    int params3[5] = { 0,                0,                0,             0,              0};
+
+
+    */
+
+    for (int j=0; j<5; j++){
+        
+        cout << j << " | " << var_fixed[j] << " (fixed) \n";
+    }
+	
+    target=params0[5];
+    range = (double)params1[3];
+	int identify_failed = 0;
+	char* filename= (char *)"movie_graph.txt";
+ 	int n_ep=params1[0], n_st=params1[1], n_rep=params1[2];
+    
+    //...............................
+    
+    for(int i=0; i<Nworkers; i++){
+        un[i] = new Universe((char *)filename,n_ep,n_st,
+                             (int)n_rep,
+                             identify_failed, target, i2string(i));
+        // CustomQueues[i] = dispatch_queue_create(i2char(i), NULL);
+    }
+    
+    //...............................
+    if(n_rep > 0){
+        
+        Results = new double[n_rep];
+        Counters = new int[n_rep];
+        
+    }else {
+        
+        Results =  NULL;
+        Counters = NULL;
+        std::cout << " Number of reruns should be positive! " <<  std::endl;
+        return 0;
+        
+    }
+    //...............................
+    srand(time(0));
+    //srandomdev();
+    
+    {
+        double r=0;
+        for (int j=0; j<100; j++){
+            
+            
+            
+            r = rand()/(double)(pow(2.,31)-1.);
+            std::cout << r << " ";
+        }
+        std::cout << "\n ";
+    }
+  	//random initiation of starting parameters
+    
+    if (range > 0.){
+        
+        for (int i=0; i < 5; i++){
+            
+            if (params0[i]==-100.){
+                
+                double r1 = (rand()/(double)(pow(2.,31)-1.));
+                double r2 = (rand()/(double)(pow(2.,31)-1.));
+                double sign = 1.;
+                
+                if(r1 > 0.5){
+                    sign=-1.;
+                }
+                
+                params0[i] = sign*r2*range;
+                
+                std::cout << par_names0[i] << ": " << params0[i] <<  std::endl;
+            }
+        }
+        
+    }
+    
+    
+    double T_start=params2[0], T_end=params2[1], Target_rejection=params2[3], starting_jump=params2[4];
+    int Annealing_repeats = (int) params2[2];
+    
+    
+    // dispatch_group_t group = dispatch_group_create();
+    
+    //.............................
+    multi_annealing( /* group, */ un, /* CustomQueues, */ T_start, T_end, Target_rejection, Annealing_repeats, 
+                    starting_jump, Results, Counters, params0, Annealing_repeats);
+    
+    //dispatch_group_wait(group, DISPATCH_TIME_FOREVER);
+    // dispatch_release(group);
+    //.............................
+    
+    
+    // stop timer
+    gettimeofday(&t2, NULL);
+    
+    // compute and print the elapsed time in millisec
+    elapsedTime = (t2.tv_sec - t1.tv_sec) * 1000.0;      // sec to ms
+    elapsedTime += (t2.tv_usec - t1.tv_usec) / 1000.0;   // us to ms
+    elapsedTime /= 1000.;
+    cout << elapsedTime << " seconds \n .....(" << elapsedTime/60. << " minutes)\n\n";
+    
+    //.....................
+    
+    for(int i=0; i<Nworkers; i++){
+        delete un[i];
+    }
+    
+    //....................
+    if(n_rep > 0){
+        
+        delete [] Results;
+        delete [] Counters;
+        
+    }
+    
+    return 0;
+	
+	
+	
+}
+


Property changes on: SwiftApps/optimizer.orig.cpp
___________________________________________________________________
Added: svn:executable
   + 

Added: SwiftApps/optimizer.sh
===================================================================
--- SwiftApps/optimizer.sh	                        (rev 0)
+++ SwiftApps/optimizer.sh	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,10 @@
+#1 /bin/sh
+bestfile=$1
+maxfile=$2
+datafile=$3
+touch best_opt_some.txt max_dist.txt multi_loss.data
+shift 3
+$(dirname $0)/Optimizer $* 2>&1
+mv best_opt_some.txt $bestfile
+mv max_dist.txt $maxfile
+mv multi_loss.data $datafile


Property changes on: SwiftApps/optimizer.sh
___________________________________________________________________
Added: svn:executable
   + 

Added: SwiftApps/optimizer.snap01.cpp
===================================================================
--- SwiftApps/optimizer.snap01.cpp	                        (rev 0)
+++ SwiftApps/optimizer.snap01.cpp	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,1722 @@
+//
+//  main.cpp
+//  optimizer
+//
+//  Created by Andrey Rzhetsky on 4/11/11.
+//  Copyright 2011 University of Chicago. All rights reserved.
+//
+
+#define Nworkers 24
+
+#include <fstream>
+#include <iostream>
+#include <stdio.h>
+#include <time.h>
+#include <ctime>    
+#include <algorithm>
+#include <string>
+
+#include <stdio.h>
+#include <sys/param.h>
+#include <sys/time.h>
+#include <sys/types.h>
+
+// #include <dispatch/dispatch.h>
+#include <fstream>
+
+
+#include <stdlib.h>
+#include <boost/numeric/ublas/io.hpp>
+#include <boost/graph/graph_traits.hpp>
+#include <boost/graph/dijkstra_shortest_paths.hpp>
+#include <boost/graph/loop_erased_random_walk.hpp>
+#include <boost/graph/random.hpp>
+#include <boost/property_map/property_map.hpp>
+#include <boost/graph/graph_concepts.hpp>
+#include <boost/graph/properties.hpp>
+
+#include <boost/graph/graph_traits.hpp>
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/adjacency_matrix.hpp>
+
+#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
+#define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
+#include <boost/graph/random.hpp>
+#include <boost/random/geometric_distribution.hpp>
+#include <boost/random/uniform_01.hpp>
+#include <boost/random.hpp>
+#include <boost/random/linear_congruential.hpp>
+#include <boost/random/uniform_int.hpp>
+#include <boost/random/uniform_real.hpp>
+#include <boost/random/variate_generator.hpp>
+#include <boost/generator_iterator.hpp>
+#include <boost/lexical_cast.hpp>
+
+#define INT_INFINITY 2147483647
+
+#define FIX_VARIABLES 1
+
+using namespace boost;
+using namespace std;
+using namespace boost::numeric::ublas;
+
+static int max_dist=0;
+
+typedef boost::adjacency_matrix<boost::directedS> Graph;
+typedef std::pair<int,int> Edge;
+typedef boost::graph_traits<Graph> GraphTraits;
+typedef boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::strict_upper> prob;
+typedef boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::strict_upper> pathlength;
+typedef boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
+
+namespace std {
+	using ::time;
+}
+
+static int var_fixed[5] = {1, 0, 1, 1, 0};
+
+typedef boost::minstd_rand base_generator_type;
+typedef adjacency_list < listS, vecS, directedS,
+no_property, property < edge_weight_t, int > > graph_t;
+typedef graph_traits < graph_t >::vertex_descriptor vertex_descriptor;
+typedef graph_traits < graph_t >::edge_descriptor edge_descriptor;
+
+
+//================================================
+string strDouble(double number)
+{
+    stringstream ss;//create a stringstream
+    ss << number;//add number to the stream
+    return ss.str();//return a string with the contents of the stream
+}
+
+//================================================
+
+double gaussian(double sigma)
+{
+    double GaussNum = 0.0;
+    int NumInSum = 10;
+    for(int i = 0; i < NumInSum; i++)
+    {
+        GaussNum += ((double)rand()/(double)RAND_MAX - 0.5);
+    }
+    GaussNum = GaussNum*sqrt((double)12/(double)NumInSum);
+    
+    
+    return GaussNum;
+    
+}
+
+
+
+//=================================================
+double diffclock(clock_t clock1,clock_t clock2)
+{
+	double diffticks=clock1-clock2;
+	double diffms=(diffticks)/CLOCKS_PER_SEC;
+	return diffms;
+}
+
+//================================================
+//================================================================
+double get_new_x(double x, double dx){
+    
+    double new_x;
+    // boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
+    double r = rand()/(double)(pow(2.,31)-1.);
+    
+    if (r > 0.5){            
+        new_x = x + rand()*dx/(double)(pow(2.,31)-1.);
+    } else {            
+        new_x = x - rand()*dx/(double)(pow(2.,31)-1.);
+    }
+    
+    return new_x;
+    
+}
+
+
+//===============================================   
+string string_wrap(string ins, int mode){
+    
+    std::ostringstream s;
+    
+    switch(mode){
+        case 0:
+            s << "\033[1;29m" << ins << "\033[0m";
+            break;
+        case 1:
+            s << "\033[1;34m" << ins << "\033[0m";
+            break;
+        case 2:
+            s << "\033[1;44m" << ins << "\033[0m";
+            break;
+        case 3:
+            s << "\033[1;35m" << ins << "\033[0m";
+            break;
+        case 4:
+            s << "\033[1;33;44m" << ins << "\033[0m";
+            break;
+        case 5:
+            s << "\033[1;47;34m" << ins << "\033[0m";
+            break;
+        case 6:
+            s << "\033[1;1;31m" << ins << "\033[0m";
+            break;
+        case 7:
+            s << "\033[1;1;33m" << ins << "\033[0m";
+            break;
+        case 8:
+            s << "\033[1;1;43;34m" << ins << "\033[0m";
+            break;
+        case 9:
+            s << "\033[1;1;37m" << ins << "\033[0m";
+            break;
+        case 10:
+            s << "\033[1;30;47m" << ins << "\033[0m";
+            break;
+        default:
+            s << ins;
+    }
+    
+    return s.str();
+}
+
+
+//===============================================
+string wrap_double(double val, int mode){
+    
+    std::ostringstream s;
+    s << string_wrap(strDouble(val),mode);
+    
+    return s.str();
+}
+
+
+
+//===============================================
+const     
+string i2string(int i){
+    
+    std::ostringstream s;
+    s << "worker" 
+    << lexical_cast<std::string>(i);
+    
+    return s.str();
+    
+}
+
+//===============================================
+char* i2char(int i){
+    
+    std::ostringstream s;
+    s << "worker" 
+    << lexical_cast<std::string>(i);
+    
+    char* a=new char[s.str().size()+1];
+    memcpy(a,s.str().c_str(), s.str().size());
+    
+    return a;
+}
+
+//================================================
+class Universe {
+	
+private:
+	
+	double alpha_i;
+	double alpha_m;
+	double beta;
+	double gamma;
+	double delta;
+	
+    double TargetNovelty;
+    double CumulativeRelativeLoss;
+    double CRLsquare;
+    string id;
+    
+    
+	int N_nodes;
+	int M_edges;
+	
+	int N_epochs;
+	int N_steps;
+	int N_repeats;
+	
+	int current_epoch;
+	double current_loss;
+	int current_repeat;
+    double current_novelty;
+	
+	int mode_identify_failed;
+    int verbose_level; // 0 is silent, higher is more
+	
+	double k_max;
+	
+	graph_t Full_g;
+	
+	double **Prob;
+	double **Tried;
+	double **Dist;
+	double **Final;
+    double **EdgeIndex;
+	double *Rank;
+	
+    base_generator_type generator;	
+    boost::uniform_real<> uni_dist;
+    boost::geometric_distribution<double> geo;
+    
+public:
+	
+    
+    
+	//======  Constructor ======
+	Universe(const std::string FileToOpen, int Epochs, int Steps, int Repeats, int identify_failed, double target, const std::string idd)
+	{
+		//typedef array_type2::index index2;
+		
+		
+		std::ifstream inFile;
+        //string line;
+        
+        //-------------------------------
+        
+        base_generator_type gene(42u);
+        generator = gene;
+        generator.seed(static_cast<unsigned int>(std::time(0)));
+        boost::uniform_real<> uni_d(0,1);
+        uni_dist = uni_d;
+        
+        //--------------------------------
+        
+		int i, k;
+        int x, y;
+		Edge* edge_array_mine;
+		int num_arcs_mine, num_nodes_mine;
+		int* weights_mine;
+		
+        TargetNovelty = target;
+        CumulativeRelativeLoss = 0.;
+        CRLsquare = 0.;
+        
+		
+		N_epochs  = Epochs;
+		N_steps   = Steps;
+		N_repeats = Repeats;
+		
+		current_epoch = 0;
+		current_loss = 0.;
+		current_repeat = 0;
+        
+        id = idd;
+        
+        verbose_level = 1;
+		
+		mode_identify_failed = identify_failed;
+        
+		
+		//-------------------------------
+		// The first pass though file with the graph
+		inFile.open(FileToOpen.c_str());
+		if (inFile.fail()) {
+			cout << "Unable to open file";
+			exit(1); // terminate with error
+		}else {
+            
+            if (verbose_level > 2){
+                std::cout <<  " Opened <" << FileToOpen << ">"<<std::endl;
+            }
+        }
+		
+		i=0;
+        std::string line;
+		//while (! inFile.eof() && ! inFile.fail()) {
+        while (1==1) {
+            
+            inFile >> x;
+            inFile >> y;
+            
+            if (verbose_level > 2){
+                std::cout << " x: " << x;
+                std::cout << " y: " << y << std::endl;
+            }
+            
+			if (i==0){
+				N_nodes=x;
+				M_edges=y;	
+                break;
+			}
+			i++;
+            
+			
+		}
+		inFile.close();
+        
+        if (verbose_level == 2){
+            std::cout << N_nodes <<  " nodes, " << M_edges << " edges"<<std::endl;
+        }
+		
+		// k_max is the longest distance possible
+		
+        //k_max = M_edges;
+		k_max = 70;
+        
+		//------------------------------------
+		// Get memory allocated for all class members
+		
+		Prob = allocate_2Dmatrix(N_nodes, N_nodes);
+		Tried = allocate_2Dmatrix(N_nodes, N_nodes);
+		Dist = allocate_2Dmatrix(N_nodes, N_nodes);
+		Final = allocate_2Dmatrix(N_nodes, N_nodes);
+        EdgeIndex = allocate_2Dmatrix(N_nodes, N_nodes);
+		Rank = allocate_1Dmatrix(N_nodes);
+		
+        //The second pass through file with the graph
+        
+		for(int i = 0; i < N_nodes; ++i) {
+			Rank[i]=0.;
+			for(int j = 0; j < N_nodes; ++j) {
+				Final[i][j] = 0.;
+				Prob[i][j]=0.;
+				Dist[i][j]=-1.;
+				Tried[i][j]=0.;
+                EdgeIndex[i][j]=-1;
+			}
+		}
+        
+		
+		// Fill in the final graph -- and we are ready to go!
+        
+	    inFile.open(FileToOpen.c_str());
+		if (!inFile) {
+            std::cout << "Unable to open file";
+			exit(1); // terminate with error
+		}
+		else {
+            
+            if (verbose_level > 2){
+                std::cout <<  " Opened <" << FileToOpen << ">"<<std::endl;
+            }
+        }
+        
+		i=0;  
+		while (inFile >> x && inFile >>y) {
+			if (i > 0) {
+				Final[x][y]=1.;
+				Final[y][x]=1.;
+                
+                
+                if (verbose_level == 2){
+                    std::cout << ".";
+                }
+			}
+			i++;
+			
+		}
+        if (verbose_level == 2){
+            std::cout << std::endl;
+        }
+		inFile.close(); 
+		
+        k=0;
+        for (int i=0; i<N_nodes-1; i++){
+            for (int j=i+1;j<N_nodes; j++){
+                if(Final[i][j] > 0.){
+                    EdgeIndex[i][j]=k;
+                    k++;
+                }
+            }
+        }
+        
+        
+		
+		//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+		// create graph -- hopefully, we can keep it, just modifying edge weights
+		
+		
+		edge_array_mine = new Edge[2*M_edges];
+		num_arcs_mine = 2*M_edges;
+		num_nodes_mine = N_nodes;
+		weights_mine = new int[2*M_edges];
+		for (int i=0; i<2*M_edges; i++){ weights_mine[i]=1;}
+		
+		k=0;
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				if (Final[i][j]>0.){
+					edge_array_mine[2*k]  =Edge(i,j);
+					edge_array_mine[2*k+1]=Edge(j,i);
+					k++;
+				}
+			}
+		}
+		graph_t g(edge_array_mine, edge_array_mine + num_arcs_mine, weights_mine, num_nodes_mine);
+		
+		Full_g = g;
+		delete edge_array_mine;
+		delete weights_mine;
+		
+		//===========================================================================
+		std::vector<edge_descriptor> p(num_edges(Full_g));
+		std::vector<int> d(num_edges(Full_g));
+		edge_descriptor s;
+		boost::graph_traits<graph_t>::vertex_descriptor u, v;
+		
+		for (int i=0; i<N_nodes-1; i++){
+			for (int j=i+1; j<N_nodes; j++){
+				if (Final[i][j] > 0.){
+					u = vertex(i, Full_g);
+					v = vertex(j, Full_g);
+					remove_edge(u,v,Full_g);
+					remove_edge(v,u,Full_g);
+					
+				}
+			}
+		}
+		
+        
+    }
+	
+	
+	//=====================================================================
+	int sample_failed_number(double pfail){
+		
+		//boost::geometric_distribution<double> geo(pfail);
+		//boost::variate_generator<base_generator_type&, geometric_distribution<double> > geom(generator, geo);
+		
+		double r, u, g;
+        
+        r=0.;
+		for(int i=0; i<N_steps; i++){
+            
+            u=(double)rand();
+            u = 1.-u /(double)(pow(2.,31)-1.);
+            g=(int)(ceil(log(u) / log(pfail)));
+            
+			//r += geom();
+            
+            r+=g;
+		}
+        
+        if (verbose_level>=3){
+            std::cout << id << " failed " << r << std::endl;
+		}
+		return r;
+		
+	}
+    
+    //=============================================
+    double get_target(void){
+        return TargetNovelty;
+    }
+	
+    //=============================================
+    void set_target(double target){
+        TargetNovelty=target;
+    }
+	
+	//=============================================
+	int sample(){
+		
+        //boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
+        // double r = uni(), Summa = 0.;
+        
+        
+        
+        double r = rand(), Summa = 0.;
+        r /= (double)(pow(2.,31)-1.);
+		int result = 0;
+		int finished = 0;
+        
+        if (verbose_level==4){
+            std::cout << id << " sampled " << r << std::endl;
+        }
+		
+		for(int i=0; i<N_nodes-1 && finished==0; i++){			
+			for( int j=i+1; j<N_nodes && finished==0; j++){
+				
+				Summa += Prob[i][j];
+				
+				if (Summa > r){
+					
+					Tried[i][j]+=1.;
+					
+					if (Final[i][j] > 0.){
+						result = 1;
+					}
+					finished = 1;										
+				}
+			}
+		}
+		
+		return result;
+		
+	}
+	
+	//===============================
+	void update_current_graph(void){
+		
+		std::vector<edge_descriptor> p(num_edges(Full_g));
+		std::vector<int> d(num_edges(Full_g));
+		edge_descriptor s;
+		boost::graph_traits<graph_t>::vertex_descriptor u, v;
+		
+		//property_map<graph_t, edge_weight_t>::type weightmap = get(edge_weight, Full_g);
+		for (int i=0; i<N_nodes-1; i++){
+			for (int j=i+1; j<N_nodes; j++){
+				if (Final[i][j] > 0. && Tried[i][j]>0){
+					//s = edge(i, j, Full_g);	
+					boost::graph_traits<graph_t>::edge_descriptor e1,e2;
+					bool found1, found2;
+					u = vertex(i, Full_g);
+					v = vertex(j, Full_g);
+					tie(e1, found1) = edge(u, v, Full_g);
+					tie(e2, found2) = edge(v, u, Full_g);
+					if (!found1 && !found2){
+						add_edge(u,v,1,Full_g);
+					    add_edge(v,u,1,Full_g);
+					}
+					
+				}
+			}
+			
+		}
+	}
+	
+	//===============================
+	void update_distances(void){
+		// put shortest paths to the *Dist[][]
+		std::vector<vertex_descriptor> p(num_vertices(Full_g));
+		std::vector<int> d(num_vertices(Full_g));
+		vertex_descriptor s;
+		
+		
+		// put shortest paths to the *Dist[][]
+		for (int j=0; j<num_vertices(Full_g); j++){
+			
+			if(Rank[j] > 0.){
+				s = vertex(j, Full_g);	 
+				dijkstra_shortest_paths(Full_g, s, predecessor_map(&p[0]).distance_map(&d[0]));
+				
+				//std::cout <<" Vertex "<< j << std::endl;
+				graph_traits < graph_t >::vertex_iterator vi, vend;
+				
+				for (boost::tie(vi, vend) = vertices(Full_g); vi != vend; ++vi) {
+					
+					if (p[*vi]!=*vi){
+						Dist[*vi][j]=d[*vi];
+						Dist[j][*vi]=d[*vi];
+                        
+                        if (Dist[*vi][j]>max_dist){
+                            max_dist=Dist[*vi][j];
+                        }
+                        
+                        
+					} else {
+						Dist[*vi][j]=-1.;
+						Dist[j][*vi]=-1.;
+					}
+				}
+			}
+			
+		}
+		
+		
+	}
+	
+	//======================================================
+	void update_ranks(void){
+		
+		for(int i=0; i<N_nodes; i++){
+			Rank[i]=0.;
+		}
+		
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				if (Tried[i][j]>0. && Final[i][j] >0.){
+					Rank[i]++;
+					Rank[j]++;
+				}
+			}
+		}
+		
+	}
+	
+	//====================================================================
+	void set_world(double a_i, double a_m, double b, double g, double d){
+		
+		alpha_i=a_i;
+		alpha_m=a_m;
+		gamma=g;
+		beta=b;
+		delta=d;
+		
+	}
+	
+	//====================================================================
+	void reset_world(){
+		
+        //====================================================
+		std::vector<edge_descriptor> p(num_edges(Full_g));
+		std::vector<int> d(num_edges(Full_g));
+		edge_descriptor s;
+		boost::graph_traits<graph_t>::vertex_descriptor u, v;
+        
+		
+		for (int i=0; i<N_nodes-1; i++){
+			for (int j=i+1; j<N_nodes; j++){
+				if (Final[i][j] > 0. && Tried[i][j] > 0){
+					u = vertex(i, Full_g);
+					v = vertex(j, Full_g);
+					remove_edge(u,v,Full_g);
+					remove_edge(v,u,Full_g);
+					
+				}
+			}
+		}
+        
+        //==================================================
+        
+		current_loss=0;
+		current_epoch=0;
+		current_repeat++;
+        current_novelty=0;
+		
+		for(int i = 0; i < N_nodes; ++i) {
+			Rank[i]=0.;
+			for(int j = 0; j < N_nodes; ++j) {
+				Prob[i][j]=0.;
+				Dist[i][j]=-1.;
+				Tried[i][j]=0.;
+			}
+		}
+	}
+	
+	
+    //==============================================
+    void show_parameters(void){
+        
+        std::cout << "Parameters: " 
+        << alpha_i << " "
+        << alpha_m << " | "
+        << beta << " "
+        << gamma << " | "
+        << delta << std::endl;
+        
+    }
+    
+    
+    
+    //===============================================
+    string file_name(){
+        
+        std::ostringstream s;
+        s << "world_" 
+        << lexical_cast<std::string>(alpha_i) << "_" 
+        << lexical_cast<std::string>(alpha_m) << "_"
+        << lexical_cast<std::string>(beta) << "_"
+        << lexical_cast<std::string>(gamma) << "_"
+        << lexical_cast<std::string>(delta) << "_"
+        << lexical_cast<std::string>(N_epochs) << "_"
+        << lexical_cast<std::string>(N_steps) << "_"
+        << lexical_cast<std::string>(N_repeats) << ".txt";
+        
+        return s.str();
+        
+    }
+    
+    
+    
+    
+    //=================================================
+    void set_verbose(int verbose){
+        
+        verbose_level = verbose;
+    }
+    
+    
+    //=============================================================
+    void update_probabilities(void){
+        
+        
+        //=========================
+		// Compute sampling probabilities
+		// first pass: \xi_i,j
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				
+				double bg = 0.;
+				
+				Prob[i][j] = alpha_i*log(min(Rank[i]+1.,Rank[j]+1.)) + 
+                alpha_m*log(max(Rank[i]+1.,Rank[j]+1.));
+				
+                if (Dist[i][j] > 0.){
+                    
+                    double k = Dist[i][j];
+                    if (k >= k_max){
+                        k = k_max-1;
+                    }
+					
+                    bg = beta * log(k/k_max) + gamma * log(1. - k/k_max);
+					
+                } else {
+                    bg = delta;
+                }
+				
+				Prob[i][j] = exp(Prob[i][j] + bg);
+			}
+		}
+        
+		
+		// second pass: sum
+		double Summa = 0.;
+		
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				Summa += Prob[i][j];
+			}
+		}
+		
+		// third pass: normalize
+		for(int i=0; i<N_nodes-1; i++){
+			for( int j=i+1; j<N_nodes; j++){
+				Prob[i][j] /= Summa;
+			}
+		}
+        
+    }
+    
+	// Now we are ready for simulations
+	//==============================================
+	void update_world(){
+		
+		int failed = 0;
+        
+		// Given current universe compute shortest paths
+		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+		
+		update_current_graph();
+		update_ranks();				
+		update_distances();
+		update_probabilities();
+		
+		//===============================
+		// sampling
+		int result;
+		double cost=0., novel=0.;
+		int publishable = 0;
+		
+		
+		//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+		if (mode_identify_failed == 1){
+			
+			while(publishable < N_steps){
+				
+		    	result = sample();
+			    publishable += result;
+			    failed += (1-result);
+				
+			}
+			
+			for(int i=0; i<N_nodes-1; i++){
+				for( int j=i+1; j<N_nodes; j++){
+					
+					cost+=Tried[i][j];
+					
+					if (Tried[i][j]>0. && Final[i][j]>0.){
+						novel+=1.;
+					}
+				}
+			}
+			
+		}
+		//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+		else {
+			
+			double pfail=0.;
+			int n_failed;
+			//, n_check = 0;
+			
+			for(int i=0; i<N_nodes-1; i++){
+				for( int j=i+1; j<N_nodes; j++){
+					if (Final[i][j] == 0.){
+						pfail += Prob[i][j];
+						Prob[i][j] = 0.;
+					}
+					
+				}
+			}
+			
+			for(int i=0; i<N_nodes-1; i++){
+				for( int j=i+1; j<N_nodes; j++){
+					Prob[i][j] /= (1.-pfail);
+				}
+				//std::cout << std::endl;
+			}			
+			
+			n_failed = sample_failed_number(pfail);
+			while(publishable < N_steps){
+				
+		    	result = sample();
+			    publishable += result;					
+			}
+            
+            
+			current_loss += (n_failed + N_steps);
+			cost = current_loss;
+			
+			for(int i=0; i<N_nodes-1; i++){
+				for( int j=i+1; j<N_nodes; j++){
+					
+					if (Tried[i][j]>0. && Final[i][j]>0.){
+						novel+=1.;
+					}
+				}
+			}
+		}
+        
+        current_novelty = novel;
+        
+        
+		//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+		if (verbose_level == 2){
+            std::cout << (current_repeat+1) << "  epoch=" << (current_epoch+1) 
+            
+		    << "  cost=" << cost 
+		    << " novel=" << novel 
+		    << " rel_loss=" << cost/novel
+		    << std::endl;
+        }
+        
+		current_epoch++;
+	}
+	
+	
+	//======  Destructor ======
+	~Universe(){
+		
+		delete_2Dmatrix(Final, N_nodes);
+		delete_2Dmatrix(Dist, N_nodes);
+		delete_2Dmatrix(Tried, N_nodes);
+		delete_2Dmatrix(Prob, N_nodes);
+        delete_2Dmatrix(EdgeIndex, N_nodes);
+		delete_1Dmatrix(Rank);		
+	}
+    
+    //================================================
+    // Allocate memory
+    double** allocate_2Dmatrix(int N, int M)
+    {
+        double **pointer;
+        
+        if (verbose_level == 2){
+            std::cout<< "["<<N<<"|"<<M<<"]"<<std::endl;
+        }
+        pointer = new double*[N];
+        for (int i = 0; i < N; ++i)
+            pointer[i] = new double[M];	
+        
+        return pointer;
+    }
+    //===================
+    double* allocate_1Dmatrix(int N)
+    {
+        double *pointer;
+        
+        if(N > 0){
+            
+            pointer = new double[N];
+            
+        }else {
+            
+            pointer = NULL;
+        }
+        
+        return pointer;
+        
+    }
+    
+    //==============================================
+    // De-Allocate memory to prevent memory leak
+    void delete_2Dmatrix(double **pointer, int N){
+        
+        if (pointer != NULL){
+            
+            for (int i = 0; i < N; ++i){
+                delete [] pointer[i];
+            }
+            delete [] pointer;
+        }
+    }
+    //====================
+    void delete_1Dmatrix(double *pointer){
+        
+        delete [] pointer;
+    }
+    
+    //===========================================
+    double get_rel_loss(){
+        
+        return CumulativeRelativeLoss ;
+    }
+    
+    //===========================================
+    double get_rel_loss_err(){
+        
+        return CRLsquare ;
+    }
+    
+    
+    
+    //==================================================================================
+    void evolve_to_target_and_save(int istart, int iend, double* storage, int* counters){
+        
+        double ALOT=100000000000.;
+        
+	// std::cout<<" evolve_to_target_and_save: istart=" << istart << "iend=" << iend << "\n";
+
+        reset_world();
+        
+        for (int k = istart; k < iend; k++){
+
+	  // std::cout<<" evolve: k=" << k << "\n";
+
+            
+            
+            for(int i=0; i< N_epochs &&  current_novelty < TargetNovelty; i++){
+	      // std::cout<<" evolve: k=" << k << " i=" << i << " cur=" << current_novelty << " Target=" << TargetNovelty << "\n";
+                update_world();
+            }
+            
+            storage[k]=current_loss/current_novelty;
+            counters[k]=1;
+            
+            
+            reset_world();
+        }
+        
+    }
+    //==============================================
+    int get_reruns(void){
+        return N_repeats;
+    }
+    
+    //==============================================
+    double get_parameter(int i){
+        
+        switch(i){
+            case 0:
+                return alpha_i;
+            case 1:
+                return alpha_m;
+            case 2:
+                return beta;
+            case 3:
+                return gamma;
+            case 4:
+                return delta;
+            default:
+                
+                std::cout << "Erroneous parameter id!!!!\n\n\n";
+                return 0.;
+        }
+    }
+    
+    
+    //==============================================
+    void evolve_to_target(){
+        
+        reset_world();
+        if (beta < -1. || gamma < -1.){
+            CumulativeRelativeLoss = 100000000000.;
+            CRLsquare = 0.;
+            return;
+        }
+        
+        
+        for (int k=0; k< N_repeats; k++){
+            
+            
+            for(int i=0; i<N_epochs &&  current_novelty < TargetNovelty; i++){
+                update_world();
+            }
+            
+            CumulativeRelativeLoss += current_loss/current_novelty;
+            CRLsquare += (current_loss/current_novelty)*(current_loss/current_novelty);
+            if(verbose_level==3){
+                std::cout <<  CumulativeRelativeLoss << " | " << CRLsquare << std::endl;
+            }
+            
+            if(verbose_level==1){
+                std::cout <<  "." ;
+                
+            }
+            else if(verbose_level==2){
+                std::cout <<  "**" << (k+1) <<  "**  curr loss " << current_loss << "; curr novelty " << current_novelty << std::endl;
+            }
+            
+            
+            reset_world();
+        }
+        
+        CumulativeRelativeLoss /= double(N_repeats);
+        CRLsquare /= double(N_repeats);
+        
+        if(verbose_level==1){
+            std::cout << std::endl;
+        }
+        
+        if(verbose_level==2){
+            std::cout <<  CumulativeRelativeLoss << " || " << CRLsquare << std::endl;
+        }
+        
+        CRLsquare = 2*sqrt((CRLsquare - CumulativeRelativeLoss*CumulativeRelativeLoss)/double(N_repeats));
+        
+    }
+    
+    
+    //================================================================    
+    int set_parameter(double value, int position){
+        
+        if (position < 0 || position > 4) {return 0;}
+        
+        else {
+            
+            switch(position){
+                case 0:
+                    alpha_i=value;
+                    return 1;
+                case 1:
+                    alpha_m=value;
+                    return 1;
+                case 2:
+                    beta=value;
+                    return 1;
+                case 3:
+                    gamma=value;
+                    return 1;
+                case 4:
+                    delta=value;
+                    return 1;
+            }
+            
+        }
+        
+        return 0;
+    }
+    
+    
+    //=================================================================
+    void try_annealing(double starting_jump, int iterations, 
+                       double temp_start, double temp_end, double target_rejection){
+        
+        double dx[5]={0.,0.,0.,0.,0};
+        double x[5]={0.,0.,0.,0.,0};
+        double rejection[5]={0., 0., 0., 0., 0.};
+        double curr_x, curr_err, x_tmp;
+        double temperature;
+        double ratio, r;
+        int cycle=10;
+        boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
+        
+        // set up parameter for annealing
+        
+        x[0]=alpha_i;
+        x[1]=alpha_m;
+        x[2]=beta;
+        x[3]=gamma;
+        x[4]=delta;
+        
+        for(int i=0;i<5;i++){
+            dx[i] = starting_jump;
+        }
+        
+        // establish the current value
+        
+        //..........................................
+        evolve_to_target();        
+        std::cout << CumulativeRelativeLoss << " +- " << CRLsquare << std::endl;
+        
+        curr_x   = CumulativeRelativeLoss;
+        curr_err = CRLsquare;
+        CumulativeRelativeLoss = 0;
+        CRLsquare = 0;
+        //...........................................
+        
+        // optimization cycle
+        for(int i=0; i<iterations; i++){
+            
+            temperature = temp_start*exp( i*(log(temp_end)-log(temp_start))/(double)iterations);
+            std::cout  << std::endl << "....T = " << wrap_double(temperature,3) << std::endl << std::endl;
+            
+            if (i % cycle == 0 && i > 0){
+                
+                for (int k=0; k<5; k++){
+                    
+                    rejection[k]/=(double)cycle;
+                    if (rejection[k] > 0){
+                        dx[k] = dx[k]/(rejection[k]/target_rejection);
+                        rejection[k]=0.;
+                    }
+                    else{
+                        dx[k]*=2.;
+                    }
+                    std::cout  << dx[k] << " ";
+                }
+                std::cout  << std::endl;
+            }
+            
+            
+            for (int j=0; j<5; j++){
+                
+                // get new value of x[j]
+                x_tmp = get_new_x(x[j],dx[j]);
+                
+                
+                
+                //.............................................
+                set_parameter(x_tmp, j);
+                
+                
+                evolve_to_target(); 
+                
+                std::cout  << std::endl << "......... " << std::endl;
+                std::cout << "Trying... " << CumulativeRelativeLoss << " +- " << CRLsquare << std::endl;
+                
+                ratio = min(1.,exp(-(CumulativeRelativeLoss-curr_x)/temperature));
+                r = uni();
+                std::cout << r << " vs " << ratio << std::endl;
+                
+                if (r > ratio){
+                    
+                    std::cout << string_wrap(id, 4) <<" "<< (i+1) << ","<< (j) 
+                    <<" "<< (i+1) << " Did not accept " 
+                    << x_tmp << "(" << j << ")" << std::endl;
+                    std::cout << alpha_i << " "<< alpha_m << " "
+                    << beta << " " << gamma << " " 
+                    << delta << " " << std::endl;
+                    set_parameter(x[j], j);
+                    CumulativeRelativeLoss = 0;
+                    CRLsquare = 0;
+                    
+                    rejection[j]+=1.;
+                }
+                
+                else {
+                    
+                    curr_x   = CumulativeRelativeLoss;
+                    curr_err = CRLsquare;
+                    x[j] = x_tmp;
+                    CumulativeRelativeLoss = 0;
+                    CRLsquare = 0;
+                    std::cout << (i+1) << string_wrap((string) " Rejection counts: ", 8) 
+                    << wrap_double(rejection[0],2) 
+                    << " "<< wrap_double(rejection[1], 7) << " "
+                    << wrap_double(rejection[2],5) << " " << wrap_double(rejection[2],9) << " " 
+                    << wrap_double(rejection[4],6) << " " 
+                    << std::endl << std::endl;
+                    
+                    std::cout << string_wrap(id, 4) <<" "<< (i+1) <<","<< (j) 
+                    <<" "
+                    << string_wrap((string) "***** Did accept! ", 3) 
+                    << wrap_double(alpha_i,2) 
+                    << " "<< wrap_double(alpha_m, 7) << " "
+                    << wrap_double(beta,5) << " " 
+                    << wrap_double(gamma,9) << " " 
+                    << wrap_double(delta,6) << " " 
+                    << std::endl << std::endl;
+                    
+                }
+                //........................................................ 
+                
+            }
+            
+        }
+        
+    }
+    
+	
+};
+
+//============================================================
+
+std::pair<double,double> multi_loss( // dispatch_group_t group, 
+                                    Universe* un[], 
+                                    // dispatch_queue_t* CustomQueues,
+                                    double* Results,
+                                    int*    Counters,
+                                    double* params){
+    
+    int N = un[0]->get_reruns();
+    int step = (int)(double)N/(double)(Nworkers);
+    int istart=0;
+    int iend = istart+step;
+    
+    double Loss=0., LossSquare=0.;
+    
+    timeval startTime, endTime;
+    double elapsedTime;
+    // start timer
+    gettimeofday(&startTime, NULL);
+
+    //err:    for(int i=0; i<Nworkers; i++){
+        
+    for(int i=0; i<Nworkers; i++){
+        for(int j=0; j<5; j++){
+            un[i]->set_parameter(params[j],j);
+        }
+    }        
+    int i;
+    #pragma omp parallel for private (i)
+    for(i=0; i<Nworkers; i++){
+
+      // dispatch_group_async(group, CustomQueues[i], ^{
+      //      std::cout<<"multi_loss: Calling evolve_to_target_and_save i=" << i << " N=" << N << " istart=" << i*step << " iend=" << (i+1)*step << "\n";
+      //un[i]->evolve_to_target_and_save(istart, iend, Results, Counters);
+      un[i]->evolve_to_target_and_save(i*step, min((i+1)*step,N), Results, Counters);
+      //});
+
+      // std::cout<<"multi_loss: Returned from evolve_to_target_and_save " << i << "\n";
+      
+
+      //istart += step;
+      //iend = min(istart+step,N);
+            
+    }
+    // err    }
+    // dispatch_group_wait(group, DISPATCH_TIME_FOREVER);
+    //dispatch_release(group);
+    
+    for (int i=0; i<N; i++){
+        
+        Loss+=Results[i]/(double)N;
+        LossSquare+=Results[i]*Results[i]/(double)N;
+        
+        std::cout<<i<<":"<< Results[i] << " ";
+    }
+    
+    std::cout<<" \n\n\n";
+    double two_std = ((LossSquare - Loss*Loss)/(double)N);
+    
+    two_std = 2.*sqrt(two_std);
+    std::pair<double,double> Res;
+    Res.first=Loss;
+    Res.second=two_std;
+
+    gettimeofday(&endTime, NULL);
+    elapsedTime = (endTime.tv_sec - startTime.tv_sec) * 1000.0;      // sec to ms
+    elapsedTime += (endTime.tv_usec - startTime.tv_usec) / 1000.0;   // us to ms
+    elapsedTime /= 1000.;
+    cout << "multi_loss(N=" << N << ") elapsed time: " << elapsedTime << " seconds " << elapsedTime/60. << " minutes\n\n";
+    
+    return Res;
+}
+//============================================================
+
+
+//============================================================
+void multi_annealing( // dispatch_group_t group, 
+                     Universe* un[], 
+                     // dispatch_queue_t* CustomQueues, 
+                     double T_start, double T_end, 
+                     double Target_rejection, 
+                     int Annealing_repeats, 
+                     double starting_jump,
+                     double* Results,
+                     int*    Counters,
+                     double* params0,
+                     double annealing_cycles){
+    //.................................
+    // re-implement annealing
+    
+    double dx[5]={0.,0.,0.,0.,0};
+    double x[5]={0.,0.,0.,0.,0};
+    double rejection[5]={0., 0., 0., 0., 0.};
+    double curr_x, curr_err, x_tmp;
+    double temperature;
+    double ratio, r;
+    int cycle=10;
+    //boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
+    
+    // set up parameter for annealing
+    
+    x[0]=params0[0];
+    x[1]=params0[1];
+    x[2]=params0[2];
+    x[3]=params0[3];
+    x[4]=params0[4];
+    
+    for(int i=0;i<5;i++){
+        dx[i] = starting_jump;
+    }
+    
+    // establish the current value
+    std::pair<double,double>Res;
+    
+    Res = multi_loss( /* group,*/ un, /* CustomQueues,*/ Results, Counters, x);
+    std::cout << Res.first << " +- " << Res.second << std::endl;
+    
+    curr_x   = Res.first;
+    curr_err = Res.second;
+    
+    // optimization cycle
+    
+    for(int i=0; i<annealing_cycles; i++){
+        
+        temperature = T_start*exp( i*(log(T_end)-log(T_start))/(double)annealing_cycles);
+        std::cout  << std::endl << "....T = " << wrap_double(temperature,3) << std::endl << std::endl;
+        
+        if (i % cycle == 0 && i > 0){
+            
+            for (int k=0; k<5; k++){
+                rejection[k]/=(double)cycle;
+                
+                if (rejection[k] > 0){
+                    dx[k] = dx[k]/(rejection[k]/Target_rejection);
+                    rejection[k]=0.;
+                }
+                else{
+                    dx[k]*=2.;
+                }
+                std::cout  << dx[k] << " ";
+            }
+            std::cout  << std::endl;
+        }
+        
+        
+        for (int j=0; j<5; j++){
+            
+            ///////////////////////////////
+            if (FIX_VARIABLES==0 || var_fixed[j]==0){
+                
+                
+                
+                // get new value of x[j]
+                double x_hold=x[j];
+                x_tmp = get_new_x(x[j],dx[j]);
+                x[j]=x_tmp;
+                
+                std::cout << wrap_double(x_tmp,10) << " " << wrap_double(j,9) << "\n\n"; 
+                //=======================================
+                //.............................................
+                for(int w=0; w<Nworkers; w++){
+                    un[w]->set_parameter(x_tmp, j);
+                }
+                
+                
+                Res = multi_loss(/* group, */ un, /* CustomQueues, */ Results, Counters, x);
+                std::cout << Res.first << " +- " << Res.second << std::endl;
+                
+                ratio = min(1.,exp(-(Res.first-curr_x)/temperature));
+                r = rand()/(double)(pow(2.,31)-1.);
+                std::cout << r << " vs " << ratio << std::endl;
+                
+                double ALOT=100000000000.;
+                
+                if (Res.first < ALOT)
+                {
+                    ofstream filestr;
+                    
+                    filestr.open ("best_opt_some.txt", ofstream::app);
+                    
+                    // >> i/o operations here <<
+                    filestr << un[0]->get_target() << "," 
+                    << Res.first 
+                    << "," << un[0]->get_parameter(0) 
+                    << "," << un[0]->get_parameter(1) 
+                    << "," << un[0]->get_parameter(2) 
+                    << "," << un[0]->get_parameter(3) 
+                    << "," << un[0]->get_parameter(4) << "," << Res.second << ",\n";
+                    
+                    filestr.close();
+                    
+                    
+                    filestr.open ("max_dist.txt", ofstream::app);
+                    
+                    // >> i/o operations here <<
+                    filestr << max_dist << ",\n";
+                    
+                    filestr.close();
+                    
+                }
+                
+                
+                if (r > ratio){
+                    
+                    std::cout << " "<< (i+1) << ","<< (j) 
+                    <<" "<< (i+1) << " Did not accept " 
+                    << x_tmp << "(" << j << ")" << std::endl;
+                    std::cout << un[0]->get_parameter(0) 
+                    << " " << un[0]->get_parameter(1) 
+                    << " " << un[0]->get_parameter(2) 
+                    << " " << un[0]->get_parameter(3) 
+                    << " " << un[0]->get_parameter(4) << " " << std::endl;
+                    
+                    x[j]=x_hold;
+                    for(int w=0; w<Nworkers; w++){
+                        un[w]->set_parameter(x[j], j);
+                    }
+                    
+                    
+                    //set_parameter(x[j], j);     
+                    rejection[j]+=1.;
+                }
+                
+                else {
+                    
+                    curr_x   = Res.first;
+                    curr_err = Res.second;
+                    x[j] = x_tmp;
+                    
+                    for(int w=0; w<Nworkers; w++){
+                        un[w]->set_parameter(x[j], j);
+                    }
+                    
+                    std::cout << (i+1) << string_wrap((string) " Rejection counts: ", 8) 
+                    << wrap_double(rejection[0],2) << " " 
+                    << wrap_double(rejection[1],7) << " "
+                    << wrap_double(rejection[2],5) << " " 
+                    << wrap_double(rejection[3],9) << " " 
+                    << wrap_double(rejection[4],6) << " " 
+                    << std::endl << std::endl;
+                    
+                    std::cout << " "<< (i+1) <<","<< (j) 
+                    <<" "
+                    << string_wrap((string) "***** Did accept! ", 3) 
+                    << wrap_double(un[0]->get_parameter(0),2) << " "
+                    << wrap_double(un[0]->get_parameter(1),7) << " "
+                    << wrap_double(un[0]->get_parameter(2),5) << " " 
+                    << wrap_double(un[0]->get_parameter(3),9) << " " 
+                    << wrap_double(un[0]->get_parameter(4),6) << " " 
+                    << std::endl << std::endl;
+                    
+                    
+                    
+                }
+                //........................................................ 
+                
+            }
+        }
+        
+    }
+    
+}
+
+
+
+//================================================
+int
+main(int argc, char* argv[])
+{
+    
+    double params0[6] = {0., 0., 0., 0., 0., 0.2}, target=50., range;
+    string par_names0[6] = {"alpha_i", "alpha_m", "beta", "gamma", "delta", "target"};
+    string par_names1[4] = {"n_epochs", "n_steps", "n_reruns", "range"};
+    string par_names2[5] = {"T_start", "T_end", "Annealing_steps","Target_rejection","Starting_jump"};
+    string par_names3[5] = {"FREEZE_alpha_i", "FREEZE_alpha_m", "FREEZE_beta", "FREEZE_gamma", "FREEZE_delta"};    
+    int params1[4] = {300, 50, 1000, 10};
+    int params3[5] = { 0, 0, 0, 0, 0};
+    
+    //          temperature_start,  temperature_end,  annealing_steps target_rejection  Starting_jump
+    double params2[5] = {1,             0.001,               100,              0.3,           1.5};
+    
+    int verbose_level = 2;
+    const std::string one="one", two="two";
+    static Universe* un[Nworkers];
+    // static dispatch_queue_t CustomQueues[Nworkers];
+    
+    static double* Results;
+    static int*    Counters;
+    
+    timeval t1, t2;
+    double elapsedTime;
+    // start timer
+    gettimeofday(&t1, NULL);
+    
+    
+    if (argc < 8) {
+        std::cout << "Usage: super_optimizer alpha_i alpha_m beta gamma delta target_innov [n_epochs n_steps n_reruns] [range] [verbose_level]\n"; 
+        std::cout << "         [T_start T_end Annealing_steps Target_rejection Starting_jump]\n"; 
+        std::cout << "         [FREEZE_alpha_i FREEZE_alpha_m FREEZE_beta FREEZE_gamma FREEZE_delta]\n"; 
+        
+        system("pwd");
+        
+        
+        return(1);
+    }
+    else {
+        for (int nArg=0; nArg < argc; nArg++){
+            //std::cout << nArg << " " << argv[nArg] << std::endl;
+            if (nArg > 0 && nArg < 7){
+                params0[nArg-1]= atof(argv[nArg]);
+                std::cout << par_names0[nArg-1] << ": " << params0[nArg-1] <<  std::endl;
+            }
+            if (nArg > 6 && nArg < 11){
+                params1[nArg-7]= atoi(argv[nArg]);
+                std::cout << par_names1[nArg-7] << ": " << params1[nArg-7] <<  std::endl;
+            }
+            if (nArg == 11){
+                verbose_level = atoi(argv[nArg]);
+                std::cout << "verbose level: " << verbose_level <<  std::endl;
+            }
+            if (nArg > 11 && nArg < 17){
+                params2[nArg-12]= atof(argv[nArg]);
+                std::cout << par_names2[nArg-12] << ": " << params2[nArg-12] <<  std::endl;
+            }
+            if (nArg > 16 && nArg < 22){
+                params3[nArg-17]= atof(argv[nArg]);
+                var_fixed[nArg-17]= atof(argv[nArg]);
+                std::cout << par_names3[nArg-17] << ": " << var_fixed[nArg-17] <<  std::endl;
+            }
+            
+            
+        }
+        
+    }
+    
+    /*
+    for target in range(58,1009,50):
+	s = ("%d" % target)
+	print s
+
+	for i in range(15):
+          # Param groups separated by "|" below for documentation. NOTE that | is not used on command line!
+	  os.system("./supe_duper_optimizer |0 0 4 50 -1 "+s+" | 40000 20 1000 2 | 1 | 2. 0.01 100 0.3 2.3 | 1 1 0 0 0")
+    */
+
+    /*   Parameters, re-iterated:
+
+                        "alpha_i", "alpha_m", "beta", "gamma", "delta", "target"};
+    double params0[6] = {0.,        0.,        0.,     0.,      0.,      0.2}, target=50., range;
+
+                    {"n_epochs", "n_steps", "n_reruns", "range"};
+    int params1[4] = {300,        50,       1000,        10};
+
+                       {"T_start",          "T_end",          "Annealing_steps","Target_rejection","Starting_jump"};
+                         temperature_start,  temperature_end,  annealing_steps target_rejection  Starting_jump
+    double params2[5] = {1,                  0.001,            100,            0.3,              1.5};
+
+                     {"FREEZE_alpha_i", "FREEZE_alpha_m", "FREEZE_beta", "FREEZE_gamma", "FREEZE_delta"};    
+    int params3[5] = { 0,                0,                0,             0,              0};
+
+
+    */
+
+    for (int j=0; j<5; j++){
+        
+        cout << j << " | " << var_fixed[j] << " (fixed) \n";
+    }
+	
+    target=params0[5];
+    range = (double)params1[3];
+	int identify_failed = 0;
+	char* filename= (char *)"movie_graph.txt";
+ 	int n_ep=params1[0], n_st=params1[1], n_rep=params1[2];
+    
+    //...............................
+    
+    for(int i=0; i<Nworkers; i++){
+        un[i] = new Universe((char *)filename,n_ep,n_st,
+                             (int)n_rep,
+                             identify_failed, target, i2string(i));
+        // CustomQueues[i] = dispatch_queue_create(i2char(i), NULL);
+    }
+    
+    //...............................
+    if(n_rep > 0){
+        
+        Results = new double[n_rep];
+        Counters = new int[n_rep];
+        
+    }else {
+        
+        Results =  NULL;
+        Counters = NULL;
+        std::cout << " Number of reruns should be positive! " <<  std::endl;
+        return 0;
+        
+    }
+    //...............................
+    srand(time(0));
+    //srandomdev();
+    
+    {
+        double r=0;
+        for (int j=0; j<100; j++){
+            
+            
+            
+            r = rand()/(double)(pow(2.,31)-1.);
+            std::cout << r << " ";
+        }
+        std::cout << "\n ";
+    }
+  	//random initiation of starting parameters
+    
+    if (range > 0.){
+        
+        for (int i=0; i < 5; i++){
+            
+            if (params0[i]==-100.){
+                
+                double r1 = (rand()/(double)(pow(2.,31)-1.));
+                double r2 = (rand()/(double)(pow(2.,31)-1.));
+                double sign = 1.;
+                
+                if(r1 > 0.5){
+                    sign=-1.;
+                }
+                
+                params0[i] = sign*r2*range;
+                
+                std::cout << par_names0[i] << ": " << params0[i] <<  std::endl;
+            }
+        }
+        
+    }
+    
+    
+    double T_start=params2[0], T_end=params2[1], Target_rejection=params2[3], starting_jump=params2[4];
+    int Annealing_repeats = (int) params2[2];
+    
+    
+    // dispatch_group_t group = dispatch_group_create();
+    
+    //.............................
+    multi_annealing( /* group, */ un, /* CustomQueues, */ T_start, T_end, Target_rejection, Annealing_repeats, 
+                    starting_jump, Results, Counters, params0, Annealing_repeats);
+    
+    //dispatch_group_wait(group, DISPATCH_TIME_FOREVER);
+    // dispatch_release(group);
+    //.............................
+    
+    
+    // stop timer
+    gettimeofday(&t2, NULL);
+    
+    // compute and print the elapsed time in millisec
+    elapsedTime = (t2.tv_sec - t1.tv_sec) * 1000.0;      // sec to ms
+    elapsedTime += (t2.tv_usec - t1.tv_usec) / 1000.0;   // us to ms
+    elapsedTime /= 1000.;
+    cout << elapsedTime << " seconds " << elapsedTime/60. << " minutes\n\n";
+    
+    //.....................
+    
+    for(int i=0; i<Nworkers; i++){
+        delete un[i];
+    }
+    
+    //....................
+    if(n_rep > 0){
+        
+        delete [] Results;
+        delete [] Counters;
+        
+    }
+    
+    return 0;
+	
+	
+	
+}
+


Property changes on: SwiftApps/optimizer.snap01.cpp
___________________________________________________________________
Added: svn:executable
   + 

Added: SwiftApps/optirun.swift
===================================================================
--- SwiftApps/optirun.swift	                        (rev 0)
+++ SwiftApps/optirun.swift	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,41 @@
+type file;
+
+app (file outfile, file best, file max) optimize ( string args[], file graph )
+{
+  optimizersh @best @max args stdout=@outfile ;
+}
+
+int minrange=58;
+int maxrange=1009;
+//int maxrange=209;
+int rangeinc=50;
+
+int nreps=2; # 15
+
+//      [alpha_i alpha_m beta gamma delta target_innov
+//      [n_epochs n_steps n_reruns] [range]
+//      [verbose_level]
+//      [T_start T_end Annealing_steps Target_rejection Starting_jump]
+//      [FREEZE_alpha_i FREEZE_alpha_m FREEZE_beta FREEZE_gamma FREEZE_delta]
+
+file graph <"movie_graph.txt">;
+
+foreach target in [minrange:maxrange:rangeinc] {
+  foreach rep in [1:nreps] {
+    file outfile <single_file_mapper; file=@strcat("output/T",target,".R",rep,".out")>;
+    // file errfile <single_file_mapper; file=@strcat("output/T",target,".R",rep,".err")>;
+    file bestfile <single_file_mapper; file=@strcat("output/T",target,".R",rep,".best_opt_some")>;
+    file maxfile <single_file_mapper; file=@strcat("output/T",target,".R",rep,".max_dist")>;
+
+    // string longargs[] = @strcat("0 0 4 50 -1 ",target," 40000 20 1000 2 1 2. 0.01 100 0.3 2.3 1 1 0 0 0");
+
+    string fastargs1[] = ["0", "0", "4", "50", "-1", @strcat(target), "40000", "20", "1000", "2", "1", "2.", "0.01", "100", "0.3", "2.3", "1", "1", "0", "0", "0"];
+    string fastargs2[] = [
+       "0", "0", "4", "50", "-1", @strcat(target),
+       "40000", "20", "1000", "2",
+       "1",
+       "2.", "0.01",  "5", "0.3", "2.3",
+       "1", "1", "0", "0", "0"];
+    (outfile, bestfile, maxfile) = optimize(fastargs2,graph);
+  }
+}

Added: SwiftApps/script-smaller-a.py
===================================================================
--- SwiftApps/script-smaller-a.py	                        (rev 0)
+++ SwiftApps/script-smaller-a.py	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,19 @@
+#! /usr/bin/env python
+import os
+
+# FULL for target in range(58,1009,50):
+# FAST for target in range(58,209,50):
+for target in range(58,209,50):
+	s = ("%d" % target)
+	print s
+
+	for i in range(15):
+# FAST	for i in range(2):
+		args="./optimizer 0 0 4 50 -1 "+s+" 40000 20 1000 2 1 2. 0.01 100 0.3 2.3 1 1 0 0 0 # > out.T"+str(target)+".i"+str(i)
+		print("\n\n **** CALLING APP: "+args+"\n\n\n")
+	        os.system(args);
+#		print("\n\n **** CALLING APP: ./optimizer 0 0 4 50 -1 "+s+" 40000 20 10 2 1 2. 0.01 100 0.3 2.3 1 1 0 0 0\n\n\n")
+#		FAST: os.system("./optimizer 0 0 4 50 -1 "+s+" 40000 20 10 2 1 2. 0.01 2 0.3 2.3 1 1 0 0 0")
+
+
+print "Done!"


Property changes on: SwiftApps/script-smaller-a.py
___________________________________________________________________
Added: svn:executable
   + 

Added: SwiftApps/sumloss.sh
===================================================================
--- SwiftApps/sumloss.sh	                        (rev 0)
+++ SwiftApps/sumloss.sh	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,31 @@
+#! /bin/awk -f
+
+BEGIN { n = ARGC-1 }
+{
+  loss += $1/n
+  loss_sq += ($1*$1)/n
+}
+END { 
+  x = (loss_sq - (loss*loss))/n
+  sdev = 2.0 * sqrt(x)
+  printf "loss sdev\n"
+  printf "%f %f\n", loss, sdev
+}
+
+# the awk script above implements this c++ logic from optimizer.cpp:
+#
+#    for (int i=0; i<N; i++){
+#        double n = double(N);
+#        Loss += Results[i]/n;
+#        LossSquare += Results[i]*Results[i]/n;
+#
+#        std::cout<<" RESULT: " << Results[i];
+#    }
+#
+#    std::cout<<" \n\n\n";
+#    double x = ((LossSquare - Loss*Loss)/(double)N);
+#    double two_std = 2.*sqrt(x);
+#
+#    std::pair<double,double> Res;
+#    Res.first=Loss;
+#    Res.second=two_std;


Property changes on: SwiftApps/sumloss.sh
___________________________________________________________________
Added: svn:executable
   + 

Added: SwiftApps/t1.py
===================================================================
--- SwiftApps/t1.py	                        (rev 0)
+++ SwiftApps/t1.py	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,22 @@
+#! /usr/bin/env python
+import os
+
+# FULL for target in range(58,1009,50):
+# FAST for target in range(58,209,50):
+for target in range(58,59,50):
+	s = ("%d" % target)
+	print s
+
+# FULL	for i in range(15):
+# FAST	for i in range(2):
+	for i in range(1):
+		args="./toptimizer 0 0 4 50 -1 "+s+" 40000 20 1000 2 1 2. 0.01 100 0.3 2.3 1 1 0 0 0 # > out.T"+str(target)+".i"+str(i)
+		args="OMP_NUM_THREADS=24 ./toptimizer 0 0 4 50 -1 "+s+" 40000 20 100 2 1 2. 0.01 2 0.3 2.3 1 1 0 0 0 # > out.T"+str(target)+".i"+str(i)
+		args="OMP_NUM_THREADS=24 ./toptimizer 0 0 4 50 -1 "+s+" 40000 20 100 2 1 2. 0.01 2 0.3 2.3 1 1 1 0 0 # > out.T"+str(target)+".i"+str(i)
+		print("\n\n **** CALLING APP: "+args+"\n\n\n")
+	        os.system(args);
+#		print("\n\n **** CALLING APP: ./optimizer 0 0 4 50 -1 "+s+" 40000 20 10 2 1 2. 0.01 100 0.3 2.3 1 1 0 0 0\n\n\n")
+#		FAST: os.system("./optimizer 0 0 4 50 -1 "+s+" 40000 20 10 2 1 2. 0.01 2 0.3 2.3 1 1 0 0 0")
+
+
+print "Done!"


Property changes on: SwiftApps/t1.py
___________________________________________________________________
Added: svn:executable
   + 

Added: SwiftApps/t2.cpp
===================================================================
--- SwiftApps/t2.cpp	                        (rev 0)
+++ SwiftApps/t2.cpp	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,44 @@
+#include <string>
+#include <sstream>
+#include <iostream>
+
+template <class T>
+bool from_string(T& t, 
+                 const std::string& s, 
+                 std::ios_base& (*f)(std::ios_base&))
+{
+  std::istringstream iss(s);
+  return !(iss >> f >> t).fail();
+}
+
+int main()
+{
+  int i;
+  float f;
+  
+  // the third parameter of from_string() should be 
+  // one of std::hex, std::dec or std::oct
+  if(from_string<int>(i, std::string("ff"), std::hex))
+    {
+      std::cout << i << std::endl;
+    }
+  else
+    {
+      std::cout << "from_string failed" << std::endl;
+    }
+  
+  if(from_string<float>(f, std::string("123.456"), std::dec))
+    {
+      std::cout << f << std::endl;
+    }
+  else
+    {
+      std::cout << "from_string failed" << std::endl;
+    }
+  return 0;
+} 
+
+/* output:
+255
+123.456
+*/

Added: SwiftApps/t2.py
===================================================================
--- SwiftApps/t2.py	                        (rev 0)
+++ SwiftApps/t2.py	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,23 @@
+#! /usr/bin/env python
+import os
+
+# FULL for target in range(58,1009,50):
+# FAST for target in range(58,209,50):
+for target in range(58,59,50):
+	s = ("%d" % target)
+	print s
+
+# FULL	for i in range(15):
+# FAST	for i in range(2):
+	for i in range(1):
+		args="./toptimizer 0 0 4 50 -1 "+s+" 40000 20 1000 2 1 2. 0.01 100 0.3 2.3 1 1 0 0 0 # > out.T"+str(target)+".i"+str(i)
+		args="OMP_NUM_THREADS=24 ./toptimizer 0 0 4 50 -1 "+s+" 40000 20 100 2 1 2. 0.01 2 0.3 2.3 1 1 0 0 0 # > out.T"+str(target)+".i"+str(i)
+		args="./toptimizer 0 0 4 50 -1 "+s+" 40000 20 75    2 1 2. 0.01 2 0.3 2.3 1 1 1 0 0 m # > out.T"+str(target)+".i"+str(i)
+		args="OMP_NUM_THREADS=24 ./toptimizer 0 0 4 50 -1 "+s+" 40000 20 96    2 1 2. 0.01 2 0.3 2.3 1 1 0 0 0 m 24 # > out.T"+str(target)+".i"+str(i)
+		print("\n\n **** CALLING APP: "+args+"\n\n\n")
+	        os.system(args);
+#		print("\n\n **** CALLING APP: ./optimizer 0 0 4 50 -1 "+s+" 40000 20 10 2 1 2. 0.01 100 0.3 2.3 1 1 0 0 0\n\n\n")
+#		FAST: os.system("./optimizer 0 0 4 50 -1 "+s+" 40000 20 10 2 1 2. 0.01 2 0.3 2.3 1 1 0 0 0")
+
+
+print "Done!"


Property changes on: SwiftApps/t2.py
___________________________________________________________________
Added: svn:executable
   + 

Added: SwiftApps/t3.py
===================================================================
--- SwiftApps/t3.py	                        (rev 0)
+++ SwiftApps/t3.py	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,22 @@
+#! /usr/bin/env python
+import os
+
+# FULL for target in range(58,1009,50):
+# FAST for target in range(58,209,50):
+for target in range(58,59,50):
+	s = ("%d" % target)
+	print s
+
+# FULL	for i in range(15):
+# FAST	for i in range(2):
+	for i in range(1):
+		args="./toptimizer 0 0 4 50 -1 "+s+" 40000 20 1000 2 1 2. 0.01 100 0.3 2.3 1 1 0 0 0 # > out.T"+str(target)+".i"+str(i)
+		args="OMP_NUM_THREADS=24 ./toptimizer 0 0 4 50 -1 "+s+" 40000 20 100 2 1 2. 0.01 2 0.3 2.3 1 1 0 0 0 # > out.T"+str(target)+".i"+str(i)
+		args="./toptimizer 0 0 4 50 -1 "+s+" 40000 20 75    2 1 2. 0.01 2 0.3 2.3 1 1 1 0 0 a # > out.T"+str(target)+".i"+str(i)
+		print("\n\n **** CALLING APP: "+args+"\n\n\n")
+	        os.system(args);
+#		print("\n\n **** CALLING APP: ./optimizer 0 0 4 50 -1 "+s+" 40000 20 10 2 1 2. 0.01 100 0.3 2.3 1 1 0 0 0\n\n\n")
+#		FAST: os.system("./optimizer 0 0 4 50 -1 "+s+" 40000 20 10 2 1 2. 0.01 2 0.3 2.3 1 1 0 0 0")
+
+
+print "Done!"


Property changes on: SwiftApps/t3.py
___________________________________________________________________
Added: svn:executable
   + 

Added: SwiftApps/tc
===================================================================
--- SwiftApps/tc	                        (rev 0)
+++ SwiftApps/tc	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,13 @@
+localhost sh /bin/sh null null null
+localhost cat /bin/cat null null null
+pbs cat /bin/cat null null null
+mcs cat /bin/cat null null null
+localhost catnap /home/wilde/swift/lab/catnap.sh null null GLOBUS::maxwalltime="00:01:00"
+beagle optimizer /home/wilde/AndreysOptimizer/Optimizer null null GLOBUS::maxwalltime="01:00:00"
+beagle optimizersh /home/wilde/AndreysOptimizer/optimizer.sh null null GLOBUS::maxwalltime="02:00:00"
+
+beagle evolve /home/wilde/AndreysOptimizer/evolve.sh null null GLOBUS::maxwalltime="02:00:00"
+localhost evolve /home/wilde/AndreysOptimizer/evolve.sh null null GLOBUS::maxwalltime="02:00:00"
+
+beagle sumloss /home/wilde/AndreysOptimizer/evolve.sh null null GLOBUS::maxwalltime="02:00:00"
+localhost sumloss /home/wilde/AndreysOptimizer/sumloss.sh null null GLOBUS::maxwalltime="02:00:00"

Added: SwiftApps/testO.sh
===================================================================
--- SwiftApps/testO.sh	                        (rev 0)
+++ SwiftApps/testO.sh	2012-01-20 21:26:24 UTC (rev 5502)
@@ -0,0 +1,21 @@
+#! /bin/bash
+
+python <<END
+import os
+
+target = 58
+s = ("%d" % target)
+
+args="0 0 4 50 -1 "+s+" 40000 20 1000 2 1 2. 0.01 100 0.3 2.3 1 1 0 0 0 "
+args="0 0 4 50 -1 "+s+" 40000 20  100 2 1 2. 0.01   2 0.3 2.3 1 1 0 0 0 "
+
+print("\n\n **** CALLING APP: "+args+"\n\n\n")
+os.system("./optimizer "+args+" >out.o.T"+s);
+os.system("./Optimizer "+args+" >out.O.T"+s);
+
+#		print("\n\n **** CALLING APP: ./optimizer 0 0 4 50 -1 "+s+" 40000 20 10 2 1 2. 0.01 100 0.3 2.3 1 1 0 0 0\n\n\n")
+#		FAST: os.system("./optimizer 0 0 4 50 -1 "+s+" 40000 20 10 2 1 2. 0.01 2 0.3 2.3 1 1 0 0 0")
+
+print "Done!"
+END
+


Property changes on: SwiftApps/testO.sh
___________________________________________________________________
Added: svn:executable
   + 




More information about the Swift-commit mailing list