compute_weights.cpp

   1 /*
   2  * Copyright (c) 2009 Mauro Iazzi
   3  *
   4  * Permission is hereby granted, free of charge, to any person
   5  * obtaining a copy of this software and associated documentation
   6  * files (the "Software"), to deal in the Software without
   7  * restriction, including without limitation the rights to use,
   8  * copy, modify, merge, publish, distribute, sublicense, and/or sell
   9  * copies of the Software, and to permit persons to whom the
  10  * Software is furnished to do so, subject to the following
  11  * conditions:
  12  *
  13  * The above copyright notice and this permission notice shall be
  14  * included in all copies or substantial portions of the Software.
  15  *
  16  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  17  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
  18  * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  19  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  20  * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  21  * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  22  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  23  * OTHER DEALINGS IN THE SOFTWARE.
  24  */
  25
  26 #include <string>
  27 #include <vector>
  28 #include <iostream>
  29 #include <fstream>
  30 #include <sstream>
  31 #include <cmath>
  32 #include <cstdlib>
  33
  34 using namespace std;
  35
  36 struct Result {
  37     std::string prepend;
  38     std::vector<std::string> states;
  39     std::vector<double> energies;
  40     std::vector< std::vector<double> > weights;
  41     std::vector<int> lengths;
  42 };
  43
  44 void createarrays (std::string filename, Result& res, const int each, double c, bool use_log = true) {
  45     const int len = *(res.lengths.end());
  46     int n = 0;
  47     double *v = new double[len];
  48     int i = 0;
  49     for (int j=0;j<len;j++) v[i] = 0;
  50     std::string line_string = "";
  51     std::string state = "";
  52     std::ifstream f(filename.c_str(), std::fstream::in);
  53     while (!f.eof()) {
  54         double energy = 0.0;
  55         double weight = 0.0;
  56         getline(f, line_string);
  57         std::stringstream line(line_string);
  58         if (line_string[0]=='#') {
  59             res.prepend.append(line_string);
  60         } else if (line_string.size()!=0) {
  61             double nw = use_log?0.0:1.0;
  62             line >> state >> energy >> weight;
  63             res.states.push_back(state);
  64             res.energies.push_back(energy);
  65             if (use_log) {
  66                 weight = std::log(weight) - c;
  67             } else {
  68                 weight = weight * std::exp(-c);
  69             }
  70             v[i] = weight;
  71             int which_l = 0;
  72             for (int j=0;j<len;j++) {
  73                 if (use_log) {
  74                     nw += v[(len+i-j)%len];
  75                 } else {
  76                     nw *= v[(len+i-j)%len];
  77                 }
  78                 if (j==res.lengths[which_l]) {
  79                     res.weights[which_l].push_back(use_log?std::exp(nw):nw);
  80                     which_l++;
  81                 }
  82             }
  83             i = (i+1) % len;
  84         }
  85     }
  86     f.close();
  87     delete v;
  88 }
  89
  90 void analyze_data_file (std::string filename, std::ostream& out, const int len, double c, bool use_log = true) {
  91     int n = 0;
  92     double *v = new double[len];
  93     int i = 0;
  94     for (int j=0;j<len;j++) v[i] = 0;
  95     std::string line_string = "";
  96     std::string state = "";
  97     std::ifstream f(filename.c_str(), std::fstream::in);
  98     while (!f.eof()) {
  99         double energy = 0.0;
 100         double weight = 0.0;
 101         getline(f, line_string);
 102         std::stringstream line(line_string);
 103         if (line_string[0]=='#') {
 104             out << line_string << "\n";
 105         } else if (line_string.size()!=0) {
 106             double nw = use_log?0.0:1.0;
 107             line >> state >> energy >> weight;
 108             if (use_log) {
 109                 weight = std::log(weight) - c;
 110             } else {
 111                 weight = weight * std::exp(-c);
 112             }
 113             out << state << " " << energy;
 114             v[i] = weight;
 115             for (int j=0;j<len;j++) {
 116                 if (use_log) {
 117                     nw += v[(len+i-j)%len];
 118                 } else {
 119                     nw *= v[(len+i-j)%len];
 120                 }
 121             }
 122             out << " " << (use_log?std::exp(nw):nw);
 123             out << "\n";
 124             i = (i+1) % len;
 125         }
 126     }
 127     f.close();
 128     delete v;
 129 }
 130
 131 int main (int argc, char **argv) {
 132     string len_s = "1000";
 133     int len = 1000;
 134     double corr = 0.0;
 135     vector<string> files;
 136     for (int i=1;i<argc;i++) {
 137         if (argv[i][0]=='-') {
 138             switch (argv[i][1]) {
 139                 case 'c':
 140                     corr = atof(argv[++i]);
 141                     break;
 142                 case 'l':
 143                     len = atoi(argv[++i]);
 144                     len_s = argv[i];
 145                     break;
 146             }
 147         } else {
 148             files.push_back(argv[i]);
 149         }
 150     }
 151     for (vector<string>::iterator iter = files.begin();iter != files.end();iter++) {
 152         string ofn("l-");
 153         ofn.append(len_s);
 154         ofn.append("-");
 155         ofn.append(*iter);
 156         fstream o(ofn.c_str(), fstream::out);
 157         o.precision(10);
 158         o << "#!prepend " << len_s << "\n";
 159         analyze_data_file(*iter, o, len, corr);
 160         o.close();
 161     }
 162     return 0;
 163 }
 164