From 09ac878bccf9787096dac14ed9f29ffe781f9f4e Mon Sep 17 00:00:00 2001 From: wyq Date: Sun, 31 Jul 2022 10:42:00 +0800 Subject: [PATCH] update curve_fit function - still is not usable --- meteoinfo-lab/milconfig.xml | 40 ++-- meteoinfo-lab/pylib/milab_debug.py | 1 + .../numeric/optimize/__init__$py.class | Bin 3150 -> 3176 bytes .../mipylib/numeric/optimize/__init__.py | 2 +- .../numeric/optimize/_lsq/__init__$py.class | Bin 2630 -> 2606 bytes .../mipylib/numeric/optimize/_lsq/__init__.py | 2 +- .../mipylib/numeric/optimize/minpack$py.class | Bin 19335 -> 18930 bytes .../pylib/mipylib/numeric/optimize/minpack.py | 63 +++++- .../pylib/mipylib/plotlib/_axes3dgl$py.class | Bin 92428 -> 83771 bytes .../pylib/mipylib/plotlib/_axes3dgl.py | 198 +----------------- .../MyParametricUnivariateFunction.java | 57 +++++ .../meteoinfo/math/optimize/OptimizeUtil.java | 30 +++ 12 files changed, 167 insertions(+), 226 deletions(-) create mode 100644 meteoinfo-math/src/main/java/org/meteoinfo/math/optimize/MyParametricUnivariateFunction.java diff --git a/meteoinfo-lab/milconfig.xml b/meteoinfo-lab/milconfig.xml index 5adb67b3..9f2de2b5 100644 --- a/meteoinfo-lab/milconfig.xml +++ b/meteoinfo-lab/milconfig.xml @@ -1,34 +1,34 @@ - - - - - - - - - - - + - - + + + + + + + + + + + + - - - + + + - - - + + + @@ -36,5 +36,5 @@
- + diff --git a/meteoinfo-lab/pylib/milab_debug.py b/meteoinfo-lab/pylib/milab_debug.py index abe944a9..306adcad 100644 --- a/meteoinfo-lab/pylib/milab_debug.py +++ b/meteoinfo-lab/pylib/milab_debug.py @@ -6,5 +6,6 @@ sys.path.append(toolbox_path) sys.path.append(os.path.join(toolbox_path, 'miml_dev')) sys.path.append(os.path.join(toolbox_path, 'EMIPS')) sys.path.append(os.path.join(toolbox_path, 'IMEP')) +sys.path.append(os.path.join(toolbox_path, 'meteoview3d')) mipylib.migl.mifolder = 'D:/MyProgram/Distribution/Java/MeteoInfo/MeteoInfo' \ No newline at end of file diff --git a/meteoinfo-lab/pylib/mipylib/numeric/optimize/__init__$py.class b/meteoinfo-lab/pylib/mipylib/numeric/optimize/__init__$py.class index dbc8761ac1b8f8b915859c209318e6c37804a75f..ac561d35a0e7e5b4ef6f9adb543bb1cabf02dd64 100644 GIT binary patch delta 1229 zcmZWp-E$LF6#w1bP4;H94rNUxrNv^Qkv6HJSWvJQsR*P(Xh8}EagCd$4dkoI)>8Z^ z;`cXr)v<~nJUF~LrHtEXN5?ltAN&J+^2Hf_&_BU?Zju=uoZZ=T?)jbbJHLDG?*1D6 zGivSz~d*`C%7xZ6MnIvI5P)PUXJf12YS@e|wKy$PEQXy_)kY_7s_FQrVn zUNKW}P4vhz4ZUKslHv*UNsjx(jFRZtihjA+#=xo^0}Mvh%T&DTY|fkF*g;{N<3TCI ztjpLX{!oVbARadG5QaojJ#jRFMum?DwE0boHtb}CY@5nO^(j#IL@G@OGTGqy%#W=oGVoPOjs$fv2!e zT-El6pE2-7DIq8U?r~k;Fb~BYG{cEq;EBRo%|IR zwGL2|gS)gg?bhnz?Zi=Pf27%7`4;0#cyf{SEBKtA{ZLjhf`Gmgjt))7XF59Q2GbE) z#jp(B2I^6KL3Ld;B^rk05q}x&^K@X1`OpV?2ULp3dY57M-h`{FKVa$*sXeQz>JqIP zQ`fejhT=_>Z&0S>p<`~=w~&>*+7ilAClEpvUK8C&+!Bqfz?raSmNDyNu7-j6-WpzR zkcIFHUTvmZ^<+9i%c-1g6g{uIXZlU2mC#}01PK$6M3{-h0D@plB=X@Pj0h%>Ahd1LNk%6`Y1N)+!cOQEvPJEGRd+(w*_X&Q!6$>gvAl-E;4Gr_X!c zcM`uQjGO=d_%nbF`0ARWvab312gEHUK65Op>Sz!?+a7-qD-0}0BZF=q*<`o3ALU4jo2;vO z6;>P2@DQ;L=87E6R5I-Po-^v2cv!Y+SS!{m8J6QymA%_O6dJud!Ky7@LdVPGpd#8&n2!4x`WcbAl7dj;RNO>AFe>}IHu zId9r67Hm6w$$k?%uv2#IVlZtx=etGEoebNaS!~OQ+iLU9XR%l6-52sY`Mgaco}eR}A8L)da(!>=+gQL^5$=sFPY9 z_$CS%lcU!}Pc+FU#dx%rPvdn1Q+PxC7Tv+$LRoUWEt0y;PKvTVz)$0hfm1ju{?JEN z2j@gTPqGVQnjeTvuKM<`zr`*7HTz6V#CpmO$!7Qv3b<5FXDlJLP>fP!lxO%HVD*U5 z1$>{@0&vSnYZcU{bsFD*)i;m&J}VJmWj1X!JrZCI%Z6CmY6-BeH2~$>vsbYzZEb8I zf;CliW0_#C`YC)&ja^?k*B$@g?=GgicRJ zxAwm`lY6!35`CZdL8Ype)I5%+u46z|f5Ol{(m$fA>Kx4(QkL1iNgh8x%OuJ|n>I_E5}_GnjPwcgpQswMKmw#q_abESW#Uk95mzWS XQsGC$IFrH05OS~wpWstDDBAKnvDZt$M{bRf|#!m1;%3*2+u!fVP@2B5~;S zd<2Ox2qn-)gUlTG4-EVOKS6OH%FxN)YwgF{j~~7jpZW9e*Y5!0c)w;Q6p>NF$skUc z@Z%(d=DLMa$#r#{CNk}o3iZYu14~)xMV|o`{S0=sHWQs+Y|K_G(PFhW70oPah=^?^ zqz~eZfdNFtZ)K=E&cOTna~6U)Z{i#-2#1dyylmkjE*Wreg~2ak?2!1(ZFWs;^A-CB zZkp78%L2nySu-q_)YM_>2_j|^zi;6V?n*Q&cGTnSq4=Zb^f8Q^$RH~wJ;VA0a#G{5 zSoXMVQfS_y{v0n%JcBE0-m=$0ks%;{c>-L+lu)#B&xcfC`KZWipIB9>`b5m3Z5gJa zK$pwXyZmw?@?P?c;tUJKs>D6yv*M#3TT-_W$aC9lqbJ+KsdU5+4Q}EHOGlIqgqt|K zN~H`nvL;DL3~;h^+rD5N2PguB!Yw4i>losE3)bjB#JYH*k)-8|^QP6zJNr6AkKCP{BL$ X$Ea|I2F#`~hc^<8;4R7$)KU2dyRvjG delta 739 zcmZuvOHUI~6#nkqna*@l35?Y=jYbfy&=zZ}*4KzH9u}nXkZMI8OKkxyls;5&!;iQN zH*VY-!!tso!KDim7ybZ$f{8IECN9KtX_s#9cb?xl=gVF8e)5{1|GxhSAcloCGp;n5 zg-XdO%*^F7Wdm)ZTM75J~yba5_{EtYeU;YA%?VnqoU5u9-2IHKa4Vz-=P;2oVQ584qoaT;etkdHK7@ZcQI zyJ6!ZgHJ@+Uh#@s?6O$l%hpv~GpWDF1BOeorcV^r#6Ie2$4Qg?ArG!2A+4L@i+YHS zi0^937{e_SqZk*1T90uXcO*wzywn_aU;NgxhJ%MD9v~woU9&D5IR?M@s`;o3>?`h5|z_OX~ULLWC|tMsbD;dC$q$2q(k~Bf6w*!=L7s zSwrhsTD1;K+(IB3vO2mpv7aSF$~%PD5nH7ehAO=#QAr!%^wNzD;u!W&1PBDTa5lJs zZqB#i848E2E1dIn@)_eT&8VSw1O01M2`n)wO@pABP|~XF=$GRFHwG}+$hPi|1H2z# zr>z@XNN!@dhE%YIyLB=*CUCD&ZS6AA?Ur4`WSxy6i${%op1y?&h}Z4|rTi+ZU-y(o z#-Z$xa;grS=@QBy$?dwiLp%IGLz+CJZUyrsFmR9xXK2*p1fHNEjRDM}D2)o{{{c!l Be;_*O>P3LN+={Dz;1*@AQ;Gl_~IFbrAVx}ZjueG*=*9yZVXz! zwJ0D+wR)|fZRG>4w$&;S4N+8(P>_cLJ}RZQwYIiZYi+40nn3!SxyfclY5Papd>rP? z%sIdFJ7>{BQ=&~b z_#+YjGKIcP%wE?Ns|f`Zx|P@hV|seFF?(ozzcdE43H25wto3AbqDpKekVoqNq5DG*hxrLKey4SS$9x+GXvO?>PO~Z1k z1s9($7!5>xK9znzZh`qjdj&EOf5cMigP~BXKNJJKO~fqzP)MPNgi_O@(2t0zrh7Qh z78MYWnh-xWhP$TB{3*8Q$HLKIW8kt#q$#4%<0{>6l6iud?yC>B1?qgwO~G)~w=@uH z4!~rKz3j4OcfK=WT?ukE^LeWML07(wwi|o1M~-?1qy%AqG}shYXor)YWph>PFkxRH zmSr5ZWE#UqC4-j`5_X^l;F zu~48)p##PPBc6yKqL-blfnHHj2YCiZP{7)yf!YSa^HogzktJb&BTDl+CefvV`dGMD zpDppm3(3q>S1WZUMQL;F3na!fnRgBlWlgoMid>QFsbC)ilMzbu9`X1idj3 zUD{Nq@KopxH_F=a!Jz0{5{QafbD&n?>8N!p5)qZgv(R2zQD-v%B`24lC0qnVf;!}8m04-AZcuboWUI*6IUpTs8%ZJIG8X9SfQlXDg-6@0prm-9;F-keL7OZjE!J$||Ic8-r`85Oy;${aq| z$!GH_7;!&Qe|W8GPZ@K9$Jgl8(8*?OXZx$%8I=Vc>^!P+r?Z_MOxebo zs&=qN)$XU*wFTSgZjwei840VdFI(N$Z?z#|wfD4H;UA&vPW;fzfbJJ`cET~gj2Mc&QXo2y>R9D<> zmuoA;nohUjv^&^}0-xgZ7Ie_+Dz}z_R#_>#^LDg~TPGM?@4-6!Uu) z&;mSnh6iDFT4&r;dVc(_4t6{#0!cLUw08ju7rYN^8nZ~OTwv`QVlnTXVQzzY@;- zjd^92`T=39NSk%{({QXGG#)9dOg$vlwG?b*(njNC*%*4+$S%Kz3}ab&K5aHOm1omo zV^_J0ju^+w=g|hkQ;|P;v)5&F;UG^?iLJ!u)>IcxdL<^nJt`3Us`s!s<^IVd#+r(; zvtE?}2O%J=#hJNmTj2HDz_J}V3UJ`ru_LlXYy70d>ay-Ty$~y}V@0ROjGl_|@i&B@ zPCqkyB-JcJC)+9B2(uAf%yybhbg&BiJ58qi8I5N>n$HB5&vg|~IUyWH0W;^KFYVq_{ zTXYi7oVVKQ&iD|_z123Io|7wfBvUJvw$t-{=e6MrhkiyR+*jVT4r6TPP^H(b!B%5N z<*YcaP~s+-eVQRDb%Rz9uR}3Dt;~tLGfrtE<*p4y10TeDPsaOD z#`}|u_kq_VxS6h>Sgrk*e+RryA89!f!A|!64Y+|afIu0}3E#KQ@{PlJv~+e*B0eI+ z84z#=1aeD;c>gR2AQVajdkeJL@*hOe>D# z{#2IWGZ~>MEMI=SI+j!4C!bGuAJkynH} z3j%=!8#v?LDTtp7}8~fY8g7e z05e_Z*&v!xsrfBdTxBgSiRYSu76K=$&c7xy;8Y(&b8@(XtL3 z(pi^jtR)%Wr`ff%{z8~H5AOr|1;w-Le2gr_I9Uh}$eOjmTfr~%#>HbH2;s5d6EL6< zASMb?ghAE$L=36&J;6C_O=o+}_~}XUpG5q`jr!N|Ma@A8G~5}p(%G@NiI$OswsR}q zPKH+#KUMl3d<`$D_~t`tkBV;(#!`4O$2XdIS=mUut71262kWZV`NiA$)DB+S$uAk# z$*(ZMQuqu$v(G$38sjrEqpOo^SDJXluj2EBfw68z_(JcN4*qq>)sSk)Hz124K8PPu z3kg7$K!RJ0Q!_4zlfWVIO@v<@nZ4|Kd=ca4z4>9xEQxRAUh(Pj-`|PsL{T|1XOARe{bC=ZbTmD?X6cpf@24`CA_fDn0v5CtVl1eO9p=yH=Ruw=8lZgwH4 zbPI@26!b!Qzd?kGMW`W4)hcT1vrw$IXl-jtKWl5P)cWGp-^|@?Rw$o;+4^xZGiT2E zoyVMWa^oqJ?pi~xSI#~43}b8{-?7b=U!>3R$Eu9FaBaW{jg5!uVu5f-p+$~SPp%XH zttU3w9Gz2WPR}2fSnOgIw1k+wDP}}t3f)0WiG}AG(WtRdp*x8=>cjEcpkJZ8A=DVI z_Xica7i@^bj6%zabqE@_E*vx82*wS|)C1r(8qMvjCWTfK%LZ>;%pa9fU_C@EVXut` zg7pfmg57W=7HACI3R^$K3Fs5Tfe=irCFYEUM_VqfS1G_=$Ye)oBZPeLXJF6LtKd-6 zEpa0X-7Um2eD#6(Ds5xUBD(D=ZMHOaA_S?S(8I)BzNo)2=#NHwmJU(r5ldtbY=JM> z3jN5mcP{Gmn6Tpuz;T~XrTyd)ppTp7o%4&IAZ9ay$jBkFsYOgr5>qYVkiR)56b@So zKQ$M3E*tj@t|!Gqh|@nN8V$q7&s5rHnSP#_?wcKG_SgF&B5>bae=y>Q$P140F$*7f zf57r#2u%d6&@T@1fy4P9&>3tO zTHA9l5sMn3rrF_Wqr%~`Dzz6xjYz~_uW*h!i5O9%u}R?>urkl6Q@Ddz7B;?=lmg#_ zwF-BkK6xA)6^zILg?9jcuoiC=Lb?^M;TUf}y5U`;(bS;uT!l^$qh^KY!?~uwoJK?8 z1!xMCNH`de`F%-VS9mCsnZmpb9@d8Ap?Y}O4fY+4{@C0wT!Ly z@gsy@Nd0`oB@)g1z$Ok23l!cLFIw-PZN!5}IrqS)kitvA6!FST^8QebM8jxtWULei z*@bZ&20=U#RXCc9M@khw6lXYcH5@VOVsRs=@ZpfJ_csNiBDo_Ge{(Dh_9(E?>!KzV zK3Zg{$?r4j>M%R_H9+wxe2f4csPOA>j3|5@K$Bg^Z-6#}8*cm>c-ri%1GL6@`_*6dKJkjG?Wd@&#TIJs}dl&hr*1WT*PO0a97q8>9 z%}Dg~F;t#T1CQ9=ib5HL@eK8dAJEOChd~(GYk}KwhSkXOrBWP_tcvGKK zJh4fqhE{g{A@*{$$6Z;{!n#nkCx@+SVafp(u0G7B)p&lwW|SPDT_lBc+MSYGd$H8s zHmOHbQX4LoI?^VU^TrpSH5qB|f{)tJBMHvBc-_dDHwk zuG4EN&#RI!{0X-606V`Qp%UL>mU>3#ufl#UXai^;f=98)@bB?F&wX3LoW`8z*$!rR zvLjBG;pFS>PQH>Jay$1(uFCN z+0L-r?5A1dX?8uXZywNxHk(TaOw8MYTL-*jtS5e=WQC9BUkB8r-Soi?N?6`ZEnSgz zGY~f@X&Kk)(y}yTjWFgv%_?AQ>>yz*G-zC!Vi^=s_Y(58LEF=AhBFdESrH`tWzE>XSC*f> znK9c2W=HJIHs~O~`gtP7d}&zEgqG970y$D8o#U6)5x=ZYJx(7AW!j|BlPWQ=c&5e6`zF(g0LOenL zMSg|6*QpK$q{4%Okz@M+zR9zB!5Lh?@B~ikDJOvfCofMqxkjRDy5$8#b0HeC0$M3# zb%MvUig6(WS(gOW37sY3VH;|BDGLt?EiDhy)gI9iC88fV--@3dnC%=pg9_hd|B_7@ z#l?G7GR<5$d}v|?_A5aT;ngt3#H-;n8E=z;bp*x&;!!*{;CkG}a7m}F()SF(5$qOr z0={@e23u$cPHFbhk3g|%d^Id!IwmcQU^ZbP!Fs|3qDz_(sDSXum3CYLY8{|3UM(ow zc9zMUqLhHNQ6LFrBlmoPHAFE3mvH5&9gUqG-m-6mO**5`O1h&eTxuQ zq|owpbQ$vmi^6-?&NA<0O37F2Z zpv2=A*!KhHk#@EN;TE9gVZnhL5ilC%>0ucK89zTe1)M$L==7`^AJspxSLo^VoK;}i zGSnrl?Ad4WR#wFnF#B0H7}F5{&eD7_91BjfkH26A@Mf16eDnDjMsw=Fj zbW5v8Vfh&7IOt6b#v2GfMGx2~ux!*QvtQfq5dF5_LHeD|_B2z+7sJ@^Z8qBhFoe@fK@f{R<0t|HcC5=U6gz`qEqMalfi{ zkZ2hajRYEd$*D+N2>P_QI76qi-eQD%-dpU@ISF)b6S`k$-KDM_DS~SdKbsKpLI(2C zg=uIBKP}|XGBT2)+ztfbXhjmjL5?i3j@5POhF_OqVFX0JX~T`Ym+QlH4W=Vd&;3s zR{eFBA)^CeewsGZQprwJnWANBo!SZx&==G=VnK~}kwJV%b`UR=9mFq}9mIRc4ic%= zTsrTGB}?bMu_RR4d0Mv4ufVZ`&ii5MsPkehn$EApk|SF~bL)HnPP)zqV#(F{U@Un$ z9|E`Cm6~D8K$B}7(!Ac(UQqB5Pg^E4KTO8C8=ufNTSDk`Ow-xbp2!hWr(8-un#0yw zY0=Vz_i7HULt9$V0!Ryn391LFM15Dtv{cHpTnlHxf*%+xd@L3`c$qgLCJVnFN8|;L ziM&){QFVSJmTWnE@G39nunm^|oHYAKV85-ulB4HR_QeapKHTuQt8&;r3(TEJLp%OB zrmFZV$KQL&X|G|%iZ5QI_N4eiVJws1#3v+gyJO;t+q}c{@qr7bH7&Z zv%oU>EN&#l-BK7I815sjd`=RM_*@wcKKnp>)ffj)h zpv9o2pk<)DK=**|1Kkf=v7ZD!i4Q9LHJ5vy&BEUe_+Ks7pAXEI`b*-x__ArcN$~ZT zy855oB;+4A|34QJ-;|(lISYx@OMHul@X;i`dpKXq#kbF8%*Ms%^we?Obn{hWSH)NJ KHDc4m*Zv1+O9&kR diff --git a/meteoinfo-lab/pylib/mipylib/numeric/optimize/minpack.py b/meteoinfo-lab/pylib/mipylib/numeric/optimize/minpack.py index 2ac45dc6..6dd73ad3 100644 --- a/meteoinfo-lab/pylib/mipylib/numeric/optimize/minpack.py +++ b/meteoinfo-lab/pylib/mipylib/numeric/optimize/minpack.py @@ -132,19 +132,62 @@ def curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, n = p0.size func = UniFunc(f) - jac_func = OptimizeUtil.getJacobianFunction(func, xdata.asarray(), func.order, 5, 0.1) - problem = LeastSquaresBuilder().start(p0.tojarray('double')). \ - model(jac_func). \ - target(ydata.tojarray('double')). \ - lazyEvaluation(False). \ - maxEvaluations(1000). \ - maxIterations(1000). \ - build() - optimum = LevenbergMarquardtOptimizer().optimize(problem) - r = tuple(optimum.getPoint().toArray()) + best = OptimizeUtil.curveFit(func, xdata.asarray(), ydata.asarray(), 5, 0.1, p0.tojarray('double')) + r = tuple(best) return r +# def curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, +# check_finite=True, bounds=(-np.inf, np.inf), method=None, +# jac=None, **kwargs): +# """ +# Use non-linear least squares to fit a function, f, to data. +# +# Assumes ``ydata = f(xdata, *params) + eps`` +# +# :param f: callable +# The model function, f(x, ...). It must take the independent +# variable as the first argument and the parameters to fit as +# separate remaining arguments. +# :param xdata: array_like or object +# The independent variable where the data is measured. +# Should usually be an M-length sequence or an (k,M)-shaped array for +# functions with k predictors, but can actually be any object. +# :param ydata: array_like +# The dependent data, a length M array - nominally ``f(xdata, ...)``. +# :param p0: array_like, optional +# Initial guess for the parameters (length N). If None, then the +# initial values will all be 1 (if the number of parameters for the +# function can be determined using introspection, otherwise a +# ValueError is raised). +# :return: +# """ +# if p0 is None: +# # determine number of parameters by inspecting the function +# from ..lib._util import getargspec_no_self as _getargspec +# args, varargs, varkw, defaults = _getargspec(f) +# if len(args) < 2: +# raise ValueError("Unable to determine number of fit parameters.") +# n = len(args) - 1 +# p0 = np.ones(n) +# else: +# p0 = np.atleast_1d(p0) +# n = p0.size +# +# func = UniFunc(f) +# jac_func = OptimizeUtil.getJacobianFunction(func, xdata.asarray(), func.order, 5, 0.1) +# problem = LeastSquaresBuilder().start(p0.tojarray('double')). \ +# model(jac_func). \ +# target(ydata.tojarray('double')). \ +# lazyEvaluation(False). \ +# maxEvaluations(1000). \ +# maxIterations(1000). \ +# build() +# optimum = LevenbergMarquardtOptimizer().optimize(problem) +# r = tuple(optimum.getPoint().toArray()) +# +# return r + def _del2(p0, p1, d): return p0 - np.square(p1 - p0) / d diff --git a/meteoinfo-lab/pylib/mipylib/plotlib/_axes3dgl$py.class b/meteoinfo-lab/pylib/mipylib/plotlib/_axes3dgl$py.class index 4673b5d3981433de28772a5aa540bbf12117f400..6981fce867ed08689eb24fe0095083070fb57d40 100644 GIT binary patch delta 20850 zcmc(H349bq_J7xO)ufY5!c01M0!d6blW?N~g4{O&4EGf#7YT+;NW$qcTpp;1*lK&N z3L@St5DZa3L0JTpOV@SPMGg-d)7y?Ryk-m6#D z9q;>qZMmO??7j5Hn~I`zWM?;oJb6RBLtRd$1e3|?14IcaEAdq3RFo9w6=>>%RUtic zY~3FwioIf{$5UCPsU}ezd7jcr&+PJSN%WXQeMmNcB2edJ@RzDJH+pOxo*_z@uTi0| z)oO?Ov~0D;Qt4~9wyIyB;bi2CMA2pzR1{4wFV5H0mx&ToF}u9b%_V!4KNM6BN3Ev* z6&&!;KYx0`;L^&1@-j``;80(X?5l`EITKBNt?F3cdx*D&C@e4EUFp`;H;59NS6<+* zEEqh~JsoVfqF{PK<-l@xSy6FbmZrW55m1m_CEv?m~}T zP)*%UlqPaMrQB(qx{#{x`VS{f-*c#YCG|c`nw3UV4^&kSe3N)T0Gk2~NmKuh!IoAQ zK~PM?Kd?r{6(0Xg{1eZ4p0fFa3&$20%+=IS9qJ*;`Cmkd5S#<#X`QZfsGrN~BlW@s zsx|c(m&RRQS)5n(kHHSL%a;&h&dlOcP5oNr-19W`8=|N=(+eExcgjCQlsnY#A;DO8 z$?Ssu<>j7oO+97A49_5EjDEwj+WnyZ7^19Je1$dOZ=jY_PD=Ldir)#toQCwhp({T6vzOroeeCQ!ghSz1tmh8fZT`dEtEI%v3@kM4V>hl6Xo+1kh z&>UGl4Rce5_62#B#h%i%&eu9Ht-=iGu>r5 zg`SdpKCro%&jQ$4k@OyPzJqq1A%(dGYLYz*H^PQj*ymNd)-(9t0;s)Zq*WPQ+!g}q`$uFMcpeuw0Q}}(hprG6Xu^)uk zIqv*?2YpD^!TDPoi1~)6D7*%V2?CrPHpajJ6Bt^9_vsZq;H1MZS*bU*{N-! zchI-7+0l0}up4H2@@JP6^wj9zj4!6H^S(=WbH;ld3X$yI)94Ni!d*TcbF&|E-aJ?4 z5AtDSAJi!K(?#4hQ}}6CR$Q3}O@NoSFQ|YSV9wy7)%;?#|CFmosaqHjoHLoH(a(|0 z*XU93C@CoPg57bHEYK(%FPmGO&$SAk*PiRCtn|QL!S`x;^RXU6i@pQVJh7G0lgQ4O z+0)1_klC}`zs|1o_&EKDn!-F!?hK9ojBJ6-{tBX8cX@7kf!BRA#5_k7TS=})&!eyD zawkOT1zw8u5=6~|Dr-!E3n{M1;~r0AHi$i2V-(38joFdR)tH7hr3*9$muxGYuQ3Nt z=V>e)XtYc}eVA)FXFgUbyTel-W{r7^hMP#!tJ$%p`fabL@KSARh(2;uFk zG5A(6^vvNKn#SN+vCP;gSX<7}Q<{@sP=e5bwdarS%KFESe2XsQ`%{I+I`h@^lp=R6 z&$)RS+zEA;@LV^j16Or+12V<=n2w&j#Dh@|?ggrOe2Zj#L=KxG>xbNYQ8qy2c-bK2 z7KpMTBFD>yaqW`E!m^Apxmo&_5HgvKaj=nW6wi$#O1#l4FEV(T#|`xfViSxHbH{3v z*yIp4kxel=@ah+`HHUWo{4rrp&$PpuVNshYd%T9l#bNbBD0G@LbGf_95K1 zEc^)WSLQx8_7*0QPmJ#hZ(@r_?bOB@y{G?$d|@1!9vlA^QS?c}X1M3L+m*OWr?(rC z%eUh0LF|O#DvBcC8ZC>)g>6#^`&7A3fw!e^8}}An^x{Z?ttyru5=4qPT&6=0^GkNQ@*;GC67V@{~E*3PoY-&x=B}G3eZ;9)6&?+v>VlTb@lZC=^j7notKmD*VXp}qzC<^Mnh6BAJ}1C z{V+iLiJ$g^hP0gYXamyYe$r!?ky5*^eih()sh*T2TS>iQYH6rWgZ;EL%un0dN-Ie_ z=+qe??c^uD<8snII=wbP+RsN?O$V-r!Vz{7zsBct7IbQ-DC$DsprX)N;7-5*J^+qE zX_A4?z#UK?xNf05_s(4qOimg0Mcv8Csy>T!LcTY~E?yFT3QXZ?E+|B zfq#5!6O?(&IOV!!;90e6*JUqKN@-03SQlM#m*p9pME5nd}%Z> zGlRcmLCPIamyF%NuEVfd5SbjL&>#kS_KneIb5fIW zbxKQ1xzS;5RLqxP7R0_&lsYBHS~q2F95iq4+7Zyas&}97+eCEns?zLgU1UPUt1(@= zxmASdkPo1FWUVqPEs~wAdULav_tNwA1AL_6J2Ei%?mOPWMAYq=iu}~K`9#cl`$0z# zALn)@(QthgMZPwgA83mn7afR-^$#PmPRX^71C?E^m7^?`eCS+~Iwh@9DW4E7*?|M; zq|Nu}GZ$9JkPd+F&VVbBZhah{VD?(NTCG&+>fL5pPm*Gcd^_wqztO1~u>4Gg@PLCa zy`qGkl)OTwpmh(yi~D`ftEX8OK~mXU$lhqKOd!$bx(E_sQ+x4;RyMV_C^QZp(n%k4 zLpn*YseKLdK{DxY&eBP|O&wtL{2+l0G!C9fAcM?;ExE}U@x(QD(hp)}8%GZ(lF`PZzo*5H z;Ujfw7;09!GZ^FECAPA$#-YzU8fX6Al#DmRK1wAMj1Ksnh|&>cvXOFR0GVPG9myb5 z%?@@FZM=LWl}s}Z9BE73#-$^zNv_fAXiJi3j6B+w}9%5K5*$c&s%k zGP3rx(A8qzKSTEsuT?8G8M^wYG38WT*bH7A5niplrK^t{%Z{xjGmZAg+mKRY((zQ{ zF}%o^8801gA2zE2<+0O`lL}+_tFgudw+FE_47SpE=1eM?W1M*_**JQ;rn=O*#u=2% zGvdBXAoGoOXPY8FYJVhIV4ii7tUa5JtAz&6t?xK&$b@( z+1BI2ww@5S^`x+^zX;oUO4!!Z!nRfk+gcs4t!IU8tr52MoY0@kY-_FGwq6jn^`fw? zmxOJ-B5dnbAhhY-_u)tsM<0xoy4Ez_xY@+uCKZt#^fO?H0DRN7&YGVOx6} z+tzzxZW`Fu`_i`dN!xlK1ee*?{wv$or2lu@S`N*-A3(sbJ^-u)9s~>{@#}Q&)5xs? zRs#r&)HT3!z*^u1;6>mi;1%FiU>&d?*Z^z<5DuxEfY*R(U^7qyybinpYy~ z4(tFB2B|xNU5b*WzKhImV2>}|ixiQN`aZDFm+ptr;Iu1^Cb$naB)g($*T76Z#azl0 z12>w_jCDV@u&MvzwviNS_q)5p^~MrxB>fZ|u6{07f4=a+uNc?<)WfD8F-FddrQmVY zSou?uO+D^2qc4OReJRZ7E5rSGJf}Hf_I3~uePi_hxwTFG&M5mi%BFrVEaoJ)7~|;A z)P^@0d^b7`i%Cu_LbEeMGkjJPknXJM4kzsm&o9{ZeiU}|ldzkgW%=P?JhOjEByp#RLn+VT~Av7zfVvC837?drRUFj==2aB&lJb^BeCvW6|#e zJ7SZ>4vt?0E9?0lnjUcuqpDWk!(!7WqE?8_WyX$3oyPg!le@SA-PTt`SU7DIoSZ@| za1vV^)rA!sp7SA*uw3*C?s9KT1>Y!R`T3-F(S9ab>9*I%_)OS=Vq>T9+fj+UKq6SJ zGOJb@8hwHAsMh%Ge6sf>D*vH$fgXPfd{H`{hbU`aa|`!6Ptj>dQR?c}ps+J#$fGP&01f2j?* z-YB~ijdxXx4rh>HW5l^=WACMq*zTempXr36Q>oZD7$+|!ksgMlE{XJ%J0e%N8;#L* zV@WUL<+{Y6-a2lesgU3KxK@53RaCF%EYM{*4rV7gk6xQ9!L9nmhgLbF#g z3AWL}4?jbaTMXf?tb5uPgxSz_7;zK7Bty-MDv5uDkk&!NL@!Nvk7T$xl#mudBScY% z;XWQ~K0pW~Bh6KrBsGqXk`HuRlfD$A8EuB!NGcg)_GnA4vC*-- zu9;bBBRo>Re?hFV{sGFmf^dSYIbGV+2@!eR1zIKjaNqa zCN;ue>o(Tk+RSb%ndcaZ zB^9R2PLfC^cX(aNY-7W_VnSUL!Euflh&&RkGat2+Az||x6rHh?Ib=b7fYO!go{Zpi zF@nDdgK)yh!DZEZv|ELkVlwN!;BDrcK`?;Z%_Bji4e=V2y-EBaK^K~BHIy$BKG5`N zBqVaNXdNXu!ntxMxWw4LKhDV6$KvQxQOi4T=n2i88c7o;_^-`V8tKy)j*ARYZi2T+ z1=4`#KntKH&gmezzAR@Fbd#3 zjRp|*Q-u8#VLwIKPbUFWfE-{NfRG<|*LVX33%B`_m-ZwaE!=m|Gid#^;Kgq%Ij?8OxCjLBP5 zZS-w3BZef~=sRXs3>l=;osC(4Y<5WIGxImNLFoz}2CH41%;&&9G0(=3rI~cM=p8Xv zm2Sh<;&V@*1>BQAl%#viRk0+}8>kCZ?d4qfrg>#1frS_67apWunJ}<>@NX8b^Yo=zcK>C)kPzlRK0H=Ds-6%tjCLqmTELq9|>ZunVLK z!l;)okg*7%YL$Vt$}n54GT^*Yh2IOxU{QLwmc-U7Q^cKhwUP-Nc^Y^ISOu&Go(0wb z&jD+J7l0Rmmw=anSAbUm9B1fXf%WE8yq5Cf^a5K@C{8bQ3$PW~2H*rk-vV|3Zv*cD zJAuCeyMW!m9$+u<9`HV}57-YJ01m>0@fz+5?tmXfx zV1_q&l37qM=@0A$n9&~vGya}|`-7kSR|xcHt4sLB=MsLkxP;&29RfXvI|La0d7shO z3ZuUujsBv~=r74z1X>3}?xHhse6ZU533FElNsSTM}#<6SQt^o?A;!* zN{l(FJxPw>Z@#SaNNv12+HBvN*n{H4uw6XjGq$gxsVrU;a68i}TszK8P8kFE;sTZ+ zt+Gi6?DC0%SjG)0t7uzE%S8=wLxV0Iu# zF~`MYl(XVV?C59|!~n5C91ssA0Es|TpqWB&l(~X=u?`rae47S~U>*5*rFC7pf^{(j z9jqgG;gz_wuHdimE(#)PIMy&kn&8D>w{c#Bxpm>%`4!}caC1o@L;2Ar#m0u2J4HI& zJR#C7Go%x?1~$^{(23WNGDnN_CUdSxv(3jlku)0{ZEhFo7#w^O98SOQW@Cj_vhnhC zgAXqDUySe0Ci2nn?VI_x??97D^_BMe%0S&@Zoo|=M9^;0J6vn*?@Ur+*i_!uTJbH$ zoc>5Ab{78*%yKx9e`CIa8H=$056n2XaSy_>sLO9V{ymxHSxFhcM?Y^q)df2z%jbba ztALdiSSdkBBe~cSMcm$KwCT8(bV-QK&crKnn!w!+w^9GeiWPFtl4|ROmh#H$Y|JD5 z%z6|~Vr60?xe@BLskB&eFri6!karzymW3uo(uAXdueO}m%AkRtz5N@pxvCwBG>3I1 zJx%ZRq?68OH=d#9=Gh<;o65uq#IoYFC8lnRvSb3ojy^$+q^?K0lB94p*V>br*_D)% zB=gCxq`l4-H0IsfJkpiKYD}I~qU&8$XC_Z3{PdvHbU7HfKrt9`8N+UC%qr6y(T%V8 z?MAZ=;{8Vp99x)|lQ|`yc9Ks_13n9-tJ~iV0SH?psNqb}k!X{>gac1En*(2o5U@hZTmy3d3QA;jqGn z16cqLD-4GfhQkWW21YApSSktejz=Dc6*du=1WX2|08@b+U>cAM%G3V47rU=~mg%m(HFa{-)X*e$@Vz-@}sliiMt7gz`^0u}>H;12O! zkS}ui23JWGbFvlFe>bK~ep0~Q;WBgN4LE?ad-z0h+tKOu=E56Dw3po*DC#DbN!$ha zboahMB{oKW#WzyeYfn7*NR|gkdP$OVAQ6@P8g3(*xE>z1m88uk&ti@@{xR<9`RIq#52{(6ZqC8U8f^;dJ}Y?Ke*U~=vhAVAwGnsh9AMu z2}FJ>L-OGO#HTV0d>T--#kV(8bvhYL1#QG`Kk3>=q@Z*gy~~oL*h`n!@i-g^3@{LXp|SF*!qZg0}v`=UcUu{<|=Y(Zn=O-n@e54IN z6lw;QQrHHd3fKwYl*2v*4)gDNFlYD*P)@RK6<5EOxGBE@Q^b$+I+Y8=_Zyw`Kv(rA zY4)1gw+~6Rv1+rZ4^Gc)vw2S+($vOk%oj!ay184_zhRz2eK6Ysc0!$_`;x91tnvAs zmV&uX%kV3<`#=cvCRg}oG2D;AlVJCpEW8eOQjEURT;G?(q_FLm)nK3GE(#kILW+6j zH+@Od?hKcJjT!i@23|^RH?7do+nmm7Tno%z{YZ+?cwJ!^ltOGWOD`zDiQQ$GO=y2N zB^Cf$O-_`68r(A|pF2gC8?L!X>{XJOi@DT6eYIb-CO?_W(L0eRKG471_VU|7?!rj50 zFt@3H@=0sBI%0D97Q-jk_r|I>NI$uZw9#_K4fzL2heCjIK0l`2pHxT+PdN_d|q*?Nje#86I zRm62~Wkk~D2!+KyTFvhW`AdL+@_x%HHTgD#8s^G@B-LCvm^i}Ocf6)y#Wqp#edCHY z-bcXGe+$P>S_mHH1C6AQL4-*91Q&#F0M#T3E&T5j*(tQ(0Wt~$rve;(>jki_mTY)X zZygUwsD&%De(wQKl@xbyK(t^XA!`-M~bCveJO z(k`6++PJKRxqUDR4rjkLu1NKY8@6iNaBdN<_?i0VcCDBcAxU6v14o(P(c-eXJs)6G zKev`85(6^l4>_kyGZpp`>}# zm-XM{2?Hop`04zF>FP+D=PB$5Je&q_FTijoz;GwPa3{ci1I_`z1LuKS0QUmyqVM?< zEC!#1ON${tO6NducWlQXVKIDng2mYF#uG0jwyu#^N2EVWaK9;|Pyfk{Vfzv<@$8!E z9Yzw&XF3vQK0S=YbPK5`_8*b`cbj%7on#bLe|!$upWR`e9fq5yFf(j8$w`C3HuMwg zg`W^NTp~zy*_)UT3@4F2U7*zIMo3>nt08R!Z|Qq6*64*YqLm9|J`=NgtrCV!m6_wa zlO!+CVn?ty0bD?YLhxrvTjX51vO$Vl6xYBuo9Xs8RdlC{4=^0VGq6YQQdM;$3J6hM zPNjzKt)lPO*GAmsY$0m#CzAf!!FFeiXsqgVUv>u7=>hyQogOlmW?`rJ(A=9vhI)50 z(S~r2B!^)uLVx}y&t_x0-0hqzxZ5q1#S|B-e#&&6l)IG0f-1;DRoa;97MF6TpwfI) zUfwMx%EAnPjM*C9>SVD?xnB@FtvxSn%=2EC@`xm}cD<}I&rHGdL6`EBB$g6*<-i|k z%=8JDvR2XsTgLNPBc}8t!E}vFStp4@EX1oD^L*8%R7UlND=h@E`?cUh>-lu<aNF0o z(|@{@Ba)BPLiP^}8?JDqgXr{#OZi$7yDY?q8}t0erJRyv5yZ>YSgy?1jhQwVOwYKK zpCvUuR`(O1YRvSUOR1G)kruLF8}X!{37&N>Vw1$pqb$@H?N0nqXquq5~9;HNJ4e`0}_W$e?k(b(_fH;M-<_=iEitM z#HrKY%z|vv&P#vCW2{bVk;LirqLB0+7wM%^u3NDdN$<8hlN$;O$HM1BjIoXj$_N+K zW`dbrXDX5)UEP92(-}h&tg|2_Avz025~{O-p-3G%3qum7vnEKwb&9!zq!g3KCB+mt zT{?3iiO_8ukmx$okwoe&3P}{CHz$t9d6-@_myRZ_u{3cOxohptl>gsyUsdcj5bEER zd+Zn~_rfu_rHw(uB;CfBGFexn`C6#)E*qg@E#d=XQe*8-y}_7RT;sJ+llWSwDPV@R zNI-(MfDvIW#A;(Ln&Js-0YbhODST<5G#AbIvE({TVuGc^ad!J3^amNIxokC3UlJ`% zrrMp+4chRH$geM0V<b3z$pf4@W(_=|| z4A)oahU!ETp{o%{pc^q}(m0&>cRHI|q`9xr00P2b0ESCVmXqAf(qyIG8TFqJz}a{J zX8ZAw>QTrD-R6QjrQndC3myT%h5TauBkR}yDf|_e6fMaj+LQmOXxw06jE&~@KNIbF zmx@i9yNP6rXwTW5&i~P@@df6~?zM_0)*HIn5{3R{w~fIQ+yue|?k4ygHmRSdYL|LY z3X@_HX46%55PAn4geF188cwkclgZu9S1 zRwjglOW=ZEHYh$)7(X9P7`~H060WO=NG!@K9nvQ*YH*ogUuhPh4qt5?$6fRTIgaL* zCPy1Pv4Blto5QkP)j@dG6Ob-gSNXOWSe%CA_*IVMycmbrm#|{eX;1})a2j${mt6Fi zl(0oWLTbOtpu9B9#e80)rKL@%-5K|X`pF6ImAGc9E9Mr7(37kf!^Yil9P z_m63%i}`~6cB0-Z2ld8P2F1OB98{GYRC^2OTdy=IJ|-OgxKHh1q4hQr7N;)eJs{~i zTId!x<|zcA%UtY3N!-aoe3#SiY|=neSqHvBvrasqV_mqWs!75*9)haiFJXrxfu>%E zM1%fyMG95zjs*7C0|}htjYzOTMU%I}Ul>y~IB`^n=2P{}`=H-&DULs$E;T;WL z!(3OVoW%FtWHOO|RjH`{Kd`KWY=1(H_kXLSKj;VlQ-}Yxv0h08@1HTLSL)~w`tklr z>-z7C;KkqW{x?PN>EC}%1n-}1VE>ne@P7H{+xXQ)@bW+2{x6E){jXZTLMQ(Zoh|-! k2$$pdJ{zLJ9;!E+nB^4zYkJ2&}Rq zqM`ySDkU1Qiz0$3&4yU9U;|XJ5`_OVGn>210sg<@@AoMWGrRNd%sb_I=Y41QczzMP zV?PUf^Yn8sDT>mbo!%IBU$4kQ6%K7NQ4%W4Jv9Z@1EaF#U+(hB}Iib zCB+VHDN&;P+Sg^CijjDvtYV@=y9>{k6i+O1Xv>`1LNYd#jM45PN|?VUhjwr6E8Pmz zt}BV6S5NVHY9=|f2Z-V<@>JA#rd18+%{jDHPVGMV_{07(p^8IWT|2s`M{V~gQ6lVz zO6{i}cWP_pQ%?jG+K)X|d$RWsGHfGJ98*fFCrzv>Lq}?f5>h>_s^` z4Aj6Wg?uQs;l{+0ngLaXm6OVf@*UcCsDK;8Clywrc!#{Qh2O*#nzgUqa%wNj7I&jX zuBb!XOOz-%o>K3u4xLEtgZ+k(?(YyKDX=7X*4}e!`(a*Osm6dl29?N?Rhsd8w4n)3(ekb3i(!TH|5w9{eAL)salgcg=pPAV+m zby1@Ce0Dki?_U3X)J&P)^UDt6Kisb@tg0z1Dle&a&`_cfB;kBk%4U~1Xaq(S^ES*= zR#7ut+=+5(KZ@>9H&NVTAW%W=+#!e5-Z4Z;X}BY+3>q04=cI;wI)NyWjb$A)sdnzL zCbVr69n=;7L5~Vg)f5M93O)0YD=ZgEWTK#`ysXken-OKA`EyvRdGv}Q-WC`o+ES_< zAl(|PxIosEgGGz}^zoE?s)iKe`8JrU>7MdwQ_$JA;D)ENu&AtN7Q~DBfV78}4Y$fm zrk9lCUPshXRXD4lsIa<*n+JLU1_G5{KBRC~A9?SJS)c&D2=g_es%+w_DQp`VJgEZ0PWQ z1w-=-`V2njyy<0uw)*6L}MB~=bO zK@RJ%D$nGB6{Q{rEq2mE*>I`({K%^ZlsV}nd3_B;3ul&8*IDe`^=2B#9w zmlV~&^k#Luz)34*$y6sDBXiZ}*sCt`PH&_I2AP#iomN=xptF!I6UK>Ni+P*jsVXln zsPR;K+Dxv*G+mD&s4RqGrd~{X$;`6q8cfQKxH}ok+d*%_Gq7I;Qwl2!N;-wRO@(^8TIXkWtf8mOAMY zd5v8qT-04)E4{4t(@`6Vce&nc`kaR4Ur<~&-APvp13>Rb2WOX5d7$=#P`jY87~}kq zEW?K4C@-t1#`?v?uhw4~A?wh4*%>@gI;|YTv<|ftl!<4bl+QlJM@Fo=4SKKUrS^lb z*cm6iM=D&0VJr7klow)oZZSKK*&KrngYPDNs|siIn^{v*fr1y!Gh@b*m(8(bbJDib zSNRfp%}Jk=HN5Vm(`4=qqQsi*^ceGnvCUGjJ7_Iijc_VU>04oR7kwN0RF{;O5+!x8 z`EzMvN>=5pu1>m#?#0ea_o1iVrg(~{m6!Bz(09!CR zt0rQwVE*-M-~jMEtUoQqrtF}w`*afD5$G2fk;<}~A}kGP?&EvqpX-OGuUr^f~Tg&0|WLOlrEZu z2@=cd4+zZ^`zk$&>@1l*h3ssZJK1Bi76;4E%Hp5>|iRgB{EAPG@-C+LRE>k z@S17|(}|*%PjD~?+L|c0W!w$pcab`wY7y49gGIuqlvNjT8|Yvzs6EZWqLEB@Fj#Oh z!@=V4OvP*mgOOG%W;s|APiHz<3eutq2TO%mk%OfpDRZz4;#Jg%6CErIx|J8s0&!S# zNR~np)|gG;Tz4>-VfY80DJ3;ks1lagQC2;zqHHR5N!A9+iJmhU;ZR{|3ou^}1_K<< z-56LM;jJ)keYr2@yKw=O@b+;q7y?x2k=rhVt;IBB>tNk@g`SFn;*xUsF{}r_EUam` z?8QBUO1`~eh%V-{>8U`jFV7Y7d;O6sEa$mPupIc>XM>O_E5>jP=65{k<-mMM&E#7y z8zyqtWZC7&%@X&n5IKHt1ah;*y{kly-y6+Wmt+qcV?I)t-*#sh8N-U4Y#b}#xe_=) zo~ns?m9uIlc`EXXJXIxmeXu_b_7q|@g|LZc#)J`$GB!DkO=8!Wb0+k4Okow^1M`?~ zO(-DMW=>I&V;Y+t#%kCM(<<_KuZ6K!R9^Y@PUT6N^Eye1%-!fDlV$EEC#jaXTb$%N znVaV%H_6;>PI8CL-R>l&%q?^hyXHHcWVO7v#7Ul#InzltBj-iKcRR@|@}}h^J7sRU zlkAnb6;AS@%&l~i&t>ibC;3k1RyoN}GWW2P{3&y*ovI>pYn^JY%su8*FOs=+PW5(~ zd(x?%#(L#_U+>hCWp1NW8!dCSPOVVpo`J6~bDQC|%N%^=wKDfSyk?nu5#FrKy$nxO z=C;Czl)2a7|H<6z@SJ4s4KLgvncWEoK<0K~ua>#DvEj+w9&DpBw+|br%)NsxQ|8{o z7AbQdz^=;NhfX?E=03uPBy*o&3zE6du=(&@{1|o!W?kO=66RUvzBUs})5tew=hDmB zjUV6q<|uP<=}Pjm*=Axw(g~s%R}Y?CIK42hys%z1-mI4!O9ZNM0qlpce$b6xkXvHnM}yfHnbJFysi9fw^mfA zI5sJ3@}E~89hLpG@$ z4U`s(=8Uy z!RfO@5eZ&7rz<&GyJU_*$e-rw8*hrNN2MW(dP>P>N$jL~!Q3MZ6lt$eg`0P5AK7mr38Rj)N2R>Vi(x*|1y+u5`yV0}Hqx1!oV!sik zt#1i@_BE7lLn-zkQM&M!KA_{iTjDV)-!+na2a@kX5*w3{9DJ(>k}uqv49Nqx25S2y zB)@_rwksjoao!+Edgt*vt$BfxKST03B(aGJ$rD1d$9yh1Vt%0H87v6=;ai)K+%%up z`4^Y0jlFFk#Z(t8%*3sm7S6-N;ft<7e*7X{Jaf^u5EOS`mQIeFGUkTD?&4PlUHznQb@(xOJ5QCbAtH|MO5 zjX#dcLf9XQaz-f#EPG;gB39n^)kC3lZQR;EQKE?_mDIm$A{Q$DU8V2TwTu=G@(whM zsaGz`ieV>fBQ|?^D+{;m;XNJuN)9@|;FZ@f5cyljA^+r7J`k^OUFrtyczE}2|lU8a&GRlCqU^D2^yOxOMt z(%H(7At@nR7v8pRw%eA(E<`I`P0#)gq?`HFYZf&VI`vIaJj@(=D$Y#b0o^Y*7j4Kg@Ax8&3^!N4--=veZvCPKxzdcceV(*T1-wi^XHvKR;MVEVa|Wk z9UkQq9c#`xoL=toU5^tQx}ag5^67HrG7L(AxD~~3sipaAm6q8qH}?Uyqij4c(@*znWDDsN;~W-IY2JZ2FL^20{nh^paXy% zM!NvG5V#2F40Hjy0^NY_Ko6iN&>OfI=mYcx`T_lcOMrojSJ4I`gY8R$pQGgiLxEwy z<-l;@3gAj$1TYe~3K#{92CfF+hH2w~0${wI7Ak~$;@l=GO85~;x-d?lxS!vS zjud9=A#>H2Yse}y|3K@IhxmP$%iMTFf-Y$Ju(|R;Yx0P><3QVp)r}rH^VPj%t@+@G zndDLPrLQx|<7SVS)6Kja9GY8OXGR>(Bu|*_4kwc*%^`=IAV2kB40#HSYaO#*G!hfY zvkm4;5Z!2gelR6&1B%b0&L%+}><6trN`|TL`Ayb2OntC552nd6IR1(@6{rHLff`^M zFddiy%mQWu*8uO3@Z0vjkWQ z+-1J|O}43h8#8Jxt{w#*2jC8APXJHa=~GD818^I)jdr>T#6Vzr4q|M@U`h9%6XVrj zV!XzQvCSsN>ozgA3u3%s6Jw_!#+x=Vc1dEqC5Z91O^n@q!u-V8BZ#qA5@Vkr#(qJJ zcO)_1vx)J(B*q7V7zcdB_)rkzpdiM_f*79&Vtnc+#^-)w9O4^|s(s;S#aBL7d@WdU zSg_(7A1l5UtoTl_;(H$}eh{oUB3W^icj#oMgmzu6|f>nr3AD{j;9I)f_XvJ ze~uShvDjV(;6-cOfY*WTz#G6$;7wo`@D}hkup8I|>;?7#`vLgX+IztJzz4tq;6vab z@GP`7b+B6aSue;p1@DqongQZvCelp)n9G6s8$Fg~j{RY31&J&e~T1^;9$q4n0QQJUOc}y-m zLE6sjd^BC9?ai0Aw^C^b^VXwT@f}4CVSG;@Sml=hk5i^A2`MBer!Q{3A5t) zF*BmK&{gh;b>;&<4rqU|cp7(V4Q-2X@NB=9J9w0z0A#y_!d%-d-0Ev4{FL6QpTE(D zg8l(D;q-*Nib5%6*GS(M}D^Y47n;7 z`-7swfy8GJO3Wz|jqKEV<r8-hpa^~;iluvgU7p* z0yE)6Gcw-naiUp7As>wB2EtXI=s_mXhP57@D7yzkkgSs%kX?>0W|Di&o+q;+R)|VLE`;f? zo0TUsBZOmY@Aa>n51wpA?l<3r`~&8(lZoU((@~$2GMe%whHM@ec<`2A17W9}nkwS$!9>+T2&)oUAbm?~ga5PNl(h&col2!v1t7kD06P zNH7<_NZq5M_v1p>21`_TYEZ-zjqXJJH64!hgCAyGK-Y_A4WZ!M(y(AJLipNXShz3U z5SQ9kYNPq_U+H9%@JpJU4kNYZx*O8OQheGRd^(*xBW=;?u$VgWaIC1%wndxGsNIRF za^b-yqHf-JV~b>daXL#JGyP(QpSgG>j1U>5Tn^JP5S zaMkH;zyjb7U?FfPummuHyMViaWq<|T11tyb1y%qnf%}066s0wN5SdlL!@wiJYG4hp z7J!3K9|Il-p0K)vlFSix$ZiIxYr#hw0o%zRriZ@#)h0 zC~9=NM7Im;AiR_}g1wX-0bUB$Wq{EXUdo#dUdk@N`FvZN?|+k zkn~c%XzZnY)!?OkEh?8@O1R%kIm``chTnXD6R5i1e18k8nAV>=7vG6`ggG>{>oFg> z1$F8ck`LqP7yn)e@#C&QhvoC@v;v}CpB4C-Q##e}C&_pwg6r1cWSlV5Uy1el4aYCOIvbN) zQIu>Y>J+)KvAZEotp?%LYEr#2sGf9Ck=1pO5F3^)WJoSd z6hi9MHTl_8mabLSqv~)$O6HVNI*SmjWRZx{lUp)uP_pWfb5!kWD0G;;Uo&zUHq6GM zFZT|zg)fcM>cIV-MPvU?Fw+<4x_Uo!Fnpd*sLEo{_vR{#<(^2g%Hp`Qkg2kGYoA79 zRhD4Qp`^W)L`jOulB_WiB$Fgtb9j=%i}F>LYVD<@pUTp#ct)}Zv-Duaf{&7Y>z6g< zqQZqf*HJRd;6$=oMUr&5NOwW=EOXuSP_CKv6C>#=Yi>pBDH%_WTfQ|AuYTtTPULolZcSPb_>3=a~M#c$I1AsgY5H5 z`UNrdj_i>0R_R|jZ8~vD9@QD@xf=I5Yb=ioHuShaB-vKyl$7NhtF+RJ% z+FMH;BvoY{#5n_wrWK`)5^;(|f&c4Hkr9nftWPSnIJZ8f3>0@isV4~yCZJBq#XBk$ z&;V*BxUqXjA{zy`05=c~7(fgV3m{}+@jwEQ1SA6pSy(EN2BZT`0fa0p3up#32U-9v zfmQ%456b~~^KF1U0Kp7HFf*JXnqeJ4c6`BMlfBX=JN!r18>8VTyNCuGZ^hpRZ^f)n zkVWV#EJ8oqBJ{T{!T@OzM5xDOlXF>wfqW_&0+(}HgiEDG7;IaFA<`n`ON%hH;puZ( zgkcSZ!Xl{bGQsO^{-ts`7vOt*!%Lm-gG_{P?_|Tds4LW5^=BfRh%viDnuRNEvoOLo z3nOi_aFvO7Ng;f5!X|GurQB4Eme%2FX&uIJiB!I~N;jZRT@MP`+uvGc5ZV~WMPm7y zuOn&xLX;HPreu6@2hMLw3PlI_R(MWRG9l2E3_%mcw(TehZlb|h=aM2KoSffqp=LU;uClFc7#D7z_*n@`0hiFyJx( zp$8id=Ym~<%$2|hU?gxAFdDcTKpeux0R_N#pb(e<6aytdDKHTz111C40Oi0GpaSp! zmB3V(6MP0{d-lDgInhh7a-wd=2+6Nhxbtv}&s4~F!t7RIE9Tj@V!mxFZsWE>WeenC zJ-b7)exc3!MUwS*a@MPCF+Z+vLY8noUl6~PPe>q(PjD8ytQ%uUOp2fEX0X_=Ru3o1 z>~WWvPtGz>4E{+ioLub>k}o{QG<4u@E?{Upe#cAiQ>Yx?;})vCj4}(I%Va6Y-D6vV z<>EBKs*EF)+-qBd6;^r#=}J~g18~1I01xn+srZNiJy~Ua9zoiZhZ;zjIa0*g0iB{M z+Kr)z%i%L~@5S#VJQCC~stZwAL|fGN&X}#X9*!g}Bi8UIBw1;gdFHFJ{H1aC=tqMq za$0+%NE{9_*kjiHQ6xS3@nDHaE|Egk@!|yQl|&MrA!2y3w;B{I34C2F!fE!TAe8UG zgFVG7bwzVCXFKU^y;bQVX=H=hdv~0mQ&T1UlyV$&T69YJL~Ob*z!SenZC_75tXKLr z?7k0ycM7{1xCOWsmwqVKCxNGc^}q&&@NXt<-jvzAS%5Kq z4dWm`*y6O=#z%#3!%ci2T+n#FumeB1V6O+0Yr7zq2#?{Si>MgsPD!peZF0RO$@R8P zuHBMcdnCE`N^Kb0mD$Mi%35$n{mQ==sU@wItVJ{>0hj`Zl;C zKe@h>~p_cn5eFK%mPIZ!7hA5YP1R%;~h;kXC zTy_Zf0{9a63iuj0Y!zmcW-SN~^lb)(evOa(7+z}O;9ZB;00!~z7{MStq=6fTsyoDr z*F&vg(WI%WhgsDk4Y!tw)M-5%O`55Cr1iPDj4U(nm(UwP~hIO~W z<;6R%G*!n7uM?;yaYl zn+7Q$_ZsW6SQ4Y^8N$DA>DzAf%wXvZZor}^@$(*ZkxeLM~k zy58E_5l1pqJ=;1W(i~oRp{nPKowp#KWDL{Wh_U4RkwFpj@(m0J54b?EilbouQcM>? zrr5I#{!Wm~L_)?cZOdzrAt-HU?TaUi4ZTC~SadLVFG`QWs9FuhNebNCnBYE=?qD4gI3vIxe&pC6WY(-kBG;;_#JA;)QqzqjwS4 zhUg8-xY}{H+G20ky9U?P)jExuRK1%yWkZ75>vUMAE?+jvN|5K3#5by>2e(1i^+_Z& zR_`fV1S2pu@@QS}#gEp9%l3Had9I^%c+jm7eJfvJ%l`j<=U1ll?RZ9Ne{u&pNp2)6Ym00M3u z!M5HM$N(~dET9?C9B2W+!_Zp;2)OkeAQxx@|-vmRtNdzk>!QM1OW%*!d|eK`b9uPF7gLNL?M3#yNQh60F~so-AVi# z0=+*l0JsG34#Z!BfJ@PnF%3N-kr3oPIh-ARe|GePbqd~1C1e|PoI%H9jq?k)4pq=a ze$BghIj$Q2Jb+Fl@<#-u81xfABG5+#NwvYR63Pua8I=gx=p8=M;vi8-daS8UNP4Ue zD-_&T6)%djUT8wX`viQ=9>+HreHefUOCJRk1C!8M{HVpAwx6UUr@V+>3`foAXOR_# z*Z8V_rPVJLpS|iMEKe$Fnxga9!}4=V;=o`Vv_uR!Hm_Q<=bYI5dW>19LjF>o=I^(~ z!E!VoNPRefuuexKKOm+=qk(V}c(9I8fSsE-FGpb7j!$8I7QIX;`QsH&7VLCP`S0cYlNb%<~A z^-?^FpY+&|cC#Pl7HG2Uy!BfeN%McauTRDU_~j3Oz-paNlD$|YGklALhurd-cZq)@ z8MH_ioV`fqo-MWT+-n3YF}S|PVvXeLOU^F3+o1frXSpb;)m3~2!a%9r zFh*5q&Y-hVy?xq4x+CB44YFQs9dC+LY`w;c&%i;kKF#W!L7J%gbn7aS&ah@=@ba0~ zeJFS7vrrvh(L4msj$UH9Q%UKK{UBiEXAO8balbX*S>b zF%!NLr72 zb1)hBotnMm0|_|~iuehMPY4v$GkWpL@=( zK+F%{u>$93IhWJ;HfEc_`Bu~BBvs6hL&0mc8TdS0Op3I>cVhsS0f@|XMCST(^a#KF zlRd(@SSa6@=)>s3YSAN3L*Ah^yd6V>d-ifp3|#@IR&1Vi^Z{rCPtm_h*K$dR{Mv@D zvU;^3nO0sa;*8WE^%0|qJIBZme=6U2L6zk5xUXB$I!hyyY)@{o^SNf9acm zTHtBm8K4f>3~T|O1)c+*x3BTrG5%{<&JH-Wiyrb>0`wSYN=!us?~f8ujf~n;d-B`Gz$JNtYcBlD>G>cjT(?Y$(H7q>Oa~ z*!HINKn{)$cUdpwkOJ>pTr|)cV8i5*1w15fHMpj4a}nW~Hh3Mlxd%55`U<$dg>J+3 zULl03z(=D&Uk}1HIYnkMF$ota5f~_(W$4?0*Kx?FZ%5`0U}9V`Yr$m8v5J7 z9$+ti&bAaU^TqErqTfL!8a$^rT5qPiTi4S4n)agN9Fl{b@?A~SHsJ;#n%kvpp$BT| z!G_Z4d9LQ7l=@4P=t$iaFP_#6`qA(lYS2%H=P-kQW_4&o`jJD{tTyCQ?|vqpA-qOX zzz`#%J>L@>hhul%;hG_8*8?QWC=HOzHb|9QnJc710;DQ}rLJ`=w+bnTE#>9SVh95& z!!Mg#p;?2>bt{X6uq&|TIl*<_?pBsbvB0Kp3$8O))Vb8H+%JWt2430m#le-Xa4T!1 zTxdXlmIqZzzZ8`|;#QuL!eIfz4+Yoxq+6+zV&Oh9Z!K9DTDl9YD(q-)93 z!Ii%3R<=vAhybw{gX{c8)Vag0?3Tijwy+rg?LpFNUm?B6t-LR#qXMdaJ6MAQZsjv6 z<_ZvdKd4U9PSp9CTRAL+-9p$N|4)J|{no7%N4N5e6voecK@T~n zM}zD9!>!ayu^6A&CAH+Ypi1d?qS7;NqDtvlpL8wxOLs*yTC4O2n^byxGSs`@73UXL zTZe=+mOUzDe;l5p8}w)1irUFdE|9In2Rzup+l>Hs^5p>qY8um zi6q46ML$F8Fw|~HLJeAvB+Q_HAqh9=86-}F5hM`?(~v|)Pr~0(2A^=3%fJUXF&#*r zmxbUm!C;|C5)BqEG+p2(y)+6F*?oWOyH01Xq+raYa=NmEo2fNr<7n zhQwho14*dCVv&RyEFMX?!3HEEg`&wwplA~$kp{)sK~svsgQgTi;Bp%*4N0`2!c4gg zgEd7GW3WslvC!TM?T9Fugd3LW_>fQ}JJ!VR@Q66n~; zT6qCUiWkG^mcxjD>3|O-R!TE#RcB&&=eyNc<;XP&(CQA|75$%Ubxx%;zWSYpMv+7q zngdCsp?SlQMhSvpZSdLTYg*VeZG~G)m#U=(sJ2{p#r|`udC_3FR0_IrrJx(v3A#Oq z#IG7(gcw9VdeL$b8gdcFG%Uh3_@_9ezG(saKJvHvdjG!}KKEGiLHw`vUGLVO_DY4* z0~CHzce(zp0p@GZw;~(#Df>f`>70lU%NFwE&d!rgg?4N&HV zzoiE#0eWBvK@V;fK@V;*p$xZ}pfFbr?BHa9GMp=5aAzb=Lwy`cgjLg-xSDWN4BO9j z#_HobM;r7yB;ZN7hb2efkA?`DDB%pubVz6o6ihr z{($as|69%5#-l6hhSlTUK5KONkfC`#7HFQ67F76*^Tf{usnj7i{Y~nW6`<4C|CLT+ zszMFzd$i}Dv9#gdb18KWosPL#C#h4j0G*BmZv_9Ak=3xwycj6$*C3?=bfPNCutn=a zC$ZCGLGd|s8`OyehF6Rl4Rf<TyKmcGu^(w?+&xW&42H-cX7 z*850Tv{iucFa()i5ranGiGR~c3bYmiG77EHgPhz?{BDMJC9jzj<2{U2BRinxiNSt1 zyf_4#Zhe$&B*!O<@h=Z{=BwO#v6Rd8$<q8uMc0T((eSgJ#X8qd1*71Fco7@)>x<_c}3g*IGhahcEpd9;Ydfd)sS zGHhL=kl^0cNN{tk7~;WhU4+1JG&=>-w$Ow8#UkE! zvGMp4Q5>r%O>wBgOJkoWV@4Tl!e&-n$4a-b^442erK~KJO=VTKxFLk`y9I-c*}`US z;3qr0xmj#Bm$2^bNlscDdXeY6b?o4VI`%Q}DeyV)1@IMccmt6yIPr1-pM70;MmeHT z{6j)S{!2oN=KDv5&hR&F|5=ImpRMlSuFd-&>iuswbe?wbWBC7c=>8EEyzKaYF?;{L z3SRX7f2jigu}Zm||3{;M`Tr+q^dC=z7o-2bRKfdC8k7HE6fpn)s6G34Yn%BGvLF8- mz$sXHmHov2LfS;p*mtn@-w$Lzup>M(h#h4=@=P`R>Hh%D_gG#4 diff --git a/meteoinfo-lab/pylib/mipylib/plotlib/_axes3dgl.py b/meteoinfo-lab/pylib/mipylib/plotlib/_axes3dgl.py index 8eaf21c2..4a7c7d9c 100644 --- a/meteoinfo-lab/pylib/mipylib/plotlib/_axes3dgl.py +++ b/meteoinfo-lab/pylib/mipylib/plotlib/_axes3dgl.py @@ -810,59 +810,7 @@ class Axes3DGL(Axes3D): :param cmap: (*string*) Color map string. :return: """ - warnings.warn("plot_slice is deprecated", DeprecationWarning) - if len(args) <= 3: - x = args[0].dimvalue(2) - y = args[0].dimvalue(1) - z = args[0].dimvalue(0) - data = args[0] - args = args[1:] - else: - x = args[0] - y = args[1] - z = args[2] - data = args[3] - args = args[4:] - if x.ndim == 3: - x = x[0,0] - if y.ndim == 3: - y = y[0,:,0] - if z.ndim == 3: - z = z[:,0,0] - - cmap = plotutil.getcolormap(**kwargs) - if len(args) > 0: - level_arg = args[0] - if isinstance(level_arg, int): - cn = level_arg - ls = LegendManage.createLegendScheme(data.min(), data.max(), cn, cmap) - else: - if isinstance(level_arg, NDArray): - level_arg = level_arg.aslist() - ls = LegendManage.createLegendScheme(data.min(), data.max(), level_arg, cmap) - else: - ls = LegendManage.createLegendScheme(data.min(), data.max(), cmap) - ls = ls.convertTo(ShapeTypes.POLYGON) - edge = kwargs.pop('edge', True) - kwargs['edge'] = edge - plotutil.setlegendscheme(ls, **kwargs) - - xslice = kwargs.pop('xslice', []) - if isinstance(xslice, numbers.Number): - xslice = [xslice] - yslice = kwargs.pop('yslice', []) - if isinstance(yslice, numbers.Number): - yslice = [yslice] - zslice = kwargs.pop('zslice', []) - if isinstance(zslice, numbers.Number): - zslice = [zslice] - graphics = JOGLUtil.slice(data.asarray(), x.asarray(), y.asarray(), z.asarray(), xslice, \ - yslice, zslice, ls) - visible = kwargs.pop('visible', True) - if visible: - for gg in graphics: - self.add_graphic(gg) - return graphics + return self.slice(*args, **kwargs) def contourslice(self, *args, **kwargs): """ @@ -1218,54 +1166,7 @@ class Axes3DGL(Axes3D): :returns: Legend """ - warnings.warn("plot_surface is deprecated", DeprecationWarning) - if len(args) <= 2: - x = args[0].dimvalue(1) - y = args[0].dimvalue(0) - x, y = np.meshgrid(x, y) - z = args[0] - args = args[1:] - else: - x = args[0] - y = args[1] - z = args[2] - args = args[3:] - - if kwargs.has_key('colors'): - cn = len(kwargs['colors']) - else: - cn = None - cmap = plotutil.getcolormap(**kwargs) - if len(args) > 0: - level_arg = args[0] - if isinstance(level_arg, int): - cn = level_arg - ls = LegendManage.createLegendScheme(z.min(), z.max(), cn, cmap) - else: - if isinstance(level_arg, NDArray): - level_arg = level_arg.aslist() - ls = LegendManage.createLegendScheme(z.min(), z.max(), level_arg, cmap) - else: - if cn is None: - ls = LegendManage.createLegendScheme(z.min(), z.max(), cmap) - else: - ls = LegendManage.createLegendScheme(z.min(), z.max(), cn, cmap) - ls = ls.convertTo(ShapeTypes.POLYGON) - facecolor = kwargs.pop('facecolor', None) - face_interp = None - if not facecolor is None: - face_interp = (facecolor == 'interp') - if not face_interp: - if not facecolor in ['flat','texturemap','none']: - kwargs['facecolor'] = facecolor - plotutil.setlegendscheme(ls, **kwargs) - graphics = JOGLUtil.surface(x.asarray(), y.asarray(), z.asarray(), ls) - if face_interp: - graphics.setFaceInterp(face_interp) - visible = kwargs.pop('visible', True) - if visible: - self.add_graphic(graphics) - return graphics + return self.surf(*args, **kwargs) def isosurface(self, *args, **kwargs): """ @@ -1343,54 +1244,7 @@ class Axes3DGL(Axes3D): :returns: Legend """ - warnings.warn("plot_isosurface is deprecated", DeprecationWarning) - if len(args) <= 3: - x = args[0].dimvalue(2) - y = args[0].dimvalue(1) - z = args[0].dimvalue(0) - data = args[0] - isovalue = args[1] - args = args[2:] - else: - x = args[0] - y = args[1] - z = args[2] - data = args[3] - isovalue = args[4] - args = args[5:] - cmap = plotutil.getcolormap(**kwargs) - cvalue = kwargs.pop('cvalue', None) - if not cvalue is None: - if len(args) > 0: - level_arg = args[0] - if isinstance(level_arg, int): - cn = level_arg - ls = LegendManage.createLegendScheme(data.min(), data.max(), cn, cmap) - else: - if isinstance(level_arg, NDArray): - level_arg = level_arg.aslist() - ls = LegendManage.createLegendScheme(data.min(), data.max(), level_arg, cmap) - else: - ls = LegendManage.createLegendScheme(data.min(), data.max(), cmap) - ls = ls.convertTo(ShapeTypes.POLYGON) - edge = kwargs.pop('edge', True) - kwargs['edge'] = edge - plotutil.setlegendscheme(ls, **kwargs) - else: - ls = plotutil.getlegendbreak('polygon', **kwargs)[0] - nthread = kwargs.pop('nthread', None) - if nthread is None: - graphics = JOGLUtil.isosurface(data.asarray(), x.asarray(), y.asarray(), z.asarray(), isovalue, ls) - else: - data = data.asarray().copyIfView() - x = x.asarray().copyIfView() - y = y.asarray().copyIfView() - z = z.asarray().copyIfView() - graphics = JOGLUtil.isosurface(data, x, y, z, isovalue, ls, nthread) - visible = kwargs.pop('visible', True) - if visible: - self.add_graphic(graphics) - return graphics + return self.isosurface(*args, **kwargs) def particles(self, *args, **kwargs): """ @@ -1473,51 +1327,7 @@ class Axes3DGL(Axes3D): :returns: Legend """ - warnings.warn("plot_particles is deprecated", DeprecationWarning) - if len(args) <= 3: - x = args[0].dimvalue(2) - y = args[0].dimvalue(1) - z = args[0].dimvalue(0) - data = args[0] - args = args[1:] - else: - x = args[0] - y = args[1] - z = args[2] - data = args[3] - args = args[4:] - cmap = plotutil.getcolormap(**kwargs) - vmin = kwargs.pop('vmin', data.min()) - vmax = kwargs.pop('vmax', data.max()) - if vmin >= vmax: - raise ValueError("Minimum value larger than maximum value") - - if len(args) > 0: - level_arg = args[0] - if isinstance(level_arg, int): - cn = level_arg - ls = LegendManage.createLegendScheme(vmin, vmax, cn, cmap) - else: - if isinstance(level_arg, NDArray): - level_arg = level_arg.aslist() - ls = LegendManage.createLegendScheme(vmin, vmax, level_arg, cmap) - else: - ls = LegendManage.createLegendScheme(vmin, vmax, cmap) - plotutil.setlegendscheme(ls, **kwargs) - alpha_min = kwargs.pop('alpha_min', 0.1) - alpha_max = kwargs.pop('alpha_max', 0.6) - density = kwargs.pop('density', 2) - graphics = JOGLUtil.particles(data.asarray(), x.asarray(), y.asarray(), z.asarray(), ls, \ - alpha_min, alpha_max, density) - s = kwargs.pop('s', None) - if s is None: - s = kwargs.pop('size', None) - if not s is None: - graphics.setPointSize(s) - visible = kwargs.pop('visible', True) - if visible: - self.add_graphic(graphics) - return graphics + return self.particles(*args, **kwargs) def volumeplot(self, *args, **kwargs): """ diff --git a/meteoinfo-math/src/main/java/org/meteoinfo/math/optimize/MyParametricUnivariateFunction.java b/meteoinfo-math/src/main/java/org/meteoinfo/math/optimize/MyParametricUnivariateFunction.java new file mode 100644 index 00000000..33874d76 --- /dev/null +++ b/meteoinfo-math/src/main/java/org/meteoinfo/math/optimize/MyParametricUnivariateFunction.java @@ -0,0 +1,57 @@ +package org.meteoinfo.math.optimize; + +import org.apache.commons.math3.analysis.ParametricUnivariateFunction; +import org.apache.commons.math3.analysis.differentiation.DerivativeStructure; +import org.apache.commons.math3.analysis.differentiation.FiniteDifferencesDifferentiator; +import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction; + +public class MyParametricUnivariateFunction implements ParametricUnivariateFunction { + + private ParamUnivariateFunction function; + private int nbPoints; + private double stepSize; + + /** + * Constructor + * @param function ParamUnivariateFunction + * @param nbPoints Number of points for difference calculation + * @param stepSize Step size for difference calculation + */ + public MyParametricUnivariateFunction(ParamUnivariateFunction function, int nbPoints, double stepSize) { + this.function = function; + this.nbPoints = nbPoints; + this.stepSize = stepSize; + } + + @Override + public double value(double v, double... parameters) { + function.setParameters(parameters); + return function.value(v); + } + + @Override + public double[] gradient(double v, double... parameters) { + function.setParameters(parameters); + + // create a differentiator + FiniteDifferencesDifferentiator differentiator = + new FiniteDifferencesDifferentiator(nbPoints, stepSize); + + // create a new function that computes both the value and the derivatives + // using DerivativeStructure + UnivariateDifferentiableFunction diffFunc = differentiator.differentiate(function); + + double y = function.value(v); + int n = parameters.length; + double[] gradients = new double[n]; + for (int i = 0; i < n; i++) { + DerivativeStructure xDS = new DerivativeStructure(n, 1, i, parameters[i]); + DerivativeStructure yDS = diffFunc.value(xDS); + int[] idx = new int[n]; + idx[i] = 1; + gradients[i] = yDS.getPartialDerivative(idx); + } + + return gradients; + } +} diff --git a/meteoinfo-math/src/main/java/org/meteoinfo/math/optimize/OptimizeUtil.java b/meteoinfo-math/src/main/java/org/meteoinfo/math/optimize/OptimizeUtil.java index b5eb7c7f..ac9f8899 100644 --- a/meteoinfo-math/src/main/java/org/meteoinfo/math/optimize/OptimizeUtil.java +++ b/meteoinfo-math/src/main/java/org/meteoinfo/math/optimize/OptimizeUtil.java @@ -1,9 +1,12 @@ package org.meteoinfo.math.optimize; +import org.apache.commons.math3.analysis.ParametricUnivariateFunction; import org.apache.commons.math3.analysis.UnivariateFunction; import org.apache.commons.math3.analysis.differentiation.DerivativeStructure; import org.apache.commons.math3.analysis.differentiation.FiniteDifferencesDifferentiator; import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction; +import org.apache.commons.math3.fitting.SimpleCurveFitter; +import org.apache.commons.math3.fitting.WeightedObservedPoints; import org.apache.commons.math3.fitting.leastsquares.MultivariateJacobianFunction; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.ArrayRealVector; @@ -104,4 +107,31 @@ public class OptimizeUtil { return jacobianFunc; } + + /** + * Get curve fitting parameters. + * @param func The uni-variate function + * @param x X values + * @param y Y values + * @param nbPoints Number of points for difference calculation + * @param stepSize Step size for difference calculation + * @param guess Guess values + * @return Curve fitting parameters + */ + public static double[] curveFit(ParamUnivariateFunction func, Array x, Array y, + int nbPoints, double stepSize, double[] guess) { + ParametricUnivariateFunction function = new MyParametricUnivariateFunction(func, nbPoints, stepSize); + SimpleCurveFitter curveFitter = SimpleCurveFitter.create(function, guess); + + IndexIterator xIter = x.getIndexIterator(); + IndexIterator yIter = y.getIndexIterator(); + WeightedObservedPoints observedPoints = new WeightedObservedPoints(); + while (xIter.hasNext()){ + observedPoints.add(xIter.getDoubleNext(), yIter.getDoubleNext()); + } + + double[] best = curveFitter.fit(observedPoints.toList()); + + return best; + } }