[Swift-commit] r5705 - SwiftApps/SciColSim
jonmon at ci.uchicago.edu
jonmon at ci.uchicago.edu
Wed Mar 7 10:06:18 CST 2012
Author: jonmon
Date: 2012-03-07 10:06:17 -0600 (Wed, 07 Mar 2012)
New Revision: 5705
Modified:
SwiftApps/SciColSim/optimizer.cpp
Log:
o added loss limiting code in the optimizer for the SciColSim app. Needs larger scale testing.
Modified: SwiftApps/SciColSim/optimizer.cpp
===================================================================
--- SwiftApps/SciColSim/optimizer.cpp 2012-03-02 00:46:25 UTC (rev 5704)
+++ SwiftApps/SciColSim/optimizer.cpp 2012-03-07 16:06:17 UTC (rev 5705)
@@ -22,7 +22,7 @@
#include <iostream>
#include <stdio.h>
#include <time.h>
-#include <ctime>
+#include <ctime>
#include <algorithm>
#include <string>
@@ -69,6 +69,8 @@
#define FIX_VARIABLES 1
+#define MAX_LOSS 160
+
using namespace boost;
using namespace std;
using namespace boost::numeric::ublas;
@@ -114,10 +116,10 @@
GaussNum += ((double)rand()/(double)RAND_MAX - 0.5);
}
GaussNum = GaussNum*sqrt((double)12/(double)NumInSum);
-
-
+
+
return GaussNum;
-
+
}
@@ -133,27 +135,27 @@
//================================================
//================================================================
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){
+
+ if (r > 0.5){
new_x = x + rand()*dx/(double)(pow(2.,31)-1.);
- } else {
+ } 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";
@@ -191,143 +193,143 @@
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
+const
string i2string(int i){
-
+
std::ostringstream s;
- s << "worker"
+ s << "worker"
<< lexical_cast<std::string>(i);
-
+
return s.str();
-
+
}
//===============================================
char* i2char(int i){
-
+
std::ostringstream s;
- s << "worker"
+ 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;
+
+ 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());
@@ -335,57 +337,57 @@
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;
+ 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) {
@@ -396,41 +398,41 @@
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;
+
+ 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();
-
+ inFile.close();
+
k=0;
for (int i=0; i<N_nodes-1; i++){
for (int j=i+1;j<N_nodes; j++){
@@ -440,19 +442,19 @@
}
}
}
-
-
-
+
+
+
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// 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++){
@@ -464,17 +466,17 @@
}
}
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.){
@@ -482,103 +484,103 @@
v = vertex(j, Full_g);
remove_edge(u,v,Full_g);
remove_edge(v,u,Full_g);
-
+
}
}
}
-
-
+
+
}
-
-
+
+
//=====================================================================
int sample_failed_number(double pfail){
-
+
//boost::geometric_distribution<double> geo(pfail);
//boost::variate_generator<base_generator_type&, geometric_distribution<double> > geom(generator, geo);
-
+
double r, u, g;
-
+
r=0.;
for(int i=0; i<N_steps; i++){
-
+
u=(double)rand();
u = 1.-u /(double)(pow(2.,31)-1.);
g=(int)(ceil(log(u) / log(pfail)));
-
+
//r += geom();
-
+
r+=g;
}
-
+
if (verbose_level>=3){
std::cout << id << " failed " << r << std::endl;
}
return (int) round(r); // FIXME: Andrey: please verify that round() is correct.
-
+
}
-
+
//=============================================
double get_target(void){
return TargetNovelty;
}
-
+
//=============================================
void set_target(double target){
TargetNovelty=target;
}
-
+
//=============================================
int sample(){
-
+
//boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
// double r = uni(), Summa = 0.;
-
-
-
+
+
+
double r = rand(), Summa = 0.;
r /= (double)(pow(2.,31)-1.);
int result = 0;
int finished = 0;
-
+
if (verbose_level==4){
std::cout << id << " sampled " << r << std::endl;
}
-
- for(int i=0; i<N_nodes-1 && finished==0; i++){
+
+ for(int 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;
+ 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);
+ //s = edge(i, j, Full_g);
boost::graph_traits<graph_t>::edge_descriptor e1,e2;
bool found1, found2;
u = vertex(i, Full_g);
@@ -589,64 +591,64 @@
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);
+ s = vertex(j, Full_g);
dijkstra_shortest_paths(Full_g, s, predecessor_map(&p[0]).distance_map(&d[0]));
-
+
//std::cout <<" Vertex "<< j << std::endl;
graph_traits < graph_t >::vertex_iterator vi, vend;
-
+
for (boost::tie(vi, vend) = vertices(Full_g); vi != vend; ++vi) {
-
+
if (p[*vi]!=*vi){
Dist[*vi][j]=d[*vi];
Dist[j][*vi]=d[*vi];
-
+
if ( (int)round(Dist[*vi][j]>max_dist)) {
// FIXME: Andrey: please verify that (int) cast is correct. Do we need to round()?
- // also, the indent on this iff statement was way off -
+ // also, the indent on this iff statement was way off -
// perhaps due to space v. tab?
max_dist=(int)round(Dist[*vi][j]);
}
-
-
+
+
} else {
Dist[*vi][j]=-1.;
Dist[j][*vi]=-1.;
}
}
}
-
+
}
-
-
+
+
}
-
+
//======================================================
void update_ranks(void){
-
+
for(int i=0; i<N_nodes; i++){
Rank[i]=0.;
}
-
+
for(int i=0; i<N_nodes-1; i++){
for( int j=i+1; j<N_nodes; j++){
if (Tried[i][j]>0. && Final[i][j] >0.){
@@ -655,30 +657,30 @@
}
}
}
-
+
}
-
+
//====================================================================
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){
@@ -686,18 +688,18 @@
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) {
@@ -707,28 +709,28 @@
}
}
}
-
-
+
+
//==============================================
void show_parameters(void){
-
- std::cout << "Parameters: "
+
+ 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) << "_"
+ s << "world_"
+ << lexical_cast<std::string>(alpha_i) << "_"
<< lexical_cast<std::string>(alpha_m) << "_"
<< lexical_cast<std::string>(beta) << "_"
<< lexical_cast<std::string>(gamma) << "_"
@@ -736,228 +738,228 @@
<< 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.)) +
+
+ 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_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;
+ 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
+ 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);
+ 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];
-
+ 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];
}
@@ -966,54 +968,98 @@
}
//====================
void delete_1Dmatrix(double *pointer){
-
+
delete [] pointer;
}
-
+
//===========================================
double get_rel_loss(){
-
+
return CumulativeRelativeLoss ;
}
-
+
//===========================================
double get_rel_loss_err(){
-
+
return CRLsquare ;
}
-
-
-
+
+
+
//==================================================================================
void evolve_to_target_and_save(int istart, int iend, double* storage, int* counters){
-
- double ALOT=100000000000.;
-
+
+ double ALOT=10000.;
+ int check = 0;
+
reset_world();
-
+
for (int k = istart; k < iend; k++){
-
-
+
+// Code for getting time to check
+
for(int i=0; i< N_epochs && current_novelty < TargetNovelty; i++){
- update_world();
+
+// Check the walltime limit
+
+ if(current_loss/current_novelty < max_loss_value((int)TargetNovelty)){
+ update_world();
+ }
+ else{
+ check=1;
+ break;
+ }
}
-
- storage[k]=current_loss/current_novelty;
+
+ if( check ){
+ storage[k]=ALOT;
+ }
+ else{
+ storage[k]=current_loss/current_novelty;
+ }
counters[k]=1;
-
-
+
reset_world();
}
-
+
}
+
//==============================================
+ double max_loss_value(int tn){
+ /* Poor check */
+ /* Do not check for equality, check for within some error */
+ if( tn == 58. ) return 128.286;
+ else if( tn == 108. ) return 131.866;
+ else if( tn == 158. ) return 135.551;
+ else if( tn == 208. ) return 139.694;
+ else if( tn == 258. ) return 144.163;
+ else if( tn == 308. ) return 148.967;
+ else if( tn == 358. ) return 154.201;
+ else if( tn == 408. ) return 159.962;
+ else if( tn == 458. ) return 166.441;
+ else if( tn == 508. ) return 173.655;
+ else if( tn == 558. ) return 181.921;
+ else if( tn == 608. ) return 191.246;
+ else if( tn == 658. ) return 202.150;
+ else if( tn == 708. ) return 215.197;
+ else if( tn == 758. ) return 202.150; // Verify with Andrey. Same as TargetNovelty of 658
+ else if( tn == 808. ) return 251.698;
+ else if( tn == 858. ) return 279.201;
+ else if( tn == 908. ) return 320.112;
+ else if( tn == 958. ) return 394.774;
+ else if( tn == 1008. ) return 1052.38; // Huge number here. Why?
+ else return 10000.; /* Pray that this does not occur. Ask Andrey what to do */
+ }
+
+
+ //==============================================
int get_reruns(void){
return N_repeats;
}
-
+
//==============================================
double get_parameter(int i){
-
+
switch(i){
case 0:
return alpha_i;
@@ -1026,72 +1072,72 @@
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;
@@ -1109,17 +1155,17 @@
delta=value;
return 1;
}
-
+
}
-
+
return 0;
}
-
-#ifdef notdef
+
+#ifdef notdef
//=================================================================
- void try_annealing(double starting_jump, int iterations,
+ 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.};
@@ -1128,41 +1174,41 @@
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();
+ 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);
@@ -1175,104 +1221,104 @@
}
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();
-
+
+
+ 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 "
+
+ 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 << " "
+ << 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)
+ 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) << " "
+ << 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)
+
+ std::cout << string_wrap(id, 4) <<" "<< (i+1) <<","<< (j)
<<" "
- << string_wrap((string) "***** Did accept! ", 3)
- << wrap_double(alpha_i,2)
+ << 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) << " "
+ << wrap_double(beta,5) << " "
+ << wrap_double(gamma,9) << " "
+ << wrap_double(delta,6) << " "
<< std::endl << std::endl;
-
+
}
- //........................................................
-
+ //........................................................
+
}
-
+
}
-
+
}
-
+
#endif
-
+
};
//============================================================
#ifdef P_DISPATCH
-std::pair<double,double> multi_loss(dispatch_group_t group,
- Universe* un[],
+std::pair<double,double> multi_loss(dispatch_group_t group,
+ Universe* un[],
dispatch_queue_t* CustomQueues,
double* Results,
int* Counters,
double* params){
#else
-std::pair<double,double> multi_loss(Universe* un[],
+std::pair<double,double> multi_loss(Universe* un[],
double* Results,
int* Counters,
double* params){
#endif
-
+
int N = un[0]->get_reruns();
int step = (int)floor((double)N/(double)(Nworkers)); // FIXME: Andrey: please check change in cast grouping and use of floor
int istart=0;
int iend = istart+step;
-
+
double Loss=0., LossSquare=0.;
-
+
timeval startTime, endTime;
double elapsedTime;
gettimeofday(&startTime, NULL);
@@ -1282,20 +1328,20 @@
un[i]->set_parameter(params[j],j);
}
}
-
-#ifdef P_DISPATCH
+
+#ifdef P_DISPATCH
for(int i=0; i<Nworkers; i++){
-
+
dispatch_group_async(group, CustomQueues[i], ^{
-
+
// un[i]->evolve_to_target_and_save(istart, iend, Results, Counters);
un[i]->evolve_to_target_and_save(i*step, min((i+1)*step,N), Results, Counters);
});
-
+
std::cout << "queued: i=" << i << " N=" << N << " istart=" << istart << " iend=" << iend << "\n";
// istart += step;
// iend = min(istart+step,N);
-
+
}
dispatch_group_wait(group, DISPATCH_TIME_FOREVER);
//dispatch_release(group);
@@ -1316,23 +1362,23 @@
un[i]->evolve_to_target_and_save(i*step, min((i+1)*step,N), Results, Counters);
}
#endif
-
+
for (int i=0; i<N; i++){
-
+
Loss+=Results[i]/(double)N;
LossSquare+=Results[i]*Results[i]/(double)N;
-
+
// std::cout<<" " << 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
@@ -1340,30 +1386,30 @@
cout << "multi_loss(N=" << N << ", target=" << un[0]->get_target() << ") elapsed time: " << elapsedTime << " seconds " << elapsedTime/60. << " minutes\n\n";
return Res;
-
-
+
+
}
//============================================================
//============================================================
#ifdef P_DISPATCH
-void multi_annealing( dispatch_group_t group,
- Universe* un[],
- dispatch_queue_t* CustomQueues,
- double T_start, double T_end,
- double Target_rejection,
- int Annealing_repeats,
+void multi_annealing( dispatch_group_t group,
+ Universe* un[],
+ dispatch_queue_t* CustomQueues,
+ double T_start, double T_end,
+ double Target_rejection,
+ int Annealing_repeats,
double starting_jump,
double* Results,
int* Counters,
double* params0,
double annealing_cycles){
#else
-void multi_annealing(Universe* un[],
- double T_start, double T_end,
- double Target_rejection,
- int Annealing_repeats,
+void multi_annealing(Universe* un[],
+ double T_start, double T_end,
+ double Target_rejection,
+ int Annealing_repeats,
double starting_jump,
double* Results,
int* Counters,
@@ -1372,7 +1418,7 @@
#endif
//.................................
// re-implement annealing
-
+
double dx[5]={0.,0.,0.,0.,0};
double x[5]={0.,0.,0.,0.,0};
double rejection[5]={0., 0., 0., 0., 0.};
@@ -1381,30 +1427,30 @@
double ratio, r;
int cycle=10;
//boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(generator, uni_dist);
-
+
// set up parameter for annealing
-
+
x[0]=params0[0];
x[1]=params0[1];
x[2]=params0[2];
x[3]=params0[3];
x[4]=params0[4];
-
+
for(int i=0;i<5;i++){
dx[i] = starting_jump;
}
-
+
// establish the current value
std::pair<double,double>Res;
-#ifdef P_DISPATCH
+#ifdef P_DISPATCH
Res = multi_loss(group, un, CustomQueues, Results, Counters, x);
#else
Res = multi_loss( un, Results, Counters, x);
#endif
std::cout << "Returned from initial multi_loss:" << std::endl;
std::cout << Res.first << " +- " << Res.second << std::endl;
-
+
if ( operation == 'm' ) {
FILE *f;
int N = un[0]->get_reruns();
@@ -1419,19 +1465,19 @@
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.;
@@ -1443,143 +1489,143 @@
}
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";
+
+ std::cout << wrap_double(x_tmp,10) << " " << wrap_double(j,9) << "\n\n";
//=======================================
//.............................................
for(int w=0; w<Nworkers; w++){
un[w]->set_parameter(x_tmp, j);
}
-
-#ifdef P_DISPATCH
+
+#ifdef P_DISPATCH
Res = multi_loss(group, un, CustomQueues, Results, Counters, x);
#else
Res = multi_loss( un, Results, Counters, x);
#endif
std::cout << Res.first << " +- " << Res.second << std::endl;
-
+
ratio = min(1.,exp(-(Res.first-curr_x)/temperature));
r = rand()/(double)(pow(2.,31)-1.);
std::cout << r << " vs " << ratio << std::endl;
-
+
double ALOT=100000000000.;
-
+
if (Res.first < ALOT)
{
ofstream filestr;
-
+
filestr.open ("best_opt_some.txt", ofstream::app);
-
+
// >> i/o operations here <<
filestr
<< "N, " << i << ", " << j << ", " << dx[j] << ", " << rejection[j] << ", |, " // FIXME: MW-DEBUGGING
-
- << un[0]->get_target() << ", "
- << Res.first
- << ", " << un[0]->get_parameter(0)
- << ", " << un[0]->get_parameter(1)
- << ", " << un[0]->get_parameter(2)
- << ", " << un[0]->get_parameter(3)
+
+ << un[0]->get_target() << ", "
+ << Res.first
+ << ", " << un[0]->get_parameter(0)
+ << ", " << un[0]->get_parameter(1)
+ << ", " << un[0]->get_parameter(2)
+ << ", " << un[0]->get_parameter(3)
<< ", " << un[0]->get_parameter(4) << ", " << Res.second << ",\n";
-
+
filestr.close();
-
-
+
+
filestr.open ("max_dist.txt", ofstream::app);
-
+
// >> i/o operations here <<
filestr << max_dist << ",\n";
-
+
filestr.close();
FILE *bf;
bf = fopen("bestdb.txt","a");
fprintf(bf, "N %2d %2d %10.5f %5.2f | %5.2f %10.5f [ %5.2f %5.2f %10.5f %10.5f %10.5f ] %10.5f\n", i, j, dx[j], rejection[j],
un[0]->get_target(),
- Res.first,
- un[0]->get_parameter(0),
- un[0]->get_parameter(1),
- un[0]->get_parameter(2),
+ Res.first,
+ un[0]->get_parameter(0),
+ un[0]->get_parameter(1),
+ un[0]->get_parameter(2),
un[0]->get_parameter(3),
un[0]->get_parameter(4),
Res.second);
fclose(bf);
}
-
-
+
+
if (r > ratio){
-
- std::cout << " "<< (i+1) << ","<< (j)
- <<" "<< (i+1) << " Did not accept "
+
+ 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)
+ 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);
+
+
+ //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) << " "
+
+ 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) << " "
+ << wrap_double(rejection[2],5) << " "
+ << wrap_double(rejection[3],9) << " "
+ << wrap_double(rejection[4],6) << " "
<< std::endl << std::endl;
-
- std::cout << " "<< (i+1) <<","<< (j)
+
+ std::cout << " "<< (i+1) <<","<< (j)
<<" "
- << string_wrap((string) "***** Did accept! ", 3)
+ << 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) << " "
+ << 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;
-
-
-
+
+
+
}
- //........................................................
-
+ //........................................................
+
}
}
-
+
}
-
+
}
@@ -1588,43 +1634,43 @@
int
main(int argc, char* argv[])
{
-
+
double params0[6] = {0., 0., 0., 0., 0., 0.2}, target=50., range;
string par_names0[6] = {"alpha_i", "alpha_m", "beta", "gamma", "delta", "target"};
string par_names1[4] = {"n_epochs", "n_steps", "n_reruns", "range"};
string par_names2[5] = {"T_start", "T_end", "Annealing_steps","Target_rejection","Starting_jump"};
- string par_names3[5] = {"FREEZE_alpha_i", "FREEZE_alpha_m", "FREEZE_beta", "FREEZE_gamma", "FREEZE_delta"};
+ string par_names3[5] = {"FREEZE_alpha_i", "FREEZE_alpha_m", "FREEZE_beta", "FREEZE_gamma", "FREEZE_delta"};
string par_names4[3] = {"Operation", "Nworkers", "initSeed"};
int params1[4] = {300, 50, 1000, 10};
int params3[5] = { 0, 0, 0, 0, 0};
-
+
// temperature_start, temperature_end, annealing_steps target_rejection Starting_jump
double params2[5] = {1, 0.001, 100, 0.3, 1.5};
-
+
int verbose_level = 2;
const std::string one="one", two="two";
static Universe* un[MAXNworkers];
#ifdef P_DISPATCH
static dispatch_queue_t CustomQueues[MAXNworkers];
-#endif
+#endif
static double* Results;
static int* Counters;
-
+
timeval t1, t2;
double elapsedTime;
// start timer
gettimeofday(&t1, NULL);
-
-
+
+
if (argc < 8) {
- std::cout << "Usage: super_optimizer alpha_i alpha_m beta gamma delta target_innov [n_epochs n_steps n_reruns] [range] [verbose_level]\n";
- std::cout << " [T_start T_end Annealing_steps Target_rejection Starting_jump]\n";
- std::cout << " [FREEZE_alpha_i FREEZE_alpha_m FREEZE_beta FREEZE_gamma FREEZE_delta]\n";
-
+ 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 {
@@ -1663,25 +1709,25 @@
initSeed = atoi(argv[nArg]);
std::cout << par_names4[2] << ": " << initSeed << std::endl;
}
-
-
+
+
}
-
+
}
-
+
for (int j=0; j<5; j++){
-
+
cout << j << " | " << var_fixed[j] << " (fixed) \n";
}
-
+
target=params0[5];
range = (double)params1[3];
int identify_failed = 0;
char* filename= (char *)"movie_graph.txt";
int n_ep=params1[0], n_st=params1[1], n_rep=params1[2];
-
+
//...............................
-
+
for(int i=0; i<Nworkers; i++){
un[i] = new Universe((char *)filename,n_ep,n_st,
(int)n_rep,
@@ -1690,25 +1736,25 @@
CustomQueues[i] = dispatch_queue_create(i2char(i), NULL);
#endif
}
-
+
//...............................
if(n_rep > 0){
-
+
Results = new double[n_rep];
Counters = new int[n_rep];
-
+
}else {
-
+
Results = NULL;
Counters = NULL;
std::cout << " Number of reruns should be positive! " << std::endl;
return 0;
-
+
}
//...............................
//srand(time(0));
//srandomdev();
-
+
if ( initSeed != 0.0 ) {
srand(initSeed);
}
@@ -1721,84 +1767,84 @@
{
double r=0;
for (int j=0; j<100; j++){
-
-
-
+
+
+
r = rand()/(double)(pow(2.,31)-1.);
std::cout << r << " ";
}
std::cout << "\n";
}
//random initiation of starting parameters
-
+
if (range > 0.){
-
+
for (int i=0; i < 5; i++){
-
+
if (params0[i]==-100.){
-
+
double r1 = (rand()/(double)(pow(2.,31)-1.));
double r2 = (rand()/(double)(pow(2.,31)-1.));
double sign = 1.;
-
+
if(r1 > 0.5){
sign=-1.;
}
-
+
params0[i] = sign*r2*range;
-
+
std::cout << par_names0[i] << ": " << params0[i] << std::endl;
}
}
-
+
}
-
-
+
+
double T_start=params2[0], T_end=params2[1], Target_rejection=params2[3], starting_jump=params2[4];
int Annealing_repeats = (int) params2[2];
-
-#ifdef P_DISPATCH
+
+#ifdef P_DISPATCH
dispatch_group_t group = dispatch_group_create();
-
+
//.............................
- multi_annealing(group, un, CustomQueues, T_start, T_end, Target_rejection, Annealing_repeats,
+ multi_annealing(group, un, CustomQueues, T_start, T_end, Target_rejection, Annealing_repeats,
starting_jump, Results, Counters, params0, Annealing_repeats);
-
+
//dispatch_group_wait(group, DISPATCH_TIME_FOREVER);
dispatch_release(group);
//.............................
#else
- multi_annealing( un, T_start, T_end, Target_rejection, Annealing_repeats,
+ multi_annealing( un, T_start, T_end, Target_rejection, Annealing_repeats,
starting_jump, Results, Counters, params0, Annealing_repeats);
#endif
-
-
+
+
// stop timer
gettimeofday(&t2, NULL);
-
+
// compute and print the elapsed time in millisec
elapsedTime = (t2.tv_sec - t1.tv_sec) * 1000.0; // sec to ms
elapsedTime += (t2.tv_usec - t1.tv_usec) / 1000.0; // us to ms
elapsedTime /= 1000.;
cout << "\n*** optimizer completed, elapsed time=" << elapsedTime << " seconds " << elapsedTime/60. << " minutes)\n\n";
-
+
//.....................
-
+
for(int i=0; i<Nworkers; i++){
delete un[i];
}
-
+
//....................
if(n_rep > 0){
-
+
delete [] Results;
delete [] Counters;
-
+
}
-
+
return 0;
-
-
-
+
+
+
}
More information about the Swift-commit
mailing list