pax_global_header00006660000000000000000000000064122223242460014511gustar00rootroot0000000000000052 comment=c3cbfe32b733c20ca20006a45d2fce0cc2149930 pyGP_PR-1.0/000077500000000000000000000000001222232424600126315ustar00rootroot00000000000000pyGP_PR-1.0/COPYRIGHT.txt000066400000000000000000000037111222232424600147440ustar00rootroot00000000000000======================================================================================= Marion Neumann [marion dot neumann at uni-bonn dot de] Daniel Marthaler [marthaler at ge dot com] Shan Huang [shan dot huang at iais dot fraunhofer dot de] Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] This file is part of pyGP_FN. The code is released under the BSD 2-Clause (FreeBSD) License. Copyright (c) 2013 Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: - Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. - Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. =======================================================================================pyGP_PR-1.0/README.md000077500000000000000000000045041222232424600141160ustar00rootroot00000000000000================================================================================ Marion Neumann [marion dot neumann at uni-bonn dot de] Daniel Marthaler [marthaler at ge dot com] Shan Huang [shan dot huang at iais dot fraunhofer dot de] Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] This file is part of pyGP_FN. The software package is released under the BSD 2-Clause (FreeBSD) License. Copyright (c) by Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 30/09/2013 ================================================================================ pyGP_FN is a library containing python code for Gaussian Process (GP) Regression and Classification. pyGP_FN is a procedural implementation of Gaussian Processes. pyGP_FN follows the structure and (a subset of) functionalities of the gpml matlab implementaion by Carl Edward Rasmussen and Hannes Nickisch (Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2013-01-21). This is a stable release. If you observe problems or bugs, please let us know. NOTE: This implementation will be maintained as to bug fixes and corrections of problems in the execution of current functionality. Future extensions will be desigend for the pyGP_OO version only which is an object-oriented implementation of GP functionality. Further it includes implementations of - minimize.py implemented in python by Roland Memisevic 2008, following minimize.m (Copyright (c) Carl Edward Rasmussen (1999-2006)) - scg.py (Copyright (c) Ian T Nabney (1996-2001)) - brentmin.py (Copyright (c) Hannes Nickisch 2010-01-10) - Mauna Loa CO2 data (Copyright (c) Pieter Tans, Aug 2012) - FITC functionality (following matlab implementations under Copyright (c) by Ed Snelson, Carl Edward Rasmussen and Hannes Nickisch, 2011-11-02) installing pyGP_FN ------------------ Download the archive and extract it to any local directory. Add the local directory to your PYTHONPATH: export PYTHONPATH=$PYTHONPATH:/path/to/local/directory/../pyGP_FN/src/ requirements ------------------ Python 2.6 or 2.7 scipy, numpy, and matplotlib: open-source packages for scientific computing using the Python programming language. acknowledgements ------------------ The following persons helped to improve this software: Roman Garnett, Maciej Kurek, Hannes Nickisch, and Zhao Xu. pyGP_PR-1.0/data/000077500000000000000000000000001222232424600135425ustar00rootroot00000000000000pyGP_PR-1.0/data/build.py000066400000000000000000000206751222232424600152250ustar00rootroot00000000000000import numpy as np def save_data_regresssion(): # n = 20 # number of labeled/training data # D = 1 # dimension of input data x = np.array([[2.083970427750732, -0.821018066101379, -0.617870699182597, -1.183822608860694,\ 0.274087442277144, 0.599441729295593, 1.768897919204435, -0.465645549031928,\ 0.588852784375935, -0.832982214438054, -0.512106527960363, 0.277883144210116,\ -0.065870426922211, -0.821412363806325, 0.185399443778088, -0.858296174995998,\ 0.370786630037059, -1.409869162416639,-0.144668412325022,-0.553299615220374]]).T y = np.array([[4.549203746331698, 0.371985574437271, 0.711307965514790, -0.013212893618430, 2.255473255338191,\ 1.009915749295733, 3.744675937965029, 0.424592771793202, 1.322833652295811, 0.278298293510020,\ 0.267229130945574, 2.200112286723833, 1.200609983308969, 0.439971697236094, 2.628580433511255,\ 0.503774817336353, 1.942525313820564, 0.579133950013327, 0.670874423968554, 0.377353755100965]]).T # TEST points # test points evenly distributed in the interval [-2, 2] xstar = np.array([np.linspace(-2,2,101)]).T np.savez('regression_data', x=x, y=y, xstar=xstar) def save_data_classification(): n1 = 80; n2 = 40 x1 = np.array([[0.089450165731417, -0.000700765006939],\ [ 1.171605560541542, 1.177765337635947],\ [ 1.404722675089394, -0.017417915887421],\ [ 0.556096196907929, -1.489370243839215],\ [ 1.213163445267992, 0.044545401368647],\ [ 0.173404742510759, -0.675668036759603],\ [ 2.225008556585363, 0.469803193769368],\ [ 1.470329290331445, 0.887642323697526],\ [ 2.715199208821485, 0.621044646503113],\ [ 0.173640760494328, -0.936054178730056],\ [ 2.038152815025167, 0.262587298316711],\ [ 1.670218375320427, -2.633186886994263],\ [ 0.270098501389591, -0.948779657473203],\ [ 1.396339236138275, -1.114992287201776],\ [-1.482070589718501, -0.654590652482805],\ [-1.493788226272929, 0.382017940248275],\ [ 1.025083846875763, -0.860344923788873],\ [ 0.750316336734172, -0.101864205602753],\ [ 0.184311310148912, -0.258523866245887],\ [ 0.221868667121623, -1.393954437105630],\ [ 2.258881477897777, -0.786806071526136],\ [ 1.211362530151533, -0.423431246029886],\ [ 1.525307406741207, -0.097975367602030],\ [ 0.978930232706465, 0.476154349549524],\ [ 1.347884229346280, -0.248408186838667],\ [ 1.205779546204216, -0.090878327349907],\ [ 0.124388644862000, 0.599612645000285],\ [ 0.784044356662233, 0.356596736271853],\ [ 1.060216683845210, -0.318474838087900],\ [ 1.678114484474938, 0.678735373910422],\ [ 0.973851135005570, 0.024880700382574],\ [ 0.016237746864886, -0.480899874254564],\ [ 0.979406721923196, 0.697708815321128],\ [ 2.217307638531248, -0.956931847027775],\ [ 2.150475558834153, 1.059031573329512],\ [ 1.050502393215048, 0.532141747419667],\ [ 1.210593098269218, -0.318123542280113],\ [ 0.426309208807901, -0.571727978045793],\ [ 0.742552105732714, -0.122112766396886],\ [ 0.757210723588679, 0.862002000781123],\ [-0.431639130160791, -0.763118261936640],\ [-0.748398486307095, -0.603667649379360],\ [ 0.975086541108249, -1.525297946453790],\ [ 0.074503762788667, -0.092155036190678],\ [-0.668889572018935, 1.305400680048752],\ [ 0.725632503186580, 0.096286255882168],\ [-1.042270707136463, 1.297009698531055],\ [ 1.943144890398260, -1.051176922438962],\ [ 1.191448645802597, 0.261349747400059],\ [ 0.778004017505022, -1.046301123377022],\ [ 0.628873970760607, 1.103926629619643],\ [ 1.295113890591403, -0.479519217798997],\ [ 1.522065175744686, 0.993476032742058],\ [ 1.100255776045601, 0.961069161713818],\ [-0.593243832838153, -0.479418953496258],\ [ 2.023196521366462, -0.275055494808503],\ [-0.788103134597041, -1.090707985778480],\ [-0.085168420896236, 1.226858390046108],\ [ 1.691706923196703, -1.153144804780540],\ [ 1.989279380395157, 1.974704317386435],\ [ 0.398799861652602, 3.051291814188982],\ [-0.707217210772927, 0.185505264874794],\ [ 0.697550136765320, 0.222287208720035],\ [ 2.186126058382323, -0.327829143438683],\ [ 1.368068331060010, 1.708138258453435],\ [ 0.883049126818189, -1.334269372314072],\ [ 1.737643116893527, 0.618452933813739],\ [ 2.002228743955222, 0.103381966018445],\ [-0.202638622737115, 0.495024938090909],\ [ 0.543309203560769, -0.802120609128192],\ [-1.796161599703804, -0.054795478648902],\ [ 1.460693782000059, 0.750052171180825],\ [ 0.133277872804608, -1.154891068006907],\ [ 0.203670382700157, -0.480336687666025],\ [-0.278985011909341, 0.030578590108392],\ [ 2.070490237052893, 2.420782751903098],\ [ 0.599023881366768, -1.673208560658818],\ [ 0.140506592147238, 0.804938444757444],\ [-0.980799204108985, -1.847987723222053],\ [-0.102350006007740, -0.822093851434857]]) x2 = np.array([[1.160257057434194, 1.544111720606185],\ [-0.458434595629321, 0.205667827100987],\ [-1.053562345687376, -0.614938261650010],\ [-1.687901005751336, -0.780028275457715],\ [-0.467035854712698, 0.561692074343868],\ [-0.703391186121452, 0.281301267639200],\ [-1.568557779993616, -0.629129013661319],\ [-2.176478596101226, -1.176211396013793],\ [ 0.768109265900499, 1.376893437232103],\ [-0.514772970064353, 0.474264363701950],\ [-1.301924381487904, -0.525179228127957],\ [-1.312024947004566, -0.049469442305628],\ [-0.623417800418214, 0.226456899059445],\ [ 0.020290591370131, 0.374055846421580],\ [-1.002901826023476, 0.076597486786743],\ [-2.553713136283273, -1.731788289864902],\ [-1.788156378743716, -0.742460481943494],\ [-1.119582270077321, -0.256154464598782],\ [-0.423084091988017, 0.395108309297119],\ [-1.645945345460644, -1.216319293733455],\ [ 0.227805611684674, 0.925948003854262],\ [-1.298719171366801, -0.965511301629466],\ [-0.618292817021891, 0.140045887498202],\ [ 0.794935039731655, 1.917830760420081],\ [-0.213709179946402, 0.617751634356751],\ [-0.474251035850546, -0.054854432018974],\ [ 0.056077816960464, 1.046282980014428],\ [ 0.887136693467512, 1.536490289895764],\ [ 1.377161915854166, 1.764872700787871],\ [-0.901195709427863, -0.340855547886558],\ [-0.783104424735034, -0.330927422324566],\ [-1.507139570543989, 0.137504213149820],\ [-0.348999111724700, 0.235931187612453],\ [-0.367309385513174, 0.655996377722041],\ [-0.050622309620072, 0.410969334468070],\ [ 1.734919039047271, 2.611080177877894],\ [-0.567413078682755, -0.458249564234885],\ [-0.622230797920433, 0.258401595566888],\ [-1.642146761593230, -1.138579130251617],\ [-0.285298076847255, 0.085451489400687]]) x = np.concatenate((x1,x2),axis=0) y = np.concatenate((-np.ones((1,n1)),np.ones((1,n2))),axis=1).T t1,t2 = np.meshgrid(np.arange(-4,4.1,0.1),np.arange(-4,4.1,0.1)) t = np.array(zip(np.reshape(t1,(np.prod(t1.shape),)),np.reshape(t2,(np.prod(t2.shape),)))) # these are the test inputs n = t.shape[0] tmm = np.zeros_like(t) S1 = np.eye(2); S2 = np.array([[1, 0.95], [0.95, 1]]) m1 = np.array([0.75, 0]); m2 = np.array([-0.75, 0]) tmm[:,0] = t[:,0] - m2[0]; tmm[:,1] = t[:,1] - m2[1] p1 = n1*np.exp( (-np.dot(tmm,S1)*tmm/2).sum(axis=1) ) S2i = np.array([[10.256410256410254,-9.743589743589741],[-9.743589743589741,10.256410256410254]]) p2 = n2*np.exp( (-np.dot(tmm,S2i)*tmm/2).sum(axis=1) ) / np.sqrt(0.0975) np.savez('classification_data', x=x, y=y, xstar=t, x1=x1,x2=x2,t1=t1,t2=t2,p1=p1,p2=p2) if __name__=='__main__': # save_data_regression() # save_data_classification() pyGP_PR-1.0/data/classification_data.npz000066400000000000000000011633601222232424600202710ustar00rootroot00000000000000PKN3Ct^XXp2.npyNUMPYF{'descr': 'QFx>Urp>M/x>,"gc>;>K>w-=|6줴=0R)=lAP=rG=jwܑ<-i<+f<'=c3&@B>#ĕ>A\>j0k> \B><>V=a=lDy=;U=?>7A=0wj<8$!X-4hT4C~3w\Trte3Rӱa2&R Gm2jR1j1m80@w#_0:;)/+bI/ w_{.Ɲ*.nmܱ-=g-Bk,4g+ϣ_4+.?e**$)e ?}Q)6R'((#W(TJvX'nU<&,%{cJ%oq$#7(<'#^l"Sf! ku 7Ä- 'n7gp(T%4^?Bvj7?d 9? ?VKVs?-Ma?n7K?xQ3?g?&б>A/@`> #>KS@@+>sNYr>.giI>4H>Vr|#=QO<=~=#:[==kB$= 5}Ɠ40#]W4EL33h3K2GXXp2I1.♟Km1L,0Ϯbka0n-\)/٦~DV|K/NOu.d':,.W2 >-<'-^l, 4îl+^nr5+6 *"d40|^)21IR)#@('z(jX'5k= #&j'%UڭJ%2E$U=(#>@&#Q1l" L !A XU, lfbo[?N20?H ?z,?| ;?{?eLfh?HbR? P:?wD ?w&S? p&C>?_>vQ>fL8x>ƂP>p%>C==`yҕ=Sb=nXG*=C<> ۵Po_>M"MY>H\G6>?*8V>M}+>=J=!=H둚g=G>{1=1«8Q !8U/7 [7+F#?6`R\X6)5!6"Y5~1-C5\M.42x]4 ?&B3lv|m3XW20s2`O}1̚@Yq1  0(6d0p&s/{>O/#C|ϖ.8(0.E L3ԛ-M"?-E#n,ر+6+Q,*7JCm&):fIR)hC(0e(&X'y&+x(%\u*I%-ɯD$| #-kZ%#ӡi"yq!!p.& x(Dl?>ojʧ?lSV?JXL`p?A?ZT?n?"[_?ހv?ma?TxH?V>.?.ʙ?=9 >NN>uW/>v黉8>ྋk\>&1>e27>`"Y=u=E`LNm=Ʒ2(5=xkEy`4XOoy31gMUp3 2l!eu2 ygs1K靗r19v0ope0Wg/SJ6gP/Ce.ja0.fG Ŝ-n-o,&ł( +a_6+KplN:ȗ*4)R) ;w(q>5h؊(P~6X'&A"&P:rT%\H%!k@w$*Z^#;h-V$#Δh"2nЪ!hE?(a?Kr?8Z4?&?+|?@x ?eg?7naW??]i~?F:g?9^O^P?Vݬϋ4?.8F?XЎg>>} J>YcAf~>砦b>4s6>{= >(1,=ۨ=1ihr=:= )Z=EjwV<= <-<9.taF<涯<7D4';._Y`Hv;[.,;j:nj:ں@:p9̛(9[ >9xRVx8K"8&9&8}'O7Ң^a7T/65 6N\(&6s527ҖH5/hNu54 b41_a&3$gq3n22([v2Y\?b1\4s1ia]0ʽsqf0;y(/8l Q//3].ʎ 1.<)q--lx"o,Jg4+h 7+ *l )AװyR)Ȳ'W((|W'dj0&BJuL%?٥.G%~$#ULE,##+g"~?o ?2Re?,4׋?S.W?Yf?.uE0?C.{?gH·?:0?u]h?ZjA˳Io?3{U?[C;?҂Hw?>]:s>w>A\>*cg>$~<>.> ='R=A'bcv=} @=|T =l!6F 6P}(6T,Z"5vG:J5яVy+4c4Gq3zRr3l,ș2LSRx2_ MB1E#ϴt1a0Iddg0vg$%/tauQ/g.`2l^p1.Pҝ-lj -7xp,/'+7^"e6+?Gk*G/)YK R)$C/("rx(GV'(&_9T %vXF%ȿX`$(Ʀ#!#92 ?3<?B706|?@+? ɸv[?-.9?>]gԜ>>oD`>q?s n>A5B>:H >ک=S6={+-|= |rD=1\ =ޫ?_>g+*>2ur>RȤF>0W> B3=_&t=\C)=7*H=[d&=ԄDh<]*:%:6GI:9f7I9LmF9#@c8z,Ə8T1 08\iD76`h7kg7-n6$U-6eʫ5Lӷ5O50jY4 LRHf48s$3(n8u3'21z2|.ʂ1="v1(T0PXAh0?/s!Q/1@# .tTR1.~pի-R=,-(W @o,Vq+ rU5+f >r"*Xi?)8z׷P) q\@();u(ɹWT'AAo}H&CN%|g̩C%1$}$Ʀ ?y*??X?G?X??By*?}$Ʀ ?PhE?/%??S}P\?1hi?v.TF? :g?9mR^L?tip!0? c ?2x>S>bmϜz>N=;w>{K>֝~O,>?1#<|=6X=,@=S߻M=[jF==<)޾Mӗ<9V-o+PjYg5+qoc0*u;L9I)vvO)*;|($D8>n''6R's 쯤&F Y%rDbB%H h?l7z?zҞ?c*? .22?y ?0NX?ʉE8?A?7}_?q??|Y[ ?N20?6X?Fy ?'*lm? FR?l?.4?+u?=9 >o]>ʁ>4 s|> $P>QefI">M\=␴1=4=,TLQ==0n<ۛ<5vqkZ;xWg:(Kh֡:|ҺP: 195<ä9!>L9*U8׭48.Ywi438@T713;l7v7>6߈6ҳ`|06{c5]>aQ5y[4%?Oh4߱n36v3b\2d{2М@7}q1Ne;hv1xm´0*h00Q/VaXnӶQ/4|7.0.@9U-q?k-<Љ>FKm5>P1=>wOxT>zp)%>lJfV7=9VO_=.=gT=}%=;.T|2wzނ1\Zv1Mg.w05 3Oh0lr h/P_OQ/-<4.f),[0.s74C-Fa-vN=;k,nv+/2+Z;4Ւ*A_$)41K)'n>(7`'˿MP'_:I&9EjN?J Ń?KYY? ?Qd@$iu@?҂Hw?YՊ>@B>z y̚>:rJ>8yoX> *>J7 =LI!=8W=< #X=<=l<' Ww<9Ua5@R5T04hui4X3pã&w3԰W2F{2N1.\v1aͤ'0(ҎCWܲ>Iz>gr\>"/<.>ގSO=tqR=8e$=Z:[=:a.!=ϩ=+iH>%݌>7C`>&1>Z` > (= 'ܞ=~_=`$=@$?z+B?D3S?7r?[е?b?g`r?~")i?U+5K?)stM*?3I_?lC>>8h >qu8\ʐ>/Ĕ)c>ڻČ4> \!n>m$=p[=a=:wv&=L9:}<6.? z: ?9>Vmka'>dŶe]U>ق`f>Z^7> > =1z.ɜ=)^vc=ͨqG(=., @{#@b]&@>&@]Ȅ$@ԝ @/@ @(!gxm@7N?*?a/?FPye?SE? ^z ?'(q?R??1?`4v4?o>E>̽>!Sp^i>ěH3]:>E}ٹ>=N%p=  f=+ p+=|H|\gٔl ?p0X H@@V6@V+@}@<$@e(@.+@;v+@ݬ\)@n-?$@19@q6v@2齀@ư?JHp=?:c۩8? pQ?T&?ߐd?Jڑlt?{U?Pem4?^?zN>T%t>0VCۘ>1yl>Go=>e͇w >%:X!=-,Q^=xh=1=N<-=4<ݴTmӰʫw?T4|!?h?)'^F? J?6^N??'^f @T^@{"@:aO*)@.@"}w1@"}w1@.@DaO*)@"@r^@ '^f @X?^N? J?d'^F?eh?4|!?>ʫw?]k X?)7?6D?x\k>}R>>bb9o>LB=D@>*G#K3>v+RM=Ibvj!=Sj=W1B)0="m1<9NÖ<>pUr3ߋ1s2}xt2͍f´1:߭n1;}0.U^0xRI//oC/HTo.dPQ!.hM- aTj,ZZ,O+?N>?v%`?ؤbx~? O?Cj?`[Uʪ?MwR?zҞ?̙ ,@dJ_@gI@ag 7&@}4k.@ժNɁ2@k4@hzsɾ4@롷2@HR0.@V&@ճ0qm&@ (-]@y@](?#z15̿?fC:?_ Ǟ,?ǽ?z?]o\? P:?[+¾?Z]/G>=Zj>&,;>ycq>4LB_*n>}_54e=Y,Aћ=ak=?Nt%1=֍??V,Y&\!?3P[D?Ld?DT?`.s?+GP?0x?}Ƽ?30v?W4@2 ;3@=Qo"@x/ڵ,@@I/2@WPT6@^$9@]:98@_{;6@1@#2n*@^\!!@}8]ZN@ad@ο?2ޭdG? %?nzh?qUw&?U<~? T_?J=?A-:D? Tb>rw`>lƠ>}Pir>|DC>2i>tT[=^ 1=;Ѐ};hP;eo ;tX:d2]:6%:mMCŰ9:U9@d8F7~18 8817X@p77H]@R6!06r:M_5 S 0N5:4 1Jb4 3T{Fo39Uv&D2)ц4q2s9@15a;4i1gf0hmTX0vT/)vm@h=@C:@ZK5@ދ~\.@U;q#@)Lq@Ł+ @6cs?}Wp ?5(RL?_?6X?r?ma?@:!l@?8?qȏy4>|_27>*>.bt>.q˥D>&U>b$= /-G0=x{*o=:!2={O;?[J>,?HP?O p?v'? ` b ?&]fx ?g?I)4?zI@m4[@݌ i @v"(L+@uQ`4@Sm';@8$T@@A@ ]A@ >@~;r8@!?1@62B&@"8@RP @>?+?RK?á?{L?ަ蓼?rEc?`:bA?Wi;?N>yQFx>Jew>u>AҮE>H >scQfDY="8=DS1p=sD`3=dI<.h>d۩R ?Z2?pyluT?8t?/l _?~(rO)?SV?L?‹I?$iu@29q@}@<$@ԃZx0@:LQp28@IK+ @@I=N>N!>iyv>pF>n\8>;0X_=_ =aDAp=a,3=\&ij٧>dwW>O?PVr6?z0qY?!y?:ח?xVWKʳ?mo?-$Ʀ ?!0v?f[ @w@YhB(@iu3@]o<@f^B@ڜ&SF@aZiG@KۖRCG@X[D@V@@`PT6@]uan_-@^\!!@G\@|=:@i(?9Q9X?{?-贷#?Ӟ^?' N6>LOuw>xG>>&y=%"=i{p= A4=JW5E~Y4E5Y M/4hkIY4;33پqid3Rn'2۹?e2 1xȸ1?^J]1whI^03=VK0${./ 1/F&j>^k5>Q)B>G=>Ie?O}@;?imP_?A,@cqfZ6@)Wz@@zE@%͈/6o~>ah+>=Zx>A yH>8?O>r4=C@=O8q=()O4=9(GK#FIp>|!c>|\uo>>u>Z5o?avA?…c?j,=?ë[?s؏Ƽ?+6V?No瓎?P?@Ey@h"@vyƎm0@"2.T:@OF C@"H@O#MM@d>)O@ pM@+ yI@WGAD@?b<@Tد32@J$@i а@ОR@;rk,?X1?S [?f?;T?ki?R7G?mE##?fe1z>`>z,>IdYy>9/oI>,tL>!=KP= q=sR'4=$><ʻQYbt>j>̯9D>%B>ܔ ? D?35Cg?U?JO?i?4?^% ?|*@\@n(Vi%@ɪNɁ2@ 3*>@wE@2L@b}^P@hs|Q@Ԟ9kP@dzL@{QbF@mam?@kYQ3@ή>&@)Lq@?FV@ sg? ?|Y[ ?4?棧?[_?SFyj?n H?G>#?5J>?>DG#[>1i> y>,7I>*fMZ>;<=/;Cө=%zH7p=<'3= EzH_'I>m5ry>\쏼Ƨ>u>r:>tZv>]쏼Ƨ>im5ry>r{H_'I>bC.>=%RS]=Sqp=a>3=O.JS<_585|8H4]L43}!cU3+ -LI2buEV27d91#$UM1yͯ=B=ZM>{BN>|7%o\>kcѬ>}\O>'}? i$+(? GPxM?2:7p?/{?$PX=߬?wn\?Qe?}jHz?W' @V@i.+@auQ8@<}>~H&>7#y>~WH>>295@=u=b'o= 2=KZu<>srR>$>2]F|q>uN=?MJ>,?-{QQ?N r?@wv?)?7r?1 ?F?޷EQ@`~0 @ɝE6/@@m';@: RE@`[7N@&F>SS@I+&OV@MD"FS>W@|Z0U@*bR@} iL@z)C@m;r8@v"(L+@UK @S @ F M?z-.9?e`v?d` b ?V~?Ӊtl?IyvI?mԇ:$?l>ܯs>W+i>p$,y>hUjPDH> .>b=~ ܺ=gQxtn=~ L1=fM 8p<¨,ne<-fFm<%8?B7u67>*-Y7\6h˃6/DI685)L&x25]?뙼4tgD4mSF3{DN3ߝ\2Ǡ{fRN2ĚB=1#=8=(t5= $>`(dV>7KY>䤔>h>?.1?J^0?o%T?/Mu?z 5蒕?I=4$?njʧ?B'?NX?Y\}i@04+"@4_7V1@ڑ$>@7OqG@b}^P@O<]U@-Q]'X@! Y@C\FlW@}3S@J.#N@qTD@q5I-9@a=6,@ó0qm&@m# @t4 ?/ԅrI?4l?\.v?| ;?sql?ͻeI?[g$?e >\"^>hY3>W4dx>QOxG>b>#z=.*=cN2l=>p{0=;"K<5ڋ;< mƀk<0%wlWZ>!>yy>5>8 ?}xQ3?{2W? V#y?U{?+7?frJ?.9K?jh"@wF8@{#@@3@TcWj@@/+7I@wPqgQ@9V@Gb8Y@g:Z@EއX@T@^O@HE@8w:@yJ;-@ @(pͦ9 @7N?lU??v/?f?%E Eٌ?L~?>>E!>p[xw>kҦ\xF>i>IF=ON=%$ k=^N/=x <&nu <ېEi< 9#X؁88j8* p7&YIiR7}R&6QE^|6Ƀt16iQ5ii$)5G:4X;4,3m8D]D3ک4<&.t =)w.K=Ud.==,i=(`pG,>:@_>$e,>aǞ>Zp>ZV!?U!6?*;bpZ?P|?(n?5bo[?V9EjN?ʹ?I @CcZ@s%@Czsɾ4@hA@sa{K@}bU-S@%Q]'X@3zu[@]n ,\@4VsZ@&U@n$HqP@pQbF@D}J;@~\.@zI@\ @jHz?i?Nx0?ZT?R[0{?Nk?xH?iyV?D|#?7>dc^>ʰl/>M3[v>JE>4TM> Pe=S==Q55. i=dDD,=+|ڀ<,iƝ&-b>;N >t_](>@<*>w牐?Y9?>]?0@?jLU?G˨?6V?2 ^?9e@w8o@EKRЩ'@fX6@;F C@&#MM@nWT@|^$}Y@uok\@V`m]@'*[@V@ LRP@jGgF@Jk<@.@}D@>+ @k?y(i?D'^F?a:#ܩ?gԋ?_ǔk?)G?Y䈉 "?>"<>G-@A>R>6ou>u00]C>Vg>|{4\=B=k\=둚g=;j53*= -=d>gґ>] h>uDc?n* 9D>]Ԣ>s>THB>vpa>*|\ 9=`'aɠ=ij=ۡd=>'=NkU'lr5>5g>+XOn>o]\>"^>1?m7VI??fb?{DT?p?^Hs0?hE?BnV?IK@v2@Yp*@l^$9@ږ56E@g&P@"1V@/D:[@%8^@h"_@!ldc\@ \W@u&ZQ@rKۖRCG@v%<@K`@.@;/KK@f[ @X?0+BX?T[f?$Z/1?Wʼn?#:h?E?.GJb ??XЎg>SR4v\>?mϜz>]o)r>isA>'m >=Q `=lon1b={g"];%=olޅ<6:rlh9cPmSg94@i19riF8Үq8"t,B8-74pjA7jK6i6uGB5\E5v5~ٰc;2uF8>]dj>F6L>R]1>Jd+v>Є@[&?쬵Z43A?0}@!d?S5(C?CWcD?JXL`p?@(?9A?ƥI@Np i@=>,@GC:@! F@9kP@2zV@Kn ,\@ŇX`_@!e_@qݛ\@C\FlW@NQ@r,G@=sd;@HR0.@S+@m+2 @P ?:~)?e`A?  ? VK~k?A[g?Q(O)D?dk?'>Χ >wvQ>Wp> u>> p >"Y=6 z4=yO`=MPC"=*S<_"G<#-b\6>Hំd6o`5~~U5b|;Ү:;a1N^;>p?s n>R6>L٣f>pL>H?IkB?3 ϖe?I>U?s p⸥??U?v8?CKI"#5 @*n@-?Y-@#m';@'F@3XTQ@1D"FS>W@JT$\@һ2d_@ջ2d_@JT$\@AD"FS>W@FXTQ@'F@Nm';@b?Y-@*n@KI"#5 @8?U?O?Ot p⸥?TI>U?4 ϖe?kB?jH?L>٣f>jR6>|q?s n>";>^>h= К9=gQ{]=v}] =jwܑ<1|<8ÆX u>>Wp>WvQ>^ >>)dk?G(O)D?A[g?UK~k? ?e`A? ~)?P ?Wm+2 @-+@ߋHR0.@ sd;@r,G@xNQ@C\FlW@qݛ\@!e_@X`_@Kn ,\@2zV@9kP@! F@WC:@=>,@ Op i@ޥI@9A?@(?JXL`p?kWcD?S5(C?0}@!d?:Z43A?:@[&?Jd+v>]1>sG6L>]dj>8>n*m>;w5+=Sc)=4BY=*b"=U%<7˙pVB`.isA>[]o)r>lϜz>Q4v\>WЎg>GJb ?E?9h?VWʼn?K$Z/1? [f?+BX?UX?[ @.KK@`@.@v%<@LKۖRCG@lu&ZQ@ \W@ ldc\@h"_@x%8^@/D:[@"1V@g&P@ږ56E@^$9@=Yp*@2@]K@onV?2hE?Hs0?p?DT?jfb?xm7VI??1?"^>]\>H,XOn>+5g>ڂr5>٣"jU->=N5=~`=ArRhV=—2=<۔|< ֤<OQ<ib HB>s>%]Ԣ>9D>Wa{_>hڣ!?WNF?3puj?h`G?`B$?tE0?]?n=R"??bn @9@8}4k.@Z+8<@E=v4G@'.8Q@W@u[@e}l^@/:]@-e{Z@[U@j:N@#]4D@M,x7@ݬ\)@ H8W\@AK@񂉆?C.W?Iuc?ni?eā?b `?}n*]gґ>Wd>j"3>ɷ2=Ls2=rc=)S="=٨Z^m<{< sqM< էI<ծ[]U;]-lp;Cq!;IA:J~:m{':V9 [u9 t:vr9,o~)8|j+&X8I ~7 yVwc7[ "%7 E: : ȷ:<;U ;*+;U5'<%c̰d<;pbfͨ<4-=/0]C>C5ou>0R>-@A>"<>㈉ "?)G?^ǔk? gԋ?aa:#ܩ?&^F?x(i?Bk?1>+ @}D@8.@k<@;GgF@LRP@V@*[@V`m]@uok\@t^$}Y@nWT@,#MM@OF C@vX6@KKRЩ'@w8o@9e@w ^?06V?˨?LU?@?E>]?݇Y9?牐?@<*>t_](>N >E&-b>Y0>s#=#-= S =8P=~"=F]XH<n<}|A;1W ak;!;:>$x:7": 9_p9 -tgl9;7!8UR8\7h:7<֢J9.CH:wo橝:,O9^n]:GK@;lEJ;'qL;;ld!핷JE>gM3[v>ʰl/>c^>7>&yV?D|#?xH?Nk?[0{?ZT?~Nx0???}jHz?\ @:I@w~\.@C}J;@LQbF@V$HqP@&U@VsZ@@n ,\@o3zu[@Q]'X@}bU-S@sa{K@hA@3zsɾ4@s%@[cZ@S @ʹ?_9EjN?wbo[?Rn?kP|?Q;bpZ?U!6?ZV!?Zp>~Ǟ>%e,>::@_> `pG,>Bi=3=d.=)w.K=2'.t =M4<+$l<:Sة@D<3R;Lu4;?k!f;9d;v-Q:ce7s:E::Է9{ʍj9=J9O#㱍<87>'M8ΜDrj7KI{9?d9]n\K:$xDѠ:r+hMu:- JB; }Z;Yd; 9#<ڐEiҦ\xF>[xw>!>?>>W7>Q $?-jI?K y><!>wlWZ>y['>';Nr=\q=)M4߂=81F=ǀ M=(n֤<$Ä<I Pt@p{0=AN2l= *=#z=b>HQOxG>V4dx>VX3>\"^> >wZg$?̻eI?sql? | ;?-v?l?ӅrI?Ct4 ?&m# @C0qm&@=6,@`q5I-9@pTD@zJ.#N@M3S@C\FlW@! Y@Q]'X@C<]U@ b}^P@7OqG@ڑ$>@4_7V1@,04+"@k\}i@@NX?B'?Eojʧ?I=4$?>{ 5蒕?Mu?jo%T?J^0?.1? h>夔>7KY>!a(dV> $>t5=m8="3#=}B=V:tR=ҹw5<"3Ȁ<5v:<(Jo;dr0;a5 \;V ; ::Y,g:}E(:p@690@_9ЫR91`&;8(Fo8@ܓE9a9$9MexjQ:qbr:Hw:-7HVF;76B6œ;F,r C;;c3&<,fFm<#¨,neTjPDH>t$,y>ˣ+i>zs>l>Tmԇ:$?oIyvI?Ӊtl?^V~?_ b ?`v?,.9?F M?0S @fUK @Dv"(L+@$;r8@F)C@6 iL@m*bR@OZ0U@%D"FS>W@kI+&OV@&F>SS@l`[7N@: RE@Im';@͝E6/@T~0 @EQ@?1 ?7r?A?pwv?N r?-{QQ?J>,?uN=?3F|q>2]0$>irR>" , >ҏv={+i=b<ny=w=$O>=2n =M(aA<3waz<P{5~WH>:7#y>7H&>|>u >v$?'DI?+RAl?uf0?'gH·?Q?l;??@L?轀@'&h@[?j&*@07@*5kB@\v6OJ@,.8Q@eR/T@p[U@fLjdT@5GwQ@BuVeWK@0kcѬ>7%o\>{BN>ZM>NB=ͯ=Rut=V-8=~jD<.fu3=Sqp=RS]=K=bC.>zH_'I>l5ry>l\쏼Ƨ>u>r:>Au>*]쏼Ƨ>im5ry>zH_'I>bC.>=>RS]=3Sqp= a>3='.JS<_<ܪp<)<:v";w,^;B&%0J;lF4:6f:gJU:U$9KorQͥ9j$7b.~.8p8ݭx̍8P M9Z~^٧9\Z:0V:qq~:3V:RVK;-U;#;2,+<%mٕhq,7I>0h> y>F#[>؎?>4J>>#?9n H?SFyj?5[_?I4?棧?V|Y[ ? ?sg?FV@ )Lq@ͮ>&@kYQ3@`m?@$QbF@WdzL@9kP@Shs|Q@a}^P@ѵ2L@\E@ 3*>@NɁ2@(Vi%@)\@{*@n% ?I4?i?_O??5Cg?D D?;ܔ ?&B>,9D>aj>FYbt>?qTD>d׸>BDXY@=&ë0= j=eKR.=<R'4=0 q=JP=9!=sL>9/oI>HdYy>z,>F`> fe1z>{mE##?;R7G?偎ki?L;T?Xf? [?X1?;rk,?2ОR@ а@ȜJ$@د32@r?b<@ GAD@t+ yI@U pM@d>)O@#MM@"H@*F C@1.T:@vyƎm0@h"@+y@P?@Uo瓎?@6V?e؏Ƽ?ë[?j,=?…c?avA?/Z5o?u>\uo>!c>#FIp>LB=D@>,qe >Am[8_=dA4П=1"P!jd=1Ud'=D]rR@ yH>=Zx>h+>?/6o~> |+g>ˆR"?cF?Lfh?A??a,@:t@iZ1!@-22?ٺ ?w{?S?mM ?1lQ)B>k5>1F&j>^T<<9>ag{_>Q v =wUΘ=_=]D@"=*S<ڬhm9ѡ*xG>&LOuw>2N6>ϹD>Y&\!?D?' ٧>2>e>ig14>J1iJr>=N5=DOm%=ؚ'RX==C=%oF>myv>;N!>M>6%B>t"H ?sC?e#dd?Du]h?'6?h?E?{R?JHp=?Iy ?"E@נI@[?j&*@ WX84@Z+8<@hA@#]4D@pTD@Rq>\%e,>ֶ r`>tUIT/>Ox=g= A=+|R=)Z>=];|ѮE>Vu>ew>~PFx>YN>Wh;?&`:bA?\rEc?"ަ蓼?+{L??RK?+?>?RP @*"8@1B&@!?1@;r8@ >@ǻ]A@A@$T@@m';@uuQ`4@Dv"(L+@v݌ i @m4[@PI@1)4?g?]fx ?_ b ?Uv'?h p?HP?>J>,?;?dm!>zʴ>·E}/>O5LdY>1yA(>o;5=2¿='Qd=sfK=24􏰐=)<`+p˥D>.bt>8*>A|_27>>qȏy4>?:!l@?ma?r?}6X?_?4(RL?Wp ?cs?[+ @E)Lq@;q#@Y~\.@K5@@C}J;@\cqfZ6@0@9g 7&@ l@r}{@0~,D?N|W?mN?tͿ?5O?)숖j̈?bSFyj?IG*ϨI?9 a&? ?O>`G~ >6i>:8i3fS>efI">E¼=Ά=l ⎀=,WuD=]DyRb=L˭J<' %<ET/B<%>;fj;/:46S 0N5:M_5!06]@R6j7,X@p7u17 88^D7~18%@d8X:U9ۆMCŰ95%:u2]:tX:o ;<hP;};q=;cHu,U{DC>t|Pir>3lƠ>qw`>^Tb>v,:D?I=? T_?T<~?Tw&?zh?U$?ޭdG?Ϳ?,ad@ 8]ZN@$^\!!@'#2n*@41@^{;6@\:98@P^$9@PT6@H/2@/ڵ,@Qo"@ ;3@W4@0v?}Ƽ?x?ٮGP?`.s?DT?Ld?3P[D?Z&\!?h^1>| $8>窍>BXD;|>8?PVM>I}骃>uVJ=(ʱ={x=?!>=p2+%=k~a<4TE<=B:<<O;gl WQd4&4:4ӬBP5p ]5)w 06$& ~6+A57Bh p7~f7C;iv98"Κ8DOw8`,qU9hTj9DO : Lq]:YK:Ec;p&!P;0cZ;TZ;IǼ+<5yq<VCpLBTycq>&,;>=Zj>]/G>+¾? P:?nn\?z?ǽ? Ǟ,?fC:?Uy15̿?(?Fy@'-]@K0qm&@GV&@HR0.@r롷2@4zsɾ4@mk4@NɁ2@Z}4k.@Og 7&@*I@dJ_@̙ ,@zҞ?@wR?Q[Uʪ?$Dj?!O?ؤbx~?*v%`?3@N>??z~A>|_27>,go>H=u>Q"E>s|>95@=1~9=Cs8q=Q3!`E6=p <<( v<߶2<~3yTe4Tn`4k׎þP5a;;5+ 16l'608yc7t?_q7:b\V7v98M>]&8r=8 3*nU9iV۰9EfOm:2;ab]:2¡: `m ;I<05VO;`2);cBv;;W*<>p<8NÖ,LB=D@>bb9o>#>,}R>]\k>86D?)7?]k X?>ʫw? 4|!?}h?&^F?| J?]N??'^f @^@A"@aO*)@8.@!}w1@!}w1@8.@aO*)@J"@&^@'^f @?^N? J?D'^F?h?J4|!?l>ʫw?x]k X?)7?`6D?ۨ\k>}R>>Lbb9o>}LB=D@>*G#K3>]+RM=bvj!={Sj=71B)0=ڞm1<79NÖ<\>p<"cs3:~w3hQ:zf4Pdqܐ4FtQ5Xa5n\D26tJɟ6^d7nKq7;r7_T:8,%k888n\عU9@hO]9Kt:pƶ\:=Y+:%oo;Q'N;ͮQpO;ͧ;t))<xXn<5TmӰ<4Go=>.1yl>/VCۘ>S%t>yN>?꯷em4?{U?ّlt?ߐd?CT&?oQ?b۩8?JHp=?Ű?轀@q6v@ɳ9@5-?$@bݬ\)@ߥv+@<.+@d(@}@<$@+@V6@p0X H@ ?~L1L?9H?ngI]>i>/VOg>Yb7> >X]=M=db=)i.c &=*$m<>SrؠaH3]:> Rp^i>˽> E>o>H4v4?v?1?R?&(q? ^z ?%SE?cOye?a/?C*?R7N?(!gxm@f @=/@ @]Ȅ$@>&@b]&@p{#@ @-  q@z,/@IL?@fʅ?[yӜ?,m?dֹ?_"W?!vo?J),k?,[L?"8"ٵ*?K?$ظk>ay>2\<'>x^fa>jDP1>\h=[%Ab=l8e$=53Z=ecC = qG(=^vc=z.ɜ=ؤ = >eY^7>sق`f>|Ķe]U>kka'>9>x: ?>.?q`pO?7)lm? .? sQ(?>`N?4(RL?i;?sg?j"B@?\ @mq6v@Np i@a阗 @/4+"@keϔP!@:t@ZaPa@dJ_@I@ءxP?5#?M ?fC:? _y?n?H?ɴ`Vd?~QxE?>#?r%? %.>l|a>T`J>TM1X>Pc(> *=(=b=$R==P o,1z2n^'fj]2TuGv3c5ᒨ3Ci4(˸P4d-S5D5268C6 A1' 7 4q7&!,7m`#981E˙8&8hޖS9t9):yc1Y:B:‚Ũ:"yCCI;nvf+;:[EP;Ft!$ٻČ4>Ô)c>t8\ʐ>Nh >lC>>dH_?qtM*?{*5K?盭")i?ȏ`r?a?[е?V7r?t2S?+B?>? I@[t@aEȻ@EdH`@*n@TUK @&=I@`}1@RP @=I@?y3<?|? hϒ? Š_?tK”?f˹Sz?!^?6XO^ ??H?Wlףj>1>Yj> q> oQ>!>=#/Ng= =>4q'J=k餪"u1}C1z2$3+2%&j(v3dX q3xqi4 P[p4HS5UDp(5m{26q 6~y7'q7+i74`eP88.ј8ō8US9v:ٴ9HnNT:AV?W:^6:2%:!9UG;|U;S H",;w O"<e<5<@$<`$=?_=L 'ܞ=0 (=Y` > &1>V7C`>%݌>+iH>ƀ< >%S? &&?NF?#dd?`|d?8sY? Ǟ,?TNx0?COQ?~L1L?ct4 ?' @ m+2 @'-]@cZ@9q@ߑ\@ ҷ=@=}{@FV@󡎀xP?g񂉆?G h?Yf?h?k٬ ?c0Ž?GE]_gs?HTV?}Y6?%ѥ?8eXA]>a>&*>%M3[v>~WH>vD_>{=d=m{={偭0u11f٥{2`֙"2pvw3k&p3.Ci4Wlqe4wR5"15ݱ߲26a6}(7Pd p7*kq7(F78vף8{8pϲWQ9gĪ9 :U:}:Q:VTE;͙t6;`B;ўd{ Ufr\>"z>CWܲ>$>4aiZ?udv"?CnPB?,Ma?) |?zÄ?*?M޻:?a/?Cfp? =up?e\=?(a@oͦ9 @Z{O4@Ғ{i@] @?M @fZW@<8@ii?-9K?[M*?l?2k ?vj7?uэ?L=yjK>3ER>Ɗdwo>-@wQcA>1rW>wnV8=RS]=V(5hM16 KBZ6 ?7zo7777i!@68 LBI8nͫR8P9)٨9T:i0(T:g9i:{:!:]q&MC;9[ ;s$:t;"t(_<9Ua<Wwb8yoX>9rJ>x y̚> @B>=Պ>тHw?>N>?/`rB\?.Vw?ߐd?qd%f?h?ԙ{u!?:?-l7z?n=R"? @ɖK@#+ @W' @T @U$iu@7Qd@ ?NKYY?* Ń?n9EjN?iV+?~dݮ?>@?'݉Ā^?`,M%qd?3( G?J6o'?:|}? K>&Rq>Q{>Ue>h8>LT>R-j=Ahr h/ 3Oh05f.w0uZZv1&zނ1_<.T|22/9v3ϙi\-3Axi4Q#4 R575R6O16e6s]A7:n7̒7]`]48ˍ3˔8#8S:eN9[F5_Ԧ9؂Ct9:UR:`z!:L'@:@QIA;(l;I?;~awOxT>P1=>{EKm5>?>,+|C>hz?GY9?js V?y΅r?gԋ?DF쟢?;K?^aߟ?iv!?Ӱm?;rk,?q?k7@s4@ОR@v9e@]Ѫ@{ocMk?"^3]?mnӖ0? J?Au#P?tUa?cf?6]P?%>ʫw?88]?Ka@?#/ ?=eH>3N>M >8>dp,~:^>Y0>Ɛ?R>;`XnӶQ/0Q/)*h0_m´0.c;hv1)@7}q1d{2a\2j6v3ɰn3c?Oh4[4>aQ5{c5N`|06=6߈6u7=13;l7g@T7HXwi438V48@*U8 >L9:ä9P19{ҺP:Jh֡:ʚWg:ҷȥ>;`۾Z;aw;u`<5vqkZ<߁ۛ</n<=TLQ=o젘= 1=M\=efI">Q$P>l s|>ʁ>o]>8 >c+u?1l?.4?@ FR?r)lm?Fy ?p6X?]N20?V|Y[ ??*? 7}_?A?ɉE8?MX?Iy ?-22?c*?zҞ?l7z?G h?6rgl? ?I=4$?(TŚ?mUK~k?(p?p|T?7?l n7?n0>{'M>>Ռg>V6T>G=Es&>Kѓla.oJVQ/qT/Չ'h0*0mi܈v1#u8A1)K{232u3Tqh3zlseg4F@}K49 P5 o[E5G/6V2ә6u.]7(I5j7D7nS 18#ӎ8l8oPmI9}q97v09Ke,K>"N=;w>lϜz>6>2x>G c ?hp!0?%lR^L?B:g?-TF?0hi?& }P\?̇?t$?hE?#Ʀ ?x*?,?X?G?JX?M?x*?$Ʀ ?hE?U$?Ƈ? }P\?1hi?-TF?:g?lR^L?$ip!0?w c ?2x> >lϜz>9N=;w>K>QR1./@# .Ľ!Q/?/MXAh0<(T0="v1gz.ʂ1z2W'2&n8u3k7s$3 LRHf4.jY4Yӷ5O5bʫ5.$U-6n6kg7S`h7PYiD7mS1 08,Ə8y@c8WJmF96I9O94GI:%:F>:iؑZ7;wzv^;:)_;?m<RS<_]*< ӄDh<9d&=(H=ZC)=2_&t=)B3=W>RȤF>}2ur>e+*>5_>HcJ'>{ ?ydO)?NF?ސZb?z?mƩ?}?@?yl??g~)?xp{?ʹ?KYY?% ?o64?N|W?Wp ?M ?Tڔo'? ?M?AqZ?.{?EC0x?BN&oe`?D? &&? L%?#>\~>]*>6n:o>=2\}-Ԩ )ؒ1.4ˤ. Q/m V]/ƴkh030]2ku1bۯ1y+y22.1\t3ԗ3"ќe4[` 40(M5'?f5O*6(y6jY7;e7lU{762dE,8^u8c؅K8]z).JD9hja9ÓlN293aF:Yܷ:i<':~T54;t+k;{0;_w5B>p?s n>=>oD`>/fԜ>kC>>?OpK?%d#?0A?U\?xs!t?4?a?v2?p?uܚ?@&>[`>Ն>lj -Pҝ-`2l^p1.깘g.auQ/e$%/Hddg0a0xD#ϴt1 ^ MB1SRx2,ș2zRr3bp3 c4Vy+4JuG:J5+Z"5P}(6 6[>6kc7\B^785$)8$98*ڨN8wzA9@CJT9=9C:X:g#:(1;Ųez;MGRԯ;VdE<_ K$~<>lbg>@\>w>z:s>>e҂Hw?3[C;?e{U?iA˳Io?u]h?0?gH·?DC.{?tE0?nYf?.W?3׋?Qe?\o ?.?I ?k9H?^6E?%V+?f9[?.|_?h _Ȓ?TQ~?5Cg?PN?D2?9nFl?F%>Q  >kcѬ>mx"o,9-)q-Ȏ 1.3].7l Q/hy(/Fɽsqf0#ia]02s1\?b1([v2Ȑ22{eq3]a&3~ b4x-hNu544ҖH5똃5:\(&65 6 R/6Ң^a7|'O7`9&8:K"8iwRVx8[ >9̛(9Bp9ں@:nj:j:[.,;]Y`Hv;4D4';_嶯<.taF< <-4s6>砦b>$WcAf~>ѡ} J>>vWЎg>..8F?ݬϋ4?B8^O^P?*:g?]i~?2maW?of?cw ?iz?&?mY4?Kr?(a?ghE?[z8?W@ ?D?[f?zGP?.7?jjb?"Wʼn?Du?nB;_?呀D?@wM)?TIO ? v>}>#ł( +o,n-C Ŝ-ha0.Ce.FH6gP/2Wg/ope0Y|9v0K靗r1ygs17i!eu2 2gMUp3 Moy3a<`4Z]^24®?F5eHk5NK#6 6 x]}6^-^7~Tr7K)%L#8Od$88t+-8wDm:93C9XMb9"R<:BoG:gxQ:6';u.zr;ÍB; ;-"%1>޾k\>u黉8>W/>tM>I8 >j.ʙ?>.?-xH?ma?ހv?[_?um? ZT?@?'JXL`p?xRV?Anjʧ?sw(Dl?w{?8Yf?[Uʪ??㢖h? " t?qZ?k _?y`|d?sql?n[T?@}@;?SX ?:F˳?%'> 6+Fر+E#n, ?- L3ԛ-8(0.B|ϖ.{>O/m&s/'6d0N 0̚@Yq1,_O}1v/s27W2iv|m3 ?&B3F/x]49L.4ѫ~1-C5Y5z(5!6 Q\X6D#?6_ [7U/7Q !8#>8K{8"6951! 99Ha7:^9$G: %J:h#;z%VVn;d 6l;Z7;i=*8V>VH\G6>w!MY>Po_>B^g/77>r:9X ?g'?c&B?ㄋY?FU#p?Pj,=?3|!?p޼?33?9X?J=G*N?&^F?paߟ? =?b?+>S,?؏Ƽ?9|?fҧ?'t(_kH?f+?Ƹ>ʫw?Z…c?_IK?:r2?ء?)=>t;.*2-6+ 6%+:5m,w-C:&JT- 62..&X*.tM/rN/ b0By0^Ip1ϭ'1Ckr2mj2"q2k3F 3s0Z4934:A5 ۵-ƂP>L8x>vQ>5_> p&C>%S?׿wD ? P:?bR?Lfh?l{?| ;?z,? ?1N20?ao[?t L<^?Si??`Y?(tͿ?>`N?$ _y?-v?SŚ?mƩ?O`|d?BiA˳Io?S?Y?lbB?a i$+(?۩R ?`40|^)͠ *nr5+ 4îl+V^l, <'-/ >-$:,. NOu.j~DV|K/+\)/-bka0I,0+♟Km1I1)FXXp22H23h3nCL3!]W4zƓ4 Ƞk>5h8o|X50B^6a4/a6f6bS7ko78{%Sw8<.`h8~Ev/9J944}?9N'0:MsZ:kT:*Ь; ԑxc;V%;%a;$02<,~r0giI>rNYr>Q@@+>" #>o.@`>б>g?XxQ3?7K?,Ma?UKVs?|?d 9?vj7?e4^?@?k ?NyV?zMͷ?7? [/`oܱ?zf?Umȡ?Ä?uo?/(w?Y|Od?xP?m7?YEV? ?}Q)*$)x-?e*̣_4+i4g+- k,w:g-iܱ-*.w_{.C)bI/;)/Fr#_0f80gj1R1!R Gm2ӱa2ɴ\Trte3B~3hT4 !!X-4Z:5*85D 66-˶H6H.6- PP77klr8v s8ye8Bi2*9#"L9*{90?*:yovz:s:NK;Ff__;O|;ITnp;r2X-\B>j0k>C@\>AÕ>?B>z~A>?I #+?؛sC?qVs{X?QAl?r|? i`G?5A{ȗ?!i?d%f?OX=߬?%?E(rO)?"hR?tT?'5F5?>@?[0{?TQ~?Wn?[QZ?[:qE?>.?a(}w(+vQ)b:t)8ҫ;>d!gc>M/x> qp>>PFx>hI>g8f ?%d#?@/ ;?γ-{QQ?I@c?-%*t?6ަ蓼?v'?h`$+z?a?TyϚ?s p⸥?s p⸥?fyϚ?a?`$+z?Ev'?Zަ蓼?^%*t?@c?-{QQ?/ ;?PKN3C0XXp1.npyNUMPYF{'descr': '?(( D?:K #BJ?۵ӕQ?}h(U?I@b[?5}a?6te?eej?Rkp?ΙhRs?>\w?-|?4 ?ƒ 5?I ?ɍ)?mB&._?t8?@Ћ6?nI"&?Rk? 9w??k粧?7E՘?LgK0ٙ?!c;I?#܉u,?)l_r?)l_r?#܉u,?!c;I?>gK0ٙ?E՘?k粧? 9w??2k?nI"&?&Ћ6?58?4B&._?ɍ)?I ? 5?  ?-|?N>\w?LΙhRs?bRkp?Feej?se?:5}a?FI@b[?n}h(U?ӕQ?J #BJ?)(( D?$AL>?p 6?o|0?vA(?y1!?{Z?zodBO?|A ?s~?x>hx4>fL>V7>wn}>aو>LN>%Z->Ƭ>=5 *>U蕴>#~*?aS>2I2?$u8?յ@?\2D|F?6N.M?lr}S?؉FY?,*:`?ˀd?u3i?S;S-p?Zs?_x?"Iy}?pRܤ? ։%?kЊ?fL?[ Cf?%"?ˬ? !H~?`"[?8ml?lqʞ?Xq!?d$?9n?T99/?(.Oģ?o^m+?^?^?j^m+?.Oģ?T99/?9n?S$?Xq!?Eqʞ?8ml?`"[? H~?lˬ?"?[ Cf?)fL?jЊ?Չ%?2I2?#~*?Pd#?c-?{l'?`k ?wpu?XP?>j^:>rD#>Pj>4sS>`E/o>=s>P[>^̒v>7Pڠ>ı3?e:?)QB?uv8H?q>P?}h(U?K\?Fb?YPg?Qn?"os?>\w?>DS}?̨} ?B?_ ?a?u:?^?Jh?El_r?cp6?bGE?m?B?qC?:bp@D?h,Zeĩ??bU1?Iej ?G!P?~9,?~9,?@!P?Iej ?bU1??yh,Zeĩ?"bp@D?tqC?B?m?bGE?cp6?l_r?xJh?^?u:?`?i_ ??˨} ?CS}?7>\w?os?Qn?mYPg?b?ʂ\?Y}h(U?/>P?u8H?v)QB?޽:?Rı3?YV,?~])#?GR%?8&?s ?U[J?*[>^> ;Y> [_> !2>]Z7>cb>ϳ>Oq?Cè>xd̢?F$?@JD|?1=?Va\?đ&ݒ?dmJ?t[*Q?c#V?|clK?+U.{?d~v?;CaYq?Dpy5k?d?;_?bq X?ᾈQ?{WJ?Wg،C?뒅dkLB>މ^^>iɥ>HXa>ܴ >Ӽ>(R>rQ(5sD?3ҿ)L?pƔ32S?FWY?5Qd(Qa?H ;f?n? Wvs?!$wx?{I?vF?x?a?H"Œ?녲?>܉u,?,?zC?h? LA?{*iì? 7??"?WAj?&i?]O=n?j?Nr/>^>tgl>*q>PE>1!>CZ*M?kT?B\sa[?Gczb?Ї{h? ûVp?[(l;gu?+U.{?'-Bҁ?Xjħ?fL?=Ñ?/?.Ě?6Um/?p`?b??IM#P?_W"?f'ǂ?)O< ?|2й?r-?n{7?Mވ?8];c"?XFuA?z3c2>?\?o? 9&?1Cڕ6_?1Cڕ6_? 9&?؏o?\?k3c2>?XFuA?&];c"?Lވ?n{7?r-?I2й?O< ?f'ǂ?6W"?L#P??b?p`?Um/?'.Ě?/?=Ñ? fL?\Xjħ?{'-Bҁ?+U.{?|[(l;gu?ûVp?0Ї{h?LGczb?ZB\sa[?ZkT?Z*M?NмD?M1=?K5?o~Y-?rbE$?s??\3?~+ ?U?>QV>od>t>SqL>6ځ>>4T?3 [\?h6Tc?Sj?K21Ioq?:w?`ڞ8~?Ҏ혃?ɍ)?6?.x!?r?dp6?zC?htzڦ?>ZSK\?%7?g`?MT? ]O=n?pL ?S?ܖ9?"Y$?"l5 ?)?lSĖ?u[X0?X5?XԿ?=jħ?- ,{MF?gpĖ?gpĖ?& ,{MF?=jħ?IԿ?X5?u[X0?QSĖ?)?l5 ?"Y$?ܖ9?S?OpL ?\O=n?T?7`?7?YSK\?ļhtzڦ?zC?cp6?r?-x!?46?ɍ)??Ҏ혃?`ڞ8~?P:w?21Ioq?.Sj?:h6Tc?Q3 [\?k>4T?;0M?E?=?U}4?,?|S#?ꬕX?4vں?Մ0?( >]t)>*>޴֎>/v>Rͷ!>V1W5\?JḂc?py5k?avc=r?_x?^) ?Xɪ?%]e? Y?a\??$?m$?$d?=Elbk?T??T?L7R䙺?R0?֮rv?3c2>?ov?#lAPz?ʠeGE(?:[;?Z?V?@Y??d=?0kZm?ͺ*?h?h? ͺ*?0kZm?d=??-Y?V?Y?[;?eGE(?kAPz?Wov?\3c2>?rv?vR0?7R䙺??޳?|T?Elbk?m$d?0$?O?$??a\?ĽY?]e?Xɪ?c^) ?H_x?vc=r?)py5k?ขc?0W5\?OZT?>cuM?FͲSD?x\:&>*2%>U9>>~C>naXc?l-snk?IFJr?-+>y?>#ހ?|gkQ?WP`#ހ?.-+>y?lIFJr?l-snk? aXc?\?bT?$?1L?(( D?b*;?i2?KO7)?l0!?-O?}\4 ?L=ʀ?*[>(>xT9>-̯G>W@lKk?(UIr?Qԟy?clK? )?vݔ?:x?`"[?ֹ'/?`>o@3?;[?m2?g'ǂ?(qP ?{TLK?=];c"?x _~?hy5j9? rN?c?.1?dtf{?8ᆎ(?-w,?6 6?X~å?!?^y|B?DK^?S&˵?~oQ?*ޛ?0?0?%ޛ?~oQ?I&˵?6K^?Ly|B?q!?0~å?6 6?-w,?ᆎ(?>tf{? 1?c?O rN?*y5j9?Jx _~? ];c"?TLK?׆qP ?f'ǂ?32?;[?>o@3?'/?N`"[?w?>vݔ?R )?bclK?jQԟy?TIr?Ψ@lKk? cmc?nK: [?ԞyBS?=:wK?íB?]9?2ؤ1?NƸ'?@C?R?`k ?п6?>(eI>Q>=nr?(Qy?кa?ހoZv?<[?\k?>܉u,?FF?qC?*iì??O?g-Vӻ?V(D;?µ?4?o5i)?}i.%?Ӽ?X+?q?i ?gP[\?7k,c?>b?l ?&k6w?0?vr?P3}W?}ٿ?Le?_,Q?-Lg?-Lg?_,Q??e?}ٿ?P3}W?vr?? k6w?l ?>b?|7k,c?4P[\?q?i ?)+?Ӽ?Ri.%?'5i)?4?~?(D;?gg-Vӻ?O?|?*iì?NqC?YFF?ۉu,?k?$<[?ZހoZv?ya?(Qy?L=nr? j? yb?"FbZ?}G˕R?KI?X'A?=,7?Zy.B/?!Dմ4%?n+?^JΞ ?ſtY?E,>z$>2za03>ǡx?-[^?J}eX?mJ?ˬ?ʤ:ܛ?ګ??uު?rc?CP?O?{TLK?uFů?ov?4q?'}?14r,>= ?h_Ć?qr`/ ?^?j%P?s:N?Rbp@D?,x?fj4?+G"?v*?"Y$?lF??Y?ٽY? 7?q?i ?L0? ē?آ<k?0?\w?6(p?>YPg?nPT`?Db'!V?ooO?EKHE? Zą?&JN?-@[o?I(?9n?dc?T?$?L7R䙺?Mވ?UG?,]~?fh[??>ᆎ(?; t?;D$?Q ??9H_?* ,?3^? $qB?466D?J?BKlFL?ع5=?RF?tD%?0q?7?~Ir?74@`ڇ@b@b@`ڇ@74@jIr?ֈ0q?7?ZD%?4F?5=?%KlFL?bJ?l466D?#qB?^? ,?H_? ??D$?F; t?(??fh[?r,]~?TG?Lވ?6R䙺?X$?MT?dc?9n?(?ȣ@[o?p%JN?Zą?-AZ2L?T@Fv?seo?5FGe?wB^?sOZT?ƎL?-?B?wu8?-Ma0?5NdO%?tګs??u?zgX>mB&._?1)?;c;I?6&?$Sp?c&3ɯ? "Ҵ?j?XV?.l5 ?s?h38*H?-J?Q?q?i ?:4Y[}?4"?KaW? @8?F_Av?qJ8?M?%"?|/?u?D&|Z^?oVr@_@vmQ@gI.@< Q~@C()@"ע.@Ld]@Jd]@ע.@C()@< Q~@gI.@fmQ@p_@ZVr@&|Z^?du?{/?i%"?M?qJ8?^Av?@8?aW?"?ӄ4Y[}?q?i ?-J?Q?,38*H?s?l5 ?XV?Xj?5 "Ҵ?wb&3ɯ?$Sp?6&?b;I?)?A&._?5h4e?xfC }?!vt? l?fc?EZ?2$CR?bv{H?$z]c?@?< Z5?i+?y1!? O?E"A ?s~?/ }mN? 8G5?"/?s2?+ɺ?g'ǂ?1U޺?N[l?ACڕ6_?up?V?:;7M?ɉU? _?~oQ??S?6Q0??F[٧?mN?`c?e7?74@@9P@" @Jn9@ X@V-,@Jsi @\r @#&A @# @_q @_q @# @#&A @\r @Jsi @vV-,@ X@0n9@" @@9P@74@le7?0c?:N?[٧??w6Q0?g?S?~oQ? _?@ɉU?s:;7M?vV?Pup?Bڕ6_?sN[l?U޺?f'ǂ?ɺ?Xs2?6"/?7G5?a/ }mN?FEA^R?+?LmÛgz?#elr?8;si?@#* aga?yW?E3O?.мD?KT;?1?ES'?d-q?g[#?;?Hzc?,?h?d!P??w-γ?ŏ|4?۞Z)? B?X5?)>9*?X+?n7~?[?,k6w?!w?` ?7૙.k?^f?Q? |/?ȟ?qVr@Gi@Υ@Rl8@ɘ @5Y @l.+ @B1%@;jB;@PnY@ i@ܪ@~$@|$@ܪ@i@EnY@;jB;@&1%@K.+ @5Y @fɘ @0l8@Υ@Gi@RVr@Eȟ?{/?^f?Q?૙.k? ?w?j6w?[?2n7~?+?=9*?EX5?+ B?Z)?O|4??w-γ?!P??d,?Hzc?5nА? Ό?>#ހ?-Z?Iw?9Ap?+OIb? 9J?XP?I_?Mt?Sȉ?+{2-?=@Hs@Xl8@:HA" @&Pێ` @9xZe@²E@@\z@veJ @@u@?^@*Ҭ@q+V@"@^[k@\[k@"@d+V@Ҭ@?^@@u@\eJ @\z@@E@9xZe@u&Pێ` @ HA" @%l8@s@=@>+{2-? ȉ?Mt?I_?WP?h 9J?؀>b?ƀ?91?/›_?^u[X0?kb?]%_?*`y?*d?eY‚?:N?b;I?\[?H4?Ү??By?HȈq?&YPg?> _?B\*7U? FK?R6^p6#B?z\\7?]GF -?W"? d?b?b=???te5?r?<ݻ?c?ͺ*?:ܞ? e[?\Jyt? ?$qB??r}?B? fb@On9@l/@%=] @,'@Fa:6@feWT@՝zϯw@UY@dt@ ӏ@mC S@ "3@ώ @M!@&Wq!@x`8!@x`8!@!Wq!@M!@ώ @"3@SC S@ӏ@dt@UY@zϯw@@eWT@Fa:6@,'@=] @l/@n9@eb@A? }?@?#qB? ?Jyt?F e[?T:ܞ?ͺ*?kc?ݻ?ar?te5??'b=?hb?EXq!? H~?>|t?1 )? ?ہu?0cm?ขc?+NZ?RlHQ?.\2D|F?Zx\x2?VQ'?>E.&?bU1?i?l jw?=E?)?gpĖ?jmY?`4fo?U@C:??Ze?h,?F?/?}T;?Ad1?s~"%?^@?N"@*$@%P@&@ q'@T)@t8+@c,@wW-@.@j@O/@~PO͟/@{PO͟/@c@O/@.@wW-@ߡc,@\8+@T)@q'@v%P@&@*$@qy>"@zHSX @AT@ ԃu@M ("@0k @;jB;@:@t0d @ l8@]f@O8 ?{/?TD;'? ;1@ ^G2@?`q2@Ĝ^ 2@Ĝ^ 2@?`q2@ ^G2@> ;1@Et+ 1@`60@$'QZ.@~x,@ސ9uH*@7zX(@ (%@px#@HtpD!@4C S@ ]@}@.ͅV@UP@V#&A @:qщ@>/s@NG?@?Ah"?fk?H_?y|B?ɉU?Td=?}^?c ?Lވ?R?b?‡!?I֡??۝?I}eX?B y?>}0u?j0)]9 l? Gczb?bq X?& N?t%C?8u8?*/?pL ?tiA ?u[X0?ٽY?t%0?@Ln9?1@?Ze?<L?%"?5rP?_@Zd]@u]N @”f@PF=@(~Ĝ@0Ҭ@]@5*b!@Q7]$@֝&@T)@k0T,@8.@q50@r2@Lk3@r4@5u4@x=5@yˀ5@176@176@tˀ5@x=5@5u4@r4@;k3@|r2@\50@x8.@k0T,@ЪT)@&@fQ7]$@ *b!@@Ҭ@~Ĝ@fPF=@f@]N @d]@H_@T5rP?<%"?L?e?N1@?@Ln9?6t%0?ٽY?#ހ?}MEmcv?|=km?#]e#c?H ǦX?ooO?h{C?M8?6Cڕ6_?2[;?9!?cPUY?DK^??" j?J?Z;D$W?`@K@#&A @[@wY@И@Dl:@x`8!@lGʬ$@\zX(@Un+@}8mw/@> ;1@jG3@kDk5@Y[7@$h9@C;@ %z=@6?@,U)@@iqU@@a̫u@@ɜ%,)A@ɜ%,)A@a̫u@@bqU@@,U)@@6?@%z=@(;@h9@Y[7@IDk5@GG3@> ;1@F}8mw/@PUn+@zX(@9Gʬ$@w`8!@l:@И@VwY@[@K#&A @ʴK@`@;D$W?(J?:" j??K^?-cPUY?89!?Z;?Bڕ6_?ޓm?oP T?Wzh"?6XE`?`d?\'D?Կ}̌?ޱR5FGe?wӐ6\?ᾈQ?{b̶F?yd0k @>gD@"@ekX @p]t¿#@<5Qm'@z8+@v@O/@#C1@r4@NT/?p6@_8@h:";@ Vd=@]?@Ab@@;bA@6PB@:B@^@C@,rC@,rC@^@C@:B@6PB@-bA@Ab@@_]?@ Vd=@h:";@_8@MT/?p6@r4@h#C1@1@O/@38+@5Qm'@>]t¿#@7kX @s"@gD@0k @1%@@gI.@38 ? ?<?q-Lg?j6w?L? v?+|Jt?eY?Oܖ9?K?]H?JY‚?(FF?Jh?îRQ?vmbI?LfC }?os?,h?d> _?bT?>KI?nk)@?up?'@|@@Z+gA@-&uB@HYC@D@H­+E@ 4?E@#E@#E@ 4?E@@­+E@D@HYC@&uB@F+gA@|@@XTp>@;@SPcQ,9@4c; 6@}s%:4@WM 1@'QZ.@9uH*@~ e&@^=?"@,=ѐ@@ v,_@bR@6# @#@`@F?yVP?BNy?tĝ~I?y|B?=-w,?@sB@@6PB@r%C@y1D@E@gFF@nɃG@ϧdG@]~;H@]~;H@ϧdG@nɃG@gFF@E@f1D@r%C@6PB@sB@@)Џ>@9;@nc?- 9@176@$P&΁3@Ʀh0@F^1-@5ܮb(@bc$@0޾!@AT@@VPF=@1%@ ߜ@Y[M@Y ?I5.?L?4}ٿ?zhK?sc?klA?fgpĖ?~dĵ?v*?r?Bȶ?@b&3ɯ??;?sCU? cB?/M5?>\w?*QHn?ϩaXc?͓fY?@i"P8 P?'( D?GfW^?<.?j@N?"??=tLj?žu?f}?`ڇ@En9@+ @( @ĥ(@ @$"3@1F#@g:'@f:. [+@B^4 0@?`q2@t$I|s5@ m8@ h4;@hK->@a̫u@@o[ B@bn?o#D@ 4?E@c* (F@rR-H@eI@`g۵I@.VgoaJ@oOJ@oOJ@'VgoaJ@Vg۵I@eI@rR-H@c* (F@ 4?E@In?o#D@o[ B@a̫u@@/K->@| h4;@ m8@t$I|s5@?`q2@^4 0@":. [+@Eg:'@k1F#@"3@@pĥ(@ @Ɩ @m9@`ڇ@}?ku?t=tLj???Aj@N?į.?fW^?'h\)?^rv?qP ?*jz?cYD%?$Xq!?a\?%%JN?Ɵ!ƿ?fY$z?Loq?H5FGe?K: [?c )@wW-@y1@hM4@. Z7@:@=@Ab@@tTB@JG=D@E@nɃG@tU7H@N>&?J@A,QK@/.&L@H\"L@MQ,zM@MQ,zM@G\"L@$.&L@A,QK@ N>&?J@^U7H@nɃG@|E@JG=D@tTB@Ab@@\=@½:@Z7@hM4@oy1@wW-@c )@Xc$@2ւ`[!@ [k@ eJ @|@/&Pێ` @.'g@'}J@u?>/?@P@}1?A@^@C@$XFE@u-8GG@D8I@vf K@H\"L@깎1N@pmO@W1P@R(սP@vWA=P@tWA=P@N(սP@W1P@ pmO@޼깎1N@G\"L@uvf K@hD8I@W-8GG@$XFE@^@C@0@}1?A@ٗ>@ :@Z7@vr4@+P~L !1@0i!,@f+$(@']t¿#@zHSX @@C@ nY@]N @Υ@ %~?7SD^?߫.k?¯e?X<k?2E'3?flA? ,{MF? B?H]%_? WAj?:*t9@176@)(8 82@7+/@iA%*@B%P@&@byƠe"@XB;H@zXl@Ǧb@@(&Pێ` @l8@&_@.:"?t[?^Av?$ 9J?U?^C ?28*H?J4?SXV?j;0?G~;?^htzڦ?cp6?k?pV?/M5?Y?Iw?Qn?g6Tc?sgX?N?X1?%D$?#H_?N?j? u?vI3D@jF@z޴(pl@MͅV@LvVV@^C S@S"@a9p&@zUn+@b70@֚3@^+7@NϺ:@]K->@]A@HYC@m F@z]5ZH@ oOJ@{1pL@\q;N@x(oP@o%aGQ@i&0Q@ER@^BR@S@S@^BR@>R@_&0Q@o%aGQ@i(oP@\q;N@V1pL@nOJ@z]5ZH@E F@HYC@v]A@K->@ Ϻ:@$+7@֚3@b70@ Un+@a9p&@S"@C S@uVV@ͅV@ ޴(pl@ F@"I3D@t??N?H_?ВD$?1?DV?lkAPz?>2?qP ?M?uު?b53?`={?eL?t+?_x?Vseo?Qd?؉FY?D3O?V:h?[? ?hz:a,?^f?Q?KL?jmQ@p_1K@Wv1@@[u@G@ġZ.#@.5Qm'@k0T,@Ʀh0@r4@ߏԹ–7@ܥ h;@r]?@wzaA@JG=D@A F@nU7H@A,QK@L_M@' DܽO@?*P@!\Q@D4spR@1 5S@qbybS@/ YS@/ YS@qbybS@* 5S@:4spR@!\Q@?*P@w' DܽO@'_M@hA,QK@GU7H@ F@JG=D@ƦwzaA@,]?@zܥ h;@Թ–7@pr4@Ʀh0@Mk0T,@4Qm'@Z.#@@[u@Z@)Wv1@ _1K@mQ@L?T^f?Q?:hz:a,? ?[? V:h?~ٽY?m?Ph?*G"?I; ?bp@D?8 n?Y?^O`ʷR@EIS@ #xS@S@S@#xS@EIS@3ʷR@Z&0Q@o7:^Q@>npbUP@H N@jNjK@3힡XI@c* (F@D@P?k.)?n e[??S?t" j?*X?}e7?9@l/@QJ@s~ @I:U@t=ѐ@px#@Kg)'@g ,@471@wg4@Y[7@/;@*@@[r8B@D@c* (F@V힡XI@NjK@H N@LnpbUP@o7:^Q@c&0Q@:ʷR@EIS@#xS@S@S@#xS@EIS@/ʷR@V&0Q@o7:^Q@:npbUP@H N@eNjK@/힡XI@c* (F@D@9r8B@*@@;@qY[7@Pwg4@471@ ,@Kg)'@lpx#@=ѐ@:U@s~ @&QJ@l/@9@d7?K*X?&" j??S? e[?Tk.)?'?,]~?juFů? ?W"?!K#?P'/?۫a\?PO)}?c?-w^y?xûVp?d?X+NZ?bw$qP?V:h?[? ?hz:a,?^f?Q?>L?amQ@c_1K@Wv1@@[u@6@Z.#@%5Qm'@k0T,@Ʀh0@r4@֏Թ–7@ܥ h;@f]?@wzaA@JG=D@7 F@dU7H@A,QK@?_M@' DܽO@?*P@!\Q@=4spR@* 5S@qbybS@/ YS@/ YS@qbybS@" 5S@24spR@!\Q@?*P@j' DܽO@_M@^A,QK@=U7H@ F@JG=D@wzaA@ ]?@pܥ h;@Թ–7@hr4@Ʀh0@Bk0T,@4Qm'@~Z.#@@y[u@P@Wv1@_1K@ mQ@L?P^f?Q?6hz:a,? ?[?V:h?yٽY?m?Ph?*G"?I; ?~bp@D?8 n?Y?PO`@]A@HYC@_ F@z]5ZH@nOJ@h1pL@\q;N@m(oP@o%aGQ@^&0Q@9R@{^BR@S@S@w^BR@2R@T&0Q@o%aGQ@^(oP@\q;N@E1pL@nOJ@z]5ZH@7 F@HYC@k]A@K->@Ϻ:@+7@֚3@b70@Un+@|a9p&@S"@C S@uVV@ͅV@ݴ(pl@F@I3D@t??lN?H_?ÒD$?1?7V?XkAPz?|>2?mqP ?M?uު?W53?`={?eL?a+?_x?7seo?Qd?؉FY?D3O?C ?U? 9J? _Av?[?F/:"?k_@3l8@}&Pێ` @Ǧb@@Xl@'YB;H@yƠe"@v%P@&@ZiA%*@(7+/@J(8 82@176@k*t9@Wް[b=@sB@@vC@$XFE@nɃG@2ƱI@k"GK@ZM@' DܽO@jWA=P@v]Q@NnQ@< o{AR@44spR@34spR@< o{AR@zNnQ@v]Q@^WA=P@g' DܽO@ZM@J"GK@2ƱI@nɃG@$XFE@vC@sB@@ް[b=@)*t9@176@(8 82@Η7+/@iA%*@0%P@&@UyƠe"@XB;H@hXl@Ǧb@@&Pێ` @l8@_@.:"?^[?x^Av? 9J?U?FC ?28*H?24?BXV?j;0?G~;?Ghtzڦ?bp6?k?pV?/M5?Y?Iw?Qn?g6Tc?sgX?N?L?(~å??/T^?u?7c?mdX@K@Ԋ_q @F@zϯw@2 @HtpD!@Ŷ%@0)@'QZ.@V22@#<5@c&=8@a<@,U)@@Sr8B@R?WD@HTzF@ݘH@nOJ@oQQL@YUMXN@iO@x&-P@o7:^Q@cQsQ@K\Q@K\Q@cQsQ@o7:^Q@x&-P@iO@YUMXN@|oQQL@nOJ@ݘH@%TzF@R?WD@2r8B@e,U)@@la<@#&=8@#<5@)22@'QZ.@s0)@ӘŶ%@HtpD!@A2 @]zϯw@@n_q @K@dX@ýc?Ru?K/T^?e?~å?ى?&}?ۼڄ?& Ғ?O< ?S?=+ ٥?n?Sԓ?jFEA^R? clK?(v?kl?Fczb?Q+PUW?Z*M?lA?E'3?<k?#e?૙.k?SD^?{%~?Υ@C]N @8nY@C@@zHSX @O]t¿#@+$(@0i!,@BP~L !1@r4@ Z7@ :@3ٗ>@;@}1?A@^@C@$XFE@X-8GG@hD8I@rvf K@G\"L@м깎1N@pmO@W1P@>(սP@bWA=P@`WA=P@:(սP@W1P@pmO@깎1N@G\"L@Uvf K@HD8I@:-8GG@$XFE@^@C@@}1?A@ٗ>@ :@Z7@^r4@P~L !1@0i!,@I+$(@]t¿#@zHSX @@C@mY@\N @mΥ@%~?SD^?߫.k?e?><k? E'3?KlA?[ ,{MF? B?!]%_?WAj?:힡XI@7J@KL@ ;ͶrM@YUMXN@\q;N@ DO@ DO@\q;N@YUMXN@;ͶrM@KL@J@#힡XI@^G@#E@`#C@&aKj B@^,U)@@4<@SPcQ,9@ Dk5@Ĝ^ 2@^4 0@m^+@<&@S"@_4]@nИ@]GI@ QJ@rVAA@~ f@|NG? ,?#qB?"R?PS"?9ܞ?.a? rN?é?|8/f?^f'ǂ?U?p`?>1=?lY?I}eX?zKh~?Zs?)ci?nc`?!OZT?b_I?#?q?i ?ē? ?=B.??/?Vu?T}J@\'g@n&Pێ` @|@JeJ @8[k@Lւ`[!@|c$@c )@wW-@~y1@hM4@ Z7@Ƚ:@c=@Ab@@tTB@JG=D@tE@nɃG@MU7H@M>&?J@fA,QK@.&L@G\"L@LQ,zM@LQ,zM@G\"L@-&L@UA,QK@M>&?J@6U7H@nɃG@ZE@JG=D@tTB@Ab@@,=@:@Z7@hM4@Ry1@wW-@κc )@8c$@ւ`[!@[k@eJ @|@&Pێ` @'g@ }J@t?>/?@a̫u@@o[ B@?n?o#D@ 4?E@c* (F@rR-H@eI@4g۵I@VgoaJ@nOJ@nOJ@UgoaJ@*g۵I@eI@rR-H@c* (F@ 4?E@&n?o#D@o[ B@qa̫u@@J->@N h4;@ m8@t$I|s5@b?`q2@]4 0@:. [+@g:'@I1F#@"3@@Mĥ(@ @ @m9@x`ڇ@}?Iu?@=tLj???%j@N?.?љfW^?g\)?>rv?TqP ? jz?cYD%?Xq!?īa\?$JN?!ƿ?vfY$z?wLoq?(5FGe?J: [?@sB@@6PB@r%C@P1D@jE@gFF@nɃG@ϧdG@]~;H@]~;H@ϧdG@nɃG@gFF@ZE@<1D@r%C@~6PB@sB@@묃Џ>@9;@\w?QHn?aXc?P͓fY? i"P8 P?'( D?eup?@|@@0+gA@&uB@HYC@D@­+E@ 4?E@#E@#E@ 4?E@­+E@֭D@HYC@%uB@+gA@|@@Tp>@;@SPcQ,9@b; 6@Ps%:4@dWM 1@'QZ.@9uH*@N~ e&@Y^=?"@=ѐ@@m v,_@<R@# @"@`@F?JVP?rBNy?8ĝ~I?y|B?,w,?@1%@0k @gD@"@=kX @>]t¿#@5Qm'@88+@,@O/@^#C1@~r4@MT/?p6@_8@jh:";@ Vd=@5]?@Ab@@bA@6PB@d:B@^@C@r,rC@p,rC@^@C@[:B@6PB@bA@Ab@@]?@q Vd=@Hh:";@k_8@MT/?p6@Xr4@=#C1@@O/@8+@4Qm'@ ]t¿#@kX @1"@gD@/k @1%@ޗ@gI.@7 ?Q ?<?*-Lg?j6w?ZL?T v?*|Jt?JeY?"ܖ9?K?a]H? Y‚?EF?ªJh?îRQ??mbI?fC }?Zos?,h?> _?TbT?KI?jnk)@?Bڕ6_?Z;?x9!?tcPUY?K^??[" j?MJ?;D$W?`@K@g#&A @[@dwY@И@l:@w`8!@9Gʬ$@"zX(@HUn+@3}8mw/@> ;1@<G3@4Dk5@Y[7@h9@;@ÿ%z=@6?@l,U)@@AqU@@za̫u@@Ȝ%,)A@Ȝ%,)A@va̫u@@:qU@@c,U)@@h6?@%z=@;@Ƈh9@hY[7@Dk5@G3@> ;1@|8mw/@ Un+@yX(@Gʬ$@w`8!@l:@hИ@%wY@[@#&A @K@`@;D$W?J?" j?m?K^?bPUY?9!?KZ;?Bڕ6_?ޓm?

UB?%|Z^?}0@,@ @W&Pێ` @V$@!ܱ@@o-@eM@ryƠe"@hc$@{+$(@O'~+@I8.@_QM1@k3@\5u4@:6@bN8@>*t9@\h:";@sČ=<@jZ=@ް[b=@6=@4=@ް[b=@_Z=@dČ=<@Hh:";@(*t9@HN8@:6@A5u4@k3@_QM1@8.@'~+@E+$(@3c$@EyƠe"@TeM@o-@!ܱ@$@%Pێ` @? @L}0@~%|Z^?TB?<?n_?`?U?ʇI?^ٽY?F}?X '?g?&"?aU1?Tʩ?<ۉu,?{ǁA?ȍ)?A>#ހ?:MEmcv?%=km?]e#c? ǦX?ooO?ph{C?SM8?Lވ?H?{lM?R0kZm?S{Ǣ? K^?K ?ĺ?,d֩?&D?`ڇ@/7@Jsi @f @Hm6@zϯw@idt@=i@@Mέ'"@td$@Kg)'@ҭa*@i-@b70@322@֚3@-5@bS6@6ugx7@ m8@SPcQ,9@t9@ʇh9@Ƈh9@{t9@SPcQ,9@ m8@(ugx7@bS6@-5@֚3@22@sb70@i-@a*@{Kg)'@ztd$@ έ'"@=i@@'dt@Ozϯw@Hm6@f @SJsi @9/7@h`ڇ@p&D?ޤd֩?ĺ?K ?K^?S{Ǣ?/kZm?lM?H?Lވ?6O?M?dc?-B?, H~?1Ux!?(Zą?,.΅.}?Br]s?8;si?]q`?W3RU?1K?)X*A?7$Hou5?#pL ?3iA ?zu[X0?ٽY?Mt%0?@Ln9?o1@?e?L?L%"?g5rP?N_@d]@]N @f@bPF=@޶~Ĝ@ޯҬ@@*b!@LQ7]$@&@T)@bk0T,@>8.@:50@Tr2@k3@jr4@O5u4@tx=5@3ˀ5@176@176@.ˀ5@mx=5@B5u4@\r4@k3@Br2@&50@8.@5k0T,@~T)@b&@$Q7]$@Ս*b!@@Ҭ@~Ĝ@*PF=@\f@\N @c]@_@5rP?%"?L?e?1@?@Ln9?s%0?cٽY?t[X0?ЮiA ?oL ? "Ҵ?C!P?B?ܨ:$?P-x!?րpV?~??'#wx?ixp? ;1@] ^G2@q?`q2@Ĝ^ 2@Ĝ^ 2@l?`q2@U ^G2@> ;1@t+ 1@`60@'QZ.@"x,@9uH*@yX(@ (%@[px#@HtpD!@C S@ ]@}@ͅV@"P@#&A @qщ@/s@|NG?@?h"? k?H_?y|B?ȉU?d=? ^?c ?Lވ?R?b?U‡!?֡?`?۝?I}eX?y?B>}0u?0)]9 l?Fczb?`bq X?o& N?t%C?u8?*/?@?'N"@*$@D%P@&@q'@T)@8+@c,@wW-@@.@@O/@PO͟/@ PO͟/@@O/@4.@wW-@xc,@8+@T)@q'@(%P@&@*$@Qqy>"@zHSX @AT@ ԃu@ ("@/k @D;jB;@?t0d @l8@f]f@7 ?j{/? D;'?/?/?h,?}e?`? U@C:?4fo?lY?gpĖ?)?=E?pjw?i?aU1?FSD?1aڜ?vsCU?`A&._?Ҏ혃?HPz?Xq%r?Dq3ZTh?@r%`?[*7U?/ʚK?|XEA?Ap 6?@,?ny1!?b?Fb=?:?te5?fr?ݻ?tc?ͺ*?L:ܞ?B e[?Jyt? ?#qB?0?}?A?eb@m9@l/@=] @},'@Fa:6@eWT@zϯw@UY@Sdt@ӏ@C S@"3@jώ @KM!@Wq!@w`8!@w`8!@ݥWq!@EM!@`ώ @"3@C S@}ӏ@1dt@dUY@]zϯw@eWT@^Fa:6@_,'@w=] @rl/@m9@eb@8A?}??#qB?G ?Iyt? e[?9ܞ?nͺ*?(c?aݻ?r?Ate5?~?a=? b?Xq!?, H~?|t?Գ )? ?įہu?0cm?pขc?J+NZ?lHQ?[2D|F?x\r9>x2?YVQ'?ʆE.&?:N?qY‚?*d?*`y?]%_?kb?^u[X0?*›_?(1?ƀ?π>b?Q 9J?WP?vI_?Mt?lj?+{2-?~=@~s@l8@GA" @,&Pێ` @9xZe@rE@\@a\z@eJ @@u@'?^@Ҭ@+V@J"@[k@[k@B"@+V@Ҭ@?^@r@u@eJ @F\z@D@WE@}9xZe@%Pێ` @GA" @l8@~s@T=@+{2-?lj?HMt?DI_?/WP? 9J?>b?~ƀ?1?_?t[X0?kkb?]%_?*`y?x*d?Y‚?:N?1b;I?]\[?G4?~??By?GȈq?YPg?= _?[*7U?FK?6^p6#B?5z\\7?GF -?W"?d?ܟ?6^?΍&ƫ?cKXپ?qP ?\ Ғ?o?DeGE(?:$?? ; t?y|B?F;?K~W'?:" j?[٧?JlFL?;D$W?#ހ?Y?Iw?=9Ap?\+OI6Q0?^?[٧?N?c?e7?h74@p@9P@" @m9@ X@V-,@rJsi @:\r @##&A @# @n_q @k_q @# @#&A @*\r @`Jsi @V-,@ X@m9@" @U@9P@J74@d7?c?M?[٧?4?5Q0??S?Y~oQ?P _?ȉU? :;7M?%V?tp?Bڕ6_? N[l?2U޺?Ef'ǂ?ɺ?r2?!/?7G5?/ }mN?BFEA^R?N+?֦mÛgz?#elr?.8;si?"* aga?yW?D3O?ƘмD?S;?1?IES'?d-q?[#?;?A&._?)?b;I?6&?$Sp?_b&3ɯ?* "Ҵ?0q?7?Hr?P74@v`ڇ@3b@1b@r`ڇ@H74@Hr?(0q?7?C%?F?(5=?JlFL?J?366D?#qB?G^?Y ,?yH_? ??D$?: t?e(??[fh[?+]~?>TG?yLވ?N6R䙺?$?S?dc?19n?((?N@[o?$JN?Zą?rAZ2L??Fv?seo?4FGe?-wB^?OZT?9ƎL?P-?B?ᣘu8?ӠMa0?MdO%?٫s???VzgX> ?^Ć?q`/ ?r^?4j%P?:N?ǽbp@D?,x?j4?*G"?2v*?"Y$?lF?V?Y?ٽY?7?Eq?i ?vK0?ē?b<k?>hz:a,?L?&,+?hMt?/p??z?TB?lj?{/?诠xh]?䯠xh]?|{/?lj?TB??z?x/p?RMt?,+?L?&hz:a,?;躽D?B@8?>\w?'p?YPg?PT`?a'!V?booO?EKHE?ǡx?-[^?J}eX?mJ?,ˬ??:ܛ?(ګ??Juު?Vrc?BP?O?TLK?uFů? ov?J4q?&}?4r,>:=nr?'Qy?ga?CހoZv?<[?k?ۉu,?:FF?qC?)iì?T?vO?g-Vӻ?(D;?F?i4?4i)?i.%?NӼ?*?8q?i ?P[\?77k,c?>b?ql ?j6w??evr?"P3}W?}ٿ?e? _,Q?9-Lg?9-Lg?_,Q?e?|ٿ?P3}W?Lvr??j6w?Yl ?>b?7k,c?P[\?q?i ?*?Ӽ?i.%?x4i)? 4??(D;?f-Vӻ?'O??i)iì?pC?EF? ۉu,?k?h<[?݀oZv?a?~'Qy?=nr?$ j?yb?EbZ? G˕R?KI?'A?Z=,7?x.B/?Dմ4%?+?JΞ ?9tY?TD,>xz$>ya03>@lKk?TIr?QQԟy??clK? )?vݔ?w? `"[?j'/?=o@3?V;[?2?~f'ǂ?qP ?TLK?\;c"?w _~?ďy5j9?ڠ rN?Ic?1?sf{?(?"-w,?-6 6?~å?!?y|B?K^?%˵?h~oQ?ݛ?'0?'0?ݛ?^~oQ?%˵?K^?y|B?!?a~å? 6 6?,w,?k(?sf{?1?#c? rN?y5j9?w _~?\;c"?TTLK?5qP ?;f'ǂ?2?;[?=o@3?/'/?_"[?^w?wuݔ? )?clK?Pԟy?PTIr?@lKk?bmc?J: [?^yBS?:wK?GB?6]9?2ؤ1?&NƸ'?M@C?R?շ`k ?`6?>O'eI>EQ>aXc?ܻl-snk?YIFJr?-+>y?}>#ހ?gkQ?O`#ހ?-+>y?HFJr?Sl-snk?aXc?\?>bT?[?1L?'( D?a*;?0i2?N7)?2l0!?O?|\4 ?<ʀ?[>(>JxT9>̯G>0W5\?ขc?py5k?uc=r?_x?L^) ?{Xɪ?r]e?Y?a\??$?$?'$d?Elbk?T???6R䙺?R0?krv?3c2>?nv?kAPz?eGE(?yZ;?Y?DV?X?H?&d=?0kZm?ͺ*? h? h?ͺ*?0kZm?d=?8?X?2V?Y?RZ;?eGE(?QkAPz?nv?3c2>?Crv?Q0?\6R䙺??f?S?BElbk?$d?$??$?a\?ZY?]e?Xɪ?^) ?_x?uc=r?zoy5k?fขc? 0W5\?OZT?=cuM?̲SD?x\:&>*2%>U9>>~C>B>4T?3 [\?h6Tc?Sj?11Ioq?":w?`ڞ8~?Ҏ혃?Dɍ)?6?-x!?[r?/cp6?czC?thtzڦ?YSK\?7?`?T?d\O=n?oL ?S?=ܖ9?"Y$?l5 ?3)?RĖ?u[X0?W5?Կ?=jħ?[ ,{MF?'gpĖ?'gpĖ?R ,{MF?=jħ?Կ?W5?t[X0?RĖ?)?rl5 ?b"Y$?ܖ9?S?oL ?+\O=n?T?`?7?3YSK\?*htzڦ?&zC?bp6?r?<-x!?e6?ȍ)?Ҏ혃?8`ڞ8~?:w?11Ioq?Sj?g6Tc?3 [\?=4T?:0M?E?$=?Ļ}4?@,?v|S#?>X?K4vں?00?4( >t)>>j֎>.v>̷!>hZ*M?(kT?$B\sa[?Gczb?χ{h?ûVp?<[(l;gu?J+U.{?4'-Bҁ?Xjħ?eL?W=Ñ?T/?-Ě?Tm/?2p`?5b??^L#P?W"?cf'ǂ?{O< ?1й?"r-?'n{7?Lވ?\;c"?XFuA?3c2>?[?@o?0 9&?Bڕ6_?Bڕ6_?+ 9&?5o?[?3c2>?pXFuA?\;c"?}Lވ?n{7?q-?1й?SO< ?@f'ǂ?W"?L#P??a?p`?Tm/?l-Ě?/?=Ñ?8eL?Wjħ?'-Bҁ?ۉ+U.{?Z(l;gu?FûVp?jχ{h?Fczb?B\sa[?kT?֤Z*M?мD?M1=? K5?@n~Y-?rbE$?r??ٕ3?}+ ?~U?>|V>nod>6t>pL>6ځ>*rQ(5sD?`ҿ)L?Ŕ32S?BFWY?Qd(Qa?;f?n? Wvs?r#wx?I?F??;`?H"Œ? 녲?ۉu,?D,?ZzC?? LA?)iì?7?9?I"?WAj? &i?J\O=n?j?N/>~^>gl>K*q>!PE>z >d2kLB>Dމ^^>iɥ>HXa>< >lӼ>R>ı3?:?Q)QB?u8H?>P?}h(U?\?b?YPg?+Qn?os?؟>\w?bCS}?˨} ??^ ?,`?au:?G^?Jh?vl_r?cp6?bGE?4m?0B?pC?bp@D?g,Zeĩ?? bU1?Hej ?i!P?}9,?}9,?a!P?Hej ?aU1??g,Zeĩ?lbp@D?pC?B?m?bGE?bp6?@l_r?ªJh? ^?4u:?_?^ ?T?L˨} ? CS}?>\w?Pos?Qn?YPg?fb?󁱤\?|h(U?>P?Mu8H?(QB?:?ı3?xU,?E~])#?GR%?8&?s ?[J?n[>>Z;Y>:Z_> 2> Z7>cb>ϳ>p?Cè>"~*?φS>2I2?Qu8?uյ@?"\2D|F?IN.M?Ɩr}S?؉FY?4,*:`?0ˀd?2i?і;S-p?Zs?_x?"Iy}?Rܤ?VՉ%?2jЊ?eL?[ Cf?"?ˬ?O H~?_"[?8ml?nqʞ?Xq!?$?H9n?ST99/?-Oģ?Ȗ^m+?^?^?–^m+?|-Oģ?IT99/?;9n?$?Wq!?Gqʞ?8ml?_"[? H~?Îˬ?j"?Z Cf?EeL?iЊ?#Չ%?Rܤ?´"Iy}?_x?Zs?;S-p?02i?ˀd?+*:`?؉FY?nr}S?N.M?[2D|F?)յ@?ȣu8?tS>2I2?I"~*?@Pd#?c-?k'?ⷆ`k ?wpu?XP?>]:>C#> Pj>sS>_E/o>i=s>[>˒v> 7Pڠ>y1!?cvA(?o|0?p 6?AL>?'( D?hJ #BJ?SӕQ? }h(U?H@b[?5}a?se?eej?Rkp?͙hRs?̟>\w?-|? ?- 5?8I ?ɍ)?A&._?w8?Ћ6?mI"&?k?3 9w??@k粧?pE՘?}fK0ٙ?Kb;I?Jۉu,?Ml_r?Ml_r?Jۉu,?Kb;I?pfK0ٙ?WE՘?4k粧? 9w??k?omI"&?Ћ6?88?QA&._?ȍ)?I ? 5? ?-|?>\w?͙hRs?ަRkp?neej??6p 6?%o|0?vA(?ey1!? {Z?!zodBO?|A ?9s~?>ahx4>4fL>^7>n}>aو>sKN>@%Z->Ƭ><5 *>U蕴>PKN3CʊXXt2.npyNUMPYF{'descr': 'gǔ?? Dq䛳?4'na IghIm|<՝2.F dпLۿOƳ\tI?7\UCA8 v.(?7K]?VW wU=_ ?›p?ENo?s\Z˿S?;D!Z޿OK;?_8? ɰ}lc?v?>?)7(U5E?d3^VֿOrs?2?9׿F?J*멿sRRM?^w:?^~@{y?(⿡SݿY;P㿊jئ?YgĮ;Fk7jRBҿ>^&?PKN3C xy.npyNUMPYF{'descr': 'I@+2jA?0x?7qߐg?Y\@?ɰE9?X'wD&#N@`:?.26?+_I9KI?'g\B^gW?. '4EUʕh俍߯{_r?Qf?}yg?`dT?sпF1f?m*Ms0@/-A)a?bs^ۿDg?[)G~eS?utPy? b ? CϿ R~J?rSC)׿?j0?H0?w {?ĥ?fIaԿ[8ܒ?ά2?3?W)?NgljSz?BM?3$D޿LW?oS?5M @o/Q,4@\N?`?@!N?F]ޖ^?"\ԿtkH?FxK⿈U6?CEBBͿ;? N8?ۿ?vkyaw>QЌ3?ўgH4?ęx̝g?4;ça8?uP7? $;=?f?^t,?ܫyJ?jh?*?t7N?e6?dWɸ?*_q޿@CaZ?6??̥?(X?m}⿓U̮޿ 3e/@Vh\ѿ#$8Cy7s`f͵#aC6?XG;?[Gs0x?c??î i@[?FTR?JAs?{/}@go'ԿeI?QȈT?\A>A?PKb*Yrb?b ]?0@ =w?mɿ"}?b?. )'rq_?v~0im?h|??̓xoz4n? &ֽ޿bѿVqO?v/]@R(X]@!4+?qVv?D?2b_*[t(3×N!Ji?N}?NVݿzRS?c0d{[bf)|0z:ZݿV7a? ;.u[5?W {T+!5mip Y?үi?_Q)yXZ?p#+D^ C TK T'?>gǔ?? Dq䛳?4'na IghIm|<՝2.F dпLۿOƳ\tI?7\UCA8 v.(?7K]?VW wU=_ ?›p?ENo?s\Z˿S?;D!Z޿OK;?_8? ɰ}lc?v?>?)7(U5E?d3^VֿOrs?2?9׿F?J*멿sRRM?^w:?^~@{y?(⿡SݿY;P㿊jئ?YgĮ;Fk7jRBҿ>^&?PKN3CFPPPx1.npyNUMPYF{'descr': 'I@+2jA?0x?7qߐg?Y\@?ɰE9?X'wD&#N@`:?.26?+_I9KI?'g\B^gW?. '4EUʕh俍߯{_r?Qf?}yg?`dT?sпF1f?m*Ms0@/-A)a?bs^ۿDg?[)G~eS?utPy? b ? CϿ R~J?rSC)׿?j0?H0?w {?ĥ?fIaԿ[8ܒ?ά2?3?W)?NgljSz?BM?3$D޿LW?oS?5M @o/Q,4@\N?`?@!N?F]ޖ^?"\ԿtkH?FxK⿈U6?CEBBͿ;? N8?ۿ?vkyaw>QЌ3?ўgH4?ęx̝g?4;ça8?uP7? $;=?f?^t,?ܫyJ?jh?*?t7N?e6?dWɸ?*_q޿@CaZ?6??̥?(X?m}⿓U̮޿ 3e/@Vh\ѿ#$8Cy7s`f͵#aC6?XG;?[Gs0x?c??î i@[?FTR?JAs?{/}@go'ԿeI?QȈT?\A>A?PKb*Yrb?b ]?0@ =w?mɿ"}?b?. )'rq_?v~0im?h|??̓xoz4n? &ֽ޿bѿVqO?v/]@R(X]@!4+?qVv?D?2b_*[t(3×NPKN3Cٛ`` xstar.npyNUMPYF{'descr': '),?7S*? J?|6H? xxԙ@Yв5?V۲ (?԰-U@A]?.hl?C?p^գw?K]&?PK 3Cqx.npyNUMPYF{'descr': ''4? k.? gM?'"ݿa?'%QʧC :-c࿑~b?OgaܰǀIA+?ux<)wgs?y~Ҏٙ~¿6YePK 3Cxx xstar.npyNUMPYF{'descr': ' ףp=(\RQGz ףp= ffffff(\Q{Gzףp= 333333(\QHzGp= ףQp= ףp(\(zG꿚RQ ףp= (\zGz233333Qῤp= ף࿸Q޿(\(ܿٿףp= ׿|GzԿQѿQοɿxGzĿQGzGzGz?Gz?Q?Gz??Q?Q?xGz?ףp= ??(\(?Q?p= ף?Q?433333?|Gz?(\? ףp= ?TQ??zG?(\(?p= ףp?Q??p= ף?HzG?Q?(\?433333?أp= ?|Gz? Q?(\?ffffff? ףp= ?Gz?RQ?(\??> ףp=?zG?Q?*\(??p= ףp?Gz?Q?\(\?@PK 3Cpb,y.npyPK 3Cqx.npyPK 3Cxx &xstar.npyPKpyGP_PR-1.0/src/000077500000000000000000000000001222232424600134205ustar00rootroot00000000000000pyGP_PR-1.0/src/Core/000077500000000000000000000000001222232424600143105ustar00rootroot00000000000000pyGP_PR-1.0/src/Core/COV/000077500000000000000000000000001222232424600147375ustar00rootroot00000000000000pyGP_PR-1.0/src/Core/COV/__init__.py000077500000000000000000000000311222232424600170450ustar00rootroot00000000000000 __all__ = ["kernels"] pyGP_PR-1.0/src/Core/COV/kernels.py000077500000000000000000001147721222232424600167730ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 30/09/2013 #================================================================================ # Created on 31/08/2009 # # covariance functions/kernels to be used by Gaussian process functions. # There are three different kinds of covariance functions: simple, composite and FITC: # # simple covariance functions: # # covSEiso - isotropic squared exponential covariance function # covSEard - squared exponential covariance function with ard (Automatic Relevance Determination) # covSEisoU - isotropic squared exponential covariance function with unit magnitude # covLIN - linear covariance function # covLINard - linear covariance function with ard # covPoly - polynomial covariance function # covPPiso - piecewise polynomial covariance function with compact support # covConst - constant covariance function # covMatern - Matern covariance function # covPeriodic - periodic covariance function # covRQiso - rational quadratic covariance function # covRQard - rational quadratic covariance function with ard # covNoise - independent covariance function (ie white noise) # # composite covariance functions (see explanation at the bottom): # # covSum - sum of (parameterized) covariance functions # covProd - product of (parametrized) covariance functions # covScale - composition of a covariance function as a scaled version of another covariance function # covMask - composition of a covariance function from another covariance function using a subset of input dimensions # # FITC: # # covFITC - covariance function to be used together with the FITC (pseudo input) approximation # # Naming convention: all covariance functions start with "cov". A trailing # "iso" means isotropic, "ard" means Automatic Relevance Determination. # # The covariance functions are written according to a special convention where # the exact behaviour depends on the number of input and output arguments # passed to the function. If you want to add new covariance functions, you # should follow this convention if you want them to work with the function gp. # There are four different ways of calling # the covariance functions: # # 1) With no input arguments: # # num_prarams = Tools.general.feval(covfunc) # # The covariance function returns a list with the number of hyperparameters it # expects. For ard covariacne function the list includes a string where "D" # is the dimension of the input space. # For example, calling "covSEiso" returns [2]; "covSEard" returns ['D + 1']. # # 2) With two input arguments: # # K = Tools.general.feval(covfunc, hyp, x) # # The function computes and returns the covariance matrix where logtheta are # the log of the hyperparameters and x is an n by D matrix of input points, where # D is the dimension of the input space. The returned covariance matrix is of # size n by n. # # 3) With three input arguments: # # Ks = Tools.general.feval(covfunc, hyp, x, z) # # The function computes test set covariances; Ks is a (n by nn) matrix of cross # covariances between training cases x and test cases z, # where z is a nn by D matrix. # # 4) With three input arguments where one is 'diag': # # kss = Tools.general.feval(covfunc, hyp, z, 'diag') # # The function computes self-variances; kss is a (nn by 1) vector. # # 5) With four input arguments: # # Deriv = Tools.general.feval(covfunc, hyp, x, z, der) or # Deriv = Tools.general.feval(covfunc, hyp, x, None, der) # # The function computes and returns the n by n resp. n by nn matrix of partial # derivatives of either the training set covariance matrix or the train-test covariance # with respect to logtheta(der), i.e. with respect to the log of hyperparameter number der, # where der in [0,1, ..., num_prarams-1]. # # About the specification of simple and composite covariance functions to be # used by the Gaussian process function gp: # # covfunc = [['kernels.covSEard']] # # Composite covariance functions can be specified as list. For example: # # covfunc = [['kernels.covSum'], [['kernels.covSEard'],['kernels.covNoise']]] # covfunc = [['kernels.covProd'],[['kernels.covPeriodic'],['kernels.covSEiso']]] # # # @author: Marion Neumann (last update 30/09/13) # Substantial updates by Daniel Marthaler Fall 2012. # Shan Huang (Sep. 2013) # # This is a python implementation of gpml functionality (Copyright (c) by # Carl Edward Rasmussen and Hannes Nickisch, 2011-02-18). # # Copyright (c) by Marion Neumann and Daniel Marthaler, 20/05/2013 import Tools import numpy as np import math import scipy.spatial.distance as spdist def covSEiso(hyp=None, x=None, z=None, der=None): # Squared Exponential covariance function with isotropic distance measure. # The covariance function is parameterized as: # # k(x^p,x^q) = sf2 * exp(-(x^p - x^q)'*inv(P)*(x^p - x^q)/2) # # where the P matrix is ell^2 times the unit matrix and # sf2 is the signal variance # # The hyperparameters of the function are: # # hyp = [ log(ell) # log(sqrt(sf2)) ] # a column vector. # Each row of x resp. z is a data point. if hyp == None: # report number of parameters return [2] ell = np.exp(hyp[0]) # characteristic length scale sf2 = np.exp(2.*hyp[1]) # signal variance n,D = x.shape if z == 'diag': A = np.zeros((n,1)) elif z == None: A = spdist.cdist(x/ell, x/ell, 'sqeuclidean') else: # compute covariance between data sets x and z A = spdist.cdist(x/ell, z/ell, 'sqeuclidean') # self covariances if der == None: # compute covariance matix for dataset x A = sf2 * np.exp(-0.5*A) else: if der == 0: # compute derivative matrix wrt 1st parameter A = sf2 * np.exp(-0.5*A) * A elif der == 1: # compute derivative matrix wrt 2nd parameter A = 2. * sf2 * np.exp(-0.5*A) else: raise Exception("Calling for a derivative in covSEiso that does not exist") return A def covSEard(hyp=None, x=None, z=None, der=None): # Squared Exponential covariance function with Automatic Relevance Detemination # (ARD) distance measure. The covariance function is parameterized as: # # k(x^p,x^q) = sf2 * exp(-(x^p - x^q)'*inv(P)*(x^p - x^q)/2) # # where the P matrix is diagonal with ARD parameters ell_1^2,...,ell_D^2, where # D is the dimension of the input space and sf2 is the signal variance. # # The hyperparameters are: # # hyp = [ log(ell_1) # log(ell_2) # . # log(ell_D) # log(sqrt(sf2)) ] if hyp == None: # report number of parameters return ['D + 1'] # USAGE: integer OR D_+_int (spaces are SIGNIFICANT) [n, D] = x.shape ell = 1/np.exp(hyp[0:D]) # characteristic length scale sf2 = np.exp(2.*hyp[D]) # signal variance if z == 'diag': A = np.zeros((n,1)) elif z == None: tmp = np.dot(np.diag(ell),x.T).T A = spdist.cdist(tmp, tmp, 'sqeuclidean') else: # compute covariance between data sets x and z A = spdist.cdist(np.dot(np.diag(ell),x.T).T, np.dot(np.diag(ell),z.T).T, 'sqeuclidean') # cross covariances A = sf2*np.exp(-0.5*A) if not der == None: if der < D: # compute derivative matrix wrt length scale parameters if z == 'diag': A = A*0. elif z == None: tmp = np.atleast_2d(x[:,der])/ell[der] A = A * spdist.cdist(tmp, tmp, 'sqeuclidean') else: A = A * spdist.cdist(np.atleast_2d(x[:,der]).T/ell[der], np.atleast_2d(z[:,der]).T/ell[der], 'sqeuclidean') elif der==D: # compute derivative matrix wrt magnitude parameter A = 2.*A else: raise Exception("Wrong derivative index in covSEard") return A def covSEisoU(hyp=None, x=None, z=None, der=None): # Squared Exponential covariance function with isotropic distance measure with # unit magnitude. The covariance function is parameterized as: # # k(x^p,x^q) = exp( -(x^p - x^q)' * inv(P) * (x^p - x^q) / 2 ) # # where the P matrix is ell^2 times the unit matrix. # # The hyperparameters of the function are: # # hyp = [ log(ell) ] if hyp == None: # report number of parameters return [1] ell = np.exp(hyp[0]) # characteristic length scale n,D = x.shape if z == 'diag': A = np.zeros((n,1)) elif z == None: A = spdist.cdist(x/ell, x/ell, 'sqeuclidean') else: # compute covariance between data sets x and z A = spdist.cdist(x/ell, z/ell, 'sqeuclidean') # self covariances if der == None: # compute covariance matix for dataset x A = np.exp(-0.5*A) else: if der == 0: # compute derivative matrix wrt 1st parameter A = np.exp(-0.5*A) * A else: raise Exception("Wrong derivative index in covSEisoU") return A def covLIN(hyp=None, x=None, z=None, der=None): # Linear Covariance function. # The covariance function is parameterized as: # k(x^p,x^q) = x^p'*x^q # # There are no hyperparameters: # # hyp = [] # # Note that there is no bias or scale term; use covConst and covScale to add these. if hyp == None: # report number of parameters return [0] n,m = x.shape if z == 'diag': A = np.reshape(np.sum(x*x,1), (n,1)) elif z == None: A = np.dot(x,x.T) + np.eye(n)*1e-16 # required for numerical accuracy else: # compute covariance between data sets x and z A = np.dot(x,z.T) # cross covariances if der: raise Exception("No derivative available in covLIN") return A def covLINard(hyp=None, x=None, z=None, der=None): # Linear covariance function with Automatic Relevance Detemination # (ARD) distance measure. The covariance function is parameterized as: # k(x^p,x^q) = x^p' * inv(P) * x^q # # where the P matrix is diagonal with ARD parameters ell_1^2,...,ell_D^2, where # D is the dimension of the input space and sf2 is the signal variance. The # hyperparameters are: # # hyp = [ log(ell_1) # log(ell_2) # . # . # log(ell_D) ] # # Note that there is no bias term; use covConst to add a bias. if hyp == None: # report number of parameters return ['D + 0'] # USAGE: integer OR D_+_int (spaces are SIGNIFICANT) n, D = x.shape ell = np.exp(hyp) # characteristic length scales x = np.dot(x,np.diag(1./ell)) if z == 'diag': A = np.reshape(np.sum(x*x,1), (n,1)) elif z == None: A = np.dot(x,x.T) else: # compute covariance between data sets x and z z = np.dot(z,np.diag(1./ell)) A = np.dot(x,z.T) # cross covariances if not der == None and der < D: if z == 'diag': A = -2.*x[:,der]*x[:,der] elif z == None: A = -2.*np.dot(x[:,der],x[:,der].T) else: A = -2.*np.dot(x[:,der],z[:,der].T) # cross covariances elif der: raise Exception("Wrong derivative index in covLINard") return A def covPoly(hyp=None, x=None, z=None,der=None): # Polynomial covariance function # The covariance function is parameterized as: # k(x^p,x^q) = sf2 * ( c + (x^p)'*(x^q) ) ** d # # The hyperparameters of the function are: # # hyp = [ log(c) # log(sqrt(sf2)) # d ] # # NOTE: d is not treated as a hyperparameter. if hyp == None: # report number of parameters return [3] c = np.exp(hyp[0]) # inhomogeneous offset sf2 = np.exp(2.*hyp[1]) # signal variance d = hyp[2] # degree of polynomical if np.abs(d-np.round(d)) < 1e-8: # remove numerical error from format of parameter d = int(round(d)) assert(d >= 1.) # only nonzero integers for d d = int(d) n, D = x.shape if z == 'diag': A = np.reshape(np.sum(x*x,1), (n,1)) elif z==None: A = np.dot(x,x.T) else: # compute covariance between data sets x and z A = np.dot(x,z.T) # cross covariances if der == None: # compute covariance matix for dataset x A = sf2 * (c + A)**d else: if der == 0: # compute derivative matrix wrt 1st parameter A = c * d * sf2 * (c+A)**(d-1) elif der == 1: # compute derivative matrix wrt 2nd parameter A = 2. * sf2 * (c + A)**d elif der == 2: # NOTE: d is not treated as a hyperparameter -> we do not want to optimize w.r.t. d A = np.zeros_like(A) else: raise Exception("Wrong derivative entry in covPoly") return A def covPPiso(hyp=None, x=None, z=None, der=None): # Piecewise polynomial covariance function with compact support # The covariance function is: # # k(x^p,x^q) = s2f * (1-r)_+.^j * f(r,j) # # where r is the distance sqrt((x^p-x^q)'*inv(P)*(x^p-x^q)), P is ell^2 times # the unit matrix and sf2 is the signal variance. # The hyperparameters are: # # hyp = [ log(ell) # log(sqrt(sf2)) # log(v) ] def ppmax(A,B): return np.maximum(A,B*np.ones_like(A)) def func(v,r,j): if v == 0: return 1 elif v == 1: return ( 1. + (j+1) * r ) elif v == 2: return ( 1. + (j+2)*r + (j*j + 4.*j+ 3)/3.*r*r ) elif v == 3: return ( 1. + (j+3)*r + (6.*j*j+36.*j+45.)/15.*r*r + (j*j*j+9.*j*j+23.*j+15.)/15.*r*r*r ) else: raise Exception (["Wrong degree in covPPiso. Should be 0,1,2 or 3, is " + str(v)]) def dfunc(v,r,j): if v == 0: return 0 elif v == 1: return ( j+1 ) elif v == 2: return ( (j+2) + 2.*(j*j+ 4.*j+ 3.)/3.*r ) elif v == 3: return ( (j+3) + 2.*(6.*j*j+36.*j+45.)/15.*r + (j*j*j+9.*j*j+23.*j+15.)/5.*r*r ) else: raise Exception (["Wrong degree in covPPiso. Should be 0,1,2 or 3, is " + str(v)]) def pp(r,j,v,func): return func(v,r,j)*(ppmax(1-r,0)**(j+v)) def dpp(r,j,v,func,dfunc): return ppmax(1-r,0)**(j+v-1) * r * ( (j+v)*func(v,r,j) - ppmax(1-r,0) * dfunc(v,r,j) ) if hyp == None: # report number of parameters return [3] ell = np.exp(hyp[0]) # characteristic length scale sf2 = np.exp(2.*hyp[1]) # signal variance v = np.exp(hyp[2]) # degree (v = 0,1,2 or 3 only) if np.abs(v-np.round(v)) < 1e-8: # remove numerical error from format of parameter v = int(round(v)) assert(int(v) in range(4)) # Only allowed degrees: 0,1,2 or 3 v = int(v) n, D = x.shape j = math.floor(0.5*D) + v + 1 if z == 'diag': A = np.zeros((n,1)) elif z == None: A = np.sqrt( spdist.cdist(x/ell, x/ell, 'sqeuclidean') ) else: # compute covariance between data sets x and z A = np.sqrt( spdist.cdist(x/ell, z/ell, 'sqeuclidean') ) # cross covariances if der == None: # compute covariance matix for dataset x A = sf2 * pp(A,j,v,func) else: if der == 0: # compute derivative matrix wrt 1st parameter A = sf2 * dpp(A,j,v,func,dfunc) elif der == 1: # compute derivative matrix wrt 2nd parameter A = 2. * sf2 * pp(A,j,v,func) elif der == 2: # wants to compute derivative wrt order A = np.zeros_like(A) else: raise Exception("Wrong derivative entry in covPPiso") return A def covConst(hyp=None, x=None, z=None, der=None): # Covariance function for a constant function. # The covariance function is parameterized as: # k(x^p,x^q) = sf2 # # The scalar hyperparameter is: # # hyp = [ log(sqrt(sf2)) ] if hyp == None: # report number of parameters return [1] sf2 = np.exp(2.*hyp[0]) # s2 n,m = x.shape if z == 'diag': A = sf2*np.ones((n,1)) elif z == None: A = sf2 * np.ones((n,n)) else: A = sf2*np.ones((n,z.shape[0])) if der == 0: # compute derivative matrix wrt sf2 A = 2. * A elif der: raise Exception("Wrong derivative entry in covConst") return A def covMatern(hyp=None, x=None, z=None, der=None): # Matern covariance function with nu = d/2 and isotropic distance measure. For d=1 # the function is also known as the exponential covariance function or the # Ornstein-Uhlenbeck covariance in 1d. The covariance function is: # # k(x^p,x^q) = s2f * f( sqrt(d)*r ) * exp(-sqrt(d)*r) # # with f(t)=1 for d=1, f(t)=1+t for d=3 and f(t)=1+t+(t*t)/3 for d=5. # Here, r is the distance sqrt( (x^p-x^q)'*inv(P)*(x^p-x^q)), # where P is ell times the unit matrix and sf2 is the signal variance. # # The hyperparameters of the function are: # # hyp = [ log(ell) # log(sqrt(sf2)) # log(d) ] def func(d,t): if d == 1: return 1 elif d == 3: return 1 + t elif d == 5: return 1 + t*(1+t/3.) else: raise Exception("Wrong value for d in covMatern") def dfunc(d,t): if d == 1: return 1 elif d == 3: return t elif d == 5: return t*(1+t/3.) else: raise Exception("Wrong value for d in covMatern") def mfunc(d,t): return func(d,t)*np.exp(-1.*t) def dmfunc(d,t): return dfunc(d,t)*t*np.exp(-1.*t) if hyp == None: # report number of parameters return [3] ell = np.exp(hyp[0]) # characteristic length scale sf2 = np.exp(2.*hyp[1]) # signal variance d = np.exp(hyp[2]) # 2 times nu if np.abs(d-np.round(d)) < 1e-8: # remove numerical error from format of parameter d = int(round(d)) try: assert(int(d) in [1,3,5]) # Check for valid values of d except AssertionError: d = 3 d = int(d) if z == 'diag': A = np.zeros((x.shape[0],1)) elif z == None: x = np.sqrt(d)*x/ell A = np.sqrt(spdist.cdist(x, x, 'sqeuclidean')) else: x = np.sqrt(d)*x/ell z = np.sqrt(d)*z/ell A = np.sqrt(spdist.cdist(x, z, 'sqeuclidean')) if der == None: # compute covariance matix for dataset x A = sf2 * mfunc(d,A) else: if der == 0: # compute derivative matrix wrt 1st parameter A = sf2 * dmfunc(d,A) elif der == 1: # compute derivative matrix wrt 2nd parameter A = 2 * sf2 * mfunc(d,A) elif der == 2: # Wants to compute derivative wrt nu A = np.zeros_like(A) # Do nothing else: raise Exception("Wrong derivative value in covMatern") return A def covPeriodic(hyp=None, x=None, z=None, der=None): # Stationary covariance function for a smooth periodic function, # with period p: # # k(x^p,x^q) = sf2 * exp( -2*sin^2( pi*||x^p - x^q)||/p )/ell**2 ) # # The hyperparameters of the function are: # hyp = [ log(ell) # log(p) # log(sqrt(sf2)) ] if hyp == None: # report number of parameters return [3] ell = np.exp(hyp[0]) # characteristic length scale p = np.exp(hyp[1]) # period sf2 = np.exp(2.*hyp[2]) # signal variance n,D = x.shape if z == 'diag': A = np.zeros((n,1)) elif z == None: A = np.sqrt(spdist.cdist(x, x, 'sqeuclidean')) else: A = np.sqrt(spdist.cdist(x, z, 'sqeuclidean')) A = np.pi*A/p if der == None: # compute covariance matix for dataset x A = np.sin(A)/ell A = A * A A = sf2 *np.exp(-2.*A) else: if der == 0: # compute derivative matrix wrt 1st parameter A = np.sin(A)/ell A = A * A A = 4. *sf2 *np.exp(-2.*A) * A elif der == 1: # compute derivative matrix wrt 2nd parameter R = np.sin(A)/ell A = 4 * sf2/ell * np.exp(-2.*R*R)*R*np.cos(A)*A elif der == 2: # compute derivative matrix wrt 3rd parameter A = np.sin(A)/ell A = A * A A = 2. * sf2 * np.exp(-2.*A) else: raise Exception("Wrong derivative index in covPeriodic") return A def covRQiso(hyp=None, x=None, z=None, der=None): # Rational Quadratic covariance function with isotropic distance measure. # The covariance function is parameterized as: # # k(x^p,x^q) = sf2 * [1 + (x^p - x^q)'*inv(P)*(x^p - x^q)/(2*alpha)]^(-alpha) # # where the P matrix is ell^2 times the unit matrix, # sf2 is the signal variance, and alpha is the shape parameter for the RQ # covariance. # # The hyperparameters of the function are: # hyp = [ log(ell) # log(sqrt(sf2)) # log(alpha) ] # # each row of x/z is a data point if hyp == None: # report number of parameters return [3] ell = np.exp(hyp[0]) # characteristic length scale sf2 = np.exp(2.*hyp[1]) # signal variance alpha = np.exp(hyp[2]) n,D = x.shape if z == 'diag': D2 = np.zeros((n,1)) elif z == None: D2 = spdist.cdist(x/ell, x/ell, 'sqeuclidean') else: D2 = spdist.cdist(x/ell, z/ell, 'sqeuclidean') if der == None: # compute covariance matix for dataset x A = sf2 * ( ( 1.0 + 0.5*D2/alpha )**(-alpha) ) else: if der == 0: # compute derivative matrix wrt 1st parameter A = sf2 * ( 1.0 + 0.5*D2/alpha )**(-alpha-1) * D2 elif der == 1: # compute derivative matrix wrt 2nd parameter A = 2.* sf2 * ( ( 1.0 + 0.5*D2/alpha )**(-alpha) ) elif der == 2: # compute derivative matrix wrt 3rd parameter K = ( 1.0 + 0.5*D2/alpha ) A = sf2 * K**(-alpha) * (0.5*D2/K - alpha*np.log(K) ) else: raise Exception("Wrong derivative index in covRQiso") return A def covRQard(hyp=None, x=None, z=None, der=None): # Rational Quadratic covariance function with Automatic Relevance Detemination # (ARD) distance measure. The covariance function is parameterized as: # # k(x^p,x^q) = sf2 * [1 + (x^p - x^q)'*inv(P)*(x^p - x^q)/(2*alpha)]^(-alpha) # # where the P matrix is diagonal with ARD parameters ell_1^2,...,ell_D^2, where # D is the dimension of the input space, sf2 is the signal variance and alpha is # the shape parameter for the RQ covariance. The hyperparameters are: # # hyp = [ log(ell_1) # log(ell_2) # ... # log(ell_D) # log(sqrt(sf2)) # log(alpha)] if hyp == None: # report number of parameters return ['D + 2'] # USAGE: integer OR D_+_int (spaces are SIGNIFICANT) [n, D] = x.shape ell = 1/np.exp(hyp[0:D]) # characteristic length scale sf2 = np.exp(2.*hyp[D]) # signal variance alpha = np.exp(hyp[D+1]) if z == 'diag': D2 = np.zeros((n,1)) elif z == None: tmp = np.dot(np.diag(ell),x.T).T D2 = spdist.cdist(tmp, tmp, 'sqeuclidean') else: D2 = spdist.cdist(np.dot(np.diag(ell),x.T).T, np.dot(np.diag(ell),z.T).T, 'sqeuclidean') if der == None: # compute covariance matix for dataset x A = sf2 * ( ( 1.0 + 0.5*D2/alpha )**(-alpha) ) else: if der < D: # compute derivative matrix wrt length scale parameters if z == 'diag': A = D2*0 elif z == None: tmp = np.atleast_2d(x[:,der])/ell[der] A = sf2 * ( 1.0 + 0.5*D2/alpha )**(-alpha-1) * spdist.cdist(tmp, tmp, 'sqeuclidean') else: A = sf2 * ( 1.0 + 0.5*D2/alpha )**(-alpha-1) * spdist.cdist(np.atleast_2d(x[:,der]).T/ell[der], np.atleast_2d(z[:,der]).T/ell[der], 'sqeuclidean') elif der==D: # compute derivative matrix wrt magnitude parameter A = 2. * sf2 * ( ( 1.0 + 0.5*D2/alpha )**(-alpha) ) elif der==(D+1): # compute derivative matrix wrt magnitude parameter K = ( 1.0 + 0.5*D2/alpha ) A = sf2 * K**(-alpha) * ( 0.5*D2/K - alpha*np.log(K) ) else: raise Exception("Wrong derivative index in covRQard") return A def covNoise(hyp=None, x=None, z=None, der=None): # Independent covariance function, ie "white noise", with specified variance. # The covariance function is specified as: # # k(x^p,x^q) = s2 * \delta(p,q) # # where s2 is the noise variance and \delta(p,q) is a Kronecker delta function # which is 1 iff p=q and zero otherwise. The hyperparameter is # # hyp = [ log(sqrt(s2)) ] tol = 1.e-9 # Tolerance for declaring two vectors "equal" if hyp == None: # report number of parameters return [1] s2 = np.exp(2.*hyp[0]) # noise variance n,D = x.shape if z == 'diag': A = np.ones((n,1)) elif z == None: A = np.eye(n) else: # compute covariance between data sets x and z M = spdist.cdist(x, z, 'sqeuclidean') A = np.zeros_like(M,dtype=np.float) A[M < tol] = 1. if der == None: A = s2*A else: # compute derivative matrix if der == 0: A = 2.*s2*A else: raise Exception("Wrong derivative index in covNoise") return A def covSum(covfunc, hyp=None, x=None, z=None, der=None): # covSum - compose a covariance function as the sum of other covariance # functions. This function doesn't actually compute very much on its own, it # merely does some bookkeeping, and calls other covariance functions to do the # actual work. def DetermineNumberOfParameters(v,no_param): if isinstance(no_param, int): v.append(no_param) elif isinstance(no_param,str): # no_param is a string pram_str = no_param.split(' ') if pram_str[0]=='D': temp = int(D) if pram_str[1]=='+': temp += int(pram_str[2]) elif pram_str[1]=='-': temp -= int(pram_str[2]) else: raise Exception(["Error: number of parameters of "+covfunc[i] +" unknown!"]) v.append(temp) elif isinstance(no_param, list): # The number of hyperparameters for this piece of covfunc is the sum # of all of them in this composition temp = [0] for jj in xrange(len(no_param)): DetermineNumberOfParameters(temp,no_param[jj]) v.append(sum(temp)) else: # This is an error, we should never be here raise Exception("Error in return of number of parameters") return v if hyp == None: # report number of parameters A = [Tools.general.feval(covfunc[0])] for i in range(1,len(covfunc)): A.append(Tools.general.feval(covfunc[i])) return A [n, D] = x.shape # SET vector v (v indicates how many parameters each covfunc has # (NOTE : v[i]=number of parameters + 1 -> this is because of the indexing of python!)) v = [0] # needed for technical reasons for ii in range(1,len(covfunc)+1): no_param = Tools.general.feval(covfunc[ii-1]) DetermineNumberOfParameters(v,no_param) AT = Tools.general.feval(covfunc[0], hyp[:v[1]], x, z) A = np.zeros_like(AT) # Allocate covariance Matrix if der == None: # compute covariance matrix for ii in range(1,len(covfunc)+1): # iteration over summand functions f = covfunc[ii-1] s = sum(v[0:ii]) A = A + Tools.general.feval(f, hyp[s:(s+v[ii])], x, z) # accumulate covariances else: # compute derivative matrices tmp = 0 for ii in range(1,len(covfunc)+1): tmp += v[ii] if der this is because of the indexing of python!)) v = [0] for ii in range(1,len(covfunc)+1): no_param = Tools.general.feval(covfunc[ii-1]) DetermineNumberOfParameters(v,no_param) AT = Tools.general.feval(covfunc[0], hyp[:v[1]], x, z) A = np.ones_like(AT) # Allocate covariance Matrix if der == None: # compute covariance matrix for ii in range(1,len(covfunc)+1): # iteration over multiplicand functions f = covfunc[ii-1] s = sum(v[0:ii]) A *= Tools.general.feval(f, hyp[s:(s+v[ii])], x, z) # accumulate covariances else: tmp = 0 flag = True for ii in range(1,len(covfunc)+1): tmp += v[ii] if der= 0) if der == None: # compute covariance matix for dataset x if z == None: A = Tools.general.feval(cov, hyp, x[:,mask]) else: # compute covariance between data sets x and z if z == 'diag': A = Tools.general.feval(cov,hyp,x[:,mask],z) else: # coumpute cross covariances A = Tools.general.feval(cov,hyp,x[:,mask],z[:,mask]) else: # compute derivatives if z == None: A = Tools.general.feval(cov, hyp, x[:,mask],z,der) elif z == 'diag': A = Tools.general.feval(cov,hyp,x[:,mask],z,der) else: A = Tools.general.feval(cov,hyp,x[:,mask],z[:,mask],der) return A def covFITC(covfunc, xu=None, hyp=None, x=None, z=None, der=None): # Covariance function to be used together with the FITC approximation. # # The function allows for more than one output argument and does not respect the # interface of a proper covariance function. In fact, it wraps a proper # covariance function such that it can be used together with infFITC.m. # Instead of outputing the full covariance, it returns cross-covariances between # the inputs x, z and the inducing inputs xu as needed by infFITC.m # # Copyright (c) by Ed Snelson, Carl Edward Rasmussen and Hannes Nickisch, 2010-12-21. # # NOTE: The first element of cov should be ['kernels.covFITC'] if hyp == None: # report number of parameters A = [Tools.general.feval(covfunc)] return A try: assert(xu.shape[1]==x.shape[1]) except AssertionError: raise Exception('Dimensionality of inducing inputs must match training inputs') n,D = x.shape if der == None: # compute covariance matrices for dataset x if z == None: K = Tools.general.feval(covfunc,hyp,x,'diag') Kuu = Tools.general.feval(covfunc,hyp,xu) Ku = Tools.general.feval(covfunc,hyp,xu,x) elif z == 'diag': K = Tools.general.feval(covfunc,hyp,x,z) return K else: K = Tools.general.feval(covfunc,hyp,xu,z) return K else: # compute derivative matrices if z == None: K = Tools.general.feval(covfunc,hyp,x,'diag',der) Kuu = Tools.general.feval(covfunc,hyp,xu,None,der) Ku = Tools.general.feval(covfunc,hyp,xu,x,der) elif z == 'diag': K = Tools.general.feval(covfunc,hyp,x,z,der) return K else: K = Tools.general.feval(covfunc,hyp,xu,z,der) return K return K, Kuu, Ku pyGP_PR-1.0/src/Core/INF/000077500000000000000000000000001222232424600147245ustar00rootroot00000000000000pyGP_PR-1.0/src/Core/INF/__init__.py000077500000000000000000000000251222232424600170350ustar00rootroot00000000000000 __all__ = ["inf"] pyGP_PR-1.0/src/Core/INF/inf.py000066400000000000000000001423331222232424600160600ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 20/05/2013 #================================================================================ # Inference methods: Compute the (approximate) posterior for a Gaussian process. # Methods currently implemented include: # # infExact Exact inference (only possible with Gaussian likelihood) # infLaplace Laplace's Approximation # infEP Expectation Propagation # infVB [NOT IMPLEMENTED!] Variational Bayes Approximation # # infFITC Large scale regression with approximate covariance matrix # infFITC_Laplace Large scale inference with approximate covariance matrix # infFITC_EP Large scale inference with approximate covariance matrix # # infMCMC [NOT IMPLEMENTED!] # Markov Chain Monte Carlo and Annealed Importance Sampling # We offer two samplers. # - hmc: Hybrid Monte Carlo # - ess: Elliptical Slice Sampling # No derivatives w.r.t. to hyperparameters are provided. # # infLOO [NOT IMPLEMENTED!] # Leave-One-Out predictive probability and Least-Squares Approxim. # # The interface to the approximation methods is the following: # # function [post nlZ dnlZ] = inf..(hyp, cov, lik, x, y) # # where: # INPUT: # hyp struct of hyperparameters # cov name of the covariance function (see covFunctions.m) # lik name of the likelihood function (see likFunctions.m) # x n by D matrix of training inputs # y 1d array (of size n) of targets # # OUTPUT: # post struct representation of the (approximate) posterior containing: # alpha 1d array containing inv(K)*m, # where K is the prior covariance matrix and m the approx posterior mean # sW 1d array containing diagonal of sqrt(W) # the approximate posterior covariance matrix is inv(inv(K)+W) # L 2d array, L = chol(sW*K*sW+identity(n)) # nlZ returned value of the negative log marginal likelihood # dnlZ 1d array of partial derivatives of the negative log marginal likelihood # w.r.t. each hyperparameter # # Usually, the approximate posterior to be returned admits the form # N(m=K*alpha, V=inv(inv(K)+W)), where alpha is a vector and W is diagonal; # if not, then L contains instead -inv(K+inv(W)), and sW is unused. # # For more information on the individual approximation methods and their # implementations, see the respective inf* function below. See also gp.py # # @author: Daniel Marthaler (Fall 2012) # Substantial updates by Shan Huang (Sep. 2013) # # This is a python implementation of gpml functionality (Copyright (c) by # Carl Edward Rasmussen and Hannes Nickisch, 2013-01-21). # # Copyright (c) by Marion Neumann and Daniel Marthaler, 20/05/2013 import numpy as np import Tools.general from ..UTIL.solve_chol import solve_chol from copy import copy, deepcopy from ..UTIL.utils import randperm, cholupdate, dnlzStruct from ..UTIL.brentmin import brentmin np.seterr(all='ignore') class postStruct: def __init__(self): self.alpha = np.array([]) self.L = np.array([]) self.sW = np.array([]) def infLaplace(hyp, meanfunc, covfunc, likfunc, x, y,nargout=1): # function [post nlZ dnlZ] = infLaplace(hyp, mean, cov, lik, x, y) # Laplace approximation to the posterior Gaussian process. # The function takes a specified covariance function (see kernels.py) and # likelihood function (see lik.py). tol = 1e-6; # tolerance for when to stop the Newton iterations smax = 2; Nline = 20; thr = 1e-4; # line search parameters maxit = 20; # max number of Newton steps in f inffunc = 'inf.infLaplace' K = Tools.general.feval(covfunc, hyp.cov, x) # evaluate the covariance$ m = Tools.general.feval(meanfunc, hyp.mean, x) # evaluate the mean vector n, D = x.shape Psi_old = np.inf # make sure while loop starts by the largest old objective val if "last_alpha" not in infLaplace.__dict__: # find a good starting point for alpha and f alpha = np.zeros((n,1)); f = np.dot(K,alpha) + m # start at mean if sizes not match vargout = Tools.general.feval(likfunc, hyp.lik, y, f, None, inffunc, None, 3) lp = vargout[0]; dlp = vargout[1]; d2lp = vargout[2] W= -d2lp; Psi_new = -lp.sum() else: alpha = last_alpha; f = np.dot(K,alpha) + m # try last one vargout = Tools.general.feval(likfunc, hyp.lik, y, f, None, inffunc, None, 3) lp = vargout[0]; dlp = vargout[1]; d2lp = vargout[2] W=-d2lp; Psi_new = np.dot(alpha.T,(f-m))/2. - lp.sum() # objective for last alpha vargout = -Tools.general.feval(likfunc, hyp.lik, y, m, None, inffunc, None, 1) Psi_def = vargout[0] # objective for default init f==m if Psi_def < Psi_new: # if default is better, we use it alpha = np.zeros((n,1)); f = np.dot(K,alpha) + m vargout = Tools.general.feval(likfunc, hyp.lik, y, f, None, inffunc, None, 3) lp = vargout[0]; dlp = vargout[1]; d2lp = vargout[2] W=-d2lp; Psi_new = -lp.sum() isWneg = np.any(W<0) # flag indicating whether we found negative values of W it = 0 # this happens for the Student's t likelihood while (Psi_old - Psi_new > tol) and it2: # do we want derivatives? dnlZ = dnlzStruct(hyp) # allocate space for derivatives if isWneg: # switch between Cholesky and LU decomposition mode Z = -post.L # inv(K+inv(W)) g = np.atleast_2d((iA*K).sum(axis=1)).T /2 # deriv. of ln|B| wrt W; g = diag(inv(inv(K)+diag(W)))/2 else: Z = np.tile(sW,(1,n))*solve_chol(post.L,np.diag(np.reshape(sW,(sW.shape[0],)))) #sW*inv(B)*sW=inv(K+inv(W)) C = np.linalg.solve(post.L.T,np.tile(sW,(1,n))*K) # deriv. of ln|B| wrt W g = np.atleast_2d((np.diag(K)-(C**2).sum(axis=0).T)).T /2. # g = diag(inv(inv(K)+W))/2 dfhat = g* d3lp # deriv. of nlZ wrt. fhat for ii in range(len(hyp.cov)): # covariance hypers dK = Tools.general.feval(covfunc, hyp.cov, x, None, ii) dnlZ.cov[ii] = (Z*dK).sum()/2. - np.dot(alpha.T,np.dot(dK,alpha))/2. # explicit part b = np.dot(dK,dlp) tmp = np.dot(dfhat.T,b-np.dot(K,np.dot(Z,b))) dnlZ.cov[ii] -= np.dot(dfhat.T,b-np.dot(K,np.dot(Z,b)))[0,0] # implicit part for ii in range(len(hyp.lik)): # likelihood hypers [lp_dhyp,dlp_dhyp,d2lp_dhyp] = Tools.general.feval(likfunc,hyp.lik,y,f,None,inffunc,ii,3) dnlZ.lik[ii] = -np.dot(g.T,d2lp_dhyp) - lp_dhyp.sum() # explicit part b = np.dot(K,dlp_dhyp) dnlZ.lik[ii] -= np.dot(dfhat.T,b-np.dot(K,np.dot(Z,b)))[0,0] # implicit part for ii in range(len(hyp.mean)): # mean hypers dm = Tools.general.feval(meanfunc, hyp.mean, x, ii) dnlZ.mean[ii] = -np.dot(alpha.T,dm) # explicit part dnlZ.mean[ii] -= np.dot(dfhat.T,dm-np.dot(K,np.dot(Z,dm)))[0,0] # implicit part vargout = [post,nlZ,dnlZ] else: vargout = [post, nlZ] return vargout # criterion Psi at alpha + s*dalpha for line search def Psi_line(s,dalpha,alpha,hyp,K,m,likfunc,y,inffunc): alpha = alpha + s*dalpha f = np.dot(K,alpha) + m [lp,dlp,d2lp] = Tools.general.feval(likfunc,hyp.lik,y,f,None,inffunc,None,3) W = -d2lp Psi = np.dot(alpha.T,(f-m))/2. - lp.sum() return Psi[0],alpha,f,dlp,W # Compute the log determinant ldA and the inverse iA of a square nxn matrix # A = eye(n) + K*diag(w) from its LU decomposition; for negative definite A, we # return ldA = Inf. We also return mwiA = -diag(w)*inv(A). def logdetA(K,w,nargout): n = K.shape[0] assert(K.shape[0] == K.shape[1]) A = np.eye(n) + K*np.tile(w.T,(n,1)) [L,U,P] = np.linalg.lu(A) # compute LU decomposition, A = P'*L*U u = np.diag(U) signU = np.prod(np.sign(u)) # sign of U detP = 1 # compute sign (and det) of the permutation matrix P p = np.dot(P,np.array(range(n)).T) for ii in range(n): if ii != p[ii]: detP = -detP j= [jj for jj,val in enumerate(p) if val == ii] p[ii,j[0]] = p[j[0],ii] if signU != detP: # log becomes complex for negative values, encoded by infinity ldA = np.inf else: # det(L) = 1 and U triangular => det(A) = det(P)*prod(diag(U)) ldA = np.log(np.abs(u)).sum() if nargout>1: iA = inv(U)*inv(L)*P; if nargout>2: mwiA = -np.tile(w,(1,n))*iA return ldA,iA,mwiA else: return ldA,iA else: return ldA def infFITC_Laplace(hyp, meanfunc, covfunc, likfunc, x, y,nargout=1): # function [post nlZ dnlZ] = infFITC_Laplace(hyp, mean, cov, lik, x, y) # # FITC-Laplace approximation to the posterior Gaussian process. The function is # equivalent to infLaplace with the covariance function: # # Kt = Q + G; G = diag(g); g = diag(K-Q); Q = Ku'*inv(Kuu + snu2*eye(nu))*Ku; # # where Ku and Kuu are covariances w.r.t. to inducing inputs xu and # snu2 = sn2/1e6 is the noise of the inducing inputs. We fixed the standard # deviation of the inducing inputs snu to be a one per mil of the measurement # noise's standard deviation sn. In case of a likelihood without noise # parameter sn2, we simply use snu2 = 1e-6. # # The implementation exploits the Woodbury matrix identity # inv(Kt) = inv(G) - inv(G)*Ku'*inv(Kuu+Ku*inv(G)*Ku')*Ku*inv(G) # in order to be applicable to large datasets. The computational complexity # is O(n nu^2) where n is the number of data points x and nu the number of # inducing inputs in xu. # The posterior N(f|h,Sigma) is given by h = m+mu with mu = nn + P'*gg and # Sigma = inv(inv(K)+diag(W)) = diag(d) + P'*R0'*R'*R*R0*P. # # The function takes a specified covariance function (see kernels.py) and # likelihood function (lik.py), and is designed to be used with # gp.py and in conjunction with covFITC. cov1 = covfunc[0] if not cov1 == ['kernels.covFITC']: raise Exception('Only covFITC supported.') # check cov tol = 1e-6; # tolerance for when to stop the Newton iterations smax = 2; Nline = 100; thr = 1e-4; # line search parameters maxit = 20; # max number of Newton steps in f inffunc = 'inf.infLaplace' diagK,Kuu,Ku = Tools.general.feval(covfunc, hyp.cov, x) # evaluate the covariance$ m = Tools.general.feval(meanfunc, hyp.mean, x) # evaluate the mean vector if hyp.lik: # hard coded inducing inputs noise sn2 = np.exp(2.*hyp.lik[-1]) snu2 = 1.e-6*sn2 # similar to infFITC else: snu2 = 1.e-6 n, D = x.shape; nu = Kuu.shape[0] rot180 = lambda A: np.rot90(np.rot90(A)) # little helper funct chol_inv = lambda A: np.linalg.solve( rot180( np.linalg.cholesky(rot180(A)) ),np.eye(nu)) # chol(inv(A)) R0 = chol_inv(Kuu+snu2*np.eye(nu)) # initial R, used for refresh O(nu^3) V = np.dot(R0,Ku); d0 = diagK - np.array([(V*V).sum(axis=0)]).T # initial d, needed Psi_old = np.inf # make sure while loop starts by the largest old objective val if "last_alpha" not in infFITC_Laplace.__dict__: # find a good starting point for alpha and f alpha = np.zeros((n,1)); f = mvmK(alpha,V,d0) + m # start at mean if sizes not match vargout = Tools.general.feval(likfunc, hyp.lik, y, f, None, inffunc, None, 3) lp = vargout[0]; dlp = vargout[1]; d2lp = vargout[2] W=-d2lp; Psi_new = -lp.sum() else: alpha = last_alpha; f = mvmK(alpha,V,d0) + m # try last one vargout = Tools.general.feval(likfunc, hyp.lik, y, f, None, inffunc, None, 3) lp = vargout[0]; dlp = vargout[1]; d2lp = vargout[2] W=-d2lp; Psi_new = np.dot(alpha.T,(f-m))/2. - lp.sum() # objective for last alpha vargout = -Tools.general.feval(likfunc, hyp.lik, y, m, None, inffunc, None, 1) Psi_def = vargout[0] # objective for default init f==m if Psi_def < Psi_new: # if default is better, we use it alpha = np.zeros((n,1)); f = mvmK(alpha,V,d0) + m vargout = Tools.general.feval(likfunc, hyp.lik, y, f, None, inffunc, None, 3) lp = vargout[0]; dlp = vargout[1]; d2lp = vargout[2] W=-d2lp; Psi_new = -lp.sum() isWneg = np.any(W<0) # flag indicating whether we found negative values of W it = 0 # this happens for the Student's t likelihood while (Psi_old - Psi_new > tol) and it2: # do we want derivatives? dnlZ = dnlzStruct(hyp) # allocate space for derivatives [d,P,R] = fitcRefresh(d0,Ku,R0,V,W) # g = diag(inv(inv(K)+W))/2 g = d/2 + 0.5*np.atleast_2d((np.dot(np.dot(R,R0),P)**2).sum(axis=0)).T t = W/(1+W*d0) dfhat = g*d3lp # deriv. of nlZ wrt. fhat: dfhat=diag(inv(inv(K)+W)).*d3lp/2 for ii in range(len(hyp.cov)): # covariance hypers ddiagK,dKuu,dKu = Tools.general.feval(covfunc, hyp.cov, x, None, ii) # eval cov derivatives dA = 2.*dKu.T-np.dot(R0tV.T,dKuu) # dQ = dA*R0tV w = np.atleast_2d((dA*R0tV.T).sum(axis=1)).T; v = ddiagK-w # w = diag(dQ); v = diag(dK)-diag(dQ); dnlZ.cov[ii] = np.dot(ddiagK.T,t) - np.dot((RVdd*RVdd).sum(axis=0),v) # explicit part dnlZ.cov[ii] -= (np.dot(RVdd,dA)*np.dot(RVdd,R0tV.T)).sum() # explicit part dnlZ.cov[ii] = 0.5*dnlZ.cov[ii] - np.dot(alpha.T,np.dot(dA,np.dot(R0tV,alpha))+v*alpha)/2. # explicit b = np.dot(dA,np.dot(R0tV,dlp)) + v*dlp # b-K*(Z*b) = inv(eye(n)+K*diag(W))*b KZb = mvmK(mvmZ(b,RVdd,t),V,d0) dnlZ.cov[ii] -= np.dot(dfhat.T,(b-KZb)) # implicit part for ii in range(len(hyp.lik)): # likelihood hypers vargout = Tools.general.feval(likfunc,hyp.lik,y,f,None,inffunc,ii,3) lp_dhyp = vargout[0]; dlp_dhyp = vargout[1]; d2lp_dhyp = vargout[2] dnlZ.lik[ii] = -np.dot(g.T,d2lp_dhyp) - lp_dhyp.sum() # explicit part b = mvmK(dlp_dhyp,V,d0) # implicit part dnlZ.lik[ii] -= np.dot(dfhat.T,b-mvmK(mvmZ(b,RVdd,t),V,d0)) if ii == len(hyp.lik)-1: # since snu2 is a fixed fraction of sn2, there is a covariance-like term in the derivative as well snu = np.sqrt(snu2); T = chol_inv(Kuu + snu2*np.eye(nu)); T = np.dot(T.T,np.dot(T,snu*Ku)); t = np.array([(T*T).sum(axis=0)]).T z = np.dot(alpha.T,np.dot(T.T,np.dot(T,alpha))-t*alpha) - np.dot(np.array([(RVdd*RVdd).sum(axis=0)]),t) z += (np.dot(RVdd,T.T)**2).sum() b = (t*dlp-np.dot(T.T,np.dot(T,dlp)))/2. KZb = mvmK(mvmZ(b,RVdd,t),V,d0) z -= np.dot(dfhat.T,b-KZb) dnlZ.lik[ii] += z for ii in range(len(hyp.mean)): # mean hypers dm = Tools.general.feval(meanfunc, hyp.mean, x, ii) dnlZ.mean[ii] = -np.dot(alpha.T,dm) # explicit part Zdm = mvmZ(dm,RVdd,t) dnlZ.mean[ii] -= np.dot(dfhat.T,(dm-mvmK(Zdm,V,d0))) # implicit part vargout = [post,nlZ[0,0],dnlZ] else: vargout = [post,nlZ[0,0]] return vargout # matrix vector multiplication with Z=inv(K+inv(W)) def mvmZ(x,RVdd,t): Zx = t*x - np.dot(RVdd.T,np.dot(RVdd,x)) return Zx # matrix vector multiplication with approximate covariance matrix def mvmK(al,V,d0): Kal = np.dot(V.T,np.dot(V,al)) + d0*al return Kal # criterion Psi at alpha + s*dalpha for line search def Psi_lineFITC(s,dalpha,alpha,hyp,V,d0,m,likfunc,y,inffunc): alpha = alpha + s*dalpha f = mvmK(alpha,V,d0) + m vargout = Tools.general.feval(likfunc,hyp.lik,y,f,None,inffunc,None,3) lp = vargout[0]; dlp = vargout[1]; d2lp = vargout[2] W = -d2lp Psi = np.dot(alpha.T,(f-m))/2. - lp.sum() return Psi[0],alpha,f,dlp,W # refresh the representation of the posterior from initial and site parameters # to prevent possible loss of numerical precision after many epfitcUpdates # effort is O(n*nu^2) provided that nu nlZ0: # if zero is better .. ttau = np.zeros((n,1)) # .. then initialize with zero instead tnu = np.zeros((n,1)) [d,P,R,nn,gg] = epfitcRefresh(d0,Ku,R0,V,ttau,tnu) # initial repres. nlZ = nlZ0 nlZ_old = np.inf; sweep = 0 # converged, max. sweeps or min. sweeps? while (np.abs(nlZ-nlZ_old) > tol and sweep < max_sweep) or (sweep < min_sweep): nlZ_old = nlZ; sweep += 1 rperm = range(n)#randperm(n) for ii in rperm: # iterate EP updates (in random order) over examples p_i = np.reshape(P[:,ii],(P.shape[0],1)); t = np.dot(R,np.dot(R0,p_i)) # temporary variables sigma_i = d[ii] + np.dot(t.T,t); mu_i = nn[ii] + np.dot(p_i.T,gg) # post moments O(nu^2) tau_ni = 1/sigma_i - ttau[ii] # first find the cavity distribution .. nu_ni = mu_i/sigma_i + m[ii]*tau_ni - tnu[ii] # .. params tau_ni and nu_ni # compute the desired derivatives of the indivdual log partition function vargout = Tools.general.feval(likfunc, hyp.lik, y[ii], nu_ni/tau_ni, 1/tau_ni, inffunc, None, 3) lZ = vargout[0]; dlZ = vargout[1]; d2lZ = vargout[2] ttau_i = -d2lZ /(1.+d2lZ/tau_ni) ttau_i = max(ttau_i,0) # enforce positivity i.e. lower bound ttau by zero tnu_i = ( dlZ + (m[ii]-nu_ni/tau_ni)*d2lZ )/(1.+d2lZ/tau_ni) [d,P[:,ii],R,nn,gg,ttau,tnu] = epfitcUpdate(d,P[:,ii],R,nn,gg,ttau,tnu,ii,ttau_i,tnu_i,m,d0,Ku,R0)# update representation # recompute since repeated rank-one updates can destroy numerical precision [d,P,R,nn,gg] = epfitcRefresh(d0,Ku,R0,V,ttau,tnu) [nlZ,nu_n,tau_n] = epfitcZ(d,P,R,nn,gg,ttau,tnu,d0,R0,Ku,y,likfunc,hyp,m,inffunc) if sweep == max_sweep: raise Exception('maximum number of sweeps reached in function infEP') infFITC_EP.last_ttau = ttau; infFITC_EP.last_tnu = tnu # remember for next call post = postStruct() post.sW = np.sqrt(ttau); # unused for FITC_EP prediction with gp.py dd = 1/(d0+1/ttau) alpha = tnu/ttau*dd; RV = np.dot(R,V); R0tV = np.dot(R0.T,V) alpha = alpha - np.dot(RV.T,np.dot(RV,alpha))*dd # long alpha vector for ordinary infEP post.alpha = np.dot(R0tV,alpha) B = R0tV*np.tile(dd.T,(nu,1)); L = np.dot(B,R0tV.T); B = np.dot(B,RV.T) post.L = np.dot(B,B.T) - L if nargout>2: # do we want derivatives? dnlZ = dnlzStruct(hyp) # allocate space for derivatives RVdd = RV*np.tile(dd.T,(nu,1)) for ii in range(len(hyp.cov)): ddiagK,dKuu,dKu = Tools.general.feval(covfunc, hyp.cov, x, None, ii) dA = 2*dKu.T - np.dot(R0tV.T,dKuu) w = np.atleast_2d ((dA*R0tV.T).sum(axis=1)).T; v = ddiagK - w z = np.dot(dd.T,(v+w)) - np.dot(np.atleast_2d((RVdd*RVdd).sum(axis=0)), v) \ - (np.dot(RVdd,dA).T * np.dot(R0tV,RVdd.T)).sum() dnlZ.cov[ii] = (z - np.dot(alpha.T,(alpha*v)) - np.dot(np.dot(alpha.T,dA),np.dot(R0tV,alpha)))/2. for ii in range(len(hyp.lik)): # likelihood hypers dlik = Tools.general.feval(likfunc, hyp.lik, y, nu_n/tau_n+m, 1/tau_n, inffunc, ii, 1) dnlZ.lik[ii] = -dlik.sum() if ii == len(hyp.lik-1): # since snu2 is a fixed fraction of sn2, there is a covariance-like term in the derivative as well v = np.atleast_2d((R0tV*R0tV).sum(axis=0)).T z = (np.dot(RVdd,R0tV.T)**2).sum() - np.dot(np.atleast_2d((RVdd*RVdd).sum(axis=0)), v) z = z + np.dot(post.alpha.T,post.alpha) - np.dot(alpha.T,(v*alpha)) dnlZ.lik[ii] += snu2*z; [junk,dlZ] = Tools.general.feval(likfunc, hyp.lik, y, nu_n/tau_n, 1/tau_n, inffunc, None, 2) # mean hyps for ii in range(len(hyp.mean)): dm = Tools.general.feval(meanfunc, hyp.mean, x, ii) dnlZ.mean[ii] = -np.dot(dlZ.T,dm) vargout = [post, nlZ[0][0], dnlZ] else: vargout = [post, nlZ[0][0]] return vargout # refresh the representation of the posterior from initial and site parameters # to prevent possible loss of numerical precision after many epfitcUpdates # effort is O(n*nu^2) provided that nu0: R = cholupdate(R,np.sqrt(r)*np.dot(R.T,v),'-') else: R = cholupdate(R,np.sqrt(-r)*np.dot(R.T,v),'+') ttemp = np.dot(R0.T,np.dot(R.T,np.dot(R,np.dot(R0,P_i)))) gg = gg + ((dbi-dwi*(hi-m[ii]))/t) * np.reshape(ttemp,(ttemp.shape[0],1)) # O(nu^2) w[ii] = w_i; b[ii] = b_i; # update site parameters O(1) P_i = P_i/t # O(nu) return d,P_i,R,nn,gg,w,b def infEP(hyp, meanfunc, covfunc, likfunc, x, y, nargout=1): # Expectation Propagation approximation to the posterior Gaussian Process. # The function takes a specified covariance function (see kernels.py) and # likelihood function (see lik.py), and is designed to be used with # gp.py. In the EP algorithm, the sites are # updated in random order, for better performance when cases are ordered # according to the targets. tol = 1e-4; max_sweep = 10; min_sweep = 2 # tolerance to stop EP iterations inffunc = 'inf.infEP' n = x.shape[0] K = Tools.general.feval(covfunc, hyp.cov, x) # evaluate the covariance matrix m = Tools.general.feval(meanfunc, hyp.mean, x) # evaluate the mean vector # A note on naming: variables are given short but descriptive names in # accordance with Rasmussen & Williams "GPs for Machine Learning" (2006): mu # and s2 are mean and variance, nu and tau are natural parameters. A leading t # means tilde, a subscript _ni means "not i" (for cavity parameters), or _n # for a vector of cavity parameters. # marginal likelihood for ttau = tnu = zeros(n,1); equals n*log(2) for likCum* nlZ0 = -Tools.general.feval(likfunc, hyp.lik, y, m, np.reshape(np.diag(K),(np.diag(K).shape[0],1)), inffunc).sum() if "last_ttau" not in infEP.__dict__: # find starting point for tilde parameters ttau = np.zeros((n,1)) # initialize to zero if we have no better guess tnu = np.zeros((n,1)) Sigma = K # initialize Sigma and mu, the parameters of .. mu = np.zeros((n,1)) # .. the Gaussian posterior approximation nlZ = nlZ0 else: ttau = infEP.last_ttau # try the tilde values from previous call tnu = infEP.last_tnu [Sigma, mu, nlZ, L] = epComputeParams(K, y, ttau, tnu, likfunc, hyp, m, inffunc) if nlZ > nlZ0: # if zero is better .. ttau = np.zeros((n,1)) # .. then initialize with zero instead tnu = np.zeros((n,1)) Sigma = K # initialize Sigma and mu, the parameters of .. mu = np.zeros((n,1)) # .. the Gaussian posterior approximation nlZ = nlZ0 nlZ_old = np.inf; sweep = 0 # converged, max. sweeps or min. sweeps? while (np.abs(nlZ-nlZ_old) > tol and sweep < max_sweep) or (sweep < min_sweep): nlZ_old = nlZ; sweep += 1 rperm = range(n) #randperm(n) for ii in rperm: # iterate EP updates (in random order) over examples tau_ni = 1/Sigma[ii,ii] - ttau[ii] # first find the cavity distribution .. nu_ni = mu[ii]/Sigma[ii,ii] + m[ii]*tau_ni - tnu[ii] # .. params tau_ni and nu_ni # compute the desired derivatives of the indivdual log partition function vargout = Tools.general.feval(likfunc, hyp.lik, y[ii], nu_ni/tau_ni, 1/tau_ni, inffunc, None, 3) lZ = vargout[0]; dlZ = vargout[1]; d2lZ = vargout[2] ttau_old = copy(ttau[ii]) # then find the new tilde parameters, keep copy of old ttau[ii] = -d2lZ /(1.+d2lZ/tau_ni) ttau[ii] = max(ttau[ii],0) # enforce positivity i.e. lower bound ttau by zero tnu[ii] = ( dlZ + (m[ii]-nu_ni/tau_ni)*d2lZ )/(1.+d2lZ/tau_ni) ds2 = ttau[ii] - ttau_old # finally rank-1 update Sigma .. si = np.reshape(Sigma[:,ii],(Sigma.shape[0],1)) Sigma = Sigma - ds2/(1.+ds2*si[ii])*np.dot(si,si.T) # takes 70# of total time mu = np.dot(Sigma,tnu) # .. and recompute mu # recompute since repeated rank-one updates can destroy numerical precision [Sigma, mu, nlZ, L] = epComputeParams(K, y, ttau, tnu, likfunc, hyp, m, inffunc) if sweep == max_sweep: raise Exception('maximum number of sweeps reached in function infEP') infEP.last_ttau = ttau; infEP.last_tnu = tnu # remember for next call sW = np.sqrt(ttau); alpha = tnu-sW*solve_chol(L,sW*np.dot(K,tnu)) post = postStruct() post.alpha = alpha # return the posterior params post.sW = sW post.L = L if nargout > 1: if nargout > 2: # do we want derivatives? dnlZ = dnlzStruct(hyp) # allocate space for derivatives ssi = np.sqrt(ttau) V = np.linalg.solve(L.T,np.tile(ssi,(1,n))*K) Sigma = K - np.dot(V.T,V) mu = np.dot(Sigma,tnu) Dsigma = np.reshape(np.diag(Sigma),(np.diag(Sigma).shape[0],1)) tau_n = 1/Dsigma-ttau # compute the log marginal likelihood nu_n = mu/Dsigma-tnu # vectors of cavity parameters F = np.dot(alpha,alpha.T) - np.tile(sW,(1,n))* \ solve_chol(L,np.diag(np.reshape(sW,(sW.shape[0],)))) # covariance hypers for ii in range(len(hyp.cov)): dK = Tools.general.feval(covfunc, hyp.cov, x, None, ii) dnlZ.cov[ii] = -(F*dK).sum()/2. for ii in range(len(hyp.lik)): dlik = Tools.general.feval(likfunc, hyp.lik, y, nu_n/tau_n, 1/tau_n, inffunc, ii, 1) dnlZ.lik[ii] = -dlik.sum() [junk,dlZ] = Tools.general.feval(likfunc, hyp.lik, y, nu_n/tau_n, 1/tau_n, inffunc, None, 2) # mean hyps for ii in range(len(hyp.mean)): dm = Tools.general.feval(meanfunc, hyp.mean, x, ii) dnlZ.mean[ii] = -np.dot(dlZ.T,dm)[0,0] vargout = [post, nlZ, dnlZ] else: vargout = [post, nlZ] else: vargout = [post] return vargout # functions to compute the parameters of the Gaussian approximation, Sigma and # mu, and the negative log marginal likelihood, nlZ, from the current site # parameters, ttau and tnu. Also returns L (useful for predictions). def epComputeParams(K, y, ttau, tnu, likfunc, hyp, m, inffunc): n = len(y) # number of training cases ssi = np.sqrt(ttau) # compute Sigma and mu L = np.linalg.cholesky(np.eye(n)+np.dot(ssi,ssi.T)*K).T # L'*L=B=eye(n)+sW*K*sW V = np.linalg.solve(L.T,np.tile(ssi,(1,n))*K) Sigma = K - np.dot(V.T,V) mu = np.dot(Sigma,tnu) Dsigma = np.reshape(np.diag(Sigma),(np.diag(Sigma).shape[0],1)) tau_n = 1/Dsigma - ttau # compute the log marginal likelihood nu_n = mu/Dsigma-tnu + m*tau_n # vectors of cavity parameters lZ = Tools.general.feval(likfunc, hyp.lik, y, nu_n/tau_n, 1/tau_n, inffunc) nlZ = np.log(np.diag(L)).sum() - lZ.sum() - np.dot(tnu.T,np.dot(Sigma,tnu))/2 \ - np.dot((nu_n-m*tau_n).T,((ttau/tau_n*(nu_n-m*tau_n)-2*tnu) / (ttau+tau_n)))/2 \ + (tnu**2/(tau_n+ttau)).sum()/2.- np.log(1.+ttau/tau_n).sum()/2. return Sigma, mu, nlZ[0], L def infExact(hyp, meanfunc, covfunc, likfunc, x, y, nargout=1): # Exact inference for a GP with Gaussian likelihood. Compute a parametrization # of the posterior, the negative log marginal likelihood and its derivatives # w.r.t. the hyperparameters. if not (likfunc[0] == 'lik.likGauss'): # NOTE: no explicit call to likGauss raise Exception ('Exact inference only possible with Gaussian likelihood') n, D = x.shape K = Tools.general.feval(covfunc, hyp.cov, x) # evaluate covariance matrix m = Tools.general.feval(meanfunc, hyp.mean, x) # evaluate mean vector sn2 = np.exp(2.*hyp.lik) # noise variance of likGauss L = np.linalg.cholesky(K/sn2+np.eye(n)).T # Cholesky factor of covariance with noise alpha = solve_chol(L,y-m)/sn2 post = postStruct() post.alpha = alpha # return the posterior parameters post.sW = np.ones((n,1))/np.sqrt(sn2) # sqrt of noise precision vector post.L = L # L = chol(eye(n)+sW*sW'.*K) if nargout>1: # do we want the marginal likelihood? nlZ = np.dot((y-m).T,alpha/2) + np.log(np.diag(L)).sum() + n*np.log(2*np.pi*sn2)/2. # -log marg lik if nargout>2: # do we want derivatives? dnlZ = dnlzStruct(hyp) # allocate space for derivatives Q = solve_chol(L,np.eye(n))/sn2 - np.dot(alpha,alpha.T) # precompute for convenience for ii in range(len(hyp.cov)): dnlZ.cov[ii] = (Q*Tools.general.feval(covfunc, hyp.cov, x, None, ii)).sum()/2. dnlZ.lik = sn2*np.trace(Q) for ii in range(len(hyp.mean)): dnlZ.mean[ii] = np.dot(-Tools.general.feval(meanfunc, hyp.mean, x, ii).T,alpha) return [post, nlZ[0][0], dnlZ] return [post, nlZ[0][0]] return [post] def infFITC(hyp, meanfunc, covfunc, likfunc, x, y, nargout=1): # FITC approximation to the posterior Gaussian process. The function is # equivalent to infExact with the covariance function: # # Kt = Q + G; G = diag(g); g = diag(K-Q); Q = Ku'*inv(Quu)*Ku; # # where Ku and Kuu are covariances w.r.t. to inducing inputs xu, snu2 = sn2/1e6 # is the noise of the inducing inputs and Quu = Kuu + snu2*eye(nu). # We fixed the standard deviation of the inducing inputs snu to be a one per mil # of the measurement noise's standard deviation sn. # The implementation exploits the Woodbury matrix identity # inv(Kt) = inv(G) - inv(G)*V'*inv(eye(nu)+V*inv(G)*V')*V*inv(G) # in order to be applicable to large datasets. The computational complexity # is O(n nu^2) where n is the number of data points x and nu the number of # inducing inputs in xu. # The function takes a specified covariance function (see kernels.py) and # likelihood function (see lik.py), and is designed to be used with # gp.py and in conjunction with covFITC and likGauss. if not (likfunc[0] == 'lik.likGauss'): # NOTE: no explicit call to likGauss raise Exception ('Exact inference only possible with Gaussian likelihood') cov1 = covfunc[0] if not cov1 == ['kernels.covFITC']: raise Exception('Only covFITC supported.') # check cov diagK,Kuu,Ku = Tools.general.feval(covfunc, hyp.cov, x) # evaluate covariance matrix m = Tools.general.feval(meanfunc, hyp.mean, x) # evaluate mean vector n, D = x.shape; nu = Kuu.shape[0] sn2 = np.exp(2*hyp.lik[0]) # noise variance of likGauss snu2 = 1.e-6*sn2 # hard coded inducing inputs noise Luu = np.linalg.cholesky(Kuu+snu2*np.eye(nu)).T # Kuu + snu2*I = Luu'*Luu V = np.linalg.solve(Luu.T,Ku) # V = inv(Luu')*Ku => V'*V = Q g_sn2 = diagK + sn2 - np.array([(V*V).sum(axis=0)]).T # g + sn2 = diag(K) + sn2 - diag(Q) Lu = np.linalg.cholesky(np.eye(nu) + np.dot(V/np.tile(g_sn2.T,(nu,1)),V.T)).T # Lu'*Lu=I+V*diag(1/g_sn2)*V' r = (y-m)/np.sqrt(g_sn2) be = np.linalg.solve(Lu.T,np.dot(V,r/np.sqrt(g_sn2))) iKuu = solve_chol(Luu,np.eye(nu)) # inv(Kuu + snu2*I) = iKuu post = postStruct() post.alpha = np.linalg.solve(Luu,np.linalg.solve(Lu,be)) # return the posterior parameters post.L = solve_chol(np.dot(Lu,Luu),np.eye(nu)) - iKuu # Sigma-inv(Kuu) post.sW = np.ones((n,1))/np.sqrt(sn2) # unused for FITC prediction with gp.py if nargout>1: # do we want the marginal likelihood nlZ = np.log(np.diag(Lu)).sum() + (np.log(g_sn2).sum() + n*np.log(2*np.pi) + np.dot(r.T,r) - np.dot(be.T,be))/2. if nargout>2: # do we want derivatives? dnlZ = dnlzStruct(hyp) # allocate space for derivatives al = r/np.sqrt(g_sn2) - np.dot(V.T,np.linalg.solve(Lu,be))/g_sn2 # al = (Kt+sn2*eye(n))\y B = np.dot(iKuu,Ku); w = np.dot(B,al) W = np.linalg.solve(Lu.T,V/np.tile(g_sn2.T,(nu,1))) for ii in range(len(hyp.cov)): [ddiagKi,dKuui,dKui] = Tools.general.feval(covfunc, hyp.cov, x, None, ii) # eval cov deriv R = 2.*dKui-np.dot(dKuui,B); v = ddiagKi - np.array([(R*B).sum(axis=0)]).T # diag part of cov deriv dnlZ.cov[ii] = ( np.dot(ddiagKi.T,1./g_sn2) + np.dot(w.T,(np.dot(dKuui,w)-2.*np.dot(dKui,al))) - np.dot(al.T,(v*al)) \ - np.dot(np.array([(W*W).sum(axis=0)]),v) - (np.dot(R,W.T)*np.dot(B,W.T)).sum() )/2. dnlZ.lik = sn2*((1./g_sn2).sum() - (np.array([(W*W).sum(axis=0)])).sum() - np.dot(al.T,al)) # since snu2 is a fixed fraction of sn2, there is a covariance-like term in the derivative as well dKuui = 2*snu2 R = -dKuui*B v = -np.array([(R*B).sum(axis=0)]).T # diag part of cov deriv dnlZ.lik += (np.dot(w.T,np.dot(dKuui,w)) -np.dot(al.T,(v*al)) \ - np.dot(np.array([(W*W).sum(axis=0)]),v) - (np.dot(R,W.T)*np.dot(B,W.T)).sum() )/2. dnlZ.lik = dnlZ.lik[0] for ii in range(len(hyp.mean)): dnlZ.mean[ii] = np.dot(-Tools.general.feval(meanfunc, hyp.mean, x, ii).T,al)[0,0] return [post, nlZ[0,0], dnlZ] return [post, nlZ[0,0]] return [post] pyGP_PR-1.0/src/Core/LIK/000077500000000000000000000000001222232424600147275ustar00rootroot00000000000000pyGP_PR-1.0/src/Core/LIK/__init__.py000077500000000000000000000000251222232424600170400ustar00rootroot00000000000000 __all__ = ["lik"] pyGP_PR-1.0/src/Core/LIK/lik.py000066400000000000000000000434301222232424600160640ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 20/05/2013 #================================================================================ # likelihood functions are provided to be used by the gp.py function: # # likErf (Error function, classification, probit regression) # likLogistic [NOT IMPLEMENTED!] (Logistic, classification, logit regression) # likUni [NOT IMPLEMENTED!] (Uniform likelihood, classification) # # likGauss (Gaussian, regression) # likLaplace [NOT IMPLEMENTED!] (Laplacian or double exponential, regression) # likSech2 [NOT IMPLEMENTED!] (Sech-square, regression) # likT [NOT IMPLEMENTED!] (Student's t, regression) # # likPoisson [NOT IMPLEMENTED!] (Poisson regression, count data) # # likMix [NOT IMPLEMENTED!] (Mixture of individual covariance functions) # # The likelihood functions have three possible modes, the mode being selected # as follows (where "lik" stands for any likelihood function defined as lik*): # # 1) With one or no input arguments: [REPORT NUMBER OF HYPERPARAMETERS] # # s = lik OR s = lik(hyp) # # The likelihood function returns a string telling how many hyperparameters it # expects, using the convention that "D" is the dimension of the input space. # For example, calling "likLogistic" returns the string '0'. # # # 2) With three or four input arguments: [PREDICTION MODE] # # lp = lik(hyp, y, mu) OR [lp, ymu, ys2] = lik(hyp, y, mu, s2) # # This allows to evaluate the predictive distribution. Let p(y_*|f_*) be the # likelihood of a test point and N(f_*|mu,s2) an approximation to the posterior # marginal p(f_*|x_*,x,y) as returned by an inference method. The predictive # distribution p(y_*|x_*,x,y) is approximated by. # q(y_*) = \int N(f_*|mu,s2) p(y_*|f_*) df_* # # lp = log( q(y) ) for a particular value of y, if s2 is [] or 0, this # corresponds to log( p(y|mu) ) # ymu and ys2 the mean and variance of the predictive marginal q(y) # note that these two numbers do not depend on a particular # value of y # All vectors have the same size. # # # 3) With five or six input arguments, the fifth being a string [INFERENCE MODE] # # [varargout] = lik(hyp, y, mu, s2, inf) OR # [varargout] = lik(hyp, y, mu, s2, inf, i) # # There are three cases for inf, namely a) infLaplace, b) infEP and c) infVB. # The last input i, refers to derivatives w.r.t. the ith hyperparameter. # # a1) [lp,dlp,d2lp,d3lp] = lik(hyp, y, f, [], 'infLaplace') # lp, dlp, d2lp and d3lp correspond to derivatives of the log likelihood # log(p(y|f)) w.r.t. to the latent location f. # lp = log( p(y|f) ) # dlp = d log( p(y|f) ) / df # d2lp = d^2 log( p(y|f) ) / df^2 # d3lp = d^3 log( p(y|f) ) / df^3 # # a2) [lp_dhyp,dlp_dhyp,d2lp_dhyp] = lik(hyp, y, f, [], 'infLaplace', i) # returns derivatives w.r.t. to the ith hyperparameter # lp_dhyp = d log( p(y|f) ) / ( dhyp_i) # dlp_dhyp = d^2 log( p(y|f) ) / (df dhyp_i) # d2lp_dhyp = d^3 log( p(y|f) ) / (df^2 dhyp_i) # # # b1) [lZ,dlZ,d2lZ] = lik(hyp, y, mu, s2, 'infEP') # let Z = \int p(y|f) N(f|mu,s2) df then # lZ = log(Z) # dlZ = d log(Z) / dmu # d2lZ = d^2 log(Z) / dmu^2 # # b2) [dlZhyp] = lik(hyp, y, mu, s2, 'infEP', i) # returns derivatives w.r.t. to the ith hyperparameter # dlZhyp = d log(Z) / dhyp_i # # # c1) [h,b,dh,db,d2h,d2b] = lik(hyp, y, [], ga, 'infVB') # ga is the variance of a Gaussian lower bound to the likelihood p(y|f). # p(y|f) \ge exp( b*f - f.^2/(2*ga) - h(ga)/2 ) \propto N(f|b*ga,ga) # The function returns the linear part b and the "scaling function" h(ga) and derivatives # dh = d h/dga # db = d b/dga # d2h = d^2 h/dga # d2b = d^2 b/dga # # c2) [dhhyp] = lik(hyp, y, [], ga, 'infVB', i) # dhhyp = dh / dhyp_i, the derivative w.r.t. the ith hyperparameter # # Cumulative likelihoods are designed for binary classification. Therefore, they # only look at the sign of the targets y; zero values are treated as +1. # # Some examples for valid likelihood functions: # lik = @likLogistic; # lik = {'likMix',{'likUni',@likErf}} # lik = {@likPoisson,'logistic'}; # # See the documentation for the individual likelihood for the computations specific # to each likelihood function. # # @author: Daniel Marthaler (Fall 2012) # Substantial updates by Shan Huang (Sep. 2013) # # This is a python implementation of gpml functionality (Copyright (c) by # Carl Edward Rasmussen and Hannes Nickisch, 2013-01-21). # # Copyright (c) by Marion Neumann and Daniel Marthaler, 20/05/2013 import numpy as np import Tools.general from scipy.special import erf def likErf(hyp=None, y=None, mu=None, s2=None, inffunc=None, der=None, nargout=None): # function [varargout] = likErf(hyp, y, mu, s2, inf, i) # likErf - Error function or cumulative Gaussian likelihood function for binary # classification or probit regression. The expression for the likelihood is # # likErf(t) = (1+erf(t/sqrt(2)))/2 = normcdf(t). # # Several modes are provided, for computing likelihoods, derivatives and moments # respectively, see lik.py for the details. In general, care is taken # to avoid numerical issues when the arguments are extreme. if mu == None: return [0] # report number of hyperparameters if not y == None: y = np.sign(y) y[y==0] = 1 else: y = 1; # allow only +/- 1 values if inffunc == None: # prediction mode if inf is not present y = y*np.ones_like(mu) # make y a vector s2zero = True; if not s2 == None: if np.linalg.norm(s2)>0: s2zero = False # s2==0 ? if s2zero: # log probability evaluation [p,lp] = cumGauss(y,mu,2) else: # prediction lp = Tools.general.feval(['lik.likErf'],hyp, y, mu, s2, 'inf.infEP',None,1) p = np.exp(lp) if nargout>1: ymu = 2*p-1 # first y moment if nargout>2: ys2 = 4*p*(1-p) # second y moment varargout = [lp,ymu,ys2] else: varargout = [lp,ymu] else: varargout = lp else: # inference mode if inffunc == 'inf.infLaplace': if der == None: # no derivative mode f = mu; yf = y*f # product latents and labels [p,lp] = cumGauss(y,f,2) if nargout>1: # derivative of log likelihood n_p = gauOverCumGauss(yf,p) dlp = y*n_p # derivative of log likelihood if nargout>2: # 2nd derivative of log likelihood d2lp = -n_p**2 - yf*n_p if nargout>3: # 3rd derivative of log likelihood d3lp = 2*y*n_p**3 + 3*f*n_p**2 + y*(f**2-1)*n_p varargout = [lp,dlp,d2lp,d3lp] else: varargout = [lp,dlp,d2lp] else: varargout = [lp,dlp] else: varargout = lp else: # derivative mode varargout = nargout*[] # derivative w.r.t. hypers if inffunc == 'inf.infEP': if der == None: # no derivative mode z = mu/np.sqrt(1+s2) [junk,lZ] = cumGauss(y,z,2) # log part function if not y == None: z = z*y if nargout>1: if y == None: y = 1 n_p = gauOverCumGauss(z,np.exp(lZ)) dlZ = y*n_p/np.sqrt(1.+s2) # 1st derivative wrt mean if nargout>2: d2lZ = -n_p*(z+n_p)/(1.+s2) # 2nd derivative wrt mean varargout = [lZ,dlZ,d2lZ] else: varargout = [lZ,dlZ] else: varargout = lZ else: # derivative mode varargout = 0 # deriv. wrt hyp.lik #end if inffunc == 'inf.infVB': if der == None: # no derivative mode # naive variational lower bound based on asymptotical properties of lik # normcdf(t) -> -(t*A_hat^2-2dt+c)/2 for t->-np.inf (tight lower bound) d = 0.158482605320942; c = -1.785873318175113; ga = s2; n = len(ga); b = d*y*np.ones((n,1)); db = np.zeros((n,1)); d2b = db h = -2.*c*np.ones((n,1)); h[ga>1] = np.inf; dh = np.zeros((n,1)); d2h = dh varargout = [h,b,dh,db,d2h,d2b] else: # derivative mode varargout = [] # deriv. wrt hyp.lik return varargout def cumGauss(y=None,f=None,nargout=1): if not y == None: yf = y*f else: yf = f # product of latents and labels p = (1. + erf(yf/np.sqrt(2.)))/2. # likelihood if nargout>1: lp = logphi(yf,p) return p,lp else: return p # safe implementation of the log of phi(x) = \int_{-\infty}^x N(f|0,1) df # logphi(z) = log(normcdf(z)) def logphi(z,p): lp = np.zeros_like(z) # initialize zmin = -6.2; zmax = -5.5; ok = z>zmax # safe evaluation for large values bd = z-5 # naive evaluation for large values of f n_p[ok] = (np.exp(-f[ok]**2/2)/np.sqrt(2*np.pi)) / p[ok] bd = f<-6 # tight upper bound evaluation n_p[bd] = np.sqrt(f[bd]**2/4+1)-f[bd]/2 interp = np.logical_and(np.logical_not(ok),np.logical_not(bd)) # linearly interpolate between both of them tmp = f[interp] lam = -5. - f[interp] n_p[interp] = (1-lam)*(np.exp(-tmp**2/2)/np.sqrt(2*np.pi))/p[interp] + \ lam *(np.sqrt(tmp**2/4+1)-tmp/2); return n_p def likGauss(hyp=None, y=None, mu=None, s2=None, inffunc=None, der=None, nargout=1): # likGauss - Gaussian likelihood function for regression. The expression for the likelihood is # # likGauss(t) = exp(-(t-y)^2/2*sn^2) / sqrt(2*pi*sn^2), # # where y is the mean and sn is the standard deviation. # # The hyperparameters are: # hyp = [ log(sn) ] # # Several modes are provided, for computing likelihoods, derivatives and moments # respectively, see lik.py for the details. In general, care is taken # to avoid numerical issues when the arguments are extreme. if mu == None: return [1] # report number of hyperparameters sn2 = np.exp(2.*hyp) if inffunc == None: # prediction mode if inffunc is not present if y == None: y = np.zeros_like(mu) s2zero = True if not (s2 == None): if np.linalg.norm(s2) > 0: s2zero = False if s2zero: # s2==0 ? lp = -(y-mu)**2 /sn2/2 - np.log(2.*np.pi*sn2)/2. # log probability s2 = 0. else: lp = Tools.general.feval(['lik.likGauss'],hyp, y, mu, s2, 'inf.infEP',None,1) # prediction if nargout>1: ymu = mu; # first y moment if nargout>2: ys2 = s2 + sn2; # second y moment varargout = [lp,ymu,ys2] else: varargout = [lp,ymu] else: varargout = lp else: if inffunc == 'inf.infLaplace': if der == None: # no derivative mode if y == None: y=0 #end ymmu = y-mu lp = -ymmu**2/(2*sn2) - np.log(2*np.pi*sn2)/2. if nargout>1: dlp = ymmu/sn2 # dlp, derivative of log likelihood if nargout>2: # d2lp, 2nd derivative of log likelihood d2lp = -np.ones_like(ymmu)/sn2 if nargout>3: # d3lp, 3rd derivative of log likelihood d3lp = np.zeros_like(ymmu) varargout = [lp,dlp,d2lp,d3lp] else: varargout = [lp,dlp,d2lp] else: varargout = [lp,dlp] else: varargout = lp else: # derivative mode lp_dhyp = (y-mu)**2/sn2 - 1 # derivative of log likelihood w.r.t. hypers dlp_dhyp = 2*(mu-y)/sn2 # first derivative, d2lp_dhyp = 2*np.ones_like(mu)/sn2 # and also of the second mu derivative varargout = [lp_dhyp,dlp_dhyp,d2lp_dhyp] elif inffunc == 'inf.infEP': if der == None: # no derivative mode lZ = -(y-mu)**2/(sn2+s2)/2. - np.log(2*np.pi*(sn2+s2))/2. # log part function if nargout>1: dlZ = (y-mu)/(sn2+s2) # 1st derivative w.r.t. mean if nargout>2: d2lZ = -1/(sn2+s2) # 2nd derivative w.r.t. mean varargout = [lZ,dlZ,d2lZ] else: varargout = [lZ,dlZ] else: varargout = lZ else: # derivative mode dlZhyp = ((y-mu)**2/(sn2+s2)-1) / (1+s2/sn2) # deriv. w.r.t. hyp.lik varargout = dlZhyp elif inffunc == 'inf.infVB': if der == None: # variational lower site bound # t(s) = exp(-(y-s)^2/2sn2)/sqrt(2*pi*sn2) # the bound has the form: b*s - s.^2/(2*ga) - h(ga)/2 with b=y/ga ga = s2; n = len(ga); b = y/ga; y = y*np.ones((n,1)) db = -y/ga**2 d2b = 2*y/ga**3 h = np.zeros((n,1)); dh = h; d2h = h # allocate memory for return args id = (ga <= sn2 + 1e-8) # OK below noise variance h[id] = y[id]**2/ga[id] + np.log(2*np.pi*sn2); h[np.logical_not(id)] = np.inf dh[id] = -y[id]**2/ga[id]**2 d2h[id] = 2*y[id]**2/ga[id]**3 id = ga < 0; h[id] = np.inf; dh[id] = 0; d2h[id] = 0 # neg. var. treatment varargout = [h,b,dh,db,d2h,d2b] else: ga = s2; n = len(ga); dhhyp = np.zeros((n,1)); dhhyp[ga<=sn2] = 2 dhhyp[ga<0] = 0 # negative variances get a special treatment varargout = dhhyp # deriv. w.r.t. hyp.lik else: raise Exception('Incorrect inference in lik.Gauss\n'); return varargout pyGP_PR-1.0/src/Core/MEAN/000077500000000000000000000000001222232424600150305ustar00rootroot00000000000000pyGP_PR-1.0/src/Core/MEAN/__init__.py000077500000000000000000000000271222232424600171430ustar00rootroot00000000000000 __all__ = ["means"] pyGP_PR-1.0/src/Core/MEAN/means.py000077500000000000000000000330561222232424600165170ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 20/05/2013 #================================================================================ # mean functions to be use by Gaussian process functions. There are two # different kinds of mean functions: simple and composite: # # simple mean functions: # # meanZero - zero mean function # meanOne - one mean function # meanConst - constant mean function # meanLinear - linear mean function # # composite covariance functions (see explanation at the bottom): # # meanScale - scaled version of a mean function # meanPow - power of a mean function # meanProd - products of mean functions # meanSum - sums of mean functions # meanMask - mask some dimensions of the data # # Naming convention: all mean functions are named "mean*". # # # 1) With no or only a single input argument: # # s = meanNAME or s = meanNAME(hyp) # # The mean function returns a string s telling how many hyperparameters hyp it # expects, using the convention that "D" is the dimension of the input space. # For example, calling "meanLinear" returns the string 'D'. # # 2) With two input arguments: # # m = meanNAME(hyp, x) # # The function computes and returns the mean vector where hyp are the # hyperparameters and x is an n by D matrix of cases, where D is the dimension # of the input space. The returned mean vector is of size n by 1. # # 3) With three input arguments: # # dm = meanNAME(hyp, x, i) # # The function computes and returns the n by 1 vector of partial derivatives # of the mean vector w.r.t. hyp(i) i.e. hyperparameter number i. # # # @author: Daniel Marthaler (Fall 2012) # Substantial updates by Shan Huang (Sep. 2013) # # This is a python implementation of gpml functionality (Copyright (c) by # Carl Edward Rasmussen and Hannes Nickisch, 2013-01-21). # # Copyright (c) by Marion Neumann and Daniel Marthaler, 20/05/2013 import Tools import numpy as np import math def meanMask(meanfunc, meanhyper=None, x=None, der=None): # meanMask - compose a mean function as another mean # function (meanfunc), but with only a subset of dimensions of x. mask = meanfunc[0] # The indicies to be masked (should be a list of integers) mean = meanfunc[1] # covariance function to $ if meanhyper == None: # report number of parameters A = Tools.general.feval(meanfunc[1]) return [A] n, D = x.shape assert(len(mask) <= D) assert(max(mask) <= D) assert(min(mask) >= 0) if der == None: # compute covariance matix for dataset x A = Tools.general.feval(mean, meanhyper, x[:,mask]) else: # compute derivatives A = Tools.general.feval(mean, meanhyper, x[:,mask], der) return A def meanConst(meanhyper=None, x=None, der=None): # Constant mean function. The mean function is parameterized as: # m(x) = c # # The hyperparameter is: # meanhyper = [ c ] if meanhyper == None: # report number of parameters return [1] n,D = x.shape c = meanhyper[0] if der == None: # evaluate mean A = c*np.ones((n,1)) elif isinstance(der, int) and der == 0: # compute derivative vector wrt c A = np.ones((n,1)) else: A = np.zeros((n,1)) return A def meanLinear(meanhyper=None, x=None, der=None): # Linear mean function. The mean function is parameterized as: # m(x) = sum_i ai * x_i # # The hyperparameter is: # meanhyper = [ a_1 # a_2 # ... # a_D ] if meanhyper == None: # report number of parameters return ['D'] n, D = x.shape c = np.array(meanhyper) c = np.reshape(c,(len(c),1)) if der == None: # evaluate mean A = np.dot(x,c) elif isinstance(der, int) and der < D: # compute derivative vector wrt meanparameters A = np.reshape(x[:,der], (len(x[:,der]),1) ) else: A = np.zeros((n,1)) return A def meanOne(meanhyper=None, x=None, der=None): # One mean function. The mean function does not have any parameters # m(x) = 1 if meanhyper == None: # report number of parameters return [0] n,D = x.shape if der == None: # evaluate mean A = np.ones((n,1)) else: A = np.zeros((n,1)) return A def meanZero(meanhyper=None, x=None, der=None): ''' Zero mean function. The mean function does not have any parameters m(x) = 1 ''' if meanhyper == None: # report number of parameters return [0] n, D = x.shape A = np.zeros((n,1)) return A def meanProd(meanfunc, meanhyper=None, x=None, der=None): # meanProd - compose a mean function as the product of other mean # functions. This function doesn't actually compute very much on its own, it # merely does some bookkeeping, and calls other mean functions to do the # actual work. def DetermineNumberOfParameters(v,no_param): if isinstance(no_param, int): v.append(no_param) elif isinstance(no_param,str): # no_param is a string pram_str = no_param.split(' ') if pram_str[0]=='D': temp = int(D) if pram_str[1]=='+': temp += int(pram_str[2]) elif pram_str[1]=='-': temp -= int(pram_str[2]) else: raise Exception(["Error: number of parameters of "+meanfunc[i] +" unknown!"]) v.append(temp) elif isinstance(no_param, list): # The number of hyperparameters for this piece of meanfunc is the sum # of all of them in this composition temp = [0] for jj in xrange(len(no_param)): DetermineNumberOfParameters(temp,no_param[jj]) v.append(sum(temp)) else: # This is an error, we should never be here raise Exception("Error in return of number of parameters") return v if meanhyper == None: # report number of parameters A = [Tools.general.feval(meanfunc[0])] for ii in range(1,len(meanfunc)): A.append(Tools.general.feval(meanfunc[ii])) return A [n, D] = x.shape # SET vector v (v indicates how many parameters each meanfunc has # (NOTE : v[i]=number of parameters + 1 -> this is because of the indexing of python!)) v = [0] for ii in range(1,len(meanfunc)+1): no_param = Tools.general.feval(meanfunc[ii-1]) DetermineNumberOfParameters(v,no_param) if der == None: # compute mean vector A = np.ones((n, 1)) # allocate space for mean vector for ii in range(1,len(meanfunc)+1): # iteration over multiplicand functions f = meanfunc[ii-1] B = Tools.general.feval(f,meanhyper[sum(v[0:ii]):sum(v[0:ii])+v[ii]], x) # accumulate means A *= B elif isinstance(der, int): # compute derivative vector tmp = 0 A = np.ones((n, 1)) # allocate space for derivative vector flag = True for ii in range(1,len(meanfunc)+1): tmp += v[ii] if der1: if pram_str[1]=='+': temp += int(pram_str[2]) elif pram_str[1]=='-': temp -= int(pram_str[2]) else: raise Exception(["Error: number of parameters of "+meanfunc[ii][0] +" unknown!"]) v.append(temp) elif isinstance(no_param, list): # The number of hyperparameters for this piece of meanfunc is the sum # of all of them in this composition temp = [0] for jj in xrange(len(no_param)): DetermineNumberOfParameters(temp,no_param[jj]) v.append(sum(temp)) else: # This is an error, we should never be here raise Exception("Error in return of number of parameters") return v if meanhyper == None: # report number of parameters A = [Tools.general.feval(meanfunc[0])] for ii in range(1,len(meanfunc)): A.append(Tools.general.feval(meanfunc[ii])) return A [n, D] = x.shape # SET vector v (v indicates how many parameters each meanfunc has # (NOTE : v[i]=number of parameters + 1 -> this is because of the indexing of python!)) v = [0] # needed for technical reasons for ii in range(1,len(meanfunc)+1): no_param = Tools.general.feval(meanfunc[ii-1]) DetermineNumberOfParameters(v,no_param) if der == None: # compute mean vector A = np.zeros((n, 1)) # allocate space for mean vector for ii in range(1,len(meanfunc)+1): # iteration over summand functions f = meanfunc[ii-1] s = sum(v[0:ii]) A = A + Tools.general.feval(f, meanhyper[s:(s+v[ii])], x) # accumulate means elif isinstance(der, int): # compute derivative vector A = np.zeros((n, 1)) # allocate space for mean vector tmp = 0 for ii in range(1,len(meanfunc)+1): tmp += v[ii] if der (tol2 - 0.5*(b-a)) ): gs = True # Is a parabolic fit possible? if abs(e) > tol1: # Yes, so fit parabola gs = False r = (xf-w)*(fx-fv) q = (xf-v)*(fx-fw) p = (xf-v)*q-(xf-w)*r q = 2.0*(q-r) if q > 0.0: p = -p q = abs(q) r = e; e = d # Is the parabola acceptable? if ( (abs(p)q*(a-xf)) and (p= xm: e = a-xf else: e = b-xf d = c*e # The function must not be evaluated too close to xf si = cmp(d,0) if (d == 0): si += 1 x = xf + si * max(abs(d),tol1) vargout = f(x,*args); fu = vargout[0][0]; varargout = vargout[1:] funccount += 1 # Update a, b, v, w, x, xm, tol1, tol2 if fu <= fx: if x >= xf: a = xf else: b = xf v = w; fv = fw w = xf; fw = fx xf = x; fx = fu else: # fu > fx if x < xf: a = x else: b = x if ( (fu <= fw) or (w == xf) ): v = w; fv = fw w = x; fw = fu elif ( (fu <= fv) or ((v == xf) or (v == w)) ): v = x; fv = fu xm = 0.5*(a+b) tol1 = seps*abs(xf) + tol/3.0; tol2 = 2.0*tol1 if funccount >= Nitmax: #print 'Warning: Maximum number of function evaluations reached (brentmin)' break # check that endpoints are less than the minimum found if ( (fa < fx) and (fa <= fb) ): xf = xlow; fx = fa elif fb < fx: xf = xupp; fx = fb fmin = fx xmin = xf vargout = [xmin,fmin,funccount] for vv in varargout: vargout.append(vv) return vargout pyGP_PR-1.0/src/Core/UTIL/solve_chol.py000066400000000000000000000022721222232424600175770ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 15/01/2010 #================================================================================ import numpy as np def solve_chol(L, B): # solve_chol - solve linear equations from the Cholesky factorization. # Solve A*X = B for X, where A is square, symmetric, positive definite. The # input to the function is R the Cholesky decomposition of A and the matrix B. # Example: X = solve_chol(chol(A),B); try: assert(L.shape[0] == L.shape[1] and L.shape[0] == B.shape[0]) except AssertionError: raise Exception('Wrong sizes of matrix arguments in solve_chol.py'); X = np.linalg.solve(L,np.linalg.solve(L.T,B)) return X pyGP_PR-1.0/src/Core/UTIL/utils.py000066400000000000000000000127471222232424600166120ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 20/05/2013 #================================================================================ # Created on 08/09/2009 # # @author: Marion Neumann (Aug 2009) # Substantial updates by Daniel Marthaler Fall 2012. import numpy as np import matplotlib.pyplot as plt import Tools.general from copy import deepcopy from random import sample class hyperParameters: def __init__(self): self.mean = np.array([]) self.cov = np.array([]) self.lik = np.array([]) class dnlzStruct: def __init__(self, hyp): self.mean = self.cov = self.lik = [] if hyp.mean != None: self.mean = np.zeros(hyp.mean.shape) if hyp.cov != None: self.cov = np.zeros(hyp.cov.shape) if hyp.lik != None: self.lik = np.zeros(hyp.lik.shape) def cholupdate(R,x,sgn='+'): # Placeholder for a python version of MATLAB's cholupdate. Now it is O(n^3) if len(x.shape) == 1: # Reshape x so that the dot product below is correct x = np.reshape(x,(x.shape[0],1)) assert(R.shape[0] == x.shape[0]) A = np.dot(R.T,R) if sgn == '+': R1 = A + np.dot(x,x.T) elif sgn == '-': R1 = A - np.dot(x,x.T) else: raise Exception('Sign needs to be + or - in cholupdate') return np.linalg.cholesky(R1).T def randperm(k): # return a random permutation of range(k) z = range(k) y = [] ii = 0 while z and ii < 2*k: n = sample(z,1)[0] y.append(n) z.remove(n) ii += 1 return y def cool_plotter(xs,ym,ys2,x,y,axisvals=None,file=None): xss = np.reshape(xs,(xs.shape[0],)) ymm = np.reshape(ym,(ym.shape[0],)) ys22 = np.reshape(ys2,(ys2.shape[0],)) fig = plt.figure() ax = fig.add_subplot(111,axisbg=[0.,1.0,0.0,0.2]) #plt.plot(xs, ym, 'g-', x, y, 'r+', linewidth = 3.0, markersize = 10.0) plt.plot(xs, ym, 'r-', x, y, 'y^', linewidth = 3.0, markersize = 5.0) plt.fill_between(xss,ymm + 1.*np.sqrt(ys22), ymm - 1.*np.sqrt(ys22), facecolor=[0.,1.0,0.0,0.8],linewidths=0.0) plt.fill_between(xss,ymm + 2.*np.sqrt(ys22), ymm - 2.*np.sqrt(ys22), facecolor=[0.,1.0,0.0,0.6],linewidths=0.0) plt.fill_between(xss,ymm + 3.*np.sqrt(ys22), ymm - 3.*np.sqrt(ys22), facecolor=[0.,1.0,0.0,0.4],linewidths=0.0) plt.grid() if axisvals: plt.axis(axisvals) #plt.xlabel('input x') #plt.ylabel('output y') #Create figure; attempt to set background to light gray fig.set_facecolor("#000000") #plt.savefig('GP.png') if file and isinstance(file,str): plt.savefig(file) plt.show() def plotter(xs,ym,ys2,x,y,axisvals=None,file=None): xss = np.reshape(xs,(xs.shape[0],)) ymm = np.reshape(ym,(ym.shape[0],)) ys22 = np.reshape(ys2,(ys2.shape[0],)) plt.plot(xs, ym, 'g-', x, y, 'r+', linewidth = 3.0, markersize = 10.0) plt.fill_between(xss,ymm + 2.*np.sqrt(ys22), ymm - 2.*np.sqrt(ys22), facecolor=[0.,1.0,0.0,0.8],linewidths=0.0) plt.grid() if axisvals: plt.axis(axisvals) plt.xlabel('input x') plt.ylabel('output y') if file and isinstance(file,str): plt.savefig(file) plt.show() def FITCplotter(u,xs,ym,ys2,x,y,axisvals=None,file=None): xss = np.reshape(xs,(xs.shape[0],)) ymm = np.reshape(ym,(ym.shape[0],)) ys22 = np.reshape(ys2,(ys2.shape[0],)) plt.plot(xs, ym, 'g-', x, y, 'r+', linewidth = 3.0, markersize = 10.0) plt.fill_between(xss,ymm + 2.*np.sqrt(ys22), ymm - 2.*np.sqrt(ys22), facecolor=[0.,1.0,0.0,0.8],linewidths=0.0) plt.grid() if axisvals: plt.axis(axisvals) plt.xlabel('input x') plt.ylabel('output y') plt.plot(u,np.ones_like(u),'kx',markersize=12) if file and isinstance(file,str): plt.savefig(file) plt.show() def convert_to_array(hyp): y = np.concatenate((np.concatenate((hyp.mean, hyp.cov),axis=0),hyp.lik),axis=0) return y def convert_to_class(x,hyp): y = deepcopy(hyp) Lm = len(hyp.mean) Lc = len(hyp.cov) y.mean = x[:Lm] y.cov = x[Lm:(Lm+Lc)] y.lik = x[(Lm+Lc):] return y def unique(x): # First flatten x y = [item for sublist in x for item in sublist] if isinstance(x,np.ndarray): n,D = x.shape assert(D == 1) y = np.array( list(set(x[:,0])) ) y = np.reshape(y, (len(y),1)) else: y = list(set(y)) return y def flatten(l, ltypes=(list, tuple)): ltype = type(l) l = list(l) i = 0 while i < len(l): while isinstance(l[i], ltypes): if not l[i]: l.pop(i) i -= 1 break else: l[i:i + 1] = l[i] i += 1 return ltype(l) def checkParameters(func,hyp,D): valt = flatten(Tools.general.feval(func)) val = 0 if isinstance(valt,list): for v in valt: if isinstance(v,str): val += eval(v) else: val += v else: val = valt res = (val == len(hyp)) return res pyGP_PR-1.0/src/Core/__init__.py000077500000000000000000000000671222232424600164270ustar00rootroot00000000000000__all__ = ["COV", "MEAN", "LIK", "INF", "UTIL", "gp"] pyGP_PR-1.0/src/Core/demoClassification.py000066400000000000000000000126351222232424600204710ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 30/09/2013 #================================================================================ from gp import gp from UTIL.solve_chol import solve_chol import Tools.general from Tools.min_wrapper import min_wrapper import Tools.nearPD import numpy as np import matplotlib.pyplot as plt import matplotlib.cm as cm import matplotlib.mlab as mlab from UTIL.utils import * from scipy.optimize import fmin_cg as cg if __name__ == '__main__': PLOT = True ## LOAD data demoData = np.load('../../data/classification_data.npz') x = demoData['x'] # training data y = demoData['y'] # training target xstar = demoData['xstar'] # test data n = xstar.shape[0] # number of test points ## DATA only needed in plotting if PLOT: x1 = demoData['x1'] x2 = demoData['x2'] t1 = demoData['t1'] t2 = demoData['t2'] p1 = demoData['p1'] p2 = demoData['p2'] ### PLOT data if PLOT: fig = plt.figure() plt.plot(x1[:,0], x1[:,1], 'b+', markersize = 12) plt.plot(x2[:,0], x2[:,1], 'r+', markersize = 12) pc = plt.contour(t1, t2, np.reshape(p2/(p1+p2), (t1.shape[0],t1.shape[1]) )) fig.colorbar(pc) plt.grid() plt.axis([-4, 4, -4, 4]) plt.show() ## DEFINE parameterized mean and covariance functions meanfunc = ['means.meanConst'] covfunc = ['kernels.covSEard'] ## DEFINE likelihood function used likfunc = ['lik.likErf'] ## SPECIFY inference method inffunc = ['inf.infLaplace'] ## SET (hyper)parameters hyp = hyperParameters() hyp.mean = np.array([-2.842117459073954]) hyp.cov = np.array([0.051885508906388,0.170633324977413,1.218386482861781]) ##----------------------------------------------------------## ## STANDARD GP (example 1) ## ##----------------------------------------------------------## print '...example 1: prediction...' ## GET negative log marginal likelihood [nlml,dnlZ,post] = gp(hyp, inffunc, meanfunc, covfunc, likfunc, x, y, None, None, True) print "nlml = ", nlml ## PREDICTION [ymu,ys2,fmu,fs2,lp,post] = gp(hyp, inffunc, meanfunc, covfunc, likfunc, x, y, xstar, np.ones((n,1)) ) ## PLOT log predictive probabilities if PLOT: fig = plt.figure() plt.plot(x1[:,0], x1[:,1], 'b+', markersize = 12) plt.plot(x2[:,0], x2[:,1], 'r+', markersize = 12) pc = plt.contour(t1, t2, np.reshape(np.exp(lp), (t1.shape[0],t1.shape[1]) )) fig.colorbar(pc) plt.grid() plt.axis([-4, 4, -4, 4]) plt.show() ##----------------------------------------------------------## ## SPARSE GP (example 2) ## ##----------------------------------------------------------## print '...example 2: FITC training and prediction...' ## SPECIFY inducing points u1,u2 = np.meshgrid(np.linspace(-2,2,5),np.linspace(-2,2,5)) u = np.array(zip(np.reshape(u2,(np.prod(u2.shape),)),np.reshape(u1,(np.prod(u1.shape),)))) del u1, u2 nu = u.shape[0] ## SPECIFY FITC covariance function covfuncF = [['kernels.covFITC'], covfunc, u] ## SPECIFY FICT inference method inffunc = ['inf.infFITC_EP'] ## GET negative log marginal likelihood [nlml,dnlZ,post] = gp(hyp, inffunc, meanfunc, covfuncF, likfunc, x, y, None, None, True) print "nlml =", nlml ## TRAINING: OPTIMIZE hyperparameters [hyp2_opt, fopt, gopt, funcCalls] = min_wrapper(hyp2,gp,'Minimize',inffunc,meanfunc,covfunc,likfunc,x,y,None,None,True) # minimize by Carl Rasmussen #[hyp2_opt, fopt, gopt, funcCalls] = min_wrapper(hyp2,gp,'CG',inffunc,meanfunc,covfunc,likfunc,x,y,None,None,True) # conjugent gradient #[hyp2_opt, fopt, gopt, funcCalls] = min_wrapper(hyp2,gp,'SCG',inffunc,meanfunc,covfunc,likfunc,x,y,None,None,True) # scaled conjugent gradient (faster than CG) #[hyp2_opt, fopt, gopt, funcCalls] = min_wrapper(hyp2,gp,'BFGS',inffunc,meanfunc,covfunc,likfunc,x,y,None,None,True) # quasi-Newton method of Broyden, Fletcher, Goldfarb, and Shanno (BFGS) print 'Optimal nlml =', fopt ## FITC PREDICTION [ymu,ys2,fmu,fs2,lp,post] = gp(hyp_opt, inffunc, meanfunc, covfuncF, likfunc, x, y, xstar, np.ones((n,1)) ) ## PLOT log predictive probabilities if PLOT: fig = plt.figure() plt.plot(x1[:,0], x1[:,1], 'b+', markersize = 12) plt.plot(x2[:,0], x2[:,1], 'r+', markersize = 12) plt.plot(u[:,0],u[:,1],'ko', markersize=12) pc = plt.contour(t1, t2, np.reshape(np.exp(lp), (t1.shape[0],t1.shape[1]) )) fig.colorbar(pc) plt.grid() plt.axis([-4, 4, -4, 4]) plt.show() pyGP_PR-1.0/src/Core/demoHousing.py000077500000000000000000000103701222232424600171470ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 30/09/2013 #================================================================================ import numpy as np from gp import gp #import Tools.general from Tools.min_wrapper import min_wrapper from UTIL.utils import convert_to_array, hyperParameters import matplotlib.pyplot as plt def HousingPlotter(y,z,ys): x = np.array(range(y.shape[0])) ym = np.reshape(y,(y.shape[0],)) zm = np.reshape(z,(z.shape[0],)) plt.plot(y, 'r-', linewidth = 3.0) plt.fill_between(x, ym + 1.*np.sqrt(zm), ym - 1.*np.sqrt(zm), facecolor=[0.,1.0,0.0,0.9],linewidths=0.0) plt.fill_between(x, ym + 2.*np.sqrt(zm), ym - 2.*np.sqrt(zm), facecolor=[0.,1.0,0.0,0.7],linewidths=0.0) plt.fill_between(x, ym + 3.*np.sqrt(zm), ym - 3.*np.sqrt(zm), facecolor=[0.,1.0,0.0,0.5],linewidths=0.0) plt.plot(ys, 'bx', linewidth = 3.0, markersize = 5.0) plt.grid() plt.show() if __name__ == '__main__': infile = '../../data/housing.txt' data = np.genfromtxt(infile) DN, DD = data.shape N = int(np.floor(DN/3.)) N = 50 # Get all data (exclude the 4th column which is binary) except the last 50 points for training x = np.concatenate((data[:-N,:4],data[:-N,5:-1]),axis=1) y = np.reshape(data[:-N,-1],(len(data[:-N,-1]),1)) # Test on the last 50 points xs = np.concatenate((data[-N:,:4],data[-N:,5:-1]),axis=1) ys = np.reshape(data[-N:,-1],(N,1)) N,D = x.shape ## DEFINE parameterized covariance function covfunc = [ ['kernels.covSum'], [ ['kernels.covSEiso'],['kernels.covNoise'] ] ] ## DEFINE parameterized mean function meanfunc = [ ['means.meanZero'] ] ## DEFINE parameterized inference and liklihood functions inffunc = ['inf.infExact'] likfunc = ['lik.likGauss'] ## SET (hyper)parameters hyp = hyperParameters() ## SET (hyper)parameters for covariance and mean hyp.cov = np.array([np.log(67.), np.log(66.), np.log(1.3)]) hyp.mean = np.array([]) sn = 0.1 hyp.lik = np.array([np.log(sn)]) print 'Initial mean = ',hyp.mean print 'Initial covariance = ',hyp.cov print 'Initial liklihood = ',hyp.lik [nlml, post] = gp(hyp,inffunc,meanfunc,covfunc,likfunc,x,y,None,None,False) print 'Initial negative log marginal likelihood = ',nlml ##----------------------------------------------------------## ## STANDARD GP (prediction) ## ##----------------------------------------------------------## vargout = gp(hyp,inffunc,meanfunc,covfunc,likfunc,x,y,xs) ym = vargout[0]; ys2 = vargout[1] m = vargout[2]; s2 = vargout[3] HousingPlotter(ym,ys2,ys) ##----------------------------------------------------------## ## STANDARD GP (training) ## ## OPTIMIZE HYPERPARAMETERS ## ##----------------------------------------------------------## ## -> parameter training using (off the shelf) conjugent gradient (CG) optimization (NOTE: SCG is faster) from time import clock t0 = clock() vargout = min_wrapper(hyp,gp,'SCG',inffunc,meanfunc,covfunc,likfunc,x,y,None,None,True) t1 = clock() hyp = vargout[0] vargout = gp(hyp,inffunc,meanfunc,covfunc,likfunc,x,y,xs) ym = vargout[0]; ys2 = vargout[1] m = vargout[2]; s2 = vargout[3] print 'Time to optimize = ',t1-t0 print 'Optimized mean = ',hyp.mean print 'Optimized covariance = ',hyp.cov print 'Optimized liklihood = ',hyp.lik [nlml, post] = gp(hyp,inffunc,meanfunc,covfunc,likfunc,x,y,None,None,False) print 'Final negative log marginal likelihood = ',nlml HousingPlotter(ym,ys2,ys) pyGP_PR-1.0/src/Core/demoMaunaLoa.py000066400000000000000000000105461222232424600172320ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 30/09/2013 #================================================================================ import numpy as np from gp import gp from Tools.min_wrapper import min_wrapper import Tools.general from UTIL.utils import convert_to_array, hyperParameters, plotter if __name__ == '__main__': # Example demo for pyGP prediction of carbon dioxide concentration using # the Mauna Loa CO2 data [Pieter Tans, Aug 2012]. # # The data is constantly updated and publically available under the link: # ftp://ftp.cmdl.noaa.gov/ccg/co2/trends/co2_mm_mlo.txt # # The used covariance function was proposed in [Gaussian Processes for # Machine Learning,Carl Edward Rasmussen and Christopher K. I. Williams, # The MIT Press, 2006. ISBN 0-262-18253-X]. # # Copyright (c) by Marion Neumann and Daniel Marthaler, 20/05/2013 ## LOAD data infile = '../../data/mauna.txt' # Note: Samples with value -99.99 were dropped. f = open(infile,'r') year = [] co2 = [] for line in f: z = line.split(' ') z1 = z[1].split('\n') if float(z1[0]) != -99.99: year.append(float(z[0])) co2.append(float(z1[0])) X = [i for (i,j) in zip(year,co2) if i < 2004] y = [j for (i,j) in zip(year,co2) if i < 2004] xx = [i for (i,j) in zip(year,co2) if i >= 2004] yy = [j for (i,j) in zip(year,co2) if i >= 2004] x = np.array(X) y = np.array(y) x = x.reshape((len(x),1)) y = y.reshape((len(y),1)) n,D = x.shape ## DEFINE parameterized covariance function covfunc = [ ['kernels.covSum'], [ ['kernels.covSEiso'],[['kernels.covProd'],[['kernels.covPeriodic'],['kernels.covSEiso']]],\ ['kernels.covRQiso'],['kernels.covSEiso'],['kernels.covNoise'] ] ] ## DEFINE parameterized mean function meanfunc = [ ['means.meanZero'] ] ## DEFINE parameterized inference and liklihood functions inffunc = ['inf.infExact'] likfunc = ['lik.likGauss'] ## SET (hyper)parameters hyp = hyperParameters() ## SET (hyper)parameters for covariance and mean hyp.cov = np.array([np.log(67.), np.log(66.), np.log(1.3), np.log(1.0), np.log(2.4), np.log(90.), np.log(2.4), \ np.log(1.2), np.log(0.66), np.log(0.78), np.log(1.6/12.), np.log(0.18), np.log(0.19)]) hyp.mean = np.array([]) sn = 0.1 hyp.lik = np.array([np.log(sn)]) ##----------------------------------------------------------## ## STANDARD GP (prediction) ## ##----------------------------------------------------------## xs = np.arange(2004+1./24.,2024-1./24.,1./12.) # TEST POINTS xs = xs.reshape(len(xs),1) vargout = gp(hyp,inffunc,meanfunc,covfunc,likfunc,x,y,xs) ym = vargout[0]; ys2 = vargout[1] m = vargout[2]; s2 = vargout[3] plotter(xs,ym,ys2,x,y) ##----------------------------------------------------------## ## STANDARD GP (training) ## ## OPTIMIZE HYPERPARAMETERS ## ##----------------------------------------------------------## ## -> parameter training using (off the shelf) conjugent gradient (CG) optimization (NOTE: SCG is faster) from time import clock t0 = clock() vargout = min_wrapper(hyp,gp,'SCG',inffunc,meanfunc,covfunc,likfunc,x,y,None,None,True) t1 = clock() hyp = vargout[0] vargout = gp(hyp,inffunc,meanfunc,covfunc,likfunc,x,y,xs) ym = vargout[0]; ys2 = vargout[1] m = vargout[2]; s2 = vargout[3] print 'Time to optimize = ',t1-t0 print 'Optimized mean = ',hyp.mean print 'Optimized covariance = ',hyp.cov print 'Optimized liklihood = ',hyp.lik plotter(xs,ym,ys2,x,y) pyGP_PR-1.0/src/Core/demoRegression.py000077500000000000000000000143231222232424600176550ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 30/09/2013 #================================================================================ import numpy as np import matplotlib.pyplot as plt from gp import gp from Tools.min_wrapper import min_wrapper import Tools.general from UTIL.utils import convert_to_array, hyperParameters, plotter, FITCplotter from time import clock if __name__ == '__main__': PLOT = True ## LOAD DATA demoData = np.load('../../data/regression_data.npz') x = demoData['x'] # training data y = demoData['y'] # training target xstar = demoData['xstar'] # test data ## PLOT data if PLOT: plt.plot(x,y,'r+',markersize=12) plt.axis([-1.9,1.9,-0.9,3.9]) plt.grid() plt.xlabel('input x') plt.ylabel('output y') plt.show() ## DEFINE parameterized mean and covariance functions covfunc = [['kernels.covPoly']] meanfunc = [ ['means.meanSum'], [ ['means.meanLinear'] , ['means.meanConst'] ] ] ## DEFINE likelihood function used likfunc = ['lik.likGauss'] ## SPECIFY inference method inffunc = ['inf.infExact'] ## SET (hyper)parameters hyp = hyperParameters() hyp.cov = np.array([np.log(0.25),np.log(1.0),1.0]) hyp.mean = np.array([0.5,1.0]) hyp.lik = np.array([np.log(0.1)]) ##----------------------------------------------------------## ## STANDARD GP (example 1) ## ##----------------------------------------------------------## print '...example 1: prediction...' ## PREDICTION t0 = clock() vargout = gp(hyp,inffunc,meanfunc,covfunc,likfunc,x,y,xstar) t1 = clock() ym = vargout[0]; ys2 = vargout[1]; m = vargout[2]; s2 = vargout[3] print 'Time for prediction =',t1-t0 ## PLOT results if PLOT: plotter(xstar,ym,s2,x,y,[-2, 2, -0.9, 3.9]) ## GET negative log marginal likelihood [nlml, post] = gp(hyp,inffunc,meanfunc,covfunc,likfunc,x,y,None,None,False) print "nlml =", nlml ##----------------------------------------------------------## ## STANDARD GP (example 2) ## ##----------------------------------------------------------## print '...example 2: prediction...' ## USE another covariance function -> for use of composite covariance functions see demoMaunaLoa.py covfunc = [ ['kernels.covSEiso'] ] ### SET (hyper)parameters hyp2 = hyperParameters() hyp2.cov = np.array([-1.0,0.0]) hyp2.mean = np.array([0.5,1.0]) hyp2.lik = np.array([np.log(0.1)]) ## PREDICTION t0 = clock() vargout = gp(hyp2,inffunc,meanfunc,covfunc,likfunc,x,y,xstar) t1 = clock() ym = vargout[0]; ys2 = vargout[1]; m = vargout[2]; s2 = vargout[3] print 'Time for prediction =',t1-t0 ## PLOT results if PLOT: plotter(xstar,ym,ys2,x,y,[-2, 2, -0.9, 3.9]) ## GET negative log marginal likelihood [nlml, post] = gp(hyp2,inffunc,meanfunc,covfunc,likfunc,x,y,None,None,False) print "nlml =", nlml ##----------------------------------------------------------## ## STANDARD GP (example 3) ## ##----------------------------------------------------------## print '...example 3: training and prediction...' ## TRAINING: OPTIMIZE HYPERPARAMETERS ## -> parameter training via off-the-shelf optimization t0 = clock() [hyp2_opt, fopt, gopt, funcCalls] = min_wrapper(hyp2,gp,'Minimize',inffunc,meanfunc,covfunc,likfunc,x,y,None,None,True) # minimize #[hyp2_opt, fopt, gopt, funcCalls] = min_wrapper(hyp2,gp,'CG',inffunc,meanfunc,covfunc,likfunc,x,y,None,None,True) # conjugent gradient #[hyp2_opt, fopt, gopt, funcCalls] = min_wrapper(hyp2,gp,'SCG',inffunc,meanfunc,covfunc,likfunc,x,y,None,None,True) # scaled conjugent gradient (faster than CG) #[hyp2_opt, fopt, gopt, funcCalls] = min_wrapper(hyp2,gp,'BFGS',inffunc,meanfunc,covfunc,likfunc,x,y,None,None,True) # quasi-Newton method of Broyden, Fletcher, Goldfarb, and Shanno (BFGS) t1 = clock() print 'Time for optimization =',t1-t0 print "Optimal nlml =", fopt ## PREDICTION vargout = gp(hyp2_opt,inffunc,meanfunc,covfunc,likfunc,x,y,xstar) ym = vargout[0]; ys2 = vargout[1]; m = vargout[2]; s2 = vargout[3] ## Plot results if PLOT: plotter(xstar,ym,ys2,x,y,[-1.9, 1.9, -0.9, 3.9]) ##----------------------------------------------------------## ## SPARSE GP (example 4) ## ##----------------------------------------------------------## print '...example 4: FITC training and prediction...' ## SPECIFY inducing points n = x.shape[0] num_u = np.fix(n/2) u = np.linspace(-1.3,1.3,num_u).T u = np.reshape(u,(num_u,1)) ## SPECIFY FITC covariance function covfunc = [['kernels.covFITC'], covfunc, u] ## SPECIFY FICT inference method inffunc = ['inf.infFITC'] ## TRAINING: OPTIMIZE hyperparameters [hyp2_opt, fopt, gopt, funcCalls] = min_wrapper(hyp2_opt,gp,'Minimize',inffunc,meanfunc,covfunc,likfunc,x,y,None,None,True) #[hyp2_opt, fopt, gopt, funcCalls] = min_wrapper(hyp2_opt,gp,'SCG',inffunc,meanfunc,covfunc,likfunc,x,y,None,None,True) print 'Optimal F =', fopt ## FITC PREDICTION vargout = gp(hyp2_opt, inffunc, meanfunc, covfunc, likfunc, x, y, xstar) ymF = vargout[0]; y2F = vargout[1]; mF = vargout[2]; s2F = vargout[3] ## Plot results if PLOT: FITCplotter(u,xstar,ymF,y2F,x,y,[-1.9, 1.9, -0.9, 3.9]) pyGP_PR-1.0/src/Core/gp.py000066400000000000000000000221051222232424600152700ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 30/09/2013 #================================================================================ import numpy as np import Tools.general from UTIL import solve_chol from UTIL.utils import checkParameters, unique def gp(hyp, inffunc, meanfunc, covfunc, likfunc, x, y, xs=None, ys=None, der=None): # Gaussian Process inference and prediction. The gp function provides a # flexible framework for Bayesian inference and prediction with Gaussian # processes for scalar targets, i.e. both regression and binary # classification. The prior is Gaussian process, defined through specification # of its mean and covariance function. The likelihood function is also # specified. Both the prior and the likelihood may have hyperparameters # associated with them. # # Two modes are possible: training or prediction: if no test cases are # supplied, then the negative log marginal likelihood and its partial # derivatives w.r.t. the hyperparameters is computed; this mode is used to fit # the hyperparameters. If test cases are given, then the test set predictive # probabilities are returned. Usage: # # training: [nlZ dnlZ post ] = gp(hyp, inf, mean, cov, lik, x, y, None, None, der) # prediction: [ymu ys2 fmu fs2 None post] = gp(hyp, inf, mean, cov, lik, x, y, xs, None, None, None) # or: [ymu ys2 fmu fs2 lp post] = gp(hyp, inf, mean, cov, lik, x, y, xs, ys, None) # # where: # # hyp column vector of hyperparameters # inffunc function specifying the inference method # covfunc prior covariance function (see below) # meanfunc prior mean function # likfunc likelihood function # x n by D matrix of training inputs # y column vector of length n of training targets # xs ns by D matrix of test inputs # ys column vector of length nn of test targets # der flag for dnlZ computation determination (when xs == None also) # # nlZ returned value of the negative log marginal likelihood # dnlZ column vector of partial derivatives of the negative # log marginal likelihood w.r.t. each hyperparameter # ymu column vector (of length ns) of predictive output means # ys2 column vector (of length ns) of predictive output variances # fmu column vector (of length ns) of predictive latent means # fs2 column vector (of length ns) of predictive latent variances # lp column vector (of length ns) of log predictive probabilities # # post struct representation of the (approximate) posterior # 3rd output in training mode and 6th output in prediction mode # # This is a python implementation of gpml functionality (Copyright (c) by # Carl Edward Rasmussen and Hannes Nickisch, 2011-02-18). # # Copyright (c) by Marion Neumann and Daniel Marthaler, 20/05/2013 if not inffunc: inffunc = ['inf.infExact'] # set default inf if not meanfunc: meanfunc = ['means.meanZero'] # set default mean if not covfunc: raise Exception('Covariance function cannot be empty') # no default covariance if covfunc[0] == 'kernels.covFITC': inffunc = ['inf.infFITC'] # only one possible inference alg for covFITC if not likfunc: likfunc = ['lik.likGauss'] # set default lik D = np.shape(x)[1] if not checkParameters(meanfunc,hyp.mean,D): raise Exception('Number of mean function hyperparameters disagree with mean function') if not checkParameters(covfunc,hyp.cov,D): raise Exception('Number of cov function hyperparameters disagree with cov function') if not checkParameters(likfunc,hyp.lik,D): raise Exception('Number of lik function hyperparameters disagree with lik function') # call the inference method # issue a warning if a classification likelihood is used in conjunction with # labels different from +1 and -1 if likfunc[0] == ['lik.likErf'] or likfunc[0] == ['lik.likLogistic']: uy = unique(y) ind = ( uy != 1 ) if any( uy[ind] != -1): raise Exception('You attempt classification using labels different from {+1,-1}\n') if not xs == None: # compute marginal likelihood and its derivatives only if needed vargout = Tools.general.feval(inffunc,hyp, meanfunc, covfunc, likfunc, x, y, 1) post = vargout[0] else: if not der: vargout = Tools.general.feval(inffunc, hyp, meanfunc, covfunc, likfunc, x, y, 2) post = vargout[0]; nlZ = vargout[1] else: vargout = Tools.general.feval(inffunc, hyp, meanfunc, covfunc, likfunc, x, y, 3) post = vargout[0]; nlZ = vargout[1]; dnlZ = vargout[2] if xs == None: # if no test cases are provided if not der: varargout = [nlZ, post] # report -log marg lik, derivatives and post else: varargout = [nlZ, dnlZ, post] # report -log marg lik, derivatives and post else: alpha = post.alpha L = post.L sW = post.sW nz = range(len(alpha[:,0])) if L == []: # in case L is not provided, we compute it K = Tools.general.feval(covfunc, hyp.cov, x[nz,:]) L = np.linalg.cholesky( (np.eye(nz) + np.dot(sW,sW.T)*K).T ) Ltril = np.all( np.tril(L,-1) == 0 ) # is L an upper triangular matrix? ns = xs.shape[0] # number of data points nperbatch = 1000 # number of data points per mini batch nact = 0 # number of already processed test data points ymu = np.zeros((ns,1)); ys2 = np.zeros((ns,1)) fmu = np.zeros((ns,1)); fs2 = np.zeros((ns,1)); lp = np.zeros((ns,1)) while nact<=ns-1: # process minibatches of test cases to save memory id = range(nact,min(nact+nperbatch,ns)) # data points to process kss = Tools.general.feval(covfunc, hyp.cov, xs[id,:], 'diag') # self-variances Ks = Tools.general.feval(covfunc, hyp.cov, x[nz,:], xs[id,:]) # cross-covariances ms = Tools.general.feval(meanfunc, hyp.mean, xs[id,:]) N = (alpha.shape)[1] # number of alphas (usually 1; more in case of sampling) Fmu = np.tile(ms,(1,N)) + np.dot(Ks.T,alpha[nz]) # conditional mean fs|f fmu[id] = np.reshape(Fmu.sum(axis=1)/N,(len(id),1)) # predictive means if Ltril: # L is triangular => use Cholesky parameters (alpha,sW,L) V = np.linalg.solve(L.T,np.tile(sW,(1,len(id)))*Ks) fs2[id] = kss - np.array([(V*V).sum(axis=0)]).T # predictive variances else: # L is not triangular => use alternative parametrization fs2[id] = kss + np.array([(Ks*np.dot(L,Ks)).sum(axis=0)]).T # predictive variances fs2[id] = np.maximum(fs2[id],0) # remove numerical noise i.e. negative variances Fs2 = np.tile(fs2[id],(1,N)) # we have multiple values in case of sampling if ys == None: [Lp, Ymu, Ys2] = Tools.general.feval(likfunc,hyp.lik,None,Fmu[:],Fs2[:],None,None,3) else: [Lp, Ymu, Ys2] = Tools.general.feval(likfunc,hyp.lik,np.tile(ys[id],(1,N)),Fmu[:],Fs2[:],None,None,3) lp[id] = np.reshape( np.reshape(Lp,(np.prod(Lp.shape),N)).sum(axis=1)/N , (len(id),1) ) # log probability; sample averaging ymu[id] = np.reshape( np.reshape(Ymu,(np.prod(Ymu.shape),N)).sum(axis=1)/N ,(len(id),1) ) # predictive mean ys|y and ... ys2[id] = np.reshape( np.reshape(Ys2,(np.prod(Ys2.shape),N)).sum(axis=1)/N , (len(id),1) ) # .. variance nact = id[-1]+1 # set counter to index of next data point if ys == None: varargout = [ymu, ys2, fmu, fs2, None, post] # assign output arguments else: varargout = [ymu, ys2, fmu, fs2, lp, post] # assign output arguments return varargout pyGP_PR-1.0/src/Tools/000077500000000000000000000000001222232424600145205ustar00rootroot00000000000000pyGP_PR-1.0/src/Tools/__init__.py000077500000000000000000000001051222232424600166300ustar00rootroot00000000000000 __all__ = ["general", "nearPD", "min_wrapper", "minimize", "scg"] pyGP_PR-1.0/src/Tools/general.py000077500000000000000000000157111222232424600165170ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 20/05/2013 #================================================================================ # Created on 08/09/2009 # # @author: Marion Neumann (last update 08/01/10) # # Substantial updates by Daniel Marthaler (July 2012). # get_nb_param() added by Marion Neumann (Aug 2012). import numpy as np from GPR.MEAN import means from GPR.COV import kernels from GPR.LIK import lik from GPR.INF import inf def feval(funcName, *args): assert(isinstance(funcName, list)) if len(funcName) > 1: if funcName[0] == ['kernels.covFITC']: if not args: return feval(funcName[0],funcName[1],funcName[2])[0] else: return feval(funcName[0],funcName[1],funcName[2], *args) # This is a composition assert(len(funcName) == 2) z = funcName[0] # Split off the module name (before the period) mod,fun = z[0].split('.') if mod == 'kernels': return getattr(kernels,fun)(funcName[1],*args) elif mod == 'means': return getattr(means,fun)(funcName[1],*args) elif mod == 'lik': return getattr(lik,fun)(funcName[1],*args) elif mod == 'inf': return getattr(inf,fun)(*args) else: raise Exception("Error in parameter function type") else: # This is either a singleton call or a call of one of the composed covariances z = funcName[0] if isinstance(z,list): # This is a singleton call mod,fun = z[0].split('.') else: # This is one called from a composition # Call the function, split off the module name (before the period) mod,fun = z.split('.') if mod == 'kernels': return getattr(kernels,fun)(*args) elif mod == 'means': return getattr(means,fun)(*args) elif mod == 'lik': return getattr(lik,fun)(*args) elif mod == 'inf': return getattr(inf,fun)(*args) else: raise Exception("Error in parameter function type") def convert_single(v,D): if isinstance(v,str): pram_str = v.split(' ') if pram_str[0]=='D': temp = int(D) if pram_str[1]=='+': temp += int(pram_str[2]) elif pram_str[1]=='-': temp -= int(pram_str[2]) else: raise Exception("ERROR: string representation is incorrect for this covariance") else: temp = v return temp def convert(v,D): if isinstance(v,list): w = [] for ii in range(len(v)): w.append(convert_single(v[ii],D)) else: w = convert_single(v,D) return w def check_hyperparameters(gp,x): check_hyperparameters_func(gp['covfunc'],gp['covtheta'],x) check_hyperparameters_func(gp['meanfunc'],gp['meantheta'],x,meanFlag=True) return True def check_hyperparameters_func(func,logtheta,x,meanFlag=False): [n,D] = x.shape if not meanFlag: if func[0][0] == 'kernels.covFITC': # This is an FITC approximation func = func[1:-1] xu = func[-1] if len(func) == 1: func = func[0] ## CHECK (hyper)parameters and mean/covariance function(s) if len(func) > 1: v = flatten(feval(func)) try: assert( sum(convert(v,D)) - len(logtheta) == 0 ) except AssertionError: if meanFlag: raise Exception('ERROR: number of hyperparameters does not match given mean function:' + str(sum(convert(v,D))) + ' hyperparameters needed (' + str(len(logtheta)) + ' given )!') else: raise Exception('ERROR: number of hyperparameters does not match given covariance function:'+ str(sum(convert(v,D))) + ' hyperparameters needed (' + str(len(logtheta)) + ' given )!') else: try: v = feval(func) assert( convert(v,D) - len(logtheta) == 0) except AssertionError: if meanFlag: raise Exception("ERROR: number of hyperparameters does not match given mean function: "+ str(convert(v,D)) + " hyperparameters needed ("+ str(len(logtheta)) + " given )!") else: raise Exception("ERROR: number of hyperparameters does not match given covariance function: "+ str(convert(v,D)) + " hyperparameters needed ("+ str(len(logtheta)) + " given )!") return True def get_nb_param(covList, dim): num_hyp = 0 for i in covList: if type(i) == int: num_hyp += i elif type(i) == str: add = True str2num = 0 str_list = i.split(' ') for j in range(0,len(str_list)): if str_list[j] == '+': add = True elif str_list[j] == '-': add = False else: try: current_number = int(str_list[j]) except ValueError: if str_list[j] == 'D': current_number = dim[1] elif str_list[j] == '+': add = True else: raise Exception( "ERROR: Not able to identify number of required hyperparameters!") return if add: str2num += current_number else: str2num -= current_number num_hyp += str2num return num_hyp def flatten(l, ltypes=(list, tuple)): ltype = type(l) l = list(l) i = 0 while i < len(l): while isinstance(l[i], ltypes): if not l[i]: l.pop(i) i -= 1 break else: l[i:i + 1] = l[i] i += 1 return ltype(l) def get_index_vec(a,b): ''' Returns vector of indices of b for the entries in a. Example: 1 2 2 a = 2 b = 3 index_vec = 0 3 1 1 Returns '-1' if entry in a is not an entry in b. ''' a = list(a) b = list(b) index_vec = [] for i in range(len(a)): try: index_vec.append( b.index(a[i]) ) except ValueError: index_vec.append(-1) return np.array(index_vec) pyGP_PR-1.0/src/Tools/min_wrapper.py000077500000000000000000000063311222232424600174230ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 20/05/2013 #================================================================================ # @author: Daniel Marthaler (Fall 2012) # update by Marion Neumann (May 20013) import numpy as np from copy import deepcopy from scipy.optimize import fmin_bfgs as bfgs from scipy.optimize import fmin_cg as cg from scg import scg from minimize import run from GPR.UTIL.utils import convert_to_array, convert_to_class def min_wrapper(hyp, F, Flag, *varargin): # Utilize scipy.optimize functions, sgc.py, or minimize.py to # minimize the negative log marginal liklihood. x = convert_to_array(hyp) # convert the hyperparameter class to an array if Flag == 'CG': aa = cg(nlml, x, dnlml, (F,hyp,varargin), maxiter=100, disp=False, full_output=True) x = aa[0]; fopt = aa[1]; funcCalls = aa[2]; gradcalls = aa[3] if aa[4] == 1: print "Maximum number of iterations exceeded." elif aa[4] == 2: print "Gradient and/or function calls not changing." gopt = dnlml(x,F,hyp,varargin) return convert_to_class(x,hyp), fopt, gopt, funcCalls elif Flag == 'BFGS': # Use BFGS aa = bfgs(nlml, x, dnlml, (F,hyp,varargin), maxiter=100, disp=False, full_output=True) x = aa[0]; fopt = aa[1]; gopt = aa[2]; Bopt = aa[3]; funcCalls = aa[4]; gradcalls = aa[5] if aa[6] == 1: print "Maximum number of iterations exceeded." elif aa[6] == 2: print "Gradient and/or function calls not changing." if isinstance(fopt, np.ndarray): fopt = fopt[0] return convert_to_class(x,hyp), fopt, gopt, funcCalls elif Flag == 'SCG': # use sgc.py aa = scg(x, nlml, dnlml, (F,hyp,varargin), niters = 100) hyp = convert_to_class(aa[0],hyp) fopt = aa[1][-1] gopt = dnlml(aa[0],F,hyp,varargin) return hyp, fopt, gopt, len(aa[1]) elif Flag == 'Minimize': # use minimize.py aa = run(x, nlml, dnlml, (F,hyp,varargin), maxnumfuneval=-100) hyp = convert_to_class(aa[0],hyp) fopt = aa[1][-1] gopt = dnlml(aa[0],F,hyp,varargin) return hyp, fopt, gopt, len(aa[1]) else: raise Exception('Incorrect usage of optimization flag in min_wrapper') def nlml(x,F,*varargin): hyp = varargin[0] temp = list(varargin[1:][0]) temp[-1] = False f = lambda z: F(z,*temp) X = convert_to_class(x,hyp) vargout = f(X) return vargout[0] def dnlml(x,F,*varargin): hyp = varargin[0] temp = list(varargin[1:][0]) temp[-1] = True f = lambda z: F(z,*temp) X = convert_to_class(x,hyp) vargout = f(X) z = convert_to_array(vargout[1]) return z pyGP_PR-1.0/src/Tools/minimize.py000077500000000000000000000212021222232424600167130ustar00rootroot00000000000000 #This program is distributed WITHOUT ANY WARRANTY; without even the implied #warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. # # #This file contains a Python version of Carl Rasmussen's Matlab-function #minimize.m # #minimize.m is copyright (C) 1999 - 2006, Carl Edward Rasmussen. #Python adaptation by Roland Memisevic 2008. # # #The following is the original copyright notice that comes with the #function minimize.m #(from http://www.kyb.tuebingen.mpg.de/bs/people/carl/code/minimize/Copyright): # # #"(C) Copyright 1999 - 2006, Carl Edward Rasmussen # #Permission is granted for anyone to copy, use, or modify these #programs and accompanying documents for purposes of research or #education, provided this copyright notice is retained, and note is #made of any changes that have been made. # #These programs and documents are distributed without any warranty, #express or implied. As the programs were written for research #purposes only, they have not been tested to the degree that would be #advisable in any important application. All use of these programs is #entirely at the user's own risk." """minimize.py This module contains a function 'minimize' that performs unconstrained gradient based optimization using nonlinear conjugate gradients. The function is a straightforward Python-translation of Carl Rasmussen's Matlab-function minimize.m """ from numpy import dot, isinf, isnan, any, sqrt, isreal, real, nan, inf def run(X, f, grad, args, maxnumlinesearch=None, maxnumfuneval=None, red=1.0, verbose=False): INT = 0.1;# don't reevaluate within 0.1 of the limit of the current bracket EXT = 3.0; # extrapolate maximum 3 times the current step-size MAX = 20; # max 20 function evaluations per line search RATIO = 10; # maximum allowed slope ratio SIG = 0.1;RHO = SIG/2;# SIG and RHO are the constants controlling the Wolfe- #Powell conditions. SIG is the maximum allowed absolute ratio between #previous and new slopes (derivatives in the search direction), thus setting #SIG to low (positive) values forces higher precision in the line-searches. #RHO is the minimum allowed fraction of the expected (from the slope at the #initial point in the linesearch). Constants must satisfy 0 < RHO < SIG < 1. #Tuning of SIG (depending on the nature of the function to be optimized) may #speed up the minimization; it is probably not worth playing much with RHO. SMALL = 10.**-16 #minimize.m uses matlab's realmin if maxnumlinesearch == None: if maxnumfuneval == None: raise "Specify maxnumlinesearch or maxnumfuneval" else: S = 'Function evaluation' length = maxnumfuneval else: if maxnumfuneval != None: raise "Specify either maxnumlinesearch or maxnumfuneval (not both)" else: S = 'Linesearch' length = maxnumlinesearch i = 0 # zero the run length counter ls_failed = 0 # no previous line search has failed f0 = f(X, *args) # get function value and gradient df0 = grad(X, *args) fX = [f0] i = i + (length<0) # count epochs?! s = -df0 d0 = -dot(s,s) # initial search direction (steepest) and slope x3 = red/(1.0-d0) # initial step is red/(|s|+1) while i < abs(length): # while not finished i = i + (length>0) # count iterations?! X0 = X; F0 = f0; dF0 = df0 # make a copy of current values if length>0: M = MAX else: M = min(MAX, -length-i) while 1: # keep extrapolating as long as necessary x2 = 0; f2 = f0; d2 = d0; f3 = f0; df3 = df0 success = 0 while (not success) and (M > 0): try: M = M - 1; i = i + (length<0) # count epochs?! f3 = f(X+x3*s, *args) df3 = grad(X+x3*s, *args) if isnan(f3) or isinf(f3) or any(isnan(df3)+isinf(df3)): print "error" return success = 1 except: # catch any error which occured in f x3 = (x2+x3)/2 # bisect and try again if f3 < F0: X0 = X+x3*s; F0 = f3; dF0 = df3 # keep best values d3 = dot(df3,s) # new slope if d3 > SIG*d0 or f3 > f0+x3*RHO*d0 or M == 0: # are we done extrapolating? break x1 = x2; f1 = f2; d1 = d2 # move point 2 to point 1 x2 = x3; f2 = f3; d2 = d3 # move point 3 to point 2 A = 6*(f1-f2)+3*(d2+d1)*(x2-x1) # make cubic extrapolation B = 3*(f2-f1)-(2*d1+d2)*(x2-x1) Z = B+sqrt(complex(B*B-A*d1*(x2-x1))) if Z != 0.0: x3 = x1-d1*(x2-x1)**2/Z # num. error possible, ok! else: x3 = inf if (not isreal(x3)) or isnan(x3) or isinf(x3) or (x3 < 0): # num prob | wrong sign? x3 = x2*EXT # extrapolate maximum amount elif x3 > x2*EXT: # new point beyond extrapolation limit? x3 = x2*EXT # extrapolate maximum amount elif x3 < x2+INT*(x2-x1): # new point too close to previous point? x3 = x2+INT*(x2-x1) x3 = real(x3) while (abs(d3) > -SIG*d0 or f3 > f0+x3*RHO*d0) and M > 0: # keep interpolating if (d3 > 0) or (f3 > f0+x3*RHO*d0): # choose subinterval x4 = x3; f4 = f3; d4 = d3 # move point 3 to point 4 else: x2 = x3; f2 = f3; d2 = d3 # move point 3 to point 2 if f4 > f0: x3 = x2-(0.5*d2*(x4-x2)**2)/(f4-f2-d2*(x4-x2)) # quadratic interpolation else: A = 6*(f2-f4)/(x4-x2)+3*(d4+d2) # cubic interpolation B = 3*(f4-f2)-(2*d2+d4)*(x4-x2) if A != 0: x3=x2+(sqrt(B*B-A*d2*(x4-x2)**2)-B)/A # num. error possible, ok! else: x3 = inf if isnan(x3) or isinf(x3): x3 = (x2+x4)/2 # if we had a numerical problem then bisect x3 = max(min(x3, x4-INT*(x4-x2)),x2+INT*(x4-x2)) # don't accept too close f3 = f(X+x3*s, *args) df3 = grad(X+x3*s, *args) if f3 < F0: X0 = X+x3*s; F0 = f3; dF0 = df3 # keep best values M = M - 1; i = i + (length<0) # count epochs?! d3 = dot(df3,s) # new slope if abs(d3) < -SIG*d0 and f3 < f0+x3*RHO*d0: # if line search succeeded X = X+x3*s; f0 = f3; fX.append(f0) # update variables #if verbose: print '%s %6i; Value %4.6e\r' % (S, i, f0) s = (dot(df3,df3)-dot(df0,df3))/dot(df0,df0)*s - df3 # Polack-Ribiere CG direction df0 = df3 # swap derivatives d3 = d0; d0 = dot(df0,s) if d0 > 0: # new slope must be negative s = -df0; d0 = -dot(s,s) # otherwise use steepest direction x3 = x3 * min(RATIO, d3/(d0-SMALL)) # slope ratio but max RATIO ls_failed = 0 # this line search did not fail else: X = X0; f0 = F0; df0 = dF0 # restore best point so far if ls_failed or (i>abs(length)):# line search failed twice in a row break # or we ran out of time, so we give up s = -df0; d0 = -dot(s,s) # try steepest x3 = 1/(1-d0) ls_failed = 1 # this line search failed if verbose: print "\n" return X, fX, i pyGP_PR-1.0/src/Tools/nearPD.py000077500000000000000000000032001222232424600162410ustar00rootroot00000000000000#================================================================================ # Marion Neumann [marion dot neumann at uni-bonn dot de] # Daniel Marthaler [marthaler at ge dot com] # Shan Huang [shan dot huang at iais dot fraunhofer dot de] # Kristian Kersting [kristian dot kersting at cs dot tu-dortmund dot de] # # This file is part of pyGP_FN. # The software package is released under the BSD 2-Clause (FreeBSD) License. # # Copyright (c) by # Marion Neumann, Daniel Marthaler, Shan Huang & Kristian Kersting, 20/05/2013 #================================================================================ import numpy as np def nearPD(A): tol = 1.e-8 # Returns the "closest" (up to tol) symmetric positive definite matrix to A. # Returns A if it is already Symmetric positive Definite count = 0; BOUND = 100 M = A.copy() while count < BOUND: M = (M+M.T)/2. eigval, Q = np.linalg.eig(M) eigval = eigval.real Q = Q.real xdiag = np.diag(np.maximum(eigval, tol)) M = np.dot(Q,np.dot(xdiag,Q.T)) try: L = np.linalg.cholesky(M) break except np.linalg.linalg.LinAlgError: count += 1 if count == BOUND: raise Exception("This matrix caused the nearPD algorithm to not converge") return M if __name__ == '__main__': A = np.array([[2.,-1,0,0.],[-1.,2.,-1,0],[0.,-1.,2.,-1.],[0.,0.,-1.,2.]]) A = np.random.random((4,4)) M = nearPD(A) try: L = np.linalg.cholesky(M) except np.linalg.linalg.LinAlgError: print "This shouldn't happen" print np.linalg.norm(M-A,'fro') pyGP_PR-1.0/src/Tools/scg.py000077500000000000000000000173211222232424600156550ustar00rootroot00000000000000from math import sqrt import numpy as np def scg(x, f, gradf, args, niters = 100, gradcheck = False, display = 0, flog = False, pointlog = False, scalelog = False, tolX = 1.0e-8, tolO = 1.0e-8, eval = None): #SCG Scaled conjugate gradient optimization. # # Description # [X, OPTIONS] = SCG(F, X, OPTIONS, GRADF) uses a scaled conjugate # gradients algorithm to find a local minimum of the function F(X) # whose gradient is given by GRADF(X). Here X is a row vector and F # returns a scalar value. The point at which F has a local minimum is # returned as X. The function value at that point is returned in # OPTIONS(8). # # [X, OPTIONS, FLOG, POINTLOG, SCALELOG] = SCG(F, X, OPTIONS, GRADF) # also returns (optionally) a log of the function values after each # cycle in FLOG, a log of the points visited in POINTLOG, and a log of # the scale values in the algorithm in SCALELOG. # # SCG(F, X, OPTIONS, GRADF, P1, P2, ...) allows additional arguments to # be passed to F() and GRADF(). The optional parameters have the # following interpretations. # # OPTIONS(1) = display is set to 1 to display error values; also logs error # values in the return argument ERRLOG, and the points visited in the # return argument POINTSLOG. If OPTIONS(1) is set to 0, then only # warning messages are displayed. If OPTIONS(1) is -1, then nothing is # displayed. # # OPTIONS(2)= tolX is a measure of the absolute precision required for the # value of X at the solution. If the absolute difference between the # values of X between two successive steps is less than OPTIONS(2), # then this condition is satisfied. # # OPTIONS(3)= tolO is a measure of the precision required of the objective # fuction at the solution. If the absolute difference between the # objective function values between two successive steps is less than # OPTIONS(3), then this condition is satisfied. Both this and the # previous condition must be satisfied for termination. # # OPTIONS(9) is set to 1 to check the user defined gradient function. # # OPTIONS(10) = funcCount returns the total number of function evaluations # (including those in any line searches). # # OPTIONS(11)= gradCount returns the total number of gradient evaluations. # # OPTIONS(14) = niters is the maximum number of iterations; default 100. # # See also # CONJGRAD, QUASINEW # # # Copyright (c) Ian T Nabney (1996-2001) if display: print '\n***** starting optimization (SCG) *****\n' nparams = len(x); # Check gradients if (gradcheck): pass #feval('gradchek', x, f, gradf, varargin{:}); eps = 1.0e-4 sigma0 = 1.0e-4; fold = f(x, *args); # Initial function value. fnow = fold; funcCount = 1; # Increment function evaluation counter. gradnew = gradf(x, *args) # Initial gradient. gradold = gradnew; gradCount = 1; # Increment gradient evaluation counter. d = -gradnew; # Initial search direction. success = 1; # Force calculation of directional derivs. nsuccess = 0; # nsuccess counts number of successes. beta = 1.0; # Initial scale parameter. betamin = 1.0e-15; # Lower bound on scale. betamax = 1.0e50; # Upper bound on scale. j = 1; # j counts number of iterations. if flog: pass #flog(j, :) = fold; if pointlog: pass #pointlog(j, :) = x; # Main optimization loop. listF = [fold] if eval is not None: evalue, timevalue = eval(x, *args) evalList = [evalue] time = [timevalue] while (j <= niters): # Calculate first and second directional derivatives. if (success == 1): mu = np.dot(d, gradnew) if (mu >= 0): d = - gradnew mu = np.dot(d, gradnew) kappa = np.dot(d, d) if (kappa < eps): print "FNEW: " , fnow #options(8) = fnow if eval is not None: return x, listF, evalList, time else: return x, listF sigma = sigma0/sqrt(kappa) xplus = x + sigma*d gplus = gradf(xplus, *args) gradCount += 1 theta = (np.dot(d, (gplus - gradnew)))/sigma; # Increase effective curvature and evaluate step size alpha. delta = theta + beta*kappa if (delta <= 0): delta = beta*kappa beta = beta - theta/kappa alpha = - mu/delta # Calculate the comparison ratio. xnew = x + alpha*d fnew = f(xnew, *args) funcCount += 1; Delta = 2*(fnew - fold)/(alpha*mu) if (Delta >= 0): success = 1; nsuccess += 1; x = xnew; fnow = fnew; listF.append(fnow) if eval is not None: evalue, timevalue = eval(x, *args) evalList.append(evalue) time.append(timevalue) else: success = 0; fnow = fold; if flog: # Store relevant variables #flog(j) = fnow; # Current function value pass if pointlog: #pointlog(j,:) = x; # Current position pass if scalelog: #scalelog(j) = beta; # Current scale parameter pass if display > 0: print('***** Cycle %4d Error %11.6f Scale %e' %( j, fnow, beta)) if (success == 1): # Test for termination # print type (alpha), type(d), type(tolX), type(fnew), type(fold) if ((max(abs(alpha*d)) < tolX) & (abs(fnew-fold) < tolO)): #options(8) = fnew; # print "FNEW: " , fnew if eval is not None: return x, listF, evalList, time else: return x, listF else: # Update variables for new position fold = fnew gradold = gradnew gradnew = gradf(x, *args) gradCount += 1 # If the gradient is zero then we are done. if (np.dot(gradnew, gradnew) == 0): # print "FNEW: " , fnew #options(8) = fnew; if eval is not None: return x, listF, evalList, time else: return x, listF # Adjust beta according to comparison ratio. if (Delta < 0.25): beta = min(4.0*beta, betamax); if (Delta > 0.75): beta = max(0.5*beta, betamin); # Update search direction using Polak-Ribiere formula, or re-start # in direction of negative gradient after nparams steps. if (nsuccess == nparams): d = -gradnew; nsuccess = 0; else: if (success == 1): gamma = np.dot((gradold - gradnew), gradnew)/(mu) d = gamma*d - gradnew; j += 1 # If we get here, then we haven't terminated in the given number of # iterations. #options(8) = fold; if (display): print "maximum number of iterations reached" if eval is not None: return x, listF, evalList, time else: return x, listF