From 579a767ba2d2522b054d0ace466cd6156ee4cd0f Mon Sep 17 00:00:00 2001 From: Buntu Ngcebetsha Date: Wed, 17 Apr 2024 17:01:32 +0200 Subject: [PATCH 1/7] make imports work --- .gitignore | 4 + ClassDynSpecMS.py | 307 ++++++++++++++++++++++---------------------- ClassSaveResults.py | 82 ++++++------ ms2dynspec.py | 28 ++-- 4 files changed, 213 insertions(+), 208 deletions(-) diff --git a/.gitignore b/.gitignore index ef25bf1..a2c2f65 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,7 @@ *.pyc *.o *.so +venv +*.txt +*.csv +*.obj diff --git a/ClassDynSpecMS.py b/ClassDynSpecMS.py index 7797ff7..b523a2c 100644 --- a/ClassDynSpecMS.py +++ b/ClassDynSpecMS.py @@ -10,7 +10,7 @@ log=logger.getLogger("DynSpecMS") from DDFacet.Array import shared_dict from DDFacet.Other import AsyncProcessPool - + from DDFacet.Other import Multiprocessing from DDFacet.Other import ModColor from DDFacet.Other.progressbar import ProgressBar @@ -21,7 +21,8 @@ import os from killMS.Other import reformat from DDFacet.Other import AsyncProcessPool -from .dynspecms_version import version +#from .dynspecms_version import version +from dynspecms_version import version import glob from astropy.io import fits from astropy.wcs import WCS @@ -67,8 +68,8 @@ def __init__(self, ColName="DATA", TChunkHours=0., ModelName="PREDICT_KMS", - UVRange=[1.,1000.], - ColWeights=None, + UVRange=[1.,1000.], + ColWeights=None, SolsName=None, FileCoords="Transient_LOTTS.csv", Radius=3., @@ -112,8 +113,8 @@ def __init__(self, if ListMSName is None: print(ModColor.Str("WORKING IN REPLOT MODE"), file=log) self.Mode="Plot" - - + + self.Radius=Radius self.ImageI = ImageI self.ImageV = ImageV @@ -123,14 +124,14 @@ def __init__(self, # identify version in logs print("DynSpecMS version %s starting up" % version(), file=log) self.FileCoords=FileCoords - + if self.Mode=="Spec": self.ListMSName = sorted(ListMSName)#[0:2] self.nMS = len(self.ListMSName) self.OutName = self.ListMSName[0].split("/")[-1].split("_")[0] self.ReadMSInfos() self.InitFromCatalog() - + elif self.Mode=="Plot": self.OutName = self.BaseDirSpecs.split("_")[-1] self.InitFromSpecs() @@ -150,10 +151,10 @@ def InitFromSpecs(self): self.NDirSelected=len(ListTargetFits) NOrig=len(ListTargetFits) - + ListOffFits=glob.glob("%s/OFF/*.fits"%self.BaseDirSpecs) NOff=len(ListOffFits) - + self.PosArray=np.zeros((self.NDirSelected+NOff,),dtype=[('Name','S200'),("ra",np.float64), ("dec",np.float64),('Type','S200'), ('iTessel',int),('iFacet',int)]) @@ -172,7 +173,7 @@ def InitFromSpecs(self): if ra<0.: ra+=2.*np.pi self.PosArray.ra[iDir]=ra dec=self.PosArray.dec[iDir]=float(F[0].header['DEC_RAD']) - + # print File,rad2hmsdms(ra,Type="ra").replace(" ",":"),rad2hmsdms(dec,Type="dec").replace(" ",":") # if self.PosArray.Type[iDir]=="Off": stop for iPol in range(4): @@ -195,7 +196,7 @@ def InitFromSpecs(self): continue self.PosArray.Type[iDir]=PosArrayTarget.Type[iS] self.PosArray.Name[iDir]=PosArrayTarget.Name[iS] - + def InitFromCatalog(self): FileCoords=self.FileCoords @@ -222,7 +223,7 @@ def InitFromCatalog(self): additional=np.array([additional],dtype=dtype) for r in additional: l.append(tuple(r)) - + self.PosArray=np.asarray(l,dtype=dtype) print("Created an array with %i records" % len(result), file=log) elif self.options.UseGaiaDB is not None: @@ -234,8 +235,8 @@ def InitFromCatalog(self): Parallax_min=1./(Dmax*1e-3) query=f"""SELECT TOP 10000 gaia_source.designation,gaia_source.source_id,gaia_source.ref_epoch,gaia_source.ra,gaia_source.dec,gaia_source.parallax, gaia_source.pmra,gaia_source.pmdec,gaia_source.phot_g_mean_mag,gaia_source.bp_rp - FROM gaiadr3.gaia_source - WHERE + FROM gaiadr3.gaia_source + WHERE CONTAINS( POINT('ICRS',gaiadr3.gaia_source.ra,gaiadr3.gaia_source.dec), CIRCLE('ICRS',{rac_deg},{decc_deg},{Radius_deg}) @@ -246,9 +247,9 @@ def InitFromCatalog(self): result = job.get_results() log.print(f"Query has returned {len(result)}") - - - + + + l=[] dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64), ("pmra",np.float64),("pmdec",np.float64),("ref_epoch",np.float64), @@ -266,7 +267,7 @@ def InitFromCatalog(self): b"Gaia DR3")) self.PosArray=np.asarray(l,dtype=dtype) - + CGGS=ClassGiveGaiaSample((rac_deg,decc_deg,Radius_deg),self.PosArray,RefCat="/data/cyril.tasse/Analyse_DataDynSpec_Jan23_TestRM/MergedCat.npz.MergeGaia.npz") indExo,indGaia=CGGS.buildRandGaiaSample() stop @@ -293,16 +294,16 @@ def InitFromCatalog(self): additional1["ra"][:]=additional["ra"][:] additional1["dec"][:]=additional["dec"][:] additional1["Type"][:]=additional["Type"][:] - + for r in additional1: l.append(tuple(r)) - + self.PosArray=np.asarray(l,dtype=dtype) log.print("Created an array with %i records" % len(l)) DoProperMotionCorr=True - + # import pylab # pylab.clf() # pylab.subplot(1,2,1) @@ -318,9 +319,9 @@ def InitFromCatalog(self): if self.PosArray.size>NMax: ind=np.int64(np.random.rand(NMax)*self.PosArray.size) self.PosArray=self.PosArray[ind] - + log.print("Created an array with %i records" % self.PosArray.size) - + elif self.options.FitsCatalog: print("Using the fits catalog: %s"%self.options.FitsCatalog, file=log) F=fits.open(self.options.FitsCatalog) @@ -345,16 +346,16 @@ def InitFromCatalog(self): additional1["ra"][:]=additional["ra"][:] additional1["dec"][:]=additional["dec"][:] additional1["Type"][:]=additional["Type"][:] - + for r in additional1: l.append(tuple(r)) - + self.PosArray=np.asarray(l,dtype=dtype) log.print("Created an array with %i records" % len(l)) DoProperMotionCorr=True else: - + #FileCoords="Transient_LOTTS.csv" if FileCoords is None: if not os.path.isfile(FileCoords): @@ -366,11 +367,11 @@ def InitFromCatalog(self): #self.PosArray=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",")[()] self.PosArray=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",") if len(self.PosArray.shape)==0: self.PosArray=self.PosArray.reshape((1,)) - + self.PosArray=self.PosArray.view(np.recarray) self.PosArray.ra*=np.pi/180. self.PosArray.dec*=np.pi/180. - + Radius=self.Radius NOrig=self.PosArray.Name.shape[0] Dist=AngDist(self.ra0,self.PosArray.ra,self.dec0,self.PosArray.dec) @@ -383,13 +384,13 @@ def InitFromCatalog(self): for iPCat,PCat in enumerate(self.PosArray): #print(iPCat,len(self.PosArray)) ra,dec=PCat["ra"],PCat["dec"] - + ra1,dec1=ProperMotionCorrection(PCat["ra"],PCat["dec"], PCat["pmra"],PCat["pmdec"],PCat["ref_epoch"],PCat["parallax"],self.tmin) if np.isnan(ra1) or np.isnan(dec1): ra1,dec1=ra,dec log.print(str((PCat["ra"],PCat["dec"],PCat["pmra"],PCat["pmdec"],PCat["ref_epoch"],PCat["parallax"],ra1,dec1))) - + Lra.append(ra1-ra) Ldec.append(dec1-dec) self.PosArray["ra"][iPCat]=ra1 @@ -397,7 +398,7 @@ def InitFromCatalog(self): print(np.array(Lra)*3600.*180/np.pi,np.array(Ldec)*3600.*180/np.pi) - + self.NDirSelected=self.PosArray.shape[0] print("Selected %i target [out of the %i in the original list]"%(self.NDirSelected,NOrig), file=log) @@ -405,9 +406,9 @@ def InitFromCatalog(self): print(ModColor.Str(" Have found no sources - returning"), file=log) self.killWorkers() return - + NOff=self.NOff - + if NOff==-1: NOff=self.PosArray.shape[0]*2 if NOff is not None: @@ -421,7 +422,7 @@ def InitFromCatalog(self): theta=np.linspace(0,2*np.pi,10) xC0=np.cos(theta)*Rrad0 yC0=np.sin(theta)*Rrad0 - + def give_iFacet_iTessel(l,m): #Plm=Polygon.Polygon(np.array([[l,m]])) Plm=Polygon.Polygon(np.array([xC0+l,yC0+m]).T) @@ -434,11 +435,11 @@ def give_iFacet_iTessel(l,m): ArrMax=Arr iFacetMax=iFacet return iFacetMax,self.DFacet[iFacetMax]["iSol"] - + self.PosArray=RecArrayOps.AppendField(self.PosArray,'iFacet',int) self.PosArray=RecArrayOps.AppendField(self.PosArray,'iTessel',int) l, m = self.CoordMachine.radec2lm(self.PosArray.ra, self.PosArray.dec) - + for iThis in range(self.PosArray.size): iFacet,iSol=give_iFacet_iTessel(l[iThis],m[iThis]) self.PosArray.iFacet[iThis]=iFacet @@ -451,7 +452,7 @@ def give_iFacet_iTessel(l,m): # pylab.scatter(self.PosArray.ra[ind],self.PosArray.dec[ind])#,c=iTessel) # pylab.draw() # pylab.show() - + self.NDir=self.PosArray.shape[0] print("For a total of %i targets"%(self.NDir), file=log) @@ -466,7 +467,7 @@ def give_iFacet_iTessel(l,m): self.DicoGrids["GridWeight2"] = np.zeros((self.NDir,self.NChan, self.NTimesGrid, 4), np.complex128) - + self.DoJonesCorr_kMS =False self.DicoJones=None @@ -486,14 +487,14 @@ def give_iFacet_iTessel(l,m): num_io_processes=1, verbose=0) self.APP=AsyncProcessPool.APP - + self.APP.registerJobHandlers(self) self.APP.startWorkers() - - - + + + def GiveOffPosArray(self,NOff): print("Making random off catalog with %i directions"%NOff, file=log) CatOff=np.zeros((NOff,),self.PosArray.dtype) @@ -512,13 +513,13 @@ def GiveOffPosArray(self,NOff): ind=np.where( (Fd.Isl_Total_flux>S0) & (Fd.Isl_Total_flux1: stop dtBin = dtBin_.flat[0] - - + + # if not np.allclose(ThisTimes, self.times): # raise ValueError("should have the same times") @@ -746,15 +747,15 @@ def ReadMSInfos(self): else: chSlice=slice(None) RevertChans=False - - + + if tmin is None: tmin=ThisTimes.min() tmax=ThisTimes.max() else: tmin=np.min([tmin,ThisTimes.min()]) tmax=np.max([tmax,ThisTimes.max()]) - + DicoMSInfos[iMS] = {"MSName": MSName, "ChanFreq": tf.getcol("CHAN_FREQ").ravel()[chSlice], # Hz "ChanWidth": np.abs(tf.getcol("CHAN_WIDTH").ravel()), # Hz @@ -767,21 +768,21 @@ def ReadMSInfos(self): "deltaTime": (ThisTimes[-1] - ThisTimes[0])/3600., # h "RevertChans": RevertChans, "Readable": True} - + if DicoMSInfos[iMS]["ChanWidth"][0] != self.ChanWidth: raise ValueError("should have the same chan width") pBAR.render(iMS+1, self.nMS) - + self.NTimesGrid=int(np.ceil((tmax-tmin)/dtBin)) self.timesGrid=tmin+np.arange(self.NTimesGrid)*dtBin self.tmin=tmin self.tmax=tmax - + for iMS in range(self.nMS): if not DicoMSInfos[iMS]["Readable"]: print(ModColor.Str("Problem reading %s"%MSName), file=log) print(ModColor.Str(" %s"%DicoMSInfos[iMS]["Exception"]), file=log) - + t.close() tf.close() @@ -795,12 +796,12 @@ def ReadMSInfos(self): #dt=np.median(dtArr) - + f0, f1 = self.Freq_minmax self.NChan = int((f1 - f0)/self.ChanWidth) + 1 # Fill properties - + self.tStart = Time(tmin/(24*3600.), format='mjd', scale='utc').isot self.tStop = Time(tmax/(24*3600.), format='mjd', scale='utc').isot self.fMin = self.Freq_minmax[0] @@ -825,24 +826,24 @@ def ReadMSInfos(self): LJob.append((iMS,iChunk)) self.LJob=LJob - + def LoadMS(self,iJob): iMS,iChunk=self.LJob[iJob] T0,T1=self.T0s[iChunk],self.T1s[iChunk] - - - if not self.DicoMSInfos[iMS]["Readable"]: + + + if not self.DicoMSInfos[iMS]["Readable"]: print("Skipping [%i/%i]: %s"%(iMS+1, self.nMS, self.ListMSName[iMS]), file=log) return "NotRead" print("Reading [%i/%i] %s:%s"%(iMS+1, self.nMS, self.ListMSName[iMS],self.ColName), file=log) MSName=self.ListMSName[iMS] - + t = table(MSName, ack=False) times = t.getcol("TIME") t0=self.tmin - + ind=np.where((times>=(T0+t0))&(times<(T1+t0)))[0] if ind.size==0: print("No Data in requested interval %f -> %f h time interval"%(T0/3600,T1/3600), file=log) @@ -852,17 +853,17 @@ def LoadMS(self,iJob): if ROW0!=0 or NROW!=t.nrows(): print(" Reading chunk in %.3f -> %.3f h"%(T0/3600,T1/3600), file=log) - + nch = self.DicoMSInfos[iMS]["ChanFreq"].size npol = self.DicoMSInfos[iMS]["npol"] #chSlice=self.DicoMSInfos[iMS]["chSlice"] RevertChans=self.DicoMSInfos[iMS]["RevertChans"] - + data = np.zeros((NROW,nch,npol),np.complex64) t.getcolnp(self.ColName,data,ROW0,NROW) if RevertChans: data=data[:,::-1,:] - + if self.ModelName: print(" Substracting %s from %s"%(self.ModelName,self.ColName), file=log) model=np.zeros((NROW,nch,npol),np.complex64) @@ -872,8 +873,8 @@ def LoadMS(self,iJob): data-=model del(model) - - + + if self.ColWeights: print(" Reading weight column %s"%(self.ColWeights), file=log) sW=t.getcol(self.ColWeights,0,1).shape @@ -884,8 +885,8 @@ def LoadMS(self,iJob): else: weights=np.zeros((NROW,nch),np.float32) t.getcolnp(self.ColWeights,weights,ROW0,NROW) - - + + if RevertChans: weights=weights[:,::-1] else: nrow,nch,_=data.shape @@ -894,13 +895,13 @@ def LoadMS(self,iJob): flag=np.zeros((NROW,nch,npol),np.bool) t.getcolnp("FLAG",flag,ROW0,NROW) if RevertChans: flag=flag[:,::-1] - + # data[:,:,:]=0 # data[:,0:300,0]=1 # flag.fill(0) # weights.fill(1) - + times = t.getcol("TIME",ROW0,NROW) A0, A1 = t.getcol("ANTENNA1",ROW0,NROW), t.getcol("ANTENNA2",ROW0,NROW) @@ -935,7 +936,7 @@ def LoadMS(self,iJob): DicoDATA["w"]=w0 DicoDATA["uniq_times"]=np.unique(DicoDATA["times"]) - + if self.DoJonesCorr_kMS or self.DoJonesCorr_Beam: self.setJones(DicoDATA) @@ -980,7 +981,7 @@ def setJones(self,DicoDATA): c0,c1=CutGainsMinMax G[G>c1]=0 G[G0)[0] #f = DicoDATA["flag"][indRow, :, :] #d = DicoDATA["data"][indRow, :, :] @@ -1189,18 +1190,18 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): d = np.array((DicoDATA["data"].flat[indArr.flat[:]]).reshape((nRowOut,nch,npol))).copy() f = np.array((DicoDATA["flag"].flat[indArr.flat[:]]).reshape((nRowOut,nch,npol))).copy() T.timeit("first") - + # for i in range(10): # d = (DicoDATA["data"].flat[indArr.flat[:]]).reshape((nRowOut,nch,npol)).copy() # f = (DicoDATA["flag"].flat[indArr.flat[:]]).reshape((nRowOut,nch,npol)).copy() # T.timeit("first %i"%i) - + nrow,nch,_=d.shape #weights = (DicoDATA["weights"][indRow, :]).reshape((nrow,nch,1)) - + indArr=nch*np.int64(indR)+np.int64(indCh) weights = np.array((DicoDATA["weights"].flat[indArr.flat[:]]).reshape((nRowOut,nch,1))).copy() - + A0s = DicoDATA["A0"][indRow].copy() A1s = DicoDATA["A1"][indRow].copy() u0 = DicoDATA["u"][indRow].reshape((-1,1,1)).copy() @@ -1214,14 +1215,14 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): # kk = np.exp( -2.*np.pi*1j* f/const.c.value *(u0*l + v0*m + w0*(n-1)) ) # Phasing term #print iTime,iDir ChanFreqs=np.array(self.DicoMSInfos[iMS]["ChanFreq"][0]).copy() - + iTimeGrid=np.argmin(np.abs(self.timesGrid-self.DicoMSInfos[iMS]["times"][iTime])) - + dcorr=d.copy() f0, _ = self.Freq_minmax ich0 = int( (ChanFreqs - f0)/self.ChanWidth ) OneMinusF=(1-f).copy() - + W=np.zeros((nRowOut,nch,npol),np.float32) for ipol in range(npol): W[:,:,ipol]=weights[:,:,0] @@ -1230,7 +1231,7 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): # weights=weights*np.ones((1,1,npol)) # W=weights - + kk=np.zeros_like(d) T.timeit("third") for iDir in range(self.NDir): @@ -1240,7 +1241,7 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): l, m = self.radec2lm(ra, dec,ra0,dec0) n = np.sqrt(1. - l**2. - m**2.) - + T.timeit("lmn") kkk = np.exp(-2.*np.pi*1j* chfreq/const.c.value *(u0*l + v0*m + w0*(n-1)) ) # Phasing term T.timeit("kkk") @@ -1258,21 +1259,21 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): # pylab.draw() # pylab.show(False) # pylab.pause(0.1) - - - + + + #DicoMSInfos = self.DicoMSInfos - + #_,nch,_=DicoDATA["data"].shape - + dcorr[:]=d[:] W=Wc.copy() dcorr*=W wdcorr=np.ones(dcorr.shape,np.float64) #kk=kk*np.ones((1,1,npol)) - + T.timeit("corr") - + if self.DoJonesCorr_kMS or self.DoJonesCorr_Beam: T1=ClassTimeIt.ClassTimeIt(" DoJonesCorr") T1.disable() @@ -1282,19 +1283,19 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): tm = DicoJones['tm'] # Time slot for the solution iTJones=np.argmin(np.abs(tm-ThisTime))#self.timesGrid[iTime])) - + #iDJones=np.argmin(AngDist(ra,DicoJones['ra'],dec,DicoJones['dec'])) lJones, mJones = self.CoordMachine.radec2lm(DicoJones['ra'], DicoJones['dec']) iDJones=np.argmin(np.sqrt((l-lJones)**2+(m-mJones)**2)) - + _,nchJones,_,_,_,_=DicoJones['G'].shape T1.timeit("argmin") - - + + for iFJones in range(nchJones): - + nu0,nu1=DicoJones['FreqDomains'][iFJones] fData=self.DicoMSInfos[iMS]["ChanFreq"].ravel() indCh=np.where((fData>=nu0) & (fData>log,"Making pdf overview: %s"%pdfname - # pBAR = ProgressBar(Title="Making pages") + # pBAR = ProgressBar(Title="Making pages") # NPages=self.DynSpecMS.NDirSelected #Selected # iDone=0 # pBAR.render(0, NPages) @@ -270,7 +270,7 @@ def PlotSpec(self,Prefix=""): # iDone+=1 # pBAR.render(iDone, NPages) - + def PlotSpecSingleDir(self, iDir=0, BoxArcSec=300.): label = ["I", "Q", "U", "V"] @@ -279,7 +279,7 @@ def PlotSpecSingleDir(self, iDir=0, BoxArcSec=300.): pylab.rc('text', usetex=True) font = {'family':'serif', 'serif': ['Times']} pylab.rc('font', **font) - + # Figure properties bigfont = 8 @@ -320,7 +320,7 @@ def PlotSpecSingleDir(self, iDir=0, BoxArcSec=300.): cmap = copy.copy(matplotlib.cm.get_cmap(cmap))#.copy() cmap.set_bad(color='black') - + ax1 = pylab.subplot(4, 1, ipol+1) spec = pylab.pcolormesh(times, freqs, Gn, cmap=cmap,#'bone_r', vmin=mean-3*sig, vmax=mean+10*sig, rasterized=True,shading='auto') @@ -343,11 +343,11 @@ def PlotSpecSingleDir(self, iDir=0, BoxArcSec=300.): sig = GiveMAD(Gn) mean = np.median(Gn) #spec = pylab.pcolormesh(times, freqs, Gn, cmap='bone_r', vmin=mean-3*sig, vmax=mean+10*sig, rasterized=True) - spec = pylab.imshow(Gn, interpolation="nearest", cmap='bone_r', vmin=mean-3*sig, vmax=mean+10*sig, extent=(times[0],times[-1],self.DynSpecMS.fMin*1.e-6,self.DynSpecMS.fMax*1.e-6),rasterized=True) + spec = pylab.imshow(Gn, interpolation="nearest", cmap='bone_r', vmin=mean-3*sig, vmax=mean+10*sig, extent=(times[0],times[-1],self.DynSpecMS.fMin*1.e-6,self.DynSpecMS.fMax*1.e-6),rasterized=True) axspec.axis('tight') cbar = pylab.colorbar(fraction=0.046, pad=0.01) cbar.ax.tick_params(labelsize=smallfont) - cbar.set_label(r'Flux density (Jy)', fontsize=8, horizontalalignment='center') + cbar.set_label(r'Flux density (Jy)', fontsize=8, horizontalalignment='center') pylab.text(times[-1]-0.02*(times[-1]-times[0]), freqs[-1]-0.1*(freqs[-1]-freqs[0]), 'I', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) pylab.text(times[0]+0.02*(times[-1]-times[0]), freqs[0]+0.1*(freqs[-1]-freqs[0]), r"$\sigma =$ %.3f Jy"%sig, horizontalalignment='left', verticalalignment='center', fontsize=bigfont+2) pylab.xlabel("Time (min since %s)"%(t0.iso), fontsize=bigfont) @@ -360,13 +360,13 @@ def PlotSpecSingleDir(self, iDir=0, BoxArcSec=300.): AG=np.abs(Gn) sig = GiveMAD(Gn) - mean = np.median(Gn) + mean = np.median(Gn) #spec = pylab.pcolormesh(times, freqs, Gn, cmap='bone_r', vmin=0, vmax=mean+10*sig, rasterized=True) - spec = pylab.imshow(Gn, interpolation="nearest", cmap='bone_r', vmin=mean-3*sig, vmax=mean+10*sig, extent=(times[0],times[-1],self.DynSpecMS.fMin*1.e-6,self.DynSpecMS.fMax*1.e-6), rasterized=True) + spec = pylab.imshow(Gn, interpolation="nearest", cmap='bone_r', vmin=mean-3*sig, vmax=mean+10*sig, extent=(times[0],times[-1],self.DynSpecMS.fMin*1.e-6,self.DynSpecMS.fMax*1.e-6), rasterized=True) axspec.axis('tight') cbar = pylab.colorbar(fraction=0.046, pad=0.01) cbar.ax.tick_params(labelsize=smallfont) - cbar.set_label(r'Flux density (Jy)', fontsize=8, horizontalalignment='center') + cbar.set_label(r'Flux density (Jy)', fontsize=8, horizontalalignment='center') pylab.text(times[-1]-0.02*(times[-1]-times[0]), freqs[-1]-0.1*(freqs[-1]-freqs[0]), 'L', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) pylab.text(times[0]+0.02*(times[-1]-times[0]), freqs[0]+0.1*(freqs[-1]-freqs[0]), r"$\sigma =$ %.3f Jy"%sig, horizontalalignment='left', verticalalignment='center', fontsize=bigfont+2) pylab.xlabel("Time (min since %s)"%(t0.iso), fontsize=bigfont) @@ -379,13 +379,13 @@ def PlotSpecSingleDir(self, iDir=0, BoxArcSec=300.): AG=np.abs(Gn) sig = GiveMAD(Gn) - mean = np.median(Gn) + mean = np.median(Gn) #spec = pylab.pcolormesh(times, freqs, Gn, cmap='bone_r', vmin=mean-3*sig, vmax=mean+10*sig, rasterized=True) - spec = pylab.imshow(Gn, interpolation="nearest", cmap='bone_r', vmin=mean-5*sig, vmax=mean+5*sig, extent=(times[0],times[-1],self.DynSpecMS.fMin*1.e-6,self.DynSpecMS.fMax*1.e-6), rasterized=True) + spec = pylab.imshow(Gn, interpolation="nearest", cmap='bone_r', vmin=mean-5*sig, vmax=mean+5*sig, extent=(times[0],times[-1],self.DynSpecMS.fMin*1.e-6,self.DynSpecMS.fMax*1.e-6), rasterized=True) axspec.axis('tight') cbar = pylab.colorbar(fraction=0.046, pad=0.01) cbar.ax.tick_params(labelsize=smallfont) - cbar.set_label(r'Flux density (Jy)', fontsize=8, horizontalalignment='center') + cbar.set_label(r'Flux density (Jy)', fontsize=8, horizontalalignment='center') pylab.text(times[-1]-0.02*(times[-1]-times[0]), freqs[-1]-0.1*(freqs[-1]-freqs[0]), 'V', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) pylab.text(times[0]+0.02*(times[-1]-times[0]), freqs[0]+0.1*(freqs[-1]-freqs[0]), r"$\sigma =$ %.3f Jy"%sig, horizontalalignment='left', verticalalignment='center', fontsize=bigfont+2) pylab.xlabel("Time (min since %s)"%(t0.iso), fontsize=bigfont) @@ -451,7 +451,7 @@ def giveBounded(x): x1=giveBounded(xc+box) y0=giveBounded(yc-box) y1=giveBounded(yc+box) - + DataBoxed=self.ImageIData[y0:y1,x0:x1] FluxI=self.ImageIData[yc,xc] sigFluxI=GiveMAD(DataBoxed) @@ -467,7 +467,7 @@ def giveBounded(x): im = pylab.imshow(DataBoxed, interpolation="nearest", cmap='bone_r', aspect="auto", vmin=vMin, vmax=vMax, origin='lower', rasterized=True) #pylab.text((ra_crop[1]-ra_crop[0])/16, (dec_crop[1]-dec_crop[0])/16, r"$\sigma =$ %.3f mJy"%rms, horizontalalignment='left', verticalalignment='center', fontsize=bigfont+2) cbar = pylab.colorbar()#(fraction=0.046*2., pad=0.01*4.) - + ax1.set_xlabel(r'RA (J2000)') raax = ax1.coords[0] raax.set_major_formatter('hh:mm:ss') @@ -483,13 +483,13 @@ def giveBounded(x): cbar.ax.tick_params(labelsize=smallfont) pylab.setp(ax1.get_xticklabels(), rotation='horizontal', fontsize=smallfont) pylab.setp(ax1.get_yticklabels(), rotation='horizontal', fontsize=smallfont) - + ra_cen, dec_cen = wcs.wcs_world2pix(np.degrees(self.DynSpecMS.PosArray.ra[iDir]), np.degrees(self.DynSpecMS.PosArray.dec[iDir]), 1) #pylab.plot(ra_cen, dec_cen, 'o', markerfacecolor='none', markeredgecolor='red', markersize=bigfont) # plot a circle at the target #pylab.plot(newra_cen, newdec_cen, 'o', markerfacecolor='none', markeredgecolor='red', markersize=bigfont) # plot a circle at the target pylab.plot(DataBoxed.shape[1]/2., DataBoxed.shape[0]/2., 'o', markerfacecolor='none', markeredgecolor='red', markersize=bigfont) # plot a circle at the target - - pylab.text(DataBoxed.shape[0]*0.9, DataBoxed.shape[1]*0.9, 'I', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) + + pylab.text(DataBoxed.shape[0]*0.9, DataBoxed.shape[1]*0.9, 'I', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) # ---- Image V ---- @@ -524,7 +524,7 @@ def giveBounded(x): self.CatFlux.Type[iDir]=self.DynSpecMS.PosArray.Type[iDir] self.CatFlux.ra[iDir]=self.DynSpecMS.PosArray.ra[iDir] self.CatFlux.dec[iDir]=self.DynSpecMS.PosArray.dec[iDir] - + newra_cen, newdec_cen = wcs.wcs_pix2world( (x1+x0)/2., (y1+y0)/2., 1) wcs.wcs.crpix = [ DataBoxed.shape[1]/2., DataBoxed.shape[0]/2. ] # update the WCS object wcs.wcs.crval = [ newra_cen, newdec_cen ] @@ -549,7 +549,7 @@ def giveBounded(x): pylab.setp(ax1.get_yticklabels(), rotation='horizontal', fontsize=smallfont) ra_cen, dec_cen = wcs.wcs_world2pix(np.degrees(self.DynSpecMS.PosArray.ra[iDir]), np.degrees(self.DynSpecMS.PosArray.dec[iDir]), 1) pylab.plot(DataBoxed.shape[1]/2., DataBoxed.shape[0]/2., 'o', markerfacecolor='none', markeredgecolor='red', markersize=bigfont) # plot a circle at the target - pylab.text(DataBoxed.shape[0]*0.9, DataBoxed.shape[1]*0.9, 'V', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) + pylab.text(DataBoxed.shape[0]*0.9, DataBoxed.shape[1]*0.9, 'V', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) #pylab.subplots_adjust(wspace=0.15, hspace=0.30) diff --git a/ms2dynspec.py b/ms2dynspec.py index 106468a..e63f84c 100755 --- a/ms2dynspec.py +++ b/ms2dynspec.py @@ -7,7 +7,7 @@ from past.builtins import cmp __author__ = "Cyril Tasse, and Alan Loh" __credits__ = ["Cyril Tasse", "Alan Loh"] -from DynSpecMS.dynspecms_version import version +from dynspecms_version import version from DDFacet.Other import logger log=logger.getLogger("ms2dynspec") from DDFacet.Other import ModColor @@ -18,8 +18,8 @@ ========================================================================= DESCRIPTION Blablabla - Modif - + Modif + Example: python ms2dynspec.py --msfile /usr/data/ --data CORRECTED --model PREDICT --sols killms.npz --srclist SRCPOS.txt --rad 10 @@ -52,10 +52,10 @@ import glob, os import pylab from DDFacet.Other import MyPickle -from DynSpecMS import logo -logo.PrintLogo(__version__) -from DynSpecMS.ClassDynSpecMS import ClassDynSpecMS -from DynSpecMS import ClassSaveResults +#from ClassDynSpecMS import logo +#logo.PrintLogo(__version__) +from ClassDynSpecMS import ClassDynSpecMS +from ClassSaveResults import ClassSaveResults from DDFacet.Data.ClassMS import expandMSList from DDFacet.Other import ModColor from DDFacet.Other import progressbar @@ -110,8 +110,8 @@ def main(args=None, messages=[]): for k in DT.keys(): MSList=DT[k] - D = ClassDynSpecMS(ListMSName=MSList, - ColName=args.data, ModelName=args.model, + D = ClassDynSpecMS(ListMSName=MSList, + ColName=args.data, ModelName=args.model, SolsName=args.sols, TChunkHours=args.TChunkHours, ColWeights=args.WeightCol, @@ -133,9 +133,9 @@ def main(args=None, messages=[]): if D.NDirSelected==0: return - + if D.Mode=="Spec": D.StackAll() - + SaveMachine=ClassSaveResults.ClassSaveResults(D,DIRNAME=args.OutDirName) if D.Mode=="Spec": SaveMachine.WriteFits() @@ -148,7 +148,7 @@ def main(args=None, messages=[]): SaveMachine.PlotSpec(Prefix="_replot") D.killWorkers() Multiprocessing.cleanupShm() - + # ========================================================================= # ========================================================================= @@ -157,7 +157,7 @@ def main(args=None, messages=[]): parser.add_argument("--ms", type=str, help="Name of MS file / directory", required=False) parser.add_argument("--data", type=str, default="CORRECTED", help="Name of DATA column", required=False) parser.add_argument("--TChunkHours", type=float, default=0., help="Chunk size in hours", required=False) - + parser.add_argument("--WeightCol", type=str, default=None, help="Name of weights column to be taken into account", required=False) parser.add_argument("--model", type=str, help="Name of MODEL column",default="")#, required=True) parser.add_argument("--sols", type=str, help="Jones solutions",default="") @@ -178,7 +178,7 @@ def main(args=None, messages=[]): parser.add_argument("--UseLoTSSDB", type=int, default=0, help="Use LoTSS DB for target list", required=False) parser.add_argument("--UseGaiaDB", type=str, default=None, help="Use Gaia DB for target list", required=False) parser.add_argument("--DoTar", type=int, default=1, help="Tar final products", required=False) - + parser.add_argument("--NCPU", type=int, default=0, help="NCPU", required=False) parser.add_argument("--BeamModel", type=str, default=None, help="Beam Model to be used", required=False) parser.add_argument("--BeamNBand", type=int, default=1, help="Number of channels in the Beam Jones matrix", required=False) From 47222e5a0a8e76b7da9dcd99235c8f0d129dda63 Mon Sep 17 00:00:00 2001 From: Buntu Ngcebetsha Date: Mon, 22 Apr 2024 17:51:00 +0200 Subject: [PATCH 2/7] Save current progress The script is almost there, I keep getting bugs. I am presently experiencing: ``` File "/data4/sirothia/proxima_btry/msdir/DynSpecMS/ClassDynSpecMS.py", line 939, in StackAll self.processJob(iJob) File "/data4/sirothia/proxima_btry/msdir/DynSpecMS/ClassDynSpecMS.py", line 987, in processJob self.delShm(iJob) File "/data4/sirothia/proxima_btry/msdir/DynSpecMS/ClassDynSpecMS.py", line 990, in delShm shared_dict.SharedDict.delete_item("DATA_%i"%(iJob)) TypeError: delete_item() missing 1 required positional argument: 'item' ``` --- ClassDynSpecMS.py | 243 ++++++++++---------------------------------- ClassSaveResults.py | 8 +- ms2dynspec.py | 124 ++++++++++++++-------- 3 files changed, 144 insertions(+), 231 deletions(-) diff --git a/ClassDynSpecMS.py b/ClassDynSpecMS.py index b523a2c..0fcdbe2 100644 --- a/ClassDynSpecMS.py +++ b/ClassDynSpecMS.py @@ -21,7 +21,6 @@ import os from killMS.Other import reformat from DDFacet.Other import AsyncProcessPool -#from .dynspecms_version import version from dynspecms_version import version import glob from astropy.io import fits @@ -36,6 +35,7 @@ #from killMS.Data import ClassJonesDomains import DDFacet.Other.ClassJonesDomains import psutil +import ClassGiveCatalog def AngDist(ra0,ra1,dec0,dec1): AC=np.arccos @@ -83,8 +83,11 @@ def __init__(self, SourceCatOff=None, SourceCatOff_FluxMean=None, SourceCatOff_dFluxMean=None, - options=None): + options=None,SubSet=None): + self.options=options + self.SubSet=SubSet + # CutGainsMinMax if BeamModel=="None": @@ -200,183 +203,14 @@ def InitFromSpecs(self): def InitFromCatalog(self): FileCoords=self.FileCoords - dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64),('Type','S200')] # should we use the surveys DB? - DoProperMotionCorr=False - if self.options.UseLoTSSDB: - print("Using the surveys database", file=log) - from surveys_db import SurveysDB - with SurveysDB() as sdb: - sdb.cur.execute('select * from transients') - result=sdb.cur.fetchall() - # convert to a list, then to ndarray, then to recarray - l=[] - for r in result: - l.append((r['id'],r['ra'],r['decl'],r['type'])) - - if FileCoords is not None and FileCoords!="": - print('Adding data from file '+FileCoords, file=log) - additional=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",")[()] - if len(additional.shape)==0: additional=additional.reshape((1,)) - if not additional.shape: - # deal with a one-line input file - additional=np.array([additional],dtype=dtype) - for r in additional: - l.append(tuple(r)) - - self.PosArray=np.asarray(l,dtype=dtype) - print("Created an array with %i records" % len(result), file=log) - elif self.options.UseGaiaDB is not None: - from astroquery.gaia import Gaia - rac_deg,decc_deg=self.ra0*180/np.pi, self.dec0*180/np.pi - Radius_deg=self.Radius - Dmax,NMax=self.options.UseGaiaDB.split(",") - Dmax,NMax=float(Dmax),int(NMax) - Parallax_min=1./(Dmax*1e-3) - query=f"""SELECT TOP 10000 gaia_source.designation,gaia_source.source_id,gaia_source.ref_epoch,gaia_source.ra,gaia_source.dec,gaia_source.parallax, - gaia_source.pmra,gaia_source.pmdec,gaia_source.phot_g_mean_mag,gaia_source.bp_rp - FROM gaiadr3.gaia_source - WHERE - CONTAINS( - POINT('ICRS',gaiadr3.gaia_source.ra,gaiadr3.gaia_source.dec), - CIRCLE('ICRS',{rac_deg},{decc_deg},{Radius_deg}) - )=1 AND (gaiadr3.gaia_source.parallax>={Parallax_min})""" - log.print(f"Sending du Gaia server") - print(f"{query}") - job = Gaia.launch_job(query) - result = job.get_results() - - log.print(f"Query has returned {len(result)}") - - - - l=[] - dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64), - ("pmra",np.float64),("pmdec",np.float64),("ref_epoch",np.float64), - ("parallax",np.float64),("GaiaDistance",np.float64), - ("g",np.float64),("G",np.float64),("b_r",np.float64), - ('Type','S200')] - for r in result: - parallax=r['parallax'] - g=r['g'] - D=1./(parallax*1e-3) - G=g+5*np.log10(parallax)-10 - l.append((r['DESIGNATION'],r['ra'],r['dec'],r['pmra'],r['pmdec'],r['ref_epoch'], - r['parallax'],D, - g,G,r['b_r'], - b"Gaia DR3")) - self.PosArray=np.asarray(l,dtype=dtype) - - - CGGS=ClassGiveGaiaSample((rac_deg,decc_deg,Radius_deg),self.PosArray,RefCat="/data/cyril.tasse/Analyse_DataDynSpec_Jan23_TestRM/MergedCat.npz.MergeGaia.npz") - indExo,indGaia=CGGS.buildRandGaiaSample() - stop - # ################################## - F=fits.open(self.options.FitsCatalog) - d=F[1].data - l=[] - dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64), - ("pmra",np.float64),("pmdec",np.float64),("ref_epoch",np.float64),("parallax",np.float64), - ('Type','S200')] - for r in d: - l.append((r['DESIGNATION'],r['ra'],r['dec'],r['pmra'],r['pmdec'],r['ref_epoch'],r['parallax'],r['Type'])) - - if FileCoords is not None and FileCoords!="": - print('Adding data from file '+FileCoords, file=log) - dtype0=[('Name','S200'),("ra",np.float64),("dec",np.float64),('Type','S200')] - additional=np.genfromtxt(FileCoords,dtype=dtype0,delimiter=",")[()] - if len(additional.shape)==0: additional=additional.reshape((1,)) - if not additional.shape: - # deal with a one-line input file - additional=np.array([additional],dtype=dtype0) - additional1=np.zeros((additional.shape[0],),dtype=dtype) - additional1["Name"][:]=additional["Name"][:] - additional1["ra"][:]=additional["ra"][:] - additional1["dec"][:]=additional["dec"][:] - additional1["Type"][:]=additional["Type"][:] - - for r in additional1: - l.append(tuple(r)) - - - self.PosArray=np.asarray(l,dtype=dtype) - log.print("Created an array with %i records" % len(l)) - DoProperMotionCorr=True - - - # import pylab - # pylab.clf() - # pylab.subplot(1,2,1) - # pylab.scatter(self.PosArray["ra"],self.PosArray["dec"]) - # pylab.scatter([rac_deg],[decc_deg]) - # pylab.subplot(1,2,2) - # D=1./(self.PosArray["parallax"]*1e-3) - # pylab.hist(D,bins=100) - # pylab.draw() - # pylab.show() - DoProperMotionCorr=True - if self.PosArray.size>NMax: - ind=np.int64(np.random.rand(NMax)*self.PosArray.size) - self.PosArray=self.PosArray[ind] - - log.print("Created an array with %i records" % self.PosArray.size) - - elif self.options.FitsCatalog: - print("Using the fits catalog: %s"%self.options.FitsCatalog, file=log) - F=fits.open(self.options.FitsCatalog) - d=F[1].data - l=[] - dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64), - ("pmra",np.float64),("pmdec",np.float64),("ref_epoch",np.float64),("parallax",np.float64), - ('Type','S200')] - for r in d: - l.append((r['DESIGNATION'],r['ra'],r['dec'],r['pmra'],r['pmdec'],r['ref_epoch'],r['parallax'],r['Type'])) - - if FileCoords is not None and FileCoords!="": - print('Adding data from file '+FileCoords, file=log) - dtype0=[('Name','S200'),("ra",np.float64),("dec",np.float64),('Type','S200')] - additional=np.genfromtxt(FileCoords,dtype=dtype0,delimiter=",")[()] - if len(additional.shape)==0: additional=additional.reshape((1,)) - if not additional.shape: - # deal with a one-line input file - additional=np.array([additional],dtype=dtype0) - additional1=np.zeros((additional.shape[0],),dtype=dtype) - additional1["Name"][:]=additional["Name"][:] - additional1["ra"][:]=additional["ra"][:] - additional1["dec"][:]=additional["dec"][:] - additional1["Type"][:]=additional["Type"][:] - - for r in additional1: - l.append(tuple(r)) - - - self.PosArray=np.asarray(l,dtype=dtype) - log.print("Created an array with %i records" % len(l)) - DoProperMotionCorr=True - else: - - #FileCoords="Transient_LOTTS.csv" - if FileCoords is None: - if not os.path.isfile(FileCoords): - ssExec="wget -q --user=anonymous ftp://ftp.strw.leidenuniv.nl/pub/tasse/%s -O %s"%(FileCoords,FileCoords) - print("Downloading %s"%FileCoords, file=log) - print(" Executing: %s"%ssExec, file=log) - os.system(ssExec) - log.print("Reading cvs file: %s"%FileCoords) - #self.PosArray=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",")[()] - self.PosArray=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",") - if len(self.PosArray.shape)==0: self.PosArray=self.PosArray.reshape((1,)) - - self.PosArray=self.PosArray.view(np.recarray) - self.PosArray.ra*=np.pi/180. - self.PosArray.dec*=np.pi/180. - - Radius=self.Radius - NOrig=self.PosArray.Name.shape[0] - Dist=AngDist(self.ra0,self.PosArray.ra,self.dec0,self.PosArray.dec) - ind=np.where(Dist<(Radius*np.pi/180))[0] - self.PosArray=self.PosArray[ind] + CGC=ClassGiveCatalog.ClassGiveCatalog(self.options, + self.ra0,self.dec0, + self.Radius, + self.FileCoords) + self.PosArray=CGC.giveCat(SubSet=self.SubSet) + DoProperMotionCorr=CGC.DoProperMotionCorr if DoProperMotionCorr: Lra,Ldec=[],[] @@ -396,12 +230,12 @@ def InitFromCatalog(self): self.PosArray["ra"][iPCat]=ra1 self.PosArray["dec"][iPCat]=dec1 - print(np.array(Lra)*3600.*180/np.pi,np.array(Ldec)*3600.*180/np.pi) + # print(np.array(Lra)*3600.*180/np.pi,np.array(Ldec)*3600.*180/np.pi) self.NDirSelected=self.PosArray.shape[0] - print("Selected %i target [out of the %i in the original list]"%(self.NDirSelected,NOrig), file=log) + print("Selected %i target [out of the %i in the original list]"%(self.NDirSelected,CGC.NOrig), file=log) if self.NDirSelected==0: print(ModColor.Str(" Have found no sources - returning"), file=log) self.killWorkers() @@ -478,7 +312,7 @@ def give_iFacet_iTessel(l,m): if self.BeamModel: self.DoJonesCorr_Beam=True - AsyncProcessPool.APP=None + #AsyncProcessPool.APP=None # AsyncProcessPool.init(ncpu=self.NCPU, # num_io_processes=1, # affinity="disable") @@ -724,7 +558,7 @@ def ReadMSInfos(self): tf = table("%s::SPECTRAL_WINDOW"%MSName, ack=False) ThisTimes = np.unique(t.getcol("TIME")) dtBin_=np.unique(t.getcol("INTERVAL")) - if dtBin_.size>1: stop + #if dtBin_.size>1: stop dtBin = dtBin_.flat[0] @@ -954,7 +788,7 @@ def setJones(self,DicoDATA): SolsName=SolsName.split(",") GD={"Beam":{"Model":self.BeamModel, "PhasedArrayMode":"A", - "At":"facet", + "At":"tessel", "DtBeamMin":5., "NBand":self.BeamNBand, "CenterNorm":1}, @@ -967,10 +801,27 @@ def setJones(self,DicoDATA): } print("Reading Jones matrices solution file:", file=log) + RADEC_ForcedBeamDirs=None + if not self.DoJonesCorr_kMS and self.DoJonesCorr_Beam: + RADEC_ForcedBeamDirs=(self.PosArray.ra,self.PosArray.dec) + + ms=ClassMS.ClassMS(self.DicoMSInfos[iMS]["MSName"],GD=GD,DoReadData=False,) JonesMachine = ClassJones.ClassJones(GD, ms, CacheMode=False) - JonesMachine.InitDDESols(DicoDATA) + JonesMachine.InitDDESols(DicoDATA,RADEC_ForcedBeamDirs=RADEC_ForcedBeamDirs) + + # print("================") + # RAJoneskillMS=DicoDATA["killMS"]["Dirs"]["ra"] + # DECJoneskillMS=DicoDATA["killMS"]["Dirs"]["dec"] + # for ThisRA,ThisDEC in zip(RAJoneskillMS,DECJoneskillMS): + # print(rad2hmsdms(ThisRA,Type="ra").replace(" ",":"),rad2hmsdms(ThisDEC,Type="dec").replace(" ",":")) + # print("================") + # RAJonesBeam=DicoDATA["Beam"]["Dirs"]["ra"] + # DECJonesBeam=DicoDATA["Beam"]["Dirs"]["dec"] + # for ThisRA,ThisDEC in zip(RAJonesBeam,DECJonesBeam): + # print(rad2hmsdms(ThisRA,Type="ra").replace(" ",":"),rad2hmsdms(ThisDEC,Type="dec").replace(" ",":")) + # print("================") CutGainsMinMax=StrToList(self.options.CutGainsMinMax) @@ -985,31 +836,40 @@ def setJones(self,DicoDATA): if self.DoJonesCorr_kMS and self.DoJonesCorr_Beam: DomainMachine=DDFacet.Other.ClassJonesDomains.ClassJonesDomains() JonesSols=DomainMachine.MergeJones(DicoDATA["killMS"]["Jones"], DicoDATA["Beam"]["Jones"]) + RAJones=JonesMachine.ClusterCat['ra'] + DECJones=JonesMachine.ClusterCat['dec'] elif self.DoJonesCorr_kMS: JonesSols=DicoDATA["killMS"]["Jones"] + RAJones=DicoDATA["killMS"]["Dirs"]["ra"] + DECJones=DicoDATA["killMS"]["Dirs"]["dec"] elif self.DoJonesCorr_Beam: JonesSols=DicoDATA["Beam"]["Jones"] + RAJones=DicoDATA["Beam"]["Dirs"]["ra"] + DECJones=DicoDATA["Beam"]["Dirs"]["dec"] else: stop + + DicoJones=shared_dict.create("DicoJones_%i"%iJob) DicoJones["G"]=np.swapaxes(JonesSols["Jones"],1,3) # Normalize Jones matrices G=DicoJones["G"] nt,nch,na,nDir,_,_=G.shape DicoJones['tm']=(JonesSols["t0"]+JonesSols["t1"])/2. - DicoJones['ra']=JonesMachine.ClusterCat['ra'] - DicoJones['dec']=JonesMachine.ClusterCat['dec'] + DicoJones['ra']=RAJones + DicoJones['dec']=DECJones DicoJones['FreqDomains']=JonesSols['FreqDomains'] DicoJones['FreqDomains_mean']=np.mean(JonesSols['FreqDomains'],axis=1) DicoJones['IDJones']=np.zeros((self.NDir,),np.int32) + lJones, mJones = self.CoordMachine.radec2lm(DicoJones['ra'], DicoJones['dec']) for iDir in range(self.NDir): ra=self.PosArray.ra[iDir] dec=self.PosArray.dec[iDir] #DicoJones['IDJones'][iDir]=np.argmin(AngDist(ra,DicoJones['ra'],dec,DicoJones['dec'])) lTarget, mTarget = self.CoordMachine.radec2lm(np.array([ra]), np.array([dec])) - lJones, mJones = self.CoordMachine.radec2lm(DicoJones['ra'], DicoJones['dec']) DicoJones['IDJones'][iDir]=np.argmin(np.sqrt((lTarget-lJones)**2+(mTarget-mJones)**2)) + # print(DicoJones['IDJones']) @@ -1127,8 +987,10 @@ def processJob(self,iJob): self.delShm(iJob) def delShm(self,iJob): - shared_dict.delDict("DATA_%i"%(iJob)) - shared_dict.delDict("DicoJones_%i"%(iJob)) + shared_dict.SharedDict.delete_item("DATA_%i"%(iJob)) + shared_dict.SharedDict.delete_item("DicoJones_%i"%(iJob)) + #shared_dict.delDict("DATA_%i"%(iJob)) + #shared_dict.delDict("DicoJones_%i"%(iJob)) def killWorkers(self): @@ -1136,7 +998,8 @@ def killWorkers(self): self.APP.terminate() self.APP.shutdown() del(self.DicoGrids) - shared_dict.delDict("Grids") + #shared_dict.delDict("Grids") + shared_dict.SharedDict.delete("Grids") #Multiprocessing.cleanupShm() @@ -1268,6 +1131,7 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): dcorr[:]=d[:] W=Wc.copy() + #W2=Wc.copy() dcorr*=W wdcorr=np.ones(dcorr.shape,np.float64) #kk=kk*np.ones((1,1,npol)) @@ -1317,6 +1181,7 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): #wdcorr[:,indCh,:] *= (np.abs(J0) * np.abs(J1))**2 #print(iDir,iFJones,np.count_nonzero(J0==0),np.count_nonzero(J1==0)) #dcorr[:,indCh,:] = 1./J0 * dcorr[:,indCh,:] * 1./J1.conj() + #W[:,indCh,:]*=(np.abs(J0) * np.abs(J1)) W[:,indCh,:]*=(np.abs(J0) * np.abs(J1))**2 T1.timeit("[%i] apply "%iFJones) diff --git a/ClassSaveResults.py b/ClassSaveResults.py index 2a6a788..13c5211 100644 --- a/ClassSaveResults.py +++ b/ClassSaveResults.py @@ -83,11 +83,17 @@ def WriteFits(self): self.CatFlux.IDFacet[:]=self.DynSpecMS.PosArray.iFacet[:] self.CatFlux.IDTessel[:]=self.DynSpecMS.PosArray.iTessel[:] + pBAR = ProgressBar(Title="Saving fits") + NJobs=self.DynSpecMS.NDir #Selected + iDone=0 + pBAR.render(0, NJobs) for iDir in range(self.DynSpecMS.NDir): self.WriteFitsThisDir(iDir) self.WriteFitsThisDir(iDir,Weight="Weight") self.WriteFitsThisDir(iDir,Weight="W2") + iDone+=1 + pBAR.render(iDone, NJobs) def SaveCatalog(self): FileName = "%s/%s.npy"%(self.DIRNAME,"Catalog") @@ -154,7 +160,7 @@ def WriteFitsThisDir(self,iDir,Weight=False): elif Weight=="W2": fitsname = "%s/%s/%s_%s_%s.W2.fits"%(self.DIRNAME,self.GiveSubDir(self.DynSpecMS.PosArray.Type[iDir],Weight=Weight),self.DynSpecMS.OutName, strRA, strDEC) - print("#%i %s %s"%(iDir,self.DynSpecMS.PosArray.Type[iDir].decode("ascii"),fitsname),file=log) + #print("#%i %s %s"%(iDir,self.DynSpecMS.PosArray.Type[iDir].decode("ascii"),fitsname),file=log) # Create the fits file prihdr = fits.Header() prihdr.set('CTYPE1', 'Time', 'Time') diff --git a/ms2dynspec.py b/ms2dynspec.py index e63f84c..d6a89f9 100755 --- a/ms2dynspec.py +++ b/ms2dynspec.py @@ -12,7 +12,9 @@ log=logger.getLogger("ms2dynspec") from DDFacet.Other import ModColor __version__ = version() +import numpy as np SaveFile = "last_dynspec.obj" +from ClassDynSpecMS import ClassGiveCatalog """ ========================================================================= @@ -52,8 +54,8 @@ import glob, os import pylab from DDFacet.Other import MyPickle -#from ClassDynSpecMS import logo -#logo.PrintLogo(__version__) +import logo +logo.PrintLogo(__version__) from ClassDynSpecMS import ClassDynSpecMS from ClassSaveResults import ClassSaveResults from DDFacet.Data.ClassMS import expandMSList @@ -86,6 +88,8 @@ def angSep(ra1, dec1, ra2, dec2): def main(args=None, messages=[]): if args is None: args = MyPickle.Load(SaveFile) + if args.UseRandomSeed!=0: + np.random.seed(args.UseRandomSeed) MSList=None if args.ms: @@ -103,50 +107,86 @@ def main(args=None, messages=[]): else: DT[T].append(MSName) t.close() + if len(DT)>1: log.print(ModColor.Str("FOUND %i time periods"%len(DT))) else: DT={0:MSList} - for k in DT.keys(): - MSList=DT[k] - D = ClassDynSpecMS(ListMSName=MSList, - ColName=args.data, ModelName=args.model, - SolsName=args.sols, - TChunkHours=args.TChunkHours, - ColWeights=args.WeightCol, - UVRange=args.uv, - FileCoords=args.srclist, - Radius=args.rad, - NOff=args.noff, - DicoFacet=args.DicoFacet, - ImageI=args.imageI, - ImageV=args.imageV, - SolsDir=args.SolsDir,NCPU=args.NCPU, - BaseDirSpecs=args.BaseDirSpecs, - BeamModel=args.BeamModel, - BeamNBand=args.BeamNBand, - SourceCatOff_FluxMean=args.SourceCatOff_FluxMean, - SourceCatOff_dFluxMean=args.SourceCatOff_dFluxMean, - SourceCatOff=args.SourceCatOff, - options=args) - - if D.NDirSelected==0: - return - - if D.Mode=="Spec": D.StackAll() - - SaveMachine=ClassSaveResults.ClassSaveResults(D,DIRNAME=args.OutDirName) - if D.Mode=="Spec": - SaveMachine.WriteFits() - if args.SavePDF: - SaveMachine.PlotSpec() - SaveMachine.SaveCatalog() - if args.DoTar: SaveMachine.tarDirectory() - else: - SaveMachine.SaveCatalog() - SaveMachine.PlotSpec(Prefix="_replot") - D.killWorkers() + L_radec=[] + for MSName in MSList: + tField = table("%s::FIELD"%MSName, ack=False) + ra0, dec0 = tField.getcol("PHASE_DIR").ravel() + if ra0<0.: ra0+=2.*np.pi + L_radec.append((ra0,dec0)) + tField.close() + L_radec=list(set(L_radec)) + if len(L_radec)>1: stop + ra0,dec0=L_radec[0] + + NChunk=1 + if args.NMaxTargets!=0: + CGC=ClassGiveCatalog.ClassGiveCatalog(args, + ra0,dec0, + args.rad, + args.srclist) + PosArray=CGC.giveCat() + NChunk=PosArray.size//args.NMaxTargets+1 + log.print(ModColor.Str("Will process the catalog in %i chunks"%NChunk,col="green")) + + + SubSet=None + for iChunk in range(NChunk): + log.print("====================================================================") + if NChunk>1: + SubSet=(iChunk,NChunk) + for ik,k in enumerate(sorted(list(DT.keys()))): + MSList=DT[k] + DIRNAME="%s"%(args.OutDirName) + if len(DT)>1: + DIRNAME+="_T%i"%ik + if NChunk>1: + DIRNAME+="_RandChunk%i"%iChunk + + + D = ClassDynSpecMS(ListMSName=MSList, + ColName=args.data, ModelName=args.model, + SolsName=args.sols, + TChunkHours=args.TChunkHours, + ColWeights=args.WeightCol, + UVRange=args.uv, + FileCoords=args.srclist, + Radius=args.rad, + NOff=args.noff, + DicoFacet=args.DicoFacet, + ImageI=args.imageI, + ImageV=args.imageV, + SolsDir=args.SolsDir,NCPU=args.NCPU, + BaseDirSpecs=args.BaseDirSpecs, + BeamModel=args.BeamModel, + BeamNBand=args.BeamNBand, + SourceCatOff_FluxMean=args.SourceCatOff_FluxMean, + SourceCatOff_dFluxMean=args.SourceCatOff_dFluxMean, + SourceCatOff=args.SourceCatOff, + options=args, + SubSet=SubSet) + + if D.NDirSelected==0: + return + + if D.Mode=="Spec": D.StackAll() + + SaveMachine=ClassSaveResults.ClassSaveResults(D,DIRNAME=DIRNAME) + if D.Mode=="Spec": + SaveMachine.WriteFits() + if args.SavePDF: + SaveMachine.PlotSpec() + SaveMachine.SaveCatalog() + if args.DoTar: SaveMachine.tarDirectory() + else: + SaveMachine.SaveCatalog() + SaveMachine.PlotSpec(Prefix="_replot") + D.killWorkers() Multiprocessing.cleanupShm() @@ -178,6 +218,7 @@ def main(args=None, messages=[]): parser.add_argument("--UseLoTSSDB", type=int, default=0, help="Use LoTSS DB for target list", required=False) parser.add_argument("--UseGaiaDB", type=str, default=None, help="Use Gaia DB for target list", required=False) parser.add_argument("--DoTar", type=int, default=1, help="Tar final products", required=False) + parser.add_argument("--UseRandomSeed", type=int, default=0, help="Use random seed", required=False) parser.add_argument("--NCPU", type=int, default=0, help="NCPU", required=False) parser.add_argument("--BeamModel", type=str, default=None, help="Beam Model to be used", required=False) @@ -187,6 +228,7 @@ def main(args=None, messages=[]): parser.add_argument("--SourceCatOff", type=str, default="", help="Read the code", required=False) parser.add_argument("--SourceCatOff_FluxMean", type=float, default=0, help="Read the code", required=False) parser.add_argument("--SourceCatOff_dFluxMean", type=float, default=0, help="Read the code", required=False) + parser.add_argument("--NMaxTargets", type=int, default=0, help="Read the code", required=False) args = parser.parse_args() From 8dd78e323b2f494c098c6454fd81653e1089e43a Mon Sep 17 00:00:00 2001 From: Buntu Ngcebetsha Date: Mon, 22 Apr 2024 17:54:32 +0200 Subject: [PATCH 3/7] class catalogue --- ClassGiveCatalog.py | 171 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 171 insertions(+) create mode 100644 ClassGiveCatalog.py diff --git a/ClassGiveCatalog.py b/ClassGiveCatalog.py new file mode 100644 index 0000000..188b8ca --- /dev/null +++ b/ClassGiveCatalog.py @@ -0,0 +1,171 @@ +import numpy as np +from DDFacet.Other import logger +log=logger.getLogger("DynSpecMS") +from astropy.io import fits +from astropy.wcs import WCS +import random + +def AngDist(ra0,ra1,dec0,dec1): + AC=np.arccos + C=np.cos + S=np.sin + D=S(dec0)*S(dec1)+C(dec0)*C(dec1)*C(ra0-ra1) + if type(D).__name__=="ndarray": + D[D>1.]=1. + D[D<-1.]=-1. + else: + if D>1.: D=1. + if D<-1.: D=-1. + return AC(D) + +class ClassGiveCatalog(): + def __init__(self,options,ra0,dec0,Radius,FileCoords): + self.options=options + self.ra0=ra0 + self.dec0=dec0 + self.Radius=Radius + self.FileCoords=FileCoords + + + def giveCat(self,SubSet=None): + FileCoords=self.FileCoords + dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64),('Type','S200')] + + self.DoProperMotionCorr=False + if self.options.UseLoTSSDB: + print("Using the surveys database", file=log) + from surveys_db import SurveysDB + with SurveysDB() as sdb: + sdb.cur.execute('select * from transients') + result=sdb.cur.fetchall() + # convert to a list, then to ndarray, then to recarray + l=[] + for r in result: + l.append((r['id'],r['ra'],r['decl'],r['type'])) + + if FileCoords is not None and FileCoords!="": + print('Adding data from file '+FileCoords, file=log) + additional=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",")[()] + if len(additional.shape)==0: additional=additional.reshape((1,)) + if not additional.shape: + # deal with a one-line input file + additional=np.array([additional],dtype=dtype) + for r in additional: + l.append(tuple(r)) + + self.PosArray=np.asarray(l,dtype=dtype) + elif self.options.UseGaiaDB is not None: + from astroquery.gaia import Gaia + rac_deg,decc_deg=self.ra0*180/np.pi, self.dec0*180/np.pi + Radius_deg=self.Radius + Dmax,NMax=self.options.UseGaiaDB.split(",") + Dmax,NMax=float(Dmax),int(NMax) + Parallax_min=1./(Dmax*1e-3) + query=f"""SELECT TOP 10000 gaia_source.designation,gaia_source.source_id,gaia_source.ref_epoch,gaia_source.ra,gaia_source.dec,gaia_source.parallax,gaia_source.pmra,gaia_source.pmdec + FROM gaiadr3.gaia_source + WHERE + CONTAINS( + POINT('ICRS',gaiadr3.gaia_source.ra,gaiadr3.gaia_source.dec), + CIRCLE('ICRS',{rac_deg},{decc_deg},{Radius_deg}) + )=1 AND (gaiadr3.gaia_source.parallax>={Parallax_min})""" + log.print(f"Sending du Gaia server") + print(f"{query}") + job = Gaia.launch_job(query) + result = job.get_results() + + log.print(f"Query has returned {len(result)}") + + + + l=[] + dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64), + ("pmra",np.float64),("pmdec",np.float64),("ref_epoch",np.float64),("parallax",np.float64), + ('Type','S200')] + for r in result: + l.append((r['DESIGNATION'],r['ra'],r['dec'],r['pmra'],r['pmdec'],r['ref_epoch'],r['parallax'],b"Gaia DR3")) + self.PosArray=np.asarray(l,dtype=dtype) + + # import pylab + # pylab.clf() + # pylab.subplot(1,2,1) + # pylab.scatter(self.PosArray["ra"],self.PosArray["dec"]) + # pylab.scatter([rac_deg],[decc_deg]) + # pylab.subplot(1,2,2) + # D=1./(self.PosArray["parallax"]*1e-3) + # pylab.hist(D,bins=100) + # pylab.draw() + # pylab.show() + + self.DoProperMotionCorr=True + if self.PosArray.size>NMax: + ind=np.int64(np.random.rand(NMax)*self.PosArray.size) + self.PosArray=self.PosArray[ind] + + + elif self.options.FitsCatalog: + print("Using the fits catalog: %s"%self.options.FitsCatalog, file=log) + F=fits.open(self.options.FitsCatalog) + d=F[1].data + l=[] + dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64), + ("pmra",np.float64),("pmdec",np.float64),("ref_epoch",np.float64),("parallax",np.float64), + ('Type','S200')] + for r in d: + l.append((r['DESIGNATION'],r['ra'],r['dec'],r['pmra'],r['pmdec'],r['ref_epoch'],r['parallax'],r['Type'])) + + if FileCoords is not None and FileCoords!="": + print('Adding data from file '+FileCoords, file=log) + dtype0=[('Name','S200'),("ra",np.float64),("dec",np.float64),('Type','S200')] + additional=np.genfromtxt(FileCoords,dtype=dtype0,delimiter=",")[()] + if len(additional.shape)==0: additional=additional.reshape((1,)) + if not additional.shape: + # deal with a one-line input file + additional=np.array([additional],dtype=dtype0) + additional1=np.zeros((additional.shape[0],),dtype=dtype) + additional1["Name"][:]=additional["Name"][:] + additional1["ra"][:]=additional["ra"][:] + additional1["dec"][:]=additional["dec"][:] + additional1["Type"][:]=additional["Type"][:] + + for r in additional1: + l.append(tuple(r)) + + + self.PosArray=np.asarray(l,dtype=dtype) + self.DoProperMotionCorr=True + else: + + #FileCoords="Transient_LOTTS.csv" + if FileCoords is None: + if not os.path.isfile(FileCoords): + ssExec="wget -q --user=anonymous ftp://ftp.strw.leidenuniv.nl/pub/tasse/%s -O %s"%(FileCoords,FileCoords) + print("Downloading %s"%FileCoords, file=log) + print(" Executing: %s"%ssExec, file=log) + os.system(ssExec) + log.print("Reading cvs file: %s"%FileCoords) + #self.PosArray=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",")[()] + self.PosArray=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",") + if len(self.PosArray.shape)==0: self.PosArray=self.PosArray.reshape((1,)) + + self.PosArray=self.PosArray.view(np.recarray) + self.PosArray.ra*=np.pi/180. + self.PosArray.dec*=np.pi/180. + + Radius=self.Radius + self.NOrig=self.PosArray.Name.shape[0] + Dist=AngDist(self.ra0,self.PosArray.ra,self.dec0,self.PosArray.dec) + ind=np.where(Dist<(Radius*np.pi/180))[0] + self.PosArray=self.PosArray[ind] + + print("Created an array with %i records" % self.PosArray.size, file=log) + + if SubSet is not None: + random.seed(42) + i,n=SubSet + ind=random.sample(range(0, self.PosArray.size), self.PosArray.size) + ii=np.int64(np.linspace(0,self.PosArray.size+1,n+1)) + L=np.sort(ind[ii[i]:ii[i+1]]) + self.PosArray=self.PosArray[L] + print("Selected %i objects for subset %i/%i" % (self.PosArray.size,i+1,n), file=log) + + return self.PosArray From c4f41b7dcf1cf39ee23e026322a75bee041c29da Mon Sep 17 00:00:00 2001 From: Buntu Ngcebetsha Date: Wed, 24 Apr 2024 15:24:46 +0200 Subject: [PATCH 4/7] hacked it to work --- ClassDynSpecMS.py | 17 +++++++++-------- ms2dynspec.py | 2 +- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/ClassDynSpecMS.py b/ClassDynSpecMS.py index 0fcdbe2..1ffd72c 100644 --- a/ClassDynSpecMS.py +++ b/ClassDynSpecMS.py @@ -969,9 +969,9 @@ def processJob(self,iJob): # print(rep) rep=self.APP.awaitJobResults("LoadMS_%i"%(iJob)) - if rep=="NotRead": - self.delShm(iJob) - return + #if rep=="NotRead": + # self.delShm(iJob) + # return # NTimes=self.NTimesGrid iMS,iChunk=self.LJob[iJob] @@ -984,14 +984,15 @@ def processJob(self,iJob): self.APP.awaitJobResults("Stack_SingleTime:%i_*"%iJob, progress="Append MS %i"%iMS) - self.delShm(iJob) + #self.delShm(iJob) def delShm(self,iJob): - shared_dict.SharedDict.delete_item("DATA_%i"%(iJob)) - shared_dict.SharedDict.delete_item("DicoJones_%i"%(iJob)) - #shared_dict.delDict("DATA_%i"%(iJob)) - #shared_dict.delDict("DicoJones_%i"%(iJob)) + shared_dict.delDict("DATA_%i"%(iJob)) + shared_dict.delDict("DicoJones_%i"%(iJob)) + #def delShm(self, iJob): + # shared_dict.SharedDict.delete("DATA_%i"%(iJob)) + # shared_dict.SharedDict.delete("DicoJones_%i"%(iJob)) def killWorkers(self): print("Killing workers", file=log) diff --git a/ms2dynspec.py b/ms2dynspec.py index d6a89f9..030bbdc 100755 --- a/ms2dynspec.py +++ b/ms2dynspec.py @@ -176,7 +176,7 @@ def main(args=None, messages=[]): if D.Mode=="Spec": D.StackAll() - SaveMachine=ClassSaveResults.ClassSaveResults(D,DIRNAME=DIRNAME) + SaveMachine=ClassSaveResults(D,DIRNAME=DIRNAME) if D.Mode=="Spec": SaveMachine.WriteFits() if args.SavePDF: From 59a8924c636c0265217353d415005c289afd4a6f Mon Sep 17 00:00:00 2001 From: Buntu Ngcebetsha Date: Tue, 14 May 2024 18:33:50 +0200 Subject: [PATCH 5/7] update for notebook --- .gitignore | 4 + dynamic_spectra_analysis.ipynb | 177 +++++++++++++++++++++++++++++++++ 2 files changed, 181 insertions(+) create mode 100644 dynamic_spectra_analysis.ipynb diff --git a/.gitignore b/.gitignore index a2c2f65..5558c1a 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,7 @@ venv *.txt *.csv *.obj +__pycache__ +.ipynb_checkpoints/ +.radiopadre +.radiopadre-session diff --git a/dynamic_spectra_analysis.ipynb b/dynamic_spectra_analysis.ipynb new file mode 100644 index 0000000..06a5b7c --- /dev/null +++ b/dynamic_spectra_analysis.ipynb @@ -0,0 +1,177 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.io import fits\n", + "from astropy.wcs import WCS\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from astropy import units as u\n", + "\n", + "from matplotlib.backends.backend_pdf import PdfPages\n", + "from matplotlib import patches\n", + "from matplotlib.patches import Ellipse\n", + "from decimal import Decimal" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.rcParams['font.size'] = 7" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "target = '~/Downloads/MSName_RandChunk0_DynSpecs_1556208112/TARGET/1556208112_14:29:42.950_-62:40:46.200.fits'\n", + "target_w = '~/Downloads/MSName_RandChunk0_DynSpecs_1556208112/TARGET_W/1556208112_14:29:42.950_-62:40:46.200.W.fits'\n", + "\n", + "target_hdu = fits.open(target)[0]\n", + "target_w_hdu = fits.open(target_w)[0]\n", + "\n", + "data = target_hdu.data\n", + "weights = target_w_hdu.data\n", + "# print(target_hdu.shape)\n", + "# print(target_w_hdu.shape)\n", + "# print(target_hdu.header)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "norm_w = np.sqrt(weights)\n", + "for i in range(4):\n", + " data[i, :, :] = data[i, :, :] * norm_w * np.cbrt(data[i, :, :])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Stokes I, Q, U, V Plotted in order below\n", + "```python\n", + "I = data[0, :, :]\n", + "Q = data[1, :, :]\n", + "U = data[2, :, :]\n", + "V = data[3, :, :]\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def plot_stokes(d, s, f=15, t0=6325, t1=6345):\n", + " \"\"\"Plot dynamic spectra for given data array and Stokes param\n", + " \n", + " Parameters\n", + " ----------\n", + " d: array\n", + " The data array extracted from fits file\n", + " s: int\n", + " Stokes parameter of interest mapped (I, Q, U, V)->(0, 1, 2, 3)\n", + " f: int\n", + " channel slice\n", + " t0: int\n", + " TBD - will be the `start_time` to zoom-in on\n", + " t1: int\n", + " TBD - will be the `end_time to` zoom-in on\n", + " \"\"\"\n", + " # Slice in time at channel 15 arbitrarily\n", + " print(\"--- Stats ---\")\n", + " print(\"Min: \", np.min(d[s, :, t0:t1]))\n", + " print(\"Max: \", np.max(d[s, :, t0:t1]))\n", + " print(\"Noise: \", np.std(d[s, :, t0:t1]))\n", + " print(\"Mean: \", np.mean(d[s, :, t0:t1]))\n", + " print(\"--- Time for max flux at this frequency ---\")\n", + " print(np.argmax(d[s, f, :]))\n", + " print(\"-------------------------------------------\")\n", + " data_sl = d[s, :, t0:t1]\n", + "\n", + " # Plot\n", + " plt.imshow(data_sl, origin='lower', cmap='viridis', aspect='auto', vmin=np.min(data_sl), vmax=np.max(data_sl))\n", + " cbar = plt.colorbar(pad=.07)\n", + " cbar.set_label(r'$Jy/B$', size=7)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_stokes(data, 0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_stokes(data, 1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_stokes(data, 2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_stokes(data, 3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.10" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From eef86187ac14dffb0e562590e324e77918c61359 Mon Sep 17 00:00:00 2001 From: Buntu Ngcebetsha Date: Tue, 14 May 2024 20:28:37 +0200 Subject: [PATCH 6/7] Dynamic spec --- dynamic_spectra_analysis.ipynb | 67 +++++++++++++++++++++++----------- 1 file changed, 45 insertions(+), 22 deletions(-) diff --git a/dynamic_spectra_analysis.ipynb b/dynamic_spectra_analysis.ipynb index 06a5b7c..b408499 100644 --- a/dynamic_spectra_analysis.ipynb +++ b/dynamic_spectra_analysis.ipynb @@ -53,8 +53,12 @@ "outputs": [], "source": [ "norm_w = np.sqrt(weights)\n", - "for i in range(4):\n", - " data[i, :, :] = data[i, :, :] * norm_w * np.cbrt(data[i, :, :])" + "for i in range(data.shape[0]):\n", + " data[i, :, :] = data[i, :, :] * norm_w * np.cbrt(data[i, :, :])\n", + "\n", + "print('--- Data has shape ---')\n", + "print('(stokes, freq, time)')\n", + "print(data.shape)" ] }, { @@ -67,6 +71,13 @@ "Q = data[1, :, :]\n", "U = data[2, :, :]\n", "V = data[3, :, :]\n", + "```\n", + "\n", + "For now, we slice the time (t0, t1):\n", + "```sh\n", + "6325, 6345\n", + "12953, 12993\n", + "19803, 19843\n", "```" ] }, @@ -76,34 +87,35 @@ "metadata": {}, "outputs": [], "source": [ - "def plot_stokes(d, s, f=15, t0=6325, t1=6345):\n", + "def plot_stokes(d, symbol, stokes, t0, t1):\n", " \"\"\"Plot dynamic spectra for given data array and Stokes param\n", " \n", " Parameters\n", " ----------\n", " d: array\n", " The data array extracted from fits file\n", - " s: int\n", - " Stokes parameter of interest mapped (I, Q, U, V)->(0, 1, 2, 3)\n", - " f: int\n", - " channel slice\n", + " symbol: str\n", + " Stokes parameter of interest mapped ('I', 'Q', 'U', 'V')\n", + " stokes: int\n", + " Stokes parameter of interest mapped ('I', 'Q', 'U', 'V')->(0, 1, 2, 3)\n", " t0: int\n", - " TBD - will be the `start_time` to zoom-in on\n", + " `start_time` to zoom-in on\n", " t1: int\n", - " TBD - will be the `end_time to` zoom-in on\n", + " `end_time to` zoom-in on\n", " \"\"\"\n", - " # Slice in time at channel 15 arbitrarily\n", - " print(\"--- Stats ---\")\n", - " print(\"Min: \", np.min(d[s, :, t0:t1]))\n", - " print(\"Max: \", np.max(d[s, :, t0:t1]))\n", - " print(\"Noise: \", np.std(d[s, :, t0:t1]))\n", - " print(\"Mean: \", np.mean(d[s, :, t0:t1]))\n", - " print(\"--- Time for max flux at this frequency ---\")\n", - " print(np.argmax(d[s, f, :]))\n", " print(\"-------------------------------------------\")\n", - " data_sl = d[s, :, t0:t1]\n", + " # d.shape -> (stokes, freq, time)\n", + " data_sl = d[stokes, :, t0:t1]\n", + " # Slice in time at channel 15 arbitrarily\n", + " print(\"--- Stats for Stokes {} ---\".format(symbol))\n", + " print(\"Min: \", np.min(data_sl))\n", + " print(\"Max: \", np.max(data_sl))\n", + " print(\"Noise: \", np.std(data_sl))\n", + " print(\"Mean: \", np.mean(data_sl))\n", + " # print(np.argmax(d[s, f, :]))\n", "\n", " # Plot\n", + " plt.title(symbol)\n", " plt.imshow(data_sl, origin='lower', cmap='viridis', aspect='auto', vmin=np.min(data_sl), vmax=np.max(data_sl))\n", " cbar = plt.colorbar(pad=.07)\n", " cbar.set_label(r'$Jy/B$', size=7)" @@ -115,7 +127,18 @@ "metadata": {}, "outputs": [], "source": [ - "plot_stokes(data, 0)" + "# probably best to split the time into bins of 40 or of 80\n", + "for i in range(0, data.shape[2], 80):\n", + " t0, t1 = i, i+80" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_stokes(data, 'Stokes I', 0, t0, t1)" ] }, { @@ -124,7 +147,7 @@ "metadata": {}, "outputs": [], "source": [ - "plot_stokes(data, 1)" + "plot_stokes(data, 'Stokes Q', 1, t0, t1)" ] }, { @@ -133,7 +156,7 @@ "metadata": {}, "outputs": [], "source": [ - "plot_stokes(data, 2)" + "plot_stokes(data, 'Stokes U', 2, t0, t1)" ] }, { @@ -142,7 +165,7 @@ "metadata": {}, "outputs": [], "source": [ - "plot_stokes(data, 3)" + "plot_stokes(data, 'Stokes V', 3, t0, t1)" ] }, { From beb5700f5584c360da8efe8143996c5dfc0cc944 Mon Sep 17 00:00:00 2001 From: Buntu Ngcebetsha Date: Thu, 30 May 2024 15:33:14 +0200 Subject: [PATCH 7/7] Save the working notebook for later analysis --- dynamic_spectra_analysis.ipynb | 168 ++++++++++++++++++++++++++++----- lst-elev.png | Bin 0 -> 126079 bytes 2 files changed, 142 insertions(+), 26 deletions(-) create mode 100644 lst-elev.png diff --git a/dynamic_spectra_analysis.ipynb b/dynamic_spectra_analysis.ipynb index b408499..b4e019f 100644 --- a/dynamic_spectra_analysis.ipynb +++ b/dynamic_spectra_analysis.ipynb @@ -35,15 +35,32 @@ "source": [ "target = '~/Downloads/MSName_RandChunk0_DynSpecs_1556208112/TARGET/1556208112_14:29:42.950_-62:40:46.200.fits'\n", "target_w = '~/Downloads/MSName_RandChunk0_DynSpecs_1556208112/TARGET_W/1556208112_14:29:42.950_-62:40:46.200.W.fits'\n", + "target_w2 = '~/Downloads/MSName_RandChunk0_DynSpecs_1556208112/TARGET_W/1556208112_14:29:42.950_-62:40:46.200.W2.fits'\n", "\n", "target_hdu = fits.open(target)[0]\n", "target_w_hdu = fits.open(target_w)[0]\n", + "target_w2_hdu = fits.open(target_w2)[0]\n", "\n", "data = target_hdu.data\n", "weights = target_w_hdu.data\n", - "# print(target_hdu.shape)\n", - "# print(target_w_hdu.shape)\n", - "# print(target_hdu.header)" + "weights2 = target_w2_hdu.data\n", + "\n", + "header = target_hdu.header" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Loop over the items and print each one\n", + "for keyword in header:\n", + "# print(f'{keyword} = {header[keyword]}')\n", + " frq_min = header['FRQ-MIN'] * 1e-6\n", + " frq_max = header['FRQ-MAX'] * 1e-6\n", + " tmin = 0\n", + " tmax = data.shape[2] * header['CDELT1']" ] }, { @@ -52,15 +69,37 @@ "metadata": {}, "outputs": [], "source": [ - "norm_w = np.sqrt(weights)\n", + "weights_data = np.where(weights==0,1, weights)\n", + "weights_data2 = np.where(weights2==0,1, weights2)\n", + "\n", + "sigma_inv = np.sqrt(np.power(weights_data, 2) / weights_data2)\n", + "\n", "for i in range(data.shape[0]):\n", - " data[i, :, :] = data[i, :, :] * norm_w * np.cbrt(data[i, :, :])\n", + " #data[i, :, :] = data[i, :, :] * norm_w * np.cbrt(data[i, :, :])\n", + " data[i, :, :] = data[i, :, :] * sigma_inv\n", "\n", "print('--- Data has shape ---')\n", "print('(stokes, freq, time)')\n", "print(data.shape)" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "stokes_idx = 0\n", + "d = data[stokes_idx, :, :]\n", + "extent=[tmin, tmax, frq_min, frq_max]\n", + "plt.imshow(d, origin='lower', cmap='viridis', aspect='auto', vmin=np.min(d), vmax=np.max(d), extent=extent)\n", + "plt.title('Dynamic Spectrum of total intensity (Stokes I)')\n", + "plt.xlabel('time (s)')\n", + "plt.ylabel(r'$\\nu$ (MHz)')\n", + "cbar_t = plt.colorbar(pad=.07)\n", + "cbar_t.set_label(r'$mJy/B$', size=7)" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -87,38 +126,79 @@ "metadata": {}, "outputs": [], "source": [ - "def plot_stokes(d, symbol, stokes, t0, t1):\n", + "def get_spectral_features(s, t0, t1):\n", + " data_sl = data[s, :, t0:t1]\n", + " t = []\n", + " if not np.max(data_sl) == 0:\n", + " t.append(t0)\n", + " print(\"Flaring from {} to {} seconds\".format(t0, t1))\n", + " return" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def plot_stokes(stokes_symbol, stokes_idx, t0, t1):\n", " \"\"\"Plot dynamic spectra for given data array and Stokes param\n", " \n", " Parameters\n", " ----------\n", - " d: array\n", - " The data array extracted from fits file\n", - " symbol: str\n", + " stokes_symbol: str\n", " Stokes parameter of interest mapped ('I', 'Q', 'U', 'V')\n", - " stokes: int\n", + " stokes_idx: int\n", " Stokes parameter of interest mapped ('I', 'Q', 'U', 'V')->(0, 1, 2, 3)\n", " t0: int\n", " `start_time` to zoom-in on\n", " t1: int\n", " `end_time to` zoom-in on\n", " \"\"\"\n", - " print(\"-------------------------------------------\")\n", - " # d.shape -> (stokes, freq, time)\n", - " data_sl = d[stokes, :, t0:t1]\n", - " # Slice in time at channel 15 arbitrarily\n", - " print(\"--- Stats for Stokes {} ---\".format(symbol))\n", + " data_sl = data[stokes_idx, :, t0:t1]\n", + " print(\"--- Stats for Stokes {} ---\".format(stokes_symbol))\n", + " print(\"----Slicing throug time axis from {} to {} seconds after start of observation\".format(t0, t1))\n", " print(\"Min: \", np.min(data_sl))\n", " print(\"Max: \", np.max(data_sl))\n", " print(\"Noise: \", np.std(data_sl))\n", " print(\"Mean: \", np.mean(data_sl))\n", - " # print(np.argmax(d[s, f, :]))\n", + " #print(np.argmax(d[stokes_idx, :, :]))\n", "\n", " # Plot\n", - " plt.title(symbol)\n", - " plt.imshow(data_sl, origin='lower', cmap='viridis', aspect='auto', vmin=np.min(data_sl), vmax=np.max(data_sl))\n", + " extent = [t0, t1, frq_min, frq_max]\n", + " plt.title(stokes_symbol)\n", + " plt.imshow(data_sl, origin='lower', cmap='viridis', aspect='auto', vmin=np.min(data_sl), vmax=np.max(data_sl), extent=extent)\n", + " plt.xlabel('time (s)')\n", + " plt.ylabel(r'$\\nu$ (MHz)')\n", " cbar = plt.colorbar(pad=.07)\n", - " cbar.set_label(r'$Jy/B$', size=7)" + " cbar.set_label(r'$mJy/B$', size=7)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for i in range(0, data.shape[2], 20):\n", + " t0, t1 = i, i+20\n", + " get_spectral_features(0, t0, t1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Flaring events every 35 minutes that repeat every 20 seconds" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(19820-17440)/60" ] }, { @@ -127,9 +207,7 @@ "metadata": {}, "outputs": [], "source": [ - "# probably best to split the time into bins of 40 or of 80\n", - "for i in range(0, data.shape[2], 80):\n", - " t0, t1 = i, i+80" + "plot_stokes('Stokes I', 0, 0, 40)" ] }, { @@ -138,7 +216,7 @@ "metadata": {}, "outputs": [], "source": [ - "plot_stokes(data, 'Stokes I', 0, t0, t1)" + "# plot_stokes(data, 'Stokes I', 0, 12953, 12993)" ] }, { @@ -147,7 +225,7 @@ "metadata": {}, "outputs": [], "source": [ - "plot_stokes(data, 'Stokes Q', 1, t0, t1)" + "# plot_stokes(data, 'Stokes I', 0, 19803, 19843)" ] }, { @@ -156,7 +234,7 @@ "metadata": {}, "outputs": [], "source": [ - "plot_stokes(data, 'Stokes U', 2, t0, t1)" + "# plot_stokes(data, 'Stokes I', 0, 19803, 19843)" ] }, { @@ -165,7 +243,45 @@ "metadata": {}, "outputs": [], "source": [ - "plot_stokes(data, 'Stokes V', 3, t0, t1)" + "# plot_stokes(data, 'Stokes Q', 1, 19803, 19843)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot_stokes(data, 'Stokes U', 2, 19803, 19843)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot_stokes(data, 'Stokes V', 3, 19803, 19843)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### During Fast Imaging meeting\n", + "Oleg advised I try this\n", + "``` sh\n", + "pip install owlcat\n", + "plot-evelation-tracks.py \n", + "```\n", + "This will help us determine if the features we see were during the duty cycle of the target source" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "![](lst-elev.png)" ] }, { diff --git a/lst-elev.png b/lst-elev.png new file mode 100644 index 0000000000000000000000000000000000000000..4012fd486a9e49144b52154817996aa902079bcb GIT binary patch literal 126079 zcmeEugOlgT<&QCa8ph3Q|gA0+&=$kQfE&l9bjlPy~?+h_p(Vh;&m)38`kbzyK2>{4HVfkDyvp8 zh?3ub{(TJd?{4-xyb7wlE7jqIH+*Z;M0z86; z_nO$-+t`Wn@>>4a3wW$;jd}l4`;?E1Fxbdy*{xc|!$yApaxGTMan&!YR#DEJRJ#^F z*zRE8*&z@MV_VC%=)uYbE7ydf;;vbd6*Y_&EUR%4B zR&-F2CNZ&j?=QbHnD$PuDftw@5oUP&_Tdo^T}>WGqpI*cr|FL7E4>|d7G>LU%m4Aq z#(*Vc&;R&4zKpqb>cIc?9V@^8|Be5-Apd_o{`ben^exY&S(#&DiujMWc;@mm`f=-y z`ZC{LLRM|UlYKR(uY7zI>`oN%mXR&(#T|ij)OShza5O=pq z8`XYs&Y?0~M0Jm-{W)4~Y!I8w#@X4~M1>H3v&NJVX`k&+CkL7i@$m38dpHh$d3rut zI?{`rfttBE^BM+?vKIBll={K7p;q8-+D0??2cmJVIXMTA?`q{JTt<&u~MWs%{HpHI;tr?$|I?#GA#Mr(^F%UlR;@|X=}F~PT*7wHR*cq z&ok*ttK~SNo3F+u^FVEWimsk;ZeKyKQDvwQPx<2xXYTf@WG>aXz7XTc^FiF|On?2g zmzS4UV&e1ZB%QpMdTygs3^Hd zhZFH0v@?vWHSF!}ljfwH$6`8Mhb0;k&YyN3?@`A)SAS_pAs=MCr!s)Aj8^>Y*|Ul1 zPEYgi@9y=t7vylOC$TFlD-Sr+8j>}WE(gZw7C6$fn)Fn=KLl;XThiCpU%h&D;@SJ# z=B8S$9dH%3$$9egude@+bnfXt3kwTfy}j{kw{n}{J9G2HC8~R@PksCLZSC5%FF(2( z+u1KKEzW(v&pXi-pqzB-_Sz3Yhk^?W3j+@*vycV-ZOz8|stB=&2i;Z?QW5?mi%&14 zU#`d&D(n0D`evw*`J+EJ?oJ?|`eL?)((rK79x=x+-On?6lX7!~E!%R>lLw-mW0iWj z)SJt7YH)0H^!$xqSJ&rWHSpanqGDVfsi>}F-t=6(#A8wIxsC`vmAqDlk@~&O9QEl2 zGK(`m{Jihoqo|kAx<5ZXy_4_qt^UT;pw(xjHF^@SlDC`gO(B z)gNUWj_4I$xcSEhmA7{`sK5C8kefYR#16L7eDfrJpO0(I?)CvS=qa9B@uX6ytf|eQDZ$8R=Ixr(N}GS&ULiy z)!`xrJC9q1Ta<{{b{(yYS9Z#;##O>BnqS6lPs_>4*}-@DjQ3(g>c!)w#TsdPbbGh5 z7|+O|FX?aVpX;{Lbj<1#{6(&Qza`HXUxuBg*qUwWiyu?Hy7iF8rxb*6tbyA+J%X{u^B)RxCIKR5P~SJfvAPl!iqBlfT=u$AZ`Ui`a%nF)^b$k95Jy*Jr zK}IvxUAK5?Zh(`U`$Ko`7MJcIO`WP5B#00%R*x^V!oilTAClV2ge$c)a_<S}6%cQ)*b`~3Nw?+yXZ3Q1KlGm34km@*@acmlH7M=RYB z{wr$fm(O-r1jt;vbg3oNG(agz@(4vc!^oFxLc~(?nR}aQ+sY%}ynUN>ZS-=11KE+t zpXyIe{Bl-C2K!of4x5Ujoh31z62KXH+&Z!)+cNC(a$THa#x9A<&Q5*bUB@+^>*Qq` zSF@N4Mk6E>kjh$euSVMr=D@kW^~4HL<1=bA8bySs~pg zIt(;&?caY2OMgDYs4~;2awl?^Kxy&(RP!C_ExKYI#!>D)wFk9EA04?EH{M%)@uQ@q zWRiOFqksPSheXrrHEYa~o4oJdJ!fy9N$LS=qbvzSad9fkVM){zQtzIg_DXT!kPAH5 zp6}2-FpwNA<7f7x%oiJ_af;rG62fNL{F0+xCQnZ8!iC47$1IcRtvMG5a1*s?=?&zu z8~W5NEzaW;c~)$XmJjAh zwCk-36~6jxZhAN>>5JbUF{5dZ;UwgpWIVw}*N94^Is-q6I4XxcHzN9wM#{M}XKwN6 z6((R+%~0RX+H$WZ`K@7RS`N{pB20tt`LApy<11DXTgg{3gKc?sjeQ*CGiOS8)Ca1> z%7s|9=ZEkals!gJwfcStc(f@(MIiJ(8-olxrTVYUztkfn(*jDQ11Q8d?4gI|Yvl7+2YU86rvH^OKW~Y*Oo4#Dz^~Or!fzB@c8B$MozOFJmT>$27b-H z*D@?yv)Si2tX)f95~rPI)^nR#@2F~=0$+h6-PSZPP%PFdatXys!!ycl?4F4Goc*(7 z>FMc42rfr9-}*$=@^_`Bc}xX|P)m*3L#|(6C1CUQM22}2j~Bf$P2bv^nb)ecySLZ4 zDoQHKqI1p0-7j2=T_^f3VQbM_QgrfafB*gW3zm@pWrp*+CEarU4ZXJ=F^jzRond^m zoHZ=ykY*^4Zhp+%?uR2;2*gX$ek>6tzlS6G!cfvK($_NcSXi~?`dgH!#w%qFI1PSD z=R9_-v4v&VF7vVnd(4qGQXI+JaXA`T`R!DHdTO<~T8s>f^E=vZQiXllRkO_0#}OxC z2ojOOL^j&9rD-J{{WOw?Ad***=-V}lQV_CEfu)wIpxp$8MFxog?5Hs&CV&C_0 zbaGONsxUp&t}CjUWm5Czi^Z|JdATD_Qm%}AS(1Ys!tsEQk-I#c3zOpyTiuS>Ik0~ zA76`A*X`T4FAmY}INbE~>C;5-t+EGE*>=*N%LQ?V?AeGIzV6C&TZsg7TSr$#K!Jl;>`}L1p=;sxG}WEc&MjT3tlwss5Z!kBs{giRF?KI>m0y>vtVH+^log zX<_DSL1Iyd%Rs7Os}sGY!R%s*hr3tl?`7M^feKr zKaQ6@CgpFQSHHL|c6WrB)AJB1&tlc`9Rfx^x_X?)9h!?S^x1cR*y?pVv}$NL1we5? zWbL1S9%bgyv7zZqHXBEAXg$2vdsWoF@9)WrFaCbbexI6c*=of6mQiZqpclKIF#^Xq^FItKhr&@W+*S489zr3QMv|*2^V8&bM5vy$cB$L?jLFlnPKIEfw38TDMA_Czv&$nHM#&8tyRBe5aJT3^WG z%jd*wHff)S09!Gwbq=W)()U+nZ#kf(rlS+io(Pl+FnRl}BMS@5P9d|vtOek~J_NwV z)-q*JP6~3l=G0(wfVf9hxJW~>yQ{g#H+wlJa1UGWx{*BYqVY^DGlD1hCE(=cR;38h zL?rsgF?_n(BM-y~QE^#boJLX!zT67<-K^LisH1~L{*3&a+6|t?A!4ZkJXH_OM^OP7 zM&=7NEIE1d#s zToe`5Z$bOM&+PMij#)gOoQV>?`X;v5b-1Hw4vB+_ryd`Wa1{eG6P3u+*^W@V&=So8Eo zV;>fYmc(({6DO{BxQ#2CH#dZr*nInE3740p?i4h60(_pkEg>eYWT_>;J;PX8JxQ(N zQo&$zhWWd@n;ORu;n9yDZx8;rvR!V0F1Bh_qRi7RJnAS`EN0rn_$Q&>G{a z@N}M?v-9%r6+DJ&i{m&`@Lj4z29_-y0RF*zc&yf9udMYdt`iM>e1R-Lzy z?o(TtNxl9*$4U{d_bAI3klb8VhZh7(AF1R7BREe)3#i%mEk+ncJ&+)Urh|dHX1#&i?z)@(*5kz%Jek=BkJ%cGaa5 z;iHaS43^x8b!7{y?5RzE3k+%7EoUwgJN+F+@({2bHRh1!a}zvZo;}q&P)+B$gY_rc z?5Y#NeULC71_uY*_Ex2$NSz@x6`A*n4<$RQug^p=N>Ue0RPEfkhcAsPxuv9}f{t9& zaG#r~8SDNS;T;bS`?A>r|7m!zO1P+tm%0wSjj;w%h>d^)7(`6mOBl`3-emo&mJ@1*5C?q0Ixwvm}9-wY4{N4pyT)E z?*2LoKnc}Ap zo~e?(_Ebpnk=Pzdw+r#gQQ41nA=_Lo!K-o|T*0_Dw*LOdAKE#n%VC$J7esA z_%5WEOF(#66Y@)H>+_!Kihu(S0Oa+n`<|_1=DCn_^PhiIEQY(guN+T0|7QQKTeY^! z0WOPeG~}>)pkDTNmK_JJI|`NP!P|yHgQvZN^*yvS)ASS_zR8;jkB8^Yi5}F+ebL@$ zU)5t__9;~8V}3g}_IfcTl_9qW;g$2N*KNDNc4Y&L_ywpxPP1Raj#<9+E<_c54)Dd% z+R@aciRzg3sBYd_tHh&-TP@*9e58vU!O$D`h{n8lakTW^yFNt3+(25%1U6v;3V;h} zZor`NKw~PmesQ78^iV2p(}x1??c-CAZRXJTc^{IyVrAHIO(Mta-p9`sa6nn3)O+h> zjSQ^ps@4>+O`%&=BR#~pf-ka?p565eMda)-GI)3ayK z(!Q78pQGn=s$riaF5AsKyWVe7Pt{hhmY5FLNY#!7KEYeM2q+T>q$7rZ5-Q#Q6>lOL zx!#2q8$5yfJOB(y2!n!x&vOt#093Ql@1;BW4evrB3jBUozrblY+q)d}%2l&Ub0svC zBmY&l*$;eUVrGcnQxNv%&6^1BG+&`6e;eAez6WT0>ZYt zz&ga8#}Xh-Bq0#JqhGy~7oL;H)?lXqqJDmNZ%gXBOMt}Oc$h4n3s)!UtxZ-;|77Mn z3_Jp?wm8CkM7#b{!7a8mOb5@Mz-lX`4!k7acQ>~O;vk5q>;oq!Xe&-MCtIC3ouG!optpXVx;%esk>EV(?b$4U!e-Y<+ zuy*RXUQrT=i`K`gw&h!vsF=~9b{t|6U+-->AZYmhFG?*TEgeO!1)Cj{VmW#11{R#t z!UT*Ts3(l{gjvTWjg+cG>d4x4<6vU%r-r7J3-5(S!?|=6eGy71#nsvmJUhdIW#W{h zu>OvO&7=HtYd#MLF^(qmBQa4pBIsx@_-7K-3O1GR%H9 zwmnlBR_&rB%|MR{v=<8|DG5sgbSek=jYta|l#ju@not+AF1`IDeyxb-;-seeyG+2x zM32Qe><7(fO(#!bH`PM6jYp>Ju8O$mCo%OHT1hk%EKZRI1G&qWFHgMNqHKc8kYoj_ zx^?^Z1d@v(&#hgz4(bFc7a)7;aSA&i2w#X-il9R8WN{tj+j{uIIRY{uW!aDSC}IO7 z;u8=?4UQAp_ZTl#0hA)DFeT@{t$9OIELMfchxu23-1+C`%}J;j1m`NXGdj&EFd~${ zF8wFVUGb^GZsMSn-e=att`y*SfS59Q{m<%+IY9tq*2S}N*w%EF*L6Uy^@YxMUai0W z`YU+~MU#ye8w2?$v)D_itv7DmAeyHj3Yz)5JJ*nzfUxETFTMG_XkqN*PBF(6Rc2vT zey1Uwf&GswjZ{Iq_hgG*8+qc;8d4*{K3@~Makuc-7m8qiM3W#}2KC9TB{PwTz$_9j z2~b9;G0?G1dp?ECA;OrLnGb-Ekn{q{+MfHt{A=V?y7YGw#1E(RxYX6FSLc@Ihf}+| zzROu|EM*qmBW&~axzF0*k?(l|b=Dv|ty}p-xY=c|b#lx%NvjN5*o}OBZ75^yJX*@^ zfHznm%KsT`W^gS>MCHc`QO}N#WViPxr^(4({B}Q?2=wh!)6>ef{fULnUL>7RR!z+O<6wM!j?kknQ{4++L?$;Fv18IHBr3-R`8M_`28U zmft1ZGr}d6B){mZHZtf9o)7=ZsSx6H#r`3r4VEkkX>t| z%$ceO)I_=CctqR9qL{Y2}| zV#n%+h^UP38(1aDlJ69;jUjjr0$XgZvK$sW`eKe%;%Tpqq2kUjQS6N^*Wq3u_n97x zE~TZVre9uM!Zrv7T#kkAVP2nb{sRQOh9z((vyQ?%NQ7y3@BT^PBSb2B`%h3ME<_Cz zgc#*9S79~r`Wi?k5vV}qsF!)lrFu}Y?ZT?}t>j!q1>H6H;*})?gG%mIHch4j$b%^e zt|V+3M7VVR#EzXilfY!op;XjEp2y>vWC0TqePQ_NvMJ9cHjacN;TqQB(Cjb9t-P_sk} zPncvhXWD!_g+i%FW!kYL5g#Q1$|<1JIly4(!cWoJIaO0rLmDxS&vlcB+6%+Qpu>Y$3Htm_3!A*E<{_WiqeU&Or#R0%#!@h$f6*sS-rewJBKy9vjG z6s?T1p<)qlfo}ul)z!-Q42paY0M(jx+mwcufDWhU#;c_`DAq-jkH{mqyllFanK>Rv zzFpe`K$Chn^y=-5-fGu&&!utJ2q{m!TKb#d0zU5sV9!bD2p2x`7LJ^vuFe4Jgb1OS zK0H%1_5Diyz$A)eBOUiAH9*5(6|gVgY4|zPJ<$rrMsgRkXuuJKMJ11j9@)mLXW{#R zy*<+9C;+M;vr41o%fp__=DsO` zpL4JNAX+Wb*KXk1kc9UXs>?**Q@qQIACcN>QTl;t3L#o~FGLT|q(X}tRP-}{i9a^gI6PzQUeqtYTJT+Pr1avc7R zsum-=2|bQd?MT$Dw_LqMGir}KW*Pc+{{hO_#Dp4>G~(%S0qwxB{kJ!(7a`P8t(dcn z%G5rt1STt!=nA2|$lz%iwL?`0n6F=IdDhQkdC@6C)IkoK8p*#Kq5d*#;rz1ox5H?n zv==xz34cFvO2n>5ElsZoOj?SDwhw4?6HvcB$loIHm1BvG6jGyhwngH@hktpw>lV2* zSh+lj@?6rQ*2gLCB64R1m4j?ZG~?nSl_VlzI9l#m{HJqkTH%{URe|-X|w-A~ufoM>us{&0o&kjua!bKLWi@w1YC3rYY6qy+sx%(v7oX zskeEnAf`9iRZAFm8)E~a%2;zW7bE{P(jk4A0f|tui#!(db$?M`$oKhc+cu;2jkpKY zU?M`f&vplkdMxA`RY#^Ks>Ty}x24#{vgZrDN${29cd7ldzT0VbFV~ztHDj9& z9`09;FIIg?Hvq9uWN$6^Yuiid=;%<|^P}o=_^s#;b_n4ylz{p?JL5&9DR5JbXldVi zXorMqBZsP^ladKZxM<|fx+Lwaes_sl`BMPOs%S)|c4IH+m6~XvQPTesnrQh$1qny8>nf2l~{eX7~`0~k@#16JNt#M@{g`|fX~ zFMv+-sie8)63-H$BrU8HVd<#HFaxBvqU2N37m+&XwmP^dB{=cM7H+MVF)NL@Fi1OE zIlBZ+9(Mg1Ho949jr+pP=$Y0n+k62hDk@|Hwk?Nr7-b})oSF_1EcoiC?@u$16%A>d zgbq99$i-|0!9^MA@#*Ot=_?d(2yPJcE()1^ex*H4i;v=L?kF!WPdvYl?4yGF;{ryX zK7+As&axAvmlI{P%iu~z^PtB!PiUm(XSQk{*3Nuy^&t_pquE#Go9vql(a)ZRi&jOb zHi+zW-w(B0xhr zJ?wFLkb~?QCI)SLKW$3Q!mj5NWh(2c@$-02FP$9>T&cwiWAe!( zyr?BmIEb8q+7gG*Igk9Sk$$-}^Y7QcUa6p#T1|wB#qXf1OwakN6fZ-fr~$t}Na=6R zh$p-(eqJfHWO?Z#(nE&y6d-3J9;p9!51RW@q2`8N#}tqWf&wLx%^sumqkFPpZN0VZPdPrr5h)5a=6s6VMhEXGG`f8#%+Ie;lK*{`AUaYF! zdB5(5c<757M?77R zi5aezCna)k#JlCOi~qK)*sI^ZJq@!-PMktZhOr+SpQ$h%#3Xh?Tfc~{qz;gJ*}xzR zU2Tc+PXdmXLxp`&z|IAUVY{870;}iSn-al=nrT|c*?{&BEtdc+XcQlIM#O2@p8Rd5Ls*#suOzPW~|2@s~$&=D%i=FeD^((8Eb`Iv_36}EoAdG-WEeI2?iqV0+5^b}hCP*sO zx*MesL;@z~gaJHw^os}8&MOvb>N%IX&|R-xGr{Kpq}6&k|E{boNjXXq77~A`3i(?T z<0gw{s>+7}KkvnL7ODdI4C;9zeB=WUep6fjfIVFLDBK^3?hx>JVXaACH3NcDu`N8- z=tK6z;~u&aHm+Q{;+ye=WO{UT^cKtE)||9QJi1Bfe9KcwZ_&4Z8_4n>8+P`gFQ*CO zp>jAEI!>%rMv3fo`R$N^C;n9&WhtULFkL$rH4hk@B&e{lIj34$2q!NLSrreK&JX6bC zvhZEa|J0lNZ>ke_QMgdEa7mU>~G>mEgQ=8Ro7vi8!j}U2ixR7xyGO` z1vaCK_xE;F`{sP-@i@Q@t0SGx%U-Z405nqa#n^t^Bjurk3uk!WCvrN>TY{(8ehwG$ z#2Tm`PzslQ>G$Nx-kd}B{dg?cpLJOy93qf48UfJa@%AYyr2QC~pNi@X02J&o+s#}V zA|L~qBpBUCDCg1GO?~$cM4m%c4??snMu;}FRM(K$(|XeNDbtj!`9(?Bsg%0gO33%w zx+m#iFqX%&uBf1R9rvXn%VFxyTRbjq3n9lHfV^Uem4(GhvoCO@Fz`{6 z=R$ovHHEf1!e7GGG5IssrNa+mjDAn{*)B0sWxBPlpomTDrC|;35C`Sb^wVR$mK^J- zEINu_sDx`a0gSGM(W3SBsODJoxka{3*!G4makYUuK0p zc)N!(*^9%A^*oSM)FGr>W!2>$fN|1lFD{zJ|@H0XCB|g_6|> z2|)F5_TAu-{g$?+s~IB1oon-TAyS=bMZ-1;8W|;L-BT9lE`O$Zi?-Z!2Qj!&)-bZz zn|TMjmf`~S&(5wRVpc8P-$f(k!c(X?#l z3e2nIYnTH+jRX3l++Xo5i$E}FYJBtvgzg0R&j8?rQns#K+RX`@UguXK)0JZJHFi_( zONR@To1nCjY8p?JV7?kkXtMiA#G`?()AmhQ0tl6%k}(T#dc&-(s%P*`c*J)g4cv zhTGnD>Yqz*?q^mhSlf<#J&97I07mQ24vb78`j0&GpIaW>y}GpE905AMdM%TzW7N#N zL+{e)<_FTSAAYT0JtHQQn&&cU<~`qb2nK_v#BC6Ly$e1*P;>{gcT}HaCZG2faBVD^ z?b#JCEksE`a?p4x^Jf3`m~aW#Myq$5Wkk9*e*gYG6Adi$JbN@{KzR1&?FLQ@fvAvWvHo9qbaJ1db6=ZL%aP=z8TmTZ4 z+C0$k^C!Yaz8k5yAA?aX-yM(A-1})4d*S2o+KCfjUz$KcssH%%U-=MY#{LK;tyQs? zS?aQ;En>^OT}gPeYT>QAx3_NJ7I8=y@@Pr9wzRawHDt$o*T;uAW3!-JVUOl&zcjo@ z?1UBm6tv80po64zgEA&TmgAWo2dKN4qgxsE>eUH!6FzhsbZkJ3rt9f5N=!TiPWp6$^u1)N-pgoNdae{j!1%CZuoE znNkytF3-bDy*#i4uRwbxVQnpy&Bewo7RR>o=mE_-&0Dx>bpHk1gUd z)h;}c%btQlaI`KVJRPI6W6Sf{zui?u^IOhA3KS^A3Y6;=cbSsXf@=)>)`RD++nwGS z2rVASBd?t)ib1MGOxuaqEsltLHr7u*GmuAd@!tLWb!h+a(59HUhP1;Ueg=?9X$3@3 z!8PmJdWi-Ok;FU?&GwD02uCIKN~l(TqUmVu8P2yOqpx&HU>7-Bkv&?c2$ZT#B_$F@ zxG%nCUk`^>8adrCFZ0<^k= z2bwb2r9iew6zyq0inlZ12*HQC`)WMEZX=y&Ih!_ZqD<}BwCSn7$871h;!yxBHUl43 zPloC(1Sq#Y3Y;-zavJV9LYqo32Rm23%APdwbr4+j$Kbd_@Ayk333GSbv$mtEfo04? z;7hs$b>#c1JPw3rOG5_Dd=wGUaU_elng*soKU5C?6$$T%LjhP~#kTiV=o9iRAv2l) zQ}UDx8XXzD z4NL+1(brb!oG*5BHaC2aMk(;bp;Grar3kZjr3>QT_3|^#Jl_s_kHT!pP64WxpAuN{ z`{(}iFY~;D4BnV|C}i@ZwQ>CQfEbASG)M+HkDw5ndeyx;bM`D33`2>Mr%#<4fFoi>Jyve37acA-ljH4bFy&<;KA2>#%PA3U z;ZzWe0;*V1r!R2HS6O+J4pfWc`7LSag?<1*&k~3G!jShuXYSS0Kz28>bkJ(-e((Rb zd;XgeP*4L5(#gV~>*`cX7RK(wC>V+jjEs{LmT!t4!cI|plcf2NQBozkH{L+r3pkvq zm2DA%HZCX4@u`k9vr0GumXau_My|@;+gTlfWV4TT(w$-+zH!6GUnP5J zZ3Ou98DoUrqo9D#L9hcM{3n^r@ac!rJe#do*= zi^u(XhCb(`B;5s#t9DazSpAUuI*!EM2|+})*DHc^yhM+vPHI*4Hn@a+)GNG)9lpGEr)^XZxS`PQwv zOTALdHE42g*{X!rTnXt9sozp9yml34f$B?900QcOzwm?~qD>vN?R1}xxY(nI=4uLU zo@o2vIsI@XWSQHcKSTOB@U-~by9&lb%M)n7JcABo3fBZ)6<5yUs(>8m@MQ|!?KwPY zK{OQ3wCZ^-I+%06#?~hy8$OceGe@b#i2$HEqa!&mkLpUt%lQ(6{E$FK0pU~kd5aSe z>6}_;iY;4@6+z8ww0%>P&))3b^C{&glo(PrUEV=VjH_azXYTKJ9_u=f%@p3!3$K$I zHq^7pQRY5&)7lt{V_mgujtrMKD?_gmG50`_<^oJY6Ilt3U{&;xkX#(l39;-bY!2;^ z$#=iu2*AsI;l6X$u`d2;22skhl6fhY>@C-6uJf#^lA z6H^8QO+e#W`y;koJ>23c<@-4*d@obC5iSJr^(tQI|>POI= zgA&3{K@&dR+l`$D!q5c99DKQiX4?1psbh3j5Na4H9_?x_lnRs_T|UiII~Zm9QyhJcxPUV`uX1uitMxOv0~#MRa5Yh=Rrpye_vgeDh5|d- zVT>eVdHnXC-RP@+Mky0=9y9EQ5J#ROaj(JggjdjB-T2Kh`UVC}G>40(BXwCvqi7Nm z5(0}l>#>i&N*a*uyw?}z6f1TuPXKs$Dr#IiiGqbrfM!`)4M z`O4bz9X=L9JuEZI;Al>(G!jTFnhL$=?QNSl9K$YtU_Gm(cyl&G^*AkP!0E>V&6{|B z%Q0MY=B1^otE5AZ!^7R0a~K1Yo3wk>CNhk)V4vg%-*hORSUhbht_-63C-bzTM}y`R zD&CgQEUIk4;DM~(shJYYamCX-*xf4g#+pq5%i_?Bt0;Sfkjl6=`dzRzIGlL!p=*du z3UAoK&tBmHzqf^4;6eY(#PuNr2B>f>?pB&tjAV5O*S7)j*rwKq;dTmTMg(R%p-UxP z5%Iq1)w^Uw9%Du0?Q(&6$&8cZJgV|8SV!1IrknpF&spd7-Jjm`I|T|?YltilqC$dS zS~ky~MkN^!QWu$|BQ_aqx%yanW-=|qHMIlO0!>}wOC-Ym)KI(mB<#J4WFYhJ*Z+_i z6|&usg!9DRl_4h*MTsP&P*|S`6hJIiptp}=s)B1YK$k=@KbXa+n12xDj3_&}oGEdWzb1-Gy7F3zZ+E5$H&s zgovpEmzu!h5*vMPik^(!L+oPOnI`^Z;AdpSh8T?yK}|43L}T{Es3r3R`70X1Pw^y3 zuR#>z4>0vX)GYF1h#Kd?iyAOjXgo%|LG&tRe*p7WvRs&e2Vk1O*Nn_ zb`1}w{VPg-20)7H;aVS1hxtj8@(36pXb0Su*bIrv1Li;sRfzx~Sg!z~_DD8|{(3Uw zgPeT0GD=DxkF^$&(Tz$BKWiV9L4so7mKtAQz&jlVv|$byKpa0yOY_#Gi-A&e4#F(@ zJaJH2h=-)VEiZ!zl6^vjfKqt5{N}nwB0S1^9y1tve4T{^Dt`r}ujjZe8PR{b>bJ;cFWi@(dvy#S_CE;+{yYS(fh~?qd=m3j{6a z%1>PTg3uaM12d{L(Y7_uO+E4zA{6{mDKa|JNW;@-G2Y| ziU0kk%Bg?f$bVjx^RGGnKYw5K|LQl|?o=8DR2tdkFQUK*Iu2eSjbW%nROAOTkd6n@ z56DX9QrIb&V%F%aT0)o~yg9`2bL~{#3S_L;L58vg#@*V-w5?FGutVf3J3 zuRK9!^fz*_Q2^lmBsDS+ZVhTThGG@zv|K!-n$me(*Cjg@F%POUZ#z_^Ip9f0vHspgqxd zu>;)+Tn-g=NWxoVr$BS1aLSNwAr!}ne~qxXi9A{YkJtflxdPaJUgCtyPhAo(oMtyr zc-iU)R{;gq?2-0A|6S`HKL0t%q%8mz*vBJw; zRwPrrH7nx3r&yU|e4U+iskQU%6TyEH$#g3kfMhfo$g=^EWpb$fB1Zf2FYfyp`+ILT z`|xpQ+$1aUU%#$*T{Cq2ufMH)@o@8h9^1d)&z1OJ{*CKv&cKgA#<`*J6VDvPHG-Z9 zH}|5a>Eu0Io-bKe#fTN@a3Z=Af#ZmrkYw-$92rx~NvQtZ05R}JQ3;pK&CO+>X4tTy z0mNe~6B9e|Cytq*p>SwkxNsrq$E)@%_jyP6v5EsO-|rnD?nJnTn=P_o7{OEwfkC8u z4MTVGoe3FR+dnbhjGB} zVpr#*d7+*a5%b@d!Q}^ScZV;tC?tw3M@YD85ZmKqyB%@CVaA-8nL&M()pw)eN{+3- zwGM#$6C4UJ55+O6K3w3No=SkH&io0NT=d3?J(RQv(P$Q4;DM87NaPPvc<1IZu=ex&u{{e;5CK&YoCEWr;8!8!Hz?i?3XkU>T?N5wf|F-^ z`@)%*nLQ-AE0>=L3`QCZP|De%K z;qfyi6SEPF&}rno=3{!o8K!@jEnTJmM!SN3qed|_oDK|?G?2xH*QZXF9~@Bt*j_f1 zwMp_J3kk2~SBrN}9r*goQFJKSD&YBQLNB3Rdl}nv7bOmW-H{HB-<&-Jibe=XnE4_S zhYuJkG%3G4M|x;oxu3T;Kj_7ADm)7YDU*on<8V(3QOR`j7J3ut4qJsdWpUGA#ZWhLljyo=~XUJ6Ktx43{&=*&%EdT~!gSyG8BaYnQ++MCSGv z78eyEcGs1cxXr#C9UsSYfQpkk2_-paNr2XOQwJS$U-t19^E9o2XVKA$U|T-y0&X)l z@dhV)Y=owvF8xg;RsmQL>@fRPhtDGD{`tT(i^cP*Uw%`-R9Kee;lrQi(w`r??Fe1*D9j^#NP$M6!Hsj6 z77NAF>;e(Cn;E%O;I86Yb!+~r&~on7jR}6V(Zv~(y|pCka`l)GN8 z@-(S1)P8>ICp3k&Sh~`)nta&D!F@#B3&D@;W%*Z*k3e78Lx4)WPxMBdN7Eq7k7O>d zCgA~s_5;(2u)}V4-?U-6r=Pj#h<2uI-5PgWMuym_Q$$`$ogM`@6)s|(^}3}J=2^BK z8=94^!iQC=gXaETioHv2$ziyoeP(Vh6emQ*$cW6R*(!?lCMg5c&<$h5=`k)}H1 zxSYYp_@3~4uBUBh%d;5UjzG-er;@?@JgjzBq&aKuF6xI8v8#)}O^Y)AlAA^jRsuyM z9V1K&W}Rs*u&2R4y1R#(3FVZcg2_M>?ZAXZL=z7MBcnI!_}*3DX7Nn9MHEb87{nUsmvNce-_7+Gu?Y5Rye0y&MH2Axl7SY~No z(KTmh3KdU3qQh(9wAo?231aiVvwp`jkf6AH)-lGAdkZ%YVa8@Z#Ra6ME>^18^`y&V zxKY_@2m=RGsVYyK(unFDHY~4Weo5CDZ{qM~RlBNpMSfVXivTkFDMSwey*_*o-*MgD ztuMWI`C&6~_K+PWuc804ChfVWzl=7N92M4;eH6&GHeV!{q(Vj)Ig$SP{)003VjNl9 z=E-aXxIG~%+;F)*qgO}_uN;teq=xygiWN2oCqn}IdWcPU2vxjSKK>3afD?2E9ui-{ z-Wta-#vP{E9dp%@*#!AbNDAWo6m}kf)Ns=P<)EJ!F}q3@dYxCj-1O#9vMA8j&f+OtU&V zH}>D(F|7DalmU^DI10K~`nAi~29esR?+~{V55{d91{ym3RFkoqlu-xu$S9g%D;)gqkHQqzr#Y z{B(mT)uW`W9P$P~Mo|}flvj)Pfi)(OLhdR_J5S~)n^0-_VtBxHUHAqC<@Enm5K&W_eI9uWb zv6a9NP3CfOqxZ3rvqP&*iEV9%0dW)rOV@o{SflB_Fj`K~GlZKbM1KRqHfzmJg?JqY z&Iobw<4&yu9HY)SNA3x>ftbbVrF|?-k~TTWH+s~1xHHY`cRqLjYIQZxljgf{Q*7kZ z;Hp)IDKm^+WWX#}VoB=yC~G(#_7Q+_-CUcKxPv$G#M}pD4 zfLZk3gKo3+iS28qivuH_Xj68q+qT6a77pQXg2p*$d2(}|Okxc4^3RpW7){R(?TZe_ z91zYw$Y(c&cE54iKyGQD1et{(Vr^D9OrBuAk`;;aSV=Fq-c4ZDCr}!2rEkImF?kvU zV*;C#ID^S)Ybee{(n2tsW5kKf$somVrepu7o;rP+gA&ZE--wf1Qh)jFk3KSsM5If! zAo_?F0&fu-)?8#9({X;ReI38IKwab2rPD44&po*g2`(D?CONIck?jlt6C{G4x#Qfr zlm+?s&BL2H#u#$UE@S1R>UB~)S`kCRR$gs zRd8Sx;CUkO5iJB$#yKr6*S;{W&B2>1*u(R;2ym;B*Y!0_d0sv?IBSlaYJs6}Dx7%< z&{)MR@aZIqaA1uKTw_}?bcfPm0>`4ux7*BQ7_Txy?0Ai=7}>#SbE4!1m0)M%oVnmB zj^&3ek~g>Y9}?(l8tuq_DKF5?kHJ@$0K!h|itB%ml|KB_Y1^C)dT?T6lZ&PT~OB>uDR& zC!^W8E)>=?B)^5;YCvXm(DSGVV%>+`ja{zQpw?bH{A!CP#&qORCPJGY=a}>hUCnung!vIm>del4< z7kT_a2!OJ8fOQI(LDoQQPnPYnHbKMEu{95?7pSrndwh&>=sJ+KB=h#|2PvX(l)R;n z4&WGt6bHbxl@rAHOLl6OC-)zf(u%i zYL9bP;^(>wm}4~HT_oOq$~~EUyFajmmywvr=|`CPeh%7xprr$91BFaRhWGPI`eBJu z??*fBUw-%(n8zkGPM@K{GXU}0Wofo|&w~q}vKS4SjK_G3{5@^fKu(*^n|aEBQL=9c zF*3$}pao99-^h78WUTPc#yv@Bq?-fmQe)7axMH&2^xm%N=SXrN(0V{87f)N^89mXt zBj95!rrznz88XO_1eB=G#~j-ze(k+POc};&lStzQbg2f8H`)a6pTlQRncw-n+_HJ| z_X$a?+9jVFe27d(x^_VczyKxRozDV$OXHiW9MWR&=g;S0>dCO%HoM07-ws^xM2U|4 zK|B@AhEWGSk3>KjIB7-(w$O|sqpFxKewZXFy!={gl=)>S(lQzYJioTAsg;Mdo1_55Ck>YM&mh_`BW=$5(XxjBpN!eyP4A@) zkFC^V5Fx|}&mNq00u@#+(;a;%Y7Ago?3Z9ZX$ZVo?1`8tR4T^=`K|&3j!e5=C1iL7 zcq=mJw;JGKcB3^1f_-3a@;gLqP{{olJ35Om9Yr_hmXOg~bc3lexRg+hqvn{9p6zT7 zwkVQP64W1P6pESEsMsm@Ek=oJ7neE&Hyr+O5?xvS-2Qf_&Ib56&>~DgbX_d@uDoH4 zk#+@(ON{|26fI)JnnLuEF5B+qIa<;ohVwFf_!8j&)xS55l;{Fi$0`Wa)iMI&^0&YnT@Ts;MkQx3m^aNV&#wc z#inx?j^*$-MBE#b^L!?WiPYZ^WzwW+2_fEt?z-auai>8{1(k+^E+Vdk6sCxJR*8y9 zZlXyN>xDgDB{GPtR?7Vg74rfamgrhk8?g1%x9Dv`;Vqswc zM)nME>6xWZU68_AJt2O_V%k>!WV;WQ`tRw5?Q)rsU7>a71? zd25Vg0zJtIeN+V;@;X@Mjz{y=r)z8&pT!8>gPxyS^Tq5>T{tRJmYgwutOS*lrimHzk z__m}Xx=dsi4o=Q8_(ImhXtId{=1J=>Yh%9@=3et8mcWQQ}4NH$>;hgp6RB zi-1t{pi$yT$MHN?;t~>Y_dBfUXCPw6u>tH4}KDNY*XXgR>n{x0yaR;?^Qyhwt|5|e?67t){k8_L0MLOPB-evJH$ zcr2CV=NC~24l#aTCKJ6m)sTFB*;fwEdwP0J9S|;&mvG3~Wvp`Vh%nWM-EL(_6Nte8 zV-Z9FBZmp0$#j2odK$V94h#GCqGX2mG&@naLdeit;KB26xOXbdoL%1*keufP`d?rP z=mj|iggw;(HXzymL*9EpRe5ghqAPJ5w;J_N6nlxfuz-kyh=>R&mW6^yF9MBSk=C>B^{qO&tbI%xOjC=37XN{4)Rak3% z-}}AqEYE!AGi}#eFS(w{b(41Jf0&N`@`exK^uBtDlwj$)E0aHO1|);fJz=y3y;Zky z6a!H13PMz-H3l~LH&>ows~t5xI3dMVzO!^}-VtS+=`w&Zj)B6P8gG;VQbWTdz^{7I>)(h#8F#2j_Wk1A!!(aJvdVjzYLg^Dnmp?F7rEnHuq7 zxsQJ>{5Mz*DPO%KutPJBC+07be}8D+SPj-B=@JJJr-llJYbFqkDImclbKr$)=XPintpeScJM+rJ z>7vBw*RJRNYQyg7tK4$B@Q3fDPTplZPYdt;8eBK}YYxKFRp5A`Oq4bs+MBfDVZ`o3 zy{nV@a7HwA+c)3=nG{GzMm{T9IXLibhx|0}(%rd^!1>=EZdd$uKo9$AC<(a?EQ+jauF6xS zc}D<}qCmO}BpS&>p38QF=tRQFLd`X$ z(@b$b;i2CJNthZMR1!K6QInuV|)K1oUM6-w+avqic^eDbTX03W#yJ#IYjqHHt7J77=lz_5}$Hz~)Vl~95hq44eMNZEBK5GS7p>zzbfkgs&0APZ4O zQm;l|O{5;s;N#g_tHIKxY}d>YVclh9IpQJ=ql*9a?oRX5R-AwItKzC{&>;iURTbe3 zab@8Wn=NM0nETd`bi8M#M++df@f?vTC|w%_O`(j^WT9~XtnjVbu*_RTLLR!l9Gpd)y*(?FkjK)@?<6SrAk(MH{M~53ItVbVan9$<%>Ih{PvE@ z7bvFtvPB?pLT$b17$%-z4TcR7I$2{Ke1iW*E}{pE`<0Trp^!l8ilgAGJL+nt5D*_$Tgv5(#Zbr zg>%5^$1fvq6;1sj-GztYkbfKQ6c6AE}Y0leKQmS=2K7<_$_@MTHr}g2csiCV6}nQ zX_3OiBUj{_Gz9_&yZn`?2-FGeps`0;JnHwUZ7!p_1oc4cPc!=^)6aHP^vOyz4>czs zy$?cQ8IIzcp|Wb>6olYW#2|X!gK<*a5w4I=@ zI+_5O6cvu91YxCcxx2qhu_Bg$_vLu`vJgagwO62P=j@jrimM7Yr_RX?BN;pHhnF+c z8px|t>jwarQrgarahsX0A4d{MnLZK3a4gRx$O76r2%#c;m?>aBhwYb{n`4bq2{Ep? z1jVtJ=mw&8Uih4DfQWUl!D*xk6RUQNOA1}2ck$#ld*Vy92q5D^)6P!nu_mzsz_4+| z<-`R%$EH<7!x01blE@euV+J@iANpv%4LL8cznID+p|!s?fmm$@l|7Uf>j}lD*2bqU zAD$n9Qiu>(vm7qYDfJFSK9b?ZxbDK+Y)emL=^9XSqCR~>c}S9hij~MC&cUkx*0cmp zi`s%2ni{T@9?&cU4d_HI z)#0mA6|fT6B}g*B_M2@)wNn0QOy<+@6cI$aPhiRc3AYioFlrRKuOVa+emu$hoffBb z%U07?tec%@9|oqJZKdWgI?rZ4QMFj#d^_oV`}WmM@h<2ayOq!Rc5n0_xBu z{LoD506?tt3DU;_@HY1khD6WKdwPt$g0NX!l)4IU7vI`H8Lp zrqEyVq<(~;`>#PjR5T+=or79i@6@t(6sqzXN z!bWpBv2(E`4XA>A059SkkhwsR+RENPdHk4AQfi>3&w+B3D%O)C5Lc+^(M%^Y=5^B) zQDXAVHpn=@3-o|HbQ11zuVo7jRyhXHi_m2fv5xG%=!p*QeH4!YO(6YHz<0=`#nLRA z)%o$Dfd7EF*cY1_`9Adegx*BUN$&^}YZY@h@cYmUrCEkL)rsK;umK;SxS~^r6R%*W zPx}a*n&0PbAbGJ^os2|)ml61CEn=vn_cjsa5|AvVADxBi%%bu#QQ~99PgbImH4JAH zhr2CFkAv@{Z2;5MgPqm@_>&&Nc~IIhHHqOFVhM9xsCf>>R*IegH&z1#cxYJiHO0Qv zqKD$hFaXo+wL`ytwRsO3hNvALjFDQny0VP%-jTrLDIU}_GAoB^pvcrSD0rR#R z?lTQy3v<8t!W0a-@jriFQr&Tp7^kH8PcT0?d^J?{0HFjGB;8U#9D8OYphuzve zn;aL)+(C(;z7L$8`!5LW&Ctn*s88$-8OXo(t_bHVT`IS(pauNK2Ho7<$UC@fDm|4V z$VuyiKBf}T{WMeos38CS>yX;e^ZIAKNNkK0;6a4rq=cYOk21a}B&L72 z)OJ8!()Ju^dN{|3omYScVaEd z$5^rE1R|pi0|gPnCBy95UTPwyDFGE~MaLnOaA4%OP5B%?)pX~PtzVwWzYl6M) z^c`%8Dtg_`vcQIc?n>C#X$#TBhuw~D+ub?Cpd|0yAP&Af_K+%yB@xX9Yp>1BWr**8 zx?Mk&n|NY+U-*x*eYIQN);_y+HD=-g=P*xG4$-bdri|+pp==qhDBnXYp3$jz;?D&$ zt(&C9P%Nfy#)z1tk2B&cK!&uO0x6qaVsMdWK1eK3pFG>`Ma&s4oPQL(3zmW80*4_8 zOtW&to;)<+2WuF`-i{!9=Q?4up5swJLwL;QPzM$wv({Q3J-a(&U6jP7 z$^y^dCiH(VtzGo=!}+e#%Xcr|I__1j=lt!lTbr(*_0)Gm89QgWw=i;40n~K}V>#br zg0t8J2LmL&4jhIY3zRHLuhkB05h5H~Hly$qe`bTm`g_Q`v}T|TOhWK!>3|dE@gHf2 zp~GoF)$M&vvEbSLp9k5@)ZM!JDqxw`nfLdGjR79M8u|F0zjV~WDul{Vi-|))%Glh! z1f9-5x=!>-lDf40OoLq=j9M*f@tkVtkf;KY&*>v^Z8o)k{WcN;$6?p(+x_Z$Vn0~r zzW)V*J7eDvUTCA&l&YnOX(&Wf9LCHsHfDx=2N#2xEg4^1I1FEu-@Y{>V(9aLpFh#~ zRRHKi#Ld1NWq(zKu4>4sDQAI<4+$(D#xOD}B*QKfOCe1kXEdZo7I`2-@M;vI1UYa# zG^}hHV(o0WnV=W;1YWgr=p3?=>FOngrzj*jgj!g z?39@Ptza{d?mgg4GJ?Rc<2dJ*2yRpdhnh}adOp&45j552pP z-@?WPRrj6r#t5Y)I4}w%IE&Logb6B9L!pbh!qUTFoXik!<1CUBMsRDGA(OX$Yo_G= z#5v5KD3o~$;XW_~9|9WhHf)CGR718dA@}3`fGr@KWwEXctaBCQ2`U?Q(0f{*R_b0m zYCqAK+@A4L;l0W&WUV!?Mo?inYVFc9y!Q)f!y}G54~hf@)@cdw!^AuRf(EUm(+LAGDtBFpDFWm{K6@(j$}hjf!d3lt z|4JN+Gr*7|DL{Z9J7&G?K0P+LSsF|T!N)t`6$+Yq2R_{%8GN%vQgprQrzVs9Wgzxo zl!Q(hsk2=S1I|qG8*gR_qEw@!jA+FL*=@Il%~F~dj%_WNvV>TG4?E;w-zA@zkHaV{ zJBVyCf$Zf`#G0h3jeQev&Iwh0^jNbwG=GqGQ*;H*;6!IE&9TtVcg{j!?%n>LNGmdU zbTprX(mL$FSZp}T1X0`mhC={W!U#gw({@kS-36l!KO?C`) z?uiMM(su>!4k?CK8EG{_F;HCW2YS|-U3M`iHB}_ns@R(U)!mz$@}f2iR_@o*2gr4f za>19GB^bO+yB$n;6NxEg$++sE>*1~i+)Xt#0ilC}nyeHR`t%URBYG>g)CB@Z2`PB% zG}5pB^4YZ~bXSr>3#AU=JFK%T0BIy*=31i9Gf*mOdaoJ7_%=FbW}0k9~a`2a@!^p!2J)mAYS6nbd^@OAJ0!hT~N` zjl!%%js=Fk(9k2wixGe7;o&HZ?2B=s848=b-KyNa>knMIL~OQO=Ysz3CvPJC^{N5q zusH}KkEcI7kJ2Mis}eehRt9l`RqGirqeV0{TJHz9FP)EzN(jX)R$JH|d024xY2nik z|5$UalTA{;WU$yAQYpc88d$7B84Lr6Xi(Fx=f4(GzKI|u3Q0gq-K82*ZYH3U5+Y}Z zBrZvkSSXy=0e!jRQUx+(q~_Y~kP*ZbyadU%OAJ92u~U0rU&1Wd4hb22q;KczZ(RNv zP(DGhvIlA&)JOw(`y@E)KGF`Bx&^MfP&0h^OO8JiHfScxA?US!(E!1G(&#t)`vly+MabYNGsq5xV95B1yB;Pe~}2o`P+gJ%5D%YoI?mkqhZjP z(01wdBGU54cGQpk3$S=B*(5Tgh>2^^x90;0ZMd&l1Byh6Kz`pg!3D`@%k;Hn=foR%z)ve zzu<~U`t-!`G!PwxdU;%nROD9#KEhBu2e3s=Xw*lJV7I3x>PD%r&5XY4=anB3z@y_N zV_7qSeAqe)DFvVfjAMzRWYTNk}jqig?9OyhGFET38R1wf-sR5_#~8~R!d&^W{z`#lxEoU zMY&feA8;dr4~DL8Nr?|l(Ai-#hZjD|-G5gM)-}BfEl^jXd24CGwdWMMd)%g-8BIJe z?4afek{WnLI%Tju$8^FAb)Q)b8ZXjKZlGFOZ^oHDw!iDy1?5UD-OT@7<2(EUrQ~~T zhi!P!Tx#B?fg9#45M}_h*DRxy=mE)+v}ETY7<@Q)fnoB1%Kxw2#fEq;7Q`~3uF){V4&zB4iR9?YQT(mR61#`s%q>Cqp( z9HwR-IJ`8_RvR6EH?S)E8)bo${~_@0rZZRPUi19*w+D~EJG1#K-<3LNJb!z#`LK-N z*=vW-O0O>b`q;B`ch@f!(D|Fwo#<`f?r{I<@6}g#4Q$9iD7)3?PQ`&Ludph+iVkfj zx6!7TAvSYu{xnXQ8>XQpZsS@fL6=)A@)MqOgi$C9%*|VD7-?RK%3Kq7dz?!xoZ}S*Cd1Re0Huug9 zS7e(W7jc5sd)Ccpjn+wo*G@ik=n^t>t#`V*nvqFxSKh?@|Dh2GY#YJTdyIKrPdLWK zQCRUrc)<@YXg$3qRMAxtlrgyIyG0WQRf^qNT5%n11Nh&yd1YH9UVKGY5h=&)1y6GT zXdV~AO-g;VSUNsAAMvWkd!%Fh4^WzC=({|d~aVI*wpqIB)bRv zE`*8|7(V2W3U+7E%SmJVMkD`C4hwc%4JR0RSqbn~8apb`i9|?0rVB02fk?ej(PE>f z=!k2Jovg@!>{jHi{`!Pw(_!Yjf+)C6wJ|kf1E^%} z!(ZsOgmo$9EG;b+Gbs8;3l%inL}k#k@_cHenj_LHZ-nVGMuxe zO3&A8ylllxu4{&9qW?`G&TLLUI;YVU9oh8Kq1CP-5&$2m%R#5PU=Ji|@NXI@ZG<4^ zRGXN%v-4}-rF`tMR0QQY#&%h2sV0K+p)&pHlwC88Rjz&*wjQNifoEXd8DP?h^nAyu z9Ej;@q@HDYo=5N7Dk*D2%d8H+VB2cWSf^w`08aS*`1l>PM0Ifurb#q^pzDb2_(`A< z${%$)0bCR(xlQ%T-)Zrc0P=bY{f1J1y3ZlLM0{%H{_~{G4{Eb>ug=R-47JSC?){vl zboc&sh4uRcwFVtKdXFs^ToN33+$Cq9i^wA=>VzPTXvgf4qkv>;3fyLB@f^{%u+WWT zZZZCP!^~|cTyB9h*$>l${srwo6H60w`UBzeqpAlQSi&b56Jh!S=D?;wcV@t9&w88k z-EnD^fM9A61F}nin^_~Tj>6=I?w18_*Z2MPt9eSvmaKN}u?53jAr5*Y*w7MY>)ls7 zfHn|?)onBmvVcUUd&`AKiyYB9k%}ZP zm`%bk`M5Vs=h_x`AHxzQ&14x*X#H5Rw|NJk?w2&-9J3#Vp(R{$!)@}lInEcA0Td68 zwI4?!kdQa9qTu5p$PQ}VPQqjGD99WEnZvSW??1PS{H9z zg;492Hd$nBY#Qg|hz%wNg+w4D+e||2)S-mwV^b+*;$b$8Ms!7#2$c~lWv0V`8v^w{ z0c$XXRif$O#i&#vv*S*_zX6qx5^u81iFHiw z55``FIpPG0A<{$nbu(BZE_8Md;O@Gc(t1}H^fteo9#L)}5RGjnp)^;2t37jCR*%%3 z?WIG6%EYC#25pSkvPJ;cg`KE?i%daEfXlBh0g&jAn4?RhU|fwlpEquG-+vl@v-q=l zYj2|WM3f?VMd`kWb$Y!8(2Q8PF$!T$7miptt)XG@%@!;xlD+R-r+W2Q^i*LzsXx9< z*9yU-6#j&FSY)*B1zGzH2BbvNudfm_JW6O>F1|ntb8Qz^f$Hy)g%m`60))M7pH)a| zzn>abpBa^jJqE90Q@rs?JLXDNAyNtkJGK*lru;=585y&!I%LYmWy0(lbR{_{U~;4l zF8-MP;2!yb+T?lbFP%s8YxQi-ebKY!+mWxLs-{E4?%F~YsZl3=5c(>h5OzyWS6(=QL)(8<-UJ)H@ zB&_ywqHS~#gBA@z)7PTy^|1xGvihZh+?F=+L)A&;0RgE{H^B@e^VXeb}&%U1HxW>Z5$F7kRo4b~aCbY7R^3 zOq0DrzC+)NJNJ`>tlm8yPS^FUX)t{QE-Tj%W)t&l4|5}c9tI#>7weUUFAOwY^l@fB zuk`4AjXa3QMbs};)Ya2lnvjbupHAM&nI3N_*qPcN=k|``0U@Vz;5ZnZs(bmtEf@5c{)S#P~&Djg91Snqw$6qIQ`oYun zIui8|(YBh6P*8m%^z?(Lcq&$xD5>1{ge8;+RU^yC4!I}bC~SY*N|}d za&jvZP--5l%(gPTzF|;#OXck_L4gvkXM3vSRk3$J?v@(7|J7s5B`orwMT_O(f*LJz z`f!ohFrr>aax0ALDU|S_0tQUIgNQ>O&hTb{F+v|MRCBf_s1C4L$@hU=hu}OMAj}ZO zP3DY1!?GOKCR9PPbq`I(Bc(BH=h%1Ut4xA#kq3&6ODu&=_ebV`AG?TR0vl+ytQa`^ zrJAI2K6<>z(F=nG&=@P%a4;h(C9%sR zO#$)TLf7W9PRcVcK3(SzS>K&dQ{vEN0)Xv+@2S|C^Og%BKLVbYyn;O#F%pyDG-*n~ zQP^V@(m(Vzh$87JrCuBg>|ni$VSRi6ua)~1B%G8hhWd6T;=l#cl2cKT?=VK8VXWVwAZJ1A99Cf< zIVJ$5{s8zBs%>)KrW}Jz`EQYlR;d+x;ecn98*u+69f_bvM5(H(##c}LIP>&6r7OwH zoyDBK^_39sa?_g=!YUlZSeu+(memm@@-jnVcKW)n<4ViQY#s|}w9^9UA^2z+roj(* zYE31ej<~CB;hlqTgN(1mBAt2hq?EO+0cC}Sudy(Rwc^)St4%wAPdDVJBMk(8Ej0EV z5vwcXJfMq095HbfAek_%U?hU%G1Sfnm5jc@GKM-n>JL?B=k_h{{uR36yr>;qss6`Z|c5YEx7)3BGuECu|VA1W(fsfGia` zeND0Re9-o}n|!?Z#lQY|O9mKvy{V6vPp`xM$kn~&-OulBt$T}4O517nPlX1hmC^(I zdEIA)L>*9HnRL$Ne$seVc2eT|$ejM{xc{&wKF>7V=?(#Fc`m4~T^N^KhU*Vtk3llo zy6JpWg6#}$z&VdGS%uox&+o_#Momi_zxKPPqUY!1FDMeYr`ThjLsiJu4?FCgSU2&B zE1rNf@E(tXgb)KU{lhdS5S}V&*ZO@IP6MmSZW*SQOkfMiyWB9QlHvj19xzkflO&Hc zAfGDKu5~IKu5l{MyW08%YljMs=fc-0k{BX@@V%a@;|FM|V3#GVBDz|2x~HL~?w`?i zId0bgEU&K%Js=u(hNp*M;udlgL3p5Z>pgp2hkYNSxNvvVxBmuaHsXf+4MT zaHv!`Y1tE0fpvPQu}rw4Als~V(8u3D3TTtZ7(wwJ;`o$rnB``@lc13eZ&S*4dW?aU zpw%Icz93bD-1yc&Fx!+VYYtP->Vyx1CWu#2E)}B>V(*u37Z%p(X##w6_1S!97X(A+ zV$AfX1((+EWEKD`evP)=s)Np#xQ4*SrO+zuuQLJWNjtV((ucQh7;^SK@({W~HtkFb zU6YcpklDbN5h({qlPeU4wOl|~kA3P}q(JE>6;Du35qKaT?Ia{5BnF6sGKuP8{G*7t zZE+xzxL|o!BUL$TjbzaX!2K#|dviBvpSJ$yJ{pFpq6OC0A6R5^&p%`0vV5qzWLA8e&s(A0M z8zy@$L)WDeTQ{d%_9#S*tkDHXO?7uK^hp+sRhYT})hp{dcY0fBN{ZOby7lu<7qPZ|FpZ_i7=zj8n;v$H-EeclG860@ znJ9_~7GG%vd!e40j&ijLMJzBOLon(lsN;D!QPNN;)fL<^qJhVRBE$hV+|L4<0q>7l zm{LgJY0fTIAPgaq36M~az&T3}>_XzCzwQnaqD2^%6M*TG%X1ce)FzJmUf|JtEL$#R zNu;bWd-B|lSt`MPvT`^C87&i_oRnbU@@t0yl@&xt;_6f0CK8=s*|SL zn$;B=gap9MV+v2B1DWFq$a@cg_j3@Hf)Z>@sQSEv{;zYua8$6hJUoxsMreeZ61#v% ze3- zDw{9IDjlrRQ5FRTm2jFK1*=G;x_}2&;5PoDa{|(uQ9cOn@c#TM)wRw1n2# zP&T0f+cFp#qL~Fi=lWAxS^>Jo?)&fLGVl+xv@F#_e3%h`-R5|mR}|m&*^Pp+i#sp9 zxpqAGpeQHhVwv5b{)9srjPEvU6a1#qYMwo2DH)oyN{YeqZy#_KlU%q$!|>H#o?bJC ziWJ&Vi1VNsc&0m-Z%>qdpbws&ESguAH`BG^U~zKB{o2z<=o16~7Q$SS=bJ!-*XnmW zLwFFHd=UTUzH(t)qKb){E+M3LA+!IUcD_EkG%&qa{onq53SgXXbmQGGC>V9G=07o| zADkSTzd6!ZYmha1qbcV5%T`&gi|nZo5gFt9fJXrTn4RdxM8djJR5+ZQHb{%sA_#dS~{?GAeckWi6Eu|m&9AKh9 z`Xej6ogXICUj%vSJN*%lfx#@9Kg*Ty%@zGkK9qOQ^ymNgzasn`bTNZdQc~2$(LP&f zS}$N`&KF-E1wD*nT(Gg1=Lj;ta$4f=v?1)8P|^v5ahK*;gG)eC_=o!c z&isbvI$D3ZkYR+fU2uKJxD<6l-Fdq;?Nqp&9tLshAxeh{nj*$fia^(E+$pr;JhtR9 zcQxk>Uv0(j%8{tqW-^On$SO!k9(|vVB!o->BUn0$T3U@l38Q40z$gXCPu~cxN!&dP z=F7r=VXjy3)sOSm7sGfampkd0KR8sm@3q%jWRAiDmsZQVI=hbMycFa>nB&G;r?(Hx z(lmJQ?UbUnr)`dyh^~l!hnqhy-oFhNG$v*GA4?UBf&P+4d&{NQ*Wz;5EMwlHbeUQ1 zyHquG9gfrQ0Af=~up$sxu)-KeNe{{#XqiV7ez5DvIy+S;gQTO)df10*0R0L463mT2ESm*rBNdKqMWU{CiC|g03@`|BJY9APyiBHQG4J_B=kd0k zhMr#N%c*_`k1>=thp-sxv5tTKI8O|=oy|dFTLZnoKywDAXGI$`8z&(d!#OA`_80@_ zL2G4Ih_ocSKv3sFT9F0(*esVJO|il+*NNJbF%!Or-==h7%ZbonHpCWBEEuI>@u=|& za$w|BW~L`Z9YVk(Pz~Gw)X8IW{|G1sU#{JE!6E|>7_a!byzkzPhg6h_&nk-#*l%YF zkS)8~EFZC-$9q2n?@f8k+OX=%$hXQR?81bI@?fddU4U3!zRKbK&@^xzom%ZzC)~zW z407S&jG8Sn5OP630Z-Cb+vLBjxsC-z7A#iXAC$nzoh--hjx+Nl|t)Wo4yE(%ue>lNPAjgPVlki`8N9r6fVSB*YZ?k(7CC1x8&x(Z3 z#Nn5`4Hus^$ipX@_b3T48@I0iy`*fV)Do}C_NvaUoO+9vZ=&3Cj=8=_Rk^r^b`Fl$ zBHnNI#7`XV{J2rn=r6A^E#+G~VvZgR%Gfq%!%FYXR;y3G9btaiUHdz_@#ww(bi=Y| zCZn=mXN{UBZ)ScaJe@Iau@w@Jo~UUJT`?G_mBSL?C1=5SGS#3H#-hFr1jmD1H+|z+83=Wg*oyi3+io zNAav0uHr2`4<$c~96mMYMczUD@iBYf+vG>pXki3S@bLNJW4aT7jYCm}-%lA|M!{5H zIb0jM2F?+lYz#URRiMnm&oWCE?F2!FV8GPvbUFZv_#lPY4L*B03?P?aJ`*ZD`D&ZN z2{}SIhj2iz)e=HfEF&wA6LTMNCT(Z(n>fYyK`kW!NnEt*0mQ7KEKe5j(eaS|5)@Zs0(k&s834U>!IMiXvZL~l z*UH^z5R0SH7Hni7kcUR{)lw{FiD%;3+xbj27!-51N=$rJ5l8xeY+%Ae}9q5_CVAi=ue-IM}F)yS;JxcAt|BRD{PZ{MhX2a=6H z9Eu4A&p%{I)ZdrS>dWqRRv6BsMQyTEATmr{XQYwXLOHNT54L%1RuhC-zldRjZ z9&5kzE37@;9P$%ClrYIo;aOe)fl4!mMzLstR!s6(@Nbu5B+^00VwV6rG>b6dZJJ{v zNI7!B`8hThh!zOS{$VfoxV}VPslw0moo9(z`_2^*9*RlqTh!8MR!}?gqWGrhc_p^6F!3r z+xSMcf&w3;?%yCVF(X`u$Qm{*q(;?BiJ~~KMN+fxjufuC_ z1^$RGt~lU4!iu+44q}jZ#d;P=%GNIKFsW@Kh^~IlBZ=nEylz z2$WpIYxXugCudqoWjvIUmUrrhevVI!Uwx!O;5lW7Py^i#eL+R8EgBX$R8ZTQ>WaXO zuy`o?WH_Ci)VZ-nEU~%A4z!=1{{_+$aL`W?-_Lr!b zzD4U}@INrAyDe+IVHH z^-)wNy(z|pp3KEpZK4Rf*(GyVoVrO>(hDxoD9i%Ol7hjt0PrM+c+p6qymW_;$^u!! zJdufnP??C3-9syMeZXz#KsVNQ`ygDrnN)V->O5St6!u#UdJ>uPBla@hm6waaj?ONH zwUd|+W)*3I0+%;p3K5)lyL3ddfE*7}{@9dE%QzG1*g58Z9HmCR(rRub1ezv!ap1}xK;{f;|j3L8;JcW_jk1fHTFy+NG z_^MiWG=oO%az)5@^#z`Bts+GU^pXfC9P&Lb!}H&MaGDUOi0sC$l)61C@NzcTI(ZkB_pyA8Mw9jL%- zTDHOy`RGsnwL*$EMg!F$)?O(5Oyc<#ym31(&cHXwDKvurCkASQd;Tzj$ zvt>@L4V6Wa<>_}fA(tHj|6~i0I;(2J2?yF_h(L={Imub&0$HPE1_$QUB`UdOAve)@ z@tSM}?Ev9;5~m`3ySwSzdFy;&pWQsF@Gva+hz85qA{V=w>Xm*@G}hb^2hu+~A3`?r zLqzz@%{>6zh{L5rvUJOF&Hikt$V%Y}*^~(SM_|2$-@ZpXov!f`{ow((8K)StvIXbx z1-d-)63P*tw5h2>(h>g=aOSmEN(>#FRd`M*sD#JkG1Fv*bfgMzasN_-#hX)~A6_Gn zgB&)xt`s<|kZ0@P(jMXbR_8g57CF;CHI2iwY+cptitbnR)ZW zT&1n2#LZ9nO1P~YQDg#K)p?)Z{LygaCm}QA%9&gSZ;Ji-(|aB=_*Cb$-Uf=JcFK7=t+&ClTBWk&~P*}7$K;xgl?<%1|+zsC0a!WTKi80XVwR&|85r4KP|HiBLi}N(nKjgvA7(lSB&$hh@555#EZW zC>g8Xrq*NW{>-zK>QWv;Z(~A%04W<|#eJ{-Q63~IO4W?msWm+38ZJb3R*>e#O1&$b zlZq&lN`7VWV_;MbUDco51mA4JZ3l&{0tnY_1kpRyCXiXHFOXUX%05vps2(A2hf8xN z42uM2I1Ky+TOMb-hc9!ZdZi-_AcWlWU=T-8iFZwDxv9c6PyYRjtk>72tHFRm@Qkn^ zAPrQxon9HR(v!%b^oV39J&L*&l` z?2uLM@Pe6t=i@IrY)JA@ zyVAQo+Slw5U`HGTg9GQDg%g3MOm_&2EAqDxOo0eN_T(&L2WL#LI0Bu+{Q^*=qAxsk z@v#79dgRd991M6PvfR$007{BbcRVrp_`N+n{Ao;vktU)SVriPcV2HSC?D-8>A7p)d zIVn^SL5Y=^LmJ{6crXe}Sc$RU$!Cz$6{Nt!B_$h!yH$py=oaEFhk#LZ&jq7XYCocx z3V@P>IC>rMe?+3NYc^<~xj5LRXS0LTyV`pStF1Q3_ZJ<%VMxaea$DCjk3RvzP6sK&vh^nGiBGjuQ6Hf&0**-I@KjRCipbN=$F|i?l z8;XHzs2Phc*eoI#68i_;xP#td?QYXuZhlR3dk}QUgR(_@5gKSoH9IQN?*V^JJD+#@ zS?YeoNHvYYothx*93xUBn=YPMfIIYZFU2(JBmDNbSwoN-L1h6)T3rygl;c@cgD=W1 zp)r>TP_xGm@^TcGIu|}pPQ(GRC*{fjKZ z64oGUBuZ=x*3HALV0t*I1^~Id9lA`IRN?jvB@Nr7g1HzWZyqC|kMhJ@WMi0~E6=5< z5XWsi(3MZ?kAJ`rV%#H%Ax5QfW$YSL@?ezPVPEYaM*{T68zc)V^Wv$X_CQ;nax@fTw_aIi zuSo%SM>esqVKRe>bwD&dI#;)LDef)7kDIXeOuc-s^It414||B?@-Gwb<=AYiM3H38dG zmv7;$t`n`A?_}<$?~-G7OPaq^tIvR*@4iW@H~*HDYw0@+n&pWSvHGE519^wmPpBpRCBgp|!Dyyz z`*D9SVqC}Rco$DlYM7VL0R|c8{pEkd0xR#r|0l}CN>5?I_tC*bmjKsWasT$0NaCrq z=L04R$aDZOwao^`K?PRV94!l*zjRCdI}(u=87?oKhW_yRgyF@k+MS;gTZBBs+>}C2 z3w`GnQ?~}dpF{)Hm~hIE1*=11{2(jbmMKi)D&^J@`)BohDZvse*P)=}!TT=yhDo>I z{GM_>-`Ehuz(DY>7v%yMikrkcEoXBe_x5Mz?pnLr zlVd|+Gr84Oq4rJZg>@(JTMz!j))aD4RFhFNTgHkOpG_3Q7dIks8^S@41fcDk2mzEK zv%Ym`D64|9CsBzY+h_uoUfeQUZ2^$Tw|l_|%uol5cH115&%xR1eDwD)dK!HL`|trn zu(4Y?iRiOOc$lf(HUvUD6`yvLLqkVjWZ*NF>6L$$Zfe(44>epjnm%Aa31%^MYS3`7 zA)T@OiPi!2<{YpQK5+u#ef}?s$aN!$XlQN6vTQ5AvA-u8?JVI6%az*Ur-M(0 z%zhwx6Sr`P4+rd2jlPJnB+Ta4Ohk+8KDZ%ZjX;p4dgD`Lr$Y2)(k&4#RRFZn?;H3V zOFnGYZX}}iVPq(3r`a=-mr<`h_=>%O#47LI9PTO^mIv#MO zRy&Hal*V~UBHcvz2hhs|5rE5caoUW~IiYXa2PVHz#S*FZn?+wapt9kH#_UYx0V;8d zV-FRzJg{RGYJ8%92=C$5IP4;O^+Dw;fgx<+1hA*dl&{m7$J%j~X>!85 z!mt0|<6y0QA|2Sz#b91}Gu6_?P>))T!?Vq0cJY;-{`YXK9&@uSY}eK$uJ3*ZnLc8h z3+Or(e0CU8K4P*%_F2`S|HXdRX~$3ShFv+9!1RL)SOM~f4m~!N$$^JY7)eFz+tt}y z#!dfch7O;gh{$%2Woi7W_ct=c|7txy=6vxG>W}P#X~l-OtrkC@$@qe$grtgDHbjY8 z1HMa!;+*tOJl06e8-5aNx0RTo%k#hKevt;gN|5uD#s zCF*?Bz)B(Jc_^aupJw-Y?WXVbLvIVJ)dmk0Mz8Lk%L?BzyIJH!#Kuw`)#{H2?k*tT z9If?f`L3F0>Up2Zfb(4T&wu;SYH&};EKt%?i7=+x9kWj}`q2KQ0?WnC`BMwbT&F@3 zC2AM3ShV@gMfes>1 zmyd--Ildkr(;5%~M(@7A?7E`FImF;vXvG5k-w&oYDcK(N^~Up2A%r_nOBpDI`oApM zp!aOo))1Z9kY3Y=S%O+g^#@RSgvV@3u(z-_%_`z^x5 zX^{1GOt?Y8@SbrXKaya%e~l{B9IE$)O}1|xdm@>C3y7(UM;c;UBQC}7m8i{}%xpt; z5)6_z!9LVaX4NAWg6{acb??8`5Ji)wrUu-SLJ#xPTUPiF z^Y^QzLV~H&Q#i30T@FFi3geMR1zT7dvis&ZHpf~~)$3bmHXA(h+MSgAA@jf{9u$U8 zSfdL;T?@skfvqLz_24jFeGuA1i%J!qMoK?u3HDY{>L{s=QQ?)25=JIwFaNXxm0Zjc zv$jy94iys+`)X)m_}^Q0338v*mOvOxW2&mIwavz75wkX7WN7yNO5SUAH590QGZ@%^ zBXVP^9GuDc6#AwGhO0;r-kt*Bn$|~niiRn2*+7eQEujVQuk|sCN5KxQYX|m!obC+- zVfE~T=;~k-cZ^6m>=Gmoz98>YD8S(o9v@&=hFYtnM6ch?+G09Iaq*J)HD^38O|)q%+h>NjwuWl=4+HYMsNSuD6R;#OtNr zKM6hlfv2N=uQ!$h6C~yaS=Cl{31SGn%z9XAHaI>Cgif(bC;^YG$H2;S-1fpbOX34= z%w7t!{mu|qcngP0ngr{3^io9}CEeSAJKC|qMb({Vf1eD1Al&3ZN}PzIVh*=Qm7Y}W z5+mqK6+=rU6d)#Xk{?_0qjFrPPq@l_)SlzS5OCPS^U$H9vOUDy+ClZM#B!hFQ%Ue$ z;ZKOaq)((PnAbzTcBr78KusvFVB_qSM9=+XjQ_uiHOCPX@w6IMN9cA5xg|yf2=%Vi zn?ltHeQ zV@TzgmuBey5bSjhA!jCqkERKm@v1IhN>CRas=XJpB&d0YXc2&HGfD~I@W!R>1~U9c zev}%n@p`95%xSQ&Y8PPR>k75kqE5}PIhA;=5*I`H#u+4U4i8&$8bPdmz#DB&<^g8t zz1;(xVp}<6Hv-70k`yDJY#W-6av_5e6CGU)3(Q6H;tpyxq~c)Dgc9$IS$>>n^WC*0 zabBrnMV&8HE3VE%C`%OoF$9Pd3>u@8m&q(}A&8hnZ1PY<-%2nW;)6^gP;D0#O|o@^ z1qSm^pPF0PB{5N`#MqNnJg4Z>^K1eg&?vEt+V*91?Q5Rrrc;Ms;Dg}sez4mB&EvVDHYtUHeXJ(}cX zYZu^COemB@HCj!}2|b5~VDt!;Z@B#Cj#0fX`pX;9F%)|$2_-ctG(9gEBEV5_Ki`R* zwm|kBuGoZT#|5t`1*;n-_B)ivY_B!&Wx-E(z|Sjf1!Y~}2};mCp@c3juMI4OnhdL#=-EZ2 zO{mq87yR2)vY~reN@jBHr9Xe|z!5*hSPkXRA1`}wW=4MT4 z1i2C2Qvn+M*d?^ZC~g2w^vX1845{>9+B^o+$VYc$^R}-ZzroQkC9Uw)OA>rY&so{=tDI7qH zw1!+;q}Q+G>N=j`b_%a{O!M=^Wc@RLCl8{$6<&8)HwlyR~r=Fe>~ zUiJY3ACgZR1;itZjnKPqY>uOd<|wtQr5Bqjl)JxA*Ri)i+WYdZO0A|97Pfz(h*Qu1z^MzP#qZ zS%gz=JJOBth{=U6t&a4KodWNzc32ma(kqTg_(|yk3qX9YNVr0F^3@We&e;Gl=?XMUWEw*8Vbrgm1rsv_3g|xAUoC?94giOHY&q-pS(?}P zr`;nOmo?afK$zpWUCzF9+U6|rd)nQmTEyS3B{0*I(cR;70Lpjd! zrXJFXW_J(GzUQ<=E%QK$!iGk{O(YdSIU5CTkn}yuC4%Hwul`IR6b|m=`?CK0X(5*l z3XvXt3F=j8w*jW{U&2aJKPgh=MZ~oa121Js^=v0X0EBOX9M~}~n_5zUS8)-6?dMPr z9Z0#GsJ8ihd^U$>smJ4fxe=IJ;)t4w7)F2{fcpwS0KqkUv`Dg?I)cFiTv2p_Qw>>w zx5#oHn|{ilHRs#yeB(oVt`TfI;t{OK;TPyZ2o%|d`>i{V4nz7p68AWh!FQ-~A{5XO zV~IGjiu=1GoZ&K$qZ5x%2~@Q{^wgO}HjEDj)QwSFwtU`e>c}H?7Zjp4qnaj<0KNuM zGzA_M@M{dn{>}Ysm90N3O19dn5Z{Z)q?9Bmg7jKMvwTH~twC+@e}J#Yytb(k%s`#v z7U_*!kdqfLl5DvN8m$!a*Vy-Z1UMU2KnCPC{ha8vR5L@jiwP>fp!zvjW6_wLt|$7Q z7-o2^ezP)s91|LPiSkX{W3*rCP*DnKwh@93U(3Y^!FXcYup*Lf-^6w36(wpo1JVn}%XS6Z48Fu{BabnPe5fT6nrW5L6FBuzdOy*b;?|E>rG)Ub7y|vw7x=w=``Ki^-y`u>M za|-Ag%U!8UyB$**^g8UQ9PXPq1$Jk@$Nhe|txZ z%<`Y9f*;vS^})1!2LMZB00;DLZ@^<0j-P?BhtA3k7ytMw3)N2lSmL-15lI21MR_hE z#ZqVp!KwtM)k?@!Y){IL7+DD(0m)G@tw%Uz#Gzee$B)$HEUG~vjr>kg;uY*+DZm^q zn^b5ZFqYQRM7(`MS}6K)4%iqaI@P{awuA`4>;RM|aaPjPUf_FTRKAVE$)?V6wfa2R zQ!P+4&B8bH?@d$cQOY1duu=`UgpzFgFgDusmdB7~b5?g^iFT z9Q!&z>rD;|`Cn85vXnJ$Mxc|r)^J#oTgIvH9w?8;*p?u3S-TsRbKZzInLs|8A$zqk#-)LYvou{~Ym=Hngwcpd@A)Laa$yN<9pQt`otC`u=0O2-~A^#snk| zJ`THL1)x{h*%&8~&bNV{nvC4`*t&0{W{#)k12iElhp7CWAm^Q!8uyhQAju%mn!M2( zw+T6=xkCDEtmNxF{34q;LrA%_yG=l9_QT*UHb>88xX!st`lXfuG4GXEHb8w3+J2_8y&aAWgrJ`ke>prz4MIp z4GED5^@btzD!N2zZKRTbTud}f3vj}N$ek5@?dX3gZ(aRnb_@FL1un%DtIClK<;quW zKvd!Jo_c=>ZS@!?O}frhRvIT{!G4lTK%A14Vu@O$7UmNzvgKnDEMF0XUI4Wn4Dw@h zP(ua`bQ`>|-6bpB)Ln^!mN9~~S?KNhG{DN8;95QxDp(S}jI0$&zpT)^bf;SC{!BO0 zzzXcVFaTFEU+}3BP-B>c;ihq<-fYmGsh89=!3LD26;FuJL0St?r-)Xrpd>YHX0pp; zjO6x^ugBL2^H4h`BUlc>~AK+s|^|42TK!^L3|g-p%~a-SrNNE&Uh zS{a)hPUQe_+XC*KXgLs!5MsUjN{nK9sotRUo?`3&EW&UeeOI}k_~QjdX-f{gD?9+~ zmcvCZzZF*~tC3TMhUQLa3vDY$uVL7~c;{e*r7QaT01i<4O%sXgTsI^e8;aH+4;imB z!(J|fC-jHQ-HC24LjpVHM}>ui2!m;{E))5VPwHUpXuu)4g`O$SfElH<| z`b1(7q}&5hyzOjrNg(}ip=>b*Zy#n{5h###jJG_j{OiOSa7K?Bhl}*O>*n{1weAfQ zoeWjfF2W`dr3I7ICvh5Lp_^yE`i1(vNQ021yeu9rQBt+*d9y#(zpsY>sVwy{1nJum zEYAMJf2Q2%@kQe3r;m4`{ujL4)=+Qa^kF|$z#u+dPjjpo>dgb4Ij2{Bg0>4ADqa2? zd+z~N<+-Pe@0=845@T{S@d%bAZhK2aAfS|JqHdLOve?uJ2NwUno&YFKB??7nyXd6uJlBc1v?*H@DVt@CIbJLvP-JW$6(?45N z5rQFkl$LVW+W=BOCY9r(Uko~sRFQNAo7azNI$GO@$`S@X!uTcMtz}c}Fk4LSvR5y9 zT{6ANd6b7Xb)N8GqpA7C<}-LjB4ZwN0pr?MMgBW_W%eEIH#PDL+aFMMCwUHLf&V(t zc4U0d6iQ}-cF@9R64l#k^5Z+ZO$WZRZr8%ol^^S9Y}%vNSdg`k)utE-Oa5bu1<4mzr&lnt5U)y@x7xb8bNxxwtAhGezO=IoLBEYK&-*tnr~ zj`(9clRU-40r>f`bq}w3A|do5BMz+T8K;k-#W^6R%3n#&sI$EKfDJjagBTH!gj;N~ zLJxYWJRB0$#E+fTgM&77Pr~>jL z^N0Jys|-ksLhB;ms;3>gP zxQBfuh=3DAEaene_))EZLmlgfTf|QcNIME-)ap>yeNE5@AoHA(g2#|I-2mKom35X! zcs$TS5QCnvAVmDN^Ko&A`&G*Cl|(5lB3|8W_4oiKz=M) zJ=Ig=KD5=pC0C83Hu-v}$mUVM4A1Y~bhQ{=Mk2Dj(C8n#tb)C0GU%C?;o=17BNC@m zuUd&?(c>A~qddGZbZ+@=eD8;X-U4M+WGhh#N(JFh;3;)j%mzMjhxoss+QvF;K%9ve z$Yz0{{D$);^GYp6Fo<}BIm;py=)PeNFp-)-WO}PS2raKgQfP$qu^HF>ZlWi&Z&^um z^o*mqhNJ>cI5^PxR;xXgDv4&ZM1UE5J{I)a^_g+6pNAmISxTqfitB)fYeD+Ej< zf=LP_f-$Wf$Nwk_8iAX)`2FI!DDx=MnT{(Y1i2dhoR}w9I6t>UN`XvPMyjhITnTDV z3;se4h|PS?`5Jszoh4$0x#w%|M-#A09vOt(PM9>1MTi0(TAcC6Zi)x_Y8R&F6i7lQ zRxsP|@eh{UA75|;R4b~BdTxFY_N~FN>B^eLwWAnvb!5I$A1RF=#6sPHrhCfekZWVgR7eY7Tj+6%CoIdfnV)cN&+Xlfp)cRE7l*64%#Ujs?@MQukB&9hBu$el)_< zBKhhwK}>}o9qEm{XcZVBba69T!|LNGz}f0ouQjuMKN=M&A@IKz7$D002|=LKV)v~x z%$e28r$;$MM1i%{>n(UU$G@!@xM+z<1bWf7E^w1@7>8tUKBDjq44XpMwBUi)q8cWL zWA%|t4#Qx>xSUfAqQcoVLx;*8ZH#w!jrKvNy>~r|T8d6N;^pmOn*+Fe=acS*C^_e# zCvNI(53X-~k8WHJ#2M6p;+6_pq?GYvWM=_PL2iBv0b%5{Aji_WAGyN~(I5(pozEJ_ z3|zNu36n)+9Q*mgNm69?#vug+kXOm80E)g$#v`8ObaLYmZ=6q_W--Qg&+9E53$oq< z;pCfN({S_Xzk@bjop>DtIJ8>U1@16p+RQY##4o()5iFqJy@*EJC43(#EtGeSG10_+ zTMm>z=3q;J-S3^ut%huO_4Nj2fI5XW8r3ga~5em%d0GB1u+{@N`-X_B`gV{7ZUV=d(?~tCz75WzI$lAq9C4w<^$z z<-cLhc^x|elUUO9%YCf%?lJ2kD&|i~fmA$=@(_yxI zV%x}jAj%7=^knu3_@0OgV-PVoC7#O37C64A$|x2l*V~tXF`yQT^OokzvujCU+=2x# z$c}D7`?iJ^J{VioDnr5yp%onU<%g!E))IuK?mgRC5%SCFqzokFRLPv->5HN$k2kQL z^36QZ@v-yxDhLfautDrQ+6ETqMy_(XhBQM;)XZKf0*3*Kodonj5>(aY+Cz@N8tBxDG@_$=?fs|$zcBTzsPw$k^s2nCZEPY*(tm6VEf zVw#MyyniK3oA%j#Tur>lp*N>`+Kk)!H0qt52ltthrkQqIgm&nA#-JjFH~_WHSpCV5 z0uVX92NgoAt9>gD-x|)$A?hM=_zhoe05h z_X0JKA>Zdhh|k>W^bH|F2mJ4;R>m~%M@9b5MRbw!fGsNZwh<$U^V$U63=W)ej4i_< zEtj$%jl65m*-n47%EB1a(K0gNUPD>0O^r~Vz_72RNE69*P`WhUCc<0Xppor>YAc`A zI%vYxP{A0O1*Sd7=|(mknKG5Cl~Eo3jgL;wz(Kevssq=wFGI{8jh&o4uhNHS`oU2* zDGmuT$8_o(vK5XhNy>ppv|(5g2G-pV9rTugD;(n8NTpO94xu<6H8?Wo;g2cjdccdx zWclII#Gp7TH5|ZhKAY#z>pmnJ-5S?~Hu1PL9*ggdEp|aNnT}H3y#59pnWFt^ANYRV zK}xiuSK-)?7h*AnL;}p;b65fb51RNbrv{NfYwN-exG=&A(E}FOw3Mpu?2V$(E@{Hr zn)xFjExps~vn^85bTC<{oLe~Cpg!kMqqUkHcigA#M&IttdAr-<65I89wXeLEC;fI# zP*lLbUi`~O3m!0Ma&pH|sb={pTwvh!l( z2j9MW_|<8PFaEa8&;Lxf_0XU0p6_UM_nu#>9j~)JtE;Q`RGZVx`Wv?k2?~wb`JHXc zJOthAAeIQy6fJBqc{U)jHY4@k*QTQ`ZvB|)b*A@EX6>70MfK%Ke_uYHBHYxAbhz;SJ2E*l7s5*t0F>?jV>}`-Pe0ZDh`D(ms}FqY`Of&3Lo{nNzX4sr zr%pM#DX{&b?j}37SqPtsi`$H+lP*Rb4W2UO?4!Gb<}u|1;uCs)9hLPHrBEQ+o7nJ? zm~U%Xgbfd)Aa>CftZbAB!3e4ep1!yscid(z`sRr{A4Im%2bPT#vLgi=zkn4wsHHaD z?D!?FRX&Qo+P?0HF_`3pSBWdJSj$t_!yEP6-tAA9W)yC(Q5y|PWDMLmTr4#>!q?VkV178mChOe0lXo2GxVz(?Tl>$W+29A9#Rr&h}PBBJwyV>%=Me2Sz89|}CL5b%E^;Ow%27SoD3icL>!@6E6oHos zn>%WOlD2Kh*G@TJ^9XZ3wW&1bFW>8e?l)A50MX6;=7r)MSTt=+TW7q!_DSXDFz2ii z|7eS!#nlrzAKyll84Lh2PgdeQljfocXr_G%O&O|e6Ab6%t6R2u$is6kyDC##EDGCu zmWLipRzy2({$OAk1(cg6-S#X2D3q#TxOUmuEe&o^&pn3pRSQKqc^53>nU~`KnyLU? zo6&wI`4hk9U@Nsd?qi#eEWTKqkQXP-EzMV{?cR~`9i&2pKLOrJ);kSB3|qkW!D5Ld zNzNPC4W|BQB1%+%=8tA%Xw%vO^}ve8&7bdRZMv=7b9TIO<3|IZByqCJd@YRKJyd2o z@ne8kyO)A4SuJ}B)YHQ-Vho_fj!0hj9}}ODM8gEXReFpZ2LQwJbi{QmN#4`H0x3`d zsiC!lUld?)9M@dgS$di(DYYzU*rEjLCB#Q6UXAUD7%VuUqYia=W`Qh6DW{LA?-zu+ z*}-+VdYn>Zs%%N?G>^^;v;FoUajxQ&K{h!o2xny#p6C@A5g))_&_(LF2?^7xMr8yF z7-oVED@8wG2?89ohpLe<=dzc`;W_|o>1fXKgrvI*SAP0pv5uFNuF8BpOg9myI>QRU ziihaofV#k)AVHMYaG;nqwY}fL3Ac9MJ5F)vdXA}?nLUOT=x!LMFIL7A4dWwj!;`f* zWaSvc9buS+!Ki!0Mh=~P(D-7m_;dn>t!P?6OO|TfSRUcPOrl^_Z{+FeNj#?@_XN8X zv|KHpy;lY&wn%|yYvSfw-#RZ9DwPFQT*#O&A_#=HMH^;GOTPYr68Xf;&nuR zjScdIkTHQ-2;t<$72@9o7{uT;=ZtY*o~}HQB^o=h^L>?N3H7|tSSgklW#V%U zXf27n_2R~*z*^p;?WI1`SK zQf&fX!AwKt1}y%R%W10BKYw|d!Yl3ok_%oOq(Mu`T;;pfJ}j7n5iu%#j8Ev z85fh7m^iEbtk+Efz(yAwo*?Qc_dOopw+92MeLl;@ z>#%*X_$*2y!AJ$;2`z-(%R0&jl00n)C6f9NL=2Wh@>-H%c19vPB&t0(`j5}; zBMij@QnGQDr?y8O0vj`(XTsM0k_T5=yi`_DIe%!W$j}Ml;z2=AqHHcjA&AA`YIwz_i~9hH za0Cl^$Y2#9lV>7KspW`{D-J23Ivqm_-GMPl_AiQKWnYN{4N2eElz2D@7Yq>iuT45z z1Mg0EbjzyQk%G6R`8++GIuOHDqsSxss#LJ6e0ou{w7wkcjEu9E9pk&-euq=1QS|Y=R@&%@*rlf9P+%+~-x?j`Pns^%=2ep6( zl}@%kjiZaUI0y$TwM)?>-~RXqo8{PEaw)4YWH*JlH}}20{nrQh%QCt8*M`*GV@1r= zCq66`iU8K=+2t#@to~ycbXw4#F?!HMIZ!bAGJa8b;nkUccaSv}DIsW@_Drh}hffMS zaNLmcDzw1;SS=H{o}ge?ec7H>>;aXSnpcBgxCB>VzP5|ld8!(|QfTxTYf6kO0fFv7 zDAj`Auq0%bm6;ZG>0FCjb-is3VA%5U3x+R*HX*!r!IXHL}bs8|q&De`4BETqavonT$YgJ>6wE zz9VqY^~Kae1bs2M%?4f5v_u`uR99MAg$XP~y^q`u`;Acc#*NnzLGJ?j4DYb|n~)$V zw@!DfR#E&24R!(0F+8H!qj<&Ql8sC1<{lF3bBTJ}M5%B~`6Uzo~EYkS-t+5~h5T%o%u=c4&;c*$wkyWOq!8n)n_z);_ zz>XF{c@5pRPF4uQ#;bTLap+gF%N+d0R(L*5@hQ4CJG$>L*2~)s z_ok9k4In3o=kV*i9&rqKJ?C?`@nFKLFNcJ*!4A5P(|wul8%HUm#*qWeAg^%L%TX?( z(Pe_JkYrY~1gPma6l4(Qe|VLrXX~OKF}Ni=!Y5Rl6H}1Y!m$}YPpBu^>D{_$Eg!USBDZ{5#^QE`yZ z_JJiB;BALvmI^j<<{+dK6E#LzI|d_CojAgg${JUP3`O#i4~i0nKaC^kJiinmC5??i~N_ z#eQw65mA`6#UJm&W)7!>?$;i@84^KVgFqU05bRhe{IY@B2>my{)kZ z5nZ=xi=gUZbnq@P7*xq$Q9fSFt?KLEH(p!*%&uVcCdyCFsSXy`g~^HJ#kBweaQ37> z<~~WMj$w<$KlE>Ok62qfhC^%WLgrZ(k00Kwp4J1Kh~H0=fsj9~ z3yt};E(YibExXyLvLGQp=iV`#gqhQw_FP|v{wa0U)x?#5f6y$_tZdhCPAkqg(Nt=#>4D~}cQE0=YFXn=^* z+HV^(yTJyzlyhbcl3}xs+pG>K$CA_$`|KSq!LxINcwulG8QG)+a*A)Zb{<4*m}$aL zYKwe$#JMLB{c%lDWZ{%u_I3&P6af9}VJZ`UoOJ<5{j7y$awv0E;WXWB;Q#K+t5#Uy zPviSBsuZlmtZ-i&V(L#|=LIRt&B#e_>O8qh{8apyP^?n69;Gt89`>vg+tSJP5RSDR z;)4X)6y%Z*R%1^!N6q$~iUv&qlsz~ZMayVQaG@2J2Fp~v%>|gmz0VP#V>WgV%WZa; zxHO4Ypog9mRKU+DDg-)5`X$u9>d3=}bBIYM; z9TsptpkrZFmEnwHuM{WFC0N2MwlxqbK``}C?=LLXs>O^V7`^El=*buoB8hy3scd|L zGpuDcp|6hkwINJ#wDE1_&PfudTkf%Rk2UbU%YNz zynKh7vP6|6&l+6Gk&m8zDjG@A)_#-wTz?Nbch**h-o;$>BR*nRAwB2IVS`C+Opf2v zJ7;d`*e`PI%TNEoy?m^HVXa=twBZNSw~n5?W7xeCPjQDr{C#$I>e8wwoBjQd%o$?S zwS1=0OR7f7m||>X)%#6|?&I$1sqa=LUUMxjnrp14yuJW$rMh+XzRgYhaA@&C z5H_Zg?fscKN~Tu&jO+I!+tsD{}{ZowmJ64-eclY}{*@5E}K#XGM`8 zIHrAe&vLUzA^TnFIXhZwCjOwjW6QqAhU?i~FI10v(dqXQ$}bA|bufvM?3}Dj^bb0x}nXMj~7g zOPG(4l5P*ozH`q*z_~fGX8oat5#Or`dUN7 zH1<-MhlBwl!cBL%0*gof<$kFb;$|vDyNl9UOUWpg#wumn*d~Dc*viL;)o^mO^MRO; zMy~$K&@LHs;E9617646VkR`e;*G_QCkT!Rsc$*^^Bg=l&v;o9*neO&X zD;&XcDN;5LIE!#XQ3xol5bc|JvN>q>8_f9)eG25t2ZwcD!Bk!&cdOVIIREwEwS#{V zW|b5g=R`vglan9+2tI2#gQ9WpaeT`|hp5LzYGm}&Y1V0weiJj(+RFU$I21S3He<(> zvN=VM0oantP63UW=64DUMMjLW>+ScI_ip=@VcLCk3UWb0T7@c%ew;`m)}}8O3+Vvu zU*hYXaw%HGGl|l97+N+X3-k$YnQ3uyiJXpP35wO|@5qavrg_TTFwd`s!Cl=ezVGAPq!-DT$t>Hep8f&p_wXfOn9jBXdLkn6#xA4U{GFG2ktPy4s32 zG-b-6>-H&MT={H|1f(Ww~nic>s5b+DbiN+5T94;zJ@^HF(AC7dO*P@V#MyCIF*WQy!!8E zBF|5?f!Fi=;1>hU^W5lp`q;lIo+rS{VdwK2K}=QtQ_CogBP{j><=&+2)zgVw#83vJ zh9?RN(V3Y+QVOw@l$s!vj1DiKe<7ZzQ;sdF9v7*L1nV-0i@P^obs)!zqjsX5v%0Hi zN27C2>I(4%pq#wO97wY;S^NNFTO`Rbun8l$iHH~~rAA=HoC=OwS|ZgH)X+LO>2gzp zGek0Tq|{`?a4(yTKT5>W`N*PEITtp510Q)Gy~W;F60&O!ILoEr@_atkZK>`X$C&rh zWkMK%K}cuo-U6|VmVBnbsPfSv^s5{oJuQ=i43oH}_sBYqYn?o2iB~-3Z--WW!t|Z) z2LKaL7_X!N-H87DwtSqD2%c9EU;o+{=yc0cE zfiL;QSB5X2Ku#ykrNdgOdVZ*;FvuAi2Ca~!(OZU_mjmu^#Or)C|LxWD_j-{a@c4f1-D_cb3jUX z`O+LFZDY_@m&sugUWHNl77UsWkk~{cn0@#`(%3nL0~mr{ixxe~w`cqkmeB@a-f3uM zKCk!jI?VCbRNn%6hhk|uq%LtV4dp;;pSoZ$3OUwiRUO*#BMhTL^-!#}C`@c%HN zbu?}iZy#zYpD9~7&s!tP9^!v$*b`sT^X}3A!0wS4lAh<;%H!t0cpoTYZYhr|9)A3uYKEeCDs} zZp4>n+V{>2^?$SMrNMH%wM6b!Z=G6l>aYLa^m>9K^H`<-rd}=E5pOBuQ?YqZL0?r% zSl1-0Y|DA6`oji$)OGPXMn<}sg$D}L>h~PfIJc$tUn*L>)-J7Y` zt@tg`pkUmZ#$4MSbEhf0eUgp#zJ?u-hSz;)^t_4QAHM%@#{eeAG+yIWpD;oQPG3g) z7{Ngo`j9i)H;9LUaXoS_pgn>ym8C`VWgrvWF!rYSppYQ|-LO1NJMGp4F@oqP{9XH< ztAyraY#xm?)(hP`(jjX=S}tTdfNtSrqIS|3il(DE%$q-dt zc;DnhHfrVT;S=FJzwg>(U^3GdRZ_q2q-&+$4qes~rft#7_pdxMT6_k*VcpKokCSG; z60sb=k74XqB2EA-nMouGz1d6)^iZ$J+6ERWQEd%0xW7bW8_^(5ovz8g98B(L0y33tkEfGz_~7>_V*eL4hLJ?ilB?U;_xF z{Al!l?e<5_3Y}w@H!Z2h2g3!9r1!bL*Xweoeahta@zqtmM-JZd58B(c+s&Fr73T7kOZTE_% zhXCJ@IGUbwr}8JC|NW=``Kxm}>j0-pUmTMPEVKC5Bmtg-L*@VCgfMAE`j7yBO+Yj* z18H|IPq+A$_Hq9=PYw4q+{xxMz&1V8zfs=FQ|aK@^*qa0o}&N7FZKBldSF5R7bk%K zAs_%xuex3TCcHf}a#@Wz7J}l6@=An97WU17YUkpi?^h1?tCn`_h55JmNSw#LhaJ5K($#^c*P$!A;Y0Mq%O(!mT~M2x)?h zxoh{GA$n!EjtJS%MhzjQF1rGQ=$|n=_sLY~(XSPwKKXhwDIuR&dS(+qh_r#4L+IEE z8i7MWI5sXWk3>*vbFG2tmY|bLv?A=^Oj5@SiH-53fscn@OSC0K3hUFhfQQF^rDGp~ zzeg51_vW+%Q+j-)7I;N7z<>A2YD-HBE&#r1E?aq8u6b(;2zUKOQ5mV&QQ+bWMzY7;JoCV~{s=N}7(o(MmOq;V$>fqrY7$(5@jb9NS>y&SIXP;_Bf7bvV zO}r)CB$)vM;E^^zwP`->jaKDg6oxo|Je8NL-yp11{TC>g}5YP=nbGo65< z;XCsBk0`q*j}@i9ppX`jCRR`g9sPupxWb# zT{9?d_GoEG_0TOWQ&$1_+}sV1?&=j`dR>cL#^*1%hB5V z5EN>w@^Km)%O)q%yorEibw5GsCNi=A+eQ<2kaVtjsiVKI&D!=o#I~^SEoZ+S5Yx*R z3ctRt%0dA|3-Y7^>Ko1NwtRa5bNbig=0Ibeo@WRc`$^HnR&{l?dY0Wkq-(6c$wV41?BPEtE8_ z09Luw$Xx{zBapGmF3%t?p#vg^t5bjJzpy1jWJr5>FDZeN&H&e+yryu%jWhzqOSI=D zs*=!778hZN%};^0F&G3Hk+H2N? zZBO2|W|sH*FLs^vGD*)7!hg*5e#B;g6uO;0l&$AHT_|gaQppm?8vvW9_^f6h3dOBQ*_US&&FY>e%mYQFr^g|PECwO^|g$Yg)0%|pchDi3KA_SH3s>gljU03n) zveI;W>yqXH3Tvd_-F^DZ@1pi>e*^pd%(4|GIoGa8xz4Trk{BdXdAL))?{t>BJ-5` zH@7CNC>)#Ia!Sjuu)tZTvR%!^?N{rFxC`G2ub^&J;T5d%*AAtm%)x^6CzzQ{1?mT# z0<4|Vu9)X7!xT{CzIR`HUvDLeRz`A5&k|G8rE316YWA4W&=(@=e;ThdN7vQVY@!Tg z-;Zm@7aR1mzx(@Zk#7n9HwYTe^G`!nuNp%+KSQ0?;?~!~j~Rv?_+fU|$Fd?28#&1f z^)P!X?KykesIPHIkdHvj5Ae>3$Ke|At6&EJ;vqS(V9(vvRt!g4!lnti0l*NsJh~n! z++DyQLVajgg8s_(_Oq4wpbYN=DHUOV8WGPpyxpJ*d=$r7R9ZH)k9bTvfX(WhNp8$O zK@YGu;j{d0hcUM^)Odhy$ zVFB(`AmC6jNx&vmk^g0v5B$6wBe3^mZWypS0L&ks|F*D00`iXwLD1F{L2MI^DH)Q< z(!oG&ANl7%g)d=*0U8aALp5403ArT*S(!;XP$+>%#0;j{%yc+Z1PPGXOQO4Jz0d(U zh;SI_xLmY(RW%XbeT3+Nt*}q(h3w)w&1-V0!42pfRuv=yC@EzZ+CRK3o6lwcvS6_@ z93wxK!kIQ%^V!|Q@(s=z?>}EqgnG1NeAVF}W_A8<9*|&}K4+Ab#|l!hLE_(@)n3cv~FuGTVGhLs1a7pK~;jwn2+JH zkcpz5s1C4Pjzy9n=82PfR)Z$8SG4AG-<;!OV zDFw2TB2y#{3buQ37WoAW2HPYJ(+Nk`j7Xvf62!w-GF66agy3b{Lm4)z$3N4*fBAJL z9L~U>TMLaf$wo$kGo!%>zX|&MlzX3c$1@$L93ESI1!<&1&_dDT=hp)AV?kBqAQlob za$zD15n+HIk@KWKIT;OInrdp+TG+U2LM4c$5F&ZMtt_=n?K@vi?xh?8{g-8OEXW!> z5~%K|_V!hPHfljKa8WNQ3m`0B4|VG2$gZB=zeS5jAOkA|DJTSU;Y^R?G(txwj^+;l zg-#)>nhZ#ag8sO7K|c39Q4!(j_5L_sPnW-pW^+o3+1e;DS04bPy9?D5>y&!WOzO{I zr@IR#LQoO``^@P!2q5(RrGP5Mlqd&5Wz_oV^2mdT-H52sK?|Z;_XO5cRJVK_7AHo(>t3RFsm5r?w!~C*+eI;hXCg1d=~PRRS(}_3#Lyb_>F2{$7hUUcd(h zk%I}SX4cfJ4}rP>mCXYmSCagKwGYh)9qL)}g-{h4#GtPl8SWh*1sk zki>?~{p$|d<0m?gKnt3U$8XSlca=7c%le_OP(YR=Yd@orv%}yER7Ozgpdb@Y0-#cF z9>qDzM}`~%plc+=bf5w7#Xr|Uw`IS&bqkQPmryn=?LyjJKp`OP3teCi815kU`!lIe zSZeub4~GIOY|f{qQ8|2P;3Fvd2A@*2btaJG3sM&F)zZ2Q>TaIQM442wQUB>n{!qUo z<%!Z{OpP+OHigqa7_~7C2B>iav~zH&>xm2R;`Cgc0=>52%}zp8h#d2kp>a+oIB{UH zuLX@S?hvJ2oC6O=>Xu=cl;*ry{vHv>CP#_c5_f>@eftarhu zP|Uy$WJPyCihZziw&YXL0gX%?BdTImG+5R}V=z3sOvXa<%VEikppQ8-ik5JJI)owwN z5a6D@QcQK+X|Ja4N3K!_3}eCo5~j2|`Oh?Otw9!DAVzzX)7cnyxIy`CqR4`E4NMxb z!{9FZ=*u7#speL7i9AU@V(K`ZLw|-P8V49;bfZyW4r|b#XZt=_Rz0qN;xXmOzXwG* zfh-ca$gwBdArQoJKnG(+i;3|YR2Nv9lrazSk(L-G_q~NH>pz*1@qn>DL_PfNQUzBH z`4SOlYKRd-a1O=>OYqD`s8IuN00e?RMkJLW^>`f&2Yvgl1`tJ64&)?-0A`=<4u-jbjl z=88|owt?H6aIyj^74cD&7kMyMV?6yOxa#d7MlS&7m_Z>}#Wrj{WM&F&j#K7X60)Sz zv57?-cOV!Za^?!LZnWG+3bZeTyXL<*!_g&};Tk{Y5On1#AnnMpC3XuN&?X00Y z{}y8W)*gc)FV%D{)~l_3=~eBd{a$k{XVw0*VE*Zxd%fDV6<%8n6)oF))9RYwdB4vi z6G3@CH)o*Yt6?Qz+mnH9aj%>6%k7A51$ecg5b)K8>2m=+MBOm`OFvY#&*K{9$yBQfi_$EaXV_=IhZy!tJI=< zc(ohL+`y|NyLX))j6djbnx;FU_#1e^9G&*$J{qeOT2M(q71c8f;LAAe-bW+8Fdj*= z5rt+2ZVvTFPnWmSWd_VBBe}ot*Wj_Zc9*rXmSovsnuEnR{qgS`mxdqcHS=qYR+M`A z{#sd?+f!H=`&)la8W0$Kb*RP{sEy;d+yV%cavJsLvcbAym$l?E{l6_jht*wXUcz?0 zX%4R$ah`C`m-AmNBjOCh{?g!XDq6 zkNLiNPY(qQu)StpJGs%RQAN!?v%UMDSEhS@I8=CW$tOl=|FY;fl82IAdb8XtZt9>4 zn%%M7I08-1AVihZbx)zfkq_J1wphj56JLBr6Bt%^GSZ%$B9eD}a0A-s@x9iVnx{E) zs#{XY_FQDEZuhMR-$3faa>s+Rpv+8<4)3NlmM2o)91!W!*s@@9*|$AieXFTQl+?Z- zXLdqe6KAi`Qk)p^dfji7_@a;;HF)GWv?fqYvnN&o=~`^UsY%r<~GzPuHbwDn0Y(D``L~RI;nD zXPA*QXDLoO3nr_4$x*=hfETrnB# zGgkkE!_Q=_|6C^nImXnc!fQyd7@0N;S$S-y4^G-;;7c$2aBAq`*Sq}fGuthT{cE)c zWWmbHhKAWrC`BxA%E7C%YEtVPr{RLsW@6aTx^iUD)cB^|kYMQec%%@ihU+imiKkKshy0@ogZL6-;&6?LU=PB-f+8;$OF1`Pf z>IpgyAkke{@|Z6up?yEjeDj_PBr+%p51ZNV$Vdjkjyx*zw%~Y#NXUMDgWjl!V%o6c zU8UDO^{MQ#CbC+k@U@8F!O0dJS8)ujIk)Yxc2g(}A~`)ei_r#n(V@^* zFds{4|3V%aB&P+cPlh&dp>Lx*rC#Lznb#iD=g+|Z@1bY{8A%D{i?w)NsLIbPQU*hOj0 z2IF|v^t3h=}L)&rE}NYkvK=us}fGuHRk+`axYi13sdV+D@fCm%VY z9I^s)&|M+yg#98A*vx!t==wrO$h5?*a5P^}Ni@WI`Fa$sqjk4$qn{sTS}9J*yae-! z^m2E~f>BzdykpPCm7PuHx;Tk72u z0Yo~T9!+x`4sKpF)?jFqDRlNQCzyp+$TIm{*W2^i8HSc$ORVA?HZ3~}N9vGS&n`fB z*J=Lr?J8JW37tM@H8DlRX>_Wf3tl3@6#9fM&*qZJb-gGn8k(tC1l`FzbdLHabtonq zE#5UIy|fEb8ha?wi)74jm|uU4k_7ECQn*CY)rU7o;T)K)9)R{s>A)^Yu6qmr!3q|Wl_(cSX z>)isTa9xkwDf8lofz?ORN`b(Io9pq^Z1eGC1 zRjDYtt5J01A#SRnf5oBZ)8EyuKJsnysk~g3ySD`w&8Y#!X;SeD>~MT9!+}US-nc@a ziMdq3KrVP^-8!8qhcPCph6w*na@5 znuvbE1C(uPKc|^^w>O3nSG3(q*m05=g9XvUib;(?m}f_Z+y@YQCqi=0rpcXyQ_ZI> zw(Kt!0J?u1pd1>da4|@kvbv#`;Z|i_+Y7m&k~<++-@@=~E*lWUoKbe(P|nN-^@z>A z+bHk#D0YNF?rH@}c@l;BQV|h@q;SzK#3nXBqPaH^=4XVqm4p_sS@8xmC9A8gw+#>< zjpyRBJ)qlC+&PY^0hO)rI}zjsPuwI zH-!_djB&dIKgD`secw6m#ig5&5|-?ub_>Sjy}t(R;m9&CCvS55>vNUIh?7W?D{0=l0SKYhH6d1&n-pStf_!sKsBO)h}IBSV=bK=@0+<7eG&JFKJG9C4v1oaQ1 zcOTI5U(@)HR<|tO$hBh)2lzP`CVIq8q^4RkMb9XpOBF11xa@)6qnQLX-yINhsP_@y zXp`aLSW=0a&ZhoAKc!_&%F&B0!@w}2G&=OIH}UH%Tb`c0Vr8;N>R7Tdb$iRBO&woE z%xazY4@pFs)w{-|t`r+BEoG;}!;_R-s-oAeFL|@#%S{dKGDB3$2lRaQC!DseN1a(Ky zXGd0X4m4&Sw2Ezwm?#*@Jy1_xRAY3fmz>@;HFm&Zq2+u|Gve0Q+i&0hW}S8JC~rk~ zpGM@XhM7JKCUy6OVFk^wc`(d=?vzn&eaj1aFC5*SS$lBD+0tus`?rko_64{Cy4?14 z9;{dl^IPlDumOPgz_~^DpJ1#>c|nTlA}2JvH=2rtZ4EaS+s4*i#CLAkd?61>ol~Xn z2>YbBwS0@eas6;a+ldohw>I6ZI%y4Qq>EXF8!JAr8E@Ru7crpXGHhqlq2d$7{eVO! zlLI|sc9t7;vW)7V>fvam_=P}ikdz;ObHnhZhc;fS4Ybg_KC|e=9BO<%`Bk0Ke>DyT+` z3|VMBQkIY>LvB;l>e6u38)C$zon~EsJ*DyjmeBiK( ztT{PKeU)|Pi<0x%U9DQaZuibNI+mPVQIc*3(cg>zs`zXmct}F0L?pf#`BT=_v#Eu1 z24LXE7gWOrt^&i#|2Tq6_^ju8Hv>5?B&(lG0*TA(xoVR5B%W`y{xLtE$A8OEaq?dr zc-}(GggrctnBlky_sxX%bfiO?Kou^9O6GMOO)TL?>-lx(=(0!s|52evbxsxEdh#_# z%gPrw&Y(&y2AJI=bXqV_(b;<2%YiI8F`@j>@T4IFZik6U>Er1Hc$^DNJC^*r%XE5z(;3z?65# ziZ^Nev98tWePOq&(5;ENmYK&i#~gZ256@6okmbkD@;UY|rPI>Bvcw?=MBJ<>)H#`NaPLupmPx*JY8j3}$fHvUe-@wzRW3)2J&5)wLsQr zDcCnb5rV!KvBk89F+Ev~^mij_Tb5OVw>ru}=lSh-z3IjA_(z&4QYg9qUN|dLX;MC2!N>@OsO$CI?*VRCxy}R^_d`JO7eT|78(C zUO!zbBbA*0$8jebUlR**9^Ftl+ zIgb3{s-uwcq8XIfc^U&rY3!Cwn6{d>k4HM6j#gSK)DSr7<2;V-P9oKybbAG}m?HGk zoP4Mf!MV{YYwK`AM|-;33t>iZkmh-;dBimLZjy2(ez8FaCw#b&a7x_Vd!&{KhK6zi zD)C+O;0wEZj#F_AY5o*>7ZPy5VVwVudF-;TBjU~3PFP_EGom$>3zIi*L!Xx&uef`d zSK6WAjLG><{5g$r&}x7DKwJFZ|1AUiOXAf%S7Z?b>idPc@Vz7Eww(Y_mAYN>dEz+#>KoaC``Igy{$|2pM^8i^o9flD1}2K4_Rbw z>uI0i72j4H@b^J>^;{E!wMl*9W;cs_#~A-j-~@012YemZWQR|w{_CgIwN@X!{7ZdS zR^fix>EHhP)mtt$$&WNGIgln>xi2VqvWxcQXq}(O_$=-zcqKxglaYZx0T!n!_xE)U zQT@pEIa|B`nf#eg^*_BzjUe>N>v`U4; zY8v6T@Z@p}Y(XyYn@TO4+9Z!(;fapn#(dZ&QzMDa;S^@(c{z5ufT-r1xbVJ5+yuR+ zzYGV~E~NFgCahOE)T9aRBLmLs~uwZuB6Dx>zQgP$?{`Uwp@_V{g<{y%p zEU%&nv}sk46rhq^h((4n^8s?Fu+F(3e}LbK8MPq#Z(fj%B-%}y3u~a07hB;wgmH?e z_9*c~i!~z><0tS}Q*h2>a2$=92rnIM2a|^Izb2fe_W=c!uZQZWdbq=%JJYw7pQ$_S z(5kUyr^k&Yg6s+#-mu@usJn&#J|G3qP2lPfo=ln+2JAtz?tCouzF?6Sw41*c&D#@df((S zVXGW^kr~Fy1$>^FB&pEMc%q<4$^5)^o>cs;kq3yty!gA&nB>e z3!~AsjRucK-9mg6lDhA4DM!$=pCqF+m{dTz0zbx{DPRs4Roab!ViUY6GgSXB-oR`6 z_f3qLXaTjWAuAFLz5bkohI>KKTQ0>Ic_T)w;_qE}&m-whvgvFQbs6!s~(!cU=j zN-R4nwp{Q?2<$lqrGLOH7SsFzGiMsZ(_RIIo*G{JRyriXUST~be3*KO8E9w6C*nmL z8$sK31J-;Yo^lA@lwkh@sLyqpnqsT774Hr@fF1A?HuusAhUEuCqhh=rJ{9jp21H;o z!S7mk0$!}U8zO$Reh(%%hjGa>F`nr$DGc9yaWsdDOXs~M{>a9odCwfq0j>YPaps@X zs|Q}m^KO;0O8&@y*nG~sJR4|7h>-VY~saaK0xwrwVxno}(dRlWcs z^Z>nahFKOsxbhM? z&uQcw0ESoYR)+7n5Zozbn7Ki+beuNg zmJb163ZF{yZ7^z!y9d+)9CYRWGNsjm@S_b#1~?A{vGi<#PdIccOuOUp6@|a;jkL3e817cRZmmIJ> z(8y>C7vilvVIyCd(sB9^q6;VjO*t?L?;7yC*~m8(Yv8jwXey7)$Lt5&tXJgtz=t@!ItdP6V& z`q`~byOCT`<0lUYwl#nLp+na(@}p`evQH2}=WeoW9<@~bKFKg+>=feWXqcmQ6*CTjUkPQ?F+-h@Elyi+QAQ=No7sEcauMVVdbu~SS{{2js;n#-@hB- z{_D`mQW-F{3EXwW8t`w4e7^$pqVb=9{g7a9&S5D!rFkT9gRL?}#%%>R0oqh@+}?uB zjaueWx0O}kh5&4ui8A!zAiEPpp8$b+h4kgPb*46*+#!k}STBLSk#<6~{FPjllq|g; zB&BIQ@KvNJDoO>SG}p}eNr%s6DW8L6b%aXxHAjFa)V;IN*j68_e2^RNIE~SUqx)>c zxrojVTv~~TvUwBvulUXP0WkM<47~aPIe#vIAbZfEr#mgEJj>+kLD{DaUe;Af z@O%&a-3qjmO;Wm)8w@W2@Q9hjisO1R3Ve7~&~u^_;MCSO@#pkkG(pe@Z!nL`6V(~J z@5C1N3*yq-bUhdxS}3Vf-Kz5}Vd40v%Tz&FPKC+$?Cfs0EKo~pQ9)cr>^^-_hw%!H zT#_1ty_Jt5@`Y&~s$udfGNQ20Nw_gWv!6*L);IPshBPx_tidZQY5!KdD&y!pgvjajKCkxn znH^uf2PKgMxjq0}@SDxfW-z0YFR<;Ub}-q5DoD@D{A1mGbh%$slCI)q1_rwNUb9(@YM^6w?S5v`O^z2n8Bda60lJUp{ zz%Lm!xZ)`5tY4g!Wo>;cS|Z>=r8${8sEN@{-*>v#eG}j)k&~>gdNLK0r>U}LaQc&4 zG0B<->)j*1|9%Y|hy#SBG|_E!cABF}o05}eCp=)x{fajbvGmv4cd1wGD7vq=z^;rh zYD#NS=vTB)rQ$VB)ZsjQd~=UrTS0`r=X|B+T(KCcoe!SHJcA+<2LZgLkp<_Y&CO|= zrl_F&o}S3?h?z(BJihzv%-I1_{2gQ!eteuJAc zt{l`d42yfL8dl~J^P2cBb#*|Z4N|m=`vbNctD3LI;!I@Wj%3q3t2r21!oAtG;;%?T z1MvEI@?K54+u@`Pf64FHHC*?^%3O+`PQ9}-)Y?v1DR<2wFQ;S|=eQTshm#r32m^xf zey|KE8v42`{fH+lmM7meV+io|sM;9sil_Z}K1Uk_>?MJ>&3 zp81HHNU`|<$QA*;oiUhcmA$GytMyAEGfhKVllD(+ z%&@-%Gxpd&T^F7pj)kO_a{Ex1eITLL{Iu}RT;%S1z{1V^>rbLA`t!3UNKE{5e7AD^ zAenMb?w`@Pka__;8o^6KezQHH9s95vB!JVrT4#BO2j zuZ1apWY`M4Q{w4xY%&VfzHak$0SiO2D!f(2pR*mY5$yu!3Sr?n~%VEq@di>J!aU*3aYK=Lbs!6oe`oTQgY)?PQb zNMFDreL|iYwFB^X#N<k~k;!nSWqY1%G`o5AcY=cLLIx zl26)|4-s6<{#thsfx2?G^0jUJgQ+awyv57JgIW@Lrg?DeB(e1!6foAmTk^5TwAr(b zz)@+z%pHP{zr=~L01)e+Jn-h99-F&Cp&Km~Rs;sQxyvLaN5c&H`YK2-k9aD37M zh%w}hQ!f5ulR!PC__Dlp+FD=I z_~^>E(si)3)n{po6HMSIM(Fy^K?YY|ro3TfxJ&PQ0PbI%>8c znnv?z4PU~R4&z9I%M|}?rMgyJ4LCt85M-Ya#7kB!1xyL?Wz_;sm{sUnztg$r*1m|f!z@RC9V6BJ3yKg+2cHzRz(~Bo7w~9{Dk&!GFze# zx-I2B{OR!*=vuX}puCVfq%!>NK#orst6rgWvM>7-3p@j6A>0Je0n#8YhbmFDyq)+D zWemJyvlMF+2z2toJJ6T;iK1IrSTirAh_6i|Tl$tZJhwmM*-0uWsVQ3cPbzx$)dW9N zeFp@a$dT^KB^H{jf@~Dz9?DL1gXeOEpmKs(kPZv&2@3Uq4|SP+F&fN9{j414l)V&O zb=^}|M1DOoPNA3qA50BWLPmmB9h_|i(mi4qW=kl&aZ*H?+#LM&068F1C4ef?Gne|) zbI_x>1&ssdGOO^O$Yw&4?#uGt4hu1DaP3?@kNx;$1?+0Cn3`J!Gj(LDPd^OSNoXok zCB)fBnTtif5Yu|~!4Z0avpq**ESWddTRVPvFi8qDxz|~C^4z)kB#AIb>MT>Hd1aT1 zA3Rg#V$P>DH+gQEq&!yt^9ug|kX$P!VO2QnWNX&k(C%w(tndjd%mPJ43)<8qCzEEy z=ALk_A3l8e&3-8$Zl{rp0jY>QbQ4Vgve*(@yUS^K%vm#5P8di!wcsNbP}HjjR)Hk( zLCb^Y#HDRz;(j!q0#x7@lHf40t+>-wi2gs7M$_*f{()rrNiw!isiYr`Zu3Js zL_OU3-y50FbfiX@a&klqa1R8l9e#J2wwwwmT-(YTctwuQU;xy zeZ||d2FB;{kCjL^<56E3)`QhQc|B4y%%t!B_D!^4H(v-p_>PATdTSfJA#2K~QcVdF z&_qM|DRJ?M({6|q56M;Ip<$(WSr(!5{P??YFY>B<$TJZP;@L$Q-}~qrr$eidUFdBI z6Z;0O-_?k*E+8=Ms^eJP(Mo4$|EE*Wg2rw(uM+JPTGZ2#I;*r+cYY=QwMyYXalf^g z^YjYc5w#bGym|cbZSjpgk_`>3`vq@uuN1_zm@R&gW65LChJs2`Yx?lQcAzsaszHIi z->&nSM}gsNfWbA=P^Ax0iKMka?$tx+`doslO}v#4da9I^vC{M>xM98Hu%|Q3|AOHv@iS2bI)6q$>MVlsMdg7VGb=^EPc{4a=1*0M? z#N&p3%9Ls3P7L(roIwkd0n{+mY$uSK3~CU8`<9?i6=ZqQaaeON9jQhX)ES6@Ca8nv znUu0^ZqM%188YUd`jp!B9#B5T12fJ(cP(}w5ua2LZ2^Im7F(mW+}4wQ57%*>LcPB2 zvnF)yL{OlbuD+eP}o0xW*XkIDB6PVwl(yRHz`Cq@{yo8 z@u~Qxk{>rOnk@d0|3}kDC9dqhCm3IhNzaAm2b&(x&B(<)@f?}=-zj^1+E;eA@Wc4L z7T+?()gS(L_`P>u`fSAnjpG}>)%oJPiq#{&Q&~5A>cvkl$bUY3QU1jj2ds{qI31H6 z6Z7^*)26N7Gv(WnpS|?*`k4BjJC}Ra3M6X3o!@w8{_b_VjyKPDt_bx#>tI)NAPv3Z zheHRqf}PPbeJr0ZMXsYn|Hg-KQ-cF4lTGWdX{$%LPW&@+A};nQC9|I ztfWTJ7)p)n0ZGcDZr{v&m9qhE?SMuK-hW}WoP9HS1RdcZp_Vl4!MK_h@3oI}#uxJV zMq9#cKT95{n5k2PSH=1G9!ermNc|)7jChj^v*(NZcJ-uJwP!8>C#G4E-K$>Nb-TR?gLeWC_qg{Wb*w=`WPoOb zJ34y%PEl5rq8`^lpJ8hyAJCR%T`&aSU)k5&E|cqZ-yOcjF)Daz7l)S)Olp}-qn?X> z!Dg)ZoikKl3fmVT%)gtdaZ-NoftihOE4dh~`q*D6Y9qog!msmxwfE*>QJ&ehxShn9 z#5Q&(3K0hqQ%(#L2SgA-lMpS411PhA5m6b9ARsDWVocHyKnno}P(%@x87UN*%m63| zDx-ozB9n+fq(o$_yVeJ!`}Do{+ynHIRkBYl%b%et8|IujGK) zYTA4&FL{ZFLLmhzaMnm)>%k(Ju2m~b|G>gLQ}MA%KNH4AtBUTg7kLNZP-_{cE$fUg ziOj!?RFw+8h>OXMV}fNPcpDqH16QtXTL$b z5cwO;4V@FEW>qkm7v1FkgU&6kh{^37ZVmRUTF;*pS~0|kXt;A9PHSB#SI!}fm^K5D z+J@Svl*VtGRTdwjD^o&EASPbK0XmzayXQ-EdUIlBG0wlEbHY+2wwuGBPP(Lt&sL6p z_Y}1wITpGdk2j7dph;s0UHO2Fo*aH5g( zH4sv4)zR36?L!eeO}HX`@=VG026hw++5}r#6O4sOEPaN9qIs2mpD_%RF4F{!ST}r%>EqRi1Sh#(|CJTT9joq*9u*Uufd=^i zI}_Lpq^*wE$UclicpL0f#1gc{=K=L_IqdP!Up{yr|HixE4D6wPacn8=c4g2|SGMg3 z`Lq}+mNrdAr#76DCC00!hoPj;PBjvv+peJ946W!*Ai^8_cJ4pLFitw78z`o1q7WXw zjFVF1qRF4amyRwCi(s{q0|m$=!Vxxo{&sr1pY%=N>yJ)=b8~P=Yd7qzEA~DGUSm7! z`#qrVY3v>VmyRsWg!J43CFgUYbIPD`-6KVK^2arFF2WxSa^{qn%GjZM$JaA(}&bAZy8({j7D4bx0Y@F$y= zhJm4;yyWhCJjBxUica^#(^uO)unua7vM!JrP9ZN%uQ^5@Et53b zR*GoDcZ$7Tp3twFPDpva*zxHI=yc%QMHNKpmUIL6SRyJpsdXEEotf8D zHL`#In8&zCj8oI;&)Ow2`f&45m23sln2DoaDa9$=lT;NX5;AfN$&;cbDy{V)cVgjy zNvCJsq@d3&E@mkrnC>EkJUyo!J0diL>Jq z(hQft+!^Q_Y~#y}Xk1;@cJHW`_eddhplLN$bofi{`MED2(eg~(#rMGRZZ)%o+-Ref z!J^jsJo!7?6OsIn>sRoqj>Vq2lFlq8&bc~i{4{H9Uz*nyWsgWVV0=+z}*K8&^5JeQBpo)~*0vgIv*qrny^Uu+1a`;?qrY^n1m0 z;aBNDzwR`V$d<$_bwxxP*?~t&JfZ8vT;FdJkK(CTU}Eua(4iKYJw6J z{b*fXe7>U`hFq6SImIC%tK%9O8F1$Z2>`VKcQ*Ogpy1u+Fqu0I3I<2`0(`mcjA}`U zf|>Q1NnO)17snPd81BrKcJDTQGFq!%2*b34Mm)vasANcIrggv3(@k_(M%7_rrH>oM zYG;<0-VPVpxPCLZYx2R4pDstB3=fBeTF#AzodEM=n>fiHW{CGuk+vDnD?a{(msE2xTKXO&+@N}_h$n!0czbfo`N-Mh9!y8zG5l+8wHkj~b ze2^~!Ib{cxnSMp5N6YH%uMiKh{q4?n#heo{Yj19Er<-8+7^t3M5z}(8;5nre@Ggk$ zHdvnF_UM9=oBx3bcMPJ^ML91hWZZ`H(! z%6j|1vR{XboMdoB)QeH5>30slTdXx2M>ii?yY#nj_9GUOHbrpCX ziOwli5{-_Eba0%jw>yFI3wfF+26mv+C5lR#A{D_P%M7*RgtvHQk5Y%{5ovP=p z&>J<9^>Z;$=SmQA1LQoZn>04nC?2r1}u$bP#LH3lTly#v6pM7j{ zTEZwuWcCJ~1SwoK*|TNWASPJd;Lzo<9n|7cTgnMqRAUIrHjKtK?bx(Fl$P}%f^9lj z@y3MT5UOsml_nA+ZQ*>|-G4jnF;rMZ8ao^_Y#+6JQy}PdSaavT@nSl>V*(@TO03cb zi-&}a#Mx*_%vG;`F@ODfMK;vuy|)kE^-aiAY<{(GX3k-9#}UU;iD1or6=5hNJsW(W zSzdjbL8~CjSJ_jj3vl&mCCNfvFh@)+m`3 z$CvN(q7erd3L7R8P6N6MH#%TVv`XL7De+f~(@=g@32n8=!{lNJ%snU# zTmatm{y5R^FyMsuLCSXnj{|L$&PLcD+p~Eo62le518RqJsEjLv&EZYal(dc-7VVbE zP!3TnI@`ft0y;r+kDCaGf=e=WrpGvxYP4(Q9=M!vbBZ!c98_)~T7QXbxv zG(E-Iz7Gqa@&wCnyDS73HuMA)4^Rt^Y!GM(#Gl45_oAaKnw=`1ZtzR-O!})rHB=Q5 zT*2Zp43~@Ty73FReDe&9MPkmF1Q4_HZ>Qf?Dh3PwGg3J=OZzxNR-Y!{YffQ!QfS4NL?z039nC?`m!VOCxj)c0Kir#<^yNVU}v;^TN^wIK7sEpFyEbc?B@Bj;Bj+o~dr_UKR zWy|UxHp4KJPB-C2K4%u(QhsXf5N7H${b?J7N<2i71k4e03`h_{!Wj=dRC1Hb`eF(` z@l8xCQt3s{E(*Y+8dEAFTuOzm(n>Q5;}l-*$P&fQIKAEWMk~ZELb}^jx>ectxiJAv z9@Jt9m1^t~(ivC1Sv4l|?_az|_pKdb+7}tEj05~OIy)oCUJEI7tfi9L9b_FF_a5B> zE+kp=V}m&tFsz#};q&jaI18IRoWxRwL}7xVXGo&GPfyckJDV93G9gf7>@<>)!G9ku zZ$%wHP*Y-?(EmRVqP%B=|lvoUaqZOZ#gWMt$L5_Y*&Daxbj3?2VW^7dWqS2(v$O4Je zrG#~0*c?&$s_!^LU)ZuA)u2hvHnB=7wy}4+`rEf~!1e z&GnFenY^Vl0Sq~;iqek@^DSwZXw>-$=9-c`h3L~ORv%Flyw-s-AQ3m%_8LVpH0XvPdTyJYXs>m&Vro% zY!N(&A5M+?%UpYq7kl9vX=&xHH7DE~ITX(o7bFaxC;XuoI~q>tkZom)$1mADFUY)e znDBAL$7>6xd3?^K)O0+im?FQ9nsaJT%%s=-o5IY=Q8u5dB|qCLX(XR|$%8L{K2(*0 z;Bp>9o!D;JYgTqUtS=aCGK5Fd`a83@fJ{^1d1wEs`gIMb-2}v3sI11CHFxOX|JG>v zU)}Nrl8{_ZB^5(Q6H2qp2?GJq`zh6@qS7%(@5fvk5VMRER8He_&X$wwTK4d)8HSc| z>O_*Pce63&kaW-&=k952g}OZL&NFEIx3Mo>)`w1jx$CXrXoIg%G_!37mfDhgn-8jZ zrVjjNq*~aqnZ83n9^FwH2=py#oEih|=F2<0^PC%>MC+FLd!Gb7tR7L(B z&jZWs&Gis{6qj}3=$DWZf;e+&VKZrspwn1d$#K9Q!^xFs+j8)bu#J(#ip0E`R=&ZV zAd>US-z4urhMoiC4e|qG{;)rcQHpA)L^{!~-2(fk@l~kEf z%SGa1ti@!tsc=KLd+W8?Gptwn+{1n5;iyu}7p?rzI=P&$V7>8V4}E(9u5WsrN)Ldt zE(lMwb{va$ouSj(S8!jgeE9uilEgsQ_4daYQz9M$DGo! zfeYnPp}G6#oFCik{a>wEw9(pG0Q*hc*O&`8@!tV|O?<`13dVLVTR18B1E7yLcZi;&! zfN?OFZdHysKxey?7KmBvQ2|V!c5>5EYnb_{)#*|ZjLg!iKG0$dY_4pqhLF8FX>J#C z+GS8ZxZc8DFndcU2eB$pw;6Q`kKvO5h#XZ>J@H_hcfG%pI+EFQZ!{AjPTbq`camy0xtzwq&(6Fva>VINi~NH--k zU?NQ;N2MW4i{vt1{Iee*udt4*UgV2Zbh5lU^wdfmcdVcRBJC%Z978wnurTeki6C3C#oK5t_kK zkU`<*XP5PnHfm~S8LOx9-NMRqsU!4V3tK7t!;lYp;hANKZsPp(fYUw4z_~USS=3R2 zXkc8?8;`QFTgMI&vduS=S$WX|@NrUYHOJh!O3U_jP?M+Q>Zt> zMOx1Q#h6(YmFhUf>h3`)vK0UwB?$2@Q&BlP_VQc+x!KsYjyQl&v}OPST8wPPh1x99 zO}rM0Zd9DnkDzQex>P#ynthlPm*IShRff;EatP9|0yM6 z&4rkGu7RXc7gP#T`FER-zdzC>PML&tIAgd_;U#=%{^|I6dNBPh+~7)D<##C(gP;SY3i7Bu^qj zoa1;)lZI3r^(R4pnmGqvAv4=3pehGS!LhT{wZf=6l5M4$tMpl4=ij%1Rf)ljP=}CJ z)#^Xt!EmssVC7G(6xzb@1le-ridOuE3l@*fLCX!25o7^?$9PGcqiT~!2?@B!4R5`( z$5RNa*UOM&#UBlyvo`YH_hBj7?>)PF^r>%mxn)n}cV&b4w;^s~(}R8AYI<3w+`78& z#<=@~0dcaj{?J8ns{=$~mz$ft*Zpv#BX7#Rdrdzsj0YW=rMkeALrGJ@a}^MqCeXBX zu-z-e(|d~6Qs$!oRql1GLy^sbdY!;9&T$|1@z5>Q=FbC_kV%t$)j8CQmHS)T>V*&P zI{p=VDoQv6s~=>z)>0jFaaVoz_En|H;CKVZapsr~#ss9N{;XASRLtkV0KN{OGT9dV z0{#=^48CM+6C^fEw#-tS*fFDPWv>6K$}Hzw$4@@!8tCuM(<*4mYw5VVI_$@^FMe6I z;nJYiMhE#fj+sy7-?rRZzIT59rRZG(A~nq$`S)#P{Pn*U+Tf^kht+>G$bzoF+iw_p zElS}FWEG{9XN;cgRC^-UM4Ne-Gh0RlE7hz1$NoAG;xIF-d18S&h>;H~$74oUfUQxgibckQ-m;j}7i2Bl>n87hW$vOG zHdo%js7ir0w;$spj4+^DK%m@DA%kf~T=11NqQ3oP9$JDddkBmO!zYNzemc_lt*FE5 z13TWb?HL>P#2~|~a4o;iGFDLizUH;XZiAiuMrYo!@CxeR-8ZxF*0JT@^ESS+!sZQ~ zu)iUI{m0(>kWa|Fs#*ru(HOe+ts}2#j6-rdn3jj{3bv?5&o_+IEL)$cTCVyx$*yPigOxP2fqXhbr@b(x0wf;>DfxnDW&128YJHt+@5y9M?owR%8CUccD> zzij63P(!xiO(BUxEr_i~#hfLRfl2R9vKerCO~`CMVVV(y0kw?UMDIE849>De(D;O7 z$2KExD5g^dI~wS8#A>s$bAzAjp{hO9`-4O>emaH3=Ka1u^H=~~W6v3xg_MDdbp1m_>ESY~L z>Q3Adm4ny*y4&h&J>YiiSes_}$}zHhg4Tn>C}`q}Ui8$rllk^zKj$(NqF|65B8HdczW! zQ_>5;4;9(=FdS_=f-vPv9He=RA=pmpsyQ~@LT`FnrB!eFhRS`L{Z_eUm6}|7J$y}@ zx2<+wZ*z0+&~#1zx}^GxuB8vowwNTn75+oL?h7yJ;qrk~`+7njEZnm(xF)yOB;lx>k4Rohs zO$I1?Hq2#jp^p-U!+#N<9SbX+={td+XM%+y!{Cob6qqI8I-rDzB#jtIz%8mjpn=aT zEj;vXT1wXZN>fyu?i`zmUvNpWpDhN9x*TZ9>&s^spFX%ZIjW?pO*>>wx){O-3~M+CBKp4I$X=unyg9DYN%Hv1rRL!sOry z?tZ*g(Ow-mnP~4kVaTJHSP2g1kGUCDYV|>@pvNz>4nj@X3UfJH)5-G|sZv09-0`-d zPW?m?ISq(ZZ(`(B$PhJ$FGSt@$>TN~r2DmmLCk08{B+@%xgP2#HE*^oYEj$$d6rqR z@BVL^esQnANa)r+V!Z*4aREay+ht#nJQEQbF&d{r2aq#YEZtf6!tJ3_^YLQ<9@pZq zoZ4u??{QtS$fv9E^@-XXy2KG3pXWG?J;DCZ~RCeRQNVB{d-Ix=I zbUFscM=^-4(ySGxB4n@uQqT}`evE&icAv%;3jy)c`e#C4&HLgR{<5F{e6>^hLx1{0 zahiL4bMJ4}`z6Osa>jmbz30K$jT6*P%?`4;9-GVBBH>1YLQ2hv$W;SLI$rJQt?6VlC&>LhoJ zvoU2_%;DbL7<_E>KA}CXVFNmNKj{Le7_VDLEd~fc;UBN19=>(KKEcHZkcTdEVI&sb z6mo<$Kg9WxiH2C)picB`{hWO}g0p_h|Ka$i4@1hmNMyPPc=;evJzptQ}l{U%*QmUW9`T8tdNwM8DSz?VBwgC9} z4b-V8IPP0)GAaLKWjpjpdK>^Th$Uzkw9siB$(W7nCMf8MtZwW*Sjxt+ZSt#5T`x~6 zH&@AnZGRnVsT&ks`ST3vB-K~O%$+D7`mvPJsVuuyUR#x6P`=51Me+l$k53n<{`tx3 zJC1KcQa*G95#>zH3y!8|85k;JX%XhjO9FOmRilO4tSpzAF2fGmHb)bo27c4*ZLZd% zCxIJpMFYq&rGc&&v8VM>CoW=6{7NC4f>!XAVVyZ# zWwBhVxr!cHc(^%Oa0XR(P|BBeq11EUu%x1>aKm;r_gmY;E=PPkiyQsd;I-1tQTQ3CIF@}YU&o!&f+*UWY9Rcnx2a^O~akTl1;4qP>DKCTO^+zCp2Sb z$bhVXBx!iCfvWiGylT3PzHU zYAX}qK7~e@V9kHH@*-4M1zC%HEs?Dn74skrvPY_&$B{uOjFp%vh zZTuN7o*Uz;$Bg%=Nos4kI@Na5<74IC4>!*<1)#Bb@kqDaVkPO4u}v^6h=cu5VB{kC z!xIeT?7-I?quQ9QQ^!FZsHJ2^U}9kl_ync!w^#%~fEGIxf6|!|-MdE6Ts@JTW*tj@ zYv}aDCAP~~xz&1p)!&~`yY168r(5x8l8} zyKqlj;oC)d%uMnM2K*kdyNFWX^~S;Oy25-D!>;V1J@u7aJ@ed^hdU z+O=xNS#A~DS#Igp_7Qj%i$4lU5|AK^l_=-K!!hj7I9?kmLPDAznmuEKbFda@(A9{B z%;&&TXk0@t9A*w)DA>q@po)2--|q)q((isohU1_uDTJzjzaMBdTuRr)*R?94tx;I_ z*oSXEL4euCX+@%&02Yq{!_~vTeqjq}DGz&-V>sd7$EhSn$YfHt_27gh7mADV#+%1B z1rIaDRli3ebum3VW!a*V&fuO70f7dKC>e<_?^s8G@s)j8`=b6WA z8?*U0&W49ip3_(Tw0;8oDXn>d>c47fE_#gI0W_aG+R?li5f1rOLd{hmd>>J)t!QP3 zRv4S>7yzo1iM2)_5%LL0v%%;Pr)mBulvoYD>yFbRx@vT}uK*N&o0f*1emz-&2`<19 zL9Pl~HfxtZ$_!kc)%0YDHEE#efmDq~$ASc(l5w`kC>8yH9bDo7fJrh6maRvOy8hEh zXYXj{kq=kRPql|FOb3ncras+6D-q<46`KnY#SZxVy`O{|5;*f8;_s=O{L^~QpC&#UdFUA9{1$n90S>!d2()Oc*Is>VP6g=1Ed zUEqSLYKIK@3E9J+kn%Yk5aPNpK0|A&Yq~M19LN#3=X6&U|7qN`!|arTYm-#Di#E=b zOpr)_8M);-L;=!Ocr<}a(1^u+o82EYX4HQ(|0TfbN|5>OBeh;ggFn2|e!@oUNj4M7 zF$cOff=GY&9xC8H*!^sym-XT74@XvY?yvlbF=n5#$hm?deEt4<0jyPv>qa%xg( z%LjE$N>K)Cz~#VgJi!6tcSmJuOX^>!h$0+p&2S*MbFh;--|>w~IOg^6!sYeUCaHZh z9JGsnz}UH@;i_c_k*vkKzRxUd0Sf|BCL39*Bd zEPZy(5W`?Ev@%gF(D2AgNo>Wp+71YM=SxMNuh%hhewZm4y@dHR_TBQ>Zl;`-AFmZW zu*vl#GgBg>a~*DnmFp(->p7~LC6u|1r)kg!aiZQQ=6NU3Bc0h^BQ=uj$6d;Naw%Tv zS3!ibUU9g`rS>lK%J6V@(8--0Kdlzn)&Sb9>WSb`RaO>|_RH?qyzA*7aa-aHx3Kbu zv9^=-5?t<>37-o5xdRSAA)y5o45_H1wjWq3^VIWx|7kpZ+7;Uds0?xR!Nlf3Pr)7N z$=cE1!hN%_>ZnSromrSHt+H}NI_o&d4@r=jZ`^n^yye#M%NB)PanT*a+s9(|Bm?8P zXE>@%e-b06_ic&t2lXSi3EmjwRCI~Qh}{U8kB(ysT&+mcr0bCLARR<0^Fc8oh@avAZGp#jK=J`Vl@uO8CxF0S!fx2}Rc)8ufinr?v zzo5tjd6R_|Fd2Vj$^ILsY9aHc%NoL?TUa|5ezHuPEMjtTkbA6=kV8j+{>gqp=%g@0 z04)9i+$w>8E8*x&?ML5#nD*oAQyaN<>5Z@rnF`6=Jc-trKuecOfy zcN^?Vi3sbbm^Ckxvfpe4Q*de3$nkWb<_WenopHJL0k;@Qet2&|%%smwa7hcX2HQ%5 zIh!y9R;ovGr3bu>=+=C`M#b@E0@Wqu^ZT0*PQaO!ELhaB?T80_ z9(g}46Efh=Isl!2F{6^u4f7w^H6cgrqiovW#joti+ z^~*YhGvxZhIGo4Ys13>0(AfO|IgzbZ9eli;JfCusEr$xX>JqPefo0OC$Rl~h!I6_B z+z3;XN^z4dt?EJkQqXSUm~G}JWql~8Xe_+tZxvLht#~-HJT_g!-C(d1-J}N{m)ZE~ z{_H)_m>`kEf9cQ~vGs4f%!;o@!cP=BK%`pCOFRv12){unTruyZf#V=2f}?w%d9H~0 ze)Noy6`TYTBe5(l8p%1xk;YHZv7z@n+UACB0=x~3Y?iHXBzg{~Ue^~kpk^wjU+r59 zeD9;1kVoy-4b(Cbf*}7rI?W^$GiYc=tKsFuL04Z)J)3U zp6fwDZh+Hy%zD+4_!au?MCDV6pw(A*f`D4e`;G;v@#PSKPVDirKDG|9Z7cc)kQ;t4 z+|lbMqLJrsPlSs7Qslt8Gz~x!8391>9C#y^_CcCyB<9dS0(9IkoUV(#GV$p!{vL@h zQ&bbvP2OJm>jCD}fxZdYCz;Rh3Em7vr}dqM8jQbCbZU4}ug_px!>{L;?>|4jV43O_ zm+*2g8~yqxO#Jw z(!o{2c(L_fUPC4~9Q|(R^o!l5=sL|aryOpD{cyuPI^n_d2RDA$wBm;2DFLDu2L2v! z|3LWuD=d2-zWX8M6qk$cqyNLVE&AK=nsaa{`oZA4V>MEzj7ciI7mh#Vkv$ocA@KA@HcjXdor+zo1i%P;f9vXSaCcukvxA#tWesHDZl51CEh?iy2 z@L}|}g0J$I{8@`4w`FJ6;cKW6654O3y{z3zvlt)s>sp~SRQ+i9z894l*LuI}BHXcG z&a!?G{yFe}B^@69?SD(z@Ra{EjrISl8>I;mu7wxNpQ1jTFiIDNm4|mf;BY(QvVvXd zmXt5QURzNE!;|;~Wo)dV4Wr`4>+YM;U~EMqVjq!u>-WchxS;&U3r~)j8xWHkxcUA^ zfYym5$9D7VfR;1J;%=HF$H2|#_eM`*8?xJ-*$|yBFTStbzZ3^?d$O$)9X#uQ*2a1hG*N@h_ykZFI2iU%5@VlNNm=mG|g%9XRyt zfA_%t`was|FN(YVYk57j_2jFt*_Z}FQXB@h2&k|APWex{UZ8yIR6iqa6pLnUVes_J zWrnH_5OtRTlDc|F`#?D#St@S&awQFfl69!9 zW}3p;Ni5lLnrKgio8iYTl$_BhoO+COw}(0)E<4f`2hE);wwfu=E92zEK?fiZQ|Y!t zbBT7o!Wc#G>D{fd$a0z>_F4_+cq(f+9EgKzz**VUsbjHM5ACuCNZd*75KsyxS%7yq zNQeUqj$=&4D=kEd@1vK=uO16H?>>Hm$>ZMYNBx}vae9Ls{RGO2Cge)~6g@FV#yv!X z@7FDrF>*6q+y$g&f*B~gf2=qXo_z8%41Xi-Stn({F4TE?2xwW zSNo~G>DfEr(Z1@F_}ks}1Ek`Y-A>R4w|;MBwA^fyXKNnmQD^dcmIH2jVtR1s;H>B| z0FVN}Jjh#Z((b&H7c|z7EfE9B7PR3!FZ~7@nA!<6)`2NcVuyAYtQ{4ByUmqM;$SSN zNT98qM@y?wMn3{^cv~uzW|PGMq>&0V;#yxaqssr=m=Z z=xZXSbpZmIPUV-VTouFt;0#V4?CVw~kvImKih<;jcV53ukOxNC3#@Z+l$<}33Tkj$-djxjfMUUyZ@}y%=%Kc@z=6=0xKra)#r}yp5-&j4g z59N?yLl(4tPto?vr*%~VYDC<3+=uy;H0o5J!3yZ^9ax4(7zhn(407Fk6rncLV3?YZ z)oQL&N~PlP4`QikWYqih&lIOZ(>Y1FNn}m0HrWj^+e^_2{5`fp2>CtRFA z*isYnty%SLtenN?`Zmc+NP|_COv19pRp2P{Oj;>Vy!xcZpS!qqc~Wan9{Po=qMYye zr)+&g=SO)bTI{B~Kw(qlYF>aQYz&T}@q-5PuhF5txY!#iG3hC-U5|c=vB~nnlT6C$ zP;pfT)vi^sS99=v2SD;_l1X&>aI6JN9Y8kwY7?55c=}CH^TnR>>FmT^+jQ&H@~m`U zH2h-&Gp@bX&S4ZWfF2FtT^0xibXW#@thM>17-JA`!XDh)OSg`L^_B{cNGvKkd!_p# zH#dMsYJe!8EsQp$9o)14<*uwY;V6fEwkB#j(QS)COvwi?#*D~RfUeReq~S_bUdqAf6VQmh(U~uxwM}Pa?bkKAA^ajZ>i#+8e0|Yj~`a^YMw! zHV22ScSWHjHju=?m@FR-(T260`Gy@31#|DBJzNmiUYy2aF|q06Vo6Lo%`meX_p6sm zIs=>%fyodx7$E9CnY_EqUo``}JUe>BCxdQA@*gFry2h`sj~2M%Sde%))sA)i1^fX! z)4bLjK(V4YfOo&>1zLT%Sf=Q|TkZ7*%f$flmO!#j4Usne9MiXxjSk5qwk|Su+UD_L zN%iX4Ruj!5#oVc%Qq`H_T|g4e*aW+Rw-!*TNLsQSUA+DN<&>d)sQLpN1DE-OotjoN z1ma^F)|zLfsZ`GnVF8_21AuV>;A%{Q5+*qf6uvahSKw||>U{l|v@h;Nuf#pW63<5f zHmvOdCj}ph{c3iS11a#KKu5BxTAmF?oKU+U&vl_zE5XjIRo}^^cIadaP$UfioQOQ+ zxwOe1FtlGvN$N$y|YzRu!J{l<0^ z?)uBFzBrqtlqL`ONviC5`RbrIn#Kyy?x?^Y#key%6b~Awr;^U&bMnnM9p8yw`49pl zv1F%jo3-bO9@Kub#BalG#tBC3+fVa=VPPDvyU(*OQW%T_L2#Hh=-zt?Z}Xr9*sw0j z>CAv{?_knRZnw4u^#pmdD@k)s%j3}Ypz^B=H^u>gm^A5NmzJPlMPxaNC}~=o2Lh}} zOxb~z@@2i8=b++j#dC)-?M`~*sCe+eBKA-NgRH z5A~823VQ6EOg6~X-|P!Hfdha8lo7>ei+K_Z-{8}?^I+~bN6cX-!#5X9x07hUEGGL9 zj_j7M@`&Un9Hxg6>Wsd|J=5GO&*Y8mlVS0DrY;@qA;k!%MG1n#Nsc?ndVTQ9Blp2^ z)KKnN0dbNaBR1@WblK+)vL7CUVyr&!41VUbdmlyg*;QdA5OngYkcmuzL^u(qHok}e zo_^{_uW)CGC8r0UZs*r_>_zWkFFHgTwE zC#;p2_TzMRKELIv9q~euTy?w_Pj$Io-F*X>1V1ym{OsE>YRqevW15ARPiA1*p4Mh!kCeF zmE}<%-*k|T&w=Wr8GP6O7L!lWct5*()=bO!B$6F>_mmL$rHta#oI-$M-Uy&yMp8M%iFEx3n;k1logFCBrx{V}EM z$?l#MjdH}^d>`eV_|lw!n27|R{*~Ja8KHr0(pK*y>YEszPVtJ`Iw#q^Yv5OYbW@&{ z?<`TNf3G&(aRa8r|LJJ2{~nCDvT07){}_*%L~*(v;Iy#_JNk0vxanldT&jtrQD)yx zzrZJFKb1d|yP#*ud|wYZBenKyB08c&n2B3&)CLOUO^PtBI)}Mirit61lpFQ!>}BV% zWv73?HCGRpeQdjQFq?iAxSh0z%1nQVjqSWrSc6j2>Y=GH5-Qb|-4fjdF=LbcM+S0Z zYiVm9I)<8L#W?05Rm{XwQhX`@s#{qW-?r~DvYYh9^{cSLG%Nb?r#_3wp1V@zz<9mR z2GI2>tZKo$rkmjJY;5yr)z3}9Um!6&jUcW|S_(X#!x{8AuY)-~_}r=_Kp@TA&~AZy zs)x=0J{uL}dXO5AhG)Smh3H*bgYg`0PL2^4YbZo^&BOsYoP_i%P+hkZEZ>S2H`zL2 zo6D`76JzBTE^JL=^+h7nTsSHUqdKDml-@;n@^sybO{PJBwJ+f+2B!7 zGS~|>&@qZDn9b2#+sWTGQNBZMZ%Dy8xy3t(kne5kPfagDkCjvuDkgc}t@Y?up70wk zASta58-pGk^gBqS?G&=r637MVU2HxjPQJNS_E=}f`?8Ma-RQDZ_5O8^to~Uej>Y6% zw-~$571_-#J!u}Ht=3QN5Y?yTRsf%l!aHc7?6`p&6^69?L&Dyq18hZ^BGsvj>6?4k ze`7Zlsce&zm&R}1iiC$vJ63?Pg;65$cW+O3^2`IUiBmU91F$w;+u-Wn#`=5&3hKkz|Z zl$RB&w*g4j;wyWFw5{_D6Dj+=5G4K?S3A!->np7KBs(N&@Ivyn+z__3Xoy9}SJr__ zPX|0Np3(E%a0+pA6JXk|sLuU!#MoOQ6|DoFi7VwR+uzVx`RXwo6l|i1E}nXk806Jn z6AQtXSnx_OSRl*rhQUF)E6ogK)53k}!H4G@io90~C`?j@9B2O;33-ONciDOvp=Hwg zN^9|RWf^A)34AQBdPqIW_xEJgmyiubgA39jLig_qtgiW1uT9v>dQUyO5`UVF8((#&v*GorE{8>=*(DaV>06||78U&`U8Kl-gVFB_h^YNmU1d_gBvyWEtyF=a?0?VVSl5)CySjRv3)vW?(ko2vevQ59SVz-ZSV9X5=a-zVVyi`<; z6jZYR*>}ILTi(5K-pll*nve6EtUt*Mk@>Y$VVl2xai#pnEh-0;q%QZ01yTeae3N1^|Kp}{2E7BN>$4;3?RWrP1zo4>{`!BlA?J8pU4{;e1QrB5zzc3 zTQD=lcY#Y-Rj=s1P~3uwxHQU|Cl)GE9JM7^pP55aYpzPuU!64ii4sClUJK+JH2oR*Yp zyek6iCRaJSc4ne^v^yzVf(LlC?|FR$GSu8BwJgitnDDGu3sroD!6z^5P>|u)0A6m} z$T|vG&0*f%g0Gu)+P?xK*iMSC1j0JN9SW?GY;YJ<`SJ;g70FhWe3nk5P+ITxw9-nGvp9&&XWj8TiCOal3T+AsUK2qJn;`I zW92x9?B5D%l8yIYz$3@(q&GnZaR5|c=gnHgCZc)dH#qgPw@&~B>?(s)RBMxXp56s>Vy5^R zHep^tL&^9U(Yz%Ocb7kcy^ktEv@)#@jX2+V>yY|VCrAto3my6t0a|>CU|^d7Iv^ek zDntSfmm(V=a8T{GThlx^hBFTH(S$U8ng=ig2N%{oSB+8&{a#}i5XO5({jLR$L#$$E zwFHNM8#-dK1^1V<#JhlxwxICW$i+wx1*BnG0npMZz#-zbV=!q6YB}EH|C#f(FO%14 ztO21G`(qWJH2!xN`t_5cg95t4sJs9R{`G#XEjO>?p&Ony;*!_}Q?QDcny`XUkr8e{ zR6$0!1(lcv?Y&}f`3B&b2J)*Hc=QexN1*6t0ortVSDoDxRlYB-I;(W#}(f4`+o|IRvj$PT;}4-9Y`g6EQ2yaf(ID%073zw~TVUx>y@ ze4kd0Gu5~o0(j||1FA$mD|e@1NVsMgt7#>$mwxRzaU3h(B)L=C1t0=vX9H2RhIeg% zjo1CyrN}D+S-w?_Bc=;_mAqoljN!6tK;Hb&-SZ1O&4s%|6;qn0Ap)cy60y5?g zPGUX8d3A^(=kEcL3lwOVh5L!K3=Q@xVhq$gR^4n0S->Vk=8iFGO!Jo|F%0+73kaD;B(B1}g5(XHadcrN6-GoVb+O zC{&%&9Z+4Q62_%z&H~7Q_6}O3cOaG16KNWN1Oup$1JRM6jCa%ExL*YjI8@FaH%3y2?Qku;gocrXecz!q@c(hyvrz;SfAnP&sZg-J@ z5Lgl}658YaXS|o(^{-%RtijQYG2pTa*F@NRu)oLdvHj}*8U7#X0(55%^{3^(sWbV| zkoA$u3^`BGN8~Q20=1k>^8_J=;Ue;9pL8)9_;sBA1+OdTdq>+^Y&{GAtS|Ua*omX& zQU6Vi)qiZ<{kL`0Hoy5EgoDFiPU@Xqxbf={73V|3qtfQts|0md>mF=Z!Qu?tX%?X} z%Fuyb&Y~F=Z5D_6JF9!86wA?h4-JcrHA_I=KrOq77zoF=IRHXj3n;rBwDd@@U^uJ8 zNw*5#N{W3}#PoHlwephpH>u(sDDzWjTROoNF_tz3R1lprF z0^$;E^h-hN0W$%|I57S*2%E{tjTsYomd#ZBb6rl)7->xga_|;-tgOM_X2Wh{Z9Rq! zif)3L%jsa-0rT>b`R(W*(~If7G}hdT8D4@34UVX7z=)GgIS`oCo;36mImC0Gt3O`n zMJMhaqVYJ^aS3%Vn~;=cB8lNNm+F^$CFihbG&``2*p&QZ`eFK?{Jr>B8CTG78} zLMun$AWJ(gqB{x$If>XrYy~v}fWe!10GE*unVsAR0J;eQ*bufNX@Jk-uvK~UIqtEU z?I|JsY={GRfpKl8yl0NTvxa=WJKZQ0dWRX7j`gr)a&*7QGlIlkf#y578Oq%kkl~bcfQo> zfbKREGTLOX2g%_1jo+`{?FpU#EO8Z%9!de}OmhFVC=Otz?Ydy2#gz+!+avH;~a*AiRE1C+h_yIrXGWbHJ)hx>%!0rTd9tcB{U>79s665T}Y%d_YU;zmIF6|lsgZLP#SgvkPMY>GuBAl8GN4SL{)Nj>9e)Z}&ct*@b zgymQ!%FMR|jb4lW=8nBrne_q+yG0HWiE(6I~bY07-7S4~B(BZf8 zg(oPzw0x%-Swmh6TvMTEL>--WPk7*JKMa%P7P{HKYy<}4dS}&lOpnD*Z|AAqf;hv- zd-fB+pF#F?b`Zt|(J<~b%$evVLDt>gYMsthPkMlc$#YO?D(dV+U#PdHIyW5K?0r^- zEY}A>BP^vf5RNYi1=_6Kxr{jtfF}bvcXZs(dI8h1?UY7d8rru8PE712NVXz1Pgs(S zM1+tK<(KS4cI){|2p@`RlTmTiA2()=gqp3&7Xon{vmH0=qT=YzeUh z=r0ean?Oz6+138i4TZR2Nb1QJ&Pt^#3RKQb=zm{Y!o#B%dWY)wKZuL4^V24h87iXY zGyW-EH6?1#-Hc_gtfv7$f({_#_@eIm$Ja526npN5h3D%YJ+fnmM^c8fiB5G94y38<>}cD(O54FZf2t=?CV} z{vY<#pt@97T&(g5iRZB{f4HT0FCfI>z7*_sm-0}a36LeMBn74+{$ z0mIw$o%OX}t~Z|Ksbm*Q=N9qL-+tcw(9{v>Qc~qhiXBK_rI7&7C$6Zr-MIwYne0j$ zv8)rjvl;V%QGO1+a~hHj)B9JjaN_WliM+Sob+$wLeMRXC6n$jL9@j zk*NLt*8_ObNCL|>pV2M)#6O9_`1+QU@?G{+?Y*L`zeMAgg`I!h+6{9@YeB!yeg$y7 z>i(rNXGB*?&GL6#(Kcm@(S_!rDh1?ZSdofS+VTKoVu_&i>yBjW3`&_b5?TpqR+&u_ zG*gtTdL~_*9oT<^BRWYz;b;ZY4O-xTlk{w@SoVr3wJFH?oSq(GWV`wGtJa>x)>@NS*Y8Dlx^S zZgsg>6j^mGHK?5=BI-vLVU*@xfE&)A8<{U<-Y)5pY-;4w|ZEPZxqW zA@UsWTe8pPA$~m^|7O%wXD|D6$B*itn08!T zDa6_f!{p-Y@6gYM<3Gd_#KsE+3n(>3J4XVL@n}KSD%6UU=QZLqy&)Lor6O!tl7D73 z-5Yws+c<>i1~E36NiW&mRHWaO57x^$<=OL+uOIDbCFihUrAl(Nm+1{^iKKxy59*PI zYBt-ER8fhY?DYd$?aL55E@ljJ7@(vPDUMT)-R+3U21f|jQ%FEpx<(w&bkb?;U}2cuR#TUf{J`OnJ>MESqX5P6v4L>X@!=n^j zZqdRrs;x&xKos8*mDk1F3e~tcllI`v)ZDxcNE38^S`Lhl#8Ee-MaKbJFugD3ZK^mo zOxorDx3C?A;J1_x30M+aaR3+dfsv4aOY`qS--Y#@4??v<-;HD>M%yB6yu>(87J>#a zHh3R00ZyuI0Zf4S<5~uJ^B)xyBoobi9b=$`aD>y*dTTo&d{%yhMwtwmIMj1RqU?#A z8-U|t4yCHk-ea8+mNnB8k&v8E7+Xf+RL~+}VTp?c(Lu5v(Q@@Vd>*^v?VJ$Q~9f?CGi0ay;*TEQ@t4%ka17u&o8_K4e(J*d}Bs zX6;L<`TZ8ga-!qN1vohfB{wZF<^~uPme<1C-9Ex0atuA%van$$<8xp?xSi`sT1gF1SBf4fAH~mADaog≫v1~0#ixCNyreNJHL+QLTsE{sUmmV5CrF#LcBKx`r5c0c%l@G%i zx`s3!pv4?(IjhGXwz{FEF0$PGm$z`Eqt9GDD-2f~zndL=vPQrqz8@X3%h^TN@pw1T zIZg@+vMF5{KR3Vyci=dZDh7aUB?C6_su+{fI6AOlt{shG(_ntlNGhX!#K%*vAv@sI z%@%AW1qz49jLz9mxB%*K@BQ0;;lOei;6&VPAMwc)+!8{&&|Wts(eviI$#1t4@ZT{O z7l>i^LJ|4!p2PZW0Fk@UIL!8T0=ub^gB;fx9lR zPMR(##Q)rN$G`uEHKz9(b;yc!GGx{$r+