<commit_before>#!/usr/bin/python

import numpy
import numpy.random
import matplotlib
import matplotlib.pyplot as plt
import scipy.integrate
from scipy.integrate import quad
from scipy.stats import geom 
import checkdist
import optparse, sys
from optparse import OptionParser

parser = OptionParser()
parser.add_option("-t", "--temperature", nargs = 2, dest="T_k", type="float",default=[0.5,1.5],
                  help="low and high temperatures, [default = %default]") 
parser.add_option("-m", "--mu", nargs = 2, dest="mu_k", type="float",default=[0.6,1.4],
                  help="low and high chemical potentials, [default = %default]") 
parser.add_option("-c", "--cuttails", dest="cuttails", type="float",default=0.01,
                  help="fraction of the tails to omit from the analysis to avoid small sample errors in binning [default = %default]")
parser.add_option("-b", "--nboots", dest="nboots", type="int",default=200,
                  help="number of bootstrap samples performed [default = %default]")
parser.add_option("-r", "--nreps", dest="nreps", type="int",default=0,
                  help="number of independent repetitions of the sampling [default = %default]")
parser.add_option("-i", "--nbins", dest="nbins",type = "int", default=30,
                  help="number of bins for bootstrapping [default = %default]") 
parser.add_option("-e", "--energytype", dest="type", default="jointEN",
                  help="the type of energy that is being analyzed [default = %default]")
parser.add_option("-v", "--verbose", dest="verbose", action="store_true",default=False,
                  help="more verbosity")
parser.add_option("-N", "--number", nargs = 2, dest="N_k", type="int",default=[50000,50000],
                  help="number of samples from the two states, [default = %default]") 
parser.add_option("-K", "--force", dest="K", type="float",default=1.0,
                  help="spring force constant prefactor [default = %default]")
parser.add_option("-g", "--figureprefix", dest="figname", default='harmonic_mu',
                  help="name prefix for the figure")
parser.add_option("-s", "--seed", dest="seed", type = "int", default=None,
                  help="random seed for generating independet or bootstrap samples")

(options, args) = parser.parse_args()

# For N 1D harmonic oscillator, the potential is given by                                
#   V(x;K) = (K/2) * (x_-x_0)**2                                                                             
# where K denotes the spring constant.                                  
#                                                
# The equilibrium distribution is given analytically by                                      
#   p(x;beta,N) = sqrt[(beta K) / (2 pi)] exp[-beta K (x-x_0)**2 / 2]                                                       
# The dimensionless free energy is therefore                                                
#   q(beta,K) = sqrt(2 pi / beta K)
#   f(beta,K) = - (1/2) * ln[ (2 pi) / (beta K) ]              
#   f_2 - f_1 = -(1/2) (ln [ 1/ beta_2 K_2] - ln [ 1/beta_1 K_1])  
#   f_2 - f_1 = (1/2) ln [beta_2/K_2 / beta_1/K_1]  
#   f_2 - f_1 = (1/2) ln [K_1 T_1/ K_2 T_2]  
#   for T_1 = 0.9, T_2 = 1.1, K_1 = K_2, df = 0.5*ln(1.1/0.9) = 0.10034
#   for T_1 = 0.8, T_2 = 1.2, K_1 = K_2, df = 0.5*ln(1.2/0.8) = 0.2027
#   for T_1 = 0.5, T_2 = 1.5, K_1 = K_2, df = 0.5*ln(1.5/0.5) = 0.5493
#   Now add number of particles.  Use the model of identical, noninteracting particles.
#   In this case, Q(beta,N) = Q^N

#   Xi(x,N) = \sum_N Q(beta,N) exp(beta mu N)
#           = \sum_N (2 pi /beta K)^(N/2) exp(beta mu N)
#           = \sum_N sqrt(2 pi /beta K)^N exp(beta mu)^N
#           = \sum_N [sqrt(2 pi /beta K)^N exp(beta mu)]^N
#      for code simplicity, we assume N=1 is the minimum 
#           set x = (sqrt(2 pi / beta K)*exp(beta mu))
#           = ((1-x)^-1)-1 = 1/(1-x) - (1-x)/(1-x) = x/(1-x)
#
#   P(x,N) propto exp(-\beta (E - mu N))
#   P(x,N) exp(-\sum_i \betaK/2 x_i^2 -\beta mu N)/[1-(sqrt(2 pi/ beta K)exp(beta mu))]^(-1)
#   
# If there are two temperatures and the same chemical, we just bin the A_i = E_i - \mu N_i 
# So we want to sample from two temperatures and two mu.

# P_1 / P_2 = exp(B_2 PV_1 - B_2 PV_2 - (B_2 - B_1) E + (B_2 mu_2 - B_2 mu_2) 3N) 

#   analytical variance = depends on the method. 

title='harmonic oscillators with number'
if (options.nreps > 0):
    reptype = 'independent'
    nreps = options.nreps
if (options.nboots > 0):
    reptype = 'bootstrap'
    nreps = options.nboots
if (options.nboots > 0 and options.nreps > 0):
    print "Can't do both bootstrap sampling and independence sampling: defaulting to bootstrap sampling"

if (options.seed):
    numpy.random.seed(options.seed) # setting the seed for independent sampling 
    print "setting random number seed for generating samples as %d" % (options.seed)<commit_msg>Fixed some typos in the explanation.<commit_after>#!/usr/bin/python

import numpy
import numpy.random
import matplotlib
import matplotlib.pyplot as plt
import scipy.integrate
from scipy.integrate import quad
from scipy.stats import geom 
import checkdist
import optparse, sys
from optparse import OptionParser

parser = OptionParser()
parser.add_option("-t", "--temperature", nargs = 2, dest="T_k", type="float",default=[0.5,1.5],
                  help="low and high temperatures, [default = %default]") 
parser.add_option("-m", "--mu", nargs = 2, dest="mu_k", type="float",default=[0.6,1.4],
                  help="low and high chemical potentials, [default = %default]") 
parser.add_option("-c", "--cuttails", dest="cuttails", type="float",default=0.01,
                  help="fraction of the tails to omit from the analysis to avoid small sample errors in binning [default = %default]")
parser.add_option("-b", "--nboots", dest="nboots", type="int",default=200,
                  help="number of bootstrap samples performed [default = %default]")
parser.add_option("-r", "--nreps", dest="nreps", type="int",default=0,
                  help="number of independent repetitions of the sampling [default = %default]")
parser.add_option("-i", "--nbins", dest="nbins",type = "int", default=30,
                  help="number of bins for bootstrapping [default = %default]") 
parser.add_option("-e", "--energytype", dest="type", default="jointEN",
                  help="the type of energy that is being analyzed [default = %default]")
parser.add_option("-v", "--verbose", dest="verbose", action="store_true",default=False,
                  help="more verbosity")
parser.add_option("-N", "--number", nargs = 2, dest="N_k", type="int",default=[50000,50000],
                  help="number of samples from the two states, [default = %default]") 
parser.add_option("-K", "--force", dest="K", type="float",default=1.0,
                  help="spring force constant prefactor [default = %default]")
parser.add_option("-g", "--figureprefix", dest="figname", default='harmonic_mu',
                  help="name prefix for the figure")
parser.add_option("-s", "--seed", dest="seed", type = "int", default=None,
                  help="random seed for generating independet or bootstrap samples")

(options, args) = parser.parse_args()

# For N 1D harmonic oscillator, the potential is given by                                
#   V(x;K) = (K/2) * (x_-x_0)**2                                                                             
# where K denotes the spring constant.                                  
#                                                
# The equilibrium distribution is given analytically by                                      
#   p(x;beta,N) = sqrt[(beta K) / (2 pi)] exp[-beta K (x-x_0)**2 / 2]                                                       
# The dimensionless free energy is therefore                                                
#   q(beta,K) = sqrt(2 pi / beta K)
#   f(beta,K) = - (1/2) * ln[ (2 pi) / (beta K) ]              
#   f_2 - f_1 = -(1/2) (ln [ 1/ beta_2 K_2] - ln [ 1/beta_1 K_1])  
#   f_2 - f_1 = (1/2) ln [beta_2/K_2 / beta_1/K_1]  
#   f_2 - f_1 = (1/2) ln [K_1 T_1/ K_2 T_2]  
#   for T_1 = 0.9, T_2 = 1.1, K_1 = K_2, df = 0.5*ln(1.1/0.9) = 0.10034
#   for T_1 = 0.8, T_2 = 1.2, K_1 = K_2, df = 0.5*ln(1.2/0.8) = 0.2027
#   for T_1 = 0.5, T_2 = 1.5, K_1 = K_2, df = 0.5*ln(1.5/0.5) = 0.5493
#   Now add number of particles.  Use the model of distinguishable noninteracting particles.
#   In this case, Q(beta,N) = Q^N

#   Xi(x,N) = \sum_N Q(beta,N) exp(beta mu N)
#           = \sum_N (2 pi /beta K)^(N/2) exp(beta mu N)
#           = \sum_N sqrt(2 pi /beta K)^N exp(beta mu)^N
#           = \sum_N [sqrt(2 pi /beta K)^N exp(beta mu)]^N
#      for code simplicity, we assume N=1 is the minimum 
#           set x = (sqrt(2 pi / beta K)*exp(beta mu))
#           = ((1-x)^-1)-1 = 1/(1-x) - (1-x)/(1-x) = x/(1-x)
#
#   P(x,N) propto exp(-\beta (E - mu N))
#   P(x,N) exp(-\sum_i \betaK/2 x_i^2 -\beta mu N)/[1-(sqrt(2 pi/ beta K)exp(beta mu))]^(-1)
#   
# If there are two temperatures and the same chemical, we just bin the A_i = E_i - \mu N_i 
# So we want to sample from two temperatures and two mu.

# P_1 / P_2 = exp(B_2 PV_1 - B_2 PV_2 - (B_2 - B_1) E + (B_2 mu_2 - B_2 mu_2) N) 

#   analytical variance = depends on the method. 

title='harmonic oscillators with number'
if (options.nreps > 0):
    reptype = 'independent'
    nreps = options.nreps
if (options.nboots > 0):
    reptype = 'bootstrap'
    nreps = options.nboots
if (options.nboots > 0 and options.nreps > 0):
    print "Can't do both bootstrap sampling and independence sampling: defaulting to bootstrap sampling"

if (options.seed):
    numpy.random.seed(options.seed) # setting the seed for independent sampling 
    print "setting random number seed for generating samples as %d" % (options.seed)
{'fraction_non_alphanumeric': 0.15348600508905852, 'fraction_numerical': 0.033282442748091605, 'mean_word_length': 2.343314052398775, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1402645', 'n_tokens_mistral': 3457, 'n_tokens_neox': 3281, 'n_words': 964}