From 225e5c703c7b55d802f3e9a6aeb92c7cd27ad386 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Tue, 5 Feb 2008 18:09:06 +0000 Subject: [PATCH] completely rewrote estimation package documentation with downloadable example and explanation diagrams git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@618726 13f79535-47bb-0310-9956-ffa450edef68 --- .../TrajectoryDeterminationProblem.java | 187 +++++++++ .../userguide/estimation-class-diagram.png | Bin 0 -> 20220 bytes .../userguide/estimation-sequence-diagram.png | Bin 0 -> 16190 bytes xdocs/userguide/estimation.xml | 394 ++++++++++++++---- 4 files changed, 509 insertions(+), 72 deletions(-) create mode 100644 src/site/resources/userguide/TrajectoryDeterminationProblem.java create mode 100644 src/site/resources/userguide/estimation-class-diagram.png create mode 100644 src/site/resources/userguide/estimation-sequence-diagram.png diff --git a/src/site/resources/userguide/TrajectoryDeterminationProblem.java b/src/site/resources/userguide/TrajectoryDeterminationProblem.java new file mode 100644 index 000000000..32a0975f1 --- /dev/null +++ b/src/site/resources/userguide/TrajectoryDeterminationProblem.java @@ -0,0 +1,187 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +import org.apache.commons.math.estimation.EstimationException; +import org.apache.commons.math.estimation.EstimatedParameter; +import org.apache.commons.math.estimation.EstimationProblem; +import org.apache.commons.math.estimation.LevenbergMarquardtEstimator; +import org.apache.commons.math.estimation.SimpleEstimationProblem; +import org.apache.commons.math.estimation.WeightedMeasurement; + +public class TrajectoryDeterminationProblem extends SimpleEstimationProblem { + + public static void main(String[] args) { + try { + TrajectoryDeterminationProblem problem = + new TrajectoryDeterminationProblem(0.0, 100.0, 800.0, 1.0, 0.0); + + double[][] distances = { + { 0.0, 806.5849 }, { 20.0, 796.8148 }, { 40.0, 791.0833 }, { 60.0, 789.6712 }, + { 80.0, 793.1334 }, { 100.0, 797.7248 }, { 120.0, 803.2785 }, { 140.0, 813.4939 }, + { 160.0, 826.9295 }, { 180.0, 844.0640 }, { 200.0, 863.3829 }, { 220.0, 883.3143 }, + { 240.0, 908.6867 }, { 260.0, 934.8561 }, { 280.0, 964.0730 }, { 300.0, 992.1033 }, + { 320.0, 1023.998 }, { 340.0, 1057.439 }, { 360.0, 1091.912 }, { 380.0, 1125.968 }, + { 400.0, 1162.789 }, { 420.0, 1201.517 }, { 440.0, 1239.176 }, { 460.0, 1279.347 } }; + for (int i = 0; i < distances.length; ++i) { + problem.addDistanceMeasurement(1.0, distances[i][0], distances[i][1]); + }; + + double[][] angles = { + { 10.0, 1.415423 }, { 30.0, 1.352643 }, { 50.0, 1.289290 }, { 70.0, 1.225249 }, + { 90.0, 1.161203 }, {110.0, 1.098538 }, {130.0, 1.036263 }, {150.0, 0.976052 }, + {170.0, 0.917921 }, {190.0, 0.861830 }, {210.0, 0.808237 }, {230.0, 0.757043 }, + {250.0, 0.708650 }, {270.0, 0.662949 }, {290.0, 0.619903 }, {310.0, 0.579160 }, + {330.0, 0.541033 }, {350.0, 0.505590 }, {370.0, 0.471746 }, {390.0, 0.440155 }, + {410.0, 0.410522 }, {430.0, 0.382701 }, {450.0, 0.356957 }, {470.0, 0.332400 } }; + for (int i = 0; i < angles.length; ++i) { + problem.addAngularMeasurement(3.0e7, angles[i][0], angles[i][1]); + }; + + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(problem); + System.out.println("initial position: " + problem.getX0() + " " + problem.getY0()); + System.out.println("velocity: " + problem.getVx0() + " " + problem.getVy0()); + + } catch (EstimationException ee) { + System.err.println(ee.getMessage()); + } + } + + public TrajectoryDeterminationProblem(double t0, + double x0Guess, double y0Guess, + double vx0Guess, double vy0Guess) { + this.t0 = t0; + x0 = new EstimatedParameter( "x0", x0Guess); + y0 = new EstimatedParameter( "y0", y0Guess); + vx0 = new EstimatedParameter("vx0", vx0Guess); + vy0 = new EstimatedParameter("vy0", vy0Guess); + + // inform the base class about the parameters + addParameter(x0); + addParameter(y0); + addParameter(vx0); + addParameter(vy0); + + } + + public double getX0() { + return x0.getEstimate(); + } + + public double getY0() { + return y0.getEstimate(); + } + + public double getVx0() { + return vx0.getEstimate(); + } + + public double getVy0() { + return vy0.getEstimate(); + } + + public void addAngularMeasurement(double wi, double ti, double ai) { + // let the base class handle the measurement + addMeasurement(new AngularMeasurement(wi, ti, ai)); + } + + public void addDistanceMeasurement(double wi, double ti, double di) { + // let the base class handle the measurement + addMeasurement(new DistanceMeasurement(wi, ti, di)); + } + + public double x(double t) { + return x0.getEstimate() + (t - t0) * vx0.getEstimate(); + } + + public double y(double t) { + return y0.getEstimate() + (t - t0) * vy0.getEstimate(); + } + + private class AngularMeasurement extends WeightedMeasurement { + + public AngularMeasurement(double weight, double t, double angle) { + super(weight, angle); + this.t = t; + } + + public double getTheoreticalValue() { + return Math.atan2(y(t), x(t)); + } + + public double getPartial(EstimatedParameter parameter) { + double xt = x(t); + double yt = y(t); + double r = Math.sqrt(xt * xt + yt * yt); + double u = yt / (r + xt); + double c = 2 * u / (1 + u * u); + if (parameter == x0) { + return -c; + } else if (parameter == vx0) { + return -c * t; + } else if (parameter == y0) { + return c * xt / yt; + } else { + return c * t * xt / yt; + } + } + + private final double t; + private static final long serialVersionUID = -5990040582592763282L; + + } + + private class DistanceMeasurement extends WeightedMeasurement { + + public DistanceMeasurement(double weight, double t, double angle) { + super(weight, angle); + this.t = t; + } + + public double getTheoreticalValue() { + double xt = x(t); + double yt = y(t); + return Math.sqrt(xt * xt + yt * yt); + } + + public double getPartial(EstimatedParameter parameter) { + double xt = x(t); + double yt = y(t); + double r = Math.sqrt(xt * xt + yt * yt); + if (parameter == x0) { + return xt / r; + } else if (parameter == vx0) { + return xt * t / r; + } else if (parameter == y0) { + return yt / r; + } else { + return yt * t / r; + } + } + + private final double t; + private static final long serialVersionUID = 3257286197740459503L; + + } + + private double t0; + private EstimatedParameter x0; + private EstimatedParameter y0; + private EstimatedParameter vx0; + private EstimatedParameter vy0; + +} diff --git a/src/site/resources/userguide/estimation-class-diagram.png b/src/site/resources/userguide/estimation-class-diagram.png new file mode 100644 index 0000000000000000000000000000000000000000..fe9b6786a0e6adc06b35cf37b60678f893840766 GIT binary patch literal 20220 zcmY(r2V7Ij7cETaO1~mVR}lpP3B7lz(nO_00tkWx5Rl#i2r5;i1?jzubg2m-H339G zI!KoiTIlsV;okpy-}|W;&zUo4X3w5IYpuQGM?D=is_V?x2?z+NG}IsK6A*xwfj=+F zNr5Z->j+i?0-bh^$BNIrCa`H_)$}iF4}H5OSN?o@s;8%TEqw9zc;q}cDLv?B7ku$P zg=YV%#h=fQ7#~T_sH-Ye2H*aJ@N%B<6}I{!K5OEfR(@C{HS%I)WctPQ2=MFj(!_VV zc(Y;prVbXe)vTMg;HpSV{J-DyMwFxxpYZ?6f*&7Xm6=ah3;a3)&Yl{v*Kq{UCZe&Y zTcJyDPB*1%Y^spw+K$E>lg>3B>%SGFrklJKti_!^VGTD@oVyQP8UxH7ee2B|@55+L zPaMk3n;T564Zm*u8Z(4ohsTm2pRsxyyEc{BaPx+HFqXOq#*xIHM`n7H1|Anpqy1@< z^VY8hZ0$B81jKkfmLtEtbt=%MnEZhq=CNd{_rIv05c?MQpvGhVhffrne2b2CS9EM1 zkH;SOxQUZmG`VTc*j#&?fg2Ix4tX$)5T`586eKBd7y4US@)m-N`;@GTu~!v3 z))OOi%cnt%=!p?63$%a<-5U6TX5+UoKHMpwD&i&O_XM`&!I_0{@9jUYzPG1|u5(f^ z-Br1c+dsysxp9lzyE7dNsK~_ zArOmtS0T^u19Fog{6=6C5H5@=qtc-4O9ckw&DvX_lUH?sHTC&p@L{C*MHHLN_s=28 znI_+$pv}M6F!ae{&NB$-!@Z$tqn^(l$G%q5!Ao5n$o_|hj|OJZ{T3Y%#@Q!Z6yO`E z$r^`N`4Zc{V%Xtk2WvMOHQ%egFQW71hg-86$7UcKO~V_D=p=#mOUUC+E9uzI7B3fR zUOfdL?c7QMv$0Odjnm!j_QQ(4lvlh^cG(lEffhr(ix{|8X5N0bO2nJnH|#p3Kn8re zh;4?(^{;Dwr_H`em12%xa2yn(ubLzcZ4VNZ8&abPW zw<7lDBeU+e{!GxY+vJtkFs*B7Vz!AE&#eR2+&25TxWMh!`>eg<&g_eH+=@2GyC2Wa z87=;Hv}*my->-GtS}C%d2!m^Er;0fe2gdxlQ_Fd)+(I^?pvvBV?Ayxt?=rVzp>tYV zy%rz$qz7NhT3?TT?PL&O@0tzVtY?=w$s7GEO+FlVpK4``PJsk5G$g+;s}cJ)^t#P> zwu<#Q92A4YvitnCAi?Y{D5tgM9}d&uHpV%R*VCly=AkpS+EQUB`wk9wmuegecup0d z7FD)icD>)0R+|++zcvt)`H{7k9yH{g&B2o4e*g-+{P>{Z4Y(jzYkocNNz$Txye-Bn zrrfYdKd3x#_hYq={O7+iw1pjFdz_)DCacht@*9yW{cJ^N1DT5F@rUeQf7qYw+rXKE zZ*yxxFHY8*r{+Ac6JZ5$*SxA{>YUYR)bG211&oDe3aF)_3cBofp7DSAeQgb$sG);| zAa!J-WSwWadVQ^qiwYWc=?|vrXWWZdf9=NUdb~hsxkJxA<#1Y_)W%6>()l6>rN$h+ z(DhRIan*p_;yQh-{-V|@3{-G&F3oGdkt{rG2+!X<+Fwo?c#|3N_jC07Y#YKVuPFUX zLd;MBvLIEM5VyC=xwua)(|hxj_tzL|KGUJUpJg|zHY?Wsj z1DRV7z`9>VBFmp)Gk%THdk^^xq+{a`?Z!2kS*{#|3UjfJ&K4KT^8$RXJ4ZrOdPqto zG`m*jWw)D9)9^){ODx3XzL-GSddWIKhcOW zLvfoP8gRZ}aO^-tcJO;RspFCkMvmE&E7MJt&;sEi1P*oxHrc3ccbKfC>wF6dmxG9W zdrqazf`szMS4Afv`6D&*8#WnaX^}-2G8WL`Q>!*vVi{xqv)QOMEGi9CP=E()!xKP_Q=@ic*o5 z-X6@NhHio8R~LOSQfc;F8}_6JN5ntM`{F#^`wL$0Fj}kms$#Ak@&~eYsYlBYw zJeyLw7hqBIT8PXo>5=1+l*&?!{7kgKK(eqtNH}PQt03eWD!V=DD0aJFTA(tNqN`HK zCU?uI78Xosdjiu;pzyMJiLp(Ss+mXApHqN6H>XTixgMdkl8Jnxvh90s|J>ae=it$k z^K{}ZcW(4z3%E~)n{-CrzU;o$sv7O5JyPXQzy5RXj*T(g`0+OyIWL-+@DXS}eT@Vo zvF*-fwGNbFeZC7tL^IVV?NQXK%@U@m$;cc4oc6|7_m|^sjSU zK@a=`>itjO$E)T46>xXTwiqVF5JL%r2&kDr@fM8RTweMyuh$UIb5w7N`RqZd{Rm}2 zFJNk-R)keX57hOs&(9<>5A&FD?CA{#fe%{aeRQ8;zdKEv$UJ1YRD;t)?;AZ^EDi`T zCP?*BqX0*Y4%b?BJT@rOm!x*|S=7w^T7+Y1O*?o$Q_yXaTUk|QxXbKX@0uqo&zio3 z-jg3TZO~$}za&QxV)l&V8w@yD(xeHi#?qx)=nx${bTZd4RV&+DYuwS%$wfx6@m?zEr$J0eTa@Wxw)a^`8 z{^MFd_~2SOiGx+M1j59MpGYx?Wx64mF!8<|X3TK=B!}I;LVgSBTb-*n!=bAFOhYpH z#_VCjw1{(%@vW;j6Cy(LDx&heBXyXi&$@E&fl~stt`Et4A#wxXCc*$1i~$vS;tE?s zSEMcvMxO`;C9BE^6TQ5OganAR=Md_|PMzG7$~eG2IW48(h~qzkLGV%{;pEB8s0Q;E z*_H~n+IEkz{-VjHNy=rtHT1aUXc8AU#|L%hNmhjC6n z8~Q<#%tJD&uW#vS;r>RRR5JH6v;u*Im??y*{f*TZHQzP#@i7nn9eTMa*daJo=c~Yb zamrV=TJPtoO(4<3*7{e=MwNyHQ$7GXuWCt6A_FV$+(4k9k5X2W>)J0Zzo=>9?9kv; z0%7_XXp8^pC#$SJZH_5%lbwjHS9F|%biaKbjBR(hfDQX8y2wm+0$swhZP*yhI^4_HDT`Mbnwe{&oR+Nw3G#Tl2tx}o$Yvf>S{A+G$;*r2S zBjz~&{=T55iwTv{kIPVoQsTpp281^&1ll-?V@mY0f=dm`ied9V<`VraB|zZ`f9{<{ zmpi0U*bXNN_T{kQ56N&fJc8VJBO|v<_I?-dXE|P(X9R#Y$5fkv(Fch<%Ck@J%}(jsCwAOr0Fdx{qp*ol#`&Zvxs7u!FxL^oPa znRo0E%TEkLlKrxnF%*o^OX&TzkzanEleIM8s;gJ0$iheE&-eTT&KAjPCT9kImnP>Q zDeTMu@p=bqpE7FU+lw=|HtPHG!g#tt)L66k<>n6CQ!WQLT?HMYhx;LRs$O)m&berw;@ZX@fTV3ipzPvKA7TVB zecG{g{8*x88hO8tszBZ&*C0ZYvF_M%7upDH&!zk7uP@~h=dg+D7T4(RcQhMQ^*152 zl{VCulerrH9hV;+)uR1GOD#WGlmEMh?Ii+47(KsW$=Ttly~P_NbX8Ju>QC=A&lEmo z8)$AE7ea(W=-1mGgahMVA02>3v5f9KZVZhdjY5^iJhV4IQr7f-LbpSmF1tr&+bp1^ zoz})Ib|xHC-68fGTyp-WHz7`Dawc!wfW~%?sH0&!#47D9W6-6rmCuhurFRW$XHK*H zMpgP=2BZ6b9^Wq#vxbkr9W02lgU$mo+g(O}6OjbYM$&Nk>KhmUbuXjUqI2z`r8?Fn zty*##TQnqmq4FXW%WETfTmAhOO!S!e&jq8@tVT{0h#tct|4hKC1Am<=2t&ng#7mI> z8o zxV6P!dUqdy8{0eYY2?PP|8S0gLD3jJy7K-7inpb&rqaIl=FRQqp<9q#;vZ zKfmE(8K3-^;|^IcD}k+BjK)WJht1LV=VH1x&)%uAeJSLOBftSDZr9?+*V}1Gzt$#Q zvgKyj&b&*hRflV7;AMZQ^85D*(Y8*9;3a7Q%!L5_f%?nyZv<|tov<0jVvFVp7{D*3 zz|Myh%k^&Zwl^Qf+xVZn{lhj;_^4ElY>gmP%44v$*|%;K=5P}uDOqB4xEW&?9fH!@ z|Er&;$^BG#q`1-Dp+v*}G+44)D4ks0?5UqwZ-su}P|Ed7kfPo_6l_V5hnK6RTPDxS znRM{$U`zHf;(1|Aco+^rkK%Kf|B*$+bF=bgkJo@enf&=ASKlNB0Qf|~?{T@|=Iw97 zo=sHg{}z^?7_I#fb#N2*CQn^^%DQ&Q;a5m%?JbN= z;4G7;k{Y^=R#cgmKKTPvm4>AlkDD@W=URj56S%?cn4SX1 zM=WVoa#!v~p@6idSkaqKSM~hy+&ydYdt8Wf74(NMcgVx35jzNa_PAYnVi)~jNxkYw zprcH?Af_0NIYFMBaEzCqk)*kP%Y2Y#Dd2@br}OY7lZYyH1x zHlJQ=ep=}Ev}AvF3!~v7ECOB-EPCmyrCCWY^C$sM_JN+9xI^Pd_VHK4p;2MojAshO}f`#kl!nfQ?}a4&N#+{f$k zNlW16QswUj7R{M;Qt*PrZSrvFa9(s1t>{8S>ulrCV+x?5syWr-X&orr)~xQ{fVq6O zl4L`^*a}o@x>4}}w0w1o`seY#5Jvv|apC~Biz9TYX!&c6r;_)2Tn`!hGl}7ycgk#w z@;@=swgv6N@?qa(iB4Il)@Z&RMHL#s9d{P)%ayZy4<(8e2d}K7Qdfktnr?#MOUYVf zwqEr0qpbZ+nxG+P>jgg|FmJtVR|;#~0!$P_2Ue$0i|3OHGjc_&uE7@KZz!N1Vu+DJ z*W2KfjE`OE#$M#z)fsa-_c&k&X#;Yn@E7G4$$?31&*nqPY)@^c1k*{KsZg-3^9URn zGKC(CiX#0-8muN2)$G1*_PfEo3U0p<89841qQo{4WmakwOFJ|sZIq=%BFT+a2ab)Y zN}FeiU!*!{dkHbX>0-&l6Q`vc`sx+lCEmD9qxeN=yMo6hOsH>SRPABcy8DbLYY(gV zO_!O2B{{9B$y`&!`Xp1MP+IFHhPpO24u;mAv*t;T6CzKZtv?I#o^8o1b(;QC`m$#p zAC@(lQ))S(Voqj;I3yv)s@ZGD3f2`JdL2=MxBL3` zP~}r*Uo-<~@l-N(q{aqhOUFOtlOe|@qc?%t#TmMJ{S&#bcRyzU09_2|K5a_?I~io# z>@e*~jUW1?RaokGdkghnx$_J!k!6 zorQh&Si;g;B2ikgA(9VX-E^3lnmm32mz4J9A}cv!2iX!fjZ$GrpdPmnIpPyZ35qJ<2}6J!mjD(tx#~&MqLnpbum$rNJr{jo;||-< zYvoh)uCipFY#D6%r4@CMr)y1wArSt`G7({1Kd6vRVXI{lD-&l>I0T$RsK|9vCt$6a zF^KBBlK688eRITci0{VXQE~wrBWhtxsB+{uP`Td`)NW_b%cB~G!}&7 zYvKHvd(#$GDyQS~wir*E+TLZ|N0o+F9#Q)H8$U+RpWA7hfUaz=Xm{HMv@A0O6-()S zWLitlYfx(Um?Q4tvwZrD7riqgwG0?~{-N%d)mof>#{OEGe$fC45J2Q4b;`Zvz6lKr z>2eRz>Hz{13z)1};FJiLD)~Qp^edBv*lK@l@TDqUm61g)rmy%#$&brde}AqRzb+PW zdwzaWW9)OnXQM|4QvITCHIzT#*V4{ymH+XPUFG#)amU$`QEOH^D?K_!j;>QGyiy}4 zp>De@!xdFVDoiz?fp%a2t#6SbiSrMbBDVgbmGf$#%V;qB*WupZ9c~~gF#`9V|6y6` zwUh1lyaKP)+QUQl3Y<9 z8oOWoAwcNpXi&F~a3^##nRsx9vzLJ-y`EWuLK)Z8n|9?T-PJo~)rkS7q@^2U-S zqcCX3VWL`59m){*;`?)9jiOew995eZzPqM@(EGQ@h=f&;)tfY!><1p-q4`^bA<2K9 zpyPi5aZ5_a-Ct4!Or~qSvhv@Sp+K=RLdATvdT(KjVSoVZHSo#@&ugwD?2|1x2|nsT zvf#=3W6wGBZNBg+H3b}OFwp!7?@`difrPQokO6E7&&5Wd_Vv!Q(X}DT_Z#2rJU>te zP69hFo3Et8ea}T5CmQ{dZHxYi0fwoy>TH6Sj+IEO67X>tywH%0`sk5r6nv-oKs%@?_T$fO~dZYPnnX zKHM0u8Y9547#aM&F86{f?_0;ZN3e0mekn(+UGWL0MA zcoC1VA`y>5X!KiVX#tx5y1?MWjyAZD;{{GkK{$Y>(B zx9o9B%#+6d3kd)w35GMp8T`1+y}9WRM6h6>SJ4lbHeW(FGn1Ohe?W0iTgp7Nly37hRk0^w*P*H8C# zSf#$%4c>+IeBLQl2m%NsPM;_>;X8vl(v3Z3vqA~93Gumn*R|Qn6EaJJaAEz_p1l-s zF16E8gQtzrfpnh?y*8@5idQ!wCIfNhZak?EkEUPlA3BN9cJ~bK0?M%jVpSk}`Q5J} z_Qs`(ncefrcI#1DF@gfp0CA1)p`?Tz&eeRBPu<#zqPrfU4+YLu)|dkrm574o@8C(J zKXfKK-tT$jSnfz3AYQbn3UtV2qB_iam~AntW~N`hp&zh0YsGGQ+u4w!fHFWl zs6Js^fd~R{%!$Uz%w5Z4%5ex8aUpB!55xlLh zXo?wxe3X$8k#tZ5@8b3uAFG1|pXcF8>e>7Dbu!zdM~Xjh=Ok%NC@<2Knfa)3Qy>BDF7ly>vU z5fPlJZ9kNi0Wd|U|Fu8@u=4koPa*j6eHWOc0rKgijILLl_>k)?h$0PnyZTwTQ(SUjm>u+6f+d zuBUxRms$$Pxm_1aN>NO}GMDG8EoNY6egm*Rk zXz4{cO@PtjNgZ(W&c8Re$sWDE2>z*aXXP>Q@SNuy&6SPRha>Hn{vA@cBEIF3IXM+^ z@a^}H`xt5Y+qaZ}u@Fo0zE-WKL{jv)=HSh|P?&tjBy3h&$%#wmCxPvXIn6CTcFJ&h zWW);Kxq95A3;z~eJgEG);v(Hj{_Hki(&XW63IF>^&)I9}<%2q85r2gLY1wb5X{o=~ zKz{h=SQJa!0q;Sdc-#LixyY}T-B>RlsZB9&F0SD<^f>;%H5Wg}me4OJZ#KTwaGd;Z z-c(qtysqxsHV&IVW6;rl`_$k&JTfgf6l*j9FQvFRkJWYxV##uKYKRVG9 z$0C9+y5Nuw@JowflTOU9$WiP0{m^6a!#B}zHaR|+Xr1!=CQ>8_yzpaQf?$~dD5fx` zu$7JH#Q&&)iiNr6qWBcYIS_SBa7ke+XksSoBZT-n;$`bD0$}~&OkJRH%4I(y@Dn{C zyUXM-Udh9&p9^*v+a7?e{p`ZCV+^_vB}wcwQUEyV<*r9J`Sk{ zz@HhoB&G_lm*Xo20E7|-Jy4RIChMISJ_|B#;&=#@fD|w(c$Y!A9a>Yb10*XNh+M?s zB8dILN_LE>gF$TC-hJxSB#Skcdr^**z`(J@RX}DUS?Eek_sh1y(YbAli217xW~*rB zR|xxO!bmaIP2M=Z7A<*1X?{h=%HN@}gJ4Nu-$SDvkn5dAYD7UW(ki(q^MfNf7sIt zMl7jBbXFF>;N7`o5?YkR=lli^4Xpv1UdbK%+9Nvqd3+5o>OEy(HVNkJQ`2p)p z_4z_f-2(jX?Fg)D6f@?b0biG9B459dkS*`3gD3*X&-Wy)|1K+Py6AbZwPXDNyt0_P zJ-{$Kl+(qTK5Af(+VG&U>4;!l8g=OVVy}5((xgiZEax?ckzzn^{|>&4t@S%5rnau= zi)0Gi_j^S$a2h%W3paU&k+ zu;weg`y=7KD=m|qvw!_&?`u?(2llP*qIilxo!#IMGK{RRT<7KT%9H2@Z(Jwn**LK+ zNR9AQOH1jMZqI&BO7>BLRFjcIj;}090HXRa5vV`-@&R84MCb-&;g5EMNpfVx!Ub+9 zj6g5`<`^0~IbDa`=U%<z$cu12kVqiqF|gA7f;FWQxy@?&K_OPIC;xy#cxj8?U_P zS8rUU8V2125H2L2LuIc{BkqJ59CK@jvDRT!sl_Z)0d#P2D&TvdVIp8!D}n;z1aNx% zZrzZN;+G0=92YXwu*ZwYGzcDlta+mTsFgX8wQpd?-1W!Kiw-HtL4T`dR)5_*%8Nqe zi2xN?ea8;ZX$q=eb3+;IVHJxeBn``0_NsztCQvnnw{VL-HZP*?|=i2%5E-)74v^JJH~9WU79 z$4P)uL}}g6RPLXCh!_C_F#mOTwHA*f0A#_;``y6fc}NbwZFa7PXPzgO6qxvvJY;=+ z313k1YoO{X|F1YRfN-;@A+Nz15_f_S`V&DYt)I_Fllp=s_U4<^N_aRxXISYm%|5dExs^ON z{F+3s!OtB1004XBI{MX-vXxJNz8rZGpWhB}W_wbo1q%a0WtFh*wb_ZexDH=5QTS~H zkG@5dm(14D<09K)&*h~u10U-YBDchUYqVlE3O)SJj(j-48M2KD``F2+ale3Ct_ki0 z?a62F@2%a=Bdwa0jr6p2BUwGmBdY@R%Dmz}{pQqzD5-*{jxQefJ&$?xqWE2E!vTq~ zVz`3Vsk$oqq#xn|PnHU#OYJzMjnC_c{DW9MdyAb&-0~gD-}g=%Bgc@@JbE?VP&1NS z4YCCQ&%q&(#3Tarx{0LYAQiI)=Bmpel{OXt2@*y8<%kPkLfxHFRJkkt?Zd_KlAd$` z2w>*G^f*7ZraDU_`Y8GCW`9=?k^G6eI~O9EYAL}Ms8(=0aanBq2kE&^ zvM;rU9xM@1qrg~j7o3K=R!aqhnnDbMB7Y#ignc&6hJjah^PpBO7!;jgO_gpl3jEA zOR{$qB>_T4v}EUg=Bnv$Y-k<5*i3*Vbpm0=kROO9B5w#8L;7*1kZ6N00fG3#(=q|z z@MGP#@kMP+om6I6jK_Xw{4Ep8p}Wk__Lt7i z;uQu2CjpC;2bpT{iu~yYohXs(knb#WycnmpV+g(x{>FcN{Dz@K zQOm~OM_#KJ5HVnRHc5`Inn*5qEo8O8*f74VBHzf@&~tD}pk$=g8C9-mhti<5eV}bd z$FXlypo?GsN_TJp0cP7FI*|qFU6KvcBmAw~L(7aV4~R{yqo*r6US#15s0$#i2snes z$&Z^l#R6&$^Vf6!z5gDMwIl|zf+Z{Ne;o_?$oo$rkl7%@NK{Bf-y^MxTSM47zBWv#yozRNEG^?w862OkF>ex++>N!mNOcC-YP zNT5=||7Z3pf4wa47y)$a4STl3Hp@<3pjR|gR}hX;MHi@eYQXO52#|Bt{TU^&WhLtu zrhlTHD)ANyK$)@6TnolQgn=h@ygs-KrEM1c)@BA1GpV4@Pu%DCZ6sY)hbh6SKCo`H%yC%Rfpy|eh2s%3<q%nsQMW>)h#98u%Fc50)e662Sb>}mYJqxhoS{Oww5b? zzA^^+Ar~={?!N)tLtnH)zim?mP+icIeISJmT%;MF2uoLC93Xs2j-fDKF6l#g5EEmuBxywo{R zcX$@Do@K^y?87tqCS~-8yTR8}h9Z7$tWx!jdMvl2ecW(uVUz||d*tG`hNHfHhO6(3 ziS5_rY1+;0lMTN%mrfO?7~jSKaic&&0SH(@#xzTCXp1!-*b~8zfohobtkgTo#vojI zYintUt2$K9j(B;tEa>>#sFM)^t}G6e$x4ptawKs3N(G&9pu#|y;l@Fg*T9hB5Mw9P zGcsr)R41V|l@K3E6hxo}aPK?@OzSW2FsQ`}2WI{y^%5sNCI)WMfF<^|`&hgugE48< z5#=E19@Jps>1|-lr9?PMVSrx&!9u>q*a_G&uNW|LgRbF>zfL|JH_3w?GGT&NV(QjO zJTIbe^6O7YAr8-q6~-wM*W@4@-}|Jh*IGv(AnQ;RrDhOt62E!$OQih^%EsWtom)4?MVV-OZ8iR7V>9wnAiq!fNy*dV&s@xr&^zCw{&-1pEYWH->|(^Y~K86vI#u;3Q3I6 z#qPmQzco2VHNeonWT>?(rrUVdjcOi>N*Z-I@XGc5DZk01A5RsyZ(Zxrs^ekU*HEgM z);RtYD}!#wLN+z=rEdyzS|){x{qu11wmz&d+8o&L`)h;Rs#l%h>&OA$WUA z?bgt$v4<3nEWqN2$HMyo%YZDjK-~cAdpm$xP4$lia;%LNO0c`mw{?wT+-B^ ztpr~4r>jSKq;7S;b+?!NaRR%@Ls7fd*0aY2NcGsYZyjGg7*pAVS9Zq^oy=D-P1ym$ zJsIvt-3l)q?t7i>&ZxzP^gA_QUH8C`QQqs%l)(_$qlnL+xuCK$7aA?;jTb_flN+0f zGHb&M#-z9fwme7N5OOgvMclp!gd2{fNcdpbBsN!`Q89GdvV{UY zN(vvNeZVn+Tf!c;Tb_Q$vG1kqe~lG*oNJ$4R^CiFf0nuNAYra9hRb`SK)gC{Eq zKvqFGt!iQQS4)k5W{%~<6!>JGHirkagKJ6j@Y}Iwmq6b>3#Y$sxyQikulu_$ljC1{ zPLFrbqfC#3jYpyn`r{@BoA*uaP%^s$!bbhs>wfj{m8R$bAbV>3#@s%|01t<0-;uIP zf(7abm!(0ir8+j{%vJA2OQ+*+J~T4&hU8X^P`}`49+U#ag#f{YSK9u@G`1DfOIXjb zB{omLU-@f)wc*FZ%N!DW-24Idhw?_NpKZo^BDW}8wd~EPZ|N?oql6|P9AGq$b&lf* z)S#fE!{OwEQSLX8 zzoPius;O}jns~(5$zqUE424EVl|eC_NGdd;4_+%1{K;HZ}L{ z5+{0Dxxffidf1_2$*MALFeW@EFSl_+49DI9rop`z2^)jnS&=gcv zC;)t!=A>%1@%8hvzN8o%i#q(53+@3#P>%R*^9Z!%)19DCcp~WMgPj%Kx5;y4h99Rg zfy)UIpyh5f6@WP~aAWKzC6KVBp^xK0Dn;+<1AV7cc)e%i%F*u6SE}}?sqzJdi<4g! z$;@Cz0fS&lE zFk)_ho|e>Qhcat2(zrP3!Pe1(57+NIo^F;{i?AWHNo)yu_*e2E2kE%;%br%ggTu$m z{L0|Nt%qr&HgUpZz4ntwLd~|trbP?7doUwwMfn#v#r1Y1%>EC1QNB2hFO;>ErVFEkxD<71GeW%IjW8do|^H0^C3 zBmI;KQ>{b2Z`bQtx{@(=IY1 zGyPl_4m+U?vXzJDoh=im6uLHX{xw^S9N@-9=B?%^b|`PY?L*lF-L0e8QOL;Y6KMM0 zsJ9GPe~vyA2C{J;yN@lNHXvHQagblR54DiuPkX@!NtPX4Ao*RPO%ot<_qD-RIUNx8 zPGduiQVmUhdtOu$VNTZ7^p?|*wVK#>QXlCei<1`k<_(hp9(D3dc4dN_1=y!0k)e#)K-^V~Y$e=bIggmq_K{8h}+ z%>rc?V591#vH*_>-^Axz;P?YhIBb6JV`4jdI8OZtN^K8~1s=y#SaSEvicxmdkA!7O zjQV2z3YA^H;{En-H{Nf}5On0E?Hl?#Re>Md%HAnyW2bK~$v;1bn?{ldf_IQ9*?=sA@GDf}EV(S$AWR3Sg_uS{MZZ zAF$uD77QDx*pR>9{l?qqw#c6ka3(jfniKUh3nD}^OFAG81e7FEjd_4`77W|op0%-J zE>@s&&NxSB&i#iOv_r@j@ou;--f~se44f)iHzS>h0z_4cfx($MqiJXN|_qS)2BJ>Z)a6yEr@;bvcbz5m^10Yv zw@zj>?Zlt-qBL=u1`U3Wic+5Fz6|z#CSN_Nk-- zlGkIKuR}flUT|Hn`^Ul#)HHl$jx`d|M}|QV9Br20q4wsUCJGS7c@4rA@53%;Y2Ar?2B4T z*~4#KxikXr(xxdP$E3?po#JQ+f$+TgnI$oxX34DsXPRmgU>EtL zj^&^vrI`IQ`bQVP;ROxv#PDyA@yA4M8?haGm3YHGbba-Og?4G)Q6$MaZ!K@lAOzHv z?ub;Y#ID|BNedtUqzzVpLDJwU8sJVUdghDjd1jv$_#db}l10(|dq?Ek)i@}h-V0R(|mLooe}ibrO|>?^BQ5nqN$RkNrQ_ z?pOuU3Ih3QyjJizAC%Kci651s*wumAM!O$=K4a5j2y>otrJdk7WB#7Nic%noPTOV` zes@5I48sKjTg~IsFsA^OtG${9^PjeE$W_BL9)Mj(d~E2B?U7$vI4@xBf$IG?c>J{H zyt3{woG}A}%N50b(*!&e&57sO_oD%h9kdYo^WqKb_klYiP%|3f+aLgu#&g_gMy1ns zmp^jvf5Kg4Z(rHsGm!ek_Xo^~RgD#!!9sNl&?zGS?Vx{LUZyK@0p;DfqhuE-K@DdD zSa0MyNO9X?Odd!S{hx*eZZJTb`+xe;o-)qX7;HqqZdXsgnSx4K0BV%Wy}a4b>|9cI z>6MDrNA?ZkfmH8iv*I6tkr00yAV zeEvCpggbJs!9ebujeG6g9}kT*+RYI*wSRNdTdiJF12+~sb=OvfM!kvg-40%L`}-Kp z6Tk@dKolNw>>K6I1m&HvlF6ojQ<(W`z?%MQ;u6b#3W@34GXlwHwEz5iyW0gj08EbO zFm~wdk?Msj;;-VN;(Y=7{eR8~q6xbwKp2VWce{;FA}%6OiR-Ka?+)HCp!}(`Knp$! zKsx(8IW>5duU88VVk>^$0`2)+~XRj6GuB ze1I~NBAa+;e2JfZQuBbfG4+`D)HMY*DT`M|ZjS53Tqz?l*yAYk@lzH+a~ZHcz4`BMo@AS^vbP^!s%PiD|e#(iX;@lT^24VswN_T10UR)?H zaXOLg(q)wR{kWD8qj&t#`E9^FOIpVF^z~`OT?fdzU&o8I^2`Opvks>O@JL*-t!>AoyEd&dFx_3v@iexXFN zOu_QvLg^?Hn2mMK=qm!uQ-Uff*{6k2OM?O*Wig+hd{l^zQyr43F{I6qI4^0n=8Zh-ePU_XW1kFX~dFy!T38>&Nrkr?BR8ZBsRDxTVKS+>l8LaGmV77Z|J$(xg8RuJ?naMH6 z4`(yO6GJ(?j#n+G#P`<46eyN+0^cP*y&*UM&YL4k!%JR@Vw`A(GMU)}&wEo%yS z=-T@xI>eHQM^ni7%^2wYN$@Hyq#Ux&p7FtZ?c36sD!*RlZ3rV{i~3ylm+lnik4E{2 z?EaY)lVIWZ;(R+NM_F!0dAl!+d*<`TzH$}Bt9m<~238{XINCIm_;O%RsE!wZBP%>; z+~+b38-2`cy`O$1#T<(JBSh`TQfKc?PaN)paLMPhK+4rZb7jnWt)2Uo zzF@4i@lJOg@)vXBDlsXN;^pQj|GxU;@&k)bbB2}o%$NvXHwtKH)4fcTolu%NeC>_)opakQbh#_R9-2olOtx#0O-X z9Di$2!&7YxjpKX{aq5lE?^ld>vbpFpG7l&SF`Dm6iK;pSy*AUd-Q4KGyoX~F1F{D7 z`S|(Ao)JtJOiy^9!_}=bX3p7=1FLRZ#ZQwQjbppAO8qjrWcj)Ef&+ZG!FPU{g%LQOUDh^zeb<&ye0IR~W&(>Y>lR}p za`cQD_@6m6IG1^m?|iF$gOP~|J{l8;s{H;a8J}$}kC%N^Em~cI2T|i2Z`P=mrre!5 z{X1Ji8INaAnX9)!xb7s|%Co#*m(?wGM94vUZcZywe8I46qS&4HHBBo|Hw!tOc5gWo zr4+D@9EHS&*j6f|oi(x8mguPmZyl#&$_H3u+k{?mjGh_MTwj)Lp%XqYTsl-357)e} z4HG_YCXu7V-I%^7|EXqEhUMqV)qJ>*gU0WR&@9W;>~QYmca1-nEs>H&XHLaV}?4r6U(m>i3= z9_j+S@`!v^`}!G?9z#Z)|z#sAmD zl}9yoW^p2GSrTb&NJIf;9F)a8WT%=~(E=?323!ITNu+=v0s{gD2#bsnRKTNDqHGlv zWD^yH1|=j=)WQRcAR$7b3?)$rYa}FL$$SZpGylEw&OP_t@80{~`M&#pzaI;GR5|6D z_%5Go9d*G39-Dx>X2f_2%fZW>`k2FTE7`w z+w;YDrFXRJcP8)Rb?Csx)_|jm5P+kr_{M9Q>uSJ2rtE#y?BHrn!2Cq=q$5g4F!|G& z9}J^9dD}8b7G4FhmwHHJloFBzIgjQ4V^<%T(A|0!Xg|+fLi%-f-S1z&ewO)0_0pg`}fFS<$ z5Z%QeI$pQ^kAh1BgFl$jy$LXP2bkyyCJG&xtXRW9230|m|7k<5)x zcfoSE9H7&Nt;B%OwmuNv+{J-VR7yPN@69- z?onU8K$4=4tc?+NC$3baK|;@^k>jKw!q1HRYW43i@|iQi_AQSunmRXhPieDZ1vsq^ z8rtzEPX30<-x);wuvg;I<-PLMJg_y?XWz`CIuA@$_qAH53=ZLPB(zkLfxrsE8Lr6u zh;7%pDxrCO&q8WZ(5Yt8ADn`qI;VM!OLIoMvH^Q4<>Y18c!$blLZI|p#@IM9zq6bi z*^dzQ8CC)~Pg*e2>S)>3*sppN=gq5q8~04{<~0Y{A&K|L4?>WiurAER^xp8!4d*MG z6tCpxQjgu}1lV!KmU~pce#sDQtF&D! zSJYZ~<@b3gxU}IILYxtANtl!9w^GnPXNyJz^VH&o)cqe{02F(C`TSq%Hm4|qSJTA}E0MrYt zN4qGP*eriC_XRKhF8x+gsw^&3SU^B^J895UUAo%DUC(NzSsPdX)uH@ZuG&3r@`RR+ z8F952sm=x|h(n>9dXF4#eoLPw^<1_RaB3vl%Vl=@!fxo&*mG()re~DS7+Z0Lsm>1g z?SoK}qf$E}9SpCS+W@j<8i+|qH(0Hib+ti7m)f^5#Bo6OKz5Yb9us80nGsjadU@|l z1U!RTGfU?}h58z4YMF!>uQ1d>IE!g$>$d$=Bcv~Mu?rW0Wo`uO2~XK6`$%B=hini3 z&0PbaJ61GS5URj^?{R^SfNgj1HNf=$6u~#2ZejZN>JqlSGyH7I5#e%*%hS9pWBc|- za3cSC+r8ahvP~~M5(pT4kWCy^t5;GW5uvYkQ5@~MgPj*Ax~KyMZnp%zb-0adX4w_A zh6=W*@I`JSAonFbVy$ut`M+#c#c>S{qk%%jlss06DJPQgr|z(e&fZ}^e<@oGSwPKc zX92%pV5J8mzYXh&0^W6)mt5O-!K%@~yA-OiliGi$cHsuw2; zy=!pMO}^w@4pbDOPy^^S-gje zE|Zl!0bs+3+^{Y$u+3Ynqgn(w>Hx}fPgkmNTFUe%WB@6WQ@3A2lU-pd?^?_|Iu}&;+^@99Sb~xH zdVIQHtv}%%kjPN=0!XqkNHWmN(2u_DE>$;p*V=22Ip&y)UCv$6kCv#hyr;f(<&c-%|#!rpR z^-Uf0t&5YS?vRk4B7w!&E(}d57sM7V{IiN9BI%a*ur5^ie-tKI;!mPmUGCL*{E0jhd5vOP{ot zu)A~Vwpu=A#FcP4zz_}dGNnh2q6Pn7J|EnXBP6X7nawAxmX2eY?=H_AXyKgV_72#pN+CvL6T}oFM&*=!`4H8+7=MzF1-sZ}oihZCu zn1zdhgI7EE*XiQqbVlIAOUGQ4?CxIfWO(u1Et8^a9@bXreH&~?&Sb5&M~3ceW+cCR zdqc)GvW&(kgLNfHCd@_(gir{cli}8`twj1A#U3|&9u(Ue{eF{boK*X30!?=1$qI`P zt70JR%N=NOcCQB}@Kv|QU31*-2JYTnt#spSRky0YvcuF;q6CUR^wmk+h}z~I6aF

H1ZuRNC(Q_qmNXCx7bfH$tJx6NjwesCLo^mZlGj8ZLWE3&`f3Cpyfh zPah`ofPO-A3)ih9 zI>?7^%FuVmtUcvO|0G1Hd3{$Ioj6UuiWO{G*QQW*qx3wx;_;+9vpg3t!DvXU<#;up zsJd=nCHSlk?<4G1w)o1+`)p6@eCwJ*)%!%Eddy2>&&AH0Ue=p+seDYgG)tpVxA9XE zn#*kY%9GlPmKLsdN{TJsPwu#t?H<;@s??a!(qN0%G*^{w!Q^vn^+#Uv-o{9uoovZ zmjr&75qH0DBz=v38j8lbYhg1KvB8{necNdYjPBx^8p0(;TQ{ad_LsJ2*!Au*-N~A5 z4oX@Jk85;antDTY779FH~ee-z24*_pDmi7P3b>e_Ks7T@;e9j|rx580#@ z{>^lyi_=k#AI}Wngcp|^VU%`l!JQr6w`a9Dq@1DWd-BGMZTdIY3k$@e6D4y@O|bdr zazcB|f=>L+F|XB-y&M31`*=sG*P zU_?GO9%Rjd`-A?uJj-FO(!A8;UbSD9O1FDNpw^PThM6D&HcxmWE(Kh=_lL9xean>z z9fQtDuy^>BKOx8duZe{c=KLSjxZi8-o=B#^NIp7oY>_dA5{(?ifYfX6L^aRTx#yA! zIqDJ!T(y;NCK*`b9WKnTD{LbIn?zvQ%g`jM3S0`GJB~)@$sD6ejyMPP1#k1Q>Sv{c z>|`}t{$~l-`xK}AXB^y?FiV$PNAJ{>b6V@h+u)sphA4A*av$=Xi~cHAyv2!4EPd_|3jW^?EdxCatO z%G7s_6KM;hh0z-0?OjYKs@dcj2Aq|#+FZNMZzi)>tHc^G0gF*l?}#gGTe-!C6Pm>- zvXxq#GblW;9@}(i@5X!6vHMLE+`il40tV9^hE=fdy3W$n!SX343)Q+!4F^+$oBi1u zviGU1d^UNx4=Ql4PZ{%*GCi&)ZZuBG?Vk1vy@NJ7u4rrAO+k^0;0JxV2OWd#gTdV1`ykt=;tG z5WVx{<-x7SX-hhlCk>hPgnX;ky9h_vk3CEwu9Lkea`UO{31QI_CP6&2uHp!-{$9t_ zwQv57_M>&ZwRz<~*yqz+P$^#SG)GO;ia2Yod-F%X;4tsU?3k{!o6_v;J=q9@8w1d{ zc6ucmJ>SfbmMpb14lNs%O9=^Du!n|QO&&V7??fNSC!S&zYtnRlA%GpVW4|&VE8z3J zo&jlo&VLss2q2VXB9G7~<6aTe=iC-VRK8xN3tx$Khcd)O6ZsB)7M2dkUbGyH?^z+B zcgDt()|s?tFxT}!>~bjxPJ8;4*K#}zxLJNvBN^l=qTQJKCYoNmhSX(iWH2PR(SAiv z26(*nG-dDa$RF-Uoe5Po-Z&jWa8N$#9yXpp7LDibaKMC45~{e_Y~3pg-kW4Q;l;1S zECQ?FMkfs;Qk%{{cHi*W%q9ppt~*ZZn^v3*1ewS)-kgrAMJsxB@2@=n@eb7HICz`V z8KP^{9ksci6Cs+4)wbNpfgXOnQ`)*VWqPn&4%N{cEATy;b%b%+FZXKV@cG`+t`DaD zlIEBjT)2B;ec`&`P44G9e<^un$2WWUbgqG_Z+9wUwAmA*q|ruZ8S4tdSGQLq?S9fP zx%Z6nTCA%>E`cYtRaY?bGD1S^{969z_qF`p?^Dkw_GjK3Szi5lLF`L`Z}7R+!fWfi zZ+JQ*@M%C_D#LGCuo-UKBnRJ#=F!Lo-+GtO_Tjd=K5Q0D`Q2{Xx{}B4P3AIzhPmRq z9t(A|`%_7|P%L^fA3^>8_G&iY$;^`!cRP1}Z_Hz`hV|0v!x0?KRp(c+8mLb?qN!IF zrmKB(QK2DQli+|EXc;K+oti})u+ajSh+@%_E6Etvpd_dcxN$ti-FAfhFuAwEv6z$jf_KP~w&RJ!?1w|n zQjs%uzPsktS7l!39&v*@gGqy4ec2PSO@n!@x?PU=1NuihmKuU9xuH3i32M|U{*i-+ zOCiH9RPGC%vrgRsO?C@2vsb^%evx0g5D(6uaZr06>8Qf5u~S=jtW$mH%HNfEUOX;o zU3`<;yavquI5F#NdL@+Og27eyhI?b)$*KL=3!nu8&3E$qbJEp2ZWo8GCp|3qVj-Aqom7f?4$DJg)M_m*VN1Hzet%nV3EXHJ1ixazNY(@LOJ>w&S?RUb_ z>W+Ja>Nipg69i5O-i0-^8r^E!aX=|YV~JUGOki`aazdbH4`*g-+hfls37q=m5u!eX zX}))Ykq^Bt*JN659O!MEtJlg3)!<1JIuyNP__9W8J*B{CAE-Td>p#!;7nU?cHy*Tn zC41HN(&7&)5^!7S=hoLfOHw&obX^gxk%OXeQ%Ia^Y^qs>L1%C0T%N>S8sEGd0V;-0 zM|M13x>(TdE9waiQz}W?#~Ed1-|lWbQ3-f$W~(~wu%EE5*UV7sA1%vIHVL@JcL^PM z{o1S5MSqf>UH5>ZVqF;t*G0deTm7?!l8c5VdH+o{KU86TYsyS6qO*! zV8%i0FuuW@;MrH2G}uqpclb$K+nxiO7l^(5Al-JbvhH%}*Hf1hU7AnA{P|Ea;k$O! zo9we^!aQgCfvG;SSNM`CQ4vrf?2Vy~Cb7(jiQ@QkXIqZ%{Vz>j-V5wfq2V)B4PM@I z&~b1?=wkUU6)6Rxl-Z9CW)GYccTi`^`Di~Ezc>ht;uv-n{od;g)A^a1eNgFsk`o#? z)%v9AC+ZJ+iJ8&~jalyXg3lYaO-4A*+kYSrZ(PNeZFZZ(uXu#4u01`Rg(H_X<}{8I zGw`4Io^N;`qQi8EwM7`5g`zcSq_sp$X*`)CA3E9kQclde>THu=y0J6pj=G1CBa=?< z7<_rHzr>L|?#(H#x3Zolc2@81GagMjBha9hDIY(Yz~w0HJAeJ{m!8)W>Ysc~lHbzm zc%`>HX!EI$_Pjx-Z7+Ijmtad;B{be>A#v6ke6_gJMN_vhbLinrV1@bUJCnYI0gYTd z*-J8hkG8lpe^-)$wFXc2jURNu2aYRiKgirPY}r8@AMwmR)!&VJon+&8dbhSU()2HV z4iU^QimbfSdORnlaDjarL7jZ;cs{=@knH8WJrN5oR-)^}wAyqc_mw?9Eb0Bc_SN|ejiQ1b2aQi%jBFAw zhA&B<=KBch4HSh^o$jmawRlsHPEtq0sh<0o_9x%@=RlsRB-c%UWQBa+bj*3fRAMW<3U7B; z_UDAkdfCJ3ulJ8wb??w_uW7XNS6Z*O8|!x~kezlX;{d|j+tkY~im%ebav&7aE8iF_ zcMZ!r>;B1r*94EC`8`;g9KLhY^k>Dq-jhN^w>lZ%7+tQSH%Y4f0en*P9&lgtYRTmj zn$|m-xJSFO)kPxP-)?7P&YsY~Q@xG}ec#bz6fGPQDw&G*J65KCnR%KIEqirfl%#8| zwC%vi%F60va~rs*Q`7!lf6+01)(*Xcy$Mo+ab^T#&z|Vdz<{mJvdCx7-TiEQ57%oJ+~BtQx`fPwMSBKjR2v|T7rZ&YrS4)r8`=ofhJGCs$?;9Y?&ZFf5b;=&JpRq$ z>dor59aKji5_9FQq+F_#Vm9r>fJhg53s?OU$KGT|G@&E(r32qK zt*nthHkL(qVYzWNURO%+3$OKbOlFXsNe_%tfa@wXjKWo_(T0MX(SA~9gJqsXw<+}l zUbj)**Hh!Q3324TdQ+u$;Z5Hf_n8G;a%|=rysoV~8<()Lv8N0pv@bV*7$X&FUn^PN zIA82HM$PfDP1BR(5qhB(>G z1en{($qWamMcDePE8W0Wgcyh;_ZrFc)qwaO3B^Or{67boxJ591w}GXDOWn_x0%{DKs4Q9ZfLL&F|Iq z;5;clo_L>o>q!fk_MucN#lx%7C5*)*-MH6e)@`M8*6U<)c?EfAtgttcXCb$sjtecP zz5Su*L9-ALF~Z0DBj8N?xHqs?=`P=21bc+x)hu_RRa4R3(Yl7CIkH+_NB!m<%EN4% zpn%*`GVCWex>l1FDev7cmhy8OZRUPmVw1DGPqX~VS4@zI^vVPmq-B|p7YQ{U%%{vZ z2XPvRcgOF`@@e-i1*NGhnb#e*x14S}Bg=z$!l*%b8S-h*kU~lApEY~pGg0qX^r9ah z_d1#_+$rmLJF}~4*$1X=BPFDoMLD=9+u?$OYSx3}B^_6PGC2I{9G@qRiaK&iSuj`9L;H~SkGYMB+;s}K!r}+}+1$&dw zMepuV%T1LJRdZu%Q>SA;zGYi#kZtd`$Xthk?)Z^v?{0fdUkWdZ75e#(77+OeYd*$) zTp@HQ+qS!0tAnEdHk{vjVUt_uWnhP# zw0G0gA-1jQc4%Il0UEp3L<^@UCaOZGy-KSzwx@+L7bTO#yP6 zSN-Q(xjVMe4VmR-r4m_za2B~%$kFKjTKX%B%7grwz+*=ur_AmcKI&Xos1IfY)o49U zwC>=p+q&gIeabMrE;vAL4uT5Xdf(ijNi52@A4o#xFqePYys8imYSLZg55%J!$jIIA zH;g9Sq5o~uxcl>~z|K_tQ3+ETjs4c1@Xab{cP{&xEI#Wcj+WQCogmdhlqAEPSo1Z5K$rX7m?<>l8L3p{HUA}aCGii##vw0uHNMD&e?@xlg-~2!+;P9w+mMabrPUQHl z*dBv%l#kdOS+7zJq3;$>zKO|66IlJI?PaHU-{_Ld>9V`(r$0);W5H?IC1Vr;@}N&v zcbL#tvpf8mi;ky8_KZ_w(RjZp6AURai#m&cRnX0nSaj(_Y#P4ZdE(lnxu?dhQSN{# z-^v}wGFbnz{{rZ`xMKk&mObvd5$O5El1QZZdxay^q)U_I==T|iheULnzUA&$dJTB< zw=MZDx3O=WvV6N$N*kd&gDi=8AAyzmB6?SY$)|H+7uE!4!^dA4dG)UTX^c=>zEV=z z(gr*#EnYsU8E-W(WWV-p-*)~xLa?sIe#=l6W~ z!|Jj9HQnf8-_j%Xd3UaVT5^xw`{ zH3p-cnK|TN>}XV^((P15xu2K+55xRtY%E(w$`2AN)_7W|ta<+KUyUx*-2?MWw!(ke zQ25ieeeh|rr~(GksE|{UTk)?fWHKur99B@T$`~e5Pm$<#0n!|B6!T*zE=*k^6i36K zhHfKe{r;V89~Xbu*-N0Hs5G>V5v$8pF;r^LuG|`|z;VndV2Si?!kHk~6TdF=d5&lp zbK-<6$~upGGDXz7E7j0*^NWFU*rZmE-L1j16V*kxMa#tjBI`JF&s=i$F>d5!0?&i~ zi-6oHTI1MHAB2W z<{$?c!Z;q}#6Kg-@hHu7Z{~M_|KiRl!yHbXaBB^Nj9x5386acnOC_PQgC>im<@&r? zW|cBCLUVTAqgA|50y{nqDzX=+(d{S`V_B){RvsK-+>qKS6Ub8R$}@UC;KcvsT!tN~ z0aBHkt@TkprytjAzFS8F4c&qCqfV7XpJc&#{51} z%=OS5DNVs)1l+Lnii|y4K|S@H-J1FwEeA5!Xp zclday&m(Tfv|CW+jHwJW#+v2(kBI%p`D&&69BU7Q#w)cDZzD?!6GG#~Nyp{9z4`UC z7&yv8IaT(liWDR3-NwnF+AEVj@f_JCmC3RtywW1AdXq*Z3=FX#%2 z=FBuIR!7~%0=3|`J+YR8$b^1(2f~M9i-#D=l__k@0FIz{_r5w|(o+2VQ&rrR0#=($ z2M%m^=oy4*Dreva*lP;ZL<8uf&VEkB%+Q_=V4o3Vkr`v3t2vXU^MTL#X=uCN2ax6f z08aXsb=T@sR%<;xy2g*uS85!rxYsfXc>Rl%#PnkwVdWL|!o50!yonBq{<2SZ{wS!g z?J_Cvg-1rws@RWX^}Cws=q$=&Y18*%w}79i>Jk1%OtvH5%m-*)W8m1NhfLTdy&gAr ztG`w+hzz5gc@5aYx)GnZkjX_X1^j#}t!x6s0UXG_h^pMtMd zY*shZCxkg(F+=t#F!8vM#)w#@gdqo2o8x2s3OQWJ*b&`!ZJrE7F^;7S(B*Pbd{z;aW!3vJfONCFwt@Z`YKWt1x4tv8+%uas%6$C-Ocp#K$=>5H>`<0QYd| zPX!bE5|k0R*>gffh?kA`V?yqGF5OvRC4JFN!K(i~LMDnPQt!E8UDicBK3(HtoD!&} zo%GtlQ3he9GGPKDkOEPrA{X(-V^*i2_u;4AalRE1jPkc9KOnBFCdcP~z(dScp!*C! zm)1nu*f*=uzQy=n4il(5*HbaHDnhO_(l+7XC?Y8|UjMi3Q9orJZdo;Qdls=(5CI8@ zBP@}kp&tY~iII6x4v?%1FPkxFLE{f5cMT$G!SpwY*ppBcfOmTq$qzpnjiPRBcw&oI zT>JOdUK%510W3pg>uVfrTS+||*55etsOdf6&BS98e%%0rZT;PWX_=&NOh83b_HL*YIG1 zUG_RF5g#BJ?WGQ`sJW*6yGm`<<=4&8>pLo$jTuQuAtQaB56oyt^6tp7!3m)(^S1u( z_VCg$h^mzP!Nv88e;3(6M{L`Hq5fWJIc(t@_nlLz@ALC3`tsKl++q3a9RLaee(#{K zno#V2?MK*qn23lYC)NJH42h+&O|T7C5`oxW(j{YeTIu-pP zrk!$I9Sji|l!VjL^ch*Lsim+X#^wR(Tklix5nKVlZb>4esZ$x0a}MZ^2~zH(J!4Aw zWwH)=D5QGz=@h4ge(jez+r-SpfG)J*XoJK@bFviAj@;gV+4s<t($3=LUh!Wlf;C2xei$*@ZvfmDw`dYHq6AgQFm zbXY$(k;$P5D!q8G81Qe;JW-!U@s` zV{O^}s!*p5aCBF?kRtVzlG9R7^?{U`rN zFCkvBX}wQe(P+^(k`YM4My666qMVG`4C!s*C3V{Nsf1K@D8U70>dNKn|Ejac6Te`^Ur0Enp}l?^E%y{dl-F$wa27< z{;}SnUk$hzrLg%f9D!v63Y+nH6{-m#uLke2CGJFX!VoJ~?09Di{XTrOhy(YrL*9A< zY&)?1I1u*O9EbN>rC6HIBO|8Lnr79)ov!~JbL%D|HmA0kBnhW{F*L?NI%8SpCpHY&!%*GJ#q=FEIQ4lNif|38y#5)S)kl2Yr4 zSn=;kGBrvV)P}V=EQXSM901uvEcZYU>AR>^ETwVB;xiIGZ{l!Y@3x9=j4PBv7X9RU z9YO?4t>>iMsqd*qTQ`|>5$X!m#{X}q{(i)+dj){eUBPsTk;4>!w6=Aet^^xl1jp-v zFzsl7dAmL-(X#)FU1wUHp1hq0@7-9i>hU5!$<}ys1m{EUVgSK4^>juK5zrj|YR(O? zT~G!6eu=B)M1}&)IJD?_ku=e~Fs+^*x&3rPW=HYDoyFt+0-6_jvdkORVp?9MpdNRm zj>xURUn}`LmUtB%o@)mXWpK8JBD2pk#3VZ@Mk2brSFE=3mR2OKVrZZcAgl~drL5rL zB)_2BuipSSQ+N;#adw=odu2AJBiG2kvTCr8GZv1@s$<(2u~8%?o_s*&ea?q8>aQJ@ zid{i6tkdb-hV*z4#pYq(Q;$sW4gXaaHH-Dl?U!@$@QAqt)~f zm%NAF11T7R<~mucLxy{*DPDSRjz0;t0KNDVv5o4c}ySuwTA6Rj1$`{_;I0itYVMhjH*Nhqbc!z7! z6XSV|8-&Z)sF5OSe>l;a+mRiv|k zz|%Gc(f=7Wt$3Mk-kwHtXPbzWT094<6yDKSPP6cr5E=(>c`< z40*0G4b&{3mQk!$qIE(%^mWrYn{J{};5j!4az~j#J0?RavrT3tas2W#5CZ_G&3_W~ z2EeGqrJOukm7x8(!BiX47$7sx`pCZj{byf$`phmsI0fBhGH6u>x5w$fA+q*_=fJI5 zZo@A} z&(t+Rtd7TuG~EfReZk51q$ZG|S(mn>h9z;{`E^JZR0uln=s3FqX;W%DH>Q*<;0-q7 z!~^M~~| z-jX-&p$I3Q&c04((7Yp3C)6__gGQ?LaLMrDfu$NN`duQYDxnx6Q!U=Doc$p&3NtEN zjj+y71FKvu6X=N3V;b6SoHTiz%>B5dDohD5DL*m4y9g#L~_m-CtgbK0Jd?KEarf>5Fjkgc_5S8@qKN&G%n?9UM} zc8ZJCMc+K3zSxJ8vxTxqzL=ywIrsgi6CK67Q`)$n;+{KC;(u6kQAs)^DW5yV*p7DS z;zaMv4Po`aggCAXQ(vsU!AKzou;jVG=@je1jp@nF8nj)xRGn0rsrRJWUKCxVe!

    ps`lCAWu@N@XwJelfHk{e#cu*+@ z4vs^5$gyVCr^cCg$-Xcq?K_86i+0!NKgTr$wS|7=6pPrA9zq}=z0X)xUT`>)oGUP1 z(j}6jm-=p2Kx?hb_GqkS5D~RNB7su?N5mB-Ttea18pZ@JAfi@c98+fQ8~So`AhH(b zs#||{6~umR7=h>Qw!iPb+hnLTqmWT#iipIs+IPAj$#gI1yDnBV-mKq0aEQYkKtq8#!_{o!9R}@QqxNt_eUlETu6M&>750 z3{tFJlwGuBiZYXFG&kmd#FJNNB76DGAa9(gIQ`D#c1F#RDxO3ey{YIgB{;Cf^%-hy-jQ7_Gyp zl%r6;%|VZISf`^Cc?B;K1Cor0%x84WOnqs8oCevRPgj1Y!B%%%NxdfDOkU=Z@3hGd z6Wu&=Dji%EX}g##uCmAIM_4ZdiCXrOJndU~o7ot3JwtoHbu zY%|0H6o;sG`U=Y*?>LbAdyb}2k~mf~6ncu4zlefgfQg67t!HdB#|mNgXzdw@Md5QG za#?v7zWSp(`lOct58DbPmJtJ5$k-Q&%fhm(IT={P4ANO=-_79UxW@_(VuQi2v}0EM z%IbiOgV@jAcIwVIDQ}sC++&d-?+RCZ3_&a`9ZRXQN$a{kP^!82-qqigllTz6#1im> zde`uS8Lys5auq<_nt%8Yd#G01@N;N2IO!bx3}}TLjTOkXJxGTaX+)so`Tl;EJCg}` zuBTE4QL0A`W(*ZV$!p_E_e~LWwn_)2%2b8t;HVx_bZHcMZ880U&u^LH>%x9e7Mdl! z)}TNSnN%`kjJd&!Eizq8&zWnKAEd>3fBfq6ra_N5E0t+ zxs|WVsO*UCOuvpJe0hrp)jySX|1op%m{J%+?2JEy|5f=n;)8khFV>i|l$^5NP%BYX zN}9$647ZXPb7?;Lt*%&`Y1%jbvX8D5Hf!>5AkP#?E~w2YHIGEJ0w;7Zh*0DsCzU3t zzhe}k-zA?W`jjIj20N{Q*EZ|yt|GEF1d(i!(@^`ll5>K~?NKPBxd$rbG#l)*lX!PD z6xlz&?Uo@$(M6H}wI6|Tr#x^r(p3kii9qWrX6CmxrjqxP*DaX%6UJrCxDWKNz- z#5G?k0U5ral@i5{z_KZSQ?Jc1HXPH-p;V{U<(74V@7WTuY!-RG$&NtxAg@2Lf3v5* zW5t?TCO!SNW(ATK_-TlhorndY4#?i+#vzq-zD*88NTh|YG5|qDhtR-R^fx~xAiN@S zWdTz%P9(EfcYbbSX62?n?F+FC7OJ&Jy&3{!c-ab!WteZNvF7IJ0R+iLq(cC{%>MpR zj31;pNi_KYaATS0<~C?x7QP_o$!$s(rlGXW!G~uULRiFt0L0{VJZIJY7V-o{bFr%9 zNkxp{J84bTB8T1DaIFOoRbESUE!k6t&w$7uqWh7A-!5byx5(6Zk^CnHntpK456u)W+@>;UvP-3s5m0f9!Kr3!WZU!jKp5ieOeYH7hGu9cq02{N00< zPj22vWocnUp3OcETWD>jIMm{4fnE7+bqDh17r9s9zm3Y8CdCrl|F~Fc56|6M`-*hy zSKVg%dM0b*c^g8BSOJuX^6OTqBkD+;E3tVlJl1g8gGl7jDgkryCn$}- z_W6Qwqiv=Y`<3Wduo8K}D%F^~&laLGhN|V#y$z5}HAecf{tYYdCjSL1_A3Lu=`a5e zup;Vxieu-s>qGO{^MjJCW`nyFUy0KwJ%Id=of;=z4{d@I6E1Vq66do4WTSHfYNL*>U;R~~Se;RA57w_v(PRN@&4;(%-$gX=*P?sy%LT?askREulc28|zn&W{^J zKn_Mp;b%Lg@uOd`(`*mflt~QM(IYk#AN*Zur-{={EExc@Voj%MV{<)!&Ow2FXP$x| z3aHK{p=?PRwEKAd80$(7$xG)zZS?^+t8u+4SL*W;p$}S+b)-FP9@!BLM+@{Ezj~?Uhw5+oD!In#gowN3SWeTK*^Vcx ztn!D#vlS1%ep3OQ=e_Z*lozWujZ_3V?o7@)+VN6Yqk@)G6hEOvygQqKN~rE8qkn!m zRpAt(a=+K}img9=?Cb@~^iNSn{)rT{Po9DUN8r|0koB+Y;WJVZn_ zoyB4t{W4-szh8wR@1vNjK6RzO1w5}BCjb!Is`-h?(dqINr#>L!V;GPNSkMxWr!L?` z<3N_|IESSZz2BYokVN`vzL~Ye8o}8p!_%bdFc9Lm`W+P8-qkycX1xG05`2;7b7ip@ z?Rzz z`R@RA3|LM+3O9h1cpnc!bi zv-#SwI3aNPKdLPP!T+n2OQgc^c{CZL&uN%RQhgAiBL7PkR<`+nDhuDgSzsOrhW%G7 zX6oskMO5h#6|BE=q=hK|XOcRp@BZrP=5!~!h`06cN#uwlTB?kDaze64h2T zsbYb!5l3M_rqQMTF6GyIerGDDB0DTY^E-ugxZ|3+OF&hF7EFA?R_49yqphbJh`frU z$5SR?Ubo(6pvPusH0$OQtDg?0wk?==r}!Y-r+i`#i;@Pr(KUyZ9i<5MV^%_q^2yVL zLm-M8G)kT{|k@b_veo(L67_Jn?7xfdahSC}V>25`7 zrE8EDj$J@G842Yuqb^cXDfmZ5H;tt5$2NrFV^9wqIzv&h{cOMywqLM4D0(05j)d_6 zr#Ys*{m#od>&j!6KCS;mER%qx_Ux_KC(jARDU$CfII>vE`FICC+(bF#w#JZldA}cA zr0izYz;Zf_lf)$3d-<^P$99H`yIwwV7i4$qoO#II$4NvVQCl@(LXIW+qCvXOdbfTL zwTp%uUojg|6#_~Sj!(AUpVrlfiM!6!^mEn`@9Eq==BL{~3KWl!*Gf9QY5_03++;9U z4C5QK>?Opy?p$}gKbUb50NALKeT?+s-!+EY+kd_7`xxHNP%$d!>xOA&qV(tZ?cb+W zTDTTdc{|0w2mVOGmr0VpN=|$0J>pfes9Ed<;OB_v$Wdvh+L@)G|Hj^uVi=UjrM!Av zji{9v<3l>(UU!!kU7-oTtIR90D0pX=fRn-#I4E%%xyDKiS#Ml&LgWOh^i4p$9&tT6 zi#I3o%y*_U*1odw?7kVEFg*7{HQ!OI&aWKoR)02K_r^caA{XpBq9&dZ}E66D zVO5HRbexx*=XlaKA~ImTdO*RM;A~cI+l9EDsy35 zud2>~!eI-_M`YAKf&@%idHuNtc5dK6oFESU8NlJTNZ&e&<3Zj=Hvw6hB`LS$svaZv zl5|wobA+G`ism#6|NT?JtR1jbVmX|wZ!_1Vj{r#AWG{f}_}h + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. +--> + @@ -29,86 +29,336 @@

    - The estimation package provides classes to fit some non-linear model - to available observations depending on it. These problems are commonly - called estimation problems. + The estimation package provides classes to fit some non-linear + model to available observations depending on it. These + problems are commonly called estimation problems.

    - The estimation problems considered here are parametric problems where - a user-provided model depends on initially unknown scalar parameters and - several measurements made on values that depend on the model are available. - As examples, one can consider the center and radius of a circle given - points approximately lying on a ring, or a satellite orbit given range, - range-rate and angular measurements from various ground stations. + The estimation problems considered here are parametric + problems where a user-provided model depends on initially + unknown scalar parameters and several measurements made on + values that depend on the model are available. As examples, + one can consider the center and radius of a circle given + points approximately lying on a ring, or a satellite orbit + given range, range-rate and angular measurements from various + ground stations.

    - One important class of estimation problems is weighted least squares problems. - They basically consist in finding the values for some parameters pk - such that a cost function J = - - sum(wi ri2) - is minimized. The various ri terms represent the deviation - ri = mesi - modi between the measurements and - the parameterized models. The wi factors are the measurements weights, - they are often chosen either all equal to 1.0 or proportional to the inverse of - the variance of the measurement type. The solver adjusts the values of the - estimated parameters pk which are not bound (i.e. the free parameters). - It does not touch the parameters which have been put in a bound state by the user. + One important class of estimation problems is weighted least + squares problems. They basically consist in finding the values + for some parameters pk such that a cost function + J = sum(wiri2) is minimized. + The various ri terms represent the deviation + ri = mesi - modi + between the measurements and the parameterized models. The + wi factors are the measurements weights, they are often + chosen either all equal to 1.0 or proportional to the inverse of the + variance of the measurement type. The solver adjusts the values of + the estimated parameters pk which are not bound (i.e. the + free parameters). It does not touch the parameters which have been + put in a bound state by the user.

    - The aim of this package is similar to the aim of the optimization package, but the - algorithms are entirely differents as: + The aim of this package is similar to the aim of the + optimization package, but the algorithms are entirely + different as:

    • - they need the partial derivatives of the measurements - with respect to the free parameters + they need the partial derivatives of the measurements with + respect to the free parameters
    • - they are residuals based instead of generic cost functions based + they are residuals based instead of generic cost functions + based

    -
    - - The - org.apache.commons.math.estimation.EstimationProblem interface is a very - simple container packing together parameters and measurements. - - + +

    - The - org.apache.commons.math.estimation.EstimatedParameter class to represent each - estimated parameter. The parameters are set up with a guessed value which will be - adjusted by the solver until a best fit is achieved. It is possible to change which - parameters are modified and which are preserved thanks to a bound property. Such - settings are often needed by expert users to handle contingency cases with very - low observability. + The problem modeling is the most important part for the + user. Understanding it is the key to proper use of the + package. One interface and two classes are provided for this + purpose: + EstimationProblem, + EstimatedParameter and + WeightedMeasurement.

    -
    -

    - The user extends the - org.apache.commons.math.estimation.WeightedMeasurement abstract class to define its - own measurements. Each measurement types should have its own implementing class, for - example in the satellite example above , the user should define three classes, one - for range measurements, one for range-rates measurements and one for angular measurements. - Each measurement would correspond to an instance of the appropriate class, set up with - the date, a reference to the ground station, the weight and the measured value. -

    -
    - -

    - The package provides two common - org.apache.commons.math.estimation.Estimator implementations to solve weighted - least squares problems. The first one is based on the - Gauss-Newton method. - The second one is based on the - Levenberg-Marquardt - method. The first one is best suited when a good approximation of the parameters is known while the second one - is more robust and can handle starting points far from the solution. + Consider the following example problem: we want to determine the + linear trajectory of a sailing ship by performing angular and + distance measurements from an observing spot on the shore. The + problem model is represented by two equations:

    -
    +

    + x(t) = x0+(t-t0)vx0
    + y(t) = y0+(t-t0)vy0 +

    +

    + These two equations depend on four parameters (x0, y0, + vx0 and vy0). We want to determine these four parameters. +

    +

    + Assuming the observing spot is located at the origin of the coordinates + system and that the angular measurements correspond to the angle between + the x axis and the line of sight, the theoretical values of the angular + measurements at ti and of the distance measurements at + tj are modeled as follows: +

    +

    + anglei,theo = atan2(y(ti), x(ti))
    + distancej,theo = sqrt(x(tj)2+y(tj)2) +

    +

    + The real observations generate a set of measurements values anglei,meas + and distancej,meas. +

    +

    + The following class diagram shows one way to solve this problem using the + estimation package. The grey elements are already provided by the package + whereas the purple elements are developed by the user. +

    + +

    + The TrajectoryDeterminationProblem class holds the linear model + equations x(t) and y(t). It delegate storage of the four parameters x0, + y0, vx0 and vy0 and of the various measurements + anglei,meas and distancej,meas to its base class + SimpleEstimationProblem. Since the theoretical values of the measurements + anglei,theo and distancej,theo depend on the linear model, + the two classes AngularMeasurement and DistanceMeasurement + are implemented as internal classes, thus having access to the equations of the + linear model and to the parameters. +

    +

    + Here are the various parts of the TrajectoryDeterminationProblem.java + source file. This example, with an additional main method is + available here. +

    +
    First, the general setup of the class: declarations, fields, constructor, setters and getters: + +public class TrajectoryDeterminationProblem extends SimpleEstimationProblem { + public TrajectoryDeterminationProblem(double t0, + double x0Guess, double y0Guess, + double vx0Guess, double vy0Guess) { + this.t0 = t0; + x0 = new EstimatedParameter( "x0", x0Guess); + y0 = new EstimatedParameter( "y0", y0Guess); + vx0 = new EstimatedParameter("vx0", vx0Guess); + vy0 = new EstimatedParameter("vy0", vy0Guess); + + // inform the base class about the parameters + addParameter(x0); + addParameter(y0); + addParameter(vx0); + addParameter(vy0); + + } + + public double getX0() { + return x0.getEstimate(); + } + + public double getY0() { + return y0.getEstimate(); + } + + public double getVx0() { + return vx0.getEstimate(); + } + + public double getVy0() { + return vy0.getEstimate(); + } + + public void addAngularMeasurement(double wi, double ti, double ai) { + // let the base class handle the measurement + addMeasurement(new AngularMeasurement(wi, ti, ai)); + } + + public void addDistanceMeasurement(double wi, double ti, double di) { + // let the base class handle the measurement + addMeasurement(new DistanceMeasurement(wi, ti, di)); + } + + public double x(double t) { + return x0.getEstimate() + (t - t0) * vx0.getEstimate(); + } + + public double y(double t) { + return y0.getEstimate() + (t - t0) * vy0.getEstimate(); + } + + // measurements internal classes go here + + private double t0; + private EstimatedParameter x0; + private EstimatedParameter y0; + private EstimatedParameter vx0; + private EstimatedParameter vy0; + +} + +
    +
    The two specialized measurements class are simple internal classes that + implement the equation for their respective measurement type, using the + enclosing class to get the parameters references and the linear models x(t) + and y(t). The serialVersionUID static fields are present because + the WeightedMeasurement class implements the + Serializable interface. + + private class AngularMeasurement extends WeightedMeasurement { + + public AngularMeasurement(double weight, double t, double angle) { + super(weight, angle); + this.t = t; + } + + public double getTheoreticalValue() { + return Math.atan2(y(t), x(t)); + } + + public double getPartial(EstimatedParameter parameter) { + double xt = x(t); + double yt = y(t); + double r = Math.sqrt(xt * xt + yt * yt); + double u = yt / (r + xt); + double c = 2 * u / (1 + u * u); + if (parameter == x0) { + return -c; + } else if (parameter == vx0) { + return -c * t; + } else if (parameter == y0) { + return c * xt / yt; + } else { + return c * t * xt / yt; + } + } + + private final double t; + private static final long serialVersionUID = -5990040582592763282L; + + } + + + private class DistanceMeasurement extends WeightedMeasurement { + + public DistanceMeasurement(double weight, double t, double angle) { + super(weight, angle); + this.t = t; + } + + public double getTheoreticalValue() { + double xt = x(t); + double yt = y(t); + return Math.sqrt(xt * xt + yt * yt); + } + + public double getPartial(EstimatedParameter parameter) { + double xt = x(t); + double yt = y(t); + double r = Math.sqrt(xt * xt + yt * yt); + if (parameter == x0) { + return xt / r; + } else if (parameter == vx0) { + return xt * t / r; + } else if (parameter == y0) { + return yt / r; + } else { + return yt * t / r; + } + } + + private final double t; + private static final long serialVersionUID = 3257286197740459503L; + + } + +
    +
    + +

    + Solving the problem is simply a matter of choosing an implementation + of the + Estimator interface and to pass the problem instance to its estimate + method. Two implementations are already provided by the library: + GaussNewtonEstimator and + LevenbergMarquardtEstimator. The first one implements a simple Gauss-Newton + algorithm, which is sufficient when the starting point (initial guess) is close + enough to the solution. The second one implements a more complex Levenberg-Marquardt + algorithm which is more robust when the initial guess is far from the solution. +

    +

    + The following sequence diagram explains roughly what occurs under the hood + in the estimate method. +

    + +

    + Basically, the estimator first retrieves the parameters and the measurements. + The estimation loop is based on the gradient of the sum of the squares of the + residuals, hence, the estimators get the various partial derivatives of all + measurements with respect to all parameters. A new state hopefully globally + reducing the residuals is estimated, and the parameters value are updated. + This estimation loops stops when either the convergence conditions are met + or the maximal number of iterations is exceeded. +

    +
    + +

    + One important tuning parameter for weighted least-squares solving is the + weight attributed to each measurement. This weights has two purposes: +

    +
      +
    • fixing unit problems when combining different types of measurements
    • +
    • adjusting the influence of good or bad measurements on the solution
    • +
    +

    + The weight is a multiplicative factor for the square of the residuals. + A common choice is to use the inverse of the variance of the measurements error + as the weighting factor for all measurements for one type. On our sailing ship + example, we may have a range measurements accuracy of about 1 meter and an angular + measurements accuracy of about 0.01 degree, or 1.7 10-4 radians. So we + would use w=1.0 for distance measurements weight and w=3 107 for + angular measurements weight. If we knew that the measurements quality is bad + at tracking start because of measurement system warm-up delay for example, then + we would reduce the weight for the first measurements and use for example + w=0.1 and w=3 106 respectively, depending on the type. +

    +

    + After a problem has been set up, it is possible to fine tune the + way it will be solved. For example, it may appear the measurements are not + sufficient to get some parameters with sufficient confidence due to observability + problems. It is possible to fix some parameters in order to prevent the solver + from changing them. This is realized by passing true to the + setBound method of the parameter. +

    +

    + It is also possible to ignore some measurements by passing true to the + setIgnored method of the measurement. A typical use is to +

      +
    1. + perform a first determination with all parameters, to check each measurement + residual after convergence (i.e. to compute the difference between the + measurement and its theoretical value as computed from the estimated parameters), +
    2. +
    3. + compute standard deviation for the measurements samples (one sample for each + measurements type) +
    4. +
    5. + ignore measurements whose residual are above some threshold (for example three + time the standard deviation on the residuals) assuming they correspond to + bad measurements, +
    6. +
    7. + perform another determination on the reduced measurements set. +
    8. +
    +

    +