staci_split.cpp.backup

#include <stdio.h>
#include <vector>
#include "Staci.h"
#include "xmlParser.h"
#include <Eigen/Dense>

#include <ga/ga.h>
#include <ga/GASimpleGA.h>
#include <ga/GAListGenome.h>

#include <igraph/igraph.h>

float Objective(GAGenome &);

float Objective2(GAGenome &);

using namespace std;
using namespace Eigen;

#define cout STD_COUT

// These global variables are coming from the datafile
string fname;
int popsize, ngen;
double pcross, pmut;
string weight_type;

Staci *wds;
int n_n, n_p, n_comm;
int fcount;
double f_ce;
vector<int> pipe_idx;
VectorXd p;
double sumW;

int n_c;
double Qmod;

bool info;

void LoadMatrices(VectorXd &p, string weight_type);

void LoadSystem(string fname);

void save_state();

void logfile_write(string msg, int debug_level);

int MyFlipMutator(GAGenome & c, float prob_mut);
int QMutator(GAGenome & c, float prob_mut);

string Load_Settings();

int global_debug_level;

vector<int> best;
double best_Q;
double obj_offset;
double nc_weight;

stringstream strstrm;
string logfilename;

void prepare_weights_vector(igraph_vector_t* weights, const igraph_t* graph) {
  int i, n = igraph_ecount(graph);
  igraph_vector_resize(weights, n);
  for (i = 0; i < n; i++) {
    VECTOR(*weights)[i] = i % 5;
  }
}

int main(int argc, char **argv) {

    best_Q = -100.;
    obj_offset = 1.0;
    info = false;
    fcount = 0;
    nc_weight = 0.5;

    string msg = Load_Settings();
    
    // Clear logfile
    ofstream ofs(logfilename.c_str(), std::ios::out | std::ios::trunc);
    ofs.close();

    logfile_write(msg,0);

    // Load system
    LoadSystem(fname);
    f_ce = pow((double) n_comm, (double) n_n);

    // Load matrices
    LoadMatrices(p, weight_type);

    // Optimize
    cout << endl << "RUNNING OPTMIZATION...";
    cout << flush;

    int *aset = new int[n_comm];
    for (int i = 0; i < n_comm; i++)
        aset[i] = i;
    GAAlleleSet<int> alleles(n_comm, aset);

    GA1DArrayAlleleGenome<int> genome(n_n, alleles, Objective);
    genome.mutator(QMutator);

    GASteadyStateGA ga(genome);

    ga.set(gaNpopulationSize, popsize);    // population size
    ga.set(gaNpCrossover, pcross);        // probability of crossover
    ga.set(gaNpMutation, pmut);        // probability of mutation
    ga.set(gaNnGenerations, ngen);        // number of generations
    ga.set(gaNscoreFrequency, 100);        // how often to record scores
    ga.set(gaNflushFrequency, 100);    // how often to dump scores to file
    ga.set(gaNselectScores,        // which scores should we track?
       GAStatistics::Maximum | GAStatistics::Minimum | GAStatistics::Mean);
    ga.set(gaNscoreFilename, "bog.dat");

    ga.evolve();
    cout << endl << "done." << endl;

    genome = ga.statistics().bestIndividual();
    //double best_obj = Objective(genome);

    strstrm.str("");
    strstrm << endl << "BEST SOLUTION FOUND:" << endl;

    for (unsigned int i = 0; i < n_n; i++)
        strstrm << genome.gene(i) << " ";
    strstrm << endl;

    info = true;
    double best_obj = Objective(genome);

    strstrm << endl << "Some more info:";
    strstrm << endl << "\t datafile   : " << fname;
    strstrm << endl << "\t # of nodes : " << n_n;
    strstrm << endl << "\t # of pipes : " << n_p;
    strstrm << endl << "\t D          : " << (wds->GetMinPipeDiameter()) << " ... " << (wds->GetMaxPipeDiameter())
    << " m";
    strstrm << endl << "\t L          : " << (wds->GetMinPipeLength()) << " ... " << (wds->GetMaxPipeLength()) << " m";
    strstrm << endl << "\t sum. L     : " << (wds->GetSumPipeLength()) << " m ";
    strstrm << endl << "\t max. cons. : " << (wds->GetMaxConsumption()) << " kg/s = "
    << 3.6 * (wds->GetMaxConsumption())
    << " m3/h" << endl;

    logfile_write(strstrm.str(), 0);

    return 0;
}

int 
MyFlipMutator(GAGenome & c, float prob_mut)
{
    GA1DArrayAlleleGenome<int> &child = (GA1DArrayAlleleGenome<int> &) c;

    register int n, i;
    if(prob_mut <= 0.0) return(0);

    float nMut = prob_mut * (float)(child.length());
  if(nMut < 1.0){       // we have to do a flip test on each bit
    nMut = 0;
    for(i=child.length()-1; i>=0; i--){
      if(GAFlipCoin(prob_mut)){
        child.gene(i, child.alleleset().allele());
        nMut++;
    }
}
}
  else{             // only flip the number of bits we need to flip
    for(n=0; n<nMut; n++){
      i = GARandomInt(0, child.length()-1); // the index of the bit to flip
      child.gene(i, child.alleleset().allele());
  }
}
return((int)nMut);
}

int QMutator(GAGenome & c, float prob_mut){
    GA1DArrayAlleleGenome<int> &child = (GA1DArrayAlleleGenome<int> &) c;
    
    // Copy of child
    vector<int> copy_of_child(child.size());
    for (unsigned int i=0; i<child.size(); i++)
        copy_of_child.at(i)=child.gene(i);

    bool mut_info = false;
    info = false;
    int mut_round=0;
    bool next_round=true;
    double obj_before=0., obj_after=1.;
    int nMut = 0;

    if(prob_mut <= 0.0) return(0);

    if (mut_info){
        cout<<endl<<endl<<"MUTATION:"<<endl<<endl<<"Initial gene:"<<flush;
        // cin.get();
    }
    
    obj_before = Objective(child);
    while (next_round && (mut_round<10)){
        if (mut_info)
            cout<<endl<<endl<<"ROUND "<<mut_round;

        int n_c = 0;
        vector<bool> cut_idx(n_n,false);
        vector<int> cut_pair(n_n,0);
        for (unsigned int i = 0; i < n_p; i++) {
            int idx_n1 = wds->agelemek.at(pipe_idx.at(i))->Get_Cspe_Index();
            int idx_n2 = wds->agelemek.at(pipe_idx.at(i))->Get_Cspv_Index();

            if (child.gene(idx_n1) != child.gene(idx_n2)){
                n_c++;
                cut_idx.at(idx_n1)=true;
                cut_idx.at(idx_n2)=true;
                cut_pair.at(idx_n1)=idx_n2;
                cut_pair.at(idx_n2)=idx_n1;
                if (mut_info){
                    cout<<endl<<"\tCut: "<<idx_n1 <<"("<<child.gene(idx_n1);
                    cout<<") --x-- "<<idx_n2<<"("<<child.gene(idx_n2)<<")";
                }
            }
        }
        if (mut_info)
            cout<<endl;

        for(unsigned int i=0; i<n_n; i++){
            if (cut_idx.at(i)){
                int other_node_community = child.gene(cut_pair.at(i));
                if (mut_info){
                    cout<<endl<<" --> idx="<<i<<", pair: "<<cut_pair.at(i)<<" with comm. "<<other_node_community;
                }


                cut_idx.at(i)=false;
                cut_idx.at(cut_pair.at(i))=false;

                child.gene(i,other_node_community);
                nMut++;
            }
        }
        if (mut_info)
            cout<<endl;

        obj_after = Objective(child);
        
        if (obj_after>obj_before){
            next_round=true;
            mut_round++;
            obj_before=obj_after;
            for (unsigned int i=0; i<child.size(); i++)
                copy_of_child.at(i)=child.gene(i);
        }
        else{
            next_round=false;
            for (unsigned int i=0; i<child.size(); i++)
                child.gene(i,copy_of_child.at(i));
        }
    }
    info = false;

    // if no improvement was made, switch to "normal" mutation
    if (mut_round==0){
        nMut = 0;
        for (int i=0; i<child.length(); i++){
            if (GAFlipCoin(prob_mut)){
                child.gene(i, child.alleleset().allele());
                nMut++;
            }
        }
        if (mut_info)
            cout<<endl<<endl<<"Done with FlipCoin mutation, nMut="<<nMut<<endl;
    }

    if (mut_info){
        cout<<endl<<"Leaving QMutator()..."<<endl;
        // cin.get();
    }
    return((int)nMut);
}


void save_state() {
    for (int i = 0; i < n_n; i++){
        wds->cspok.at(i)->Set_dprop("concentration", best.at(i));
    }

    for (int i = 0; i < wds->agelemek.size(); i++) {
        int idx_n1 = wds->agelemek.at(i)->Get_Cspe_Index();
        int idx_n2 = wds->agelemek.at(i)->Get_Cspv_Index();
        if (idx_n1<0)
            idx_n1=idx_n2;
        if (idx_n2<0)
            idx_n2=idx_n1;
        double comm1 = (double) best.at(idx_n1);
        double comm2 = (double) best.at(idx_n2);
        double comm = (comm1+comm2)/2.;
        // cout<<endl<<"edge #"<<i<<": "<<wds->agelemek.at(i)->Get_nev()<<", community # "<<comm;
        wds->agelemek.at(i)->Set_dprop("concentration", comm);
    }

    wds->save_mod_prop_all_elements("concentration");
}

float
Objective(GAGenome &c) {
    GA1DArrayAlleleGenome<int> &genome = (GA1DArrayAlleleGenome<int> &) c;

    fcount++;

    // Extract genome to more convenient form
    vector<int> tmp(n_n);
    //tmp.at(0) = 0;
    for (int i = 0; i < n_n; i++)
        tmp.at(i) = genome.gene(i);

    // Compute number of cuts
    n_c = 0;
    for (unsigned int i = 0; i < n_p; i++) {
        int idx_n1 = wds->agelemek.at(pipe_idx.at(i))->Get_Cspe_Index();
        int idx_n2 = wds->agelemek.at(pipe_idx.at(i))->Get_Cspv_Index();

        if (genome.gene(idx_n1) != genome.gene(idx_n2))
            n_c++;
    }

        // Modularity with weight W, p = abs(Apn.transpose()) * w
    Qmod = 0.;
    double Qtmp;
    for (int m = 0; m < n_comm; m++) {
        Qtmp = 0.;
        for (int i = 0; i < n_n; i++) {
            if (genome.gene(i) == m) {
                Qtmp += p(i);
            }
        }
        Qmod += (Qtmp / 2 / sumW) * (Qtmp / 2 / sumW);
    }

    double Q = 1. - nc_weight*((double) n_c / (double) n_p) - (1.-nc_weight) *Qmod;

    if ((Q > best_Q) || (info)) {
        best = tmp;
        best_Q = Q;

        save_state();

        strstrm.str("");
        if (fcount < 1000)
            strstrm << endl << "fcount : " << fcount << " (" << ((double) fcount) / f_ce * 100.
        << " % of compl. enum.)";
        if (fcount < 1e6)
            strstrm << endl << "fcount : " << ((double) fcount) / 1000. << "k (" << ((double) fcount) / f_ce * 100.
        << " % of compl. enum.)";
        else
            strstrm << endl << "fcount : " << ((double) fcount) / 1000000. << "M (" << ((double) fcount) / f_ce * 100.
        << " % of compl. enum.)";
        strstrm << endl << "n_c/np : " << ((double) n_c / (double) n_p) << " (=n_c/np with n_c = " << n_c << ", np = "
        << n_p
        << ")";
        strstrm << endl << "Qmod   : " << Qmod;
        strstrm << endl << "Q      : " << Q << endl << "----------------------------------------------------";
        logfile_write(strstrm.str(), 0);

    }

    return (obj_offset + Q);

}

void LoadMatrices(VectorXd &p, string weight_type) {

    MatrixXd A = MatrixXd::Zero(n_n, n_n);
    MatrixXd Apn = MatrixXd::Zero(n_p, n_n);
    MatrixXd abs_Apn = MatrixXd::Zero(n_p, n_n);
    VectorXd k = VectorXd::Zero(n_n);
    p = VectorXd::Zero(n_n);
    VectorXd W = VectorXd::Zero(n_p);

    // Topological incidence matrix
    for (unsigned int j = 0; j < n_p; j++) {
        int idx_n1 = wds->agelemek.at(pipe_idx.at(j))->Get_Cspe_Index();
        int idx_n2 = wds->agelemek.at(pipe_idx.at(j))->Get_Cspv_Index();
        Apn(j, idx_n1)++;
        Apn(j, idx_n2)--;
        abs_Apn(j, idx_n1)++;
        abs_Apn(j, idx_n2)++;
    }

    A = Apn.transpose() * Apn;
    for (unsigned i = 0; i < n_n; i++) {
        A(i, i) = 0;
        for (unsigned int j = 0; j < n_n; j++) {
            if (A(i, j) != 0) {
                A(i, j) = 1.;
                k(i)++;
            }
        }
    }

    if (0 == strcmp(weight_type.c_str(), "topology")) {
        for (unsigned int i = 0; i < W.size(); i++)
            W(i) = 1;
    }

    if (0 == strcmp(weight_type.c_str(), "dp")) {
        double dp;
        for (unsigned int i = 0; i < W.size(); i++){
            dp = wds->agelemek.at(pipe_idx.at(i))->Get_dprop("headloss");
            W(i) = abs(dp);
            cout<<endl<<wds->agelemek.at(pipe_idx.at(i))->Get_nev()<<" dp="<<dp;
        }
            // cin.get();
    }

    // Final computations
    p = abs_Apn.transpose() * W;
    sumW = W.sum();

}

void LoadSystem(string fname) {
    wds = new Staci(fname.c_str());
    wds->build_system();
    wds->ini();
    wds->solve_system();
    wds->set_res_file(wds->get_def_file());
    wds->save_results(true);

    n_n = wds->cspok.size();

    for (unsigned int i = 0; i < wds->agelemek.size(); i++) {
        string tipus = wds->agelemek.at(i)->Get_Tipus();
        if (0 == strcmp(tipus.c_str(), "Cso"))
            pipe_idx.push_back(i);
    }
    n_p = pipe_idx.size();
}

void logfile_write(string msg, int debug_level) {
    if (global_debug_level >= debug_level) {
        ofstream logfile;
        logfile.open(logfilename.c_str(), ios::app);
        logfile << msg;
        logfile.close();
        cout << msg;
    }
}

string Load_Settings() {

    XMLNode xMainNode = XMLNode::openFileHelper("staci_split_settings.xml", "settings");
    global_debug_level = atoi(xMainNode.getChildNode("global_debug_level").getText());
    n_comm = atoi(xMainNode.getChildNode("n_comm").getText());
    weight_type = xMainNode.getChildNode("weight_type").getText();
    fname = xMainNode.getChildNode("fname").getText();
    logfilename = xMainNode.getChildNode("logfilename").getText();

    popsize = atoi(xMainNode.getChildNode("popsize").getText());
    ngen = atoi(xMainNode.getChildNode("ngen").getText());
    pmut = atof(xMainNode.getChildNode("pmut").getText());
    pcross = atof(xMainNode.getChildNode("pcross").getText());

    stringstream msg;
    msg.str("");
    msg << "Settings:" << endl << endl;
    msg << "\t global_debug_level     : " << global_debug_level << endl;
    msg << "\t n_comm                 : " << n_comm << endl;
    msg << "\t weight_type            : " << weight_type;
    if ((weight_type=="topology") || (weight_type=="dp")){
     msg<< " (ok)"<<endl;
 }
 else{
    cout<<endl<<"!!!! UNKNOWN WEIGHT_TYPE!!!"<<endl<<"Exiting..."<<endl;
    exit(-1);
}
msg << "\t fname                  : " << fname << endl;
msg << "\t logfilename            : " << logfilename << endl<<endl;
msg << "\t popsize : " << popsize << endl;
msg << "\t ngen    : " << ngen << endl;
msg << "\t pmut    : " << pmut << endl;
msg << "\t pcross  : " << pcross << endl << endl;

return msg.str();
}