[Swift-commit] r5997 - SwiftApps/SciColSim

wilde at ci.uchicago.edu wilde at ci.uchicago.edu
Sat Oct 27 09:31:40 CDT 2012


Author: wilde
Date: 2012-10-27 09:31:40 -0500 (Sat, 27 Oct 2012)
New Revision: 5997

Added:
   SwiftApps/SciColSim/gengraphs.sh
   SwiftApps/SciColSim/graphsim.cpp
   SwiftApps/SciColSim/testgraph.py
Modified:
   SwiftApps/SciColSim/Makefile
Log:
Changes to make test graphs and single-graph simulation.

Modified: SwiftApps/SciColSim/Makefile
===================================================================
--- SwiftApps/SciColSim/Makefile	2012-10-27 13:57:02 UTC (rev 5996)
+++ SwiftApps/SciColSim/Makefile	2012-10-27 14:31:40 UTC (rev 5997)
@@ -2,11 +2,14 @@
 
 ifeq ($(UNAME), Linux)
 
-all:	openmp-optimizer
+all:	openmp-optimizer graphsim
 
 openmp-optimizer: optimizer.cpp
 	g++ -DP_OPENMP -static -O -fopenmp -I boost_1_47_0 -o openmp-optimizer optimizer.cpp
 
+graphsim: graphsim.cpp
+	g++ -DP_OPENMP -static -O -fopenmp -I boost_1_47_0 -o graphsim graphsim.cpp
+
 clean:
 	@rm -rvf openmp-optimizer
 

Added: SwiftApps/SciColSim/gengraphs.sh
===================================================================
--- SwiftApps/SciColSim/gengraphs.sh	                        (rev 0)
+++ SwiftApps/SciColSim/gengraphs.sh	2012-10-27 14:31:40 UTC (rev 5997)
@@ -0,0 +1,9 @@
+#! /bin/sh
+
+mkdir -p graph
+
+for s in $(seq 1000 1000 10000); do
+  ./trimgraph.sh $s 1.0 > graph/$s.100
+  ./trimgraph.sh $s 0.1 > graph/$s.10
+  ./trimgraph.sh  $s 0.5 > graph/$s.50
+done
\ No newline at end of file


Property changes on: SwiftApps/SciColSim/gengraphs.sh
___________________________________________________________________
Added: svn:executable
   + *

Added: SwiftApps/SciColSim/graphsim.cpp
===================================================================
--- SwiftApps/SciColSim/graphsim.cpp	                        (rev 0)
+++ SwiftApps/SciColSim/graphsim.cpp	2012-10-27 14:31:40 UTC (rev 5997)
@@ -0,0 +1,1866 @@
+//
+//  main.cpp
+//  optimizer
+//
+//  Created by Andrey Rzhetsky on 4/11/11.
+//  Copyright 2011 University of Chicago. All rights reserved.
+//
+
+// Select OpenMP or Grand Central Dispatch for multithreading:
+
+#define MAXNworkers 1
+int Nworkers=MAXNworkers;
+
+#define NUniverse 1
+
+// Add operation code to enable existing code to be used at lower level from Swift scripts:
+
+char operation = 'n'; // n: normal; m: do one multi_loss (with n_reruns).
+                      // Not used: a: analyze and generate next annealing parameter set. g: tbd
+
+unsigned initSeed = 0;
+
+#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>
+
+#ifdef P_DISPATCH
+#include <dispatch/dispatch.h>
+#endif
+
+#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
+
+#define MAX_LOSS 160
+
+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 (int) round(r);  // FIXME: Andrey: please verify that round() is correct.
+
+	}
+
+    //=============================================
+    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 ( (int)round(Dist[*vi][j]>max_dist)) {
+						  // FIXME: Andrey: please verify that (int) cast is correct. Do we need to round()?
+						  // also, the indent on this iff statement was way off -
+						  // perhaps due to space v. tab?
+						  max_dist=(int)round(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.;
+        double ratio = 0.0;
+        int check = 0;
+
+        reset_world();
+
+        for (int k = istart; k < iend; k++){
+
+            for(int i=0; i< N_epochs &&  current_novelty < TargetNovelty; i++){
+
+                if(ratio < max_loss_value((int)TargetNovelty)){
+                    update_world();
+                }
+                else{
+                    check = 1;
+                    break;
+                }
+                ratio = current_loss/current_novelty;
+            }
+            if( check ){
+                storage[k]=ALOT;
+            }
+            else{
+                storage[k]=current_loss/current_novelty;
+            }
+            counters[k]=1;
+
+            reset_world();
+        }
+
+    }
+
+    //==============================================
+    double Omax_loss_value(int tn){
+
+        if( tn == 58 )           return 128.286;
+        else if( tn == 108 )     return 131.866;
+        else if( tn == 158 )     return 135.551;
+        else if( tn == 208 )     return 139.694;
+        else if( tn == 258 )     return 144.163;
+        else if( tn == 308 )     return 148.967;
+        else if( tn == 358 )     return 154.201;
+        else if( tn == 408 )     return 159.962;
+        else if( tn == 458 )     return 166.441;
+        else if( tn == 508 )     return 173.655;
+        else if( tn == 558 )     return 181.921;
+        else if( tn == 608 )     return 191.246;
+        else if( tn == 658 )     return 202.150;
+        else if( tn == 708 )     return 215.197;
+        else if( tn == 758 )     return 202.150;  // Verify with Andrey. Same as TargetNovelty of 658
+        else if( tn == 808 )     return 251.698;
+        else if( tn == 858 )     return 279.201;
+        else if( tn == 908 )     return 320.112;
+        else if( tn == 958 )     return 394.774;
+        else if( tn == 1008 )    return 1052.38;  // Huge number here. Why?
+        else                      return 10000.;   /* Pray that this does not occur. Ask Andrey what to do */
+    }
+    double max_loss_value(int tn){
+
+      switch(tn) {
+        case 58:      return 128.286;
+        case 108:     return 131.866;
+        case 158:     return 135.551;
+        case 208:     return 139.694;
+        case 258:     return 144.163;
+        case 308:     return 148.967;
+        case 358:     return 154.201;
+        case 408:     return 159.962;
+        case 458:     return 166.441;
+        case 508:     return 173.655;
+        case 558:     return 181.921;
+        case 608:     return 191.246;
+        case 658:     return 202.150;
+        case 708:     return 215.197;
+        case 758:     return 202.150;  // Verify with Andrey. Same as TargetNovelty of 658
+        case 808:     return 251.698;
+        case 858:     return 279.201;
+        case 908:     return 320.112;
+        case 958:     return 394.774;
+        case 1008:    return 1052.38;  // Huge number here. Why?
+        default:      return 10000.;   /* Pray that this does not occur. Ask Andrey what to do */
+      }
+    }
+
+
+    //==============================================
+    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
+
+};
+
+//============================================================
+
+#ifdef P_DISPATCH
+std::pair<double,double> multi_loss(dispatch_group_t group,
+                                    Universe* un[],
+                                    dispatch_queue_t* CustomQueues,
+                                    double* Results,
+                                    int*    Counters,
+                                    double* params){
+#else
+std::pair<double,double> multi_loss(Universe* un[],
+                                    double* Results,
+                                    int*    Counters,
+                                    double* params){
+
+#endif
+
+    int N = un[0]->get_reruns();
+    int step = (int)floor((double)N/(double)(Nworkers)); // FIXME: Andrey: please check change in cast grouping and use of floor
+    int istart=0;
+    int iend = istart+step;
+
+    double Loss=0., LossSquare=0.;
+
+    timeval startTime, endTime;
+    double elapsedTime;
+    gettimeofday(&startTime, NULL);
+
+    if ( N > NUniverse ) {
+      std::cerr << "Error: Number of reruns=" << N << " is greater than max NUniverse=" << NUniverse << "\n";
+      exit(1);
+    }
+
+    for(int i=0; i<N; i++){
+        for(int j=0; j<5; j++){
+            un[i]->set_parameter(params[j],j);
+        }
+    }
+
+#ifdef P_DISPATCH
+    for(int i=0; i<Nworkers; i++){
+
+        dispatch_group_async(group, CustomQueues[i], ^{
+
+            // 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 << "queued: i=" << i << " N=" << N << " istart=" << istart << " iend=" << iend << "\n";
+        // istart += step;
+        // iend = min(istart+step,N);
+
+    }
+    dispatch_group_wait(group, DISPATCH_TIME_FOREVER);
+    //dispatch_release(group);
+#else
+    int i;
+
+    // Execute rerun loop in parallel
+
+    #pragma omp parallel for private (i)
+    for(i=0; i<N; i++){
+        un[i]->evolve_to_target_and_save(i, i+1, Results, Counters);
+     }
+#endif
+
+    for (int i=0; i<N; i++){
+        Loss+=Results[i]/(double)N;
+        LossSquare+=Results[i]*Results[i]/(double)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 << ", target=" << un[0]->get_target() << ") elapsed time: " << elapsedTime << " seconds " << elapsedTime/60. << " minutes\n\n";
+
+    return Res;
+}
+//============================================================
+
+
+//============================================================
+#ifdef P_DISPATCH
+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){
+#else
+void multi_annealing(Universe* un[],
+                     double T_start, double T_end,
+                     double Target_rejection,
+                     int Annealing_repeats,
+                     double starting_jump,
+                     double* Results,
+                     int*    Counters,
+                     double* params0,
+                     double annealing_cycles){
+#endif
+    //.................................
+    // 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;
+
+#ifdef P_DISPATCH
+    Res = multi_loss(group, un, CustomQueues, Results, Counters, x);
+#else
+    Res = multi_loss(       un,               Results, Counters, x);
+#endif
+    std::cout << "Returned from initial multi_loss:" << std::endl;
+    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<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);
+                }
+
+#ifdef P_DISPATCH
+                Res = multi_loss(group, un, CustomQueues, Results, Counters, x);
+#else
+                Res = multi_loss(       un,               Results, Counters, x);
+#endif
+                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
+
+		    << "N, " << i << ", " << j << ", " << dx[j] << ", " << rejection[j] << ", |, " // FIXME: MW-DEBUGGING
+
+                    << 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();
+
+                    FILE *bf;
+                    bf = fopen("bestdb.txt","a");
+                    fprintf(bf, "N %2d %2d %10.5f %5.2f | %5.2f %10.5f [ %5.2f %5.2f %10.5f %10.5f %10.5f ] %10.5f\n", i, j, dx[j], rejection[j],
+			    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);
+                    fclose(bf);
+                }
+
+
+                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[3] = {"Operation", "Nworkers", "initSeed"};
+
+    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[NUniverse];
+#ifdef P_DISPATCH
+    static dispatch_queue_t CustomQueues[MAXNworkers];
+#endif
+    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]= (int)round(atof(argv[nArg]));   // FIXME: Andrey: please verify that round() is correct.
+	        var_fixed[nArg-17]= (int)round(atof(argv[nArg])); // FIXME: ditto
+                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;
+            }
+            if (nArg == 24 ){
+                initSeed = atoi(argv[nArg]);
+                std::cout << par_names4[2] << ": " << initSeed <<  std::endl;
+            }
+
+
+        }
+
+    }
+
+    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<NUniverse; i++){
+        un[i] = new Universe((char *)filename,n_ep,n_st,
+                             (int)n_rep,
+                             identify_failed, target, i2string(i));
+#ifdef P_DISPATCH
+        CustomQueues[i] = dispatch_queue_create(i2char(i), NULL);
+#endif
+    }
+
+    //...............................
+    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();
+
+    if ( initSeed != 0.0 ) {
+      srand(initSeed);
+    }
+    else {
+        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];
+
+#ifdef P_DISPATCH
+    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);
+    //.............................
+#else
+    multi_annealing(       un,                T_start, T_end, Target_rejection, Annealing_repeats,
+                    starting_jump, Results, Counters, params0, Annealing_repeats);
+#endif
+
+
+    // 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 << "\n*** optimizer completed, elapsed time=" << 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;
+
+
+
+}
+

Added: SwiftApps/SciColSim/testgraph.py
===================================================================
--- SwiftApps/SciColSim/testgraph.py	                        (rev 0)
+++ SwiftApps/SciColSim/testgraph.py	2012-10-27 14:31:40 UTC (rev 5997)
@@ -0,0 +1,223 @@
+#! /usr/bin/env python 
+
+import os, sys
+
+unbuffered = os.fdopen(sys.stdout.fileno(), 'w', 0)
+sys.stdout = unbuffered
+
+app              = "./orig-optimizer";     # For Mac only: original code (+1 loop fix) using Grnd Central Dispatch
+app              = "./dispatch-optimizer"; # For Mac only: sing Grand Central Dispatch
+app              = "./openmp-optimizer";   # For Mac or Linux: Using OpenMP (Default)
+
+#paramset="default"
+paramset="mw3"
+#paramset="VERYFASTTEST"
+#paramset="VERYFASTTEST"
+
+startTarget = int(sys.argv[1])
+print("Running with startTarget " + str(startTarget) + "\n\n");
+
+print("Running with parameter set " + paramset + "\n\n");
+
+if paramset == "beta":
+  startTarget = 58
+  endTarget = 59
+  incrTarget = 50
+  optimizerRepeats = 100
+  evolveReruns = 32
+  annealingSteps = 100
+  NWorkers = "8"
+  seed = ""
+  openmp = "OMP_NUM_THREADS=" + NWorkers
+  operation = "n"
+  
+
+elif paramset == "default":
+  
+  # Parameter names and optimizer command line positions
+  
+  startTarget      = 58    # target (target innovation)
+  endTarget        = 1009  # "
+  incrTarget       = 50    # "
+  optimizerRepeats = 15    # wrapper control var (# of repeats of optimizer to do for each target loss)
+  evolveReruns     = 1000  # n_reruns = 
+  annealingSteps   = 100
+
+# Following params are for modified optimizer openmp-optimizer
+
+# NWorkers         = "2"
+# operation        = "n"  # n=normal, m=manual (runs 1 multi_loss call)
+# seed             = "1234567"
+# openmp           = "OMP_NUM_THREADS=24"
+
+  NWorkers         = ""
+  operation        = ""
+  seed             = ""
+  openmp           = ""
+
+elif paramset == "FULLSCALE":  # FULL INITIAL TEST PARAMETERS from Andrey
+  
+  startTarget      = 58
+  endTarget        = 1009
+  incrTarget       = 50
+  optimizerRepeats = 15
+  evolveReruns     = 1000
+  annealingSteps   = 100
+
+elif paramset == "FASTTEST":  # FAST TEST PARAMETERS by Mike
+  
+  startTarget      = 58
+  endTarget        = 209
+  incrTarget       = 50
+  optimizerRepeats = 2
+  evolveReruns     = 10
+  annealingSteps   = 2
+  NWorkers         = "2"
+  openmp           = "OMP_NUM_THREADS=" + NWorkers
+  operation        = "n"  # n=normal, m=manual (runs 1 multi_loss call)
+  seed             = "1234567"
+
+elif paramset == "VERYFASTTEST":  #VERY FAST TEST PARAMETERS by Ketan
+  
+  startTarget      = 58
+  endTarget        = 209
+  incrTarget       = 50
+  optimizerRepeats = 2
+  evolveReruns     = 1
+  annealingSteps   = 2
+  NWorkers         = "1"
+  openmp           = "OMP_NUM_THREADS=" + NWorkers
+  operation        = "n"  # n=normal, m=manual (runs 1 multi_loss call)
+  seed             = "1234567"
+
+elif paramset == "MANUALTEST":  # for manual (one multi_loss call) mode
+  
+  startTarget      = 258  # runtime increases with greater target
+  endTarget        = 259
+  incrTarget       = 50
+  optimizerRepeats = 1
+  evolveReruns     = 60
+  annealingSteps   = 5
+  NWorkers         = "4"
+  openmp           = "OMP_NUM_THREADS=" + NWorkers
+  operation        = "m"  # n=normal, m=manual (runs 1 multi_loss call)
+  seed             = "1234567"
+
+elif paramset == "mw":  # FAST TEST PARAMETERS by Mike
+  
+  startTarget      = 58
+  endTarget        = 59
+  incrTarget       = 50
+  optimizerRepeats = 1
+  evolveReruns     = 100
+  annealingSteps   = 5
+  NWorkers         = "4"
+  openmp           = "OMP_NUM_THREADS=" + NWorkers
+  operation        = "n"  # n=normal, m=manual (runs 1 multi_loss call)
+  seed             = "1234567"
+  app              = "./openmp-optimizer";
+
+elif paramset == "mw2":  # Timing TEST PARAMETERS by Mike
+  
+  startTarget      = 58
+  endTarget        = 209
+  incrTarget       = 50
+  optimizerRepeats = 1
+  evolveReruns     = 240
+  annealingSteps   = 1
+  NWorkers         = "8"
+  openmp           = "OMP_NUM_THREADS=" + NWorkers
+  operation        = "m"  # n=normal, m=manual (runs 1 multi_loss call)
+  seed             = "123456"
+  app              = "./openmp-optimizer";
+
+elif paramset == "mw3":  # Timing TEST PARAMETERS by Mike
+  
+  #  startTarget      = 9800
+  endTarget        = startTarget+1
+  incrTarget       = 50
+  optimizerRepeats = 1
+  evolveReruns     = 1
+  annealingSteps   = 1
+  NWorkers         = "8"
+  openmp           = "OMP_NUM_THREADS=" + NWorkers
+  operation        = "m"  # n=normal, m=manual (runs 1 multi_loss call)
+  seed             = "123456"
+  app              = "./openmp-optimizer";
+
+# Ensure we dont pass new parameters to original optimizer versions
+# (but they would be ignored)
+
+if app == "./orig-optimizer":
+  NWorkers = ""
+  operation = ""
+  seed = ""
+
+# Run the test
+
+for target in range(startTarget,endTarget,incrTarget):
+  for i in range(optimizerRepeats):
+    args = "rm -f bestdb.txt; " + \
+           openmp + " " + app + " 0 0 0 0 0 " + str(target) + " 40000 20 " + str(evolveReruns) + \
+           " 2 1 2. 0.01 " + str(annealingSteps) + " 0.3 2.3 0 0 0 0 0 " + operation + " " + NWorkers + " " + seed + \
+           " # >& out.T"+str(target)+".i"+str(i) + "; #mv bestdb.txt best.T" + str(target) + ".i" + str(i) 
+    print("\n**** Calling optimizer: "+args+"\n")
+    os.system(args);
+
+print sys.argv[0] + " Done!"
+
+#            openmp + " " + app + " 0 0 4 50 -1 " + str(target) + " 40000 20 " + str(evolveReruns) + \
+
+###   Argument worksheet:
+
+#       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"];
+# 
+#
+#   long running args = "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 24? 12345
+# 
+#   [ alpha_i alpha_m beta gamma delta target_innov ]
+#     0       0       4    50    -1    target
+#         
+#   [ n_epochs n_steps n_reruns ]  [ range ] [ verbose_level]
+#     40000    20      1000          2         1
+#         
+#   [ T_start T_end Annealing_steps Target_rejection Starting_jump ]
+#     2.      0.01  100             0.3              2.3 
+# 
+#   [ FREEZE_alpha_i FREEZE_alpha_m FREEZE_beta FREEZE_gamma FREEZE_delta ]
+#     1              1              0           0            0
+# 
+#   [ operation-code:(n,m) Nworkers seed ]
+#     n                    24       12345 
+# 
+
+# From older test, using separate file per run:
+#    args = app + " 0 0 4 50 -1 " + str(target) + " 40000 20 " + str(evolveReruns) + \
+#      " 2 1 2. 0.01 " + str(annealingSteps) + " 0.3 2.3 1 1 0 0 0 > out.T"+str(target)+".i"+str(i)
+
+### Optimization loop logic:
+
+# a) 20 targets (parallel)
+# b)   15 repeats (parallel) of optimizer (optimizer == multi_annealing):
+#        1 multi_loss to initialize
+# c)     100 Annealing_cycles (serial)
+# d)        3 repeats (1 per non-fixed param, serial) of multi_loss:
+# e)          1000 to 10000 evolve_reruns (parallel) == multi_loss
+# f)            evolve()


Property changes on: SwiftApps/SciColSim/testgraph.py
___________________________________________________________________
Added: svn:executable
   + *




More information about the Swift-commit mailing list