based on demo_JPAC_fit.ipynb
import PyPWA as pwa
import pandas
import numpy as npy
from IPython.display import display
import warnings
warnings.filterwarnings('ignore')
import matplotlib.pyplot as plt
#import ROOT
import pyroot
import emcee
import corner
from scipy import stats
import os
pre-processed by process_data.py
outputpath = "example/"
plotpath = outputpath+"/plots_mcmc/"
data = pwa.read("../TUTORIAL/simdata_JPAC-np.csv")
mc = pwa.read("../TUTORIAL_FILES/etapi_acc.csv")
initial = []
Vs = dict()
with open(f'{outputpath}/wave_set.csv', 'r') as file:
for line in file:
tokens = line.split(',')
if line.find('errordef') and line.find('limit') and line.find('error') and line.find('fix'):
initial.append(tokens[0])
Vs.update({tokens[0]: tokens[1].strip("\n")})
print(initial)
print(Vs)
['r.1.0.0', 'r.1.1.0', 'r.1.1.1', 'r.1.2.0', 'r.1.2.1', 'r.1.2.2', 'i.1.0.0', 'i.1.1.0', 'i.1.1.1', 'i.1.2.0', 'i.1.2.1', 'i.1.2.2'] {'errordef': '1', 'r.1.0.0': '20', 'limit_r.1.0.0': '[-500', 'error_r.1.0.0': '0.1', 'r.1.1.0': '10', 'limit_r.1.1.0': '[-500', 'error_r.1.1.0': '0.1', 'r.1.1.1': '10', 'limit_r.1.1.1': '[-500', 'error_r.1.1.1': '0.1', 'r.1.2.0': '10', 'limit_r.1.2.0': '[-500', 'error_r.1.2.0': '0.1', 'r.1.2.1': '10', 'limit_r.1.2.1': '[-500', 'error_r.1.2.1': '0.1', 'r.1.2.2': '10', 'limit_r.1.2.2': '[-500', 'error_r.1.2.2': '0.1', 'i.1.0.0': '0', 'limit_i.1.0.0': '[-500', 'error_i.1.0.0': '0.1', 'i.1.1.0': '10', 'limit_i.1.1.0': '[-500', 'error_i.1.1.0': '0.1', 'i.1.1.1': '10', 'limit_i.1.1.1': '[-500', 'error_i.1.1.1': '0.1', 'i.1.2.0': '10', 'limit_i.1.2.0': '[-500', 'error_i.1.2.0': '0.1', 'i.1.2.1': '10', 'limit_i.1.2.1': '[-500', 'error_i.1.2.1': '0.1', 'i.1.2.2': '10', 'limit_i.1.2.2': '[-500', 'error_i.1.2.2': '0.1', 'fix_i.1.0.0': 'True'}
To use as start parameters for MCMC
minuitresults = pandas.read_csv(f"{outputpath}/minuit.csv")
minuitresults
r.1.0.0 | r.1.1.0 | r.1.1.1 | r.1.2.0 | r.1.2.1 | r.1.2.2 | i.1.0.0 | i.1.1.0 | i.1.1.1 | i.1.2.0 | i.1.2.1 | i.1.2.2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 154.600390 | -0.087148 | 0.728922 | 75.851323 | 74.024774 | 70.926984 | 0.0 | 13.230514 | -0.835655 | 31.386883 | -21.400500 | -19.979625 |
1 | 194.342544 | -1.138565 | -1.479713 | 82.987672 | 77.896378 | 84.728050 | 0.0 | 35.397619 | 19.381017 | 1.177333 | 33.451012 | 4.224923 |
2 | 281.450602 | 1.672613 | 0.627708 | 89.276097 | 85.600139 | 89.558392 | 0.0 | -9.862055 | 11.428773 | 43.070308 | -5.905347 | 5.226695 |
3 | 500.845930 | -0.551314 | -7.070445 | 86.068007 | 82.920974 | 91.906067 | 0.0 | 32.661097 | -12.264725 | 33.943332 | 92.798853 | 70.652099 |
4 | 423.063231 | 472.461037 | -109.284008 | -220.114224 | 109.897267 | 13.009496 | 0.0 | 25.224835 | -401.768350 | 218.521003 | 328.147971 | 64.844637 |
5 | 506.378309 | 2.658741 | -1.158987 | -82.175805 | -82.963506 | -83.671884 | 0.0 | -2.891133 | 13.541788 | -85.975670 | -97.063776 | -80.395342 |
6 | 238.013167 | -3.645763 | -2.833371 | -151.660664 | -151.462880 | -140.581859 | 0.0 | 91.905955 | 70.987310 | -0.444448 | -2.569730 | 4.098666 |
7 | 165.091849 | 6.079976 | 5.965817 | -164.410058 | -168.327166 | -162.503420 | 0.0 | 2.468134 | 0.568305 | 108.579034 | 103.145066 | 121.461943 |
8 | 87.954769 | -11.963408 | -39.614286 | -276.333138 | -251.956378 | -277.150383 | 0.0 | 34.292389 | 62.774841 | -79.652286 | -146.368104 | -50.158416 |
9 | -12.548544 | -17.245086 | -12.611267 | 340.928511 | 340.258386 | 341.078252 | 0.0 | 101.395859 | 79.270563 | 42.019583 | 49.062637 | 38.301091 |
10 | -42.958850 | -29.775779 | -1.092006 | -160.050995 | -172.109519 | -178.292795 | 0.0 | -54.245874 | -11.594841 | 120.000661 | 105.580083 | 98.206337 |
11 | -17.028557 | -12.610656 | -12.411596 | -91.120698 | -91.671758 | -83.495385 | 0.0 | 65.026192 | 49.760398 | -22.011373 | -19.851340 | -26.043280 |
12 | 11.355555 | -10.285448 | -6.366640 | -77.454752 | -74.817511 | -68.622674 | 0.0 | 49.653160 | 36.833179 | -15.637405 | -15.169712 | -14.031923 |
13 | 9.942147 | -4.822680 | 8.399374 | -99.923514 | -99.767960 | -106.759770 | 0.0 | 22.136604 | -29.272696 | -15.316771 | -22.484958 | -20.448491 |
14 | 17.363018 | 26.569997 | -2.358907 | -102.144946 | -92.556346 | -83.513940 | 0.0 | 19.837680 | -2.366917 | 110.185077 | 120.484261 | 124.393057 |
15 | 47.359134 | -1.488056 | -21.988237 | 43.459463 | 61.501180 | 17.440002 | 0.0 | 0.792152 | 5.608713 | 193.763131 | 184.189201 | 192.893229 |
16 | -35.009986 | 6.411683 | -14.683955 | -90.871617 | -106.763795 | -77.705572 | 0.0 | -0.045683 | -13.577677 | 177.777245 | 167.599435 | 187.078731 |
17 | -24.231628 | -1.587288 | 1.203721 | -152.260419 | -153.805295 | -158.533883 | 0.0 | -8.353227 | 3.760762 | 64.707286 | 70.265285 | 60.804298 |
18 | 28.364028 | -6.513438 | -8.596102 | 108.529432 | 114.269032 | 106.140341 | 0.0 | 16.489773 | 17.921101 | 77.379549 | 62.262829 | 73.061897 |
19 | -24.711095 | -6.448676 | -7.051386 | -98.962835 | -98.935334 | -96.915881 | 0.0 | -9.155978 | -13.035711 | 35.139867 | 32.170966 | 41.962102 |
apinitial = []
for val in initial:
if val.startswith('r'):
apinitial.append(val.replace('r.','a.'))
else:
apinitial.append(val.replace('i.','p.'))
print(initial)
print(apinitial)
['r.1.0.0', 'r.1.1.0', 'r.1.1.1', 'r.1.2.0', 'r.1.2.1', 'r.1.2.2', 'i.1.0.0', 'i.1.1.0', 'i.1.1.1', 'i.1.2.0', 'i.1.2.1', 'i.1.2.2'] ['a.1.0.0', 'a.1.1.0', 'a.1.1.1', 'a.1.2.0', 'a.1.2.1', 'a.1.2.2', 'p.1.0.0', 'p.1.1.0', 'p.1.1.1', 'p.1.2.0', 'p.1.2.1', 'p.1.2.2']
apminuitresults = minuitresults.copy()
for i in range(len(minuitresults)):
# amp = sqrt(Re^2+Im^2)
apminuitresults.loc[i][0] = npy.sqrt(minuitresults.loc[i][0]**2+minuitresults.loc[i][6]**2)
apminuitresults.loc[i][1] = npy.sqrt(minuitresults.loc[i][1]**2+minuitresults.loc[i][7]**2)
apminuitresults.loc[i][2] = npy.sqrt(minuitresults.loc[i][2]**2+minuitresults.loc[i][8]**2)
apminuitresults.loc[i][3] = npy.sqrt(minuitresults.loc[i][3]**2+minuitresults.loc[i][9]**2)
apminuitresults.loc[i][4] = npy.sqrt(minuitresults.loc[i][4]**2+minuitresults.loc[i][10]**2)
apminuitresults.loc[i][5] = npy.sqrt(minuitresults.loc[i][5]**2+minuitresults.loc[i][11]**2)
# # phase = atan2(Im,Re)
apminuitresults.loc[i][6] = npy.arctan2(minuitresults.loc[i][6],minuitresults.loc[i][0])
apminuitresults.loc[i][7] = npy.arctan2(minuitresults.loc[i][7],minuitresults.loc[i][1])
apminuitresults.loc[i][8] = npy.arctan2(minuitresults.loc[i][8],minuitresults.loc[i][2])
apminuitresults.loc[i][9] = npy.arctan2(minuitresults.loc[i][9],minuitresults.loc[i][3])
apminuitresults.loc[i][10] = npy.arctan2(minuitresults.loc[i][10],minuitresults.loc[i][4])
apminuitresults.loc[i][11] = npy.arctan2(minuitresults.loc[i][11],minuitresults.loc[i][5])
apminuitresults.columns = apinitial
print(apminuitresults)
a.1.0.0 a.1.1.0 a.1.1.1 a.1.2.0 a.1.2.1 a.1.2.2 \ 0 154.600390 13.230801 1.108894 82.088730 77.056139 73.687329 1 194.342544 35.415926 19.437422 82.996023 84.775090 84.833321 2 281.450602 10.002888 11.445998 99.122515 85.803596 89.710780 3 500.845930 32.665749 14.156789 92.519466 124.448845 115.924304 4 423.063231 473.133938 416.366187 310.163990 346.061411 66.136781 5 506.378309 3.927792 13.591294 118.931404 127.688371 116.036180 6 238.013167 91.978237 71.043833 151.661315 151.484678 140.641594 7 165.091849 6.561844 5.992824 197.028104 197.415652 202.880175 8 87.954769 36.319293 74.229187 287.583883 291.385721 281.652626 9 12.548544 102.851899 80.267467 343.508217 343.777416 343.222009 10 42.958850 61.880626 11.646150 200.041195 201.912953 203.550498 11 17.028557 66.237710 51.284939 93.741571 93.796519 87.462745 12 11.355555 50.707266 37.379369 79.017510 76.339899 70.042603 13 9.942147 22.655849 30.453903 101.090613 102.270324 108.700456 14 17.363018 33.158683 3.341667 150.247600 151.931347 149.827270 15 47.359134 1.685768 22.692294 198.577128 194.185625 193.680023 16 35.009986 6.411846 19.999296 199.655703 198.716075 202.574943 17 24.231628 8.502698 3.948705 165.439621 169.095473 169.794448 18 28.364028 17.729566 19.876087 133.290031 130.130978 128.855783 19 24.711095 11.198989 14.820654 105.016441 104.034472 105.610160 p.1.0.0 p.1.1.0 p.1.1.1 p.1.2.0 p.1.2.1 p.1.2.2 0 0.000000 1.577383 -0.853511 0.392342 -0.281426 -0.274578 1 0.000000 1.602950 1.646997 0.014186 0.405617 0.049823 2 0.000000 -1.402794 1.515928 0.449501 -0.068878 0.058295 3 0.000000 1.587675 -2.093747 0.375650 0.841553 0.655389 4 0.000000 0.053340 -1.836378 2.359827 1.247635 1.372799 5 0.000000 -0.827247 1.656174 -2.333600 -2.278031 -2.376163 6 0.000000 1.610444 1.610689 -3.138662 -3.124628 3.112446 7 0.000000 0.385621 0.094974 2.557930 2.591840 2.499730 8 0.000000 1.906460 2.133737 -2.860953 -2.615316 -2.962552 9 3.141593 1.739261 1.728566 0.122632 0.143205 0.111826 10 3.141593 -2.072798 -1.664700 2.498242 2.591344 2.638124 11 3.141593 1.762350 1.815236 -2.904570 -2.928337 -2.839243 12 0.000000 1.775053 1.741956 -2.942380 -2.941548 -2.939894 13 0.000000 1.785305 -1.291368 -2.989492 -2.919924 -2.952347 14 0.000000 0.641334 -2.354499 2.318346 2.225846 2.162049 15 0.000000 2.652409 2.891841 1.350156 1.248534 1.480629 16 3.141593 -0.007125 -2.395319 2.043328 2.137991 1.964476 17 3.141593 -1.758578 1.261026 2.739741 2.713060 2.775355 18 0.000000 1.946984 2.018043 0.619386 0.498903 0.602866 19 3.141593 -2.184411 -2.066648 2.800398 2.827206 2.732987
# import AmplitudeJPACfit
# amplitude = AmplitudeJPACfit.FitAmplitude(initial)
import AmplitudeJPACfitAngles
amplitude = AmplitudeJPACfitAngles.FitAmplitude(apinitial)
Here the user defines number of bins, variable to be binned and range
#Define number of bins
nbins = len(minuitresults)
binsda = pwa.bin_by_range(data, "mass", nbins, .6, 2.0)
binsma = pwa.bin_by_range(mc, mc["mass"], nbins, .6, 2.0)
Plot events per bin for visualisation and to make sure there are enough events per bin
counts=[]
centers=[]
for bin in binsda:
counts.append(len(bin))
centers.append(npy.mean(bin["mass"]))
# Add yerr to argment list when we have errors
yerr = npy.empty(len(binsda))
yerr = npy.sqrt(counts)
plt.errorbar(centers,counts, yerr, linestyle="dashed", fmt="o",label='data')
plt.legend(loc='upper right')
plt.show()
See if we can improve results by that
mcmcresults = []
numwalker=100
numsteps=1000
ncpu=6
for index,dbin in enumerate(binsda):
with pwa.LogLikelihood(
amplitude, dbin, binsma[index], generated_length=len(binsma[index]), num_of_processes=ncpu, is_minimizer=False) as Likelihood:
print(f"Find random start parameters for fit #{index}...")
startpars = npy.zeros((numwalker,len(apinitial)))
for i in range(numwalker):
a=0
while True:
a+=1
#add some randomization here for moves to work better
startvals=apminuitresults.loc[index].copy()
for p in range(0,6):
startvals[p] = startvals[p] + 2*npy.random.rand()
for p in range(7,12): #initialise the phases with an overall offset to go round on the circle
startvals[p] = (startvals[p] + 0.1*npy.random.randn()) % (2*npy.pi)
# startvals[p] = (startvals[p]+i*2*npy.pi/numwalker) % (2*npy.pi)
startvals[6]=0 #fix phase for s wave
nll = Likelihood(startvals.to_dict())
if not npy.any(npy.isnan(nll)):
break
# print(f"After {a} steps found nll={nll}.") #debugging
startpars[i] = list(startvals)
cov = 1 #npy.array([0.002,0.002,0.002,0.005,0.005,0.005,0.005,0.005,0.005,1E5])
# parlimits = [(-500,1500) for i in range(len(initial))]
# parlimits[6] = (0,0)
parlimits = [(0,1500),(0,1500),(0,1500),(0,1500),(0,1500),(0,1500),
(0,0),(0,2*npy.pi),(0,2*npy.pi),(0,2*npy.pi),(0,2*npy.pi),(0,2*npy.pi)] #fix phase
# parlimits = [(minuitresults.loc[index][0],minuitresults.loc[index][0]),(minuitresults.loc[index][1],minuitresults.loc[index][1]),(minuitresults.loc[index][2],minuitresults.loc[index][2]),(minuitresults.loc[index][3],minuitresults.loc[index][3]),(minuitresults.loc[index][4],minuitresults.loc[index][4]),(minuitresults.loc[index][5],minuitresults.loc[index][5]),
# (0,2*npy.pi),(0,2*npy.pi),(0,2*npy.pi),(0,2*npy.pi),(0,2*npy.pi),(0,2*npy.pi)] #fix i.1.0.0 to 10 to fix phase
print(f"Start chain for bin {index} ...")
# print(startpars)
optimizer = pwa.mcmc(apinitial,Likelihood,
nsteps=numsteps,
startpars=startpars,
parlimits=parlimits,
nwalker=numwalker,
emceemoves=[(emcee.moves.StretchMove(), 0.5), (emcee.moves.DEMove(), 0.3), (emcee.moves.DESnookerMove(), 0.2)])
# emceemoves=emcee.moves.StretchMove())
# emceemoves=emcee.moves.GaussianMove(cov,'vector'))
mcmcresults.append(optimizer)
print("Mean acceptance fraction: {0:.3f}".format(npy.mean(optimizer.acceptance_fraction)))
try:
print("Mean autocorrelation time: {0:.3f} steps".format(npy.mean(optimizer.get_autocorr_time())))
except:
print("Couldn't get autocorrelation time.")
# break #only one bin
Find random start parameters for fit #0... Start chain for bin 0 ...
100%|██████████| 1000/1000 [02:29<00:00, 6.70it/s]
Mean acceptance fraction: 0.164 Couldn't get autocorrelation time.
Find random start parameters for fit #1... Start chain for bin 1 ...
100%|██████████| 1000/1000 [01:49<00:00, 9.12it/s]
Mean acceptance fraction: 0.065 Couldn't get autocorrelation time.
Find random start parameters for fit #2... Start chain for bin 2 ...
100%|██████████| 1000/1000 [01:45<00:00, 9.50it/s]
Mean acceptance fraction: 0.056 Couldn't get autocorrelation time.
Find random start parameters for fit #3... Start chain for bin 3 ...
100%|██████████| 1000/1000 [03:57<00:00, 4.22it/s]
Mean acceptance fraction: 0.182 Couldn't get autocorrelation time. Find random start parameters for fit #4... Start chain for bin 4 ...
100%|██████████| 1000/1000 [05:54<00:00, 2.82it/s]
Mean acceptance fraction: 0.120 Couldn't get autocorrelation time.
Find random start parameters for fit #5... Start chain for bin 5 ...
100%|██████████| 1000/1000 [03:22<00:00, 4.94it/s]
Mean acceptance fraction: 0.167 Couldn't get autocorrelation time. Find random start parameters for fit #6... Start chain for bin 6 ...
100%|██████████| 1000/1000 [04:28<00:00, 3.73it/s]
Mean acceptance fraction: 0.207 Couldn't get autocorrelation time.
Find random start parameters for fit #7... Start chain for bin 7 ...
100%|██████████| 1000/1000 [02:43<00:00, 6.10it/s]
Mean acceptance fraction: 0.160 Couldn't get autocorrelation time. Find random start parameters for fit #8... Start chain for bin 8 ...
100%|██████████| 1000/1000 [04:47<00:00, 3.48it/s]
Mean acceptance fraction: 0.198 Couldn't get autocorrelation time. Find random start parameters for fit #9... Start chain for bin 9 ...
100%|██████████| 1000/1000 [03:21<00:00, 4.96it/s]
Mean acceptance fraction: 0.106 Couldn't get autocorrelation time.
Find random start parameters for fit #10... Start chain for bin 10 ...
100%|██████████| 1000/1000 [02:47<00:00, 5.96it/s]
Mean acceptance fraction: 0.128 Couldn't get autocorrelation time. Find random start parameters for fit #11... Start chain for bin 11 ...
100%|██████████| 1000/1000 [03:03<00:00, 5.45it/s]
Mean acceptance fraction: 0.173 Couldn't get autocorrelation time.
Find random start parameters for fit #12... Start chain for bin 12 ...
100%|██████████| 1000/1000 [03:27<00:00, 4.82it/s]
Mean acceptance fraction: 0.194 Couldn't get autocorrelation time.
Find random start parameters for fit #13... Start chain for bin 13 ...
100%|██████████| 1000/1000 [03:00<00:00, 5.53it/s]
Mean acceptance fraction: 0.178 Couldn't get autocorrelation time.
Find random start parameters for fit #14... Start chain for bin 14 ...
100%|██████████| 1000/1000 [02:44<00:00, 6.08it/s]
Mean acceptance fraction: 0.157 Couldn't get autocorrelation time.
Find random start parameters for fit #15... Start chain for bin 15 ...
100%|██████████| 1000/1000 [02:48<00:00, 5.94it/s]
Mean acceptance fraction: 0.141 Couldn't get autocorrelation time. Find random start parameters for fit #16... Start chain for bin 16 ...
100%|██████████| 1000/1000 [03:19<00:00, 5.02it/s]
Mean acceptance fraction: 0.165 Couldn't get autocorrelation time.
Find random start parameters for fit #17... Start chain for bin 17 ...
100%|██████████| 1000/1000 [02:42<00:00, 6.16it/s]
Mean acceptance fraction: 0.141 Couldn't get autocorrelation time.
Find random start parameters for fit #18... Start chain for bin 18 ...
100%|██████████| 1000/1000 [02:53<00:00, 5.78it/s]
Mean acceptance fraction: 0.166 Couldn't get autocorrelation time.
Find random start parameters for fit #19... Start chain for bin 19 ...
100%|██████████| 1000/1000 [02:08<00:00, 7.81it/s]
Mean acceptance fraction: 0.090 Couldn't get autocorrelation time.
for i in range(len(binsda)):
samples = mcmcresults[i].get_chain(discard=0)
pwa.write(f"{outputpath}/mcmc_mix_{numwalker}_{numsteps}_bin{i}.npy",samples)
# just save one bin
pwa.write(f"{outputpath}/mcmc_mix_50k_bin0.npy",mcmcresults[0].get_chain(discard=0))
burnin=500
for index, dbin in enumerate(binsda):
# print(samples[:,:,0])
print(f"Show results for bin #{index}")
samples = mcmcresults[index].get_chain(discard=0)
numsteps = len(samples)
numpars = samples.shape[2]
fig, axes = plt.subplots(len(initial), figsize=(10, 70), sharex=True)
for i in range(len(initial)):
ax = axes[i]
ax.plot(samples[:,:,i])
ax.set_ylabel(apinitial[i])
ax.set_xlim(0, numsteps)
axes[-1].set_xlabel("step number");
fig.show()
samples = mcmcresults[index].get_chain(discard=burnin)
fig2, axes2 = plt.subplots(len(initial), figsize=(5, 70), sharex=False)
for i in range(len(initial)):
ax = axes2[i]
ax.hist(samples[:,:,i], 100,histtype='barstacked')
ax.set_ylabel(apinitial[i])
fig2.show()
flatsample = mcmcresults[index].get_chain(discard=burnin,flat=True)
varrange = []
for i in range(len(initial)):
varrange.append((flatsample[:,i].min(),flatsample[:,i].max()))
fig3 = corner.corner(flatsample,
color='royalblue',
bins=50,
# truths=npy.mean(flatsample,axis=0),
# quantiles=minuitresults.loc[index],
range=varrange,
labels=apinitial,
fill_contours=True,
truth_color='red',
label_kwargs={'fontsize':20, 'labelpad':20},
hist_kwargs = {'histtype':'stepfilled','alpha':1})
fig3.show()
# Draw mean and minuit values (from emcee documentation)
# Extract the axes
axes = npy.array(fig3.axes).reshape((numpars, numpars))
# Loop over the diagonal
for i in range(numpars):
ax = axes[i, i]
# ax.axvline(npy.mean(flatsample,axis=0)[i], color="k")
ax.axvline(apminuitresults.loc[index][i], color="r")
# Loop over the histograms
for yi in range(numpars):
for xi in range(yi):
ax = axes[yi, xi]
# ax.axvline(npy.mean(flatsample,axis=0)[xi], color="k")
ax.axvline(apminuitresults.loc[index][xi], color="r")
# ax.axhline(npy.mean(flatsample,axis=0)[yi], color="k")
ax.axhline(apminuitresults.loc[index][yi], color="r")
# ax.plot(npy.mean(flatsample,axis=0)[xi], npy.mean(flatsample,axis=0)[yi], "sk")
ax.plot(apminuitresults.loc[index][xi], apminuitresults.loc[index][yi], "sr")
break #only one bin for now
Show results for bin #0
WARNING:root:Too few points to create valid contours WARNING:root:Too few points to create valid contours