-
Notifications
You must be signed in to change notification settings - Fork 0
/
RatioCorrelators.py
executable file
·3928 lines (3337 loc) · 169 KB
/
RatioCorrelators.py
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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/env python
#IMPORT THIS FIRST, put in base import file
# import matplotlib
# matplotlib.use('Agg') # Must be before importing matplotlib.pyplot or pylab!
# import matplotlib.pyplot as pl
##
from Params import nboot,outputdir,Debug_Mode,Wipe_All_Fits,this_dir
from FileIO import ReadFuns,WriteFuns
# from PredefFitFuns import ConstantFitFun,RF_Prop_Fit
from PredefFitFuns import ConstantFitFun
import SetsOfFits as sff
from NullPlotData import null_series
from TwoPtCorrelators import NNQCorr
# from TwoPtCorrelators import TwoPointCorr,NNQCorr,NNQFullCorr
# from FlowOpps import FlowOp
from MiscFuns import ODNested,mkdir_p,CheckClass,CombineCfgs,op_str_list,MultiSplit,op_dict
from MiscFuns import flip_fun_2arg,GenSinkFun,CreateNameAndLL,GetInfoList,Series_fix_key
from XmlFormatting import untstr,tstr,unxmlfitr,AvgStdToFormat,AvgFormat
from XmlFormatting import rm_tsumfitr,unxmlfitr_int,KeyForamtting
# from XmlFormatting import rm_tsumfitr,tflowTOLeg,tsumTOLeg,tsinkTOLeg
from FileIO import WriteXml,WritePickle,ReadPickleWrap
from FileIO import WriteExcel,Construct_File_Object
from copy import deepcopy
# import cPickle as pickle
# import dill as pickle
import pandas as pa
# import warnings
# import itertools
import operator as op
import os
import numpy as np
from SetsOfTwoPtCorrs import SetOfTwoPt,SetOfNNQ
from ThreePtCorrelators import ThreePointCorr,NJNQCorr,NJNQFullCorr
# import FitFunctions as ff
from collections import OrderedDict
# from FlowOpps import FlowOp
"""
How to use:
Ratio(cfglist=[], Info={})
creates class using:
configuration list cfglist
Information Dictionary Info (See below for what can be in Info)
.LoadPickle()
Does everything basically.
"""
## function to center the x-axis around zero
def CenterT(tlist):
return tlist - np.mean(tlist),-tlist[0]+np.mean(tlist)+1
## for some reason, CROMA ouptut of the three-point correlator is shifted by 3 in time compared to the two-point correlator
twopttshift = 1
class RatioCorr(object):
"""
Rat(q,tau,pp,t) uses bootstrap class
pp = sink momentum
t = source sink separation time
q = current insersion momentum
tau = current insersion time
Info is a dictionary containing information for the correlator:
see comments below to see what is required.
missing elements just sets field to default
"""
def Construct_Rat_File(this_file):
return Construct_File_Object(this_file,RatioCorr)
## TODO: create variational flag and interface with this function.
## array of all classes that inherit this class. Append this when adding new classes please!
parentlist = ['SetOfCorrs.SetsOfRatios','RatioCorrelators.RatioFOCorr','RatioCorrelators.RatioFOFullCorr']
def __init__(self, thisnboot=nboot, cfglist={}, Info={},name='',
filename='',thissym='Not Set',thiscol='Not Set',thisshift=0.0,man_load_cfgs=False):
self.current = 0
## Initialising plotting parameters
self.thisshift = thisshift
self.thiscol = thiscol
self.thissym = thissym
## the actual shift ammount needs to be scaled by the x axis size
## see GetShift
self.thisshiftscale = 'Not Set'
self.thischecklist = ['C3file','Sym']
self.C3class = ThreePointCorr(thisnboot=thisnboot, cfglist=cfglist, Info=Info,name=name,filename=filename)
self.C3file = self.C3class.PickleFile
self.checkmomlist = self.C3class.checkmomlist
self.latparams = self.C3class.latparams
if 'pmom' not in list(Info.keys()):
Info['pmom'] = [self.latparams.GetAvgmomstr(ip) for ip in self.C3class.qmom]
else:
if isinstance(Info['pmom'],str):
if Info['pmom'] == 'All':
Info['pmom'] = self.latparams.GetMomList(list_type='str_Avg',list_order='squared')
else:
Info['pmom'] = [Info['pmom']]
Info['pmom'] = list(Info['pmom'])
for ip in self.C3class.qmom:
fmt_p = self.latparams.GetAvgmomstr(ip)
if fmt_p not in Info['pmom']:
Info['pmom'].append(fmt_p)
Info['pmom'] = self.latparams.MomSqrdSort(Info['pmom'])
if '0 0 0' not in Info['pmom']:
Info['pmom'] = ['0 0 0'] + Info['pmom']
if 'AllC2p' in list(Info.keys()) and Info['AllC2p']:
Info['pmom'] = self.latparams.GetMomList(list_type='str_Avg',list_order='squared')
self.nt = self.latparams.nt
self.nxyz = self.latparams.nxyz
self.dim_label = 'RC'+str(self.nxyz)+'x'+str(self.nt)
self.dim_ll = str(self.nxyz)+r'^{3}\times ' + str(self.nt)
# Info['pmom'] = self.latparams.momstrlistAvg
## both current and sink momentum are requred for the two point correlator construction when constructing the ratio function
## but the two point correlator has been averaged over, so we need to get the averaged momentum anaolgy
## bah, this was too hard to keep unique files for the two-point correlator.
## easier just to read and write the whole 2 point mom file.
## e.g. q = [-1,1,-1] corresponds to p = [1 1 1] in the two-point correlator
# ppforp = self.latparams.GetAvgmomstr(Info['ppmom'])
# if ppforp not in Info['pmom']: Info['pmom'].append(ppforp)
if 'show_cfgs' in list(Info.keys()): self.show_cfgs = Info['show_cfgs']
else: self.show_cfgs = True
## used to pick out different configuraiton/source lists.
## for example, if there are less flowed operators than 2 point correlators, passing 'TopCharge' will read that result
if 'fo_for_cfgs' in list(Info.keys()): self.fo_for_cfgs = Info['fo_for_cfgs']
else: self.fo_for_cfgs = False
## Does symmetrising over smearings
if 'Symmetric' in list(Info.keys()): self.Sym = Info['Symmetric']
else: self.Sym = False ## default to no symmetrizing, (usually full matrix isn't computed nowerdays.
self.sinklab,self.jsmcoeffs,self.sm2ptList,InfoList = GetInfoList(deepcopy(Info))
# self.C2filelist = [ikey.PickleFile for ikey in self.C2classlist.itervalues()]
self.C2Set = SetOfTwoPt(cfglist=cfglist,InfoDict = InfoList)
self.C2class = 'Not Set'
self.defgammamomlist = OrderedDict()
for igamma in self.C3class.gammalist:
self.defgammamomlist[igamma] = self.C3class.qform
self.Rat_Col_Names = ['current_operator','momentum','source_sink_separation']
self.Rat_Stats = pa.DataFrame()
'''
self.Rat_Stats DataFrame:
columns = boot, Avg, Std
rows = current_operator, momenta, time separation
'''
self.Rat_Fit_Col_Names = ['current_operator','momentum']
self.Rat_Fit_Stats = pa.DataFrame()
'''
self.Rat_Fit_Stats DataFrame:
columns = boot
rows = current_operator, momenta
Avg/Std elements are pa.Series instances with:
index = fit parameters
values = Avg/Std values
boot elements are SetOfFitFuns objects
'''
# # Ratboot [ igamma , ip , it ] bs
# self.Ratboot = ODNested()
# self.RatAvg = ODNested()
# self.RatStd = ODNested()
#
# ## self.FitDict = { igamma , ip , ifit } = thisff
# self.FitDict = ODNested()
# self.FitDictAvg = ODNested()
# self.FitDictStd = ODNested()
self.nboot = thisnboot
self.centert = 'Not Set'
self.cfglist = cfglist
## Info can contain fit ranges to be calculated
## must bein form Info['FitsRF'] = [ 'fit#-#' , ...]
if 'FitsRF' in list(Info.keys()): self.PredefFits = Info['FitsRF']
else: self.PredefFits = []
if 'RFFitFun' in list(Info.keys()): self.Fun,self.npar = Info['RFFitFun']
else: self.Fun,self.npar = (ConstantFitFun,1)
## Info can contain initial guess for fits, make sure it array of correct length
if 'iGuessRF' in list(Info.keys()): self.iGuess = np.array(Info['iGuessRF'])
else: self.iGuess = 'None' ## Defaults to value set for function in FitFunctions.py
## is needed?
self.Info = Info
self.SqrtFact = 'Not Set'
self.SetCustomName(string=name,stringfn=filename)
if not man_load_cfgs:
self.LoadCfgs(cfglist,name,filename)
def LoadCfgs(self,cfglist,name='',filename=''):
if len(cfglist) == 0:
self.CombineCfgLists()
else:
self.ImportCfgList(cfglist)
self.SetCustomName(string=name,stringfn=filename)
## just for consistancy with other classes
def GetCfgList(self):
self.CombineCfgLists()
def CombineCfgLists(self):
if isinstance(self.C2class,str):
thiscfglist = CombineCfgs(self.C3class.C3_cfgs,self.C2Set.cfglist)
else:
thiscfglist = CombineCfgs(self.C3class.C3_cfgs,self.C2class.C2_cfgs)
self.ImportCfgList(thiscfglist)
def SetCustomName(self,string='',stringfn='',stringLL=''):
if string == '':
self.name = 'Rat_'+self.C3class.name
else:
self.name = string
if stringfn == '':
self.filename = 'Rat_'+self.C3class.filename
else:
self.filename = stringfn
if stringLL == '':
# self.LegLab = '$'+'\ '.join([self.C3class.DS,self.C3class.ism,self.C3class.jsm])+'$' ## customise this how you like
# self.LegLab = '$'+'\ '.join([self.latparams.GetPionMassLab(),self.C3class.DS,self.C3class.ism,self.C3class.jsm,self.C3class.Proj,self.C3class.gammalist[0]])+'$' ## customise this how you like
self.LegLab = '$'+'\ '.join([ self.dim_ll,self.latparams.GetPionMassLab(),self.C3class.stream_name,
self.C3class.DS,self.C3class.Proj])+'$' ## customise this how you like
else:
self.LegLab = stringLL
self.SubDir = '_'.join([self.C3class.Proj,self.C3class.DS,self.C3class.ppform])
if isinstance(self.fo_for_cfgs,str):
self.Ratdir = outputdir + '/'+self.dim_label+self.C3class.kappafolder+'/Rat/'+self.fo_for_cfgs+'/'+self.SubDir+'/'
else:
self.Ratdir = outputdir + '/'+self.dim_label+self.C3class.kappafolder+'/Rat/'+self.SubDir+'/'
mkdir_p(self.Ratdir+'/Pickle/')
mkdir_p(self.Ratdir+'/Excel/')
self.HumanFile = self.Ratdir+self.filename+'.xml'
self.ExcelFile = self.Ratdir+'/Excel/'+self.filename+'.xlsx'
self.PickleFile = self.Ratdir+'/Pickle/'+self.filename+'.py3p'
def ImportCfgList(self,cfglist):
self.cfglist = cfglist
self.C2Set.ImportCfgList(cfglist)
if self.C2class != 'Not Set': self.C2class.ImportCfgList(cfglist)
self.C3class.ImportCfgList(cfglist)
self.stream_name = self.C3class.stream_name
def Stats(self):
if hasattr(self.C2class,'Stats'):
self.C2class.Stats()
if hasattr(self.C3class,'Stats'):
self.C3class.Stats()
if 'boot' not in self.Rat_Stats: return
lAvg,lStd = [],[]
for (igamma,ip,it),iRat in self.items():
iRat.Stats()
lAvg.append(iRat.Avg)
lStd.append(iRat.Std)
indicies = pa.MultiIndex.from_tuples(self.Rat_Stats.index,names=self.Rat_Col_Names)
self.Rat_Stats.loc[:,'Avg'] = pa.Series(lAvg,index=indicies)
self.Rat_Stats.loc[:,'Std'] = pa.Series(lStd,index=indicies)
def CombineWithAlpha(self,alpha_in,flow_time,this_fun=op.truediv):
tsink_str = self.C3class.tsink.replace('tsink','t')
if isinstance(this_fun,str):
this_fun = op_dict[this_fun.replace('alpha_rat_','').replace('_corr','')]
if not isinstance(alpha_in,NNQCorr):
print(type(alpha_in))
print(alpha_in)
raise AttributeError('Please pass in NNQCorr class instance')
if 'boot' not in self.Rat_Stats:
raise AttributeError('please compute Rat_Stats before combining with alpha')
if 'Alphaboot' not in alpha_in.NNQ_Stats :
raise AttributeError('please compute NNQ_Stats in alpha before combining with alpha')
if ('p000',tsink_str,flow_time) not in alpha_in.NNQ_Stats['Alphaboot']:
print(alpha_in.NNQ_Stats['Alphaboot'])
raise AttributeError('alpha ratio needs p000,'+tsink_str+','+flow_time+' result')
alpha_p0 = alpha_in.NNQ_Stats.loc[('p000',tsink_str,flow_time),'Alphaboot']
vall,ilist,corrl = [],[],[]
for (igamma,ip,it),irat in self.Rat_Stats['boot'].items():
ilist.append((igamma,ip,flow_time,it))
output = this_fun(irat,alpha_p0[('p000',tsink_str,flow_time)])
output.Stats()
vall.append(output)
corrl.append(irat.corr(alpha_p0[('p000',tsink_str,flow_time)]))
if len(ilist) > 0:
indicies = pa.MultiIndex.from_tuples(ilist,names=self.Rat_Col_Names)
self.Rat_Stats.loc[:,'alpha_rat_'+this_fun.__name__] = pa.Series(vall,index=indicies)
self.Rat_Stats.loc[:,'alpha_rat_'+this_fun.__name__+'_corr'] = pa.Series(corrl,index=indicies)
self.Write()
def RemoveVals(self):
# for iC2 in self.C2classlist.itervalues():
# iC2.RemoveVals()
self.C2class.RemoveVals()
self.C3class.RemoveVals()
## self.Rat = [ igamma , ip, it , iconf ]
## Configs must be the same, so we must check for them.
def Read(self,WipeData=True,DefWipe=False,CheckMom=True):
if isinstance(self.C2class,str): self.C2Set.LoadPickle(DefWipe=DefWipe,CheckCfgs=True,CheckMom=CheckMom,NoFit=True)
self.C3class.LoadPickle(WipeData=WipeData,DefWipe=DefWipe,CheckCfgs=True,CheckMom=CheckMom)
def CreateC2Sink(self):
if isinstance(self.C2class,str):
thisname,thisleglab,jsmlab = CreateNameAndLL(self.C2Set.SetC2,self.sinklab)
thisfun = GenSinkFun(self.jsmcoeffs)
self.C2class = self.C2Set.CombineCorrs(thisfun,filename=thisname,LegLab=thisleglab,jsmlab=jsmlab)
self.C2Set = 'Complete'
def Bootstrap(self,WipeData=True,DefWipe=False):
if isinstance(self.C2class,str) :
self.Read(WipeData=WipeData,DefWipe=DefWipe)
self.CreateC2Sink()
# self.C2class.Bootstrap(WipeData=WipeData)
self.C3class.Bootstrap(WipeData=WipeData,DefWipe=DefWipe)
def CalcRatio(self):
self.CalcSqrtFac()
## if DS has a divide in it, we ignore the sqrt factor as they cancel
div_bool = 'divide' in self.C3class.DS
## double the square root factor if we multiply two ds types
mult_bool = 'multiply' in self.C3class.DS
lRat,lRatAvg,lRatStd = [],[],[]
ilist = []
tsinkint = int(self.C3class.tsink.replace('tsink',''))
for (igamma,ip,it),tC3 in self.C3class.C3_Stats['boot'].items():
ip2pt = self.latparams.GetAvgmomform(ip,pref='p')
ict = int(untstr(it))
if ict > tsinkint: continue
if (ip2pt,it) not in self.SqrtFact.index: continue
# warnings.warn('DEBUG, please check if tsink for 2-pt corr is correct via g4.')
# print igamma, ip, it, self.C3class.tsink
# print 'SqrtFact'
# print self.SqrtFact
# print list(self.SqrtFact.index), ip2pt,it
# print type(self.SqrtFact[ip2pt,it])
# print self.SqrtFact[ip2pt,it]
# print self.SqrtFact[ip2pt,it].Avg
if div_bool:
this_boot = tC3
elif mult_bool:
this_boot = tC3*(self.SqrtFact[ip2pt,it]**2)
else:
this_boot = tC3*self.SqrtFact[ip2pt,it]
this_boot.Stats()
ilist.append((igamma,ip,it))
lRat.append(this_boot)
lRatAvg.append(this_boot.Avg)
lRatStd.append(this_boot.Std)
# self.Ratboot[igamma][ip][it] = tC3/self.C2class[ip.replace('q','p'),'t'+str(int(self.C3class.tsink.replace('tsink',''))+1)]
if len(ilist) > 0:
indicies = pa.MultiIndex.from_tuples(ilist,names=self.Rat_Col_Names)
self.Rat_Stats.loc[:,'boot'] = pa.Series(lRat,index=indicies)
self.Rat_Stats.loc[:,'Avg'] = pa.Series(lRatAvg,index=indicies)
self.Rat_Stats.loc[:,'Std'] = pa.Series(lRatStd,index=indicies)
def CalcSqrtFac(self):
SqrtFact_hold = []
ilist = []
ipp = self.latparams.GetAvgmomform(self.C3class.ppform,pref='p')
if Debug_Mode:
print()
print('Outputting sqrt factor debug info:')
print()
for iq in self.C3class.qform:
ip = self.latparams.GetAvgmomform(iq,pref='p')
sqrt_fac = self.CalcSqrtFacMom(self.C2class.C2_Stats['boot'][ipp],self.C2class.C2_Stats['boot'][ip],self.C3class.tsink)
for it,tdata in enumerate(sqrt_fac):
if (ip,tstr(it+1)) not in ilist:
if Debug_Mode:
for ierr in tdata.GetNaNIndex():
print(iq, tstr(it+1), ' iboot'+str(ierr).zfill(4))
SqrtFact_hold.append(tdata)
ilist.append((ip,tstr(it+1)))
# if neg_sqrt: warnings.warn('negative square root in square root factor for '+ipp+iq)
indicies = pa.MultiIndex.from_tuples(ilist,names=['momentum','source_sink_separation'])
self.SqrtFact = pa.Series(SqrtFact_hold,index=indicies)
def CheckNegForSqrt(self,thisboot,negsqrt):
for ict in range(len(thisboot.bootvals)):
if thisboot.bootvals[ict] < 0.0:
negsqrt = True
# if G2_abs:
# thisboot.bootvals[ict] = abs(thisboot.bootvals[ict])
# else:
thisboot.bootvals[ict] = float('NaN')
# thisboot.bootvals[ict] = 0.0
return thisboot,negsqrt
## taking squareroot of each element before combinding makes the numbers not too small
def CalcSqrtFacMom(self,data2ptpp,data2ptp,thistsink):
dictout = []
thistsink = thistsink.replace('tsink','t')
thistsinkp1 = 't'+str(int(thistsink.replace('t',''))+twopttshift)
inttsink = int(thistsink.replace('t',''))
# NegSQRTRatFac = False
# data2ptpp[thistsinkp1],NegSQRTRatFac = self.CheckNegForSqrt(data2ptpp[thistsinkp1],NegSQRTRatFac)
# data2ptp[thistsinkp1],NegSQRTRatFac = self.CheckNegForSqrt(data2ptp[thistsinkp1],NegSQRTRatFac)
data2ptpp = data2ptpp.apply(lambda x : x.Sqrt())
data2ptp = data2ptp.apply(lambda x : x.Sqrt())
# NegSQRTRatFac = any(data2ptpp.isnull().values) or any(data2ptp.isnull().values)
two = data2ptpp[thistsinkp1]*data2ptp[thistsinkp1]
# two = two**(0.5)
thisshift = twopttshift
for it in range(thisshift,inttsink+thisshift):
thiststr = tstr(it)
tflip = tstr(inttsink-it+twopttshift+thisshift)
# data2ptpp[thiststr],NegSQRTRatFac = self.CheckNegForSqrt(data2ptpp[thiststr],NegSQRTRatFac)
# data2ptp[thiststr],NegSQRTRatFac = self.CheckNegForSqrt(data2ptp[thiststr],NegSQRTRatFac)
# data2ptpp[tflip],NegSQRTRatFac = self.CheckNegForSqrt(data2ptpp[tflip],NegSQRTRatFac)
# data2ptp[tflip],NegSQRTRatFac = self.CheckNegForSqrt(data2ptp[tflip],NegSQRTRatFac)
one = data2ptpp[thiststr]/data2ptp[thiststr]
three = data2ptp[tflip]/data2ptpp[tflip]
inter = one*three
ott = inter/two
# ott.Sqrt()
ott.Stats()
dictout.append(ott)
return dictout
def ImportFitRanges(self,fit_range):
if 'PreDef' == fit_range:
self.fit_range = self.PredefFits
else:
self.fit_range = fit_range
if not isinstance(self.fit_range,str):
print('SetsOfFits implementation now accepts 1 fit range to scan between. Selecting first one')
self.fit_range = self.fit_range[0]
self.fit_range = rm_tsumfitr(self.fit_range)
## NB, iGuess here is not implemented yet, needs to be passed into the function
def Fit(self,fit_range='PreDef',iGuess='PreDef',EstDir=False,WipeFit=True,show_timer=False,min_fitr=2):
self.ImportFitRanges(fit_range)
if iGuess != 'PreDef': self.iGuess = iGuess
lFit,ilist = [],[]
for (igamma,ip),pdata in self.Rat_Stats['boot'].groupby(level=('current_operator','momentum')):
this_key = (igamma,ip)
if this_key in self.Rat_Fit_Stats.index:
if WipeFit:
self.Rat_Fit_Stats.drop([(igamma,ip)],inplace=True)
else:
continue
if isinstance(self.fit_range,str):
ifitmin,ifitmax = np.array(unxmlfitr(self.fit_range),dtype=np.float64)
elif isinstance(self.fit_range,(tuple,list,np.ndarray)):
ifitmin,ifitmax = np.array(self.fit_range,dtype=np.float64)
else:
print(self.fit_range)
raise IOError('incorrect format for fitr#-# ' )
# ifitmin,ifitmax = np.array(unxmlfitr(ifit),dtype=np.float64)
# ##debugging
# print
# print 'FitData:'
# for it,iy in zip(tdatarange[0],ydata):
# print it, iy
# print
fit_info = {}
if EstDir:
fit_info['Funs'] = (self.Fun,self.npar,'Estimate')
else:
fit_info['Funs'] = (self.Fun,self.npar)
fit_info['iGuess'] = self.iGuess
fit_info['name'] = '_'.join(this_key)
thisff = sff.SetOfFitFuns(data=pdata[this_key],name=fit_info['name'])
thisff.ScanRange_fitr(self.fit_range,fit_info=fit_info,min_fit_len=min_fitr)
lFit.append(thisff)
ilist.append(this_key)
if len(ilist) > 0:
indicies = pa.MultiIndex.from_tuples(ilist,names=self.Rat_Fit_Col_Names)
if 'boot' in self.Rat_Fit_Stats.columns:
this_df = pa.DataFrame()
this_df.loc[:,'boot'] = pa.Series(lFit,index=indicies)
self.Rat_Fit_Stats = self.Rat_Fit_Stats.append(this_df)
else:
self.Rat_Fit_Stats.loc[:,'boot'] = pa.Series(lFit,index=indicies)
def DoFits(val):
val.DoFits(show_timer=show_timer)
return val
self.Rat_Fit_Stats['boot'] = self.Rat_Fit_Stats['boot'].apply(DoFits)
self.Write()
def GetFuns(self):
if self.Rat_Fit_Stats is not None and 'boot' in self.Rat_Fit_Stats.columns:
for ival in self.Rat_Fit_Stats['boot'].values:
ival.GetFuns()
if hasattr(self.C2class,'GetFuns'):
self.C2class.GetFuns()
if hasattr(self.C2Set,'GetFuns'):
self.C2Set.GetFuns()
if not hasattr(self,'Fun'):
self.Fun = ReadFuns(self.Fun_name)[0]
def RemoveFuns(self):
if self.Rat_Fit_Stats is not None and 'boot' in self.Rat_Fit_Stats.columns:
for ival in self.Rat_Fit_Stats['boot'].values:
ival.RemoveFuns()
if hasattr(self.C2class,'GetFuns'):
self.C2class.RemoveFuns()
if hasattr(self.C2Set,'GetFuns'):
self.C2Set.RemoveFuns()
if hasattr(self,'Fun'):
self.Fun_name = self.Fun.__name__
WriteFuns(self.Fun)
del self.Fun
def Write(self):
def FixDictArray(thislist,ikey):
outDict = OrderedDict()
for ic,ilist in enumerate(thislist):
outDict[ikey+str(ic+1)] = ilist
return outDict
# C2name = self.C2class.name.replace(self.C2class.jsm,self.C3class.jsm)
# C2LL = self.C2class.LegLab.replace(self.C2class.jsm,self.C3class.jsm)
# self.C2class.SetCustomName(C2name,C2LL)
if hasattr(self.C3class,'Write'):
self.C3class.Write()
if not hasattr(self,'show_cfgs'):
self.show_cfgs = True
outDict = ODNested()
outDict['name'] = self.name
if hasattr(self.C2class,'Write'):
self.C2class.Write()
outDict['C2file'] = self.C2class.HumanFile
outDict['C3file'] = self.C3class.HumanFile
outDict['kud'] = self.C3class.kud
outDict['ks'] = self.C3class.ks
outDict['nt'] = self.nt
outDict['nxyz'] = self.nxyz
for istream,incfg in zip(self.C3class.stream_list,self.C3class.ncfg_list):
outDict['ncfg_'+istream] = incfg
outDict['nxsrc_avg'] = self.C3class.xsrcAvg
outDict['nMeas'] = self.C3class.nMeas
outDict['nboot'] = self.nboot
outDict['tsrc'] = self.C3class.tsrc
outDict['tsink'] = self.C3class.tsink
# outDict['xsrc'] = FixDictArray(self.xsrc,'xsrc')
outDict['ism'] = self.C3class.ism
outDict['sink_type'] = self.C3class.jsm
outDict['Projector'] = self.C3class.Proj
outDict['Doub_Sing'] = self.C3class.DS
outDict['ppmom'] = self.C3class.ppmom
for icg,igamma in enumerate(self.C3class.gammalist):
outDict['gamma'+str(icg+1)] = igamma
for icq,iqmom in enumerate(self.C3class.qmom):
outDict['qmom'+str(icq+1)] = iqmom
excel_params = pa.Series(deepcopy(outDict))
outDict['gamma_list'] = FixDictArray(self.C3class.gammalist,'gamma')
outDict['qmom_list'] = FixDictArray(self.C3class.qmom,'qmom')
for icg,igamma in enumerate(self.C3class.gammalist):
del outDict['gamma'+str(icg+1)]
for icq,iqmom in enumerate(self.C3class.qmom):
del outDict['qmom'+str(icq+1)]
for (igamma,ip),fitdata in self.Rat_Fit_Stats['boot'].items():
outDict['Fits'][igamma][ip] = fitdata.GetOutputDict()
for col_key,col_data in self.Rat_Stats.items():
if 'Avg' in col_key or 'Std' in col_key: continue
for (igamma,ip,it),idata in col_data.items():
outDict[col_key][igamma][ip][it] = AvgStdToFormat(idata.Avg,idata.Std)
y = getattr(self, "SqrtFact", None)
if y is not None and hasattr(self.SqrtFact,'iteritems'):
for (ip,it),idata in self.SqrtFact.items():
if np.isnan(idata.Avg): idata.Stats()
outDict['SqrtFact'][ip][it] = AvgStdToFormat(idata.Avg,idata.Std)
for (ip,it),idata in self.C2class.C2_Stats['boot'].items():
if np.isnan(idata.Avg): idata.Stats()
outDict['G2'][ip][it] = AvgStdToFormat(idata.Avg,idata.Std)
if any(np.array(self.C3class.ncfg_list) > 0) and self.show_cfgs:
for ((istream,iccfg),icfg),ixsrc_list in zip(iter(self.C3class.C3_cfgs['configs'].items()),self.C3class.C3_cfgs['xsrc_list'].values):
outDict['cfglist'][istream][icfg] = 'n_xsrc_'+str(len(ixsrc_list))
## Write human readable file with data
WriteXml(self.HumanFile,{'Results':outDict})
if any(np.array(self.C3class.ncfg_list) > 0):
for ((istream,iccfg),icfg),ixsrc_list in zip(iter(self.C3class.C3_cfgs['configs'].items()),self.C3class.C3_cfgs['xsrc_list'].values):
outDict['cfglist'][istream][icfg] = 'n_xsrc_'+str(len(ixsrc_list))
## TODO, include fitting in the excell file
WriteExcel(self.ExcelFile,{'Rat_results':deepcopy(self.Rat_Stats)},cfglist=self.C3class.C3_cfgs,params=excel_params)
## pickles rest of data for reading
self.RemoveFuns()
WritePickle(self.PickleFile,self.__dict__)
self.GetFuns()
##################################
def ReadAndWrite(self,WipeData=True,DefWipe=False,CheckCfgs=False):
self.Read(WipeData=WipeData,DefWipe=DefWipe)
self.CreateC2Sink()
self.Bootstrap(WipeData=WipeData,DefWipe=DefWipe)
self.CalcRatio()
self.Fit()
self.Write()
## add things you want to update to the depreciated pickled files here.
def FixDict(self,thisdict):
thisdict['LegLab'] = self.LegLab
return thisdict
def GetTflowList(self,load_dict):
return load_dict['Rat_Stats'].index.levels[2]
def FixRead(self,readdict):
rewrite = False
if not isinstance(readdict['C2class'],str):
if not isinstance(readdict['C2Set'],str):
readdict['C2Set'] = 'Complete'
rewrite = True
return readdict,rewrite
def LoadPickle(self,DefWipe=False,WipeData = True,CheckCfgs=False,CheckMom=True):
# print
# print 'NJN'
# print self.HumanFile
# print self.PickleFile
if os.path.isfile(self.PickleFile) and not DefWipe:
print('Loading Pickle file ' , self.PickleFile)
loadeddict = ReadPickleWrap(self.PickleFile)
loadeddict,rewrite = self.FixRead(loadeddict)
checklist = self.thischecklist
if CheckCfgs: checklist = checklist + ['cfglist']
if CheckMom: checklist = checklist + ['checkmomlist']
if CheckClass(self.__dict__,loadeddict,checklist):
loadeddict = self.FixDict(loadeddict)
self.__dict__.update(loadeddict)
if rewrite: self.Write()
if Wipe_All_Fits:
self.Rat_Fit_Stats = pa.DataFrame()
elif len(self.Rat_Fit_Stats.values)> 0 and not isinstance(self.Rat_Fit_Stats.values[0][0],sff.SetOfFitFuns):
print('Legacy file, fits need to be recomputed')
self.Rat_Fit_Stats = pa.DataFrame()
else:
if os.path.isfile(self.PickleFile+'.bak'):
self.PickleFile = self.PickleFile+'.bak'
self.HumanFile = self.HumanFile+'.bak'
self.ExcelFile = self.ExcelFile.replace('.xlsx','.bak.xlsx')
print('using backupfile for '+self.HumanFile)
self.LoadPickle(DefWipe=DefWipe,WipeData=WipeData,CheckCfgs=CheckCfgs,CheckMom=CheckMom)
return
print('Warning, file ' , self.HumanFile , ' has different parameters than this instance, reading and writing over file now')
print()
self.ReadAndWrite(WipeData=WipeData,DefWipe=DefWipe)
# self.Write()
else:
if not os.path.isfile(self.PickleFile):
print()
print('pickle file not found: \n', self.PickleFile)
self.ReadAndWrite(WipeData=WipeData,DefWipe=DefWipe)
## Routines for Plotting ###############################################
def ImportPlotParams(self,thiscol,thissym,thisshift):
self.CheckCol(thiscol)
self.CheckSym(thissym)
self.thisshift = thisshift
def CheckCol(self,thiscol):
if 'PreDefine' == thiscol :
if 'Not Set' == self.thiscol:
raise IOError('Pass in color to initialize it')
else:
self.thiscol = thiscol
def CheckSym(self,thissym):
if 'PreDefine' == thissym :
if 'Not Set' == self.thissym:
raise IOError('Pass in symbol to initialize it')
else:
self.thissym = thissym
def GetShift(self,xlims,thisshift):
if thisshift != 'PreDefine': self.thisshift = thisshift
if self.thisshiftscale == 'Not Set':
xlen = np.abs(xlims[1]-xlims[0])
self.thisshiftscale = self.thisshift*xlen
return self.thisshiftscale
else:
return self.thisshiftscale
## make sure all xlims are the same, because this sets the window of xlim (needed for calculating shifts)
## gammamomlist is a dictionary where the keys are the gamma matricies, and the values are lists of momenta
## TODO: SHIFTS NEED TO BE FIXED!! LIKE NOW. A better way needs to be implemented.
def Plot(self,plot_class,xlims='tsink',gammamomlist={},thiscol='PreDefine',thissym='PreDefine',thisshift='PreDefine'):
if xlims == 'tsink': xlims = [0,int(self.C3class.tsink.replace('tsink',''))+1]
self.CheckCol(thiscol)
self.CheckSym(thissym)
self.thisshiftscale = 'Not Set'
thisshift = self.GetShift(xlims,thisshift)
if 'boot' not in self.Rat_Stats.columns:
raise EnvironmentError('no data to plot, maybe try to read in and bootstrap before plotting')
if len(list(gammamomlist.keys())) == 0:
gammamomlist = self.defgammamomlist
# for (igamma,ip),Rat_data in self.Rat_Stats['boot'].groupby(level=('current_operator','momentum')):
# if igamma not in gammamomlist.keys(): continue
# if ip not in gammamomlist[igamma]: continue
igamma = list(gammamomlist.keys())[0]
ip = self.latparams.TOpform(list(gammamomlist.values())[0][0],pref='q')
ip = ip.replace('qq','q')
if (igamma,ip) not in self.Rat_Stats['boot'].index:
print(', '.join([igamma,ip]) ,' not computed')
return plot_class
# tlist,databoot = zip(*list(self.Rat_Stats['boot'][igamma,ip].iteritems()))
this_key = (igamma,ip,slice(None))
databoot = self.Rat_Stats['boot']
dataavg = databoot.apply(lambda x : x.Avg)
dataerr = databoot.apply(lambda x : x.Std)
tlist = list(databoot[this_key].keys())
# dataavg,dataerr = [],[]
# for idata in databoot[xlims[0]:xlims[1]]:
# dataavg.append(idata.Avg)
# dataerr.append(idata.Std)
tlist = tlist[xlims[0]:xlims[1]]
# print 'plotting ' , igamma,ip
# pleg,gammaleg = '',''
# if len(gammamomlist[igamma]) > 1:
# pleg = ' $'+ip+'$'
# if '_cmplx' in igamma:
# gammaleg = ' $i'+igamma.replace('_cmplx','')+'$'
# else:
# gammaleg = ' $'+igamma+'$'
# if len(gammamomlist.keys()) > 1:
tlistplot,self.centert = CenterT(np.array(list(map(untstr,tlist)))-twopttshift+thisshift)
hold_series = pa.Series()
# hold_series['x_data'] = tlistplot
hold_series['x_data'] = 'from_keys'
hold_series['y_data'] = dataavg
hold_series['yerr_data'] = dataerr
hold_series['xerr_data'] = None
hold_series['key_select'] = this_key
hold_series['type'] = 'error_bar_vary'
hold_series['fit_class'] = None
hold_series['label'] = self.LegLab#+pleg+gammaleg
hold_series['symbol'] = self.thissym
hold_series['color'] = self.thiscol
hold_series['shift'] = thisshift
hold_series['xdatarange'] = None
hold_series['Phys'] = None
hold_series['ShowPar'] = None
hold_series['fmt_class'] = KeyForamtting(self.latparams)
plot_class.AppendData(hold_series)
return plot_class
## state is 1 or 2 state fit
## fitr is fit range e.g. fitr5-10
## fit range of 'Data' is just the original fit range
## tflowlist is just for consistancy with other function in NJNQ, does nothing
def GetBestFitr(self,igamma,ip):
if 'boot' not in self.Rat_Fit_Stats.columns or (igamma,ip) not in self.Rat_Fit_Stats.index:
errstr = ' '.join(igamma,ip)+ ' not foun in ratio fit list, cannot find best'
raise EnvironmentError(errstr)
chi_str,cut_fit_result = self.Rat_Fit_Stats.at[(igamma,ip),'boot'].GetBestChi()
return r'$\chi^{2}_{pdf}={:.3f}$'.format(chi_str),cut_fit_result
def Get_Extrapolation(self,keep_keys=slice(None),fmted=True):
if 'boot' in self.Rat_Fit_Stats:
return sff.PullSOFSeries_Wpar(self.Rat_Fit_Stats['boot'],keep_keys=keep_keys,fmted=fmted)
else:
return pa.Series()
def FitPlot(self,plot_class,fitr,xlims = 'Data',gammamomlist={},
thiscol='PreDefine',thisshift='PreDefine',tflowlist=[]):
self.CheckCol(thiscol)
thisshift = self.GetShift(xlims,thisshift)
if self.centert == 'Not Set':
print('print centert not set, defaulting to 0')
self.centert = 0.9
if 'boot' not in self.Rat_Fit_Stats.columns:
raise EnvironmentError('No fitting has been done, maybe try to Fit before plotting fit')
if len(list(gammamomlist.keys())) == 0:
gammamomlist = self.defgammamomlist
# for igamma,momlist in gammamomlist.iteritems():
# for ip in momlist:
igamma,momlist = next(iter(gammamomlist.items()))
ip = self.latparams.TOpform(momlist[0],pref='q')
ip = ip.replace('qq','q')
if fitr == 'Best':
chi_legend,fitr = self.GetBestFitr(igamma,ip)
if 'boot' not in self.Rat_Fit_Stats.columns:
print('no fits done, attempting fit now')
self.Fit(fit_range=fitr)
fit_data = sff.PullSOFSeries(self.Rat_Fit_Stats['boot'])
this_fun = fit_data.index[0][2]
this_key = (igamma,ip,this_fun,fitr)
if this_key not in fit_data:
print(this_key , 'not present in Rat Fit Stats, fitting now')
self.Fit(fit_range=fitr)
fit_data = sff.PullSOFSeries(self.Rat_Fit_Stats['boot'])
this_fun = fit_data.index[0][2]
this_key = (igamma,ip,this_fun,fitr)
hold_series = null_series
hold_series['key_select'] = this_key
hold_series['type'] = 'fit_vary'
hold_series['fit_class'] = fit_data
# hold_series['label'] = fit_data[this_key].name +' ' + chi_legend
hold_series['label'] = 'ACTUAL_VALUE_fit'
hold_series['symbol'] = self.thissym
hold_series['color'] = self.thiscol
hold_series['shift'] = thisshift
hold_series['xdatarange'] = xlims
hold_series['ShowPar'] = r'Par1'
hold_series['fmt_class'] = KeyForamtting(self.latparams)
plot_class.AppendData(hold_series)
return plot_class
## ##################### ###############################################
## Comparisons
def __str__(self):
return '_'.join([self.SubDir,self.name])
def __iter__(self):
return self
def iteritems(self):
return iter(self.Rat_Stats['boot'].items())
def items(self):
return list(self.Rat_Stats['boot'].items())
def itemsAvgStd(self):
outlist = []
for ivals in self.Rat_Stats.itertuples():
outlist.append((ivals[0],ivals[2],ivals[3]))
return outlist
def values(self):
return self.Rat_Stats['boot'].values
def itervalues(self):
return iter(self.Rat_Stats['boot'].values)
def valuesAvgStd(self):
outlist = []
for ivals in self.Rat_Stats.itertuples():
outlist.append((ivals[2],ivals[3]))
return outlist
def keys(self):
return self.Rat_Stats.index
def __setitem__(self, key, value ):
self.Rat_Stats['boot'][key] = value
def __getitem__(self, key):
return self.Rat_Stats['boot'][key]
def __reversed__(self):
self.Rat_Stats = self.Rat_Stats.iiloc[::-1]
return self
def __next__(self):
if self.current >= len(list(self.keys())):
self.current = 0
raise StopIteration
else:
self.current += 1
if self.current >= len(list(self.keys()))+1:
return 0.0
else:
thiskey = list(self.keys())[self.current-1]
return self[thiskey]
## unitary arithmetic operations
def __neg__(self):
for ival in list(self.keys()):
self[ival] = -self[ival]
return self
def __pos__(self):
return self
def __abs__(self):
for ival in list(self.keys()):
self[ival] = abs(self[ival])
return self
## Casting, just casts the bootstrap values, ignore return value (only needed since return has to be same type).
def __complex__(self):
for ival in list(self.keys()):
complex(self[ival])
self.Stats()
return 1+0j
def __int__(self):
for ival in list(self.keys()):
int(self[ival])
self.Stats()
return 1
def __float__(self):
for ival in list(self.keys()):
np.float64(self[ival])
self.Stats()
return 1.0
## BLOCK FOR OVERLOADING ARITHMETIC OPERATIONS ####
def ReduceIndicies(self):
indexl = []
for (igamma,ip,it) in list(self.Rat_Stats['boot'].keys()):
first_gamma = MultiSplit(igamma,op_str_list+['FO','FOFull'])[0]
this_gamma = first_gamma.replace('FO','')
this_gamma = this_gamma.replace('FOFull','')
this_gamma = this_gamma.replace('__','_')
if this_gamma[-1] == '_':
this_gamma = this_gamma[:-1]