From cae723eb0e82ea97a1adc232d60e0ab66ab7a52b Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 24 Jan 2014 09:36:01 -0700 Subject: [PATCH] Moved parameters to separate file. Documentation updates. --- README.md | 18 ++- SimPEG/Data.py | 3 +- SimPEG/Model.py | 2 +- SimPEG/ObjFunction.py | 85 +---------- SimPEG/Optimization.py | 3 +- SimPEG/Parameters.py | 164 +++++++++++++++++++++ SimPEG/Problem.py | 4 +- SimPEG/Regularization.py | 4 +- SimPEG/Utils/__init__.py | 77 ---------- SimPEG/__init__.py | 1 + docs/{api_Optimize.rst => api_Inverse.rst} | 16 +- docs/api_Parameters.rst | 11 ++ docs/index.rst | 14 +- docs/simpeg-framework.png | Bin 0 -> 16020 bytes 14 files changed, 218 insertions(+), 184 deletions(-) create mode 100644 SimPEG/Parameters.py rename docs/{api_Optimize.rst => api_Inverse.rst} (99%) create mode 100644 docs/api_Parameters.rst create mode 100644 docs/simpeg-framework.png diff --git a/README.md b/README.md index c2eafa26..37e0a743 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ ![SimPEG](https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png) -Simulation and Parameter Estimation in Geoscience - A python package for simulation and gradient based parameter estimation in the context of geophysical applications. +Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications. The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features: @@ -10,4 +10,20 @@ The vision is to create a package for finite volume simulation with applications * supports 1D, 2D and 3D problems * designed for large-scale inversions +Documentation: +[http://simpeg.readthedocs.org/en/latest/](http://simpeg.readthedocs.org/en/latest/) + +Code: +[https://github.com/simpeg/simpeg](https://github.com/simpeg/simpeg) + +Tests: +[https://travis-ci.org/simpeg/simpeg](https://travis-ci.org/simpeg/simpeg) + +Build Status: [![Build Status](https://travis-ci.org/simpeg/simpeg.png)](https://travis-ci.org/simpeg/simpeg) + +Bugs & Issues: +[https://github.com/simpeg/simpeg/issues](https://github.com/simpeg/simpeg/issues) + +Code Snippets & Tutorials: +[http://www.row1.ca/simpeg](http://www.row1.ca/simpeg) diff --git a/SimPEG/Data.py b/SimPEG/Data.py index 6c94811a..948b1c83 100644 --- a/SimPEG/Data.py +++ b/SimPEG/Data.py @@ -1,5 +1,4 @@ -import Utils -import numpy as np +import Utils, numpy as np class BaseData(object): diff --git a/SimPEG/Model.py b/SimPEG/Model.py index 8482edde..33cfeef8 100644 --- a/SimPEG/Model.py +++ b/SimPEG/Model.py @@ -1,4 +1,4 @@ -from SimPEG import Utils, np, sp +import Utils, Parameters, numpy as np, scipy.sparse as sp class BaseModel(object): diff --git a/SimPEG/ObjFunction.py b/SimPEG/ObjFunction.py index fb2fe2f4..ea2fe831 100644 --- a/SimPEG/ObjFunction.py +++ b/SimPEG/ObjFunction.py @@ -1,11 +1,11 @@ -from SimPEG import Utils, np, sp +import Utils, Parameters, numpy as np, scipy.sparse as sp class BaseObjFunction(object): """BaseObjFunction(data, reg, **kwargs)""" __metaclass__ = Utils.Save.Savable - beta = Utils.ParameterProperty('beta', default=1, doc='Regularization trade-off parameter') + beta = Parameters.ParameterProperty('beta', default=1, doc='Regularization trade-off parameter') debug = False #: Print debugging information counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things @@ -213,84 +213,3 @@ class BaseObjFunction(object): dmisfit = self.data.prob.Jt_approx(m, self.data.Wd * self.data.Wd * self.data.prob.J_approx(m, v, u=u), u=u) return dmisfit - - - -class BetaSchedule(Utils.Parameter): - """BetaSchedule""" - - beta0 = 'guess' #: The initial Beta (regularization parameter) - beta0_ratio = 0.1 #: When beta0 is set to 'guess', estimateBeta0 is used with this ratio - - coolingFactor = 2. - coolingRate = 3 - - beta = None #: Beta parameter - - def __init__(self, *args, **kwargs): - Utils.Parameter.__init__(self, *args, **kwargs) - Utils.setKwargs(self, **kwargs) - - def initialize(self): - self.beta = self.beta0 - - @Utils.requires('parent') - def nextIter(self): - if self.beta is 'guess': - if self.debug: print 'BetaSchedule is estimating Beta0.' - self.beta = self.estimateBeta0() - - opt = self.parent.parent.opt - if opt.iter > 0 and opt.iter % self.coolingRate == 0: - if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % opt.iter - self.beta /= self.coolingFactor - - return self.beta - - @Utils.requires('parent') - def estimateBeta0(self): - """estimateBeta0(u=None) - - The initial beta is calculated by comparing the estimated - eigenvalues of JtJ and WtW. - - To estimate the eigenvector of **A**, we will use one iteration - of the *Power Method*: - - .. math:: - - \mathbf{x_1 = A x_0} - - Given this (very course) approximation of the eigenvector, - we can use the *Rayleigh quotient* to approximate the largest eigenvalue. - - .. math:: - - \lambda_0 = \\frac{\mathbf{x^\\top A x}}{\mathbf{x^\\top x}} - - We will approximate the largest eigenvalue for both JtJ and WtW, and - use some ratio of the quotient to estimate beta0. - - .. math:: - - \\beta_0 = \gamma \\frac{\mathbf{x^\\top J^\\top J x}}{\mathbf{x^\\top W^\\top W x}} - - - :param numpy.array u: fields - :param float ratio: desired ratio of the eigenvalues, default is 0.1 - :rtype: float - :return: beta0 - """ - objFunc = self.parent - data = objFunc.data - - m = objFunc.m_current - u = objFunc.u_current - - if u is None: - u = data.prob.field(m) - - x0 = np.random.rand(*m.shape) - t = x0.dot(objFunc.dataObj2Deriv(m,x0,u=u)) - b = x0.dot(objFunc.reg.modelObj2Deriv()*x0) - return self.beta0_ratio*(t/b) diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 2362dd20..4e5963f6 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -1,4 +1,5 @@ -from SimPEG import Solver, Utils, sp, np +import Utils, numpy as np, scipy.sparse as sp +from Solver import Solver norm = np.linalg.norm diff --git a/SimPEG/Parameters.py b/SimPEG/Parameters.py new file mode 100644 index 00000000..2a07aa69 --- /dev/null +++ b/SimPEG/Parameters.py @@ -0,0 +1,164 @@ +import Utils, numpy as np + + +class Parameter(object): + """Parameter""" + + debug = False #: Print debugging information + + current = None #: This hold + currentIter = 0 + + def __init__(self, **kwargs): + Utils.setKwargs(self, **kwargs) + + @property + def parent(self): + """This is the parent of the Parameter instance.""" + return getattr(self,'_parent',None) + @parent.setter + def parent(self, p): + startupName = '_startup_paramProperty_'+self._propertyName + if getattr(self,'_parent',None) is not None: + delattr(self._parent,startupName) + print 'Warning: Parameter %s has switched to a new parent.' % self._propertyName + if self.debug: print '%s function has been deleted' % startupName + self._parent = p + + prop = self + def _startup_paramProperty(self, *args): + if prop.debug: print 'initializing %s' % prop._propertyName + prop.initialize() + + Utils.hook(self._parent, _startup_paramProperty, name=startupName, overwrite=True) + + @property + def inv(self): return self.parent.inv + @property + def objFunc(self): return self.parent.objFunc + @property + def opt(self): return self.parent.opt + @property + def reg(self): return self.parent.reg + @property + def data(self): return self.parent.data + @property + def prob(self): return self.parent.prob + @property + def model(self): return self.parent.model + @property + def mesh(self): return self.parent.mesh + + def initialize(self): + pass + + def get(self): + if (self.current is None or + not self.opt.iter == self.currentIter): + self.current = self.nextIter() + self.currentIter = self.opt.iter + return self.current + + def nextIter(self): + raise NotImplementedError('Getting the Parameter is not yet implemented.') + + +def ParameterProperty(name, default=None, doc=""): + def getter(self): + out = getattr(self,'_'+name,default) + if isinstance(out, Parameter): + out = out.get() + return out + def setter(self, value): + if isinstance(value, Parameter): + value._propertyName = name + value.parent = self + setattr(self, '_'+name, value) + + return property(fget=getter, fset=setter, doc=doc) + + +class BetaEstimate(Parameter): + """BetaEstimate""" + + beta0 = 'guess' #: The initial Beta (regularization parameter) + beta0_ratio = 0.1 #: When beta0 is set to 'guess', estimateBeta0 is used with this ratio + + beta = None #: Beta parameter + + def __init__(self, **kwargs): + Parameter.__init__(self, **kwargs) + + def initialize(self): + self.beta = self.beta0 + + @Utils.requires('parent') + def nextIter(self): + if self.beta is 'guess': + if self.debug: print 'BetaSchedule is estimating Beta0.' + self.beta = self.estimateBeta0() + return self.beta + + @Utils.requires('parent') + def estimateBeta0(self): + """estimateBeta0(u=None) + + The initial beta is calculated by comparing the estimated + eigenvalues of JtJ and WtW. + + To estimate the eigenvector of **A**, we will use one iteration + of the *Power Method*: + + .. math:: + + \mathbf{x_1 = A x_0} + + Given this (very course) approximation of the eigenvector, + we can use the *Rayleigh quotient* to approximate the largest eigenvalue. + + .. math:: + + \lambda_0 = \\frac{\mathbf{x^\\top A x}}{\mathbf{x^\\top x}} + + We will approximate the largest eigenvalue for both JtJ and WtW, and + use some ratio of the quotient to estimate beta0. + + .. math:: + + \\beta_0 = \gamma \\frac{\mathbf{x^\\top J^\\top J x}}{\mathbf{x^\\top W^\\top W x}} + + :rtype: float + :return: beta0 + """ + objFunc = self.parent + data = objFunc.data + + m = objFunc.m_current + u = objFunc.u_current + + if u is None: + u = data.prob.field(m) + + x0 = np.random.rand(*m.shape) + t = x0.dot(objFunc.dataObj2Deriv(m,x0,u=u)) + b = x0.dot(objFunc.reg.modelObj2Deriv()*x0) + return self.beta0_ratio*(t/b) + + +class BetaSchedule(BetaEstimate): + """BetaSchedule""" + + coolingFactor = 2. + coolingRate = 3 + + @Utils.requires('parent') + def nextIter(self): + if self.beta is 'guess': + if self.debug: print 'BetaSchedule is estimating Beta0.' + self.beta = self.estimateBeta0() + + if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: + if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter + self.beta /= self.coolingFactor + + return self.beta diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index 501b1c45..e53c4675 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -1,6 +1,4 @@ -import Utils, Data -import scipy.sparse as sp -import numpy as np +import Utils, Data, numpy as np, scipy.sparse as sp class BaseProblem(object): diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index bc754592..aa62449a 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -1,4 +1,4 @@ -from SimPEG import Utils, Model, np, sp +import Utils, Model, Parameters, numpy as np, scipy.sparse as sp class BaseRegularization(object): """ @@ -23,7 +23,7 @@ class BaseRegularization(object): assert isinstance(model, self.modelPair), "Incorrect model for this regularization" self.model = model - mref = Utils.ParameterProperty('mref', default=None, doc='Reference model.') + mref = Parameters.ParameterProperty('mref', default=None, doc='Reference model.') @property def parent(self): diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 8853500e..437f1387 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -267,83 +267,6 @@ def timeIt(f): return wrapper -class Parameter(object): - """Parameter""" - - debug = False #: Print debugging information - - current = None #: This hold - currentIter = 0 - - def __init__(self, *args, **kwargs): - pass - - @property - def parent(self): - """This is the parent of the Parameter instance.""" - return getattr(self,'_parent',None) - @parent.setter - def parent(self, p): - startupName = '_startup_paramProperty_'+self._propertyName - if getattr(self,'_parent',None) is not None: - delattr(self._parent,startupName) - print 'Warning: Parameter %s has switched to a new parent.' % self._propertyName - if self.debug: print '%s function has been deleted' % startupName - self._parent = p - - prop = self - def _startup_paramProperty(self, *args): - if prop.debug: print 'initializing %s' % prop._propertyName - prop.initialize() - - hook(self._parent, _startup_paramProperty, name=startupName, overwrite=True) - - @property - def inv(self): return self.parent.inv - @property - def objFunc(self): return self.parent.objFunc - @property - def opt(self): return self.parent.opt - @property - def reg(self): return self.parent.reg - @property - def data(self): return self.parent.data - @property - def prob(self): return self.parent.prob - @property - def model(self): return self.parent.model - @property - def mesh(self): return self.parent.mesh - - def initialize(self): - pass - - def get(self): - if (self.current is None or - not self.opt.iter == self.currentIter): - self.current = self.nextIter() - self.currentIter = self.opt.iter - return self.current - - def nextIter(self): - raise NotImplementedError('Getting the Parameter is not yet implemented.') - - -def ParameterProperty(name, default=None, doc=""): - def getter(self): - out = getattr(self,'_'+name,default) - if isinstance(out, Parameter): - out = out.get() - return out - def setter(self, value): - if isinstance(value, Parameter): - value._propertyName = name - value.parent = self - setattr(self, '_'+name, value) - - return property(fget=getter, fset=setter, doc=doc) - - if __name__ == '__main__': class MyClass(object): def __init__(self, url): diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 10d6d4e8..4967608e 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -10,6 +10,7 @@ import Regularization import ObjFunction import Optimization import Inversion +import Parameters import Examples import Tests diff --git a/docs/api_Optimize.rst b/docs/api_Inverse.rst similarity index 99% rename from docs/api_Optimize.rst rename to docs/api_Inverse.rst index 462ae0be..21f4386e 100644 --- a/docs/api_Optimize.rst +++ b/docs/api_Inverse.rst @@ -1,6 +1,15 @@ .. _api_Inverse: +Regularization +************** + +.. automodule:: SimPEG.Regularization + :show-inheritance: + :members: + :undoc-members: + + Objective Function ****************** @@ -27,10 +36,3 @@ Inversion :undoc-members: -Regularization -************** - -.. automodule:: SimPEG.Regularization - :show-inheritance: - :members: - :undoc-members: diff --git a/docs/api_Parameters.rst b/docs/api_Parameters.rst new file mode 100644 index 00000000..e90b58a5 --- /dev/null +++ b/docs/api_Parameters.rst @@ -0,0 +1,11 @@ +.. _api_Parameters: + + +Parameters +********** + +.. automodule:: SimPEG.Parameters + :show-inheritance: + :members: + :undoc-members: + :inherited-members: diff --git a/docs/index.rst b/docs/index.rst index b2882cf0..5ba56947 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -3,7 +3,7 @@ :alt: SimPEG :align: center -SimPEG (Simulation and Parameter Estimation in Geoscience) is a python +SimPEG (Simulation and Parameter Estimation in Geophysics) is a python package for simulation and gradient based parameter estimation in the context of geoscience applications. @@ -16,11 +16,10 @@ these goals, this package has the following features: * provides a framework for geophysical and hydrogeologic problems * supports 1D, 2D and 3D problems - -.. raw:: html - - - +.. image:: simpeg-framework.png + :width: 400 px + :alt: Framework + :align: center Meshing & Operators =================== @@ -49,7 +48,8 @@ Inversion .. toctree:: :maxdepth: 2 - api_Optimize + api_Inverse + api_Parameters Testing SimPEG ============== diff --git a/docs/simpeg-framework.png b/docs/simpeg-framework.png new file mode 100644 index 0000000000000000000000000000000000000000..ad8e5c1c0d848cd29af422bc8e316944b0264359 GIT binary patch literal 16020 zcmc(Gc|6o#8|YYyN|q>OqG?fN9c5oaQjxXEzHcpdp~)C3QY3!KGFe7RLP%mT7$y6l zvKvg6nW?c%W^BVS_bdJ0-}~PC-ut=tud7cp-#OoNp7WgNJllEB^E|RJGv2%V;BF8I zwD;nLb5pkMN_4Z9ERE-4&;2?Baq>` zQ?brs?YbBb(8|THu77AWQKe^WA;SEp`XSs()_O#!KecbE#Ok7Kdwx+54}aAq%;t|A zYYJj!G{QWe8$jy+;V;FrJU>e;Z&1xTHKk}5==*Kv)zrE3HEQ0F+YVu5zpfOosR z=}0Nk?#39uaV+T|fE8$a0IIdBc#iu&6Qp<{OFNh>Ml*|n;D72;6kIA>5Ypz@MQCP+ z2mg3XFO;p&caB*(@AjqGQ}I@>!D2@pb4of$V6aVmvPLE*;qcmr+q1_Z!oUh2;8kag&c5{{ z6Y_jTFvb=L^nI^Uw@O5=hT{N2F+39h{&>Nf<<;1F!^W+4Xnziv|8U-TKqRb)< zyzy6%m$Uc44YER*?|D-f7-QfgGYJTV${766G7QN&Wn{Yp^j*wO>KAqxyy#DvquO;@ z7lwQ_|LMt}SXu4*!*n4errn|D^PIAwIK`*CvhG3koc*|Jm7VqFf9-W+HM7|~DtID2 z>et}lek04}E1PH>(oA6Y*hlJ8OAxwa zCqx+7nyruRFY8PpD2E_4*7UCa|1{$FcoDG})4#j^`#@{gtuX+|7~tbSzaW1=3PFEE z{?}KfIK=wjP@1JVBQINC-RT=%JLEr??9zD=5+X;yUq#9X6hl>9u@XiO;!6y@;nt1P zBunf*C`ZQg;_*mo z@~+0dSiBmH`E0Gn!YFW6&Qr2lLD7bE6=rs*H(A_ox-(kfm{0H_g0)?)BAFZVjW>F`BUp37HhRhYDCqX->${NoIyK?Y>8SiCyrJ&~`!9t<9EF!1Jn zHE0wVPd^iv`0Wmo6^zXgce)9eU$L*MSGfye5R&hPj2e(h98wqDTh@{2Wccw=pQahQ za@6;e*(PJCFSWiOV(=)Iql*Ngp+m}Qi~e|JUX(IuCG*8v@_iQ-&u+3sc}kP}sG|2N zDe|2E^qkrB;a_!PtaN(hQPb`-&L-O9`G^R-dwnMcEnX41(h7H(&DWk0hVty_4%&kX ziPTNb>MUbBnzT7~>IIMlPP=Ig`h7^OU{*1~=kyLh!Zg^d5+UdeBSwm>s+xlKIYDbe zDf6L9{K#wUwTK?c%>eck2DhPRu6^+3-V*RTIZ9+uGPJCfM-aNQ4@&TU)4iXgOdMJ^ zjV}A3Jzsp}UFq~(san+|LcoEKA$dJhCHYM{qPK9m0zzT}_Rd#HNiA}ouXUv^%u@vT z3y*s?w67P7o0*Z4nR0A|4l$rIzvLODg1^>5zKn5WO1VYGb90h+1Ro?Lm7rQld=GD= zAHnqBf3{weI9K31J-C2*oZI*+4qqi+1XVw6NbY4Z>^ptir&!nuc-r+7dSTkcd|egD z(oK#26Rx9@0D*#}?st#{@6&Ldp@F(_)j^r<>Lt4c&8j{AJ~S zgqr8*_++?psE=}{UzvZwN~@CR=3DtUtVG-zv7v^g$6g_tu}N(2n-47*2<{2$21nBv z4+Q>-#6EqGGer{Gl;Zy2hBQy1R-^%0>F)J>-guf74W8+cC|-=TNl&`j1}15s^a8j$ zgM3xI+OoPnzS}%&)qPvdknDaC;XIe(PUzx@q7@(UEEV{=tlFaFsd=K}rotiCyN%&n zzj}ET3hL3zR%xfMvsF-`Qfn)Txs9h|@cr%&s0eDd-K;2-XD?@$c%izKBLAa?3uAj! z70j-;vAkeNocTNOO0u}5sMLicz75x}+jzAhzEFZ(R1L@!T|)i!T`0AkAL@k-q&?=t z4)$2z1wj}V2HqNm{e1jg$;`u{0*}<3+6c|Q%~t$rC)Q$GjjctXW%-*bX|3{}15y_h zrr*Jf&D*lt{ivax9^KRH4SpXA$uh#B*tdq{VUCZX5E@0rK-xh1u;&x=Hohf3m2K6S=xkG;j2;+dlCIWOcnBXG8Bptno_jCUUIaQcwE9BzT1-CN%TT zsE?9DnL?S3s_hM52eg3p6VXjx6fdgrN*w;q1KEyt4v5NVC&uVYb!OC)K=HUwGfm%x z&dyYxMw6eM@y?~@-!R*aYK3~H-JOS_!(du})pqXJX%A&;ho+x-e)GC@89l5gIAipB z44yQ(Yz!xqAc!aAxoVOhCGX-z8UMJ%&rNzUJB)B(7$n}=XbMRr95rodiz2Lpk7W|; zpy^Vy%i{|CNxHs#EWfhjawLLHEMAK0tl;P8UTKb%?V{*c`f~{#DJxj+-8sOZd0^Fj z7{17PvhS+9zcGzoB?K)yBN+m8-M6KKAt~jj{kAo}An^OQLI!?rsp(iixo_JFK-+Jt z@)*;thAndes1?YTp8tEmZ=3NC1nV}D-zH#tz;Ao=4}`z%21NLu0l&v@HO$EXmSB5^ z|2p{=+$~%2&rpVOMsY`7vW(~aM*|pu+Sc14_Spr9!08l6T^o`OjG0ozn&4(<76VUu{7&3*>`3k4zylZP~-D2_Ae;oE-$B2LP}tSk1Cen6f}A_!H{uhb=dY4X<}Gm3-J#v3sNJGP^w+%U%1v`3yCZoF zq+N5xP|Zc!GM~Oj1kyi4pPzdXWs|HPCZteYS& z>`^G&6Tg<25K&st-K_gVOAyMz)%@aSPM)YuKDx@`Ns0QPYIn`X^E>7V8JL@khj>xZ zcPHO8&kiK^Idw>Lb~N?bWo1rHoP*65%xhp=7(G;r?OBl38FP1%Yb#uZGy{&u*|rWl zN*?R1RPc~yJafvUoo(ah6Pc3cyad;Eag+-l_7^`DG+~zIVOI+9x*ZQ}8^)8nseRbl z#54O_4L{7$wMXQaUT%+UJ=2FcY0%X>9XOlxZiIIu;Ay48(l%}x#W>2H6HCsvlSQN^ z%+ouAwFvav(8BF;!qTNzk|(u2P}ORK7FB z)}y2|;d%U@*Wp3^aOU1KL9n$CksA}2V9xpD#S4+Jg6-FVXP3hLqRi?=F{W#&y+W^D zUX5GxlpZhO-`?{bNB4_zjG2s`IUF8*Lht&X$F271BzPhE#3F5bg}Yl6VstLos~}P= ze3Vtzb9uim*{7Gj8~nZe-A^nWEjKimoOvsP`Pfg$^vk`$n$g=6D-BE@$q2rProikp zca~}KN417y06cCw z?%mMiddAE<16NZ1&_X*qAO>sOSI6A{Sb`ZY!E-x@_YGT4uX}EC*>>hVi~Lgk9P+r` zjMX2g%Ol-)5l=N$GAxm>(T0Ye)~K{QgO#j{R44E2Lru$nDAeH+v!8cKCz}6ZC$W>S zyBQoDK^Y8{!L`v9z1Os!u5;VC@QGL-?rvr+MJLw~8iPuCQ#LMT#r;&O*BKJ@vPffQTYjpCa_-@R%$)v+@V#BsKjh^W?)xKu zhcn)oY^xXpCH+5g7wdZOkBsg9_*uaseT^Sbhr3*J1x>e}w`2>Fk^4vZ|22wFDeqtFJ?y6%5k11OS_&V?EB8J1UPmtprf`720+y*ap?Mv8WZLUZOdG3?|Tg1 zG8X?1g|WLas%xDZqQRHuT^y_$*B0w!;O&de>YSA{AUiggYqzizqKLl8>`vyHK{8{b z$JF8v9e1=!d+K%(vSg}o4yh$9c7R4>>3C_U&W%wOCxjAE0ov(PRkQ2CW z-iA|fY(QZ|!_R!gdLCl^&_4_R`T$=tib_9jFudNET_cA{vIu`)?Wf99yHdE4e5?97 zEd)3)NU2q6p_*nBbRmp$n;YJp#6ZQR?ka!a6rffWk=dXpbxWG9qxf%FOF-9CfdiJ$ zK?2TSyM>8}AOstR6lyiuJ-{>M2@P%f@8vy9^XJ``o6pzu-QLvn%<_gm?lLw(WZB8|wc*c<2<-{5EEKWY@sOPo> z?$A5w=|3B^UjG2E&O=`Uw(U_`vK^WM&pR$S=ykE}kH!PS6*-rQXsn$UDr)7e5LB_s zcU_a`Ti-sYO+xaph)tQxg|XtHzt`i+Y2ckzAuhUsxI`>VgzFJ0!ACuH1AJ}JgUom4 zZI45=^eBL(;Z>b4W)F9+W%Z!Wx3QiDB+a=qbF&dd zVIc4UC(g0JGtD8-UcmlF4jHKQC@?f6i9n<=u%y0v1MQs;!zXhCIIP}rB^7%(p z{}Io-t}hk7ZkwBr+<#m2f8<=y6i_3EM2IJ+*vWrX%Kr)o#P#ng<-Y&|l*gmmKkD|s z0XnN>nZB*^|0|$>G>f{dRIoZ~nqDr8E~d8%Z-zH~%i87~&ilL8_!sv?-@WHC>n3Ts znd5YKd)-R10Ss4q)bRd4HYM+!MXC$s<1}&T(BCbACj^d#YIvpRZo4$#&%UtD4?@s@ zVfjBjc!Da0l=AAz9lmT6i+|F|r0rH4cw1myi3-cRGWbz{e6U$UCkLb#60)s!@9Ai% z9Cp?8{9^*%DM}c|n8xA1{&r+tb6w5-$sDf~2yi7PrP(iCfO!UQLy9T-GmrG3x1_1$ zp#@Nt+u*j&-(f&LDSy?gz4VWl5QDGE6`Q+$Ce-^lEx$Yb#UQyTScmPOr)*m?z248( zb4kTAmA$Y0;Q^(3h{~r4GguBKS&T!uo0ipui?U4hfz>4Q6q&$@u&Ja6CzaO1&8@inML~BI#duklU&nZ?RQkb%&lBD+{#Cl;RkCt zeb0-){=-)>cw)C>TF$F(K5se8#d28zD@q{KfSQvj^S9I(kQsWaQkr$E5wl{OmS;@>;txs`j_Atko_SlrP4`;6AOFw zI}?em5JqK@V+SGYWR;!iSSYlU{bz4@QMkRRKb&^=-e0kJmBNr00^3HsusfA3-idbE?TPF78vcd#OPn}+#sdQkUq0e3zo{*7j6E=amZ;wu*E16aSa?^x zO(V$*n~wn5xZfAi#uA00FeU)L>kG#Oq@#$XK-a5MP|T_((3LoUjA;3|;l6WF+vwyB za6}Poet(D{?5s;kW?tXDc{SnUjFy-g7)pF60UWj_hPPaid|6EIEl0bCl8aL74IT$| zR})Q`hU<+yTM@{8eUuN2hp2*)!eS@u7A(vw%6dbnl2SMOhyXxqQ+O>m`BEGa?c&)irs1@ zzTnr#?u7IL<77PVHSCHp4M6l2G#DOw8Nm5(JzCB}{+T8DXJXP?U-%ot#mla_zgmnV zsak!IYlC7Zf3?`4c~Q9moWJ31;qfE0@$AH(-d{?UExsp~=?X#rV@5ua@zpqrH%!&O zqp1jp?SJCUnTi3x1Av_FCu=wbu)2y65*A@a@57>70e>MKKpNt?Y)fTKlqz6~-fN%l zzkLK(dg>?wAHq~lfAYOjYAJ-8dA4n_MbW^!*AX@l(#Ei(%Iz-rni0D2bp(iN!iUNJs`Q1Svr1{`q^IdbVwK*_fKXA9XGy!c}^-ZBJd6 z*;vzD&<~o^W|YoR#*w%VtwI&+JLR(G^hk6OzyZlj2b(qDkfTw##d^8T>&0_*KRXlW zoTmli7FGreaE<}2tb{(6-*sg`GVVmQA3nGKQn+t*?fw@rnrPKN6o*j6!1Uu zs<}OmmG5ytUZ{6wZCYDWZpudm#B{_Z3e#fl1#`U3HaP> zj#v#o>N-$RMUNi|>MGb;*5^$DD$j)pnbnP#uuat|vPjfGAU&x~A-eD!*lup$TzsKP(%Q`ymzfQk-~uTa@A zAm5QhR9Vs76am*%`7)%UuC`A18S6MsP2u!>$f-!Z3OlSbYu=F>r9eqMU5mmGrwL>( zFt$8WyriMGj>+x9p!Ih>TBdDp4Sfl*aRL%Ag0&M&&&urs@*U3L)F2+fpkzbF^5 zTiauSgrGN?{3ji@WBVS$K4FPM^~KCnzaM(=g7bz@(m#G zKNE47WVkF($NQsChSM`59e=1q@7;|Wv|@*R1wveFPh|K}m1~Zl(hFCLern=J%>lpE zHeCe^M0SN@@u)th&?C!~-a0>ggRX{aTV2|kOS`W1q8%+9?hFy;69H~1;_$o~OYo;l zFX(GoZ%WeU(2KiJW$|d*v|qQ{e1vAGKuxiEWb5fTJZdsBy#6?%DpL8=qe(*|J-W3a zb)@H_@6LJ>aL|eTi68z3l-rxqG=faYI(5jC2c|8iT z3y{%PkvROb`+rd>K|hY7?r0wKQO9Jw-%GE z*O^u(h2fB8MP`uM>v?z@34boKVQ9)d_vGNK2y`vwsb9guTDSM1nq=lt1MfFJ3NEL| z%R;}@CF@n|TqH_4HzXTU!}w#=7;HUNg8_G{g3zrqM4J9W&YUQ8rg2@9t%r<~!KH^E zd~q`)S;M8S3R%GQRy8NQriYOuXK=HA4shwak~PxTYs=tz&2u-6*j3V7bo1O4+jE`m zvP&anXVjtmEqxo&mUVgTc$|VLEi)of;JOyG3YA8Acywn_%y&bo!m;tNhqaCq6PJa) z>~vZU$y?CWB@qPAc=Lx_!v_x?+!3`U3~& zYx5b%U;S!sxq@4R_(N}dQp>`9AH13RRbu~-LG&nczMy^28%%FV;qo30Gb--;H%q-P z6!&FHyBUc8IS;CGQg+B&5xIBn{$qsig(cqQ3_vUaA@z|n+r}Fr{J*ifzlA4hO$4g{ zKLObG2?YPs>i%!+{Qs>dP;4Hv>8y z{hNWbq20SZA8ral#bYrrAR!>T=0}geyi(pf$wewDLbTI1l3;cf z4uC(F7p5QpRUZ2MqTLI7j`)I>VY{%2Os}m-Pbz{2Y`J&*CV1C@xGr~ZMiCHId8a0Q ziMyT0sQPhkjr-Q?Ck9{xtPP)m#C6Ry4FjU`&&kOKf#$e2Ef{cn55PB`v!XhNN^*V6 zEXsNh#D~BpCyMEBSHFwdoK_r`2hLDZKW|ZUF#B)baQ?9M`qUsj0)rHAsn~KBr;hNa z>RleQOIH8vjXE01GbQgW>$0Cl$Y|HTjsYE&I*Huu#Vw!QxyNWYpyu{GE`qLwOF{1E z8gOCN6F1nKL?ZiC)F2WyN@t)qr%`|%ZPrM`SchMHug=peDi6wAVd$;4_lunr@e(1i zH{;m77~=XK_ELb)`OR;iXvp;f_U(2giWRJHe#iIit8uSKld{2``#9pjr#;SFTb{Q- zuHh6n?&GR|VR7{ap=8`zcG=)&I(q^=q?lq?)G+l_a?$qTVdlvt1ryI~2hfvzR3CQ0 zSuVs=MnwJ2!d{)F^pl^}Rz@#%qjtJ-I-1pYg5^h!2OgWOU#?rt6?J*(YIg)AeMv54 zQ8zUGLJ+?cI42{&&KXgiurj(IyCeGI{!Vazf=(-}1T(i!tVX3U-0r2Vz3)eV4%dUD zm&Wp~osZd6eLDG>*&Uv+G7v^@<^`F)V8EDH_#XJo?oiPa_xvvX*JE zn+~NqFWk@%7|Ikbthiwqc}-j4ZH1f==vDWQ-2w^v(1C>tywI$|Rln6MhznShoC!CG z_{a{wc4;VMY;mLgi@$Pcw^ylNM!jj;4lclvufxx22e$Ba!Li(%2>|MUX;!2yMSPZ# z1}VL&{#KWh|7eex|M2nT?$_PSw53_nRG)^FFlZ=`P>2=jWO$;EzP5IRAQb`fO(*by zzDK(&l_>f?B-zP@d_G;1QCpCN=yK~d3IZt|NhrkPEgi*f4J){PYIgf|qA=nbt5U!g z6sYdLOtw`1p?e|ih8wgwRvVol08{%A3Q|YdSgmx7tonO*f|X^=U|toghxbz;iX-x% zz+02Dgf;ETq;{w=I6&D1ehNU~R=Yn(wA9gVyHn{Jwb#=F`I-my8ZM6wcdp*KP<#^9 zC~!dKEGB_t&5kO!!a1cd0{E#;X+NB{u7E&disksP6!&+!e(6drfM*K{^_id@0~o1s zC_X`&8esV7RB(20OZB6KzMR+d`1Mm#4}m^2{*&`fJ`6Q zl7Dznon3>*6ef7HOtSZUp$>P(wZA_HQVN=^!)P@}OA&eeP)GTj5d6x(W*tu^^UN`GpX3aGW2;}F!^M*<&KyUB#FYfuWw1hDFVY_?~XviqbC0sF7 zAWHYZ;Wl@>m-OOGQ=gPo#zBoapZm!qsByvU)55SqS+o@(5Si^p?x03p`TTLOCs*g5 z7FO7%iWX``)coY)0V(Ycdc(RX1f7()k9WP?@McP4ON2tPJ3wNUyx^tO&kKisoS1$3 zDD<^#nKcbwLzMujpFUL>j>@YDYd`nMZYH8R?$KP4<2wgXbe>=**e`#cex_mA3+mYR z;(h60M0*h@s8J|^u>@SOP4w-(J3SQN^qpjr&~VZF4ZuxuX_dn=pW0tK1;70qtG)F8 zoMn>Q(}Y@3w6YM}(i6aW&!iLux=+JpB_7H^>#-{DSUdbgd>E|1dG{O^WTg z4U^J>jOIFDfLv1woTiap!T2Y3%if630DE4nr+h7v0(~D+N7r_OoBKNyjxfvJ)0jx1 zO=3Z})fFHgth>Fa?$Q{orola(B`9j`qtgp84feP4j-Hb1W=AB~#`#F?Qa6?!fJWNS zuJf{d0ETvb`WdRR1T|iO#U_NczlgV+>M6v^eXX26A_sCcXqFpU@2h_~3^mTfJgX&s z79qi@SydpA(pPfGufiUII=}SCpH0WpfdvPB?SDQ7WNT$QY!MZX^&|q@c?ZH;tFu`b z24&cn<1IHrAC%+&O0Zg%zUK+IgmHI;0IYPDg7`RC7;(3w)(QGdUBp?zL?o`msZs!B z@i1`lWT?-yq_mBJz-|)BCcNeSIg_msQOXH+F4c6+8rfQzG)PhU@%X5+k>emx(<25H zunDL{?=$2&oXn4|*3Kq~Dp)r$hWT@N7O#J@LC4)j3v*;#zt`5(MH2!SdRc zSCL0WHs($NE4Q4z&$A^=prgXX-B-|YV#EWt?tEg8^FAdRSp-=795Mxg&ij_j23P^} zbPBJ9J6{?LHJsvnqkg#&whRKro$NLTj(p+`Ji;Q4N6zUfZdP@90=hr?Xy|c(|Km@U zCNW?GLQl&Ix@Uds<-YER@qvDM0M4D)6W!O{M$pzSz_7@DO)1KEgvFMQA#p03Cd^L_3s2JZJKRKk3wSmOgzm)5g=XDTH0|=kox(_FCiKpx1-MH zYRv}p00FR&^EeV)u?WC?YUT#j*;e?bnZjwLRHvk$)ZuI0vvh1m$=h_hcsqP3j1#0J zK3OcL_i7$$Oo=DVi0v6n5EQFO zELi?^TcvqQ=cIp9BV2DwFOJXifHE9~_-F?{()CDDQK_{FxdJm`3+k3vttcJK!FDe* z4rNyoI6$wPZV1WY&fFu@1y}UYhSOyU>+JPOb`bBz<1s#vlGCQyQP(uDo9XKSl^Czt8Aba!jgY;PX!?vo1} zlqJr&``t}%O|;Kv^@w+hI#oq*fPRRod+d5i*Ty67(myx><~z^U%Ev*xh~);`WP^`Q zP3E8rx?>G+Z#sI4=&j3NzmN)5-(I7DQv3)CL+Ke4j_)K|UXERF`r*QH>)v6{m)TrD zxb@q({e;h+c$mCv%?6ZxcqCi(aTtHDFgHZ_!)4&Nf%Nar*TP4XMl_<$1OdNQ6&MUe zm!fuj*|h}>fUl_!xW@$UFS~$v+2(!#SXjZ{ysh{5k{|s6wpa%A`k>_}wsGn=ln)^Q zZ~%%Q=MQe+*udVMD|`;y``~Lg06cEq5$1oDo9rw*Z8HQ$ZnOb<$QmWUD`2Igt#dq zRUBMH(p|tNLL?)O1~@qUWN3Sd$LT|MGI`{q7_Y^AI z1n<1W0xsRnceIj?c!8jB`y58q*80>hDMyrBrz*me?eLf$QXmr|aCUG7$W%LBRcuI3 zV6m&-YfFc|AKsS{a3kD9C)q_t7u?LIFy@r+BTr$LvOlmq?O-s3yK( zlR`G1No0-_GempP9VSU0WI#i`U4867Es-WeTX9Pop-XSg@L`d@+ED~ka(mHeS5gHE zqd|DjRtsQ8=*&ZX#*#&*Z7s}tr<%_uayv0w)YMCng=am*#Q$vHYKOl%^@~1 z{t1Kes@9mVqKSiqtoFAAGn`d}+RHhjj5E8F#dy?0WXiJ-Nod?!$!ynls_X6=HP5;M zj~a}+$<(b*sjX-;kZu$SyO%_dt^jUTg*tR8RuAxN7HbT!l3sqdoFzQS6c+-G+l5-H z!mQ6$bhY$2XoxkmC*iVb>^BP-?Jf+igW$KHAV1x@@svH@btGD~r`@refM`!*s_R9~ z-Y;$r!z8S?(9qMA(v=Ao_Xli{SR0WXi6iilhU;kPGb6(@j%b14LcME!AE_6qIWu-~ zUx1+4iIAfKD+;vq_wwVNxQJ9KT1dl&fsjb&O6wh=(GI#yKzPY&W*B`+r+uZ2Fhm%g z?-40m^4Mg@_94y4_JEz6>Msi|`04_r7+SGQk{FDBnj%El^^(@dJB@Z??u3$-s-;+w z^mxAm!PnL@YAfP^a8;)En=_qDTw4jsSAz(?LfYg6RC)k>Ab;@~5&61j3P|cR&Dpz_ z<%75pqVX0X|LLMQlyXqTBWT?V5^{_xUOR?AT~$J`>myl}J4IgG7>*XWya7Z#+-Fuk zBKWU$sor@lBr6)XDi{XAl~S6gaBWgy2P6o(-AS0yA#>6ed;=s9W{qrWXKij?sY9#A zx1LH8(*`m_vnMzAyA~Oe*MZyzucfoKwtQFoRQ9VKR8fhjuSKi_*?lH!$OAL}Q=tK5 zjt@0sg{%?tz2Z=twU~Z-c=$m`1zer(xK?hnWai#BI@(l4Fzh4sw!}TF2lC={rTV;i z>gLCQaM6o#LJLpYs)Rz7FD|D+&%WHDtxofMllcaAsfBo_wuj*?c|Q!eFfFm#7#IhF7a;BT+w2GA5z_~55i#R8hZ164zRDv z>{|D5YfHA1QuW0?WuDs{6NWG@0J%#xi7VDFPRPQq-6PcZx&h000Y20MM6pza{d~m6 z{$bU+O!}V0|W>PeBE=RD(4Cga`h9^=4MU}m?s|g zRLEa$qMQ6&7jJh1Leo#i4A6~{CAaHJINWlPdy^b;7?-BXWAJ(9$$ALjRvZ1?LJK{^ zj-)!jfxfxi03;%J9FUYqZx%si}=iuF=eO8sONd+nRo!e1)oj z6;??JI_a-VMCyd5U#b~_Dh*ViDvYYvp4Y}f1P(%d621gh@rIYt8M7xm2UL_DjJ`a^ zEjb}e*Q(T*l#0TLO$U`B-5$)CoPlvn|61QYIE2ww@_PQ8mRfMcBg3u&MqFyy&vnD@ z&`g$n%*)}?s?3Q}e!KJp0joH1=uDN$4=_2FH9|b#=Pw=A#vItU+%WA2Bur<{1xkmS z>s0sgK5=10c%0u#OpFt{kMHb5R<5k~A%jdJ)mJ^Jr9CV1dflTB4ZOu2I|mol<34Z1 z;Gc0Z();Q^RXaD35#?QDlT69wU8sfdhBgdM%?; zLIlXv1F#B~<#~6fzC)i#t9!~*zyxy#s2bSg^sng>!>=bnR_&2b?$~9h6s}~IllnHZ zDG5yGbJ6>9(;P~h1D%DD`il&BCTTLjU(4K2gI1QQk=(<_Jm{UU0hl1Op}L%zRcf>< z#)|vpot=6eKXnP{ZbIn+zo7ip4XYcthd)fws`gdI!!7F`8tlki;03iMI+?H1XPHon{jlNR<|HF|v?0K(hSFW)3 zz7XKPQz8bxt2~2ndPKVM0;28Ei}x$^8d+RX$BA|K9#X2RwknBul?3>0nMJHe8_HF4 z;Ak28ds$*I+$q%y8{=`rTdoX7a{BrJcZHPZ#08p~0P;#-^>Zyv70C_BKQ7EK(!P*w zw_4kEe!bor|0V298{84xJoh-lFR#S%Cd)f&q{*C0uo?w zd%j6paqvg~37$s@-%2TAwvzjzY*Tvy2w*)*F}+G`$YbZvXJem)G_vncr@!3WVsujN zX}=(FkW4U4R!+ErUSF=14c|)#yGJ-&t{2VgmkLxE(E%yJC*9+XH9T8#+2yh`VSA4H zx&&>~!i40evc6YW`?Swc<)eR^RO-!GrILgOWA(xTTv&Hp4U3p5}XA6O}S>|7$VH=-p%05WSjQ zKF~EGjtMwbXY~5{c>65(1@$G#eMzuArP4Rwae@r1lkechX>h3FnUW*gEi(Jd;+F+Q z6Ln0o^UFb?LBC>|;i*Dmj|0|c&Ng3YLgMM2}<2(SF(Nl)UcT`AY8 zI>$COL(T370?=0{CgT~Rh2a_!HsGWJ-XxolL;XwpMMIjOo@@BI3nXSzKK*t5NZn_- z;g@6311}}#^XFwbmLLnZrelyGkn1cEF35xpUmUZ0ua4=N7?cuab;d5u&??opv$7AHjEvc+R-btcpl^Cly@wV`b)x_HFcy3=t%!bji>`rE zDxIY)9+0%vTb!)|d{HU1)Qs|7;q~~QvX^ylJM6cJlEEn?OUoZWo>aQqV0}Ej-0=$h zMSnVYyi0@&2sHA~G?><#NR|xA8~)}nceWppc-xsNz7zE3?dx$=^=j#?5tgK-0(VLM zsC@~11b1x(1VSjLsm50x>^z?2f80bew%YcD!pN|~J+Ikqw!vk!25VWLvbO6;q7_q=d|~VzQ<2n)TiDNc+_e@6+OKN6#bREw3^DY zrYgWi;vSv|I-q!YJAv53A9s$FUnl*eroJvO<9Zvc!a=O z(U95?`p)}#f==(Cs%IhR#WoKKNh%g}hsO;@nQ;H`foYUI=KvuetNWDKebe;jG?%EOX8Ykc}c{ zJ@l#38!QEdA^7tx62)E&x|_5$z3|~Q?K>G;wmBsJX)Flbg=^Pykt(KPhUAw8{SFzY~<9FjE-fOvViz zKe@^1!+x+NOf2X(ufzkKTFCoX&{BAoig>43q(66a#?s@au_mx&K$=6CBB*imGdzwU ze0B3==)ilRVr5f1RAG`Ct`R+-&zjBH7&AZP19ElrHUCIdFHqnB0eNG`L5G3F)UC(R zU4I^dw6$FzUg@yAwcMbsyf#oYklf0>-TB`h|NFSDkzj8&WXLdrh0RR5#Xyh{*DR%M%s-qXvT6 z{RlR*UE0EHfe|~;SE?|q^7kZ98y=BdjuF& zK@s3@NK_QklDlCIj;2KAA-XvjUP$>K#)b=t7$DN~Rr`dBDZo_2OAMU z72Vl3uzuU)#!9PEWG%0oOx%Yajd!U)6%VmP&!t08&(}Bn3=vteP-m%d0F&+KYZEcT zo7$BOS79S?K{yBIG9iNXFEEI&@8>dGX5RRmG(UhD3m83+gi9GoKM+6M9iB)T^YyLB z3W&f)H;cRxu{{o9IBWdQwwN*zIH6Bf+b z6pJ95G}C56PO0%c6UJ)O9QofbX`~Y41vLo!pyvFy-;%34(-vg4um1JRbmu^`#9eo2c=fcnYYqH6SavzYfe_ zg{!jb`g{naG=k60(SoYvChxcA@H|qGcR2j9wZ`d|jj{ro09ui>u9kiYsDvpcuREF4 zcgoML(6;J(A=_g1e~=M