diff --git a/examples/Data/randomGrid3D.xml b/examples/Data/randomGrid3D.xml new file mode 100644 index 000000000..6a82ce31c --- /dev/null +++ b/examples/Data/randomGrid3D.xml @@ -0,0 +1,3414 @@ + + + + + + + 32 + 1 + + + + + + 2 + 0 + 0 + 1 + + + + + + 9 + 7 + + 3.20776022311033415e+01 + -3.21030367555546334e+01 + -4.10921918858809718e+01 + -2.97297166857905353e+01 + 5.12273135695288744e+01 + 5.32024088136580744e+01 + 7.91786461660617107e+01 + 6.02804880595302208e+01 + -3.20503372222358784e+00 + -4.49785699465720157e+00 + 9.06604158034029126e+01 + 5.65784169718906416e+00 + 3.12765298180226239e+01 + 2.69796747523965372e+01 + 6.45574939874757661e+01 + -3.41086208590283135e+01 + 3.16243899688174857e+01 + -6.07548743284260979e+01 + 5.22759578856884062e+01 + 2.73875651690404389e+01 + -6.70630095670253041e+01 + 8.09592862312814248e+01 + -2.90562646005293850e+01 + 1.97217242898365228e+01 + 1.92063557124931243e+01 + -3.42017680808024593e+01 + 5.18983240742203904e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -0.00000000000000000e+00 + -8.81008786560076906e+00 + 9.93170954106413291e+01 + 7.64832733308121160e+00 + 9.74088247837876793e+01 + 6.98423466470551624e+00 + 2.15114230210294117e+01 + 8.81008786560076906e+00 + -9.93170954106413291e+01 + -7.64832733308121160e+00 + -0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 2.08303433971152465e+01 + 9.34532104006540187e+00 + -9.73589326595991906e+01 + -9.74088247837876793e+01 + -6.98423466470551624e+00 + -2.15114230210294117e+01 + -2.08303433971152465e+01 + -9.34532104006540187e+00 + 9.73589326595991906e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -9.96847734946016431e+01 + -9.41956956915240795e+00 + 3.58611803229205250e+01 + 1.37202268725982691e+02 + 8.31655615217660085e+01 + -2.96050013164699770e+01 + 5.56978428429618404e+01 + -1.64096733189184050e+02 + -6.76911341847846018e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 1 + 2 + + + + + + 9 + 7 + + -5.23677913323444741e+01 + 2.23054746283748031e+01 + -9.06323879014932388e+00 + 5.52805622243718773e+01 + 7.94539068269289146e+01 + 1.94541508392102571e+01 + -6.13821741431047201e+01 + 5.56974973159738553e+01 + -7.78665862090978944e+00 + -3.50785425178010168e+00 + -8.59643527732278869e+01 + 4.37949256925254957e+01 + 1.46448346762102268e+01 + -9.84525875940883211e+00 + -6.31187688571179351e+01 + -1.68743787040301001e+01 + 4.96429783126596220e+01 + 6.02933235638437850e+01 + -8.50324227920073668e+01 + -8.14408668091479448e+00 + 7.80614049299559465e+00 + -3.12168328896096590e+01 + -8.16862254520412101e+00 + 7.01147846710003364e+01 + 4.14413157080187986e+01 + -1.74987783157454646e+00 + 7.04591203724149580e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -0.00000000000000000e+00 + -2.60374020501451531e+00 + -5.42127786775135831e+01 + 8.39892562474970532e+01 + -9.97226081614345645e+01 + 7.26909685044996312e+00 + 1.60051631037066500e+00 + 2.60374020501451531e+00 + 5.42127786775135831e+01 + -8.39892562474970532e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -6.97294474564191535e+00 + -8.37146036187337330e+01 + -5.42516652512392241e+01 + 9.97226081614345645e+01 + -7.26909685044996312e+00 + -1.60051631037066500e+00 + 6.97294474564191535e+00 + 8.37146036187337330e+01 + 5.42516652512392241e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.08200476037924709e+01 + 1.28950357649989712e+02 + 1.43352339174516089e+02 + -2.38602247867744630e+01 + -5.20893463997723600e+01 + 2.84608516229487805e+01 + 6.76429765595175212e+01 + 1.20773294728176111e+02 + -1.20589752280510609e+02 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 2 + 3 + + + + + + 9 + 7 + + -1.96993691143317093e+01 + -2.72709220997507202e+01 + 4.46131383257576886e+01 + 9.30649268771392428e+01 + 3.22308546546604946e+00 + -1.69350730671129561e+01 + -3.00372292195983448e+01 + 4.73051200126789553e+01 + -6.91348627543176519e+01 + 8.77985890195871921e-01 + 2.95766095687473936e+01 + -7.96244285285575017e+01 + 7.31430977105361357e+00 + 9.52212779955403192e+01 + 2.32608409272044945e+01 + -1.06369613385613393e+00 + -2.30114188850419632e+00 + -5.58240668503478901e+01 + 8.63825257313862096e+01 + 3.81456041339674954e+01 + 2.72998827723357493e+01 + 3.09037953339841458e+01 + -1.77447464898715950e+01 + -1.04957062263738177e+01 + 3.97018822008379715e+01 + -7.27393475106917862e+01 + -4.57204758371371796e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -4.59114589009974594e+01 + 8.09444187944725400e+01 + 3.66079090854293696e+01 + -7.62990082722604797e+01 + -1.48208463505055157e+01 + -6.29190261377608593e+01 + 4.59114589009974594e+01 + -8.09444187944725400e+01 + -3.66079090854293696e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -4.55038380606685564e+01 + -5.68185144075392898e+01 + 6.85642555729200751e+01 + 7.62990082722604797e+01 + 1.48208463505055157e+01 + 6.29190261377608593e+01 + 4.55038380606685564e+01 + 5.68185144075392898e+01 + -6.85642555729200751e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -8.55922050845758786e-01 + 1.40043738778271234e+02 + -3.81613152068798769e+01 + -5.81227180496655009e+01 + -3.94713178138164338e+01 + 3.22753989340650662e+01 + 1.32631681079235051e+02 + -3.12863815406522541e+01 + -4.03072840384836510e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 3 + 4 + + + + + + 9 + 7 + + 8.66048172555662319e+01 + -1.52328612569721304e+01 + 4.48534445123967629e+01 + 2.02981055828854124e+01 + -4.93746572823205057e+01 + -2.74451016941264996e+01 + 4.12224960002507945e+00 + -6.40449795031895519e+01 + -5.03376341440287618e+01 + 4.92248012617162161e+01 + 3.93527186553859920e+01 + -6.89623533878909143e+01 + -4.08605085624496382e+01 + 5.95700970246627293e+01 + 3.55176671176527492e+01 + 9.26276585311693879e+00 + -6.84332455006064038e+01 + 4.88721436512163976e+00 + -4.03054194818879719e+00 + -1.78980733327891492e+01 + 3.45198923956961394e+01 + -8.80228615455582002e+01 + -4.55655125601066899e+01 + -9.64228954249413128e+00 + 1.13078887631363152e+01 + -3.24259783798813004e+01 + 8.58793556515557270e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -0.00000000000000000e+00 + -8.08548316492842929e+01 + 5.81668844249688917e+01 + 8.89436648998715995e+00 + 4.55868302729299728e+01 + 5.23625800953892124e+01 + 7.19722245829668879e+01 + 8.08548316492842929e+01 + -5.81668844249688917e+01 + -8.89436648998715995e+00 + -0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 3.72066809139564114e+01 + 6.22476807764455131e+01 + -6.88541148612105047e+01 + -4.55868302729299728e+01 + -5.23625800953892124e+01 + -7.19722245829668879e+01 + -3.72066809139564114e+01 + -6.22476807764455131e+01 + 6.88541148612105047e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -4.49706199136485765e+01 + 2.66285359051983370e+01 + 1.14028652024075640e+02 + 5.86398529753103119e+01 + 8.34746189227890056e+00 + 1.60810773504427289e+02 + 1.79694799386287997e+02 + -4.53728709677613509e+01 + -1.70780906392258842e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 4 + 5 + + + + + + 9 + 7 + + 6.43003424386522084e+01 + 1.19796180071396279e+01 + 4.55761542738693137e+01 + -6.52039759607597773e+01 + 5.89636417155482349e+01 + 1.80719252805888360e+01 + -1.51516354244859457e+01 + -5.00500232613797280e+01 + -5.34352949392747476e+01 + 7.42583468150113362e+01 + 7.45315542640838657e+00 + -5.24573072801267415e+01 + 4.09274245687232323e+01 + -1.68412174720425938e+01 + -1.36705786744805522e+01 + 2.70644906205976383e+01 + -8.31804320954470597e+01 + 4.33744783589924054e+01 + -1.22788002717957223e+01 + -6.11064170046755617e+00 + -7.16592171005377452e+01 + -5.30800076837930419e+01 + -7.82610172559931669e+01 + 2.92797087670289500e+01 + -2.22315273959740871e+01 + -2.13187310855278938e+01 + -6.26781275987031279e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 6.86610695050494400e+01 + 1.26119199040840346e+01 + 7.16002584545316836e+01 + -7.26277336162735736e+01 + 1.63646703530006263e+01 + 6.67638365734117087e+01 + -6.86610695050494400e+01 + -1.26119199040840346e+01 + -7.16002584545316836e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -3.29694467444833483e+00 + -9.78424092128285139e+01 + 2.03959092771819215e+01 + 7.26277336162735736e+01 + -1.63646703530006263e+01 + -6.67638365734117087e+01 + 3.29694467444833483e+00 + 9.78424092128285139e+01 + -2.03959092771819215e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.74531092467603024e+01 + -9.68528064078233442e-01 + -1.43955992851650620e+01 + -3.71890171910416072e+01 + 2.70827825121624883e+01 + 1.59657441575989878e+02 + -1.61093211332558155e+02 + -2.36356649401689936e+01 + -3.50656643781457547e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 5 + 6 + + + + + + 9 + 7 + + -6.40637543398504619e+01 + -1.35919926214779476e+01 + -7.55586787939104454e+01 + 7.49591051894634859e+01 + -1.51203828101076443e+01 + -6.04179319959996945e+01 + 1.63126149967003826e+01 + -3.67352169146467711e+00 + -1.49855747686980401e+01 + -7.64033153697704415e-01 + 5.83204655642280656e+00 + -2.26064773215979020e+00 + -1.21801830224353154e+01 + 9.23282518419687364e+01 + -3.58808093806811215e+01 + 6.76481543691339766e+01 + 3.42590164955414807e+01 + 6.51923982648987277e+01 + -2.26801968652302843e+00 + -9.78783400225560598e+01 + 1.94155506146273105e+01 + 1.64832141235940597e+01 + 2.49817026599325320e+00 + -1.62466150125265507e+01 + -7.07730069025061255e+01 + 1.59540577889256383e+01 + 6.50575730811547004e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -0.00000000000000000e+00 + -2.60428350122129437e+01 + -6.97694761102178944e+01 + 6.67382270354284231e+01 + -9.26426416957655618e+01 + 3.75218214997053465e+01 + 3.07471152699694139e+00 + 2.60428350122129437e+01 + 6.97694761102178944e+01 + -6.67382270354284231e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -2.71866085445878198e+01 + -6.10273144964611944e+01 + -7.44080318325472234e+01 + 9.26426416957655618e+01 + -3.75218214997053465e+01 + -3.07471152699694139e+00 + 2.71866085445878198e+01 + 6.10273144964611944e+01 + 7.44080318325472234e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.03933661487596595e+02 + 4.68471191704998802e+01 + 1.18875204735324314e+01 + -2.97055153316788463e+01 + 7.27443939906840598e+01 + 7.23444094058463492e+01 + 3.82415703375670404e+01 + -2.27415768204949451e+01 + -3.07813391103898937e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 6 + 7 + + + + + + 9 + 7 + + 1.85630374271583953e+01 + -5.29313735130414926e+01 + 3.66303674355102942e+01 + 1.04645716047590440e+01 + 5.73485858047936077e+01 + -5.08591115829850366e+01 + 9.38445241664948213e+01 + -1.36042785531324792e+01 + -2.31303222840427054e+01 + -9.56612084104992455e+00 + -2.98701442261074703e+01 + 7.62840716168778670e+01 + -1.92910486636706899e+01 + -5.04699434611332478e+01 + 3.46538201809487560e+01 + 3.08818026153629788e+01 + 7.17821170583860777e+01 + 5.45013702935509130e+01 + -5.34080385938858839e+01 + 6.23631140795577465e+01 + 4.43545654716690905e+01 + -7.95485560661077074e+01 + -2.30201810552558968e+01 + -5.51074857191607634e+01 + 1.18903119380668851e+01 + -2.67045605220385767e+01 + -2.32776001751549053e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.79872266803328138e+01 + 8.64209977399809475e+01 + -4.69879859748907549e+01 + 9.82455450157955568e+01 + -1.33902787243874393e+01 + 1.29812680517931049e+01 + -1.79872266803328138e+01 + -8.64209977399809475e+01 + 4.69879859748907549e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 4.92671908064704667e+00 + -4.84985730234346519e+01 + -8.73133200250165942e+01 + -9.82455450157955568e+01 + 1.33902787243874393e+01 + -1.29812680517931049e+01 + -4.92671908064704667e+00 + 4.84985730234346519e+01 + 8.73133200250165942e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -8.68503697620646733e+01 + -6.63265206223623283e-01 + 5.77899517416623070e+01 + 1.54728954524297080e+02 + -7.36729123490597857e+01 + -4.33714485898333280e+01 + -8.08065933920789270e+00 + -1.49263004123200204e+02 + 1.24142002075972883e+02 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 7 + 8 + + + + + + 9 + 7 + + 4.74062871911341759e+01 + 5.71222009314007693e+01 + -2.90538462259461063e+01 + 4.09276699085499303e+01 + 3.99974137866047528e+01 + -1.99073676332006535e+01 + -7.78029399881401531e+01 + 5.27745236895511312e+01 + -3.37039333639074385e+01 + -9.12848490867447815e+00 + -6.92290536396862421e+01 + -2.32360574632871630e+00 + 1.58993825041621264e+01 + 5.31010116114126660e+01 + -6.55564733827918502e+01 + -3.41955054720470475e+00 + 4.69927897084473187e+01 + 7.43534572740371971e+01 + -6.28440645552295862e+01 + 4.27829075558255525e+01 + 5.46299457034644860e+01 + 7.56316908352858093e+01 + 1.42185369226067397e+01 + 5.51031175481571864e+01 + 2.63284733063378651e+00 + 2.23353419791525312e+01 + 4.20873574259948739e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.80585028832817862e+01 + -9.52829540657210430e+01 + -2.43936290478516966e+01 + -9.74947022342857252e+01 + -2.06160114466111395e+01 + 8.35243127964282905e+00 + -1.80585028832817862e+01 + 9.52829540657210430e+01 + 2.43936290478516966e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -1.29874366163019737e+01 + 2.22741719608807998e+01 + -9.66187753679099757e+01 + 9.74947022342857252e+01 + 2.06160114466111395e+01 + -8.35243127964282905e+00 + 1.29874366163019737e+01 + -2.22741719608807998e+01 + 9.66187753679099757e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 7.39787443001634273e+01 + -3.29520218737291373e+01 + 1.75142945360573606e+02 + -4.89974435425370416e+01 + -9.39530685963477481e+01 + -3.92888663032043226e+01 + 4.46158732311312889e+01 + 1.62422063134678211e+02 + -1.51571321993876627e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 8 + 9 + + + + + + 9 + 7 + + 9.41130157063151671e+01 + 1.30096002577236085e+01 + -9.65142385029295902e+00 + -2.97023439769390656e+01 + 3.10157410489098062e-01 + 4.06407354840771795e+00 + -9.58311576289570155e+00 + 9.66152903465384156e+01 + 2.35615250878309794e+01 + 1.77277266667175297e+01 + -4.94856229862016690e+01 + 8.47852845508831336e+01 + 1.34310978254269742e+01 + -8.43883523083151488e+01 + -5.15167787166412623e+01 + 3.41073399110957187e-01 + 2.18139939555895035e+00 + 9.36548913434189423e+00 + -2.80240621096084510e+01 + -7.66577929709818839e+00 + 9.19688865464188510e+00 + -9.38736839860108603e+01 + -1.81862826844923511e+01 + 1.54239226370280336e+00 + -8.73344590429244683e+00 + 2.30576115914954443e+01 + -9.65297250180462356e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -9.99929178544290664e+01 + 6.62091181731783163e-01 + -9.88946016968861308e-01 + -6.23564336177379053e-01 + 4.16306145486879302e+01 + 9.09203118110426516e+01 + 9.99929178544290664e+01 + -6.62091181731783163e-01 + 9.88946016968861308e-01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.01367967132287018e+00 + 9.09200394168724984e+01 + -4.16235376434812636e+01 + 6.23564336177379053e-01 + -4.16306145486879302e+01 + -9.09203118110426516e+01 + -1.01367967132287018e+00 + -9.09200394168724984e+01 + 4.16235376434812636e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 5.55282235936336566e+00 + -1.76495493086246057e+02 + -9.69626169859593645e+00 + -1.18086675146736617e+01 + 9.21060305839002780e+01 + 5.32168712759193419e+00 + 1.99148216329777483e+02 + 1.06989143861811957e+01 + -5.26829527913709494e+00 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 9 + 10 + + + + + + 9 + 7 + + -8.63087793249897430e+01 + -1.92823330270961684e+01 + -4.42211315553459627e+01 + 2.21505354932965801e+01 + -2.63320996942306031e+01 + -5.69788664986630025e+01 + -4.53775288609686314e+01 + 2.58526272351369890e+01 + 5.53658756361871127e+01 + -1.54315640909267664e+01 + -4.32654755327350813e+01 + 1.97045712499729326e+01 + 7.39853470327932428e+01 + -5.76879953686396973e+01 + 2.69216116215563304e+01 + 6.54801163688819514e+01 + 5.53812651216026666e+01 + -2.49089008945160231e+01 + -9.71189061792836705e-01 + -7.25245368327360609e+01 + 4.96429238183381116e+01 + 6.85930292060070612e-01 + -1.12841816293779988e+01 + -7.69094846357958062e+01 + -5.16757922033419734e-02 + -6.79147932987382035e+01 + -4.02434226181160639e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 7.83449748282976799e+01 + 6.13201754047114846e+01 + 1.00946028891096997e+01 + -1.98027611608243852e+01 + 9.23644092983748521e+00 + 9.75834966034573341e+01 + -7.83449748282976799e+01 + -6.13201754047114846e+01 + -1.00946028891096997e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 5.89059892503364324e+01 + -7.84507759508154265e+01 + 1.93793752003053505e+01 + 1.98027611608243852e+01 + -9.23644092983748521e+00 + -9.75834966034573341e+01 + -5.89059892503364324e+01 + 7.84507759508154265e+01 + -1.93793752003053505e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -1.06986312175636783e+02 + 1.28415573927117009e+02 + 5.26747761738054123e+01 + -8.33248123208498015e+01 + -6.72578872464476092e+01 + 1.08240234163658471e+02 + -1.79030185090507175e+01 + -1.01966657528716965e+02 + 5.84227812301104876e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 10 + 11 + + + + + + 9 + 7 + + 1.17302291731448682e+00 + 2.23736436785092039e+01 + 8.70065012231049764e+01 + 7.67707833816357947e+01 + 5.73162887867323292e+01 + -1.32881532104113664e+00 + -2.51982778355239070e+01 + -7.72901736529863381e+00 + 4.53049731782644329e+01 + 8.54937155952362531e+01 + -4.11223504877864485e+01 + 2.12416634370425328e+01 + -2.90259802009985357e+01 + -8.76138896821732693e+00 + -5.06100435935214055e+00 + -3.81597643193563130e+01 + -9.02600433698408864e+01 + 7.99447271884073940e-01 + -2.24916749275050378e+01 + 2.24299023473424732e+01 + 3.83080111263447165e+01 + -5.03856743318980023e+01 + 7.20086537445330919e+01 + -4.63731853719149143e+01 + 2.85914723565250490e+01 + -2.34837711903294171e+01 + -7.92472197577816644e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -2.63903093070078647e+01 + -2.67183596992454575e+01 + 9.26805310173727577e+01 + -8.93547271763020774e+01 + 4.29554589851731663e+01 + -1.30599109729635039e+01 + 2.63903093070078647e+01 + 2.67183596992454575e+01 + -9.26805310173727577e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -3.63219534982506076e+01 + -8.62609865371063478e+01 + -3.52101959056759028e+01 + 8.93547271763020774e+01 + -4.29554589851731663e+01 + 1.30599109729635039e+01 + 3.63219534982506076e+01 + 8.62609865371063478e+01 + 3.52101959056759028e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -1.04067559180773088e+01 + 7.90108728350369383e-01 + 5.78189525430708571e+01 + -1.16280923178694167e+02 + 8.10714898076704600e+01 + 7.53829806570444134e+01 + -5.78201683295110769e+01 + 6.19589188064910701e+01 + -1.33505720595716895e+02 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 11 + 12 + + + + + + 9 + 7 + + 3.68344254692879192e+00 + 4.86691203808377892e+00 + -1.57231292941544805e+01 + 3.72290369400325147e+01 + 2.06331817836064921e+01 + -8.89534944072540839e+01 + -8.53296212948073389e+01 + 4.59306709423131920e+01 + -2.45010409318042690e+01 + 2.57025147685637876e+01 + 2.56000347344376031e+01 + -9.16529948262070491e+01 + 2.63288801817711970e-01 + 6.59382843634221683e+00 + 1.98864164108950128e+01 + -3.68109153857383404e+01 + -8.63480489179755466e+01 + -3.21511205713796500e+01 + 8.19562719339757138e+01 + -5.68704775649636858e+01 + 6.94175446640136773e+00 + -5.03482257069523769e+01 + -7.72158009035328661e+01 + -3.65965187643956114e+01 + -7.19547325119217618e+00 + -9.46558534641098959e+00 + -4.77765582369805042e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -7.29991262018860425e+01 + 6.43134564832916169e+01 + -2.31280541536218962e+01 + -5.24720785637083011e+01 + -7.44213220616246360e+01 + -4.13297446617319864e+01 + 7.29991262018860425e+01 + -6.43134564832916169e+01 + 2.31280541536218962e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -4.37927910159324512e+01 + -1.80345817187895037e+01 + 8.80735222258265793e+01 + 5.24720785637083011e+01 + 7.44213220616246360e+01 + 4.13297446617319864e+01 + 4.37927910159324512e+01 + 1.80345817187895037e+01 + -8.80735222258265793e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 9.48709591330575535e+01 + 9.60518679589807789e+01 + -5.19581436253384297e+01 + 2.54957713470891250e+01 + -1.34158356561916207e+02 + -2.25547332522903439e+01 + 1.09221138500000791e+02 + -4.57450052559468858e+01 + -6.82129290695598485e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 12 + 13 + + + + + + 9 + 7 + + 1.93941690967042213e+01 + -3.90998618366849584e+01 + 1.63698620253740224e+00 + 1.66123034210985914e+01 + 8.85555525420248557e+01 + 2.63934202043697148e+01 + -8.61979615512473032e+01 + 1.74318161103396001e+01 + -3.92907393347052789e+01 + -7.36585675537442341e+01 + 5.45643959900256945e+01 + -4.73596788713465688e+00 + 1.23230795012408407e+01 + 3.75671217770436385e+01 + -8.32914041112130690e+01 + -4.08563087068878872e+01 + 8.53843323058891990e+00 + 3.64960868565984242e+01 + -6.45808132649689668e+01 + -7.41043072657268738e+01 + -2.22937919750649804e+00 + -1.10797951271104100e+00 + -1.86389067743133623e+01 + -4.83331185591061043e+01 + 2.23431972515867372e+01 + -4.99220528272985842e+00 + -8.44039814228390526e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -0.00000000000000000e+00 + 5.32207681440972067e+01 + -3.06680603945050478e+01 + 7.89114688102511224e+01 + -5.56684880542022569e+01 + -8.29014582098623265e+01 + 5.32613034365356164e+00 + -5.32207681440972067e+01 + 3.06680603945050478e+01 + -7.89114688102511224e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 6.37853374680371630e+01 + -4.67634290692784660e+01 + -6.11932383991170568e+01 + 5.56684880542022569e+01 + 8.29014582098623265e+01 + -5.32613034365356164e+00 + -6.37853374680371630e+01 + 4.67634290692784660e+01 + 6.11932383991170568e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -6.90389218165145877e+01 + 4.83379444724982577e+01 + 1.61042728694415416e+02 + -1.53499694264126390e+02 + -6.29210660700233788e+01 + -1.36328982272427429e-01 + -7.32570951739669738e+01 + -6.73028752573054874e+01 + -7.84207879822722163e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 13 + 14 + + + + + + 9 + 7 + + 6.06804300720567937e+01 + 6.78223389316968621e+01 + 1.14214673416742691e+01 + 2.85045083366314387e+01 + 2.67345752695483867e+01 + 9.28156907783291274e+00 + -3.79052729172203868e+01 + 4.99857449838935750e+01 + -7.77096553357260120e+01 + -7.82277784247720263e+01 + 5.40246939742984793e+01 + 2.32633737948870305e+01 + 2.47956201341556586e+01 + -1.05317240685613065e+01 + 9.52360095902563302e+01 + -1.29078256939381077e+01 + 9.88755369559818398e+00 + 1.89881719971782701e+01 + 9.09983932556493968e+00 + 4.26215664457260388e+01 + 1.04271324229186870e+01 + -6.00988012459201499e+01 + -6.92347566083080466e+01 + 1.36255057217382767e+01 + 5.89311254734215737e+01 + -4.94916274785966195e+01 + -5.89796648826666967e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -2.89030917873941853e+01 + 8.70749520313458305e+01 + -3.97814531392241904e+01 + -9.56630827209552450e+01 + -2.78467078053438328e+01 + 8.55192777842975183e+00 + 2.89030917873941853e+01 + -8.70749520313458305e+01 + 3.97814531392241904e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -3.63123800557649812e+00 + 4.05279359595652906e+01 + 9.13471429077244466e+01 + 9.56630827209552450e+01 + 2.78467078053438328e+01 + -8.55192777842975183e+00 + 3.63123800557649812e+00 + -4.05279359595652906e+01 + -9.13471429077244466e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.43780450624784670e+01 + -6.63129258801445189e+01 + 4.67182022743824987e+00 + -2.52301003257763057e+01 + -9.40350333946463337e+01 + -1.71054540098886392e+01 + 9.90721820600476235e+01 + -1.66885836907954435e+01 + 5.08294565998758543e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 14 + 15 + + + + + + 9 + 7 + + -2.76355843137293995e+01 + -7.18186699571359952e+01 + 3.07651775764515065e+01 + -6.40640399324160086e+01 + -2.68919032866484322e+01 + 9.38443754617159698e+00 + -7.15464009647350991e+01 + 4.97633461308791567e+01 + -2.49225935590679519e+01 + 1.43932895637872207e+01 + -2.88438285303894943e+01 + 6.70843814351506182e+01 + -2.31361564910724091e+01 + 9.02877315480068603e+01 + 3.51897296575981642e+01 + 1.03017517985169054e+01 + -1.55151545070516619e+01 + 6.47165615401468699e+01 + 5.11107200952757168e+01 + 4.26513944039867638e+01 + 5.62224381230135037e+01 + -6.85270728206629229e+01 + -1.63576632701424352e+01 + 1.39839419419296167e+01 + 4.30853333607984084e+01 + -1.81999495152501289e+01 + -7.11401633803200042e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -8.15059339771809164e+01 + -5.54598536753216109e+01 + -1.67626774955357369e+01 + -1.04808841952049807e+01 + 4.25686446848321935e+01 + -8.98780371179923492e+01 + 8.15059339771809164e+01 + 5.54598536753216109e+01 + 1.67626774955357369e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 5.69818724946287603e+01 + -7.14990567770542356e+01 + -4.05086544703406517e+01 + 1.04808841952049807e+01 + -4.25686446848321935e+01 + 8.98780371179923492e+01 + -5.69818724946287603e+01 + 7.14990567770542356e+01 + 4.05086544703406517e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -0.00000000000000000e+00 + 2.31230716791378725e+01 + 2.46912491883496621e+01 + 3.19609408881894597e+00 + 1.52981265769411010e+01 + 1.32876273309615431e+01 + 2.19836964146123970e+00 + 2.74814898754146313e+01 + -2.79167095337952915e+01 + 5.37372498127196874e+00 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 15 + 16 + + + + + + 9 + 7 + + -5.10421334268497020e+01 + -3.97471595810313332e+01 + -1.64389363316277826e+00 + 2.54089618298921245e+01 + -3.15327493370988705e+01 + 9.13980507822707153e+01 + -5.91789668556331776e+01 + -4.81012757100621684e+01 + 1.65563384387473489e+00 + -6.78890977332195860e+01 + -3.50650650222022904e+01 + -8.69105008709039062e+00 + -5.80472933502369273e+01 + 6.88255140605843536e+01 + 3.93666342323419940e+01 + 4.37233513120340049e+01 + 4.71317258377345070e+01 + 1.72241668652535651e+01 + -4.44141172746378885e+01 + 7.81566561420669927e+01 + -4.26269969902096406e+01 + 1.02994128500864690e+01 + -1.21754201166700025e+01 + -9.79654030213998794e+00 + -3.57503904294827279e+01 + 2.58360247551576450e+01 + 8.83514468333149239e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -8.11187605223852870e+01 + -5.83599337897738124e+01 + 3.72354929672284474e+00 + 4.13366904132922599e+01 + -5.27200270517053724e+01 + 7.42420148793328991e+01 + 8.11187605223852870e+01 + 5.83599337897738124e+01 + -3.72354929672284474e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -4.13645345312568224e+01 + 6.17633943021323262e+01 + 6.68898976474907556e+01 + -4.13366904132922599e+01 + 5.27200270517053724e+01 + -7.42420148793328991e+01 + 4.13645345312568224e+01 + -6.17633943021323262e+01 + -6.68898976474907556e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 6.98804353612945732e+01 + -9.46593060564561881e+01 + -1.56916014852066439e+02 + -3.53367495425767473e+01 + -1.16911080827186709e+02 + 7.34112761651981316e+01 + 2.36033204137547230e+01 + 1.27622949564988858e+02 + -4.72506285651327431e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 16 + 17 + + + + + + 9 + 7 + + -5.85435867217505290e+01 + 2.33746012476837990e+01 + -5.77863483723881544e+01 + -5.63329646716383934e+01 + -7.47638541726440593e+01 + 3.41977282497739097e+01 + -4.10879444630993547e+01 + 7.01966859196727366e+00 + -3.18903501380708931e+01 + 7.49027535040175820e+01 + -1.88037812404022127e+01 + -5.50573390958028739e+01 + -1.05105126900153234e+01 + 3.04370847103721403e+01 + 2.74657650627104459e+01 + -5.06989372188159564e+01 + -7.70752893724180126e+01 + -2.65036593071036783e+01 + -8.31624785461406368e+00 + 7.75118867469858053e+00 + -5.96862946091213260e+01 + 3.40815941189988152e+01 + -5.16283223754141716e+01 + -6.67021661737265248e+01 + 4.67870803895238083e+01 + -6.21976850537821235e+01 + 4.44353424269327562e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -9.98522291990096846e+01 + -2.17344032867269910e+00 + 4.98081129196222250e+00 + -1.69510675143138556e+00 + -7.46243396331473150e+01 + -6.65457327513788073e+01 + 9.98522291990096846e+01 + 2.17344032867269910e+00 + -4.98081129196222250e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 5.16322932762926534e+00 + -6.65318276575533361e+01 + 7.44772245149880234e+01 + 1.69510675143138556e+00 + 7.46243396331473150e+01 + 6.65457327513788073e+01 + -5.16322932762926534e+00 + 6.65318276575533361e+01 + -7.44772245149880234e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -3.50474015754750923e+01 + -2.85468905697825797e+01 + -8.26608268609589487e+01 + 7.28338697761986396e+01 + -1.03233140690458100e+02 + -6.32144076335772187e+00 + 4.02504937228882511e+01 + 1.40718088483059116e+01 + 7.44299004116042227e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 17 + 18 + + + + + + 9 + 7 + + -8.82393493416751795e+01 + -6.12481728289156635e+00 + 6.45898243198292366e-01 + 3.86841248560220095e+01 + -4.52281770219861770e+01 + -4.30849606407744687e+01 + -2.62811256127213326e+01 + -5.90797411700203270e+01 + -5.09550094510600857e+01 + 4.68685179819303457e+01 + -2.73625956906887646e+00 + 1.87456681275038339e+00 + 7.09242852066529537e+01 + -3.96201491012311280e+01 + 4.26233675831482230e+01 + -5.26328438238840803e+01 + -5.34225198350244597e+01 + 6.12486839444391720e+01 + -3.76709521788862034e+00 + 9.35400534400023815e+01 + -3.47030072986988287e+01 + 3.20432748703922687e+00 + 2.41487264929615755e+01 + 7.47759826419143678e+01 + -2.25246616168526836e+00 + -2.57366654666450962e+01 + -5.63865719112063317e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 9.03046090553630876e+01 + 2.33581936525576914e+01 + -3.60481951371722431e+01 + 1.31196185024049097e+01 + -9.49105504326207523e+01 + -2.86332503730943557e+01 + -9.03046090553630876e+01 + -2.33581936525576914e+01 + 3.60481951371722431e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -4.09017504968844747e+01 + 2.11277591302666750e+01 + -8.87731074167399186e+01 + -1.31196185024049097e+01 + 9.49105504326207523e+01 + 2.86332503730943557e+01 + 4.09017504968844747e+01 + -2.11277591302666750e+01 + 8.87731074167399186e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -0.00000000000000000e+00 + 4.26251548517658989e+01 + -5.58446106332349572e+01 + -4.99134201995332649e+00 + 7.19670254876823492e+01 + -1.87430011171921045e+01 + -5.57531304221817550e+01 + -9.47135556308831283e+00 + -7.80687492596315167e+01 + 5.77908380638228465e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 18 + 19 + + + + + + 9 + 7 + + 6.21423423624779545e+01 + 1.02374321697767794e+01 + -2.61952819865811541e+01 + 1.64810134236028816e+01 + -9.33026901138422886e+01 + -2.97993293115291316e+01 + 6.47591857034336442e+01 + 2.72657426895055011e+01 + -2.33482872492159856e+01 + -6.84639855508772825e+01 + -1.58264297171629540e+01 + 2.10899487969025330e+01 + -3.01026938649666285e+01 + 2.36828763275138563e+01 + -9.22944706954510821e+01 + 6.50721593326321113e+01 + 1.34399733570446998e+01 + -1.47336875153178202e+01 + 1.81996160352434408e+01 + 5.50471142576061965e+01 + 8.12609439867173649e+01 + 3.09138639218085887e+00 + -1.18364488818499112e+01 + 7.45772273447166401e-02 + 3.87852402909780736e+01 + -7.92723675773540748e+01 + 4.57949875755143552e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -7.14353163221378225e+01 + -6.97621833385175307e+01 + 5.49848686450887225e+00 + 4.54649938592219129e+01 + -4.02948505639650421e+01 + 7.94308463470498509e+01 + 7.14353163221378225e+01 + 6.97621833385175307e+01 + -5.49848686450887225e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -5.31970855906319784e+01 + 5.92415630606653920e+01 + 6.05021263328172836e+01 + -4.54649938592219129e+01 + 4.02948505639650421e+01 + -7.94308463470498509e+01 + 5.31970855906319784e+01 + -5.92415630606653920e+01 + -6.05021263328172836e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 8.66602846584249278e+01 + -1.40571348898168708e+02 + -1.29029667681538402e+01 + -4.84107069228954998e+01 + -6.46593270163667029e+01 + 1.03797317950936275e+02 + 6.32154043923236131e+01 + 1.69243370521808174e+01 + -9.00005277714220426e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 19 + 20 + + + + + + 9 + 7 + + 5.40325063921370585e+01 + 5.12983919091620777e+00 + 4.24856501928409855e+01 + -6.01366381168498805e+01 + -6.66348312341534665e+01 + 1.23707663998909716e+01 + 2.51821654645419386e+01 + -4.49879933093015580e+01 + 6.61923323512717872e+01 + -8.16765559697753787e+01 + -1.72226454079742588e+01 + 3.80134460532522880e+01 + -3.51922609512735107e+01 + -2.39876166901732191e+01 + 7.97471671594631992e+00 + -4.83518374187914013e+00 + 8.66273522162050824e+01 + 4.67131848907273408e+01 + 1.17456223558336159e+01 + -5.59162876565777989e+00 + 8.16530126090797523e+01 + 6.71811935630704937e+01 + -7.04411581080473894e+01 + -2.09183319463410733e+01 + -1.56755267575532145e+01 + 1.18327312907135145e+01 + -5.35254607305620027e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -0.00000000000000000e+00 + 1.72947781738215802e+01 + -7.37458142645349284e+01 + 6.52874071041194242e+01 + -9.78660995310737718e+01 + -5.39934004805816503e+00 + 1.98260860892763304e+01 + -1.72947781738215802e+01 + 7.37458142645349284e+01 + -6.52874071041194242e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -1.10958195052130399e+01 + -6.73231164274660614e+01 + -7.31059558722925402e+01 + 9.78660995310737718e+01 + 5.39934004805816503e+00 + -1.98260860892763304e+01 + 1.10958195052130399e+01 + 6.73231164274660614e+01 + 7.31059558722925402e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -5.38048032045895752e+00 + 1.65536142881929180e+02 + 8.22017033006773943e+01 + -7.27309788858385105e+01 + -1.01358249880792499e+01 + 1.16499706864703938e+02 + 7.80821322304104655e+01 + 9.19602704161340796e+01 + -8.91929751904364565e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 20 + 21 + + + + + + 9 + 7 + + 6.79396504738739537e+01 + 4.82037868237773779e+01 + -5.47024548652807283e+01 + -7.00832895591236564e+01 + 5.14309272143008727e+01 + -3.67255818969308834e+01 + 1.86928047352814559e+01 + -2.23350198922307115e+01 + 1.77311293144151279e+01 + 7.44105813075135192e+00 + 5.99087969415584212e+01 + 5.30460372996668283e+01 + -1.22158131494639797e+01 + -1.50072994727732638e+01 + -6.43870652455934760e+01 + -1.47523300546991756e+01 + 7.84637907473050689e+01 + -5.14524303561770466e+01 + 1.48130048314395335e+00 + -3.96193494527523526e+01 + -4.51446929817367320e+01 + 2.41623483145540554e+01 + -4.80504534286710054e+01 + -6.06763607845337773e+01 + 9.54718801387472098e+01 + 2.61345793187026096e+01 + 4.22977195343479639e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -4.31638784413171734e+01 + 8.58566719453584426e+01 + 2.76678781327787284e+01 + -6.66070179567136051e+01 + -9.65155343163612223e+00 + -7.39618325575448523e+01 + 4.31638784413171734e+01 + -8.58566719453584426e+01 + -2.76678781327787284e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -6.08307879023215037e+01 + -5.03535440542508113e+01 + 6.13525536906117850e+01 + 6.66070179567136051e+01 + 9.65155343163612223e+00 + 7.39618325575448523e+01 + 6.08307879023215037e+01 + 5.03535440542508113e+01 + -6.13525536906117850e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.33814642652391399e+02 + 1.77445244940366287e-01 + -1.49228346303260491e+01 + -6.12411970511728132e-01 + 5.97803936488639280e+01 + -1.02665793382055853e+02 + 2.53293581625552058e+01 + -3.42663008134674527e+01 + 5.61197444293587395e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 21 + 22 + + + + + + 9 + 7 + + -6.74405298173737435e+01 + 7.32759993494661614e+01 + -6.72758148030300518e+00 + 2.84952832715093294e+01 + 3.71106482095568992e+01 + 5.87105666468560798e+01 + 3.06497832963929113e+01 + 2.67945705145493136e+01 + 5.23800862467740629e+01 + -2.34193589858697777e+01 + -2.67958038635302600e+01 + 2.45872057837964846e+01 + 1.16252649553510334e+01 + 8.99424898505579051e+01 + -2.22061693363508326e+01 + 6.98562239082597500e+01 + -3.31539117547968303e+01 + -5.85832970823880999e+01 + 5.89049685210982901e+01 + 4.59874265404548481e+01 + -5.08083954437895216e+01 + -2.22191873368612889e+01 + 1.16539331220696187e+01 + -7.08270509506904347e+01 + 6.46018011815524602e+01 + 1.94415681363770645e+01 + 4.04539860386468320e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 9.78145711723120286e+01 + 1.76232674770143696e+01 + -1.10331368980137903e+01 + -1.40616464377447681e+01 + 1.69815554662718142e+01 + -9.75392068524539155e+01 + -9.78145711723120286e+01 + -1.76232674770143696e+01 + 1.10331368980137903e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -1.53159970565603700e+01 + 9.69589976091932186e+01 + 1.90885572211179095e+01 + 1.40616464377447681e+01 + -1.69815554662718142e+01 + 9.75392068524539155e+01 + 1.53159970565603700e+01 + -9.69589976091932186e+01 + -1.90885572211179095e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -2.25457773322750015e+01 + -1.39359117049240723e+02 + -1.01360767749937452e+02 + -1.06577486560693373e+02 + 3.69148207609075527e+01 + -6.52361857149522422e+01 + -9.51115531141346224e+01 + 7.07220482565039248e+01 + -3.57407782543718326e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 22 + 23 + + + + + + 9 + 7 + + 9.41929589378244003e+01 + 8.94097901372627391e+00 + 3.01620867599012286e+01 + 2.46922319059140740e+01 + 1.35465056881596215e+01 + -4.89963411381141611e+01 + -1.68344538935716344e+01 + -1.83167626535859007e+01 + 7.95313381974681306e+01 + -5.18828216803281350e+00 + -4.19655627938833504e+01 + -6.55823195034947481e+00 + 7.35562337314085681e+01 + 5.97403143502788865e+01 + -2.03387857576245445e-01 + 6.60668686562018337e+01 + -6.61937508072900442e+01 + -2.10487850516630495e+01 + 1.54947760148805180e+01 + -8.96286058833042887e+01 + -5.79472767879112993e+00 + -6.70795355915261116e+00 + -1.69646265262306706e+01 + -8.65671984840189168e+01 + -1.21336256996353509e+01 + 3.69909891189522924e+01 + -4.84851501469951671e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -0.00000000000000000e+00 + 8.51395966446430350e+01 + -5.19584580090174413e+01 + 7.18106708732509347e+00 + -5.24346283807042397e+01 + -8.39535649250525182e+01 + 1.42270405549584336e+01 + -8.51395966446430350e+01 + 5.19584580090174413e+01 + -7.18106708732509347e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -1.36338907320492919e+00 + -1.58782107839693651e+01 + -9.87219509153607788e+01 + 5.24346283807042397e+01 + 8.39535649250525182e+01 + -1.42270405549584336e+01 + 1.36338907320492919e+00 + 1.58782107839693651e+01 + 9.87219509153607788e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 3.07007154914777551e+01 + 2.70881185233511914e+01 + 3.78166527462940882e+00 + -1.15162420127751133e+02 + -6.75390943713141123e+00 + 3.95885323960980529e+00 + -1.29970649317129876e+01 + 1.14524822858325365e+02 + 2.24989739182521227e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 23 + 24 + + + + + + 9 + 7 + + -1.83421298069983294e+01 + -8.40178887379047552e+01 + 1.37665553988889080e+01 + 3.46271545199608042e+01 + 4.05922653756383482e+01 + -7.00698095538007770e+00 + 9.19927590625682683e+01 + -3.22667773518094236e+01 + 3.91925773350330608e+00 + 4.07207141973007651e+01 + -3.85719388369921674e+01 + -7.75337815930943464e+01 + 7.24865255322073381e+01 + -3.53413936296266158e+01 + 5.86044610642822335e+01 + -1.83884977004002117e+01 + 1.46625501311897892e+01 + 1.87515674799128185e+01 + -2.53722481762123628e+01 + -3.66484706705800178e+01 + 3.56472023092901580e+01 + -4.40875208799477321e+01 + -6.95087704456724254e+01 + 1.97777311636326232e+01 + 1.27979930474377834e+01 + 3.35520358632383733e+01 + 9.10304934695741821e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -2.81762066375282245e+01 + 9.59328037613137070e+01 + -1.73162987170684679e+00 + -8.31278175295505264e+01 + -2.53085961788951153e+01 + -4.94898061445724977e+01 + 2.81762066375282245e+01 + -9.59328037613137070e+01 + 1.73162987170684679e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -4.79152098220707003e+01 + -1.25048839237672293e+01 + 8.68778484181097070e+01 + 8.31278175295505264e+01 + 2.53085961788951153e+01 + 4.94898061445724977e+01 + 4.79152098220707003e+01 + 1.25048839237672293e+01 + -8.68778484181097070e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -3.78847831737257366e+01 + 2.29972416224487830e+01 + -1.37159770597291157e+02 + -4.30909041970757869e+01 + -7.29718405033518991e+01 + -1.27753851859623239e+02 + -4.00153271264082644e+00 + -1.83214678716897197e+02 + 3.84258732608477231e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 24 + 25 + + + + + + 9 + 7 + + -1.42022591781322873e+01 + -3.56662961537267762e+01 + -1.29665633112741379e+01 + -4.42547548664715933e+01 + 8.53484085500931968e+01 + -1.21228131477824768e+01 + 7.44576769541145751e+01 + 3.58775048496875613e+01 + 4.62292669635554603e+01 + 1.27026107710752729e+01 + -9.11117710856913021e+01 + 2.52571349253896464e+01 + -3.19164772916096373e+00 + -2.63470878340458832e+01 + 1.88346366483448016e+01 + 5.40205809637968812e+01 + -1.90775864349462978e+01 + -8.09820522299282857e+01 + 6.16280452912953720e+01 + 1.93876135715390099e+01 + 7.12990641830701719e+01 + 6.58414093447228908e+01 + 1.80759831933341069e+01 + -6.99233759491107918e+01 + 3.27566373625739402e+01 + 9.70961890249048842e+00 + 4.91628253265810322e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 9.34197183252336458e+01 + 2.76233972145661149e+01 + 2.25766284985226449e+01 + -2.39170403241009595e+01 + 1.53702168542746431e+00 + 9.70855949483436831e+01 + -9.34197183252336458e+01 + -2.76233972145661149e+01 + -2.25766284985226449e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 2.64713318548450189e+01 + -9.60967506769340645e+01 + 8.04258037981846208e+00 + 2.39170403241009595e+01 + -1.53702168542746431e+00 + -9.70855949483436831e+01 + -2.64713318548450189e+01 + 9.60967506769340645e+01 + -8.04258037981846208e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 4.99415394471549448e+01 + 8.89774612805008616e+01 + -7.21549206200566857e+01 + -8.47137695131298614e+01 + -3.96302269514171712e+01 + 2.91958404364136328e+01 + -1.14979727399682560e+02 + 6.32312875848932592e+01 + -5.83619937525406982e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 25 + 26 + + + + + + 9 + 7 + + 4.75279810851286726e+01 + 1.27698906927141866e+01 + -3.85248964418578721e+00 + 8.38265009034030300e+01 + 2.34263445942633055e+01 + 4.00325151654156208e+00 + -3.67133537961906153e+00 + 7.31784749045807636e+00 + 9.95211070535517166e+01 + 7.72966987607180016e+01 + 3.83523138848579492e+01 + 1.84348658933893006e+01 + -5.00282378767150675e+01 + 1.62934700342947103e+00 + -1.23502536597363957e+01 + 3.08280910111897235e+01 + -9.23259671614170685e+01 + 9.05617058795025009e+00 + -3.32567607414421573e+01 + 4.65581163028938505e+01 + 8.06472888898161528e+01 + -1.09449242551094059e+01 + 8.15866317932763252e+01 + -5.45355209481987444e+01 + -6.85221032691171050e+00 + 2.03202505385987919e+01 + 3.49109728776780415e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -0.00000000000000000e+00 + 5.75486223700384372e+01 + -3.56976774925786202e+01 + 7.35787461495951334e+01 + -7.77278254550152923e+01 + 4.09823526552186745e+00 + 6.27820803871861983e+01 + -5.75486223700384372e+01 + 3.56976774925786202e+01 + -7.35787461495951334e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -2.54271747023806824e+01 + -9.33213817372218273e+01 + -2.53885505161348775e+01 + 7.77278254550152923e+01 + -4.09823526552186745e+00 + -6.27820803871861983e+01 + 2.54271747023806824e+01 + 9.33213817372218273e+01 + 2.53885505161348775e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 5.11234651369603128e+01 + 1.45916426283882128e+01 + 8.14361101436695947e+01 + -9.64495575814352151e+01 + 5.69392664765279974e+01 + 1.45591038998347216e+02 + -1.52360008121073719e+02 + 3.92577375478538881e+00 + -7.47401324352360774e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 7 + 10 + + + + + + 9 + 7 + + -5.40766118135344485e+01 + -6.55945677789347137e+00 + 5.92166444746720977e-01 + 8.33981325801565987e+01 + -1.15536936327200639e+01 + 1.11233523720273890e+01 + -9.80808079024424018e+00 + -8.64371084278765380e+01 + 4.74584491239462309e+01 + 1.25135196935038806e+01 + 9.77748374151099995e+01 + 6.14999725352584115e+00 + 1.36807989121553284e+01 + 4.43989849701706429e+00 + -9.89599830977766572e+01 + -2.01505456436379626e+00 + -1.55263904030224698e+01 + -7.65145679760131925e-01 + 8.31600037168701789e+01 + -1.88305696684158264e+01 + -2.79941305069657931e+00 + 5.18811705043763141e+01 + -7.17174238817171439e+00 + 6.66998670151083850e+00 + -8.31970807941664781e+00 + -4.61019176373218968e+01 + -8.79761519740697651e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 7.83449748282976799e+01 + 6.13201754047114846e+01 + 1.00946028891096997e+01 + -1.98027611608243852e+01 + 9.23644092983748521e+00 + 9.75834966034573341e+01 + -7.83449748282976799e+01 + -6.13201754047114846e+01 + -1.00946028891096997e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 5.89059892503364324e+01 + -7.84507759508154265e+01 + 1.93793752003053505e+01 + 1.98027611608243852e+01 + -9.23644092983748521e+00 + -9.75834966034573341e+01 + -5.89059892503364324e+01 + 7.84507759508154265e+01 + -1.93793752003053505e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -6.07889604414481113e+01 + 8.49672366649234476e+01 + -1.19149060367207383e+02 + -3.26841768394553398e+01 + -8.97350957633814517e+01 + 9.13622785379021281e+01 + -1.77493968880864855e+02 + -4.85855715361152676e+01 + -7.39158595946127228e+00 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 9 + 14 + + + + + + 9 + 7 + + 4.43615564152903019e+01 + -3.09502672537490433e+01 + -6.65256070847414378e+01 + 4.96695435786911474e+01 + 2.78172662190901541e+01 + 6.19720498005987466e+01 + -7.45919697090238145e+01 + 1.35179597700267684e+00 + 1.13620760877569249e+00 + 5.07131897166802617e+01 + -6.73003419642306540e+01 + 3.11633063652204427e+01 + 5.47056695452382584e+01 + 6.11187671649742512e+01 + -2.76098847773204135e+01 + 6.65963623398869657e+01 + 1.43297991149407444e+00 + -1.99391285797603790e-01 + 7.65967691022320984e-01 + 1.84943157433388130e+01 + -6.50656556956834748e+01 + 8.44674932355512387e-01 + 2.30555803308049860e+01 + -6.99700625824593629e+01 + -3.40876610014427706e-01 + -9.55299072188111609e+01 + -2.94902469873441930e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -2.89030917873941853e+01 + 8.70749520313458305e+01 + -3.97814531392241904e+01 + -9.56630827209552450e+01 + -2.78467078053438328e+01 + 8.55192777842975183e+00 + 2.89030917873941853e+01 + -8.70749520313458305e+01 + 3.97814531392241904e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -3.63123800557649812e+00 + 4.05279359595652906e+01 + 9.13471429077244466e+01 + 9.56630827209552450e+01 + 2.78467078053438328e+01 + -8.55192777842975183e+00 + 3.63123800557649812e+00 + -4.05279359595652906e+01 + -9.13471429077244466e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -7.02581401004128878e+01 + -1.05108950320020767e+02 + -1.10573546354685661e+02 + -2.82871715200345584e+01 + -9.82678476821147342e+01 + -1.38389269150812240e+01 + 2.98238876860104902e+01 + -1.16573384649075678e+02 + 1.35327233223631623e+02 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 10 + 25 + + + + + + 9 + 7 + + 3.95766223652825175e+01 + 2.36448299245490254e+01 + -2.45799893899839930e+01 + -7.00228044208398188e+01 + -4.93985238797199173e+01 + 1.28698697829822883e+01 + 7.78345939301136092e+00 + 2.90892349948341078e+01 + 9.40990441114148695e+01 + -8.42616415682074660e+01 + -2.90334966132265606e+01 + -1.43512687624874875e+01 + -2.62231810185785967e+01 + -4.77560189769582593e+01 + -8.01321381020614343e-01 + 4.26634609899314583e+01 + -8.24081707666526171e+01 + 1.64578293010562966e+01 + -8.45857606147821528e+00 + 2.58771125144124774e+01 + -9.15441805767733570e+01 + 5.92769857321209486e+01 + -7.03919737513036381e+01 + -3.24535577243972782e+01 + -1.66359962010244793e+01 + 2.48028458121687763e+01 + -2.15529183777999371e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 9.34197183252336458e+01 + 2.76233972145661149e+01 + 2.25766284985226449e+01 + -2.39170403241009595e+01 + 1.53702168542746431e+00 + 9.70855949483436831e+01 + -9.34197183252336458e+01 + -2.76233972145661149e+01 + -2.25766284985226449e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 2.64713318548450189e+01 + -9.60967506769340645e+01 + 8.04258037981846208e+00 + 2.39170403241009595e+01 + -1.53702168542746431e+00 + -9.70855949483436831e+01 + -2.64713318548450189e+01 + 9.60967506769340645e+01 + -8.04258037981846208e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -6.19976450903283478e+01 + 1.85138049290959714e+02 + 2.04096238885501400e+01 + -5.38401179783636792e+01 + 1.95403129638118465e+01 + 3.38132128740523674e+00 + -1.81977553628765463e+02 + -6.94268461743989320e+01 + -2.32863678513998718e+00 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 13 + 22 + + + + + + 9 + 7 + + 1.61359876726555882e+01 + 1.08711432986265653e+01 + 8.51184696478907732e+00 + 6.31122752697449414e+01 + -1.20771880436756209e+01 + 7.50149628126125521e+01 + -4.10840851091186110e+01 + -8.68858010361213076e+01 + 2.35731234098279678e+01 + 4.08177879353251924e+01 + -2.17654374458011439e+01 + 8.78457730201184290e+01 + -2.17842465512247614e+01 + 1.46672234096332783e+01 + 6.07972286227193015e-01 + 6.89954693463132855e+01 + -4.93810369540159471e+01 + -4.74803457645439977e+01 + -6.20341282977744939e+01 + -7.54126088063309794e+01 + 1.25298075512294904e+01 + -5.28542109714301915e+01 + 5.80638264069757923e+01 + 5.82126617561848008e+01 + -2.29292298928228675e+01 + -1.00045379878072604e+00 + 1.50754924253946800e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 9.78145711723120286e+01 + 1.76232674770143696e+01 + -1.10331368980137903e+01 + -1.40616464377447681e+01 + 1.69815554662718142e+01 + -9.75392068524539155e+01 + -9.78145711723120286e+01 + -1.76232674770143696e+01 + 1.10331368980137903e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -1.53159970565603700e+01 + 9.69589976091932186e+01 + 1.90885572211179095e+01 + 1.40616464377447681e+01 + -1.69815554662718142e+01 + 9.75392068524539155e+01 + 1.53159970565603700e+01 + -9.69589976091932186e+01 + -1.90885572211179095e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 8.03183434009733901e+01 + -1.57959181406014523e+02 + -6.44060214101481705e+01 + -6.64898917378079801e+01 + -6.21866664385735888e+01 + -6.61756460634393022e+01 + -1.52823432141161049e+02 + -1.42511660560011357e+01 + -7.24093903591090111e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 20 + 15 + + + + + + 9 + 7 + + 4.71343070271561331e+01 + -6.72091990499062604e+01 + 5.47388912655424136e+01 + 8.35500845890857278e+01 + 2.23841738875626213e+01 + -3.32945376880477681e+01 + -2.59734415135038539e+01 + -2.12006909898996589e+01 + 2.34658474313030787e+01 + 6.32652845860905977e+00 + 9.30238234940692976e+00 + 3.37625583972403689e+01 + 3.44554933454442613e+00 + 8.79955033997688929e+01 + 4.09722885299349642e+01 + -1.92439510927881194e+01 + 4.62720459594682580e+01 + -8.23910435376807584e+01 + -5.55999617765192156e+00 + 6.20368695338309522e+01 + 7.15470037167418553e+01 + 3.41059120969163914e+01 + -1.39920442875337869e+01 + -2.48395357584865195e+01 + 9.22321060162172301e+01 + 2.27421056426905714e+01 + 1.25600461547882558e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -8.15059339771809164e+01 + -5.54598536753216109e+01 + -1.67626774955357369e+01 + -1.04808841952049807e+01 + 4.25686446848321935e+01 + -8.98780371179923492e+01 + 8.15059339771809164e+01 + 5.54598536753216109e+01 + 1.67626774955357369e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 5.69818724946287603e+01 + -7.14990567770542356e+01 + -4.05086544703406517e+01 + 1.04808841952049807e+01 + -4.25686446848321935e+01 + 8.98780371179923492e+01 + -5.69818724946287603e+01 + 7.14990567770542356e+01 + 4.05086544703406517e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -0.00000000000000000e+00 + -1.44773765128431023e+02 + 3.21567907629271019e+01 + 6.77990836285083986e+01 + 3.24643872643349098e+01 + 3.07230481384375675e+00 + -8.66387073171390654e+00 + 6.03331707496851593e+01 + 1.38479259212844596e+02 + 6.83329111783498746e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 27 + 0 + 0 + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + + + + + + 3 + 82 + + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + + + + 29 + 0 + 0 + 3 + 6 + 9 + 12 + 15 + 18 + 21 + 24 + 27 + 30 + 33 + 36 + 39 + 42 + 45 + 48 + 51 + 54 + 57 + 60 + 63 + 66 + 69 + 72 + 75 + 78 + 81 + 82 + + 0 + 3 + 0 + + + + + + + + + + + diff --git a/examples/Data/toy3D.xml b/examples/Data/toy3D.xml new file mode 100644 index 000000000..13dbcbe6c --- /dev/null +++ b/examples/Data/toy3D.xml @@ -0,0 +1,169 @@ + + + + + + + 2 + 1 + + + + + + 2 + 0 + 0 + 1 + + + + + + 9 + 7 + + 1.22427709071730959e+01 + 1.51514613104920048e+01 + 3.60934366857813060e+00 + -1.18407259026383116e+01 + 7.84826117220921216e+00 + 1.23509165494819051e+01 + -6.09875015991639735e+00 + 6.16547190708139126e-01 + 3.94972084922329048e+00 + -4.89208482920378174e+00 + 3.02091647632478866e+00 + -8.95328692238917512e+00 + 7.89831607220345955e+00 + -2.36793602009719084e+00 + 1.48517612051941725e+01 + -3.97284286249233731e-01 + -1.95744531643153863e+01 + -3.85954855417462017e+00 + 4.79268277145419042e+00 + -9.01707953629520453e+00 + 1.37848069005841385e+01 + 1.04829326688375950e+01 + -5.00630568442241675e+00 + 4.70463561852773182e+00 + -1.59179134598689274e+01 + -2.04767784956723942e+00 + 9.54135497908261954e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + -0.00000000000000000e+00 + -1.76201757312015372e+00 + 1.98634190821282672e+01 + 1.52966546661624236e+00 + 1.94817649567575373e+01 + 1.39684693294110307e+00 + 4.30228460420588288e+00 + 1.76201757312015372e+00 + -1.98634190821282672e+01 + -1.52966546661624236e+00 + -0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 4.16606867942304948e+00 + 1.86906420801308037e+00 + -1.94717865319198360e+01 + -1.94817649567575373e+01 + -1.39684693294110307e+00 + -4.30228460420588288e+00 + -4.16606867942304948e+00 + -1.86906420801308037e+00 + 1.94717865319198360e+01 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00891462904330531e+01 + -1.08132497987816869e+01 + 8.66487736568128497e+00 + 2.88370015604634311e+01 + 1.89391698948574643e+01 + 2.12398960190661290e+00 + 1.22150946112124039e+01 + -2.33658532501596561e+01 + 1.51576204760307363e+01 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 9 + 0 + + + + + + + + + + + + 2 + 0 + 0 + 1 + + + + + + 3 + 7 + + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 1.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + 0.00000000000000000e+00 + + + + 4 + 0 + 0 + 3 + 6 + 7 + + 0 + 3 + 0 + + + + + + + + + + + diff --git a/gtsam/base/numericalDerivative.h b/gtsam/base/numericalDerivative.h index cc1cbdb51..a9a088108 100644 --- a/gtsam/base/numericalDerivative.h +++ b/gtsam/base/numericalDerivative.h @@ -67,28 +67,29 @@ struct FixedSizeMatrix { } /** - * Numerically compute gradient of scalar function + * @brief Numerically compute gradient of scalar function + * @return n-dimensional gradient computed via central differencing * Class X is the input argument * The class X needs to have dim, expmap, logmap + * int N is the dimension of the X input value if variable dimension type but known at test time */ -template -typename internal::FixedSizeMatrix::type numericalGradient(boost::function h, const X& x, - double delta = 1e-5) { + +template ::dimension> +typename Eigen::Matrix numericalGradient( + boost::function h, const X& x, double delta = 1e-5) { double factor = 1.0 / (2.0 * delta); BOOST_STATIC_ASSERT_MSG( (boost::is_base_of::structure_category>::value), "Template argument X must be a manifold type."); - static const int N = traits::dimension; - BOOST_STATIC_ASSERT_MSG(N>0, "Template argument X must be fixed-size type."); - - typedef typename traits::TangentVector TangentX; + BOOST_STATIC_ASSERT_MSG(N>0, "Template argument X must be fixed-size type or N must be specified."); // Prepare a tangent vector to perturb x with, only works for fixed size - TangentX d; + typename traits::TangentVector d; d.setZero(); - Eigen::Matrix g; g.setZero(); // Can be fixed size + Eigen::Matrix g; + g.setZero(); for (int j = 0; j < N; j++) { d(j) = delta; double hxplus = h(traits::Retract(x, d)); @@ -108,37 +109,34 @@ typename internal::FixedSizeMatrix::type numericalGradient(boost::function +template ::dimension> // TODO Should compute fixed-size matrix -typename internal::FixedSizeMatrix::type numericalDerivative11(boost::function h, const X& x, - double delta = 1e-5) { - +typename internal::FixedSizeMatrix::type numericalDerivative11( + boost::function h, const X& x, double delta = 1e-5) { typedef typename internal::FixedSizeMatrix::type Matrix; BOOST_STATIC_ASSERT_MSG( (boost::is_base_of::structure_category>::value), "Template argument Y must be a manifold type."); typedef traits TraitsY; - typedef typename TraitsY::TangentVector TangentY; BOOST_STATIC_ASSERT_MSG( (boost::is_base_of::structure_category>::value), "Template argument X must be a manifold type."); - static const int N = traits::dimension; - BOOST_STATIC_ASSERT_MSG(N>0, "Template argument X must be fixed-size type."); + BOOST_STATIC_ASSERT_MSG(N>0, "Template argument X must be fixed-size type or N must be specified."); typedef traits TraitsX; - typedef typename TraitsX::TangentVector TangentX; // get value at x, and corresponding chart const Y hx = h(x); // Bit of a hack for now to find number of rows - const TangentY zeroY = TraitsY::Local(hx, hx); + const typename TraitsY::TangentVector zeroY = TraitsY::Local(hx, hx); const size_t m = zeroY.size(); // Prepare a tangent vector to perturb x with, only works for fixed size - TangentX dx; + Eigen::Matrix dx; dx.setZero(); // Fill in Jacobian H @@ -146,9 +144,9 @@ typename internal::FixedSizeMatrix::type numericalDerivative11(boost::funct const double factor = 1.0 / (2.0 * delta); for (int j = 0; j < N; j++) { dx(j) = delta; - const TangentY dy1 = TraitsY::Local(hx, h(TraitsX::Retract(x, dx))); + const auto dy1 = TraitsY::Local(hx, h(TraitsX::Retract(x, dx))); dx(j) = -delta; - const TangentY dy2 = TraitsY::Local(hx, h(TraitsX::Retract(x, dx))); + const auto dy2 = TraitsY::Local(hx, h(TraitsX::Retract(x, dx))); dx(j) = 0; H.col(j) << (dy1 - dy2) * factor; } diff --git a/gtsam/linear/Preconditioner.h b/gtsam/linear/Preconditioner.h index 623b29863..31901db3f 100644 --- a/gtsam/linear/Preconditioner.h +++ b/gtsam/linear/Preconditioner.h @@ -70,21 +70,20 @@ public: Preconditioner() {} virtual ~Preconditioner() {} - /* Computation Interfaces */ + /* + * Abstract interface for raw vectors. VectorValues is a speed bottleneck + * and Yong-Dian has profiled preconditioners (outside GTSAM) with the the + * three methods below. In GTSAM, unfortunately, we are still using the + * VectorValues methods called in iterative-inl.h + */ - /* implement x = L^{-1} y */ + /// implement x = L^{-1} y virtual void solve(const Vector& y, Vector &x) const = 0; -// virtual void solve(const VectorValues& y, VectorValues &x) const = 0; - /* implement x = L^{-T} y */ + /// implement x = L^{-T} y virtual void transposeSolve(const Vector& y, Vector& x) const = 0; -// virtual void transposeSolve(const VectorValues& y, VectorValues &x) const = 0; -// /* implement x = L^{-1} L^{-T} y */ -// virtual void fullSolve(const Vector& y, Vector &x) const = 0; -// virtual void fullSolve(const VectorValues& y, VectorValues &x) const = 0; - - /* build/factorize the preconditioner */ + /// build/factorize the preconditioner virtual void build( const GaussianFactorGraph &gfg, const KeyInfo &info, @@ -113,14 +112,7 @@ public: /* Computation Interfaces for raw vector */ virtual void solve(const Vector& y, Vector &x) const { x = y; } -// virtual void solve(const VectorValues& y, VectorValues& x) const { x = y; } - virtual void transposeSolve(const Vector& y, Vector& x) const { x = y; } -// virtual void transposeSolve(const VectorValues& y, VectorValues& x) const { x = y; } - -// virtual void fullSolve(const Vector& y, Vector &x) const { x = y; } -// virtual void fullSolve(const VectorValues& y, VectorValues& x) const { x = y; } - virtual void build( const GaussianFactorGraph &gfg, const KeyInfo &info, @@ -145,8 +137,6 @@ public: /* Computation Interfaces for raw vector */ virtual void solve(const Vector& y, Vector &x) const; virtual void transposeSolve(const Vector& y, Vector& x) const ; -// virtual void fullSolve(const Vector& y, Vector &x) const ; - virtual void build( const GaussianFactorGraph &gfg, const KeyInfo &info, diff --git a/gtsam/linear/SubgraphPreconditioner.cpp b/gtsam/linear/SubgraphPreconditioner.cpp index d796e28b7..657846f95 100644 --- a/gtsam/linear/SubgraphPreconditioner.cpp +++ b/gtsam/linear/SubgraphPreconditioner.cpp @@ -52,17 +52,21 @@ #include #include -using namespace std; +using std::cout; +using std::endl; +using std::vector; +using std::ostream; namespace gtsam { /* ************************************************************************* */ +// Convert any non-Jacobian factors to Jacobians (e.g. Hessian -> Jacobian with Cholesky) static GaussianFactorGraph::shared_ptr convertToJacobianFactors(const GaussianFactorGraph &gfg) { - GaussianFactorGraph::shared_ptr result(new GaussianFactorGraph()); - for(const GaussianFactor::shared_ptr &gf: gfg) { - JacobianFactor::shared_ptr jf = boost::dynamic_pointer_cast(gf); + auto result = boost::make_shared(); + for (const auto &factor : gfg) { + auto jf = boost::dynamic_pointer_cast(factor); if( !jf ) { - jf = boost::make_shared(*gf); // Convert any non-Jacobian factors to Jacobians (e.g. Hessian -> Jacobian with Cholesky) + jf = boost::make_shared(*factor); } result->push_back(jf); } @@ -70,20 +74,23 @@ static GaussianFactorGraph::shared_ptr convertToJacobianFactors(const GaussianFa } /*****************************************************************************/ -static std::vector iidSampler(const vector &weight, const size_t n) { +static vector iidSampler(const vector &weight, const size_t n) { /* compute the sum of the weights */ const double sum = std::accumulate(weight.begin(), weight.end(), 0.0); + if (sum==0.0) { + throw std::runtime_error("Weight vector has no non-zero weights"); + } /* make a normalized and accumulated version of the weight vector */ const size_t m = weight.size(); - vector w; w.reserve(m); + vector cdf; cdf.reserve(m); for ( size_t i = 0 ; i < m ; ++i ) { - w.push_back(weight[i]/sum); + cdf.push_back(weight[i]/sum); } vector acc(m); - std::partial_sum(w.begin(),w.end(),acc.begin()); + std::partial_sum(cdf.begin(),cdf.end(),acc.begin()); /* iid sample n times */ vector result; result.reserve(n); @@ -91,26 +98,26 @@ static std::vector iidSampler(const vector &weight, const size_t for ( size_t i = 0 ; i < n ; ++i ) { const double value = rand() / denominator; /* binary search the interval containing "value" */ - vector::iterator it = std::lower_bound(acc.begin(), acc.end(), value); - size_t idx = it - acc.begin(); + const auto it = std::lower_bound(acc.begin(), acc.end(), value); + const size_t idx = it - acc.begin(); result.push_back(idx); } return result; } /*****************************************************************************/ -vector uniqueSampler(const vector &weight, const size_t n) { +static vector UniqueSampler(const vector &weight, const size_t n) { const size_t m = weight.size(); - if ( n > m ) throw std::invalid_argument("uniqueSampler: invalid input size"); + if ( n > m ) throw std::invalid_argument("UniqueSampler: invalid input size"); vector result; size_t count = 0; - std::vector touched(m, false); + vector touched(m, false); while ( count < n ) { - std::vector localIndices; localIndices.reserve(n-count); - std::vector localWeights; localWeights.reserve(n-count); + vector localIndices; localIndices.reserve(n-count); + vector localWeights; localWeights.reserve(n-count); /* collect data */ for ( size_t i = 0 ; i < m ; ++i ) { @@ -134,16 +141,16 @@ vector uniqueSampler(const vector &weight, const size_t n) { } /****************************************************************************/ -Subgraph::Subgraph(const std::vector &indices) { +Subgraph::Subgraph(const vector &indices) { edges_.reserve(indices.size()); for ( const size_t &idx: indices ) { - edges_.push_back(SubgraphEdge(idx, 1.0)); + edges_.emplace_back(idx, 1.0); } } /****************************************************************************/ -std::vector Subgraph::edgeIndices() const { - std::vector eid; eid.reserve(size()); +vector Subgraph::edgeIndices() const { + vector eid; eid.reserve(size()); for ( const SubgraphEdge &edge: edges_ ) { eid.push_back(edge.index_); } @@ -169,7 +176,7 @@ Subgraph::shared_ptr Subgraph::load(const std::string &fn) { } /****************************************************************************/ -std::ostream &operator<<(std::ostream &os, const SubgraphEdge &edge) { +ostream &operator<<(ostream &os, const SubgraphEdge &edge) { if ( edge.weight() != 1.0 ) os << edge.index() << "(" << std::setprecision(2) << edge.weight() << ")"; else @@ -178,7 +185,7 @@ std::ostream &operator<<(std::ostream &os, const SubgraphEdge &edge) { } /****************************************************************************/ -std::ostream &operator<<(std::ostream &os, const Subgraph &subgraph) { +ostream &operator<<(ostream &os, const Subgraph &subgraph) { os << "Subgraph" << endl; for ( const SubgraphEdge &e: subgraph.edges() ) { os << e << ", " ; @@ -212,16 +219,16 @@ SubgraphBuilderParameters::Skeleton SubgraphBuilderParameters::skeletonTranslato if (s == "NATURALCHAIN") return NATURALCHAIN; else if (s == "BFS") return BFS; else if (s == "KRUSKAL") return KRUSKAL; - throw invalid_argument("SubgraphBuilderParameters::skeletonTranslator undefined string " + s); + throw std::invalid_argument("SubgraphBuilderParameters::skeletonTranslator undefined string " + s); return KRUSKAL; } /****************************************************************/ -std::string SubgraphBuilderParameters::skeletonTranslator(Skeleton w) { - if ( w == NATURALCHAIN )return "NATURALCHAIN"; - else if ( w == BFS ) return "BFS"; - else if ( w == KRUSKAL )return "KRUSKAL"; - else return "UNKNOWN"; +std::string SubgraphBuilderParameters::skeletonTranslator(Skeleton s) { + if ( s == NATURALCHAIN ) return "NATURALCHAIN"; + else if ( s == BFS ) return "BFS"; + else if ( s == KRUSKAL ) return "KRUSKAL"; + else return "UNKNOWN"; } /****************************************************************/ @@ -231,7 +238,7 @@ SubgraphBuilderParameters::SkeletonWeight SubgraphBuilderParameters::skeletonWei else if (s == "RHS") return RHS_2NORM; else if (s == "LHS") return LHS_FNORM; else if (s == "RANDOM") return RANDOM; - throw invalid_argument("SubgraphBuilderParameters::skeletonWeightTranslator undefined string " + s); + throw std::invalid_argument("SubgraphBuilderParameters::skeletonWeightTranslator undefined string " + s); return EQUAL; } @@ -245,12 +252,14 @@ std::string SubgraphBuilderParameters::skeletonWeightTranslator(SkeletonWeight w } /****************************************************************/ -SubgraphBuilderParameters::AugmentationWeight SubgraphBuilderParameters::augmentationWeightTranslator(const std::string &src) { +SubgraphBuilderParameters::AugmentationWeight +SubgraphBuilderParameters::augmentationWeightTranslator( + const std::string &src) { std::string s = src; boost::algorithm::to_upper(s); if (s == "SKELETON") return SKELETON; // else if (s == "STRETCH") return STRETCH; // else if (s == "GENERALIZED_STRETCH") return GENERALIZED_STRETCH; - throw invalid_argument("SubgraphBuilder::Parameters::augmentationWeightTranslator undefined string " + s); + throw std::invalid_argument("SubgraphBuilder::Parameters::augmentationWeightTranslator undefined string " + s); return SKELETON; } @@ -263,7 +272,9 @@ std::string SubgraphBuilderParameters::augmentationWeightTranslator(Augmentation } /****************************************************************/ -std::vector SubgraphBuilder::buildTree(const GaussianFactorGraph &gfg, const FastMap &ordering, const std::vector &w) const { +vector SubgraphBuilder::buildTree(const GaussianFactorGraph &gfg, + const FastMap &ordering, + const vector &weights) const { const SubgraphBuilderParameters &p = parameters_; switch (p.skeleton_) { case SubgraphBuilderParameters::NATURALCHAIN: @@ -273,21 +284,21 @@ std::vector SubgraphBuilder::buildTree(const GaussianFactorGraph &gfg, c return bfs(gfg); break; case SubgraphBuilderParameters::KRUSKAL: - return kruskal(gfg, ordering, w); + return kruskal(gfg, ordering, weights); break; default: - cerr << "SubgraphBuilder::buildTree undefined skeleton type" << endl; + std::cerr << "SubgraphBuilder::buildTree undefined skeleton type" << endl; break; } return vector(); } /****************************************************************/ -std::vector SubgraphBuilder::unary(const GaussianFactorGraph &gfg) const { - std::vector result ; +vector SubgraphBuilder::unary(const GaussianFactorGraph &gfg) const { + vector result ; size_t idx = 0; - for ( const GaussianFactor::shared_ptr &gf: gfg ) { - if ( gf->size() == 1 ) { + for (const auto &factor : gfg) { + if ( factor->size() == 1 ) { result.push_back(idx); } idx++; @@ -296,8 +307,8 @@ std::vector SubgraphBuilder::unary(const GaussianFactorGraph &gfg) const } /****************************************************************/ -std::vector SubgraphBuilder::natural_chain(const GaussianFactorGraph &gfg) const { - std::vector result ; +vector SubgraphBuilder::natural_chain(const GaussianFactorGraph &gfg) const { + vector result ; size_t idx = 0; for ( const GaussianFactor::shared_ptr &gf: gfg ) { if ( gf->size() == 2 ) { @@ -311,7 +322,7 @@ std::vector SubgraphBuilder::natural_chain(const GaussianFactorGraph &gf } /****************************************************************/ -std::vector SubgraphBuilder::bfs(const GaussianFactorGraph &gfg) const { +vector SubgraphBuilder::bfs(const GaussianFactorGraph &gfg) const { const VariableIndex variableIndex(gfg); /* start from the first key of the first factor */ Key seed = gfg[0]->keys()[0]; @@ -319,7 +330,7 @@ std::vector SubgraphBuilder::bfs(const GaussianFactorGraph &gfg) const { const size_t n = variableIndex.size(); /* each vertex has self as the predecessor */ - std::vector result; + vector result; result.reserve(n-1); /* Initialize */ @@ -347,30 +358,31 @@ std::vector SubgraphBuilder::bfs(const GaussianFactorGraph &gfg) const { } /****************************************************************/ -std::vector SubgraphBuilder::kruskal(const GaussianFactorGraph &gfg, const FastMap &ordering, const std::vector &w) const { +vector SubgraphBuilder::kruskal(const GaussianFactorGraph &gfg, + const FastMap &ordering, + const vector &weights) const { const VariableIndex variableIndex(gfg); const size_t n = variableIndex.size(); - const vector idx = sort_idx(w) ; + const vector idx = sort_idx(weights) ; /* initialize buffer */ vector result; result.reserve(n-1); // container for acsendingly sorted edges - DSFVector D(n) ; + DSFVector dsf(n); size_t count = 0 ; double sum = 0.0 ; for (const size_t id: idx) { const GaussianFactor &gf = *gfg[id]; - if ( gf.keys().size() != 2 ) continue; - const size_t u = ordering.find(gf.keys()[0])->second, - u_root = D.find(u), - v = ordering.find(gf.keys()[1])->second, - v_root = D.find(v) ; - if ( u_root != v_root ) { - D.merge(u_root, v_root) ; + const auto keys = gf.keys(); + if ( keys.size() != 2 ) continue; + const size_t u = ordering.find(keys[0])->second, + v = ordering.find(keys[1])->second; + if ( dsf.find(u) != dsf.find(v) ) { + dsf.merge(u, v) ; result.push_back(id) ; - sum += w[id] ; + sum += weights[id] ; if ( ++count == n-1 ) break ; } } @@ -378,33 +390,40 @@ std::vector SubgraphBuilder::kruskal(const GaussianFactorGraph &gfg, con } /****************************************************************/ -std::vector SubgraphBuilder::sample(const std::vector &weights, const size_t t) const { - return uniqueSampler(weights, t); +vector SubgraphBuilder::sample(const vector &weights, const size_t t) const { + return UniqueSampler(weights, t); } /****************************************************************/ Subgraph::shared_ptr SubgraphBuilder::operator() (const GaussianFactorGraph &gfg) const { - const SubgraphBuilderParameters &p = parameters_; - const Ordering inverse_ordering = Ordering::Natural(gfg); + const auto &p = parameters_; + const auto inverse_ordering = Ordering::Natural(gfg); const FastMap forward_ordering = inverse_ordering.invert(); - const size_t n = inverse_ordering.size(), t = n * p.complexity_ ; + const size_t n = inverse_ordering.size(), m = gfg.size(); - vector w = weights(gfg); - const vector tree = buildTree(gfg, forward_ordering, w); + // Make sure the subgraph preconditioner does not include more than half of + // the edges beyond the spanning tree, or we might as well solve the whole thing. + size_t numExtraEdges = n * p.complexity_; + const size_t numRemaining = m - (n - 1); + numExtraEdges = std::min(numExtraEdges, numRemaining/2); - /* sanity check */ + // Calculate weights + vector weights = this->weights(gfg); + + // Build spanning tree. + const vector tree = buildTree(gfg, forward_ordering, weights); if ( tree.size() != n-1 ) { - throw runtime_error("SubgraphBuilder::operator() tree.size() != n-1 failed "); + throw std::runtime_error("SubgraphBuilder::operator() failure: tree.size() != n-1"); } - /* down weight the tree edges to zero */ + // Downweight the tree edges to zero. for ( const size_t id: tree ) { - w[id] = 0.0; + weights[id] = 0.0; } /* decide how many edges to augment */ - std::vector offTree = sample(w, t); + vector offTree = sample(weights, numExtraEdges); vector subgraph = unary(gfg); subgraph.insert(subgraph.end(), tree.begin(), tree.end()); @@ -450,7 +469,7 @@ SubgraphBuilder::Weights SubgraphBuilder::weights(const GaussianFactorGraph &gfg break; default: - throw invalid_argument("SubgraphBuilder::weights: undefined weight scheme "); + throw std::invalid_argument("SubgraphBuilder::weights: undefined weight scheme "); break; } } @@ -484,21 +503,20 @@ double SubgraphPreconditioner::error(const VectorValues& y) const { /* ************************************************************************* */ // gradient is y + inv(R1')*A2'*(A2*inv(R1)*y-b2bar), -VectorValues SubgraphPreconditioner::gradient(const VectorValues& y) const { +VectorValues SubgraphPreconditioner::gradient(const VectorValues &y) const { VectorValues x = Rc1()->backSubstitute(y); /* inv(R1)*y */ - Errors e = (*Ab2()*x - *b2bar()); /* (A2*inv(R1)*y-b2bar) */ + Errors e = (*Ab2() * x - *b2bar()); /* (A2*inv(R1)*y-b2bar) */ VectorValues v = VectorValues::Zero(x); - Ab2()->transposeMultiplyAdd(1.0, e, v); /* A2'*(A2*inv(R1)*y-b2bar) */ + Ab2()->transposeMultiplyAdd(1.0, e, v); /* A2'*(A2*inv(R1)*y-b2bar) */ return y + Rc1()->backSubstituteTranspose(v); } /* ************************************************************************* */ // Apply operator A, A*y = [I;A2*inv(R1)]*y = [y; A2*inv(R1)*y] Errors SubgraphPreconditioner::operator*(const VectorValues& y) const { - Errors e(y); VectorValues x = Rc1()->backSubstitute(y); /* x=inv(R1)*y */ - Errors e2 = *Ab2() * x; /* A2*x */ + Errors e2 = *Ab2() * x; /* A2*x */ e.splice(e.end(), e2); return e; } @@ -508,8 +526,10 @@ Errors SubgraphPreconditioner::operator*(const VectorValues& y) const { void SubgraphPreconditioner::multiplyInPlace(const VectorValues& y, Errors& e) const { Errors::iterator ei = e.begin(); - for (size_t i = 0; i < y.size(); ++i, ++ei) - *ei = y[i]; + for(const auto& key_value: y) { + *ei = key_value.second; + ++ei; + } // Add A2 contribution VectorValues x = Rc1()->backSubstitute(y); // x=inv(R1)*y @@ -522,8 +542,10 @@ VectorValues SubgraphPreconditioner::operator^(const Errors& e) const { Errors::const_iterator it = e.begin(); VectorValues y = zero(); - for (size_t i = 0; i < y.size(); ++i, ++it) - y[i] = *it; + for(auto& key_value: y) { + key_value.second = *it; + ++it; + } transposeMultiplyAdd2(1.0, it, e.end(), y); return y; } @@ -534,9 +556,10 @@ void SubgraphPreconditioner::transposeMultiplyAdd (double alpha, const Errors& e, VectorValues& y) const { Errors::const_iterator it = e.begin(); - for (size_t i = 0; i < y.size(); ++i, ++it) { + for(auto& key_value: y) { const Vector& ei = *it; - axpy(alpha, ei, y[i]); + axpy(alpha, ei, key_value.second); + ++it; } transposeMultiplyAdd2(alpha, it, e.end(), y); } @@ -563,47 +586,52 @@ void SubgraphPreconditioner::print(const std::string& s) const { } /*****************************************************************************/ -void SubgraphPreconditioner::solve(const Vector& y, Vector &x) const -{ - /* copy first */ - std::copy(y.data(), y.data() + y.rows(), x.data()); +void SubgraphPreconditioner::solve(const Vector &y, Vector &x) const { + assert(x.size() == y.size()); - /* in place back substitute */ - for (auto cg: boost::adaptors::reverse(*Rc1_)) { + /* back substitute */ + for (const auto &cg : boost::adaptors::reverse(*Rc1_)) { /* collect a subvector of x that consists of the parents of cg (S) */ - const Vector xParent = getSubvector(x, keyInfo_, KeyVector(cg->beginParents(), cg->endParents())); - const Vector rhsFrontal = getSubvector(x, keyInfo_, KeyVector(cg->beginFrontals(), cg->endFrontals())); + const KeyVector parentKeys(cg->beginParents(), cg->endParents()); + const KeyVector frontalKeys(cg->beginFrontals(), cg->endFrontals()); + const Vector xParent = getSubvector(x, keyInfo_, parentKeys); + const Vector rhsFrontal = getSubvector(y, keyInfo_, frontalKeys); /* compute the solution for the current pivot */ - const Vector solFrontal = cg->get_R().triangularView().solve(rhsFrontal - cg->get_S() * xParent); + const Vector solFrontal = cg->get_R().triangularView().solve( + rhsFrontal - cg->get_S() * xParent); /* assign subvector of sol to the frontal variables */ - setSubvector(solFrontal, keyInfo_, KeyVector(cg->beginFrontals(), cg->endFrontals()), x); + setSubvector(solFrontal, keyInfo_, frontalKeys, x); } } /*****************************************************************************/ -void SubgraphPreconditioner::transposeSolve(const Vector& y, Vector& x) const -{ +void SubgraphPreconditioner::transposeSolve(const Vector &y, Vector &x) const { /* copy first */ + assert(x.size() == y.size()); std::copy(y.data(), y.data() + y.rows(), x.data()); /* in place back substitute */ - for(const boost::shared_ptr & cg: *Rc1_) { - const Vector rhsFrontal = getSubvector(x, keyInfo_, KeyVector(cg->beginFrontals(), cg->endFrontals())); -// const Vector solFrontal = cg->get_R().triangularView().transpose().solve(rhsFrontal); - const Vector solFrontal = cg->get_R().transpose().triangularView().solve(rhsFrontal); + for (const auto &cg : *Rc1_) { + const Vector rhsFrontal = getSubvector( + x, keyInfo_, KeyVector(cg->beginFrontals(), cg->endFrontals())); + const Vector solFrontal = + cg->get_R().transpose().triangularView().solve( + rhsFrontal); // Check for indeterminant solution - if ( solFrontal.hasNaN()) throw IndeterminantLinearSystemException(cg->keys().front()); + if (solFrontal.hasNaN()) + throw IndeterminantLinearSystemException(cg->keys().front()); /* assign subvector of sol to the frontal variables */ - setSubvector(solFrontal, keyInfo_, KeyVector(cg->beginFrontals(), cg->endFrontals()), x); + setSubvector(solFrontal, keyInfo_, + KeyVector(cg->beginFrontals(), cg->endFrontals()), x); /* substract from parent variables */ - for (GaussianConditional::const_iterator it = cg->beginParents(); it != cg->endParents(); it++) { - KeyInfo::const_iterator it2 = keyInfo_.find(*it); - Eigen::Map rhsParent(x.data()+it2->second.colstart(), it2->second.dim(), 1); + for (auto it = cg->beginParents(); it != cg->endParents(); it++) { + const KeyInfoEntry &info = keyInfo_.find(*it)->second; + Eigen::Map rhsParent(x.data() + info.colstart(), info.dim(), 1); rhsParent -= Matrix(cg->getA(it)).transpose() * solFrontal; } } @@ -626,25 +654,24 @@ void SubgraphPreconditioner::build(const GaussianFactorGraph &gfg, const KeyInfo } /*****************************************************************************/ -Vector getSubvector(const Vector &src, const KeyInfo &keyInfo, const KeyVector &keys) { - +Vector getSubvector(const Vector &src, const KeyInfo &keyInfo, + const KeyVector &keys) { /* a cache of starting index and dim */ - typedef vector > Cache; - Cache cache; + vector > cache; /* figure out dimension by traversing the keys */ - size_t d = 0; - for ( const Key &key: keys ) { + size_t dim = 0; + for (const Key &key : keys) { const KeyInfoEntry &entry = keyInfo.find(key)->second; - cache.push_back(make_pair(entry.colstart(), entry.dim())); - d += entry.dim(); + cache.emplace_back(entry.colstart(), entry.dim()); + dim += entry.dim(); } /* use the cache to fill the result */ - Vector result = Vector::Zero(d, 1); + Vector result = Vector::Zero(dim); size_t idx = 0; - for ( const Cache::value_type &p: cache ) { - result.segment(idx, p.second) = src.segment(p.first, p.second) ; + for (const auto &p : cache) { + result.segment(idx, p.second) = src.segment(p.first, p.second); idx += p.second; } @@ -653,7 +680,6 @@ Vector getSubvector(const Vector &src, const KeyInfo &keyInfo, const KeyVector & /*****************************************************************************/ void setSubvector(const Vector &src, const KeyInfo &keyInfo, const KeyVector &keys, Vector &dst) { - /* use the cache */ size_t idx = 0; for ( const Key &key: keys ) { const KeyInfoEntry &entry = keyInfo.find(key)->second; @@ -663,10 +689,10 @@ void setSubvector(const Vector &src, const KeyInfo &keyInfo, const KeyVector &ke } /*****************************************************************************/ -boost::shared_ptr -buildFactorSubgraph(const GaussianFactorGraph &gfg, const Subgraph &subgraph, const bool clone) { - - GaussianFactorGraph::shared_ptr result(new GaussianFactorGraph()); +GaussianFactorGraph::shared_ptr buildFactorSubgraph( + const GaussianFactorGraph &gfg, const Subgraph &subgraph, + const bool clone) { + auto result = boost::make_shared(); result->reserve(subgraph.size()); for ( const SubgraphEdge &e: subgraph ) { const size_t idx = e.index(); diff --git a/gtsam/linear/SubgraphPreconditioner.h b/gtsam/linear/SubgraphPreconditioner.h index e440f32e4..a9746d104 100644 --- a/gtsam/linear/SubgraphPreconditioner.h +++ b/gtsam/linear/SubgraphPreconditioner.h @@ -129,11 +129,11 @@ namespace gtsam { enum AugmentationWeight { /* how to weigh the graph edges */ SKELETON = 0, /* use the same weights in building the skeleton */ -// STRETCH, /* stretch in the laplacian sense */ -// GENERALIZED_STRETCH /* the generalized stretch defined in jian2013iros */ +// STRETCH, /* stretch in the laplacian sense */ +// GENERALIZED_STRETCH /* the generalized stretch defined in jian2013iros */ } augmentationWeight_ ; - double complexity_; + double complexity_; /* factor multiplied with n, yields number of extra edges. */ SubgraphBuilderParameters() : skeleton_(KRUSKAL), skeletonWeight_(RANDOM), augmentationWeight_(SKELETON), complexity_(1.0) {} @@ -145,7 +145,7 @@ namespace gtsam { friend std::ostream& operator<<(std::ostream &os, const PreconditionerParameters &p); static Skeleton skeletonTranslator(const std::string &s); - static std::string skeletonTranslator(Skeleton w); + static std::string skeletonTranslator(Skeleton s); static SkeletonWeight skeletonWeightTranslator(const std::string &s); static std::string skeletonWeightTranslator(SkeletonWeight w); static AugmentationWeight augmentationWeightTranslator(const std::string &s); @@ -170,7 +170,7 @@ namespace gtsam { std::vector unary(const GaussianFactorGraph &gfg) const ; std::vector natural_chain(const GaussianFactorGraph &gfg) const ; std::vector bfs(const GaussianFactorGraph &gfg) const ; - std::vector kruskal(const GaussianFactorGraph &gfg, const FastMap &ordering, const std::vector &w) const ; + std::vector kruskal(const GaussianFactorGraph &gfg, const FastMap &ordering, const std::vector &weights) const ; std::vector sample(const std::vector &weights, const size_t t) const ; Weights weights(const GaussianFactorGraph &gfg) const; SubgraphBuilderParameters parameters_; @@ -249,8 +249,8 @@ namespace gtsam { /* A zero VectorValues with the structure of xbar */ VectorValues zero() const { - VectorValues V(VectorValues::Zero(*xbar_)); - return V ; + assert(xbar_); + return VectorValues::Zero(*xbar_); } /** @@ -285,15 +285,18 @@ namespace gtsam { /*****************************************************************************/ /* implement virtual functions of Preconditioner */ - /* Computation Interfaces for Vector */ - virtual void solve(const Vector& y, Vector &x) const; - virtual void transposeSolve(const Vector& y, Vector& x) const ; + /// implement x = R^{-1} y + void solve(const Vector& y, Vector &x) const override; - virtual void build( + /// implement x = R^{-T} y + void transposeSolve(const Vector& y, Vector& x) const override; + + /// build/factorize the preconditioner + void build( const GaussianFactorGraph &gfg, const KeyInfo &info, const std::map &lambda - ) ; + ) override; /*****************************************************************************/ }; @@ -310,9 +313,9 @@ namespace gtsam { /* sort the container and return permutation index with default comparator */ - template - std::vector sort_idx(const Container &src) - { + template + std::vector sort_idx(const Container &src) + { typedef typename Container::value_type T; const size_t n = src.size() ; std::vector > tmp; @@ -329,6 +332,6 @@ namespace gtsam { idx.push_back(tmp[i].first) ; } return idx; - } + } } // namespace gtsam diff --git a/gtsam/linear/tests/testGaussianBayesNet.cpp b/gtsam/linear/tests/testGaussianBayesNet.cpp index 74f07b92e..184933732 100644 --- a/gtsam/linear/tests/testGaussianBayesNet.cpp +++ b/gtsam/linear/tests/testGaussianBayesNet.cpp @@ -21,6 +21,8 @@ #include #include +#include +#include #include #include // for operator += using namespace boost::assign; @@ -28,13 +30,11 @@ using namespace boost::assign; // STL/C++ #include #include -#include -#include using namespace std; using namespace gtsam; -static const Key _x_=0, _y_=1; +static const Key _x_ = 11, _y_ = 22, _z_ = 33; static GaussianBayesNet smallBayesNet = list_of(GaussianConditional(_x_, Vector1::Constant(9), I_1x1, _y_, I_1x1))( @@ -42,9 +42,9 @@ static GaussianBayesNet smallBayesNet = static GaussianBayesNet noisyBayesNet = list_of(GaussianConditional(_x_, Vector1::Constant(9), I_1x1, _y_, I_1x1, - noiseModel::Diagonal::Sigmas(Vector1::Constant(2))))( + noiseModel::Isotropic::Sigma(1, 2.0)))( GaussianConditional(_y_, Vector1::Constant(5), I_1x1, - noiseModel::Diagonal::Sigmas(Vector1::Constant(3)))); + noiseModel::Isotropic::Sigma(1, 3.0))); /* ************************************************************************* */ TEST( GaussianBayesNet, Matrix ) @@ -140,11 +140,33 @@ TEST( GaussianBayesNet, optimize3 ) TEST(GaussianBayesNet, ordering) { Ordering expected; - expected += 0, 1; + expected += _x_, _y_; const auto actual = noisyBayesNet.ordering(); EXPECT(assert_equal(expected, actual)); } +/* ************************************************************************* */ +TEST( GaussianBayesNet, MatrixStress ) +{ + GaussianBayesNet bn; + using GC = GaussianConditional; + bn.emplace_shared(_x_, Vector2(1, 2), 1 * I_2x2, _y_, 2 * I_2x2, _z_, 3 * I_2x2); + bn.emplace_shared(_y_, Vector2(3, 4), 4 * I_2x2, _z_, 5 * I_2x2); + bn.emplace_shared(_z_, Vector2(5, 6), 6 * I_2x2); + + const VectorValues expected = bn.optimize(); + for (const auto keys : + {KeyVector({_x_, _y_, _z_}), KeyVector({_x_, _z_, _y_}), + KeyVector({_y_, _x_, _z_}), KeyVector({_y_, _z_, _x_}), + KeyVector({_z_, _x_, _y_}), KeyVector({_z_, _y_, _x_})}) { + const Ordering ordering(keys); + Matrix R; + Vector d; + boost::tie(R, d) = bn.matrix(ordering); + EXPECT(assert_equal(expected.vector(ordering), R.inverse() * d)); + } +} + /* ************************************************************************* */ TEST( GaussianBayesNet, backSubstituteTranspose ) { diff --git a/gtsam/slam/dataset.cpp b/gtsam/slam/dataset.cpp index c32a049e2..b66f22ced 100644 --- a/gtsam/slam/dataset.cpp +++ b/gtsam/slam/dataset.cpp @@ -66,6 +66,7 @@ string findExampleDataFile(const string& name) { namesToSearch.push_back(name + ".graph"); namesToSearch.push_back(name + ".txt"); namesToSearch.push_back(name + ".out"); + namesToSearch.push_back(name + ".xml"); // Find first name that exists for(const fs::path& root: rootsToSearch) { diff --git a/tests/testSubgraphPreconditioner.cpp b/tests/testSubgraphPreconditioner.cpp index accf9a65e..d9b07184b 100644 --- a/tests/testSubgraphPreconditioner.cpp +++ b/tests/testSubgraphPreconditioner.cpp @@ -1,6 +1,6 @@ /* ---------------------------------------------------------------------------- - * GTSAM Copyright 2010, Georgia Tech Research Corporation, + * GTSAM Copyright 2010, Georgia Tech Research Corporation, * Atlanta, Georgia 30332-0415 * All Rights Reserved * Authors: Frank Dellaert, et al. (see THANKS for the full author list) @@ -15,210 +15,302 @@ * @author Frank Dellaert **/ -#include - -#if 0 - #include + +#include +#include #include -#include +#include #include #include -#include -#include +#include +#include +#include -#include +#include + +#include #include +#include +#include +#include using namespace boost::assign; +#include + using namespace std; using namespace gtsam; using namespace example; // define keys // Create key for simulated planar graph -Symbol key(int x, int y) { - return symbol_shorthand::X(1000*x+y); -} +Symbol key(int x, int y) { return symbol_shorthand::X(1000 * x + y); } /* ************************************************************************* */ -TEST( SubgraphPreconditioner, planarOrdering ) { +TEST(SubgraphPreconditioner, planarOrdering) { // Check canonical ordering Ordering expected, ordering = planarOrdering(3); expected += key(3, 3), key(2, 3), key(1, 3), key(3, 2), key(2, 2), key(1, 2), key(3, 1), key(2, 1), key(1, 1); - CHECK(assert_equal(expected,ordering)); + EXPECT(assert_equal(expected, ordering)); } /* ************************************************************************* */ /** unnormalized error */ static double error(const GaussianFactorGraph& fg, const VectorValues& x) { double total_error = 0.; - for(const GaussianFactor::shared_ptr& factor: fg) + for (const GaussianFactor::shared_ptr& factor : fg) total_error += factor->error(x); return total_error; } /* ************************************************************************* */ -TEST( SubgraphPreconditioner, planarGraph ) - { +TEST(SubgraphPreconditioner, planarGraph) { // Check planar graph construction GaussianFactorGraph A; VectorValues xtrue; boost::tie(A, xtrue) = planarGraph(3); - LONGS_EQUAL(13,A.size()); - LONGS_EQUAL(9,xtrue.size()); - DOUBLES_EQUAL(0,error(A,xtrue),1e-9); // check zero error for xtrue + LONGS_EQUAL(13, A.size()); + LONGS_EQUAL(9, xtrue.size()); + DOUBLES_EQUAL(0, error(A, xtrue), 1e-9); // check zero error for xtrue // Check that xtrue is optimal - GaussianBayesNet::shared_ptr R1 = GaussianSequentialSolver(A).eliminate(); - VectorValues actual = optimize(*R1); - CHECK(assert_equal(xtrue,actual)); + GaussianBayesNet::shared_ptr R1 = A.eliminateSequential(); + VectorValues actual = R1->optimize(); + EXPECT(assert_equal(xtrue, actual)); } /* ************************************************************************* */ -TEST( SubgraphPreconditioner, splitOffPlanarTree ) -{ +TEST(SubgraphPreconditioner, splitOffPlanarTree) { // Build a planar graph GaussianFactorGraph A; VectorValues xtrue; boost::tie(A, xtrue) = planarGraph(3); // Get the spanning tree and constraints, and check their sizes - GaussianFactorGraph T, C; + GaussianFactorGraph::shared_ptr T, C; boost::tie(T, C) = splitOffPlanarTree(3, A); - LONGS_EQUAL(9,T.size()); - LONGS_EQUAL(4,C.size()); + LONGS_EQUAL(9, T->size()); + LONGS_EQUAL(4, C->size()); // Check that the tree can be solved to give the ground xtrue - GaussianBayesNet::shared_ptr R1 = GaussianSequentialSolver(T).eliminate(); - VectorValues xbar = optimize(*R1); - CHECK(assert_equal(xtrue,xbar)); + GaussianBayesNet::shared_ptr R1 = T->eliminateSequential(); + VectorValues xbar = R1->optimize(); + EXPECT(assert_equal(xtrue, xbar)); } /* ************************************************************************* */ - -TEST( SubgraphPreconditioner, system ) -{ +TEST(SubgraphPreconditioner, system) { // Build a planar graph GaussianFactorGraph Ab; VectorValues xtrue; size_t N = 3; - boost::tie(Ab, xtrue) = planarGraph(N); // A*x-b + boost::tie(Ab, xtrue) = planarGraph(N); // A*x-b - // Get the spanning tree and corresponding ordering - GaussianFactorGraph Ab1_, Ab2_; // A1*x-b1 and A2*x-b2 - boost::tie(Ab1_, Ab2_) = splitOffPlanarTree(N, Ab); - SubgraphPreconditioner::sharedFG Ab1(new GaussianFactorGraph(Ab1_)); - SubgraphPreconditioner::sharedFG Ab2(new GaussianFactorGraph(Ab2_)); + // Get the spanning tree and remaining graph + GaussianFactorGraph::shared_ptr Ab1, Ab2; // A1*x-b1 and A2*x-b2 + boost::tie(Ab1, Ab2) = splitOffPlanarTree(N, Ab); // Eliminate the spanning tree to build a prior - SubgraphPreconditioner::sharedBayesNet Rc1 = GaussianSequentialSolver(Ab1_).eliminate(); // R1*x-c1 - VectorValues xbar = optimize(*Rc1); // xbar = inv(R1)*c1 + const Ordering ord = planarOrdering(N); + auto Rc1 = Ab1->eliminateSequential(ord); // R1*x-c1 + VectorValues xbar = Rc1->optimize(); // xbar = inv(R1)*c1 // Create Subgraph-preconditioned system - VectorValues::shared_ptr xbarShared(new VectorValues(xbar)); // TODO: horrible - SubgraphPreconditioner system(Ab2, Rc1, xbarShared); + VectorValues::shared_ptr xbarShared( + new VectorValues(xbar)); // TODO: horrible + const SubgraphPreconditioner system(Ab2, Rc1, xbarShared); - // Create zero config - VectorValues zeros = VectorValues::Zero(xbar); + // Get corresponding matrices for tests. Add dummy factors to Ab2 to make + // sure it works with the ordering. + Ordering ordering = Rc1->ordering(); // not ord in general! + Ab2->add(key(1, 1), Z_2x2, Z_2x1); + Ab2->add(key(1, 2), Z_2x2, Z_2x1); + Ab2->add(key(1, 3), Z_2x2, Z_2x1); + Matrix A, A1, A2; + Vector b, b1, b2; + std::tie(A, b) = Ab.jacobian(ordering); + std::tie(A1, b1) = Ab1->jacobian(ordering); + std::tie(A2, b2) = Ab2->jacobian(ordering); + Matrix R1 = Rc1->matrix(ordering).first; + Matrix Abar(13 * 2, 9 * 2); + Abar.topRows(9 * 2) = Matrix::Identity(9 * 2, 9 * 2); + Abar.bottomRows(8) = A2.topRows(8) * R1.inverse(); + + // Helper function to vectorize in correct order, which is the order in which + // we eliminated the spanning tree. + auto vec = [ordering](const VectorValues& x) { return x.vector(ordering); }; // Set up y0 as all zeros - VectorValues y0 = zeros; + const VectorValues y0 = system.zero(); // y1 = perturbed y0 - VectorValues y1 = zeros; - y1[1] = Vector2(1.0, -1.0); + VectorValues y1 = system.zero(); + y1[key(3, 3)] = Vector2(1.0, -1.0); - // Check corresponding x values - VectorValues expected_x1 = xtrue, x1 = system.x(y1); - expected_x1[1] = Vector2(2.01, 2.99); - expected_x1[0] = Vector2(3.01, 2.99); - CHECK(assert_equal(xtrue, system.x(y0))); - CHECK(assert_equal(expected_x1,system.x(y1))); + // Check backSubstituteTranspose works with R1 + VectorValues actual = Rc1->backSubstituteTranspose(y1); + Vector expected = R1.transpose().inverse() * vec(y1); + EXPECT(assert_equal(expected, vec(actual))); + + // Check corresponding x values + // for y = 0, we get xbar: + EXPECT(assert_equal(xbar, system.x(y0))); + // for non-zero y, answer is x = xbar + inv(R1)*y + const Vector expected_x1 = vec(xbar) + R1.inverse() * vec(y1); + const VectorValues x1 = system.x(y1); + EXPECT(assert_equal(expected_x1, vec(x1))); // Check errors - DOUBLES_EQUAL(0,error(Ab,xtrue),1e-9); - DOUBLES_EQUAL(3,error(Ab,x1),1e-9); - DOUBLES_EQUAL(0,error(system,y0),1e-9); - DOUBLES_EQUAL(3,error(system,y1),1e-9); + DOUBLES_EQUAL(0, error(Ab, xbar), 1e-9); + DOUBLES_EQUAL(0, system.error(y0), 1e-9); + DOUBLES_EQUAL(2, error(Ab, x1), 1e-9); + DOUBLES_EQUAL(2, system.error(y1), 1e-9); - // Test gradient in x - VectorValues expected_gx0 = zeros; - VectorValues expected_gx1 = zeros; - CHECK(assert_equal(expected_gx0,gradient(Ab,xtrue))); - expected_gx1[2] = Vector2(-100., 100.); - expected_gx1[4] = Vector2(-100., 100.); - expected_gx1[1] = Vector2(200., -200.); - expected_gx1[3] = Vector2(-100., 100.); - expected_gx1[0] = Vector2(100., -100.); - CHECK(assert_equal(expected_gx1,gradient(Ab,x1))); + // Check that transposeMultiplyAdd <=> y += alpha * Abar' * e + // We check for e1 =[1;0] and e2=[0;1] corresponding to T and C + const double alpha = 0.5; + Errors e1, e2; + for (size_t i = 0; i < 13; i++) { + e1 += i < 9 ? Vector2(1, 1) : Vector2(0, 0); + e2 += i >= 9 ? Vector2(1, 1) : Vector2(0, 0); + } + Vector ee1(13 * 2), ee2(13 * 2); + ee1 << Vector::Ones(9 * 2), Vector::Zero(4 * 2); + ee2 << Vector::Zero(9 * 2), Vector::Ones(4 * 2); + + // Check transposeMultiplyAdd for e1 + VectorValues y = system.zero(); + system.transposeMultiplyAdd(alpha, e1, y); + Vector expected_y = alpha * Abar.transpose() * ee1; + EXPECT(assert_equal(expected_y, vec(y))); + + // Check transposeMultiplyAdd for e2 + y = system.zero(); + system.transposeMultiplyAdd(alpha, e2, y); + expected_y = alpha * Abar.transpose() * ee2; + EXPECT(assert_equal(expected_y, vec(y))); // Test gradient in y - VectorValues expected_gy0 = zeros; - VectorValues expected_gy1 = zeros; - expected_gy1[2] = Vector2(2., -2.); - expected_gy1[4] = Vector2(-2., 2.); - expected_gy1[1] = Vector2(3., -3.); - expected_gy1[3] = Vector2(-1., 1.); - expected_gy1[0] = Vector2(1., -1.); - CHECK(assert_equal(expected_gy0,gradient(system,y0))); - CHECK(assert_equal(expected_gy1,gradient(system,y1))); - - // Check it numerically for good measure - // TODO use boost::bind(&SubgraphPreconditioner::error,&system,_1) - // Vector numerical_g1 = numericalGradient (error, y1, 0.001); - // Vector expected_g1 = (Vector(18) << 0., 0., 0., 0., 2., -2., 0., 0., -2., 2., - // 3., -3., 0., 0., -1., 1., 1., -1.); - // CHECK(assert_equal(expected_g1,numerical_g1)); + auto g = system.gradient(y0); + Vector expected_g = Vector::Zero(18); + EXPECT(assert_equal(expected_g, vec(g))); } /* ************************************************************************* */ -TEST( SubgraphPreconditioner, conjugateGradients ) -{ +BOOST_CLASS_EXPORT_GUID(gtsam::JacobianFactor, "JacobianFactor"); + +// Read from XML file +static GaussianFactorGraph read(const string& name) { + auto inputFile = findExampleDataFile(name); + ifstream is(inputFile); + if (!is.is_open()) throw runtime_error("Cannot find file " + inputFile); + boost::archive::xml_iarchive in_archive(is); + GaussianFactorGraph Ab; + in_archive >> boost::serialization::make_nvp("graph", Ab); + return Ab; +} + +TEST(SubgraphSolver, Solves) { + // Create preconditioner + SubgraphPreconditioner system; + + // We test on three different graphs + const auto Ab1 = planarGraph(3).get<0>(); + const auto Ab2 = read("toy3D"); + const auto Ab3 = read("randomGrid3D"); + + // For all graphs, test solve and solveTranspose + for (const auto& Ab : {Ab1, Ab2, Ab3}) { + // Call build, a non-const method needed to make solve work :-( + KeyInfo keyInfo(Ab); + std::map lambda; + system.build(Ab, keyInfo, lambda); + + // Create a perturbed (non-zero) RHS + const auto xbar = system.Rc1()->optimize(); // merely for use in zero below + auto values_y = VectorValues::Zero(xbar); + auto it = values_y.begin(); + it->second.setConstant(100); + ++it; + it->second.setConstant(-100); + + // Solve the VectorValues way + auto values_x = system.Rc1()->backSubstitute(values_y); + + // Solve the matrix way, this really just checks BN::backSubstitute + // This only works with Rc1 ordering, not with keyInfo ! + // TODO(frank): why does this not work with an arbitrary ordering? + const auto ord = system.Rc1()->ordering(); + const Matrix R1 = system.Rc1()->matrix(ord).first; + auto ord_y = values_y.vector(ord); + auto vector_x = R1.inverse() * ord_y; + EXPECT(assert_equal(vector_x, values_x.vector(ord))); + + // Test that 'solve' does implement x = R^{-1} y + // We do this by asserting it gives same answer as backSubstitute + // Only works with keyInfo ordering: + const auto ordering = keyInfo.ordering(); + auto vector_y = values_y.vector(ordering); + const size_t N = R1.cols(); + Vector solve_x = Vector::Zero(N); + system.solve(vector_y, solve_x); + EXPECT(assert_equal(values_x.vector(ordering), solve_x)); + + // Test that transposeSolve does implement x = R^{-T} y + // We do this by asserting it gives same answer as backSubstituteTranspose + auto values_x2 = system.Rc1()->backSubstituteTranspose(values_y); + Vector solveT_x = Vector::Zero(N); + system.transposeSolve(vector_y, solveT_x); + EXPECT(assert_equal(values_x2.vector(ordering), solveT_x)); + } +} + +/* ************************************************************************* */ +TEST(SubgraphPreconditioner, conjugateGradients) { // Build a planar graph GaussianFactorGraph Ab; VectorValues xtrue; size_t N = 3; - boost::tie(Ab, xtrue) = planarGraph(N); // A*x-b + boost::tie(Ab, xtrue) = planarGraph(N); // A*x-b - // Get the spanning tree and corresponding ordering - GaussianFactorGraph Ab1_, Ab2_; // A1*x-b1 and A2*x-b2 - boost::tie(Ab1_, Ab2_) = splitOffPlanarTree(N, Ab); - SubgraphPreconditioner::sharedFG Ab1(new GaussianFactorGraph(Ab1_)); - SubgraphPreconditioner::sharedFG Ab2(new GaussianFactorGraph(Ab2_)); + // Get the spanning tree + GaussianFactorGraph::shared_ptr Ab1, Ab2; // A1*x-b1 and A2*x-b2 + boost::tie(Ab1, Ab2) = splitOffPlanarTree(N, Ab); // Eliminate the spanning tree to build a prior - Ordering ordering = planarOrdering(N); - SubgraphPreconditioner::sharedBayesNet Rc1 = GaussianSequentialSolver(Ab1_).eliminate(); // R1*x-c1 - VectorValues xbar = optimize(*Rc1); // xbar = inv(R1)*c1 + SubgraphPreconditioner::sharedBayesNet Rc1 = + Ab1->eliminateSequential(); // R1*x-c1 + VectorValues xbar = Rc1->optimize(); // xbar = inv(R1)*c1 // Create Subgraph-preconditioned system - VectorValues::shared_ptr xbarShared(new VectorValues(xbar)); // TODO: horrible + VectorValues::shared_ptr xbarShared( + new VectorValues(xbar)); // TODO: horrible SubgraphPreconditioner system(Ab2, Rc1, xbarShared); // Create zero config y0 and perturbed config y1 VectorValues y0 = VectorValues::Zero(xbar); VectorValues y1 = y0; - y1[1] = Vector2(1.0, -1.0); + y1[key(2, 2)] = Vector2(1.0, -1.0); VectorValues x1 = system.x(y1); // Solve for the remaining constraints using PCG ConjugateGradientParameters parameters; VectorValues actual = conjugateGradients(system, y1, parameters); - CHECK(assert_equal(y0,actual)); + EXPECT(assert_equal(y0,actual)); // Compare with non preconditioned version: VectorValues actual2 = conjugateGradientDescent(Ab, x1, parameters); - CHECK(assert_equal(xtrue,actual2,1e-4)); + EXPECT(assert_equal(xtrue, actual2, 1e-4)); } -#endif - /* ************************************************************************* */ -int main() { TestResult tr; return TestRegistry::runAllTests(tr); } +int main() { + TestResult tr; + return TestRegistry::runAllTests(tr); +} /* ************************************************************************* */