# _____________________________________________________________________________ # /home/jetsu/Jaavso/progs/dcm.py # =============================== # DCM = Discrete Chi-square Method # _____________________________________________________________________________ # Version 2: The most important revisions are marked. # _____________________________________________________________________________ import os import numpy as np import pylab as pl from scipy import optimize os.system('clear') # _____________________________________________________________________________ # - Read input file 'dcm.dat' variables. # _____________________________________________________________________________ def dcmInput(Val): file0='dcm.dat' file = open(file0, 'r') rivi = file.readline() while (len(rivi) > 0): osat=rivi.split("=") luku=np.float(osat[0]) loppu=osat[2].strip() if ((Val == 1.) & (luku == 1.)): NY=loppu # Tag: String! if ((Val == 2.) & (luku == 2.)): NY=np.float(loppu) # RealData if ((Val == 3.) & (luku == 3.)): NY=loppu # file1: String if ((Val == 4.) & (luku == 4.)): NY=np.float(loppu) # Dummy if ((Val == 5.) & (luku == 5.)): NY=np.int(np.float(loppu)) # K1: Integer! if ((Val == 6.) & (luku == 6.)): NY=np.int(np.float(loppu)) # K2: Integer! if ((Val == 7.) & (luku == 7.)): NY=np.int(np.float(loppu)) # K3: Integer! if ((Val == 8.) & (luku == 8.)): NY=np.int(np.float(loppu)) # nL: Integer! if ((Val == 9.) & (luku == 9.)): NY=np.int(np.float(loppu)) # nS: Integer! if ((Val == 10.) & (luku == 10.)): NY=np.float(loppu) # c if ((Val == 11.) & (luku == 11.)): NY=np.float(loppu) # TestStat if ((Val == 12.) & (luku == 12.)): NY=np.float(loppu) # PMIN if ((Val == 13.) & (luku == 13.)): NY=np.float(loppu) # PMAX if ((Val == 14.) & (luku == 14.)): NY=np.int(np.float(loppu)) # Rounds: Integer! if ((Val == 15.) & (luku == 15.)): NY=np.float(loppu) # NonLinear if ((Val == 16.) & (luku == 16.)): NY=np.float(loppu) # SimT if ((Val == 17.) & (luku == 17.)): NY=np.int(np.float(loppu)) # SimN: Integer! if ((Val == 18.) & (luku == 18.)): NY=np.float(loppu) # SimSN if ((Val == 19.) & (luku == 19.)): NY=np.float(loppu) # SimDT if ((Val == 20.) & (luku == 20.)): NY=np.int(np.float(loppu)) # SimMany: Integer! if ((Val == 21.) & (luku == 21.)): NY=np.int(np.float(loppu)) # SimRounds: Integer! if ((Val == 22.) & (luku == 22.)): NY=np.float(loppu) # SimDF if ((Val == 23.) & (luku == 23.)): NY=np.float(loppu) # SimDA if ((Val == 24.) & (luku == 24.)): NY=np.float(loppu) # PrintScreen rivi = file.readline() # Read next line: Important in while! file.close() return NY # _____________________________________________________________________________ # - Indices for sorting vector in increasing order. # _____________________________________________________________________________ def IndexSortOrder(y): k=sorted(range(len(y)), key=lambda i:y[i]) return k # _____________________________________________________________________________ # - Frequencies for simulated data in decreasing order. # _____________________________________________________________________________ def SimulatedALLF(): K1 =dcmInput(5) PMIN=dcmInput(12) PMAX=dcmInput(13) FMIN=1./PMAX FMAX=1./PMIN ALLF=FMIN+np.random.uniform(0,1,K1)*(FMAX-FMIN) k=IndexSortOrder(-1.*ALLF) ALLF=ALLF[k] return ALLF # _____________________________________________________________________________ # - Storing/writing tested frequencies to ALLF.dat # _____________________________________________________________________________ def WriteALLF(ALLF): code1=open('ALLF.dat',"w") if (np.size(ALLF) < 1.5): code1.write("%15.15e\n" %(ALLF)) else: for i in range(np.size(ALLF)): code1.write("%15.15e\n" %(ALLF[i])) code1.close() return # _____________________________________________________________________________ # - Reading tested frequencies from ALLF.dat # _____________________________________________________________________________ def ReadALLF(): ALLF=-1. file = open('ALLF.dat', 'r') i=0 rivi = file.readline() while (len(rivi) > 0): rivi=rivi.rstrip('\n') ALLF=np.append(ALLF,np.float(rivi)) rivi = file.readline() ALLF=ALLF[1:] return ALLF # _____________________________________________________________________________ # - Simulated uniform random B_II free parameters between -0.5 and 0.5. # _____________________________________________________________________________ def LinSimulatedBeta(): K1=dcmInput(5) K2=dcmInput(6) K3=dcmInput(7) p=K1*2*K2+K3+1 BETA=0.5-np.random.uniform(0,1,p) return BETA # _____________________________________________________________________________ # - Linear model. # _____________________________________________________________________________ def LinearModel(T,BETA,Check): if (Check != 0): print('input BETA',BETA) K1=dcmInput(5) K2=dcmInput(6) K3=dcmInput(7) ALLF=ReadALLF() if (Check != 0): print('input ALLF',ALLF) H =0. H1=0. H2=0. H3=0. H4=0. H5=0. H6=0. PL=0. I=0 if (Check != 0): print('ALLF =',ALLF) print('1./ALLF =',1./ALLF) for i in range(1,K1+1): F=ALLF[i-1] X=2.*np.pi*T*F for j in range(1,K2+1): I=I+1 I1=I I=I+1 I2=I if (Check != 0): print('i',i,'j',j,'I1',I1,'I2',I2,'Beta',BETA[I1-1],BETA[I2-1]) if (i==1): H1=H1+BETA[I1-1]*np.cos(1.*j*X)+BETA[I2-1]*np.sin(1.*j*X) if (i==2): H2=H2+BETA[I1-1]*np.cos(1.*j*X)+BETA[I2-1]*np.sin(1.*j*X) if (i==3): H3=H3+BETA[I1-1]*np.cos(1.*j*X)+BETA[I2-1]*np.sin(1.*j*X) if (i==4): H4=H4+BETA[I1-1]*np.cos(1.*j*X)+BETA[I2-1]*np.sin(1.*j*X) if (i==5): H5=H5+BETA[I1-1]*np.cos(1.*j*X)+BETA[I2-1]*np.sin(1.*j*X) if (i==6): H6=H6+BETA[I1-1]*np.cos(1.*j*X)+BETA[I2-1]*np.sin(1.*j*X) for k in range(0,K3+1): X=(T-np.min(T))/((np.max(T)-np.min(T))/2.) I=I+1 I3=I if (Check != 0): print('k',k,'Beta',BETA[I3-1]) PL=PL+BETA[I3-1]*(X**(1.*k)) G=H1+H2+H3+H4+H5+H6+PL return G,H1,H2,H3,H4,H5,H6,PL # _____________________________________________________________________________ # - Nonlinear model. # _____________________________________________________________________________ def NonLinearModel(T,BETA,Check): if (Check != 0): print('Input BETA =',BETA) K1=dcmInput(5) K2=dcmInput(6) K3=dcmInput(7) H =0. H1=0. H2=0. H3=0. H4=0. H5=0. H6=0. PL=0. ALLF=BETA[0:K1] if (Check != 0): print('ALLF =',ALLF) print('1./ALLF =',1./ALLF) I=K1 for i in range(1,K1+1): F=ALLF[i-1] X=2.*np.pi*T*F for j in range(1,K2+1): I=I+1 I1=I I=I+1 I2=I if (Check != 0): print('i',i,'j',j,'I1',I1,'I2',I2,'Beta',BETA[I1-1],BETA[I2-1]) if (i==1): H1=H1+BETA[I1-1]*np.cos(1.*j*X)+BETA[I2-1]*np.sin(1.*j*X) if (i==2): H2=H2+BETA[I1-1]*np.cos(1.*j*X)+BETA[I2-1]*np.sin(1.*j*X) if (i==3): H3=H3+BETA[I1-1]*np.cos(1.*j*X)+BETA[I2-1]*np.sin(1.*j*X) if (i==4): H4=H4+BETA[I1-1]*np.cos(1.*j*X)+BETA[I2-1]*np.sin(1.*j*X) if (i==5): H5=H5+BETA[I1-1]*np.cos(1.*j*X)+BETA[I2-1]*np.sin(1.*j*X) if (i==6): H6=H6+BETA[I1-1]*np.cos(1.*j*X)+BETA[I2-1]*np.sin(1.*j*X) for k in range(0,K3+1): X=(T-np.min(T))/((np.max(T)-np.min(T))/2.) I=I+1 I3=I if (Check != 0): print('k',k,'Beta',BETA[I3-1]) PL=PL+BETA[I3-1]*(X**(1.*k)) G=H1+H2+H3+H4+H5+H6+PL return G,H1,H2,H3,H4,H5,H6,PL # _____________________________________________________________________________ # - Linear model subroutine for optimize.leastsq in LinearLSF. # _____________________________________________________________________________ def Lfunct(BETA,T,Y,EY): G,H1,H2,H3,H4,H5,H6,PL=LinearModel(T,BETA,0.) # Check=0=No check prints return (Y-G)/EY # _____________________________________________________________________________ # - Nonlinear model subroutine for optimize.leastsq in NonLinearLSF. # _____________________________________________________________________________ def Nfunct(BETA,T,Y,EY): G,H1,H2,H3,H4,H5,H6,PL=NonLinearModel(T,BETA,0.) # Check=0=No check prints return (Y-G)/EY # _____________________________________________________________________________ # - Linear least squares fit. # _____________________________________________________________________________ def LinearLSF(T,Y,EY): K1 =dcmInput(5) K2 =dcmInput(6) K3 =dcmInput(7) TestStat=dcmInput(11) # Version 2 p=K1*2*K2+K3+1 if (TestStat == 1): EY=EY # Version 2 else: EY=np.mean(EY)+0.*EY # Version 2 # Linear model: Can use smaller ftol and xtol input values. ny=optimize.leastsq(Lfunct,np.ones(p),args=(T,Y,EY),ftol=1.0e-6,xtol=1.0e-6) BETA=ny[0] G,H1,H2,H3,H4,H5,H6,PL=LinearModel(T,BETA,0.) EPSILON=Y-G LinZ=TestStatistic(EPSILON,EY,0.) # Check=0=No check print return BETA,LinZ # _____________________________________________________________________________ # - Nonlinear least squares fit. # _____________________________________________________________________________ def NonLinearLSF(T,Y,EY,INITIALBETA,Check): TestStat=dcmInput(11) # Version 2 if (TestStat == 1): EY=EY # Version 2 else: EY=np.mean(EY)+0.*EY # Version 2 # Nonlinear model INITIALBETA is accurate: Can use larger ftol and xtol. if (Check != 0.): print('Beta_initial = ',INITIALBETA) ny=optimize.leastsq(Nfunct,INITIALBETA,args=(T,Y,EY),ftol=1.0e-4,xtol=1.0e-4) FINALBETA=ny[0] if (Check != 0.): print('Beta_final = ',FINALBETA) G,H1,H2,H3,H4,H5,H6,PL=NonLinearModel(T,FINALBETA,0.) EPSILON=Y-G NonLinZ=TestStatistic(EPSILON,EY,0.) # Check=0=No check print return FINALBETA,NonLinZ # _____________________________________________________________________________ # - Test statistic Chi^2 or R. # _____________________________________________________________________________ def TestStatistic(EPSILON,EY,Check): TestStat=dcmInput(11) if (TestStat == 1): Z=np.sqrt(np.sum((EPSILON*EPSILON)/(EY*EY))/np.size(EPSILON)) if (Check != 0.): print('z from Chi^2') else: Z=np.sqrt(np.sum((EPSILON*EPSILON))/np.size(EPSILON)) if (Check != 0.): print('z from R') return Z # _____________________________________________________________________________ # - Simulating data. # _____________________________________________________________________________ def SIMdata(Check,Tfile): Tag =dcmInput(1) SimT =dcmInput(16) SimN =dcmInput(17) SimSN =dcmInput(18) SimDT =dcmInput(19) PrintScreen=dcmInput(24) if (SimT == 1): # Simulated time points T=SimDT*np.random.uniform(0,1,SimN) # -"- k=IndexSortOrder(T) # T=T[k] # T0=T[0] # Version 2 T=T-T0 # Version 2 if (SimT != 1): # Real data time pts from file1 T0,T,ny1,ny2=DataRead(Tfile) # -"- T=T+T0 # -"- ALLF=SimulatedALLF() if (PrintScreen == 1): print('Simulated P = ',1./ALLF) BETA=LinSimulatedBeta() if (Check != 0): print('Simulated BETA = ',BETA) print("Input SN",SimSN) WriteALLF(ALLF) G,H1,H2,H3,H4,H5,H6,PL=LinearModel(T,BETA,Check) HARMONICPART=H1+H2+H3+H4+H5+H6 # Error from all h_i(t) only! ERROR=(2.**(5./2.))*np.std(HARMONICPART)/SimSN EY=ERROR*np.random.normal(0.,1.,np.size(T)) # Both +/- Y=G+EY EY=np.abs(EY) # Only + if (Check != 0): print('SN here = ',((2.**(5./2.))*np.std(Y-PL))/np.mean(EY)) print('sy here = ',np.std(Y)) print('Error here = ',np.mean(EY)) filenow=Tag+'SimulatedData.dat' code1=open(filenow,"w") for i in range(np.size(T)): code1.write("%20.10f %20.10f %20.10f \n"%(T[i],Y[i],EY[i])) code1.close() kuva=Tag+'gsim.eps' SolveBeta=0 # Simulated BETA and linear model used, because SolveBeta=0 CheckPlot(T0,T,Y,EY,BETA,kuva,SolveBeta) return BETA,filenow # _____________________________________________________________________________ # - Reading analysed simulated or real data. # _____________________________________________________________________________ def DataRead(file1): T=np.loadtxt(file1,skiprows=0,usecols=(0,)) Y=np.loadtxt(file1,skiprows=0,usecols=(1,)) E=np.loadtxt(file1,skiprows=0,usecols=(2,)) T0=np.min(T) T=T-T0 return T0,T,Y,E # _____________________________________________________________________________ # - Writing tested frequencies into file. # _____________________________________________________________________________ def WriteTestedFrequencies(j,f): ny='%5i' %(j) file='frequency'+ny+'a.dat' file=file.replace(' ','') code1=open(file,"w") for i in range(np.size(f)): code1.write("%15.15e \n" %(f[i])) code1.close() return # _____________________________________________________________________________ # - Writing tested frequencies for long search into separate files. # _____________________________________________________________________________ def LongFrequencyInterval(): dummy=dcmInput(4) K1 =dcmInput(5) nL =dcmInput(8) PMIN =dcmInput(12) PMAX =dcmInput(13) f1=dummy f2=dummy f3=dummy f4=dummy f5=dummy f6=dummy FMIN=1./PMAX FMAX=1./PMIN FSTEP=(FMAX-FMIN)/(nL+0.) # "0." = Division with integer would fail. F=FMIN+np.arange(nL+1)*FSTEP MM=np.size(F) if (K1 == 1): f1=F if (K1 == 2): f1=F[1:MM-0] f2=F[0:MM-1] if (K1 == 3): f1=F[2:MM-0] f2=F[1:MM-1] f3=F[0:MM-2] if (K1 == 4): f1=F[3:MM-0] f2=F[2:MM-1] f3=F[1:MM-2] f4=F[0:MM-3] if (K1 == 5): f1=F[4:MM-0] f2=F[3:MM-1] f3=F[2:MM-2] f4=F[1:MM-3] f5=F[0:MM-4] if (K1 == 6): f1=F[5:MM-0] f2=F[4:MM-1] f3=F[3:MM-2] f4=F[2:MM-3] f5=F[1:MM-4] f6=F[0:MM-5] if (K1 >=1): f1=f1[0:] WriteTestedFrequencies(1,f1) if (K1 >= 2): f2=f2[0:] WriteTestedFrequencies(2,f2) if (K1 >= 3): f3=f3[0:] WriteTestedFrequencies(3,f3) if (K1 >= 4): f4=f4[0:] WriteTestedFrequencies(4,f4) if (K1 >= 5): f5=f5[0:] WriteTestedFrequencies(5,f5) if (K1 >= 6): f6=f6[0:] WriteTestedFrequencies(6,f6) return # _____________________________________________________________________________ # - Writing tested frequencies for short search into separate files. # _____________________________________________________________________________ def ShortFrequencyInterval(ALLF): nS =dcmInput(9) c =dcmInput(10) PMIN =dcmInput(12) PMAX =dcmInput(13) for i in range(np.size(ALLF)): LIM=c*((1./PMIN)-(1./PMAX))/2. FMID=ALLF[i] FMIN=FMID-LIM FMAX=FMID+LIM FSTEP=(FMAX-FMIN)/(nS+0.) # nS to float, not integer! ny=np.floor((FMAX-FMIN)/FSTEP)+1 F=FMIN+(np.arange(ny))*FSTEP k=(F > 0.).nonzero() # No negative frequencies F=F[k] # close to zero ny='%5i' %(i+1) file='frequency'+ny+'b.dat' file=file.replace(' ','') code1=open(file,"w") for i in range(np.size(F)): code1.write("%0.15e \n" %(F[i])) code1.close() return # _____________________________________________________________________________ # - Reading tested frequencies from the file. # _____________________________________________________________________________ def Frequencyintervalread(file): f=np.loadtxt(file,skiprows=0,usecols=(0,)) return f # _____________________________________________________________________________ # New version 2 routine # - Pmax and Pmin can change in short search, and some plots outside original # limit may not be seen. This correct plot limits. # _____________________________________________________________________________ def NewFrequencyLimits(): K1=dcmInput(5) f1=Frequencyintervalread('frequency1a.dat') Longf=f1 # Long search frequencies if (K1 >= 1.5): f2=Frequencyintervalread('frequency2a.dat') Longf=np.append(Longf,f2) if (K1 >= 2.5): f3=Frequencyintervalread('frequency3a.dat') Longf=np.append(Longf,f3) if (K1 >= 3.5): f4=Frequencyintervalread('frequency4a.dat') Longf=np.append(Longf,f4) if (K1 >= 4.5): f5=Frequencyintervalread('frequency5a.dat') Longf=np.append(Longf,f5) if (K1 >= 5.5): f6=Frequencyintervalread('frequency6a.dat') Longf=np.append(Longf,f6) # _______________________________________________ f1=Frequencyintervalread('frequency1b.dat') Shortf=f1 # Short search frequencies if (K1 >= 1.5): f2=Frequencyintervalread('frequency2b.dat') Shortf=np.append(Shortf,f2) if (K1 >= 2.5): f3=Frequencyintervalread('frequency3b.dat') Shortf=np.append(Shortf,f3) if (K1 >= 3.5): f4=Frequencyintervalread('frequency4b.dat') Shortf=np.append(Shortf,f4) if (K1 >= 4.5): f5=Frequencyintervalread('frequency5b.dat') Shortf=np.append(Shortf,f5) if (K1 >= 5.5): f6=Frequencyintervalread('frequency6b.dat') Shortf=np.append(Shortf,f6) # _______________________________________________ Bothf=np.append(Longf,Shortf) x1=np.min(Bothf) x2=np.max(Bothf) x1=x1-0.03*(x2-x1) x2=x2+0.03*(x2-x1) return x1,x2 # _____________________________________________________________________________ # - Periodogram computation uses linear LinearLSF, because tested freqs fixed. # _____________________________________________________________________________ def Zgram(T,Y,EY,stage): dummy =dcmInput(4) K1 =dcmInput(5) if (stage < 0.5): # Long search f1=Frequencyintervalread('frequency1a.dat') if (K1 >= 1.5): f2=Frequencyintervalread('frequency2a.dat') if (K1 >= 2.5): f3=Frequencyintervalread('frequency3a.dat') if (K1 >= 3.5): f4=Frequencyintervalread('frequency4a.dat') if (K1 >= 4.5): f5=Frequencyintervalread('frequency5a.dat') if (K1 >= 5.5): f6=Frequencyintervalread('frequency6a.dat') else: # Short search f1=Frequencyintervalread('frequency1b.dat') if (K1 >= 1.5): f2=Frequencyintervalread('frequency2b.dat') if (K1 >= 2.5): f3=Frequencyintervalread('frequency3b.dat') if (K1 >= 3.5): f4=Frequencyintervalread('frequency4b.dat') if (K1 >= 4.5): f5=Frequencyintervalread('frequency5b.dat') if (K1 >= 5.5): f6=Frequencyintervalread('frequency6b.dat') F1=0.+dummy F2=0.+dummy F3=0.+dummy F4=0.+dummy F5=0.+dummy F6=0.+dummy Z =0.+dummy if (K1 == 1): for i1 in range (np.size(f1)): WriteALLF(np.vstack([f1[i1]])) BETA,Linz=LinearLSF(T,Y,EY) F1=np.append(F1,f1[i1]) Z=np.append(Z,Linz) if (K1 == 2): for i1 in range (np.size(f1)): for i2 in range (np.size(f2)): if (f2[i2] < f1[i1]): WriteALLF(np.vstack([f1[i1],f2[i2]])) BETA,Linz=LinearLSF(T,Y,EY) F1=np.append(F1,f1[i1]) F2=np.append(F2,f2[i2]) Z=np.append(Z,Linz) if (K1 == 3): for i1 in range (np.size(f1)): for i2 in range (np.size(f2)): if (f2[i2] < f1[i1]): for i3 in range(np.size(f3)): if (f3[i3] < f2[i2]): WriteALLF(np.vstack([f1[i1],f2[i2],f3[i3]])) BETA,Linz=LinearLSF(T,Y,EY) F1=np.append(F1,f1[i1]) F2=np.append(F2,f2[i2]) F3=np.append(F3,f3[i3]) Z=np.append(Z,Linz) if (K1 == 4): for i1 in range (np.size(f1)): for i2 in range (np.size(f2)): if (f2[i2] < f1[i1]): for i3 in range(np.size(f3)): if (f3[i3] < f2[i2]): for i4 in range(np.size(f4)): if (f4[i4] < f3[i3]): WriteALLF(np.vstack([f1[i1],f2[i2],f3[i3],f4[i4]])) BETA,Linz=LinearLSF(T,Y,EY) F1=np.append(F1,f1[i1]) F2=np.append(F2,f2[i2]) F3=np.append(F3,f3[i3]) F4=np.append(F4,f4[i4]) Z=np.append(Z,Linz) if (K1 == 5): for i1 in range (np.size(f1)): for i2 in range (np.size(f2)): if (f2[i2] < f1[i1]): for i3 in range(np.size(f3)): if (f3[i3] < f2[i2]): for i4 in range(np.size(f4)): if (f4[i4] < f3[i3]): for i5 in range(np.size(f5)): if (f5[i5] < f4[i4]): WriteALLF(np.vstack([f1[i1],f2[i2],f3[i3],f4[i4],f5[i5]])) BETA,Linz=LinearLSF(T,Y,EY) F1=np.append(F1,f1[i1]) F2=np.append(F2,f2[i2]) F3=np.append(F3,f3[i3]) F4=np.append(F4,f4[i4]) F5=np.append(F5,f5[i5]) Z=np.append(Z,Linz) if (K1 == 6): for i1 in range (np.size(f1)): for i2 in range (np.size(f2)): if (f2[i2] < f1[i1]): for i3 in range(np.size(f3)): if (f3[i3] < f2[i2]): for i4 in range(np.size(f4)): if (f4[i4] < f3[i3]): for i5 in range(np.size(f5)): if (f5[i5] < f4[i4]): for i6 in range(np.size(f6)): if (f6[i6] < f5[i5]): WriteALLF(np.vstack([f1[i1],f2[i2],f3[i3],f4[i4],f5[i5],f6[i6]])) BETA,Linz=LinearLSF(T,Y,EY) F1=np.append(F1,f1[i1]) F2=np.append(F2,f2[i2]) F3=np.append(F3,f3[i3]) F4=np.append(F4,f4[i4]) F5=np.append(F5,f5[i5]) F6=np.append(F6,f6[i6]) Z=np.append(Z,Linz) if (np.size(F1) != 1): F1=F1[1:] if (np.size(F2) != 1): F2=F2[1:] if (np.size(F3) != 1): F3=F3[1:] if (np.size(F4) != 1): F4=F4[1:] if (np.size(F5) != 1): F5=F5[1:] if (np.size(F6) != 1): F6=F6[1:] Z =Z[1:] # Z[0]=dummy for all K1 values # Best frequencies ___________________________ j=(Z == np.min(Z)).nonzero() if (K1 == 1): ALLF=F1[j] if (K1 == 2): ALLF=np.vstack([F1[j],F2[j]]) if (K1 == 3): ALLF=np.vstack([F1[j],F2[j],F3[j]]) if (K1 == 4): ALLF=np.vstack([F1[j],F2[j],F3[j],F4[j]]) if (K1 == 5): ALLF=np.vstack([F1[j],F2[j],F3[j],F4[j],F5[j]]) if (K1 == 6): ALLF=np.vstack([F1[j],F2[j],F3[j],F4[j],F5[j],F6[j]]) return F1,F2,F3,F4,F5,F6,Z,ALLF # _____________________________________________________________________________ # - Print detected periods to screen. # _____________________________________________________________________________ def PCheck(F1,F2,F3,F4,F5,F6,Z,txt): K1 =dcmInput(5) PrintScreen=dcmInput(24) j=(Z == np.min(Z)).nonzero() ZMIN=Z[j[0]] if (PrintScreen == 1): if (K1 == 1): print(txt,'Detected P=',1./F1[j]) if (K1 == 2): print(txt,'Detected P=',1./F1[j],1./F2[j]) if (K1 == 3): print(txt,'Detected P=',1./F1[j],1./F2[j],1./F3[j]) if (K1 == 4): print(txt,'Detected P=',1./F1[j],1./F2[j],1./F3[j],1./F4[j]) if (K1 == 5): print(txt,'Detected P=',1./F1[j],1./F2[j],1./F3[j],1./F4[j],1./F5[j]) if (K1 == 6): print(txt,'Detected P=',1./F1[j],1./F2[j],1./F3[j],1./F4[j],1./F5[j],1./F6[j]) print(txt,'ZMIN=',ZMIN) return # _____________________________________________________________________________ # - Storing periodograms into the file. # _____________________________________________________________________________ def WritePeriodogram(F1,F2,F3,F4,F5,F6,Z,OUTfile): dummy=dcmInput(4) K1 =dcmInput(5) Zcode=open(OUTfile,"w") for i in range(np.size(Z)): if (K1 == 1): Zcode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e\n" \ %(F1[i],dummy,dummy,dummy,dummy,dummy,Z[i])) if (K1 == 2): Zcode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e\n" \ %(F1[i],F2[i],dummy,dummy,dummy,dummy,Z[i])) if (K1 == 3): Zcode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e\n" \ %(F1[i],F2[i],F3[i],dummy,dummy,dummy,Z[i])) if (K1 == 4): Zcode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e\n" \ %(F1[i],F2[i],F3[i],F4[i],dummy,dummy,Z[i])) if (K1 == 5): Zcode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e\n" \ %(F1[i],F2[i],F3[i],F4[i],F5[i],dummy,Z[i])) if (K1 == 6): Zcode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e\n" \ %(F1[i],F2[i],F3[i],F4[i],F5[i],F6[i],Z[i])) Zcode.close() return # _____________________________________________________________________________ # - Reading periodograms from the file. # _____________________________________________________________________________ def ReadPeriodogram(file1): dummy=dcmInput(4) F1=dummy F2=dummy F3=dummy F4=dummy F5=dummy F6=dummy Z =dummy file = open(file1, 'r') rivi = file.readline() while (len(rivi) > 0): rivi=rivi.rstrip() osat=rivi.split() F1=np.append(F1,np.array(float(osat[0]))) F2=np.append(F2,np.array(float(osat[1]))) F3=np.append(F3,np.array(float(osat[2]))) F4=np.append(F4,np.array(float(osat[3]))) F5=np.append(F5,np.array(float(osat[4]))) F6=np.append(F6,np.array(float(osat[5]))) Z =np.append( Z,np.array(float(osat[6]))) rivi = file.readline() F1=F1[1:] F2=F2[1:] F3=F3[1:] F4=F4[1:] F5=F5[1:] F6=F6[1:] Z = Z[1:] return F1,F2,F3,F4,F5,F6,Z # _____________________________________________________________________________ # - Separating periodogram slices at global minimum. # _____________________________________________________________________________ def SeparatePeriodograms(F1,F2,F3,F4,F5,F6,Z): dummy=dcmInput(4) K1 =dcmInput(5) FF1 =dummy FF2 =dummy FF3 =dummy FF4 =dummy FF5 =dummy FF6 =dummy ZZ1 =dummy ZZ2 =dummy ZZ3 =dummy ZZ4 =dummy ZZ5 =dummy ZZ6 =dummy j=(Z==np.min(Z)).nonzero() if (K1==1): FF1=F1 ZZ1=Z if (K1==2): k=((F2==F2[j])).nonzero() FF1=F1[k] ZZ1=Z[k] k=((F1==F1[j])).nonzero() FF2=F2[k] ZZ2=Z[k] if (K1==3): k=((F2==F2[j]) & (F3==F3[j])).nonzero() FF1=F1[k] ZZ1=Z[k] k=((F1==F1[j]) & (F3==F3[j])).nonzero() FF2=F2[k] ZZ2=Z[k] k=((F1==F1[j]) & (F2==F2[j])).nonzero() FF3=F3[k] ZZ3=Z[k] if (K1==4): k=((F2==F2[j]) & (F3==F3[j]) & (F4==F4[j])).nonzero() FF1=F1[k] ZZ1=Z[k] k=((F1==F1[j]) & (F3==F3[j]) & (F4==F4[j])).nonzero() FF2=F2[k] ZZ2=Z[k] k=((F1==F1[j]) & (F2==F2[j]) & (F4==F4[j])).nonzero() FF3=F3[k] ZZ3=Z[k] k=((F1==F1[j]) & (F2==F2[j]) & (F3==F3[j])).nonzero() FF4=F4[k] ZZ4=Z[k] if (K1==5): k=((F2==F2[j]) & (F3==F3[j]) & (F4==F4[j]) & (F5==F5[j])).nonzero() FF1=F1[k] ZZ1=Z[k] k=((F1==F1[j]) & (F3==F3[j]) & (F4==F4[j]) & (F5==F5[j])).nonzero() FF2=F2[k] ZZ2=Z[k] k=((F1==F1[j]) & (F2==F2[j]) & (F4==F4[j]) & (F5==F5[j])).nonzero() FF3=F3[k] ZZ3=Z[k] k=((F1==F1[j]) & (F2==F2[j]) & (F3==F3[j]) & (F5==F5[j])).nonzero() FF4=F4[k] ZZ4=Z[k] k=((F1==F1[j]) & (F2==F2[j]) & (F3==F3[j]) & (F4==F4[j])).nonzero() FF5=F5[k] ZZ5=Z[k] if (K1==6): k=((F2==F2[j]) & (F3==F3[j]) & (F4==F4[j]) & (F5==F5[j]) & (F6==F6[j])).nonzero() FF1=F1[k] ZZ1=Z[k] k=((F1==F1[j]) & (F3==F3[j]) & (F4==F4[j]) & (F5==F5[j]) & (F6==F6[j])).nonzero() FF2=F2[k] ZZ2=Z[k] k=((F1==F1[j]) & (F2==F2[j]) & (F4==F4[j]) & (F5==F5[j]) & (F6==F6[j])).nonzero() FF3=F3[k] ZZ3=Z[k] k=((F1==F1[j]) & (F2==F2[j]) & (F3==F3[j]) & (F5==F5[j]) & (F6==F6[j])).nonzero() FF4=F4[k] ZZ4=Z[k] k=((F1==F1[j]) & (F2==F2[j]) & (F3==F3[j]) & (F4==F4[j]) & (F6==F6[j])).nonzero() FF5=F5[k] ZZ5=Z[k] k=((F1==F1[j]) & (F2==F2[j]) & (F3==F3[j]) & (F4==F4[j]) & (F5==F5[j])).nonzero() FF6=F6[k] ZZ6=Z[k] return FF1,FF2,FF3,FF4,FF5,FF6,ZZ1,ZZ2,ZZ3,ZZ4,ZZ5,ZZ6 # _____________________________________________________________________________ # - Plotting separate periodograms: red, blue,green, yellow, magenta, and cyan. # _____________________________________________________________________________ def PlotSeparate(FF1,FF2,FF3,FF4,FF5,FF6,ZZ1,ZZ2,ZZ3,ZZ4,ZZ5,ZZ6,txt): dummy=dcmInput(4) koko=3 x1,x2=NewFrequencyLimits() # Version 2 pl.xlim([x1,x2]) if (np.min(FF1) != dummy): pl.plot(FF1,ZZ1,'or',ms=koko) pl.plot(FF1,ZZ1,'r',ms=koko) MarkMinima(FF1,ZZ1,1,txt) # Version 2 if (np.min(FF2) != dummy): pl.plot(FF2,ZZ2,'ob',ms=koko) pl.plot(FF2,ZZ2,'b',ms=koko) MarkMinima(FF2,ZZ2,2,txt) # Version 2 if (np.min(FF3) != dummy): pl.plot(FF3,ZZ3,'og',ms=koko) pl.plot(FF3,ZZ3,'g',ms=koko) MarkMinima(FF3,ZZ3,3,txt) # Version 2 if (np.min(FF4) != dummy): pl.plot(FF4,ZZ4,'oy',ms=koko) pl.plot(FF4,ZZ4,'y',ms=koko) MarkMinima(FF4,ZZ4,4,txt) # Version 2 if (np.min(FF5) != dummy): pl.plot(FF5,ZZ5,'om',ms=koko) pl.plot(FF5,ZZ5,'m',ms=koko) MarkMinima(FF5,ZZ5,5,txt) # Version 2 if (np.min(FF6) != dummy): pl.plot(FF6,ZZ6,'oc',ms=koko) pl.plot(FF6,ZZ6,'c',ms=koko) MarkMinima(FF6,ZZ6,6,txt) # Version 2 return # _____________________________________________________________________________ # - New version 2: Marks periodogram minima. # _____________________________________________________________________________ def MarkMinima(X,Y,m,txt): PrintScreen=dcmInput(24) j=(Y == np.min(Y)).nonzero() pl.plot(X[j],Y[j],'Dk',ms=10,markerfacecolor="None") if ((X[j] == np.min(X)) | (X[j] == np.max(X))): ny='%1i' %(m) txt=txt+ny+' period at edge' if (PrintScreen == 1.): print(txt) txt=txt+'\n' Notes(txt) # Version 2 return # _____________________________________________________________________________ # - Plotting all periodograms. # _____________________________________________________________________________ def PlotPeriodograms(file1,file2): Tag =dcmInput(1) F1,F2,F3,F4,F5,F6,Z=ReadPeriodogram(file1) # Long FF1,FF2,FF3,FF4,FF5,FF6,ZZ1,ZZ2,ZZ3,ZZ4,ZZ5,ZZ6\ =SeparatePeriodograms(F1,F2,F3,F4,F5,F6,Z) pl.axes([0.17,0.56,0.80,0.37]) txt='Long search: ' PlotSeparate(FF1,FF2,FF3,FF4,FF5,FF6,ZZ1,ZZ2,ZZ3,ZZ4,ZZ5,ZZ6,txt) pl.tick_params(labelbottom='off') pl.title('Long search',fontsize=14) pl.ylabel('$z_1(f_1), ..., z_{K_1}(f_{K_1})$',fontsize=14) pl.xticks(fontsize=10) pl.yticks(fontsize=10) # _________________________________________________________________________ F1,F2,F3,F4,F5,F6,Z=ReadPeriodogram(file2) # Short FF1,FF2,FF3,FF4,FF5,FF6,ZZ1,ZZ2,ZZ3,ZZ4,ZZ5,ZZ6\ =SeparatePeriodograms(F1,F2,F3,F4,F5,F6,Z) pl.axes([0.17,0.12,0.80,0.37]) txt='Short search: ' PlotSeparate(FF1,FF2,FF3,FF4,FF5,FF6,ZZ1,ZZ2,ZZ3,ZZ4,ZZ5,ZZ6,txt) pl.title('Short search',fontsize=16) pl.ylabel('$z_1(f_1), ..., z_{K_1}(f_{K_1})$',fontsize=14) pl.xlabel('$f_1, ..., f_{K_1}$',fontsize=14) pl.xticks(fontsize=10) pl.yticks(fontsize=10) kuva=Tag+'z.eps' pl.savefig(kuva) pl.clf() return # _____________________________________________________________________________ # - Plotting modelling results: red, blue, green, yellow, magenta, cyan. # _____________________________________________________________________________ def CheckPlot(T0,T,Y,EY,BETAFIXED,kuva,SolveBeta): NonLinear=dcmInput(15) TT=np.min(T)+(np.max(T)-np.min(T))*np.arange(10001)/10000. if (SolveBeta == 0): # Simulated BETAFIXED is given, not solved! G,H1,H2,H3,H4,H5,H6,PL=LinearModel(T,BETAFIXED,0.) EPSILON=Y-G GG,HH1,HH2,HH3,HH4,HH5,HH6,PPL=LinearModel(TT,BETAFIXED,0.) if (SolveBeta != 0): # Linear and/or nonlinear BETA is solved. LINBETA,LinZ=LinearLSF(T,Y,EY) if (NonLinear !=1): G,H1,H2,H3,H4,H5,H6,PL=LinearModel(T,LINBETA,0.) EPSILON=Y-G GG,HH1,HH2,HH3,HH4,HH5,HH6,PPL=LinearModel(TT,LINBETA,0.) if (NonLinear ==1): ALLF=ReadALLF() # This comes from short search. NONLINBETA=np.hstack((ALLF,LINBETA)) FINALBETA,NonLinZ=NonLinearLSF(T,Y,EY,NONLINBETA,0.) G,H1,H2,H3,H4,H5,H6,PL=NonLinearModel(T,FINALBETA,0.) EPSILON=Y-G GG,HH1,HH2,HH3,HH4,HH5,HH6,PPL=NonLinearModel(TT,FINALBETA,0.) koko=2.0 pl.axes([0.15,0.56,0.82,0.42]) x1=np.min(T0+T) x2=np.max(T0+T) pl.xlim([x1,x2]) y1=np.min([np.min(Y-EY),np.min(GG)]) y2=np.max([np.max(Y+EY),np.max(GG)]) pl.ylim([y1,y2]) pl.errorbar(T0+T,Y,EY,fmt='ok',ms=koko,\ capsize=0,ecolor='k',markeredgecolor='k') pl.plot(T0+TT,PPL,':k') pl.plot(T0+TT,GG,'k') pl.xticks(fontsize=10) pl.yticks(fontsize=10) pl.ylabel('$y$',fontsize=14) pl.text(x1+0.95*(x2-x1),y1+0.90*(y2-y1),'a',fontsize=16) pl.tick_params(labelbottom='off') pl.axes([0.15,0.11,0.82,0.42]) Dy=np.mean(GG-PPL)-3.0*np.std(GG-PPL)-3.0*np.std(EPSILON) # Version 2 y1=np.min(Dy-1.10*(EPSILON+EY)) # Version 2 y2=np.max([np.max(Y+2*EY-PL),np.max(GG-PPL)]) pl.ylim([y1,y2]) x1=np.min(T0+T) x2=np.max(T0+T) pl.xlim([x1,x2]) pl.plot(T0+TT,GG-PPL,'k',zorder=9) pl.plot(T0+TT,HH1,'r',zorder=3) if (np.size(HH2) > 1): pl.plot(T0+TT,HH2,'b',zorder=4) if (np.size(HH3) > 1): pl.plot(T0+TT,HH3,'g',zorder=5) if (np.size(HH4) > 1): pl.plot(T0+TT,HH4,'y',zorder=6) if (np.size(HH5) > 1): pl.plot(T0+TT,HH5,'m',zorder=7) if (np.size(HH6) > 1): pl.plot(T0+TT,HH6,'c',zorder=8) pl.errorbar(T0+T,Y-PL,EY,fmt='ok',ms=koko,\ capsize=0,ecolor='k',markeredgecolor='k',zorder=2) pl.plot(T0+T,0.*T+Dy,':b') pl.errorbar(T0+T,Dy+EPSILON,EY,fmt='ob',ms=koko,\ capsize=0,ecolor='b',markeredgecolor='b',zorder=1) pl.xticks(fontsize=10) pl.yticks(fontsize=10) pl.xlabel('t',fontsize=14) pl.ylabel('$y$',fontsize=14) pl.text(x1+0.95*(x2-x1),y1+0.90*(y2-y1),'b',fontsize=16) pl.savefig(kuva) pl.clf() return # _____________________________________________________________________________ # - Storing modelling results T0, DELTAT and BETA into OUTfile. # _____________________________________________________________________________ def WriteBeta(T0,T,ALLF,BETA,OUTfile,Style,Check): dummy=dcmInput(4) K1 =dcmInput(5) K2 =dcmInput(6) K3 =dcmInput(7) FF=np.zeros(6)+dummy # 6 BB=np.zeros(24)+dummy # 6 x 4 MM=np.zeros(7)+dummy # 7 if (Check != 0): print('K1,K2,K3=',K1,K2,K3) print('ALLF = ',ALLF) print('Beta =',BETA) I=0 J=0 for i in range(1,K1+1): FF[i-1]=ALLF[i-1] for j in range(1,K2+1): I=I+1 I1=I I=I+1 I2=I if (i==1): if (Check != 0): print('i',i,'j',j,'I1',I1,'I2',I2,'Beta',\ BETA[I1-1],BETA[I2-1]) BB[J]=BETA[I1-1] ; J=J+1 BB[J]=BETA[I2-1] ; J=J+1 if (i==2): if (Check != 0): print('i',i,'j',j,'I1',I1,'I2',I2,'Beta',\ BETA[I1-1],BETA[I2-1]) BB[J]=BETA[I1-1] ; J=J+1 BB[J]=BETA[I2-1] ; J=J+1 if (i==3): if (Check != 0): print('i',i,'j',j,'I1',I1,'I2',I2,'Beta',\ BETA[I1-1],BETA[I2-1]) BB[J]=BETA[I1-1] ; J=J+1 BB[J]=BETA[I2-1] ; J=J+1 if (i==4): if (Check != 0): print('i',i,'j',j,'I1',I1,'I2',I2,'Beta',\ BETA[I1-1],BETA[I2-1]) BB[J]=BETA[I1-1] ; J=J+1 BB[J]=BETA[I2-1] ; J=J+1 if (i==5): if (Check != 0): print('i',i,'j',j,'I1',I1,'I2',I2,'Beta',\ BETA[I1-1],BETA[I2-1]) BB[J]=BETA[I1-1] ; J=J+1 BB[J]=BETA[I2-1] ; J=J+1 if (i==6): if (Check != 0): print('i',i,'j',j,'I1',I1,'I2',I2,'Beta',\ BETA[I1-1],BETA[I2-1]) BB[J]=BETA[I1-1] ; J=J+1 BB[J]=BETA[I2-1] ; J=J+1 if (K2 == 1): J=J+2 for k in range(0,K3+1): I=I+1 I3=I if (Check != 0): print('k',k,'Beta',BETA[I3-1]) MM[k]=BETA[I3-1] if (Check != 0): print('FF=',FF) print('BB=',BB) print('MM=',MM) AllBeta=np.hstack((FF,BB,MM)) col='%0.20e %0.20e '%(T0,np.max(T)-np.min(T)) for i in range(np.size(AllBeta)): col=col+'%0.8e '%(AllBeta[i]) col=col+'\n' # This string contains all results. if (Check != 0): print(col) if (Style == 1): Resultcode=open(OUTfile,"w") # Original model results written first. if (Style == 2): Resultcode=open(OUTfile,"a") # Bootstrap results appended later. Resultcode.write(col) Resultcode.close() return # ______________________________________________________________________________ # - New version 2: Error notes. # ______________________________________________________________________________ def NotesLineOne(): tag=dcmInput(1) file1=tag+'Notes.dat' ErrorNote=open(file1,"w") ErrorNote.write('') ErrorNote.close() return # ______________________________________________________________________________ # - New version 2: Error notes. # ______________________________________________________________________________ def Notes(ny): tag=dcmInput(1) file1=tag+'Notes.dat' ErrorNote=open(file1,"a") ErrorNote.write(ny) ErrorNote.close() return # _____________________________________________________________________________ # - Modelling original data. # _____________________________________________________________________________ def dcmModel(T0,T,Y,EY,BetaFile,Style): K1 =dcmInput(5) NonLinear =dcmInput(15) PrintScreen=dcmInput(24) #Style=1 # "w" in WriteBeta(...): Writes to OUTfile if (NonLinear == 1): LINBETA,LinZ=LinearLSF(T,Y,EY) ALLF=ReadALLF() # Value stored after short search NONLINBETA=np.hstack((ALLF,LINBETA)) FINALBETA,NonLinZ=NonLinearLSF(T,Y,EY,NONLINBETA,0.) if (PrintScreen == 1): print('Final ZMIN =',NonLinZ) BETA_I=FINALBETA[0:K1] BETA_II=FINALBETA[K1:] WriteBeta(T0,T,BETA_I,BETA_II,BetaFile,Style,0.) OrigG,H1,H2,H3,H4,H5,H6,PL=NonLinearModel(T,FINALBETA,0.) if (NonLinear != 1): LINBETA,LinZ=LinearLSF(T,Y,EY) ALLF=ReadALLF() # Value stored after short search! if (PrintScreen ==1): print('Final ZMIN =',LinZ) WriteBeta(T0,T,ALLF,LINBETA,BetaFile,Style,0.) OrigG,H1,H2,H3,H4,H5,H6,PL=LinearModel(T,LINBETA,0.) FINALBETA=np.hstack((ALLF,LINBETA)) return OrigG,FINALBETA # _____________________________________________________________________________ # Version 2: Signal plots, Points file and Curves file. # _____________________________________________________________________________ def PhasePlots(T0,T,Y,EY): Tag =dcmInput(1) dummy =dcmInput(4) K1 =dcmInput(5) NonLinear =dcmInput(15) TT=np.min(T)+(np.max(T)-np.min(T))*(np.arange(5001)/5000.) if (NonLinear == 1): LINBETA,LinZ=LinearLSF(T,Y,EY) ALLF=ReadALLF() # Value stored after short search NONLINBETA=np.hstack((ALLF,LINBETA)) FINALBETA,NonLinZ=NonLinearLSF(T,Y,EY,NONLINBETA,0.) BETA_I=FINALBETA[0:K1] BETA_II=FINALBETA[K1:] OrigG,H1,H2,H3,H4,H5,H6,PL=NonLinearModel(T,FINALBETA,0.) OrigGG,HH1,HH2,HH3,HH4,HH5,HH6,PPL=NonLinearModel(TT,FINALBETA,0.) if (NonLinear != 1): LINBETA,LinZ=LinearLSF(T,Y,EY) ALLF=ReadALLF() # Value stored after short search! OrigG,H1,H2,H3,H4,H5,H6,PL=LinearModel(T,LINBETA,0.) OrigGG,HH1,HH2,HH3,HH4,HH5,HH6,PPL=LinearModel(TT,LINBETA,0.) FINALBETA=np.hstack((ALLF,LINBETA)) # Plots begin ________________________________________________________ x1=np.min(T+T0)-0.0*(np.max(T+T0)-np.min(T+T0)) x2=np.max(T+T0)+0.0*(np.max(T+T0)-np.min(T+T0)) if (K1 >= 1): pl.axes([0.15,0.84,0.595,0.128]) pl.plot(TT+T0,HH1,'r') pl.errorbar(T0+T,Y-(OrigG-H1),EY,fmt='ok',ms=1) # suppress scientific notations ... pl.ticklabel_format(useOffset=False, style='plain') pl.xlim([x1,x2]) pl.xticks(fontsize=7) pl.yticks(fontsize=7) if (K1 > 1): pl.tick_params(labelbottom='off') if (K1==1): pl.xlabel('$t$',fontsize=12) pl.ylabel('$h_1(t)$',fontsize=12) phi=T*FINALBETA[0] phi=phi-np.floor(phi) PHI=TT*FINALBETA[0] PHI=PHI-np.floor(PHI) YY=HH1 k=IndexSortOrder(PHI) PHI=PHI[k] YY=YY[k] pl.axes([0.78,0.84,0.195,0.128]) pl.errorbar(phi,Y-(OrigG-H1),EY,fmt='ok',ms=1,zorder=1) pl.plot(PHI,YY,'r',linestyle='-',zorder=2) pl.xlim([0.,1.]) pl.xticks(fontsize=7) pl.yticks(fontsize=7) pl.tick_params(labelleft='off') if (K1 > 1): pl.tick_params(labelbottom='off') if (K1==1): pl.xlabel('$\phi$',fontsize=12) if (K1 >= 2): pl.axes([0.15,0.69,0.595,0.128]) pl.plot(TT+T0,HH2,'b') pl.errorbar(T0+T,Y-(OrigG-H2),EY,fmt='ok',ms=1) # suppress scientific notations ... pl.ticklabel_format(useOffset=False, style='plain') pl.xlim([x1,x2]) pl.xticks(fontsize=7) pl.yticks(fontsize=7) if (K1 > 2): pl.tick_params(labelbottom='off') if (K1==2): pl.xlabel('$t$',fontsize=12) pl.ylabel('$h_2(t)$',fontsize=12) phi=T*FINALBETA[1] phi=phi-np.floor(phi) PHI=TT*FINALBETA[1] PHI=PHI-np.floor(PHI) YY=HH2 k=IndexSortOrder(PHI) PHI=PHI[k] YY=YY[k] pl.axes([0.78,0.69,0.195,0.128]) pl.errorbar(phi,Y-(OrigG-H2),EY,fmt='ok',ms=1,zorder=1) pl.plot(PHI,YY,'b',zorder=2) pl.xlim([0.,1.]) pl.xticks(fontsize=7) pl.yticks(fontsize=7) pl.tick_params(labelleft='off') if (K1 > 2): pl.tick_params(labelbottom='off') if (K1==2): pl.xlabel('$\phi$',fontsize=12) if (K1 >= 3): pl.axes([0.15,0.54,0.595,0.128]) pl.plot(TT+T0,HH3,'g') pl.errorbar(T0+T,Y-(OrigG-H3),EY,fmt='ok',ms=1) # suppress scientific notations ... pl.ticklabel_format(useOffset=False, style='plain') pl.xlim([x1,x2]) pl.xticks(fontsize=7) pl.yticks(fontsize=7) if (K1 > 3): pl.tick_params(labelbottom='off') if (K1==3): pl.xlabel('$t$',fontsize=12) pl.ylabel('$h_3(t)$',fontsize=12) phi=T*FINALBETA[2] phi=phi-np.floor(phi) PHI=TT*FINALBETA[2] PHI=PHI-np.floor(PHI) YY=HH3 k=IndexSortOrder(PHI) PHI=PHI[k] YY=YY[k] pl.axes([0.78,0.54,0.195,0.128]) pl.errorbar(phi,Y-(OrigG-H3),EY,fmt='ok',ms=1,zorder=1) pl.plot(PHI,YY,'g',zorder=2) pl.xlim([0.,1.]) pl.xticks(fontsize=7) pl.yticks(fontsize=7) pl.tick_params(labelleft='off') if (K1 > 3): pl.tick_params(labelbottom='off') if (K1==3): pl.xlabel('$\phi$',fontsize=12) if (K1 >= 4): pl.axes([0.15,0.39,0.595,0.128]) pl.plot(TT+T0,HH4,'y') pl.errorbar(T0+T,Y-(OrigG-H4),EY,fmt='ok',ms=1) # suppress scientific notations ... pl.ticklabel_format(useOffset=False, style='plain') pl.xlim([x1,x2]) pl.xticks(fontsize=7) pl.yticks(fontsize=7) if (K1 > 4): pl.tick_params(labelbottom='off') if (K1==4): pl.xlabel('$t$',fontsize=12) pl.ylabel('$h_4(t)$',fontsize=12) phi=T*FINALBETA[3] phi=phi-np.floor(phi) PHI=TT*FINALBETA[3] PHI=PHI-np.floor(PHI) YY=HH4 k=IndexSortOrder(PHI) PHI=PHI[k] YY=YY[k] pl.axes([0.78,0.39,0.195,0.128]) pl.errorbar(phi,Y-(OrigG-H4),EY,fmt='ok',ms=1,zorder=1) pl.plot(PHI,YY,'y',zorder=2) pl.xlim([0.,1.]) pl.xticks(fontsize=7) pl.yticks(fontsize=7) pl.tick_params(labelleft='off') if (K1 > 4): pl.tick_params(labelbottom='off') if (K1==4): pl.xlabel('$\phi$',fontsize=12) if (K1 >= 5): pl.axes([0.15,0.24,0.595,0.128]) pl.plot(TT+T0,HH5,'m') pl.errorbar(T0+T,Y-(OrigG-H5),EY,fmt='ok',ms=1) # suppress scientific notations ... pl.ticklabel_format(useOffset=False, style='plain') pl.xlim([x1,x2]) pl.xticks(fontsize=7) pl.yticks(fontsize=7) if (K1 > 5): pl.tick_params(labelbottom='off') if (K1==5): pl.xlabel('$t$',fontsize=12) pl.ylabel('$h_5(t)$',fontsize=12) phi=T*FINALBETA[4] phi=phi-np.floor(phi) PHI=TT*FINALBETA[4] PHI=PHI-np.floor(PHI) YY=HH5 k=IndexSortOrder(PHI) PHI=PHI[k] YY=YY[k] pl.axes([0.78,0.24,0.195,0.128]) pl.errorbar(phi,Y-(OrigG-H5),EY,fmt='ok',ms=1,zorder=1) pl.plot(PHI,YY,'m',zorder=2) pl.xlim([0.,1.]) pl.xticks(fontsize=7) pl.yticks(fontsize=7) pl.tick_params(labelleft='off') if (K1 > 5): pl.tick_params(labelbottom='off') if (K1==5): pl.xlabel('$\phi$',fontsize=12) if (K1 >= 6): pl.axes([0.15,0.09,0.595,0.128]) pl.plot(TT+T0,HH6,'c') pl.errorbar(T0+T,Y-(OrigG-H6),EY,fmt='ok',ms=1) # suppress scientific notations ... pl.ticklabel_format(useOffset=False, style='plain') pl.xlim([x1,x2]) pl.xticks(fontsize=7) pl.yticks(fontsize=7) pl.tick_params(labelbottom='on') pl.xlabel('$t$',fontsize=12) pl.ylabel('$h_6(t)$',fontsize=12) phi=T*FINALBETA[5] phi=phi-np.floor(phi) PHI=TT*FINALBETA[5] PHI=PHI-np.floor(PHI) YY=HH6 k=IndexSortOrder(PHI) PHI=PHI[k] YY=YY[k] pl.axes([0.78,0.09,0.195,0.128]) pl.errorbar(phi,Y-(OrigG-H6),EY,fmt='ok',ms=1,zorder=1) pl.plot(PHI,YY,'c',zorder=2) pl.xlim([0.,1.]) pl.xticks(fontsize=7) pl.yticks(fontsize=7) pl.tick_params(labelleft='off') pl.xlabel('$\phi$',fontsize=12) Figure=Tag+'Signals.eps' pl.savefig(Figure) pl.clf() OUTfile=Tag+'Points.dat' Ecode=open(OUTfile,"w") for i in range(np.size(T)): if (K1==1): Ecode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e \n" %(T[i]+T0,Y[i],EY[i],OrigG[i],H1[i],dummy,dummy,dummy,dummy,dummy,PL[i])) if (K1==2): Ecode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e \n" %(T[i]+T0,Y[i],EY[i],OrigG[i],H1[i],H2[i],dummy,dummy,dummy,dummy,PL[i])) if (K1==3): Ecode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e \n" %(T[i]+T0,Y[i],EY[i],OrigG[i],H1[i],H2[i],H3[i],dummy,dummy,dummy,PL[i])) if (K1==4): Ecode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e \n" %(T[i]+T0,Y[i],EY[i],OrigG[i],H1[i],H2[i],H3[i],H4[i],dummy,dummy,PL[i])) if (K1==5): Ecode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e \n" %(T[i]+T0,Y[i],EY[i],OrigG[i],H1[i],H2[i],H3[i],H4[i],H5[i],dummy,PL[i])) if (K1==6): Ecode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e \n" %(T[i]+T0,Y[i],EY[i],OrigG[i],H1[i],H2[i],H3[i],H4[i],H5[i],H6[i],PL[i])) Ecode.close() OUTfile=Tag+'Curves.dat' Ecode=open(OUTfile,"w") for i in range(np.size(TT)): if (K1==1): Ecode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e \n" %(TT[i]+T0,OrigGG[i],HH1[i],dummy,dummy,dummy,dummy,dummy,PPL[i])) if (K1==2): Ecode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e \n" %(TT[i]+T0,OrigGG[i],HH1[i],HH2[i],dummy,dummy,dummy,dummy,PPL[i])) if (K1==3): Ecode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e \n" %(TT[i]+T0,OrigGG[i],HH1[i],HH2[i],HH3[i],dummy,dummy,dummy,PPL[i])) if (K1==4): Ecode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e \n" %(TT[i]+T0,OrigGG[i],HH1[i],HH2[i],HH3[i],HH4[i],dummy,dummy,PPL[i])) if (K1==5): Ecode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e \n" %(TT[i]+T0,OrigGG[i],HH1[i],HH2[i],HH3[i],HH4[i],HH5[i],dummy,PPL[i])) if (K1==6): Ecode.write("%0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e %0.10e \n" %(TT[i]+T0,OrigGG[i],HH1[i],HH2[i],HH3[i],HH4[i],HH5[i],HH6[i],PPL[i])) Ecode.close() return # _____________________________________________________________________________ # - Bootstrap estimates for free parameters. # _____________________________________________________________________________ def BootStrap(T0,T,Y,EY,OrigG,OUTfile,Check): K1 =dcmInput(5) Rounds =dcmInput(14) NonLinear =dcmInput(15) PrintScreen=dcmInput(24) OrigEPSILON=Y-OrigG OrigE=EY for i in range(0,Rounds): if (PrintScreen == 1): print('Bootsrap round ',i+1,' out of ',Rounds) l=np.floor(np.size(T)*np.random.uniform(0,1,np.size(T))) l=l.astype(int) StarEPSILON=OrigEPSILON[l] StarY=OrigG+StarEPSILON StarEY=OrigE[l] stage=1. # Short frequency intervals F1,F2,F3,F4,F5,F6,Z,BootALLF=Zgram(T,StarY,StarEY,stage) WriteALLF(BootALLF) Style=2. # "a" in WriteBeta(...): Appends to OUTfile if (NonLinear == 1): LINBETA,LinZ=LinearLSF(T,StarY,StarEY) ALLF=ReadALLF() NONLINBETA=np.hstack((ALLF,LINBETA)) FINALBETA,NonLinZ=NonLinearLSF(T,StarY,StarEY,NONLINBETA,0.) BETA_I=FINALBETA[0:K1] BETA_II=FINALBETA[K1:] WriteBeta(T0,T,BETA_I,BETA_II,OUTfile,Style,Check) else: LINBETA,LinZ=LinearLSF(T,StarY,StarEY) ALLF=ReadALLF() WriteBeta(T0,T,ALLF,LINBETA,OUTfile,Style,Check) return # _____________________________________________________________________________ # - Analysing results for the chosen harmonic h_i(t), where VAL=i. # _____________________________________________________________________________ def AnalyseHarmonic(Val,INfile,Check): dummy =dcmInput(4) K1 =dcmInput(5) HH=dummy i=0 # Line counter AnalyseCode=open(INfile,"r") rivi = AnalyseCode.readline() while (len(rivi) > 0): # while begins _____________________________ rivi=rivi.rstrip(' \n') rivi = np.fromstring(rivi, dtype=float, sep=' ') j=(rivi != dummy).nonzero() rivi=rivi[j] T0=rivi[0] TSPAN=rivi[1] ALLF=rivi[2:2+K1] BETA=rivi[2+K1:] WriteALLF(ALLF) FNOW=ALLF[Val-1] TT=(1./FNOW)*np.arange(1000)/1000. # Note: 1001 out! G,H1,H2,H3,H4,H5,H6,PL=LinearModel(TT,BETA,0.) # - Checking that linear and nonlinear model give same results! if (Check != 0): NONLINBETA=np.hstack((ALLF,BETA)) GCHECK,ny1,ny2,ny3,ny4,ny5,ny6,ny7=NonLinearModel(TT,NONLINBETA,0.) print('G_linear/G_nonlinear =',G/GCHECK) if (Val == 1): HH=H1 if (Val == 2): HH=H2 if (Val == 3): HH=H3 if (Val == 4): HH=H4 if (Val == 5): HH=H5 if (Val == 6): HH=H6 if (i == 0): ColF=FNOW ColA=np.max(HH)-np.min(HH) Which=0. # Which=1 gives maxima. Any other value gives minima. AMin,BMin,NMin\ =HarmonicEpochs(Which,i,TT,HH,0.) ColTMIN1=AMin ColTMIN2=BMin Which=1. # Which=1 gives maxima. Any other value gives minima. AMax,BMax,NMax\ =HarmonicEpochs(Which,i,TT,HH,0.) ColTMAX1=AMax ColTMAX2=BMax if (i != 0): ColF=np.append(ColF,FNOW) ColA=np.append(ColA,np.max(HH)-np.min(HH)) if (Check != 0): print('Round ',i,' minima') Which=0. # Which=1 gives maxima. Any other value gives minima. aMin,bMin,nMin\ =HarmonicEpochs(Which,i,TT,HH,0.) cMin,dMin=EpochSelect(AMin,BMin,NMin,aMin,bMin,nMin,FNOW,0.) ColTMIN1=np.append(ColTMIN1,cMin) ColTMIN2=np.append(ColTMIN2,dMin) if (Check != 0): print('Round ',i,' maxima') Which=1. # Which=1 gives maxima. Any other value gives minima. aMax,bMax,nMax\ =HarmonicEpochs(Which,i,TT,HH,0.) cMax,dMax=EpochSelect(AMax,BMax,NMax,aMax,bMax,nMax,FNOW,0.) ColTMAX1=np.append(ColTMAX1,cMax) ColTMAX2=np.append(ColTMAX2,dMax) rivi = AnalyseCode.readline() # while ends __________________________ i=i+1 # Next line AnalyseCode.close() # Epoch correction, where dummies remain dummies! j=(ColTMIN1 != dummy).nonzero() if (np.size(j) > 0): ColTMIN1[j]=ColTMIN1[j]+T0 j=(ColTMIN2 != dummy).nonzero() if (np.size(j) > 0): ColTMIN2[j]=ColTMIN2[j]+T0 j=(ColTMAX1 != dummy).nonzero() if (np.size(j) > 0): ColTMAX1[j]=ColTMAX1[j]+T0 j=(ColTMAX2 != dummy).nonzero() if (np.size(j) > 0): ColTMAX2[j]=ColTMAX2[j]+T0 if (Check != 0): # Program check print(' =========================================================') print('1/ColF =',1./ColF[0],'+/-',np.std(1./ColF[1:])) print('ColF =',ColF[0],'+/-',np.std(ColF[1:])) print('ColA =',ColA[0],'+/-',np.std(ColA[1:])) print('ColTMIN1=',ColTMIN1[0],'+/-',np.std(ColTMIN1[1:])) ny=(ColTMIN2[1:] != dummy).nonzero() if (np.size(ny) > 1): print('ColTMIN2=',ColTMIN2[0],'+/-',np.std(ColTMIN2[1:])) else: print('ColTMIN2=',' ... ','+/-',' ... ') print('ColTMAX1=',ColTMAX1[0],'+/-',np.std(ColTMAX1[1:])) ny=(ColTMAX2[1:] != dummy).nonzero() if (np.size(ny) > 1): print('ColTMAX2=',ColTMAX2[0],'+/-',np.std(ColTMAX2[1:])) else: print('ColTMAX2=',' ... ','+/-',' ... ') return ColF,ColA,ColTMIN1,ColTMIN2,ColTMAX1,ColTMAX2 # __________________________________________________________________________ # - Selecting bootstrap primary and secondary minimum/maximum epochs closest # to primary and secondary minimum/maximum epochs for original data. # __________________________________________________________________________ def EpochSelect(A,B,N,a,b,n,FNOW,Check): dummy =dcmInput(4) if ((N == 2) & (n == 2)): # Original: 2 epochs, Bootstrap: 2 epochs Case=1 ny1=a-1.*(1./FNOW) ny2=a ny3=a+1.*(1./FNOW) ny4=b-1.*(1./FNOW) ny5=b ny6=b+1.*(1./FNOW) Candidates=np.hstack((ny1,ny2,ny3,ny4,ny5,ny6)) if (np.min(np.abs(Candidates-A)) <= np.min(np.abs(Candidates-B))): ny=np.abs(Candidates-A) ny=(ny == np.min(ny)).nonzero() c=Candidates[ny] ny=np.abs(Candidates-B) ny=(ny == np.min(ny)).nonzero() d=Candidates[ny] if (np.min(np.abs(Candidates-A)) > np.min(np.abs(Candidates-B))): ny=np.abs(Candidates-B) ny=(ny == np.min(ny)).nonzero() d=Candidates[ny] ny=np.abs(Candidates-A) ny=(ny == np.min(ny)).nonzero() c=Candidates[ny] if ((N == 2) & (n == 1)): # Original: 2 epochs, Bootstrap: 1 epoch Case=2 ny1=a-1.*(1./FNOW) ny2=a ny3=a+1.*(1./FNOW) Candidates=np.hstack((ny1,ny2,ny3)) if (np.min(np.abs(Candidates-A)) <= np.min(np.abs(Candidates-B))): ny=np.abs(Candidates-A) ny=(ny == np.min(ny)).nonzero() c=Candidates[ny] d=dummy if (np.min(np.abs(Candidates-A)) > np.min(np.abs(Candidates-B))): ny=np.abs(Candidates-B) ny=(ny == np.min(ny)).nonzero() c=dummy d=Candidates[ny] if ((N == 1) & (n == 2)): # Original: 1 epoch, Bootstrap: 2 epochs Case=3 ny1=a-1.*(1./FNOW) ny2=a ny3=a+1.*(1./FNOW) ny4=b-1.*(1./FNOW) ny5=b ny6=b+1.*(1./FNOW) Candidates=np.hstack((ny1,ny2,ny3,ny4,ny5,ny6)) ny=np.abs(Candidates-A) ny=(ny == np.min(ny)).nonzero() c=Candidates[ny] d=dummy if ((N == 1) & (n == 1)): # Original: 1 epoch, Bootstrap: 1 epoch Case=4 ny1=a-1.*(1./FNOW) ny2=a ny3=a+1.*(1./FNOW) Candidates=np.hstack((ny1,ny2,ny3)) ny=np.abs(Candidates-A) ny=(ny == np.min(ny)).nonzero() c=Candidates[ny] d=dummy if (Check != 0.): print('CASE =', Case ,' _______________________________________') print(A,B) print(a,b) print('Candidates =',Candidates) print(c,d) return c,d # __________________________________________________________________________ # - Solving minima or maxima epochs of harmonic HH=HH(TT). # "Which = 1" gives maxima. Any other "Which" value gives minima. # __________________________________________________________________________ def HarmonicEpochs(Which,i,TT,HH,Check): dummy=dcmInput(4) if (Which == 1): ny=((np.roll(HH, 1) < HH) & (np.roll(HH, -1) < HH)).nonzero() else: ny=((np.roll(HH, 1) > HH) & (np.roll(HH, -1) > HH)).nonzero() HowMany=np.size(ny) if (HowMany == 1): # One epoch only E=np.append(TT[ny],dummy) Y=np.append(HH[ny],dummy) if (HowMany == 2): # Two epochs E=TT[ny] Y=HH[ny] if (Which == 1): k=IndexSortOrder(-1.*Y) # Maxima to descreasing order if (Which != 1): k=IndexSortOrder(+1*Y) # Minima to increasing order E=E[k] Y=Y[k] if (Check != 0.): print('i =',i,'E[0]=',E[0],'E[1]=',E[1]) return E[0],E[1],HowMany # __________________________________________________________________________ # - Plotting parameter pairs. # __________________________________________________________________________ def ParameterPair(Place,X,Y,txt1,txt2,Diagonal,AboveZero,ShowLimits,TakeIt): pl.axes(Place) x1=np.min(X) x2=np.max(X) lev=x2-x1 x1=x1-0.1*lev x2=x2+0.1*lev if ((AboveZero != 0) and (x1 < 0)): # Amplitude never negative! x1=0 pl.xlim([x1,x2]) y1=np.min(Y) y2=np.max(Y) lev=y2-y1 y1=y1-0.1*lev y2=y2+0.1*lev if ((AboveZero != 0) and (y1 < 0)): # Amplitude never negative! y1=0 pl.ylim([y1,y2]) # suppress scientific notations ... pl.ticklabel_format(useOffset=False, style='plain') xticksnyt=[x1,x2] pl.xticks(xticksnyt,fontsize=7) yticksnyt=[y1,y2] pl.yticks(yticksnyt,fontsize=7) if (Diagonal == 1): pl.plot([x1,x2],[x1,x2],'g',linewidth=4) if (ShowLimits != 0): ErrorLimits(X,Y,TakeIt) pl.plot(X[1:],Y[1:],'or',ms=3) pl.plot(X[0],Y[0],'ob',ms=5) pl.xlabel(txt1,fontsize=12) pl.ylabel(txt2,fontsize=12) return # __________________________________________________________________________ # - Plotting 1,2,3 sigma error ellipses. # _____________________________________________________________________________ def ErrorLimits(X,Y,TakeIt): zz=2.*np.pi*np.arange(501)/500. a1=np.std(X[1:]) a2=np.std(Y[1:]) xx=X[0]+1.*a1*np.cos(zz) yy=Y[0]+1.*a2*np.sin(zz) pl.plot(xx,yy,':k') xx=X[0]+2.*a1*np.cos(zz) yy=Y[0]+2.*a2*np.sin(zz) pl.plot(xx,yy,':k') xx=X[0]+3.*a1*np.cos(zz) yy=Y[0]+3.*a2*np.sin(zz) pl.plot(xx,yy,':k') i= (yy >= xx).nonzero() # Version 2: all lines below .... if ((np.size(i) > 0) & (TakeIt != 0)): ny1='%1i' %(TakeIt) ny2='%1i' %(TakeIt+1) txt=ny1+' and '+ny2+' frequencies intersect'+'\n' Notes(txt) return # _____________________________________________________________________________ # - Graphical check of frequencies and amplitudes of harmonics. # __________________________________________________________________________ def CheckFreqAmp(file1,Check): Tag =dcmInput(1) K1 =dcmInput(5) if (K1 >= 2.): Val=1 F1,A1,ColTMIN1,ColTMIN2,ColTMAX1,ColTMAX2\ =AnalyseHarmonic(Val,file1,Check) Val=2 F2,A2,ColTMIN1,ColTMIN2,ColTMAX1,ColTMAX2\ =AnalyseHarmonic(Val,file1,Check) Place=[0.15+0.00,0.87,0.15,0.12] ParameterPair(Place,F1,F2,'$f_1$','$f_2$',1,0,1,1) Place=[0.15+0.27,0.87,0.15,0.12] ParameterPair(Place,A1,A2,'$A_1$','$A_2$',0,1,1,0) if (K1 >= 3.): Val=2 F2,A2,ColTMIN1,ColTMIN2,ColTMAX1,ColTMAX2\ =AnalyseHarmonic(Val,file1,Check) Val=3 F3,A3,ColTMIN1,ColTMIN2,ColTMAX1,ColTMAX2\ =AnalyseHarmonic(Val,file1,Check) Place=[0.15+0.00,0.67,0.15,0.12] ParameterPair(Place,F2,F3,'$f_2$','$f_3$',1,0,1,2) Place=[0.15+0.27,0.67,0.15,0.12] ParameterPair(Place,A2,A3,'$A_2$','$A_3$',0,1,1,0) if (K1 >= 4.): Val=3 F3,A3,ColTMIN1,ColTMIN2,ColTMAX1,ColTMAX2\ =AnalyseHarmonic(Val,file1,Check) Val=4 F4,A4,ColTMIN1,ColTMIN2,ColTMAX1,ColTMAX2\ =AnalyseHarmonic(Val,file1,Check) Place=[0.15+0.00,0.47,0.15,0.12] ParameterPair(Place,F3,F4,'$f_3$','$f_4$',1,0,1,3) Place=[0.15+0.27,0.47,0.15,0.12] ParameterPair(Place,A3,A4,'$A_3$','$A_4$',0,1,1,0) if (K1 >= 5.): Val=4 F4,A4,ColTMIN1,ColTMIN2,ColTMAX1,ColTMAX2\ =AnalyseHarmonic(Val,file1,Check) Val=5 F5,A5,ColTMIN1,ColTMIN2,ColTMAX1,ColTMAX2\ =AnalyseHarmonic(Val,file1,Check) Place=[0.15+0.00,0.27,0.15,0.12] ParameterPair(Place,F4,F5,'$f_4$','$f_5$',1,0,1,4) Place=[0.15+0.27,0.27,0.15,0.12] ParameterPair(Place,A4,A5,'$A_4$','$A_5$',0,1,1,0) if (K1 >= 6.): Val=5 F5,A5,ColTMIN1,ColTMIN2,ColTMAX1,ColTMAX2\ =AnalyseHarmonic(Val,file1,Check) Val=6 F6,A6,ColTMIN1,ColTMIN2,ColTMAX1,ColTMAX2\ =AnalyseHarmonic(Val,file1,Check) Place=[0.15+0.00,0.07,0.15,0.12] ParameterPair(Place,F5,F6,'$f_5$','$f_6$',1,0,1,5) Place=[0.15+0.27,0.07,0.15,0.12] ParameterPair(Place,A5,A6,'$A_5$','$A_6$',0,1,1,0) kuva=Tag+'fA.eps' pl.savefig(kuva) pl.clf() return # __________________________________________________________________________ # - Storing main results. # __________________________________________________________________________ def StoreResults(T0,T,Y,FINALBETA,file7): Tag =dcmInput(1) dummy =dcmInput(4) K1 =dcmInput(5) K2 =dcmInput(6) K3 =dcmInput(7) nL =dcmInput(8) nS =dcmInput(9) PMIN =dcmInput(12) PMAX =dcmInput(13) Rounds =dcmInput(14) PrintScreen=dcmInput(24) filenow1=Tag+'Params.dat' filenow2=Tag+'Notes.dat' Storefile=open(filenow1,"w") txt="%10s %16i"%('n',np.size(T)) Storefile.write("%60s\n" %(txt)) txt="%10s %0.10e"%('T1',T0) Storefile.write("%60s\n" %(txt)) txt="%10s %0.10e"%('DT',np.max(T)-np.min(T)) Storefile.write("%60s\n" %(txt)) txt="%10s %0.10e"%('my',np.mean(Y)) Storefile.write("%60s\n" %(txt)) txt="%10s %0.10e"%('sy',np.std(Y)) Storefile.write("%60s\n" %(txt)) txt="%10s %0.10e"%('sigma',np.mean(EY)) Storefile.write("%60s\n" %(txt)) txt="%10s %16i"%('K1',K1) Storefile.write("%60s\n" %(txt)) txt="%10s %16i"%('K2',K2) Storefile.write("%60s\n" %(txt)) txt="%10s %16i"%('K3',K3) Storefile.write("%60s\n" %(txt)) txt="%10s %16i"%('p',K1*(2*K2+1)+(K3+1)) Storefile.write("%60s\n" %(txt)) txt="%10s %0.10e"%('PMIN',PMIN) Storefile.write("%60s\n" %(txt)) txt="%10s %0.10e"%('PMAX',PMAX) Storefile.write("%60s\n" %(txt)) txt="%10s %16i"%('nL',nL) Storefile.write("%60s\n" %(txt)) txt="%10s %16i"%('nS',nS) Storefile.write("%60s\n" %(txt)) G,H1,H2,H3,H4,H5,H6,PL=NonLinearModel(T,FINALBETA,0.) ny1=np.sum(((Y-G)*(Y-G))/(EY*EY)) ny2=np.sqrt(ny1/(1.*np.size(T))) txt="%10s %0.10e %12s %0.10e"%('CHI2',ny1,' gives ZMIN',ny2) Storefile.write("%60s\n" %(txt)) ny1=np.sum((Y-G)*(Y-G)) ny2=np.sqrt(ny1/(1.*np.size(T))) txt="%10s %0.10e %12s %0.10e"%('R',ny1,'gives ZMIN',ny2) Storefile.write("%60s\n" %(txt)) Check=0 for Val in range(1,K1+1): Har='%1i' %(Val) F,A,TMIN1,TMIN2,TMAX1,TMAX2\ =AnalyseHarmonic(Val,file7,Check) txt='..........................................................' Storefile.write("%60s\n" %(txt)) txt="%10s %0.10e %3s %0.10e"%('F'+Har,F[0],'+/-',np.std(F[1:])) Storefile.write("%60s\n" %(txt)) txt="%10s %0.10e %3s %0.10e"%('P'+Har,1./F[0],'+/-',np.std(1./F[1:])) Storefile.write("%60s\n" %(txt)) NowSN=(A[0])/np.mean(EY) txt="%10s %0.10e %3s %0.10e %6s %10.2f"%\ ('A'+Har,A[0],'+/-',np.std(A[1:]),'SN',NowSN) Storefile.write("%78s\n" %(txt)) txt="%10s %0.10e %3s %0.10e"%\ ('T'+Har+'MIN1',TMIN1[0],'+/-',np.std(TMIN1[1:])) Storefile.write("%60s\n" %(txt)) ny=(TMIN2[1:] != dummy).nonzero() if (np.size(ny) > 1): ny1='%10i' %(np.size(ny)) ny2='%10i' %(np.floor(Rounds)) ny3='('+ny1+' / '+ny2+')' ny3=ny3.replace(" ", "") txt="%10s %0.10e %3s %0.10e %10s"%\ ('T'+Har+'MIN2',TMIN2[0],'+/-',np.std(TMIN2[1:]),ny3) Storefile.write("%71s\n" %(txt)) else: txt="%10s %16s %3s %15s"%('T'+Har+'MIN2','...','+/-','...') Storefile.write("%59s\n" %(txt)) txt="%10s %0.10e %3s %0.10e"%\ ('T'+Har+'MAX1',TMAX1[0],'+/-',np.std(TMAX1[1:])) Storefile.write("%60s\n" %(txt)) ny=(TMAX2[1:] != dummy).nonzero() if (np.size(ny) > 1): ny1='%10i' %(np.size(ny)) ny2='%10i' %(np.floor(Rounds)) ny3='('+ny1+' / '+ny2+')' ny3=ny3.replace(" ", "") txt="%10s %0.10e %3s %0.10e %10s"%\ ('T'+Har+'MAX2',TMAX2[0],'+/-',np.std(TMAX2[1:]),ny3) Storefile.write("%71s\n" %(txt)) else: txt="%10s %16s %3s %15s"%('T'+Har+'MAX2','...','+/-','...') Storefile.write("%59s\n" %(txt)) NN=39 Info=['T1','DeltaT',\ 'F1','F2','F3','F4','F5','F6',\ 'B11','C11','B12','C12',\ 'B21','C21','B22','C22',\ 'B31','C31','B32','C32',\ 'B41','C41','B42','C42',\ 'B51','C51','B52','C52',\ 'B61','C61','B62','C62',\ 'M0','M1','M2','M3','M4','M5','M6'] I=1 txt='...........................................................' Storefile.write("%59s\n" %(txt)) txt=' i BETA[i]' Storefile.write("%43s\n" %(txt)) for i in range(2,NN-1): # From "2" = T1 and DT given elsewhere ny=np.loadtxt(file7,skiprows=0,usecols=(i, )) if (ny[0] != dummy): if (ny[0] < 0): # minus sign = Yes txt="%4i %10s %0.5e %0.5e "%(I,Info[i],ny[0],np.std(ny[1:])) if (ny[0] >= 0): # minus sign = No txt="%4i %10s %0.5e %0.5e "%(I,Info[i],ny[0],np.std(ny[1:])) Storefile.write("%59s\n" %(txt)) I=I+1 Storefile.close() if (PrintScreen == 1): # Print results to screen PrintIt='cat '+filenow1 # -"- os.system(PrintIt) # -"- PrintIt='cat '+filenow2 # Version 2 os.system(PrintIt) # Version 2 return # __________________________________________________________________________ # - Separating simulated and detected results. # __________________________________________________________________________ def SeparateXY(X,Y): ny=np.arange(np.size(X)) ny=ny/2. ny=ny-np.floor(ny) j=(ny == 0).nonzero() Xsim=X[j] Ysim=Y[j] j=(ny != 0).nonzero() Xdet=X[j] Ydet=Y[j] return Xsim,Xdet,Ysim,Ydet # __________________________________________________________________________ # - Frequency information. # __________________________________________________________________________ def FrequencyInfo(Txt,Fdet,Fsim,Fapu,Aapu,SimDF,SimDA,EEcode): ny1=np.mean(np.abs(Fdet-Fsim)/Fsim) ny2=ny1 # May, or may not, change! ny3=ny1 # -"- j1=(Fapu != 1).nonzero() # No Fapu=1 if (np.size(j1) > 0): ny2=np.mean(np.abs(Fdet[j1]-Fsim[j1])/Fsim[j1]) j2=((Fapu == 0) & (Aapu == 0.)).nonzero() # No Fapu=1 OR Apu=1 if (np.size(j2) > 0): ny3=np.mean(np.abs(Fdet[j2]-Fsim[j2])/Fsim[j2]) EEcode.write("%4s, %5s %5i, %26s %0.3e %8s\n"\ %(Txt,'m=',np.size(Fdet),'mean relative accuracy = ',\ ny1,' for all')) if (np.size(j1) > 0): EEcode.write("%4s, %5s %5i, %26s %0.3e %10s %0.3e\n" \ %(Txt,'m=',np.size(j1),'mean relative accuracy = ',\ ny2,'for DF > ',SimDF)) else: EEcode.write("%10s \n" %('m=0 (No result)',)) if (np.size(j2) > 0): EEcode.write("%4s, %5s %5i, %26s %0.3e %10s %0.3e %8s %0.3e %10s\n" \ %(Txt,'m=',np.size(j2),'mean relative accuracy = ',\ ny3,'for DF < ',SimDF,' or DA < ',SimDA,' excluded')) else: EEcode.write("%10s \n" %('m=0 (No result)',)) return # _____________________________________________________________________________ # - Plot simulated versus detected frequencies and amplitudes. # _____________________________________________________________________________ def fAplot(Fsim,Fdet,Asim,Adet,Fapu,Aapu,paikka1,paikka2,\ FMIN,FMAX,txt1,txt2,txt3,txt4,Sym1): xx=[FMIN,FMAX] yy=[FMIN,FMAX] pl.axes(paikka1) # Frequency pl.xlim([FMIN,FMAX]) # No pl.ylim !!! pl.xticks(fontsize=7) pl.yticks(fontsize=7) pl.xlabel(txt1,fontsize=9) pl.ylabel(txt2,fontsize=9) pl.plot(Fsim,Fdet,Sym1,ms=6) j=(Fapu != 0).nonzero() if (np.size(j) > 0): pl.plot(Fsim[j],Fdet[j],'Dk',ms=10,markerfacecolor="None") j=(Aapu != 0).nonzero() if (np.size(j) > 0): pl.plot(Fsim[j],Fdet[j],'ok',ms=10,markerfacecolor="None") pl.plot(xx,yy,'k') pl.axes(paikka2) # Amplitude pp=[0.,1.05*np.max(Asim)] qq=pp pl.xlim([pp[0],pp[1]]) # No pl.ylim !!! pl.xticks(fontsize=7) pl.yticks(fontsize=7) pl.xlabel(txt3,fontsize=9) pl.ylabel(txt4,fontsize=9) pl.plot(Asim,Adet,Sym1,ms=6) j=(Fapu != 0).nonzero() if (np.size(j) > 0): pl.plot(Asim[j],Adet[j],'Dk',ms=10,markerfacecolor="None") j=(Aapu != 0).nonzero() if (np.size(j) > 0): pl.plot(Asim[j],Adet[j],'ok',ms=10,markerfacecolor="None") pl.plot(pp,qq,'k') return # _____________________________________________________________________________ # - Graphical check of many model simulations: even one error cause damage # _____________________________________________________________________________ def SimManyCheck(Betafile): Tag =dcmInput(1) dummy =dcmInput(4) K1 =dcmInput(5) PMIN =dcmInput(12) PMAX =dcmInput(13) SimDF =dcmInput(22) SimDA =dcmInput(23) PrintScreen=dcmInput(24) FMAX=1./PMIN FMIN=1./PMAX DF=FMAX-FMIN Val=1 X,Y,ColTMIN1,ColTMIN2,\ ColTMAX1,ColTMAX2=AnalyseHarmonic(Val,BetaFile,0.) F1sim,F1det,A1sim,A1det=SeparateXY(X,Y) if (K1 >= 2): Val=2 X,Y,ColTMIN1,ColTMIN2,\ ColTMAX1,ColTMAX2=AnalyseHarmonic(Val,BetaFile,0.) F2sim,F2det,A2sim,A2det=SeparateXY(X,Y) if (K1 >= 3): Val=3 X,Y,ColTMIN1,ColTMIN2,\ ColTMAX1,ColTMAX2=AnalyseHarmonic(Val,BetaFile,0.) F3sim,F3det,A3sim,A3det=SeparateXY(X,Y) if (K1 >= 4): Val=4 X,Y,ColTMIN1,ColTMIN2,\ ColTMAX1,ColTMAX2=AnalyseHarmonic(Val,BetaFile,0.) F4sim,F4det,A4sim,A4det=SeparateXY(X,Y) if (K1 >= 5): Val=5 X,Y,ColTMIN1,ColTMIN2,\ ColTMAX1,ColTMAX2=AnalyseHarmonic(Val,BetaFile,0.) F5sim,F5det,A5sim,A5det=SeparateXY(X,Y) if (K1 >= 6): Val=6 X,Y,ColTMIN1,ColTMIN2,\ ColTMAX1,ColTMAX2=AnalyseHarmonic(Val,BetaFile,0.) F6sim,F6det,A6sim,A6det=SeparateXY(X,Y) Fapu=0.*F1sim Aapu=0.*A1sim for i in range(np.size(F1sim)): allF=F1sim[i] allA=A1sim[i] if (K1 >= 2): allF=np.append(allF,F2sim[i]) allA=np.append(allA,A2sim[i]) if (K1 >= 3): allF=np.append(allF,F3sim[i]) allA=np.append(allA,A3sim[i]) if (K1 >= 4): allF=np.append(allF,F4sim[i]) allA=np.append(allA,A4sim[i]) if (K1 >= 5): allF=np.append(allF,F5sim[i]) allA=np.append(allA,A5sim[i]) if (K1 >= 6): allF=np.append(allF,F6sim[i]) allA=np.append(allA,A6sim[i]) if (K1 >= 2): # Some frequencies to compare. DallF=np.abs(np.roll(allF,1)-allF) DallF=DallF[1:] else: # No frequencies to compare. DallF=FMAX-FMIN if (np.min(DallF/DF) <= SimDF): Fapu[i]=1. if ((np.min(allA)/np.max(allA) < SimDA) & (K1 >=2)): # Some amps to compare Aapu[i]=1. Check=0 # Program check if (Check != 0): print('.....................................................') print('FREQUENCIES') print('all f ........... = ',allF) print('differences ..... = ',DallF) print('differences in DF = ',DallF/DF) print('Fapu[i] ......... = ',Fapu[i]) print('AMPLITUDES') print('all a ........... = ',allA) print('allA/a_max ...... = ',allA/np.max(allA)) print('Aapu[i] ......... = ',Aapu[i]) Lev=0.20 Kor=0.09 filenow=Tag+'ManyfA.dat' EEcode=open(filenow,"w") if (K1 >= 1): # ___________________________________________________________ FrequencyInfo('f1',F1det,F1sim,Fapu,Aapu,SimDF,SimDA,EEcode) paikka1=[0.15,0.85,Lev,Kor] paikka2=[0.45,0.85,Lev,Kor] fAplot(F1sim,F1det,A1sim,A1det,Fapu,Aapu,paikka1,paikka2,FMIN,FMAX,\ '$f_{\mathrm{1,sim}}$','$f_{\mathrm{1,det}}$',\ '$A_{\mathrm{1,sim}}$','$A_{\mathrm{1,det}}$','or') if (K1 >= 2): # ___________________________________________________________ FrequencyInfo('f2',F2det,F2sim,Fapu,Aapu,SimDF,SimDA,EEcode) paikka1=[0.15,0.69,Lev,Kor] paikka2=[0.45,0.69,Lev,Kor] fAplot(F2sim,F2det,A2sim,A2det,Fapu,Aapu,paikka1,paikka2,FMIN,FMAX,\ '$f_{\mathrm{2,sim}}$','$f_{\mathrm{2,det}}$',\ '$A_{\mathrm{2,sim}}$','$A_{\mathrm{2,det}}$','ob') if (K1 >= 3): # ___________________________________________________________ FrequencyInfo('f3',F3det,F3sim,Fapu,Aapu,SimDF,SimDA,EEcode) paikka1=[0.15,0.53,Lev,Kor] paikka2=[0.45,0.53,Lev,Kor] fAplot(F3sim,F3det,A3sim,A3det,Fapu,Aapu,paikka1,paikka2,FMIN,FMAX,\ '$f_{\mathrm{3,sim}}$','$f_{\mathrm{3,det}}$',\ '$A_{\mathrm{3,sim}}$','$A_{\mathrm{3,det}}$','og') if (K1 >= 4): # ___________________________________________________________ FrequencyInfo('f4',F4det,F4sim,Fapu,Aapu,SimDF,SimDA,EEcode) paikka1=[0.15,0.37,Lev,Kor] paikka2=[0.45,0.37,Lev,Kor] fAplot(F4sim,F4det,A4sim,A4det,Fapu,Aapu,paikka1,paikka2,FMIN,FMAX,\ '$f_{\mathrm{4,sim}}$','$f_{\mathrm{4,det}}$',\ '$A_{\mathrm{4,sim}}$','$A_{\mathrm{4,det}}$','oy') if (K1 >= 5): # ___________________________________________________________ FrequencyInfo('f5',F5det,F5sim,Fapu,Aapu,SimDF,SimDA,EEcode) paikka1=[0.15,0.21,Lev,Kor] paikka2=[0.45,0.21,Lev,Kor] fAplot(F5sim,F5det,A5sim,A5det,Fapu,Aapu,paikka1,paikka2,FMIN,FMAX,\ '$f_{\mathrm{5,sim}}$','$f_{\mathrm{5,det}}$',\ '$A_{\mathrm{5,sim}}$','$A_{\mathrm{5,det}}$','om') if (K1 >= 6): # ___________________________________________________________ FrequencyInfo('f6',F6det,F6sim,Fapu,Aapu,SimDF,SimDA,EEcode) paikka1=[0.15,0.05,Lev,Kor] paikka2=[0.45,0.05,Lev,Kor] fAplot(F6sim,F6det,A6sim,A6det,Fapu,Aapu,paikka1,paikka2,FMIN,FMAX,\ '$f_{\mathrm{6,sim}}$','$f_{\mathrm{6,det}}$',\ '$A_{\mathrm{6,sim}}$','$A_{\mathrm{6,det}}$','oc') EEcode.close() if (PrintScreen == 1): PrintIt='cat '+filenow os.system(PrintIt) kuva=Tag+'Many.eps' pl.savefig(kuva) return # _____________________________________________________________________________ # - Saving residuals. # _____________________________________________________________________________ def SaveResiduals(T0,T,Y,EY,G): Tag=dcmInput(1) OUTfile=Tag+'Residuals.dat' Ecode=open(OUTfile,"w") for i in range(np.size(T)): Ecode.write("%0.10e %0.10e %0.10e \n"\ %(T[i]+T0,Y[i]-G[i],EY[i])) Ecode.close() return # _____________________________________________________________________________ # - Saving model. # _____________________________________________________________________________ def SaveModel(T0,T,Y,EY,G): Tag=dcmInput(1) OUTfile=Tag+'Model.dat' Ecode=open(OUTfile,"w") for i in range(np.size(T)): Ecode.write("%0.10e %0.10e %0.10e %0.10e \n"\ %(T[i]+T0,Y[i],EY[i],G[i])) Ecode.close() return # _____________________________________________________________________________ # # MAIN PROGRAM # ============= file1=dcmInput(3) # Real data file Tfile=file1 # avoids input=output RealData=dcmInput(2) # Modes 1, 2 and 3 SimMany =dcmInput(20) # -"- # _____________________________________________________________________________ if (SimMany != 1): # One sample of data NotesLineOne() # Empty error notes file if (RealData == 1.): # Mode 1: Real data T0,T,Y,EY=DataRead(file1) # -"- if (RealData != 1.): # Mode 2: Simulates one BETA,file1=SIMdata(0.,Tfile) # Simulated data T0,T,Y,EY=DataRead(file1) # -"- LongFrequencyInterval() # Long search freqs F1,F2,F3,F4,F5,F6,Z,ALLF=Zgram(T,Y,EY,0.) # Long search WritePeriodogram(F1,F2,F3,F4,F5,F6,Z,'z1.dat') # Store periodograms PCheck(F1,F2,F3,F4,F5,F6,Z,'Long: ') # Prints long interval .. ShortFrequencyInterval(ALLF) # Short search freqs F1,F2,F3,F4,F5,F6,Z,ALLF=Zgram(T,Y,EY,1.) # Short search WriteALLF(ALLF) # ALLF from short search WritePeriodogram(F1,F2,F3,F4,F5,F6,Z,'z2.dat') # Store periodograms PCheck(F1,F2,F3,F4,F5,F6,Z,'Short: ') # Prints short interval .. PlotPeriodograms('z1.dat','z2.dat') # Plots periodograms BETANOW=dcmInput(4) # BETA can be dummy, SolveBeta=1 # because it is solved. kuva=dcmInput(1)+'gdet.eps' # CheckPlot(T0,T,Y,EY,BETANOW,kuva,SolveBeta) # ALLF from short search BetaFile=dcmInput(1)+'AllBeta.dat' # Free parameter file Style=1 # "w"=writes OrigG,FINALBETA=dcmModel(T0,T,Y,EY,BetaFile,Style)# Original sample PhasePlots(T0,T,Y,EY) BootStrap(T0,T,Y,EY,OrigG,BetaFile,0.) # Bootstrap samples CheckFreqAmp(BetaFile,0.) # Check freqs & amps StoreResults(T0,T,Y,FINALBETA,BetaFile) # Store results SaveResiduals(T0,T,Y,EY,OrigG) # Store residuals SaveModel(T0,T,Y,EY,OrigG) # Store data & model # _____________________________________________________________________________ if (SimMany == 1): # Mode 3: Simulates many NotesLineOne() # Empty error notes file BetaFile=dcmInput(1)+'AllBeta.dat' # Free parameter output SimRounds=dcmInput(21) # Number of samples for i in range(SimRounds): # Model loop begins if (dcmInput(24) == 1): # Information to screen print('Model =',i+1,'out of ',SimRounds) # -"- BETA,file1=SIMdata(0.,Tfile) # Simulated data T0,T,Y,EY=DataRead(file1) # -"- ALLF=ReadALLF() # SIMdata stores ALLF if (i == 0): Style=1 # 1st round: w"=writes else: Style=2 # Others: "a"=appends WriteBeta(T0,T,ALLF,BETA,BetaFile,Style,0.) # Store simulated Beta LongFrequencyInterval() # Long search freqs F1,F2,F3,F4,F5,F6,Z,ALLF=Zgram(T,Y,EY,0.) # Long search WritePeriodogram(F1,F2,F3,F4,F5,F6,Z,'z1.dat')# Store periodograms PCheck(F1,F2,F3,F4,F5,F6,Z,'Long :') # Long search results ... ShortFrequencyInterval(ALLF) # Short search freqs F1,F2,F3,F4,F5,F6,Z,ALLF=Zgram(T,Y,EY,1.) # Short search WritePeriodogram(F1,F2,F3,F4,F5,F6,Z,'z2.dat')# Store periodiograms PCheck(F1,F2,F3,F4,F5,F6,Z,'Short :') # Short search results ... PlotPeriodograms('z1.dat','z2.dat') # Plots periodogram WriteALLF(ALLF) # ALLF from short search BETANOW=dcmInput(4) # BETA can be dummy, SolveBeta=1 # because BETA is solved. kuva=dcmInput(1)+'gdet.eps' # CheckPlot(T0,T,Y,EY,BETANOW,kuva,SolveBeta) # ALLF from short search Style=2. # All rounds:"a"=appends OrigG,FINALBETA=dcmModel(T0,T,Y,EY,BetaFile,Style) # Model loop ends SimManyCheck(BetaFile) # Results checked. # _____________________________________________________________________________