[Swift-commit] r6110 - SwiftApps/SciColSim

wilde at ci.uchicago.edu wilde at ci.uchicago.edu
Sat Dec 22 22:33:23 CST 2012


Author: wilde
Date: 2012-12-22 22:33:21 -0600 (Sat, 22 Dec 2012)
New Revision: 6110

Added:
   SwiftApps/SciColSim/maxloss.cpp
   SwiftApps/SciColSim/maxloss.py
   SwiftApps/SciColSim/maxloss_ar.py
Modified:
   SwiftApps/SciColSim/Makefile
   SwiftApps/SciColSim/TODO
   SwiftApps/SciColSim/graphsim.cpp
   SwiftApps/SciColSim/samplegraph.sh
   SwiftApps/SciColSim/testgraph.py
Log:
Development snapshot. Test programs for dynamic max-loss; fine-grained OpenMP of graphsim; revised graph sampler.

Modified: SwiftApps/SciColSim/Makefile
===================================================================
--- SwiftApps/SciColSim/Makefile	2012-12-19 14:06:35 UTC (rev 6109)
+++ SwiftApps/SciColSim/Makefile	2012-12-23 04:33:21 UTC (rev 6110)
@@ -7,8 +7,11 @@
 openmp-optimizer: optimizer.cpp
 	g++ -DP_OPENMP -static -O -fopenmp -I boost_1_47_0 -o openmp-optimizer optimizer.cpp
 
+#t: t.cpp
+#	g++ -DP_OPENMP -static -O -fopenmp -I boost_1_47_0 -o t t.cpp
+
 graphsim: graphsim.cpp
-	g++ -DP_OPENMP -static -pg -O -fopenmp -I boost_1_47_0 -o graphsim graphsim.cpp
+	g++ -DP_OPENMP -static -pg -O2 -fopenmp -I boost_1_47_0 -o graphsim graphsim.cpp
 
 openmptest: openmptest.cpp
 	g++ -DP_OPENMP -static -O -fopenmp -o openmptest openmptest.cpp

Modified: SwiftApps/SciColSim/TODO
===================================================================
--- SwiftApps/SciColSim/TODO	2012-12-19 14:06:35 UTC (rev 6109)
+++ SwiftApps/SciColSim/TODO	2012-12-23 04:33:21 UTC (rev 6110)
@@ -1,6 +1,28 @@
 
+--- Questions and Open Issues ---
 
+When sampling the graph, we compress the node ids. Is this OK?
 
+Same 40K epoch count?  Is this really feasible?
+
+What values for otehr loop parameters?
+
+What sampling?
+
+Limits?
+-- placed correctly?
+-- size
+-- anomalous value for tn 758
+-- adjust all for new graph?
+
+What is range of TI to sweep over?
+
+Algorithmic improvements?
+
+Graph partitioning? => for greater total parallelism?
+
+
+
 x = done, - = in progress
 
 

Modified: SwiftApps/SciColSim/graphsim.cpp
===================================================================
--- SwiftApps/SciColSim/graphsim.cpp	2012-12-19 14:06:35 UTC (rev 6109)
+++ SwiftApps/SciColSim/graphsim.cpp	2012-12-23 04:33:21 UTC (rev 6110)
@@ -6,8 +6,6 @@
 //  Copyright 2011 University of Chicago. All rights reserved.
 //
 
-// Select OpenMP or Grand Central Dispatch for multithreading:
-
 #define MAXNworkers 1
 int Nworkers=MAXNworkers;
 
@@ -34,13 +32,8 @@
 #include <sys/types.h>
 #include "unistd.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>
@@ -88,157 +81,122 @@
 typedef boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 
 namespace std {
-	using ::time;
+  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;
+			 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
+  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 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 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.);
 
-    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;
-
+  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();
+  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();
+  std::ostringstream s;
+  s << string_wrap(strDouble(val),mode);
+  return s.str();
 }
 
-
-
-//===============================================
 const
 string i2string(int i){
-
-    std::ostringstream s;
-    s << "worker"
+  std::ostringstream s;
+  s << "worker"
     << lexical_cast<std::string>(i);
-
-    return s.str();
-
+  return s.str();
 }
 
-//===============================================
 char* i2char(int i){
-
-    std::ostringstream s;
-    s << "worker"
+  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;
+  char* a=new char[s.str().size()+1];
+  memcpy(a,s.str().c_str(), s.str().size());
+  return a;
 }
 
-//================================================
-
 typedef struct {
   unsigned long size,resident,share,text,lib,data,dt;
 } statm_t;
@@ -253,13 +211,12 @@
     perror(statm_path);
     exit(99);
   }
-
   if(7 != fscanf(f,"%ld %ld %ld %ld %ld %ld %ld",
 		 &sp->size,&sp->resident,&sp->share,&sp->text,&sp->lib,&sp->data,&sp->dt))
-    {
-      perror(statm_path);
-      exit(99);
-    }
+  {
+    perror(statm_path);
+    exit(99);
+  }
   fclose(f);
 }
 
@@ -277,529 +234,522 @@
 	 MB(size), MB(resident), MB(share), MB(text), MB(lib), MB(data), MB(dt) );
 }
 
+timeval startTime, endTime;
+double elapsedTime;
+char timenow[100];
+
+char *now()
+{
+  timeval endTime;
+  double elapsedTime;
+
+  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.;
+  sprintf(timenow, "%10.3f", elapsedTime);
+  return(timenow);
+}
+
 //================================================
 
 class Universe {
 
 private:
 
-	double alpha_i;
-	double alpha_m;
-	double beta;
-	double gamma;
-	double delta;
+  double alpha_i;
+  double alpha_m;
+  double beta;
+  double gamma;
+  double delta;
 
-    double TargetNovelty;
-    double CumulativeRelativeLoss;
-    double CRLsquare;
-    string id;
+  double TargetNovelty;
+  double CumulativeRelativeLoss;
+  double CRLsquare;
+  string id;
 
 
-	int N_nodes;
-	int M_edges;
+  int N_nodes;
+  int M_edges;
 
-	int N_epochs;
-	int N_steps;
-	int N_repeats;
+  int N_epochs;
+  int N_steps;
+  int N_repeats;
 
-	int current_epoch;
-	double current_loss;
-	int current_repeat;
-    double current_novelty;
+  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
+  int mode_identify_failed;
+  int verbose_level; // 0 is silent, higher is more
 
-	double k_max;
+  double k_max;
 
-	graph_t Full_g;
+  graph_t Full_g;
 
 #define GraphRes float // was double
 
-	GraphRes **Prob;
-	GraphRes **Tried;
-	GraphRes **Dist;
-	GraphRes **Final;
-        GraphRes **EdgeIndex;
-	GraphRes *Rank;
+  GraphRes **Prob;
+  GraphRes **Tried;
+  GraphRes **Dist;
+  GraphRes **Final;
+  GraphRes **EdgeIndex;
+  GraphRes *Rank;
 
-    base_generator_type generator;
-    boost::uniform_real<> uni_dist;
-    boost::geometric_distribution<double> geo;
+  double max_loss;
 
+  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;
 
-	//======  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;
+    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;
 
-		std::ifstream inFile;
-        //string line;
+    TargetNovelty = target;
+    CumulativeRelativeLoss = 0.;
+    CRLsquare = 0.;
 
-        //-------------------------------
+    N_epochs  = Epochs;
+    N_steps   = Steps;
+    N_repeats = Repeats;
 
-        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;
+    current_epoch = 0;
+    current_loss = 0.;
+    current_repeat = 0;
 
-        //--------------------------------
+    id = idd;
 
-		int i, k;
-        int x, y;
-		Edge* edge_array_mine;
-		int num_arcs_mine, num_nodes_mine;
-		int* weights_mine;
+    verbose_level = 1;
+    mode_identify_failed = identify_failed;
 
-        TargetNovelty = target;
-        CumulativeRelativeLoss = 0.;
-        CRLsquare = 0.;
+    // 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;
+      }
+    }
 
-		N_epochs  = Epochs;
-		N_steps   = Steps;
-		N_repeats = Repeats;
+    i=0;
+    std::string line;
+    //while (! inFile.eof() && ! inFile.fail()) {
+    while (1==1) {
 
-		current_epoch = 0;
-		current_loss = 0.;
-		current_repeat = 0;
+      inFile >> x;
+      inFile >> y;
 
-        id = idd;
+      if (verbose_level > 2){
+	std::cout << " x: " << x;
+	std::cout << " y: " << y << std::endl;
+      }
 
-        verbose_level = 1;
+      if (i==0){
+	N_nodes=x;
+	M_edges=y;
+	break;
+      }
+      i++;
+    }
+    inFile.close();
 
-		mode_identify_failed = identify_failed;
+    if (verbose_level == 2){
+      std::cout << N_nodes <<  " nodes, " << M_edges << " edges"<<std::endl;
+    }
 
+    // k_max is the longest distance possible
 
-		//-------------------------------
-		// 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 {
+    //k_max = M_edges;
+    k_max = 70;
 
-            if (verbose_level > 2){
-                std::cout <<  " Opened <" << FileToOpen << ">"<<std::endl;
-            }
-        }
+    //------------------------------------
+    // Get memory allocated for all class members
 
-		i=0;
-        std::string line;
-		//while (! inFile.eof() && ! inFile.fail()) {
-        while (1==1) {
+    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);
 
-            inFile >> x;
-            inFile >> y;
+    //The second pass through file with the graph
 
-            if (verbose_level > 2){
-                std::cout << " x: " << x;
-                std::cout << " y: " << y << std::endl;
-            }
+    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;
+      }
+    }
 
-			if (i==0){
-				N_nodes=x;
-				M_edges=y;
-                break;
-			}
-			i++;
+    // 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++;
+	}
+      }
+    }
 
-		}
-		inFile.close();
+    // create graph -- hopefully, we can keep it, just modifying edge weights
 
-        if (verbose_level == 2){
-            std::cout << N_nodes <<  " nodes, " << M_edges << " edges"<<std::endl;
-        }
+    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_max is the longest distance possible
+    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);
 
-        //k_max = M_edges;
-		k_max = 70;
+    Full_g = g;
+    delete edge_array_mine;
+    delete weights_mine;
 
-		//------------------------------------
-		// Get memory allocated for all class members
+    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;
 
-		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);
+    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);
 
-        //The second pass through file with the graph
+	}
+      }
+    }
+    set_max_loss();
+  }
 
-		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;
-			}
-		}
+  void set_max_loss() {
+    // Based on formulas from random_analytical.pdf
+    // (9): E[RLt] = (1/t) * SUM(i=0 to T of: (V*(V-1) / (2 * (E-i))
+    // (12) Var[RLt] = 1 / (T**2) * SUM( i=0 to T of: ((V*(V-1))/2) - E + 1 ) / ((E-i)**2)
+    // (9) + 3sigma = E[RLt] + 3 * SQRT(Var)
 
+    int E = M_edges;
+    int V = N_nodes;
+    int T = TargetNovelty;
+    double sum, ERLt, VarRLt;
 
-		// Fill in the final graph -- and we are ready to go!
+    sum = 0;
+    for(int i=0;i<=T;i++) {
+      sum += (V*(V-1)) / (2*(E-i));
+    }
+    ERLt = (1.0/T) * sum;
 
-	    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);
-
-				}
-			}
-		}
-
-
+    sum = 0;
+    for(int i=0;i<=T;i++) {
+      sum += (((V*(V-1))/2.0) - E + 1.0 ) / ((E-i)*(E-i));
     }
+    VarRLt = (1.0 / (T*T)) * sum;
+    max_loss = ERLt + (3.0 * sqrt(VarRLt));
+    cout << "set_max_loss: V=" << V << " E=" << E << " T=" << T << "  ERLt=" << ERLt << " VarRLt=" << VarRLt << " max_loss=" << max_loss << endl;
+  }
 
+  int sample_failed_number(double pfail) {
 
-	//=====================================================================
-	int sample_failed_number(double pfail){
+    //boost::geometric_distribution<double> geo(pfail);
+    //boost::variate_generator<base_generator_type&, geometric_distribution<double> > geom(generator, geo);
 
-		//boost::geometric_distribution<double> geo(pfail);
-		//boost::variate_generator<base_generator_type&, geometric_distribution<double> > geom(generator, geo);
+    double r, u, g;
 
-		double r, u, g;
+    r=0.;
+    for(int i=0; i<N_steps; i++){
 
-        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)));
 
-            u=(double)rand();
-            u = 1.-u /(double)(pow(2.,31)-1.);
-            g=(int)(ceil(log(u) / log(pfail)));
+      //r += geom();
 
-			//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;
+      r+=g;
     }
 
-    //=============================================
-    void set_target(double target){
-        TargetNovelty=target;
+    if (verbose_level>=3){
+      std::cout << id << " failed " << r << std::endl;
     }
+    return (int) round(r);  // FIXME: Andrey: please verify that round() is correct.
 
-	//=============================================
-	int sample(){
+  }
 
-        //boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
-        // double r = uni(), Summa = 0.;
+  double get_target(void){
+    return TargetNovelty;
+  }
 
+  void set_target(double target){
+    TargetNovelty=target;
+  }
 
+  int sample(){
 
-        double r = rand(), Summa = 0.;
-        r /= (double)(pow(2.,31)-1.);
-		int result = 0;
-		int finished = 0;
+    //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;
-        }
+    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;
-
+    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;
+  void update_current_graph(void){
 
-		//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);
-					}
+    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;
+    }
+  }
 
+  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.;
-					}
-				}
-			}
-
-		}
-
-
+    // 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.;
-		}
+  void update_ranks(void){
 
-		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]++;
-				}
-			}
-		}
+    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 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;
+  //====================================================================
+  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);
 
-				}
-			}
-		}
+    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.;
-			}
-		}
-	}
+    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;
+  //==============================================
+  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";
+  //===============================================
+  string file_name(){
 
-        return s.str();
+    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 set_verbose(int verbose){
 
+    verbose_level = verbose;
+  }
 
-    //=============================================================
+
+  //=============================================================
   void update_probabilities(void){
 
 
@@ -807,13 +757,11 @@
     // Compute sampling probabilities
     // first pass: \xi_i,j
     int i, j;
-    //#pragma omp parallel for private (j)
-    //#pragma omp parallel for 
 
-    //    #pragma omp parallel for default(none) shared( Prob, alpha_i, alpha_m, beta, k_max, gamma, delta, Rank, Dist) private(i,j)
-    #pragma omp parallel for private(i,j)
+    std::cout << now() << " pass1 update_probabilities" << std::endl;
 
     for( i=0; i<N_nodes-1; i++){
+#pragma omp parallel for private(j)
       for(  j=i+1; j<N_nodes; j++){
 
 	double bg = 0.;
@@ -838,686 +786,441 @@
       }
     }
 
+    std::cout << now() << " pass2 update_probabilities" << std::endl;
 
     // second pass: sum
     double Summa = 0.;
 
-    #pragma omp parallel
+// pragma omp parallel
     for(i=0; i<N_nodes-1; i++){
       for(j=i+1; j<N_nodes; j++){
 	Summa += Prob[i][j];
       }
     }
 
+    std::cout << now() << " pass3 update_probabilities" << std::endl;
+
     // third pass: normalize
-    #pragma omp parallel 
     for( i=0; i<N_nodes-1; i++){
-      for(  j=i+1; j<N_nodes; j++){
+#pragma omp parallel for private(j)
+      for( j=i+1; j<N_nodes; j++){
 	Prob[i][j] /= Summa;
       }
     }
 
+    std::cout << now() << " exit  update_probabilities" << std::endl;
   }
 
-	// Now we are ready for simulations
-	//==============================================
-	void update_world(){
+  // Now we are ready for simulations
+  //==============================================
+  void update_world(){
 
-		int failed = 0;
+    int failed = 0;
 
-		// Given current universe compute shortest paths
-		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Given current universe compute shortest paths
+    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-std::cout << "Before updates" << std::endl;
-		update_current_graph();
-std::cout << "after current" << std::endl;
-		update_ranks();
-std::cout << "after ranks" << std::endl;
-		update_distances();
-std::cout << "after distances" << std::endl;
-		update_probabilities();
-std::cout << "after probs" << std::endl;
+    std::cout << now() << " Before updates" << std::endl;
+    update_current_graph();
+    std::cout << now() << " after current" << std::endl;
+    update_ranks();
+    std::cout << now() << " after ranks" << std::endl;
+    update_distances();
+    std::cout << now() << " after distances" << std::endl;
+    update_probabilities();
+    std::cout << now() << " after probs" << std::endl;
 
-		//===============================
-		// sampling
-		int result;
-		double cost=0., novel=0.;
-		int publishable = 0;
+    // sampling
+    int result;
+    double cost=0., novel=0.;
+    int publishable = 0;
 
+    if (mode_identify_failed == 1){
+      while(publishable < N_steps){
+	std::cout << now() << " failed:1 before sample\n";
+	result = sample();
+	std::cout << now() << " failed:1 after sample\n";
+	publishable += result;
+	failed += (1-result);
+      }
 
-		//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-		if (mode_identify_failed == 1){
+      std::cout << now() << " failed:1 before for\n";
+      for(int i=0; i<N_nodes-1; i++){
+	int j;
+        #pragma omp parallel for private(j)
+	for( j=i+1; j<N_nodes; j++){
 
-			while(publishable < N_steps){
+	  cost+=Tried[i][j];
 
-std::cout << "failed:1 before sample\n";
-		    	result = sample();
-std::cout << "failed:1 after sample\n";
-			    publishable += result;
-			    failed += (1-result);
+	  if (Tried[i][j]>0. && Final[i][j]>0.){
+	    novel+=1.;
+	    std::cout << now() << " failed:1 novel: " << novel << std::endl;
+	  }
+	}
+      }
+      std::cout << now() << " failed:1 after for\n";
+    }
+    else {
+      double pfail=0.;
+      int n_failed;
+      //, n_check = 0;
 
-			}
-
-std::cout << "failed:1 before for\n";
-			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.;
-std::cout << "failed:1 novel: " << novel << std::endl;
-					}
-				}
-			}
-std::cout << "failed:1 after for\n";
-
-		}
-		//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-		else {
-
-			double pfail=0.;
-			int n_failed;
-			//, n_check = 0;
-
-std::cout << "failed:0 before for\n";
-			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;
-			}
-
-std::cout << "failed:0 after for\n";
-			n_failed = sample_failed_number(pfail);
-
-
-std::cout << "failed:0 before while\n";
-			while(publishable < N_steps){
-
-		    	result = sample();
-			    publishable += result;
-			}
-std::cout << "failed:0 after while\n";
-
-
-			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.;
-std::cout << "failed:0 novel: " << novel << std::endl;
-
-					}
-				}
-			}
-		}
-
-        current_novelty = novel;
-
-
-		//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-verbose_level = 2;
-		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++;
+      std::cout << now() << " failed:0 before for\n";
+      for (int i=0; i<N_nodes-1; i++) {
+	int j;
+	// #pragma omp parallel for private(j)
+	for (int j=i+1; j<N_nodes; j++) {
+	  if (Final[i][j] == 0.) {
+	    pfail += Prob[i][j];
+	    Prob[i][j] = 0.;
+	  }
 	}
+      }
 
-
-	//======  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);
+      for(int i=0; i<N_nodes-1; i++){
+	int j;
+	// #pragma omp parallel for private(j)
+	for(j=i+1; j<N_nodes; j++){
+	  Prob[i][j] /= (1.-pfail);
 	}
+	//std::cout << std::endl;
+      }
 
-    //================================================
-    // Allocate memory
-    double** allocate_2Dmatrix_double(int N, int M)
-    {
-        double **pointer;
+      std::cout << now() << " failed:0 after for\n";
+      n_failed = sample_failed_number(pfail);
+      std::cout << now() << " failed:0 before while\n";
+      while(publishable < N_steps){
+	result = sample();
+	publishable += result;
+      }
+      std::cout << now() << " failed:0 after while\n";
 
-        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];
+      current_loss += (n_failed + N_steps);
+      cost = current_loss;
 
-        return pointer;
+      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.;
+	    std::cout << now() << " failed:0 novel: " << novel << std::endl;
+	  }
+	}
+      }
     }
-    // Allocate memory
-    GraphRes** allocate_2Dmatrix(int N, int M)
-    {
-        GraphRes **pointer;
+    current_novelty = novel;
 
-        if (verbose_level == 2){
-            std::cout<< "["<<N<<"|"<<M<<"]"<<std::endl;
-        }
-        pointer = new GraphRes*[N];
-        for (int i = 0; i < N; ++i)
-            pointer[i] = new GraphRes[M];
-
-        return pointer;
+    verbose_level = 2;
+    if (verbose_level == 2) {
+      std::cout << (current_repeat+1) << "  epoch=" << (current_epoch+1)
+		<< "  cost=" << cost
+		<< " novel=" << novel
+		<< " rel_loss=" << cost/novel
+		<< std::endl;
     }
-    //===================
-    GraphRes* allocate_1Dmatrix(int N)
-    {
-        GraphRes *pointer;
+    current_epoch++;
+  }
 
-        if(N > 0){
+  //======  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);
+  }
 
-            pointer = new GraphRes[N];
-
-        }else {
-
-            pointer = NULL;
-        }
-
-        return pointer;
-
+  // Allocate memory
+  double** allocate_2Dmatrix_double(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;
+  }
+  
+  // Allocate memory
+  GraphRes** allocate_2Dmatrix(int N, int M) {
+    GraphRes **pointer;
 
-    //==============================================
-    // De-Allocate memory to prevent memory leak
-    void delete_2Dmatrix_double(double **pointer, int N){
-
-        if (pointer != NULL){
-
-            for (int i = 0; i < N; ++i){
-                delete [] pointer[i];
-            }
-            delete [] pointer;
-        }
+    if (verbose_level == 2){
+      std::cout<< "["<<N<<"|"<<M<<"]"<<std::endl;
     }
-    void delete_2Dmatrix(GraphRes **pointer, int N){
+    pointer = new GraphRes*[N];
+    for (int i = 0; i < N; ++i)
+      pointer[i] = new GraphRes[M];
 
-        if (pointer != NULL){
+    return pointer;
+  }
 
-            for (int i = 0; i < N; ++i){
-                delete [] pointer[i];
-            }
-            delete [] pointer;
-        }
+  GraphRes* allocate_1Dmatrix(int N) {
+    GraphRes *pointer;
+    if(N > 0){
+      pointer = new GraphRes[N];
+    }else {
+      pointer = NULL;
     }
-    //====================
-    void delete_1Dmatrix(GraphRes *pointer){
+    return pointer;
+  }
 
-        delete [] pointer;
-    }
+  // De-Allocate memory to prevent memory leak
 
-    //===========================================
-    double get_rel_loss(){
-
-        return CumulativeRelativeLoss ;
+  void delete_2Dmatrix_double(double **pointer, int N) {
+    if (pointer != NULL){
+      for (int i = 0; i < N; ++i){
+	delete [] pointer[i];
+      }
+      delete [] pointer;
     }
+  }
 
-    //===========================================
-    double get_rel_loss_err(){
+  void delete_2Dmatrix(GraphRes **pointer, int N) {
+    if (pointer != NULL){
 
-        return CRLsquare ;
+      for (int i = 0; i < N; ++i){
+	delete [] pointer[i];
+      }
+      delete [] pointer;
     }
+  }
 
+  void delete_1Dmatrix(GraphRes *pointer) {
+    delete [] pointer;
+  }
 
+  double get_rel_loss() {
+    return CumulativeRelativeLoss ;
+  }
 
-    //==================================================================================
-    void evolve_to_target_and_save(int istart, int iend, double* storage, int* counters){
+  double get_rel_loss_err() {
+    return CRLsquare ;
+  }
 
-        double ALOT=100000000000.;
-        double ratio = 0.0;
-        int check = 0;
+  void evolve_to_target_and_save(int istart, int iend, double* storage, int* counters) {
+    double ALOT=100000000000.;
+    double ratio = 0.0;
 
-        reset_world();
+    reset_world();
+    for (int k = istart; k < iend; k++){
+      for(int i=0; i< N_epochs &&  current_novelty < TargetNovelty; i++){
+	if(ratio < max_loss) {
+	  update_world();
+	}
+	else{
+	  break;
+	}
+	ratio = current_loss/current_novelty;
+      }
+      if (ratio < max_loss) {
+	storage[k]=current_loss/current_novelty;
+      }
+      else{
+	storage[k]=ALOT;
+      }
+      counters[k]=1;
+      reset_world();
+    }
+  }
 
-        for (int k = istart; k < iend; k++){
+  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 */
+  }
 
-            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 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 */
     }
+  }
 
-    //==============================================
-    double Omax_loss_value(int tn){
+  int get_reruns(void) {
+    return N_repeats;
+  }
 
-        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 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.;
     }
-    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 */
+  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);
 
-    //==============================================
-    int get_reruns(void){
-        return N_repeats;
+    if(verbose_level==1){
+      std::cout << std::endl;
     }
-
-    //==============================================
-    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.;
-        }
+    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){
 
-    //==============================================
-    void evolve_to_target(){
+    if (position < 0 || position > 4) {return 0;}
 
-        reset_world();
-        if (beta < -1. || gamma < -1.){
-            CumulativeRelativeLoss = 100000000000.;
-            CRLsquare = 0.;
-            return;
-        }
+    else {
 
-
-        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));
-
+      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;
+      }
     }
-
-
-    //================================================================
-    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
-
+    return 0;
+  }
 };
 
-//============================================================
-
-#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;
 
-    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.;
 
-    double Loss=0., LossSquare=0.;
+  // timeval startTime, endTime;
+  // double elapsedTime;
+  gettimeofday(&startTime, NULL);
 
-    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);
+  }
 
-    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);
     }
+  }
 
-    for(int i=0; i<N; i++){
-        for(int j=0; j<5; j++){
-            un[i]->set_parameter(params[j],j);
-        }
-    }
+  int i;
 
-#ifdef P_DISPATCH
-    for(int i=0; i<Nworkers; i++){
+  // Execute rerun loop in parallel
 
-        dispatch_group_async(group, CustomQueues[i], ^{
+  //#pragma omp parallel for private (i)
+  for(i=0; i<N; i++){
+    un[i]->evolve_to_target_and_save(i, i+1, Results, Counters);
+  }
 
-            // 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);
-        });
+  for (int i=0; i<N; i++){
+    Loss+=Results[i]/(double)N;
+    LossSquare+=Results[i]*Results[i]/(double)N;
+  }
 
-        std::cout << "queued: i=" << i << " N=" << N << " istart=" << istart << " iend=" << iend << "\n";
-        // istart += step;
-        // iend = min(istart+step,N);
+  double two_std = ((LossSquare - Loss*Loss)/(double)N);
 
-    }
-    dispatch_group_wait(group, DISPATCH_TIME_FOREVER);
-    //dispatch_release(group);
-#else
-    int i;
+  two_std = 2.*sqrt(two_std);
+  std::pair<double,double> Res;
+  Res.first=Loss;
+  Res.second=two_std;
 
-    // Execute rerun loop in parallel
+  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";
 
-    #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;
+  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,
@@ -1527,190 +1230,173 @@
                      int*    Counters,
                      double* params0,
                      double annealing_cycles){
-#endif
-    //.................................
-    // re-implement annealing
+  //.................................
+  // 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);
+  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
+  // 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];
+  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;
-    }
+  for(int i=0;i<5;i++){
+    dx[i] = starting_jump;
+  }
 
-    // establish the current value
-    std::pair<double,double>Res;
+  // 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;
+  Res = multi_loss(       un,               Results, Counters, x);
 
-    if ( operation == 'm' ) {
-      FILE *f;
-      int N = un[0]->get_reruns();
+  std::cout << "Returned from initial multi_loss:" << std::endl;
+  std::cout << Res.first << " +- " << Res.second << std::endl;
 
-      f = fopen("multi_loss.data","w");
-      for(int i=0; i<N; i++) {
-        fprintf(f,"%.20e\n",Results[i]);
-      }
-      fclose(f);
-      exit(0);
+  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;
+  curr_x   = Res.first;
+  curr_err = Res.second;
 
-    // optimization cycle
+  // optimization cycle
 
-    for(int i=0; i<annealing_cycles; i++){
+  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;
+    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){
+    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 k=0; k<5; k++){
-                rejection[k]/=(double)cycle;
+    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;
 
-                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;
-        }
+	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(un, Results, Counters, x);
+	std::cout << Res.first << " +- " << Res.second << std::endl;
 
-        for (int j=0; j<5; j++){
+	ratio = min(1.,exp(-(Res.first-curr_x)/temperature));
+	r = rand()/(double)(pow(2.,31)-1.);
+	std::cout << r << " vs " << ratio << std::endl;
 
-            ///////////////////////////////
-            if (FIX_VARIABLES==0 || var_fixed[j]==0){
+	double ALOT=100000000000.;
 
+	if (Res.first < ALOT)
+	{
+	  ofstream filestr;
 
+	  filestr.open ("best_opt_some.txt", ofstream::app);
 
-                // get new value of x[j]
-                double x_hold=x[j];
-                x_tmp = get_new_x(x[j],dx[j]);
-                x[j]=x_tmp;
+	  // >> i/o operations here <<
+	  filestr
 
-                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);
-                }
+	    << "N, " << i << ", " << j << ", " << dx[j] << ", " << rejection[j] << ", |, " // FIXME: MW-DEBUGGING
 
-#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;
+	    << 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";
 
-                ratio = min(1.,exp(-(Res.first-curr_x)/temperature));
-                r = rand()/(double)(pow(2.,31)-1.);
-                std::cout << r << " vs " << ratio << std::endl;
+	  filestr.close();
 
-                double ALOT=100000000000.;
 
-                if (Res.first < ALOT)
-                {
-                    ofstream filestr;
+	  filestr.open ("max_dist.txt", ofstream::app);
 
-                    filestr.open ("best_opt_some.txt", ofstream::app);
+	  // >> i/o operations here <<
+	  filestr << max_dist << ",\n";
 
-                    // >> i/o operations here <<
-                    filestr
+	  filestr.close();
 
-		    << "N, " << i << ", " << j << ", " << dx[j] << ", " << rejection[j] << ", |, " // FIXME: MW-DEBUGGING
+	  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);
+	}
 
-                    << 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";
+	if (r > ratio){
 
-                    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)
+	  std::cout << " "<< (i+1) << ","<< (j)
                     <<" "<< (i+1) << " Did not accept "
                     << x_tmp << "(" << j << ")" << std::endl;
-                    std::cout << un[0]->get_parameter(0)
+	  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);
-                    }
+	  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.;
-                }
+	  //set_parameter(x[j], j);
+	  rejection[j]+=1.;
+	}
 
-                else {
+	else {
 
-                    curr_x   = Res.first;
-                    curr_err = Res.second;
-                    x[j] = x_tmp;
+	  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);
-                    }
+	  for(int w=0; w<Nworkers; w++){
+	    un[w]->set_parameter(x[j], j);
+	  }
 
-                    std::cout << (i+1) << string_wrap((string) " Rejection counts: ", 8)
+	  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) << " "
@@ -1718,7 +1404,7 @@
                     << wrap_double(rejection[4],6) << " "
                     << std::endl << std::endl;
 
-                    std::cout << " "<< (i+1) <<","<< (j)
+	  std::cout << " "<< (i+1) <<","<< (j)
                     <<" "
                     << string_wrap((string) "***** Did accept! ", 3)
                     << wrap_double(un[0]->get_parameter(0),2) << " "
@@ -1727,238 +1413,176 @@
                     << 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"};
 
-    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};
 
-    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};
 
-    //          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];
+  static double* Results;
+  static int*    Counters;
 
-    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;
 
-    timeval t1, t2;
-    double elapsedTime;
-    // start timer
-    gettimeofday(&t1, NULL);
+  gettimeofday(&t1, NULL);   // start timer
 
-
-    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);
+  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;
+      }
     }
-    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));
+  }
 
-    for (int j=0; j<5; j++){
+  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;
+  }
 
-        cout << j << " | " << var_fixed[j] << " (fixed) \n";
-    }
+  //srand(time(0));
+  //srandomdev();
 
-    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];
+  if ( initSeed != 0.0 ) {
+    srand(initSeed);
+  }
+  else {
+    timeval t;
+    gettimeofday(&t, NULL);
+    srand(t.tv_usec);
+  }
 
-    //...............................
-
-    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
+  {
+    double r=0;
+    for (int j=0; j<100; j++){
+      r = rand()/(double)(pow(2.,31)-1.);
+      std::cout << r << " ";
     }
+    std::cout << "\n";
+  }
 
-    //...............................
-    if(n_rep > 0){
+  //random initiation of starting parameters
 
-        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;
-
+  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;
+      }
     }
-    //...............................
-    //srand(time(0));
-    //srandomdev();
+  }
 
-    if ( initSeed != 0.0 ) {
-      srand(initSeed);
-    }
-    else {
-        timeval t;
-        gettimeofday(&t, NULL);
-        srand(t.tv_usec);
-    }
+  double T_start=params2[0], T_end=params2[1], Target_rejection=params2[3], starting_jump=params2[4];
+  int Annealing_repeats = (int) params2[2];
 
-    {
-        double r=0;
-        for (int j=0; j<100; j++){
+  // print_memory_status();
 
+  multi_annealing(un, T_start, T_end, Target_rejection, Annealing_repeats,
+		  starting_jump, Results, Counters, params0, Annealing_repeats);
 
+  // stop timer
 
-            r = rand()/(double)(pow(2.,31)-1.);
-            std::cout << r << " ";
-        }
-        std::cout << "\n";
-    }
-  	//random initiation of starting parameters
+  gettimeofday(&t2, NULL);
 
-    if (range > 0.){
+  // compute and print the elapsed time in millisec
 
-        for (int i=0; i < 5; i++){
+  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";
 
-            if (params0[i]==-100.){
+  for(int i=0; i<Nworkers; i++){
+    delete un[i];
+  }
 
-                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];
-
-    // print_memory_status();
-
-#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;
-
-
-
+  if(n_rep > 0){
+    delete [] Results;
+    delete [] Counters;
+  }
+  return 0;
 }
 

Added: SwiftApps/SciColSim/maxloss.cpp
===================================================================
--- SwiftApps/SciColSim/maxloss.cpp	                        (rev 0)
+++ SwiftApps/SciColSim/maxloss.cpp	2012-12-23 04:33:21 UTC (rev 6110)
@@ -0,0 +1,101 @@
+#include <fstream>
+#include <iostream>
+#include <stdio.h>
+#include <time.h>
+#include <ctime>
+#include <algorithm>
+#include <string>
+
+#include <sys/param.h>
+#include <sys/time.h>
+#include <sys/types.h>
+#include <unistd.h>
+#include <stdlib.h>
+
+#include <math.h>
+
+using namespace std;
+
+void set_max_loss(long long v, long long e, long long t) {
+
+  // Based on formulas from random_analytical.pdf
+
+  // t-=100;
+  double E = e;
+  double V = v;
+  double T = t;
+
+  double max_loss;
+  double sum, ERLt, VarRLt, absVarRLt;
+
+  // (9): E[RLt] = (1/t) * SUM(i=0 to T of: (V*(V-1) / (2 * (E-i))
+
+  sum = 0;
+  for(int i=0;i<t;i++) {
+    sum += ((double)(v*(v-1))) / ((double)(2*(e-i)));
+  }
+  ERLt = (1.0/(double) t) * sum;
+
+  // (12): Var[RLt] = 1 / (T**2) * SUM( i=0 to T of: ((V*(V-1))/2) - E + 1 ) / ((E-i)**2)
+
+  sum = 0;
+  for(int i=0;i<t;i++) {
+    sum +=  
+      ( (((double)(v*(v-1)))/ 2.0) - (double)e + (double)i )
+      /  
+      (double)((e-i)*(e-i))
+      ;
+    if ( (double)((e-i)*(e-i)) == 0.0 ) {
+      printf( "zero denom: %d %d\n", e, i);
+    }
+  }
+  int tmp = t*t;
+  VarRLt = (1.0/(double)((t) * t)) * sum;
+  if (VarRLt < 0.0) absVarRLt=0.0;
+  else absVarRLt = VarRLt;
+
+  // (9) + 3sigma = E[RLt] + 3 * SQRT(Var)
+  max_loss = ERLt + (3.0 * sqrt(absVarRLt));
+  printf("%12d %12.3f %12.8e %12.3f\n", t, ERLt, VarRLt, max_loss);
+}
+
+struct Graph {
+  const char *n;
+  int v;
+  int e;
+} G[] =
+
+{
+  "movie",    500,    1008,
+  "big.001", 1857,    1312,
+  "big.002", 3094,    2683,
+  "big.005", 5405,    6776,
+  "big.01",  7765,   13421,
+  "big.10", 19281,  134594,
+  "big.20", 23258,  269133,
+  "big.30", 25484,  403401,
+  "big.40", 26821,  537106,
+  "big.50", 27774,  671179,
+  "big.60", 28482,  804774,
+  "big.70", 29055,  938061,
+  "big.80", 29487, 1072415,
+  "big.90", 29779, 1205741,
+  "big",    30060, 1338753,   
+  "",        0,     0
+};
+
+int main()
+{
+  int E, V, T;
+  int i=0;
+  while (G[i].v > 0) {
+    cout << "\nGraph " << G[i].n << " V=" << G[i].v << " E=" << G[i].e << endl;
+    printf("Pct          T          ERLt       VarRLt       max_loss\n");
+    for(double p = 0.1; p<=1.0; p+=0.1) {
+      printf("%.2f",p);
+      set_max_loss(G[i].v, G[i].e, (int)((G[i].e * p)+0.5));
+    }
+    i++;
+  }
+}
+

Added: SwiftApps/SciColSim/maxloss.py
===================================================================
--- SwiftApps/SciColSim/maxloss.py	                        (rev 0)
+++ SwiftApps/SciColSim/maxloss.py	2012-12-23 04:33:21 UTC (rev 6110)
@@ -0,0 +1,39 @@
+#! /usr/bin/env python
+
+n = 0
+v = 1
+e = 2
+
+graphs = [ 
+     [ "movie",    500,    1008 ],
+     [ "big.001", 1857,    1312 ],
+     [ "big.002", 3094,    2683 ],
+     [ "big.005", 5405,    6776 ],
+     [ "big.01",  7765,   13421 ],
+     [ "big.10", 19281,  134594 ],
+     [ "big.20", 23258,  269133 ],
+     [ "big.30", 25484,  403401 ],
+     [ "big.40", 26821,  537106 ],
+     [ "big.50", 27774,  671179 ],
+     [ "big.60", 28482,  804774 ],
+     [ "big.70", 29055,  938061 ],
+     [ "big.80", 29487, 1072415 ],
+     [ "big.90", 29779, 1205741 ],
+     [ "big",    30060, 1338753 ], 
+     ]
+
+for G in graphs:
+	N = G[n]
+	V = G[v]
+	E = G[e]
+	print "\nGraph %s: V=%d E=%d\n" % (N,V,E)
+	for pct in range(1,11):
+		Tmax = int(E * (pct/10.0))
+		tl = 0.
+		vl = 0.
+		for t in range(Tmax):
+			tl = tl + float(V*(V-1.))/float(2.*(E-t))
+			vl = vl + float(V*(V-1.)/2. - E + t)/float((E-t)*(E-t))
+		tl /= float(Tmax)
+		vl /= float(Tmax)*float(Tmax)
+		print pct, Tmax, tl, vl


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

Added: SwiftApps/SciColSim/maxloss_ar.py
===================================================================
--- SwiftApps/SciColSim/maxloss_ar.py	                        (rev 0)
+++ SwiftApps/SciColSim/maxloss_ar.py	2012-12-23 04:33:21 UTC (rev 6110)
@@ -0,0 +1,15 @@
+V = [500]
+E = [1008]
+
+curr = 0
+
+
+for Tmax in range(101,1009,101):
+        tl = 0.
+        vl = 0.
+        for t in range(Tmax):
+                tl = tl + float(V[0]*(V[0]-1.))/float(2.*(E[0]-t))
+                vl = vl + float(V[0]*(V[0]-1.)/2. - E[0] + t)/float((E[0]-t)*(E[0]-t))
+        tl /= float(Tmax)
+        vl /= float(Tmax)*float(Tmax)
+        print Tmax, tl, vl

Modified: SwiftApps/SciColSim/samplegraph.sh
===================================================================
--- SwiftApps/SciColSim/samplegraph.sh	2012-12-19 14:06:35 UTC (rev 6109)
+++ SwiftApps/SciColSim/samplegraph.sh	2012-12-23 04:33:21 UTC (rev 6110)
@@ -2,7 +2,7 @@
 
 fraction=${1:-1.0}
 
-tmp=$(mktemp -t samplegraph)
+tmp=$(mktemp -t samplegraph.XXXX)
 
 awk '
 
@@ -16,6 +16,7 @@
      
 BEGIN { 
   getmsec="perl -MTime::HiRes=gettimeofday -e \"print int(1000*gettimeofday()).qq(\\n);\""
+  getmsec="date +%N"
   getmsec | getline msec
   close(getmsec)
   _cliff_seed = (msec % 10000) / 10000

Modified: SwiftApps/SciColSim/testgraph.py
===================================================================
--- SwiftApps/SciColSim/testgraph.py	2012-12-19 14:06:35 UTC (rev 6109)
+++ SwiftApps/SciColSim/testgraph.py	2012-12-23 04:33:21 UTC (rev 6110)
@@ -138,12 +138,14 @@
   endTarget        = startTarget+1
   incrTarget       = 50
   optimizerRepeats = 1
+  # evolveReruns     = 1
   evolveReruns     = 1
   annealingSteps   = 1
-  NWorkers         = "2"
+  NWorkers         = "8"
   openmp           = "OMP_NUM_THREADS=" + NWorkers
   operation        = "m"  # n=normal, m=manual (runs 1 multi_loss call)
   seed             = "" # "123456"
+  app              = "./openmp-optimizer";
   app              = "./graphsim";
 
 # Ensure we dont pass new parameters to original optimizer versions
@@ -158,9 +160,7 @@
 
 for target in range(startTarget,endTarget,incrTarget):
   for i in range(optimizerRepeats):
-    args = "rm -f bestdb.txt; " + \
-           openmp + " " + "/usr/bin/time -l " + \
-                          app + " 0 0 0 0 0 " + str(target) + " 40000 20 " + str(evolveReruns) + \
+    args = "rm -f bestdb.txt;" +  openmp +" " + "/usr/bin/time --verbose " + app + " 0 0 0 0 0 " + str(target) + " 40000 20 " + str(evolveReruns) + \
            " 2 2 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")




More information about the Swift-commit mailing list