SSABase.hpp 1.8 KB
Newer Older
wadaniel's avatar
wadaniel committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73
#pragma once

#include <stdio.h>
#include <stdlib.h>
#include <cmath>
#include <assert.h>
#include <algorithm>
#include <vector>
#include <fstream>
#include <iomanip>
#include <iostream>

class SSABase
{
protected:
  /** Multiplier for number of molecules. */
  const int omega;
  /** Number of samples of stochastic algo to run. */
  const int numSamples;
  /** End time for algorithm. */
  const double tend;
  /** Width of trajectory bin. */
  const double bin_dt;

  /** Trajectories (binned) for S1*omega and S2*omega (trajS(i) = avg of S for t in [i,i+1]*bin_dt). */
  std::vector<double> trajS1;
  std::vector<double> trajS2;
  std::vector<int> ntraj;

  /** Store timing result here. */
  double time;

public:
  SSABase(int _omega, int num_samples, double T, double bin_dt)
    : omega(_omega)
    , numSamples(num_samples)
    , tend(T)
    , bin_dt(bin_dt)
    , trajS1(static_cast<unsigned int>(std::ceil(T/bin_dt)))
    , trajS2(static_cast<unsigned int>(std::ceil(T/bin_dt)))
    , ntraj(static_cast<unsigned int>(std::ceil(T/bin_dt)))
    , time(0.0)
  { }

  /** Return time taken by start/stopTiming. */
  double getTime() {
    return time;
  }
  /** Get omega. */
  double getOmega() {
    return omega;
  }

  /**
   * Dump trajectories to file (nbins rows and 2 columns).
   */
  void dumpTrajToFile(std::string filename) {
    // Open file stream
    std::ofstream outfile(filename.c_str());

    // Dump it
    const int nbins = trajS1.size();
    int totalevals = 0;
    for (int i = 0; i < nbins; ++i) {
      // must rescale
      outfile << i*bin_dt+bin_dt/2 << " "
              << (trajS1[i] / omega) << " "
              << (trajS2[i] / omega) << std::endl;
      totalevals += ntraj[i];
    }
    std::cout << "Average number of time steps per sample: " << double(totalevals)/numSamples << std::endl;
  }
};