[Swift-commit] r7764 - SwiftApps/SimulatedAnnealing
wilde at ci.uchicago.edu
wilde at ci.uchicago.edu
Fri Apr 11 16:43:06 CDT 2014
Author: wilde
Date: 2014-04-11 16:43:06 -0500 (Fri, 11 Apr 2014)
New Revision: 7764
Added:
SwiftApps/SimulatedAnnealing/sa.swift
Log:
New simpler kinder gentler annealer.
Added: SwiftApps/SimulatedAnnealing/sa.swift
===================================================================
--- SwiftApps/SimulatedAnnealing/sa.swift (rev 0)
+++ SwiftApps/SimulatedAnnealing/sa.swift 2014-04-11 21:43:06 UTC (rev 7764)
@@ -0,0 +1,388 @@
+import "math";
+
+type file;
+
+type Res
+{
+ float result;
+/*
+ float sdev;
+ float tavg;
+ float tsdev;
+*/
+}
+
+type Checkpoint
+{
+ int i;
+ int j;
+ float dx0, dx1, dx2, dx3, dx4;
+ float rej0, rej1, rej2, rej3, rej4;
+ float result;
+ float alpha_i;
+ float alpha_m;
+ float beta;
+ float gamma;
+ float delta;
+}
+
+global boolean restart = (@arg("restart","false") == "true");
+tracef("Restart flag is %b\n", restart);
+
+global int iterate_adjust = @toint(@arg("iterate_adjust","1")); // FIXME!!! : 1 for 0.93, 0 for 0.94
+
+global boolean FIX_VARIABLES = true;
+global int var_fixed[] = [0,0,0,0,0];
+global int Nworkers = @toint(@arg("nworkers","4"));
+global int rerunsPerApp= @toint(@arg("rerunsperapp", "100"));
+
+( float nx ) newx( float x, float dx )
+{
+ float r = (random());
+
+ if (r > 0.5) {
+ nx = x + (random())*dx; //Java already returns a float between [0-1]
+ }
+ else {
+ nx = x - (random())*dx;
+ }
+ // tracef("newx(%f,%f)=%f\n",x,dx,nx);
+}
+
+/*
+function y = simple_objective(x)
+ y = (4 - 2.1*x(1)^2 + x(1)^4/3)*x(1)^2 + x(1)*x(2) + (-4 + 4*x(2)^2)*x(2)^2;
+*/
+
+(float y) objective ( float x[] )
+{
+ float x0 = x[0];
+ float x1 = x[1];
+ float x0s = x0 * x0;
+ float x1s = x1 * x1;
+ float r = random() * .1;
+ tracef("r=%f\n",r);
+ y = ( (4.0 - 2.1*x0s + pow(x[0],4.0/3.0))*x0s + x0*x1 + (-4.0 + 4.0*x1s)*x1s ) + r;
+}
+
+(float s) sum (float a[])
+{
+ string t[] = system("echo " + strjoin(a,"+") + "|bc");
+ s = toFloat(t[0]);
+}
+
+(float a) avg (float f[])
+{
+ float n = toFloat( length(f) );
+ float t1 = sum(f);
+ float t2 = t1 / n;
+ a = t2;
+}
+
+/* Program structure:
+
+ main
+ optimizer_sweep() - do N complete annealing optimizations
+ multi_annealing()
+ N_objective()
+ objective()
+ sumobjective()
+*/
+
+(file bestfile, file maxfile) multi_annealing (string best_filename, # app
+ float T_start, # ann
+ float T_end, # ann
+ float Target_rejection,# ann
+ int obj_runs, # app
+ float starting_jump, # ann
+ float params0[], # app f(x)
+ int annealing_cycles) # ann
+{
+ int cycle = 10; // const
+ int NEVOPARAMS = 5; // const - 5 params: alpha_i, alpha_m, beta, gamma, delta
+
+ float rejection[][]; // [i][j] where i is cycle and j is evolve-parameter (alpha_i, alpha_m, beta, gamma, delta)
+ float trejection[][];
+
+ float x[][], dx[][], curr_result[];
+
+ Res objres[][];
+
+ file ckptFile <single_file_mapper; file=@strcat("ckpt")>;
+ Checkpoint ckpt;
+ int restartIndex;
+
+ if (restart) {
+ ckpt = readData(ckptFile);
+ restartIndex = ckpt.i;
+ }
+ else {
+ restartIndex = -1;
+/*
+ mlres[0][0] = multi_loss(T_start, T_end, annealing_cycles, Target_rejection,
+ starting_jump, 0, 0, params0, target_innov, evolve_reruns ); // FIXME: serves for all evolve-params ???
+*/
+ objres[0][0] = N_objective(params0, obj_runs ); // FIXME: serves for all evolve-params ???
+
+ tracef( "multi_annealing: initial: %f\n", objres[0][0].result );
+
+ foreach j in [0:NEVOPARAMS-1] {
+ x[0][j] = params0[j];
+ dx[0][j] = starting_jump;
+ rejection[0][j] = 0.0;
+ curr_result[j] = objres[0][0].result;
+ }
+ }
+
+ iterate iter_i { // number of annealing cycles
+ int i = iter_i + 1; // i ranges [1:n] in the swift script so that [0] can be the initial condition
+
+ // set new temperature, rejection threshold, and dx values for this cycle
+ float temperature = T_start*exp( @tofloat( i-1 ) * ( jlog( T_end ) - jlog( T_start ) ) / @tofloat( annealing_cycles ) );
+
+ tracef("multi_annealing: i=%i ....T =%f\n", i, temperature );
+
+ if ( i < restartIndex ) {
+ tracef( "skipping index %i - less than restart index %i\n", i, restartIndex );
+ }
+ else { if ( i == restartIndex ) { // reset variables from the restart file
+ dx[i] = [ckpt.dx0, ckpt.dx1, ckpt.dx2, ckpt.dx3, ckpt.dx4];
+ rejection[i] = [ckpt.rej0, ckpt.rej1, ckpt.rej2, ckpt.rej3, ckpt.rej4];
+ x[i][0] = ckpt.alpha_i;
+ x[i][1] = ckpt.alpha_m;
+ x[i][2] = ckpt.beta;
+ x[i][3] = ckpt.gamma;
+ x[i][4] = ckpt.delta;
+ curr_result[((i+1)*NEVOPARAMS)-1] = ckpt.result;
+ }
+ else { // i > restartIndex: proceed as normal whether restarting or not
+
+ // 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) {
+ tracef( "multi_annealing: new cycle at i=%i\n", i );
+ tracef( "multi_annealing: New cycle at %i: prev dx[0-4]=[%f %f %f %f %f]\n",
+ i, dx[i-1][0], dx[i-1][1], dx[i-1][2], dx[i-1][3], dx[i-1][4] );
+ 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 );
+ // FIXME: re-enable: rejection[i][k]=0.0;
+ trejection[i][k]=0.0;
+ }
+ else {
+ dx[i][k] = dx[i-1][k] * 2.0;
+ // FIXME: re-enable: rejection[i][k]=rejection[i-1][k];
+ trejection[i][k]=newrejection;
+ }
+ // FIXME: HANGS? : tracef("Recomputed rejection: i=%d k=%d dx[i][k]=%f\n", i, k, dx[i][k]);
+ }
+ tracef( "multi_annealing: New cycle at %i: dx[0-4]=[%f %f %f %f %f]\n",
+ i, dx[i][0], dx[i][1], dx[i][2], dx[i][3], dx[i][4] );
+ }
+ else { // If not new cycle, set dx[i][*] from previous dx ([i-1]). rejection[i][j] is set later.
+ foreach k in [0:NEVOPARAMS-1] {
+ dx[i][k] = dx[i-1][k];
+ }
+ }
+ iterate j { // Try a new value for each non-fixed param; then write results and accept or reject
+ int curr = ( i * NEVOPARAMS ) + j;
+ int prev = curr-1;
+
+ if ( /* (!FIX_VARIABLES) || */ ( var_fixed[j] == 0 ) ) { // Adjustable vars
+ // fixed=1,1,0,0,0: FIXME: FIX_VARIABLES flag has faulty logic but OK when TRUE
+ float try_x[];
+ foreach k in [0:NEVOPARAMS-1] { // 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-1][j], dx[i-1][j]); // permute x[i-1][j]
+ }
+ else { // k > j
+ try_x[k] = x[i-1][k]; // use x[i-1][k] (from prior cycle)
+ }
+ }
+ }
+ tracef( "multi_annealing: %f %i\n", try_x[j], j );
+
+ # Call objective function here:
+
+# mlres[i][j] = multi_loss( T_start, T_end, annealing_cycles, Target_rejection, starting_jump,
+# i, j, try_x, target_innov, evolve_reruns ); // do the N evolve()'s, N=evolve_rerusn
+ objres[i][j] = N_objective( try_x, obj_runs );
+
+ tracef( "multi_annealing: %f\n", objres[i][j].result );
+ // Beyond this point, x[] and dx[] are being set for this i,j
+ fprintf( best_filename, "%i, %i, %f, %f, | %f [ %f, %f, %f, %f, %f ]\n",
+ i-1, j, dx[i][j], rejection[i][j], objres[i][j].result,
+ try_x[0], try_x[1], try_x[2], try_x[3], try_x[4] );
+ // Note i-1 field: print that way to match with C++ output
+
+ // fprintf( "max_dist_swift.txt", color( Red,"multi_annealing: AF: max_dist.txt - tbd\n" ) );
+ // FIXME: max_dist is global set in evolve()
+ float ratio = min( 1.0, exp( -( objres[i][j].result - curr_result[prev] ) / temperature ) );
+ float r = (random()); //Java already returns a random float between [0.0-1.0]
+ tracef("multi_annealing: AR: %f vs %f\n", r, ratio);
+ if (r > ratio) { // Reject new parameter
+ x[i][j] = x[i-1][j];
+ if (i %% cycle == 1 && i > 1) {
+ rejection[i][j] = trejection[i][j] + 1.0; // FIXME: triple-check this!
+ }
+ else {
+ rejection[i][j] = rejection[i-1][j] + 1.0; // FIXME: AR: Is this correct? incr rejection?
+ }
+ curr_result[curr] = curr_result[prev];
+ // FIXME: AR: the following prints seem to replicate values in the .cpp version - please clarify.
+ tracef( "multi_annealing: %i,%i %i Did not accept: %f (%i)\n", i, j, i, try_x[j], j );
+ tracef( "multi_annealing: %f %f %f %f %f\n", try_x[0],try_x[1],try_x[2],try_x[3],try_x[4] );
+ }
+ else { // Accept new parameter
+ tracef( "multi_annealing: Accepting try_x[j], i=%i j=%i\n",i,j );
+ x[i][j] = try_x[j];
+ if (i %% cycle == 1 && i > 1) {
+ rejection[i][j] = trejection[i][j]; // FIXME: triple-check this!
+ }
+ else {
+ rejection[i][j] = rejection[i-1][j]; // FIXME: AR: Is this correct? no incr of rejection?
+ }
+ curr_result[curr] = objres[i][j].result;
+ tracef( "multi_annealing: Accepting try_x[j], i=%i j=%i try_x[j]=%f\n", i, j, try_x[j] );
+ float rj[];
+ foreach k in [0:NEVOPARAMS-1] {
+ if (k <= j) {
+ rj[k] = rejection[i][k]; // Was either set from previous j or just set for this j
+ }
+ else {
+ rj[k] = rejection[i-1][k]; // Not yet set, use previous
+ }
+ }
+ tracef("multi_annealing: [%i][%i] Rejection counts: %f %f %f %f %f\n\n",
+ i, j, rj[0], rj[1], rj[2], rj[3], rj[4]);
+ tracef("multi_annealing: %i ***** Did accept! %f %f %f %f %f\n\n",
+ i, try_x[0], try_x[1], try_x[2], try_x[3], try_x[4]);
+ }
+ }
+ else {// Fixed Vars
+ x[i][j] = x[i-1][j];
+ rejection[i][j] = rejection[i-1][j];
+ curr_result[curr] = curr_result[prev];
+ // dx[i][j] not set for fixed vars
+ }
+ } until( j == (NEVOPARAMS-iterate_adjust) );
+ }} // of if/else
+ } until( iter_i == (annealing_cycles-iterate_adjust) );
+}
+
+(Res r) N_objective(float x[], int obj_runs)
+{
+ tracef("%q\n", x);
+ float result[];
+
+ tracef( "N_objective: entered: obj_runs=%i x=%q\n", obj_runs, x );
+
+ foreach i in [1:obj_runs] { // repeats of the objective() - same as n_reruns
+ result[i] = objective(x);
+ }
+
+ r.result = avg(result);
+ tracef("multi_loss: returning: r.result=%f\n", r.result);
+}
+
+optimizer_serial_sweep() // Implements logic of python driver script
+{
+
+ int minrange = @toint(@arg("minrange", "58"));
+ int maxrange = @toint(@arg("maxrange", "59"));
+ int rangeinc = @toint(@arg("rangeinc", "50"));
+
+ // FIXME: add provision for random priming and random param values when x[i] == -100 (see optimizer.cpp main())
+
+ int nreps = @toint(@arg("nreps", "1"));
+
+//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] {
+# iterate i {
+# int target_innov = minrange + (i*rangeinc);
+ foreach rep in [1:nreps] {
+
+ string best_filename = @strcat( "best.R", @strcat(rep), ".txt" );
+ file outfile; // <single_file_mapper; file=@strcat("output/R",rep,".out")>;
+ file lossfile; // <single_file_mapper; file=@strcat("output/R",rep,".loss_data")>;
+
+ (outfile,lossfile) = multi_annealing(best_filename,
+ @tofloat(@arg("tstart", "2.0")),
+ @tofloat(@arg("tend", "0.01")),
+ @tofloat(@arg("trejection", "0.3")),
+ @toint(@arg("evoreruns", "100")),
+ @tofloat(@arg("startingjump", "2.3")),
+ [ @tofloat(@arg("alphai", "0.0")),
+ @tofloat(@arg("alpham", "0.0")),
+ @tofloat(@arg("beta", "4.0")),
+ @tofloat(@arg("gamma", "50.0")),
+ @tofloat(@arg("delta", "-1.0"))
+ ],
+ @toint(@arg("annealingcycles", "50")) );
+ }
+ # } until(target_innov >= (maxrange-rangeinc));
+}
+
+optimizer_sweep() // Implements logic of python driver script
+{
+
+ int minrange = @toint(@arg("minrange", "58"));
+ int maxrange = @toint(@arg("maxrange", "59"));
+ int rangeinc = @toint(@arg("rangeinc", "50"));
+
+ // FIXME: add provision for random priming and random param values when x[i] == -100 (see optimizer.cpp main())
+
+ int nreps = @toint(@arg("nreps", "1"));
+
+//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] {
+
+ string best_filename = @strcat( "best.R", @strcat(rep), ".txt" );
+ file outfile; // <single_file_mapper; file=@strcat("output/R",rep,".out")>;
+ file lossfile; // <single_file_mapper; file=@strcat("output/R",rep,".loss_data")>;
+
+ (outfile,lossfile) = multi_annealing(best_filename,
+ @tofloat(@arg("tstart", "2.0")),
+ @tofloat(@arg("tend", "0.01")),
+ @tofloat(@arg("trejection", "0.3")),
+ @toint(@arg("evoreruns", "100")),
+ @tofloat(@arg("startingjump", "2.3")),
+ [ @tofloat(@arg("alphai", "0.0")),
+ @tofloat(@arg("alpham", "0.0")),
+ @tofloat(@arg("beta", "4.0")),
+ @tofloat(@arg("gamma", "50.0")),
+ @tofloat(@arg("delta", "-1.0"))
+ ],
+ @toint(@arg("annealingcycles", "50")) );
+ }
+# }
+}
+
+main()
+{
+ # optimizer_sweep();
+
+ float a[][];
+ a[0] = [1.1,1.2,1.3,1.4];
+ float y[];
+ y[0] = objective(a[0]);
+
+ iterate c {
+ float na[] = [ newx(a[0],0.5), newx(a[1],0.5), newx(a[2],0.5), newx(a[3],0.5) ];
+ float ny = objective(na);
+ trace(c,ny);
+ } until (c==10);
+
+}
+
+main();
More information about the Swift-commit
mailing list