1D ising model - Helix-coil transistion

[Päällikkö]

Great! If you comment out everything dealing with neighbors in the localEnergy() function (the loop), and multiply the variable ener by .5f, you should get a curve that goes on top of the green one.

Finally, there are a couple of issues in the code you should be aware of. Most critically the way the random generator is instantiated: Each time a particle is created, it constructs a new (Mersenne-twister) random generator, std::mt19937, with default arguments. This means that the seed is always the same, which is a bad thing (each temperature is going to get exactly the same random numbers). To avoid this, one could make the variable static: "static std::mt19937 mGen;". This requires an instantiation, "std::mt19937 Particle::mGen = std::mt19937{};" somewhere in the code (usually in the .cpp file where you define the class, i.e. does not have to be inside a function). What static does is "share" the variable among all variables that are of type Particle (so mGen will be destructed only at program exit rather than when each Particle is destructed).

There might also be a correlation between the random number generator inside the particle and the random number generator that is used to sample the Monte Carlo energy condition. One might avoid this by seeding one of the generators from system time or something else (like a Singleton class for a global random generator etc).

These are issues you don't need to worry about for simple coursework, but if you start doing serious research and/or start noticing curious non-random patterns, these might turn out to be the culprits.

This might have been a bit difficult to follow especially if you're not very comfortable with C++.

 

[Matt atkinson]

Once again, thank you so much you really helped me, helping me with my program and pointing out some of the pitfall's with it! :D
Now to calculate errors!

 

[Matt atkinson]

Just a question I am trying to work out errors with your code using a statistics class,

class statistics{
      private:
              int n;
              double sum;
              double sumsq;
      public:
             statistics();
             int getNumber() const;
             double getAverage() const;
             double getSqAverage() const;
             void add(double x);
};

int statistics::getNumber()             const{
                                              return n;
                                              }
double statistics::getAverage()         const {
                                              if(n==0) return 0.0;
                                              return sum/n;
                                              }
double statistics::getSqAverage()       const{
                                              if(n==0) return -1;
                                              return sumsq/n;
                                              }
void statistics::add(double x)                {
                                               n++;
                                               sum += x;
                                               sumsq += x*x;
       }

But I am not sure where to implement it in the main i tried doing something like;

res.add();
error=sqrt((res.getSqAverage()-res.getAverage()*res.getAverage())/(res.getNumber()));

Im pretty sure that res isn't the thing i should be using I am just strugglign to figure out what i should.

 

[Päällikkö]

Well, you're not really specifying what the errors might be, nor have you supplied a minimal example (i.e. please do not include Particle in the example: These objects are supposed to work independently) of using the Statistics class (I prefer user defined type names to be written with the first letter capitalized).

You did not write out the constructor: I hope you are initializing n, sum, and sumsq to zero. I also hope that res.add() is not an actual use case, but rather you are adding something (you did declare Statistics res; and didn't just try to use the std::vector<float> res; I had used, right? This part of your text was a bit ambiguous.).

If you are using the STL's sqrt from <cmath>, it might be better to be explicit about it: std::sqrt. I don't think the standard mandates the STL implementation to declare the function in the global namespace (but allows it to), whereas it must be present in std:: (don't quote me on this, though).

 

[Matt atkinson]

Oh okay, and my interpretation of the errors was that each value was an average of many measurements so each of the points i plot would have a spread of values that could be represented as an error bar?

Also no I didn't declare statistics res, I realize that was a mistake now.

Also sorry i forgot to paste the constructor,

statistics::statistics()                     {
                                              n=0;
                                              sum=0.0;
                                              sumsq=0.0;
                                              }

 

[Päällikkö]

So medres in my code collects several data points at each temperature, say 10000. res, on the other hand, stores an average of these 10000 values to represent the result at each particular temperature. What you want is not only the mean of medres (which is res), but also the standard deviation, sigma (which is an error bar estimate, although quite often people use 2*sigma as their error estimate, or at LHC/CERN, 5*sigma when they announced they had found the Higgs boson). From what I understand, you are trying to achieve a universal way of going about this with your statistics class, which is a good idea and ought to clarify the code.

If you paste your main, or better yet, a minimal example of how you want to use your statistics class (that is, one without interaction with the Particle class etc.), things would probably be easier to debug (if there sitll are bug). This might also clarify to yourself what the limits of your clss design might be (for example, it is not immediately clear to me how you suppose to use the class for multiple temperatures, maybe it is just a dropin for medres, rather than res?).

Finally, I'd recommend using

statistics::statistics() : n(0), sum(0.0), sumsq(0.0) {}

Here it doesn't really matter, but suppose n were const int rather than an int. Only the command I listed would then work (you can initialize with a value, but cannot assign one; same is true if you tried const int x; x = 1; versus const int x = 1; The latter is cheaper too when not talking about consts).

 

[Matt atkinson]

well for example i want to calculate $\sigma$ for each mean value of medres, for each temperature. So inside the loop for temperature, i could put something like this ;

A=medres.getAverage();
B=medres.getSqAverage();
Cnum=medres.getNumber();
Sigma=std::sqrt((B-A*A)/(Cnum));

Not sure if I am understanding it right.

 

[Päällikkö]

That'd be correct, yes, if you define "statistics medres;" and rather than "medres.push_back(particle.mean());" you use "medres.add(particle.mean());"
Sigma should be an array like res (res stores the means, Sigma should store the sigmas at each temperature), or rather res could be replaced by A should A also be made into an array. Given that you compute A, B, C using the methods of statistics, you'd do well to make Sigma a method of statistics, too (say, getError(), or whatever you consider appropriate).

 

[Matt atkinson]

Oh okay thankyou! ill try that when i get back to my code.
Define sigma, in the statistics class then have something like;

sigma=std::sqrt((sumsq-sum*sum)/(n))

i believe this should work
also with

class statistics{
      private:
              int n;
              double sum;
              double sumsq;
              double sigma;
              double medres;
      public:
             statistics();
             int getNumber() const;
             double getAverage() const;
             double getSqAverage() const;
             void add(double x);
};

If i have understood correctly

 

[Matt atkinson]

thankyou for all your help, just checked and all the errors i get are the exact same (~0.00042) have i done something incorrectly?
I think its my definition of medres2.

int main() {
  std::uniform_real_distribution<float> dist{0.f, 1.f};
  std::mt19937                          gen;

  std::vector<float> Ts;
  std::ofstream ofile;
  std::array<float, N>   A,B,C;
  statistics medres2;
  float error;
  ofile.open ("helixfrac.txt");
  for(float T = 420.f; T < 434.f; T+=.025f)
    Ts.push_back(T);
  std::vector<float> res;
  for(auto T : Ts) {
    Particle particle(T);
    std::vector<float> medres;
    medres.reserve(10000);
    for(int i = 0; i < 100000; ++i) {
      particle.choose();
      auto Eb = particle.localEnergy();
      particle.perturb();
      auto Ef = particle.localEnergy();
      if(std::exp(-(Ef-Eb)/T) < dist(gen))
        particle.perturb();
      if(!(i%10) && i > 1000){
        medres.push_back(particle.mean());
        medres2.add(particle.mean());
                            }
    }
    float xx = 0.f;
    for(const auto x : medres)
      xx += x;
    res.push_back(xx/medres.size());
  }
  for(unsigned int i = 0; i < Ts.size(); ++i){
    A[i]=medres2.getSqAverage();
    B[i]=medres2.getAverage();
    C[i]=1/(double) medres2.getNumber();
    std::cout << Ts[i] << "," << (1+res[i])*.5 << std::endl;
    ofile << Ts[i] << "," << (1+res[i])*.5 << "," << std::sqrt((A[i]*C[i]-B[i]*B[i]*C[i])) << std::endl;
    }
    ofile.close();
  std::cout << std::endl;
}

 

[Päällikkö]

It sounds more or less correct. There are caveats in the way that you compute the error (it gives an underestimate), but it is not really your fault, it's a simplification made in your course materials. The reason is that the measurements are highly correlated with each other (the length of the polymer is, perhaps, N = 512, while the !(i%10) means that you are storing an average after every 10 attempts at flipping: Obviously even if all the flips were successful, you'd only move the average computer by Particle::mean() a few percentages at best). If you are interested in better methods, look up block-averaging and other autocorrelation based calculations (the idea is that you wait for as long as there is no correlation in the mean data to take averages). Vary N and the number of flips per temperature (i < 10000) and see how they affect the error estimates.

Finally, there's no reason for you to keep the medres's in your code. I cleaned it up a bit (still not pretty, but now it doesn't have anything that's not needed). I also use a couple of different loop styles, so you can compare and see what you feel comfortable with. For example, stats is an array (std::vector) of Statistics objects, so that each temperature gets its own Statistics. Now I initialize stat as the iterator (if you're not familiar with iterators, just think of it as a pointer) to the first Statistics object in the array, and then increment the iterator at each temperature. When I print out, I use a more traditional loop, where everything is indexed (note, though, that many data types, such as std::list, do not permit direct indexing, i.e. random access, but rather you must traverse through them).

int main() {
  auto dist = std::uniform_real_distribution<float>{0.f, 1.f};
  auto gen  = std::mt19937{};

  auto Ts = std::vector<float>{};
  for(float T = 426.f; T < 428.f; T+=.05f)
    Ts.push_back(T);

  std::vector<Statistics> stats(Ts.size());

  auto stat = stats.begin();
  for(auto T : Ts) {
    auto particle = Particle(T);
    for(int i = 0; i < 10000; ++i) {
      particle.choose();
      auto Eb = particle.localEnergy(); 
      particle.perturb();
      auto Ef = particle.localEnergy();
      if(std::exp(-(Ef-Eb)/T) < dist(gen))
        particle.perturb();
      if(!(i%100) && i > 5000)
        stat->add(particle.mean());
    }
    ++stat;
  }
  for(unsigned int i = 0; i < Ts.size(); ++i) {
    auto mean  = stats[i].getAverage();
    auto stdev = stats[i].getStdDeviation();
    std::cout << Ts[i] << " " 
              << (1+mean)*.5 << " " 
              << (1+mean-2*stdev)*.5 << " " 
              << (1+mean+2*stdev)*.5 << std::endl;
  }
  std::cout << std::endl;
}

You'll want to add the file outstream (i.e. std:fstream) instead of std::cout to write directly to a file. The reason I'm using 2*stdev is that I consider this a better estimate of the error than stdev (and this is indeed more commonly used). You might want to change this.

Finally, as you may notice, I made some modifications to your statistics class:

class Statistics {
public:
  int    getNumber() const;
  float  getAverage() const;
  float  getStdDeviation() const;
  void   add(double x);
private:
  int    n     = 0;
  double sum   = 0.f;
  double sumsq = 0.f;
};

int Statistics::getNumber() const {
  return n;
}

float Statistics::getAverage() const {
  if(n==0) 
    return 0.0;
  return sum/n;
}

float Statistics::getStdDeviation() const {
  if(n==0) 
    return 0.0;
  return std::sqrt(sumsq - sum*sum/n)/n;
}

void Statistics::add(double x) {
   n++;
   sum += x;
   sumsq += x*x;
}

[/code]

 

[Matt atkinson]

Wow thank you, your help has been greatly appreciated with this. I learned quite a lot about coding and physics too, while also doing the problem.
A lot of the code you used, i didn't understand at first but your explanations were great and anything else I looked up.

Thanks again!