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