text:
struct TanYam7ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  c7::T2
  a21::T
  a31::T
  a32::T
  a41::T
  a43::T
  a51::T
  a53::T
  a54::T
  a61::T
  a63::T
  a64::T
  a65::T
  a71::T
  a73::T
  a74::T
  a75::T
  a76::T
  a81::T
  a83::T
  a84::T
  a85::T
  a86::T
  a87::T
  a91::T
  a93::T
  a94::T
  a95::T
  a96::T
  a97::T
  a98::T
  a101::T
  a103::T
  a104::T
  a105::T
  a106::T
  a107::T
  a108::T
  b1::T
  b4::T
  b5::T
  b6::T
  b7::T
  b8::T
  b9::T
  btilde1::T
  btilde4::T
  btilde5::T
  btilde6::T
  btilde7::T
  btilde8::T
  btilde9::T
  btilde10::T
end

"""
On the Optimization of Some Nine-Stage Seventh-order Runge-Kutta Method, by M. Tanaka, S. Muramatsu and S. Yamashita,
Information Processing Society of Japan, Vol. 33, No. 12 (1992) pages 1512-1526.
"""
function TanYam7ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  c1    =T2(0.07816646510113846)
  c2    =T2(0.1172496976517077)
  c3    =T2(0.17587454647756157)
  c4    =T2(0.4987401101913988)
  c5    =T2(0.772121690184484)
  c6    =T2(0.9911856696047768)
  c7    =T2(0.9995019582097662)
  a21   =T(0.07816646510113846)
  a31   =T(0.029312424412926925)
  a32   =T(0.08793727323878078)
  a41   =T(0.04396863661939039)
  a43   =T(0.13190590985817116)
  a51   =T(0.7361834837738382)
  a53   =T(-2.8337999624233303)
  a54   =T(2.5963565888408913)
  a61   =T(-12.062819391370866)
  a63   =T(48.20838100175243)
  a64   =T(-38.05863046463434)
  a65   =T(2.6851905444372632)
  a71   =T(105.21957276320198)
  a73   =T(-417.92888626241256)
  a74   =T(332.3155504499333)
  a75   =T(-19.827591183572938)
  a76   =T(1.2125399024549859)
  a81   =T(114.67755718631742)
  a83   =T(-455.5612169896097)
  a84   =T(362.24095553923144)
  a85   =T(-21.67190442182809)
  a86   =T(1.3189132007137807)
  a87   =T(-0.0048025566150346555)
  a91   =T(115.21334870553768)
  a93   =T(-457.69356568613233)
  a94   =T(363.93688218862735)
  a95   =T(-21.776682078900294)
  a96   =T(1.3250670887878468)
  a97   =T(-0.004518190986768983)
  a98   =T(-0.0005320269334859959)
  a101  =T(115.18928245800194)
  a103  =T(-457.598022271643)
  a104  =T(363.8610256312148)
  a105  =T(-21.77212754027556)
  a106  =T(1.3248804645074317)
  a107  =T(-0.0045057252106918315)
  a108  =T(-0.0005330165949429136)
  b1    =T(0.05126014249744686)
  b4    =T(0.27521638456212627)
  b5    =T(0.33696650340710543)
  b6    =T(0.18986072244906577)
  b7    =T(8.461099418514403)
  b8    =T(-130.15941672640542)
  b9    =T(121.84501355497527)
  # bhat1 =T(0.051002417590377186)
  # bhat4 =T(0.27613929504666546)
  # bhat5 =T(0.333788602069686)
  # bhat6 =T(0.20196531139081392)
  # bhat7 =T(5.755075459041811)
  # bhat8 =T(-85.61797108513936)
  # bhat10=T(80)
  btilde1 =T(-0.0002577249070696835)
  btilde4 =T(0.0009229104845391819)
  btilde5 =T(-0.0031779013374194105)
  btilde6 =T(0.01210458894174817)
  btilde7 =T(-2.706023959472591)
  btilde8 =T(44.541445641266066)
  btilde9 =T(-121.84501355497527)
  btilde10=T(80)

  TanYam7ConstantCache(c1,c2,c3,c4,c5,c6,c7,a21,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,a91,a93,a94,a95,a96,a97,a98,a101,a103,a104,a105,a106,a107,a108,b1,b4,b5,b6,b7,b8,b9,btilde1,btilde4,btilde5,btilde6,btilde7,btilde8,btilde9,btilde10)
end

"""
On the Optimization of Some Nine-Stage Seventh-order Runge-Kutta Method, by M. Tanaka, S. Muramatsu and S. Yamashita,
Information Processing Society of Japan, Vol. 33, No. 12 (1992) pages 1512-1526.
"""
function TanYam7ConstantCache(T::Type,T2::Type)
  c1    =T2(36259//463869)
  c2    =T2(36259//309246)
  c3    =T2(36259//206164)
  c4    =T2(76401//153188)
  c5    =T2(5164663400901569152//6688924124988083687)
  c6    =T2(3486//3517)
  c7    =T2(44151//44173)
  a21   =T(36259//463869)
  a31   =T(36259//1236984)
  a32   =T(36259//412328)
  a41   =T(36259//824656)
  a43   =T(108777//824656)
  a51   =T(17751533712975975593187//24112920357813127230992)
  a53   =T(-68331192803887602162951//24112920357813127230992)
  a54   =T(7825717455900471140481//3014115044726640903874)
  a61   =T(-parse(BigInt,"2425518501234340256175806929031336393991001205323654593685210322030691047097621496102266496")//parse(BigInt,"201073929944556265242953373967503382318096046559546854970564286270157897072030532387737241"))
  a63   =T(parse(BigInt,"126875939114499086848675646731069753055308007638565564293214808307459627250976287910912")//parse(BigInt,"2631823273838775215546306644775636213113650954300949659959480717139276934490785884841"))
  a64   =T(-parse(BigInt,"18238165682427587123600563411903599919711680222699338744428834349094610403849667513626245575680")//parse(BigInt,"479212348415302218688412787744011607018072627155280851781784515530195350673833210517995172867"))
  a65   =T(parse(BigInt,"74777425357290689294313120787550134201356775453356604582280658347816977407509825814840320")//parse(BigInt,"27848089034948481594251542168496834020714916243735255636715495143105044359669539013058107"))
  a71   =T(parse(BigInt,"42210784012026021620512889337138957588173072058924928398799062235")//parse(BigInt,"401168555464694196502745570125544252560955194769351196028554688"))
  a73   =T(-parse(BigInt,"53537582181289418572806048482253962781541488")//parse(BigInt,"128102133978061070595749084326726258918069"))
  a74   =T(parse(BigInt,"6373437319382536771018620806214785516542915567996760353063349991182871200304")//parse(BigInt,"19178871740288180724887392022898914045213833131528843480576173243533301485"))
  a75   =T(-parse(BigInt,"836513109281956728811652083904588515347012294160401579661057793958992")//parse(BigInt,"42189346226535262916910956145917457264775063492307360825161811325023"))
  a76   =T(parse(BigInt,"10038768138260655813133796321688310283082351149893792474426644227234755871856831386997923013888351")//parse(BigInt,"8279123943002224665888560854425725483235895533066047643118716510648226939201056966728652698557760"))
  a81   =T(parse(BigInt,"1454976871505621321312348899226731229297985195430097820532172928754404221419640982320963761")//parse(BigInt,"12687546780768188413911065021432924447284583965992535848754097389537051103097048673168256"))
  a83   =T(-parse(BigInt,"1452249436938195913836212549773886207822959770792")//parse(BigInt,"3187825000852340545619892931005470986913487349"))
  a84   =T(parse(BigInt,"3193785703967379485471835519262043520640585789136428552340853315619929163223926155626278646291801931779256")//parse(BigInt,"8816743814108800069900425523882492176796603795861854625575345408990649746129323017714575203134405597571"))
  a85   =T(-parse(BigInt,"314398569508916946629277462588835135011587938712337655816458752800894863689255534896547161759213480")//parse(BigInt,"14507196201560052990013371105817112064769849230048646555812475120383456376679192045076337148816813"))
  a86   =T(parse(BigInt,"5021633516852870452803558794670341128133410978274753232000155240629688617274518068065484524425884625107263111090060721584249881611265924113")//parse(BigInt,"3807402575192378287101053794016079417728266285278436439472658972755893033722804748992796724254152818232996309281540415603729279478920107136"))
  a87   =T(-parse(BigInt,"894451839895008223904010765658125850176064186717638397881061173697811879745")//parse(BigInt,"186244934020117483847289332768639722211239803963523669807238114327710091115676"))
  a91   =T(parse(BigInt,"152015786770038627019906826956584678402371493198250158080970494807155603994339")//parse(BigInt,"1319428594672311986480108760138089275639618425553698631283119461253421932416"))
  a93   =T(-parse(BigInt,"19887569115365707672105043997835466942389220328")//parse(BigInt,"43451712251082409470704235239058276887205131"))
  a94   =T(parse(BigInt,"6298831527954572673520838478029639446424615570453903300371170696118960335541193275024146681623960")//parse(BigInt,"17307483347318198085207889427954666589398911583434527253470846782562794571553580157056644256313"))
  a95   =T(-parse(BigInt,"16267621644623777942279856217571823792451732234540266142050307930357537283432611648312520")//parse(BigInt,"747020211145282116967827947968990352912884924402891384654470989583659988117513448655559"))
  a96   =T(parse(BigInt,"491920517345271821393960134665582163547632868347911487496995665146055538579545277983570189994492481977206720065882583432234119698425636137169515")//parse(BigInt,"371241970695441505578374965290296000309261530083026613438333515399198575394818137422626328203755084156959422247928840402063855870066548878130304"))
  a97   =T(-parse(BigInt,"17535891839112183607157943692398769696531153719141528498448224128785868799210475")//parse(BigInt,"3881175428498724209649715816699297677268154716152409333146177577349474565697791732"))
  a98   =T(-parse(BigInt,"31140449219386755112730831706895080247696102690585728771850210691242594436100540310")//parse(BigInt,"58531715707220748822628340615174217489020037018063169180406742622693159384762890406389"))
  a101  =T(parse(BigInt,"24861126512935523838485032295435745281790804119672244200744512677831357181363")//parse(BigInt,"215828469302253893975010055544246846578750854407392771457340001283636121600"))
  a103  =T(-76626859319946149305867456329803//167454524692981091214376557800)
  a104  =T(parse(BigInt,"257532657386915224604779230484778835596042580268896440943054087972106955277512448850995064336363")//parse(BigInt,"707777528357579864776572552477247532276956780876653359042572831013312547307465249178438602200"))
  a105  =T(-parse(BigInt,"103092665221253777021612043042409780416654274677686197534469014507504059634284484983141143")//parse(BigInt,"4735075386204034224907103653335170134874540866215348781137359896717512695961598377363000"))
  a106  =T(parse(BigInt,"1318945254307068672853031172410281620677291556423152759282406612372948205789241763483098989903852936890735513699395545618802215742952753372919")//parse(BigInt,"995520191927224509158660659519643916330847017611189618002256023928790665495276022949114110343406997764203331763292012060684160018593393766400"))
  a107  =T(-parse(BigInt,"2175691361381933486174620849991740173349017185199505364607841")//parse(BigInt,"482872625303278742130341621563226511344221688759361797916327450"))
  a108  =T(-parse(BigInt,"11327601987184122343710458559595782081610122892585097")//parse(BigInt,"21251874884678431935286330856983429378055579208005268000"))
  b1    =T(parse(BigInt,"677260699094873524061210073954310211")//parse(BigInt,"13212228177645157882237395248920447488"))
  b4    =T(parse(BigInt,"5627843976805934592544586970647029617399366281651959837492864")//parse(BigInt,"20448796992082885248862284273169726631726393791864145954479875"))
  b5    =T(parse(BigInt,"1359735671458057021603668186882234273947181034928034734244224")//parse(BigInt,"4035225037829041960922838374222759264846456609494840689395475"))
  b6    =T(parse(BigInt,"3575764371063841994042920363615768888383369782579963896064642431626191680598750790399139608006651160426580137040859330533720256407")//parse(BigInt,"18833618269956378326078572170759846509476617594300797062242096554507068838086062412372695473217373611870290738365243380652826304000"))
  b7    =T(parse(BigInt,"14322850798205614664394883796805489119964080948503151")//parse(BigInt,"1692788382425178679633337406927131793062126418747780"))
  b8    =T(-parse(BigInt,"16735096417960349589058935251250023138290806176584545269411")//parse(BigInt,"128573843052304513208482301684749747737236254208431871400"))
  b9    =T(33050288141543277444692395096256051//271248590133163812341791503489000)
  # bhat1 =T(parse(BigInt,"962650826879437817605721930727384851")//parse(BigInt,"18874611682350225546053421784172067840"))
  # bhat4 =T(parse(BigInt,"99703652969826806275610089806158069716600653757297413344")//parse(BigInt,"361062893830367886445877712954352019629670588714825566425"))
  # bhat5 =T(parse(BigInt,"17540887447270394964911517553576959050951784592644178144")//parse(BigInt,"52550888012671962193116521992300249584518949945887203425"))
  # bhat6 =T(parse(BigInt,"101855668513773837712956593596043266148479443244790887636953159551191054134940671472736229702711787350735239179")//parse(BigInt,"504322587935299170723833764883183242017770187561624249681119708768991642691172146267201689787026963930014131200"))
  # bhat7 =T(parse(BigInt,"179578338747395946570172802104016572846366090083599")//parse(BigInt,"31203472487100067827342625012481692038011546889360"))
  # bhat8 =T(-parse(BigInt,"500374162579884236288722085953024481890963958534161489781")//parse(BigInt,"5844265593286568782203740985670443078965284282201448700"))
  # bhat10=T(80)
  btilde1 =T(parse(BigInt,"-11350400930890172457349074817136051")//parse(BigInt,"44040760592150526274124650829734824960"))
  btilde4 =T(parse(BigInt,"18872409140206580874590465524732661000311743892579167244576")//parse(BigInt,"20448796992082885248862284273169726631726393791864145954479875"))
  btilde5 =T(parse(BigInt,"-4274515681501734477669162831906773100582117137555409033632")//parse(BigInt,"1345075012609680653640946124740919754948818869831613563131825"))
  btilde6 =T(parse(BigInt,"151982138295746861476872192192436808638630079892291260212523545728864842939698890632387350810275352451114165347163409431707619557")//parse(BigInt,"12555745513304252217385714780506564339651078396200531374828064369671379225390708274915130315478249074580193825576828920435217536000"))
  btilde7 =T(parse(BigInt,"-6107634561545846083950679043550120057398294081957207")//parse(BigInt,"2257051176566904906177783209236175724082835224997040"))
  btilde8 =T(parse(BigInt,"1908954947067632130235683120094494845563199696277664164743")//parse(BigInt,"42857947684101504402827433894916582579078751402810623800"))
  btilde9 =T(-33050288141543277444692395096256051//271248590133163812341791503489000)
  btilde10=T(80)

  TanYam7ConstantCache(c1,c2,c3,c4,c5,c6,c7,a21,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,a91,a93,a94,a95,a96,a97,a98,a101,a103,a104,a105,a106,a107,a108,b1,b4,b5,b6,b7,b8,b9,btilde1,btilde4,btilde5,btilde6,btilde7,btilde8,btilde9,btilde10)
end

struct TsitPap8ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  c7::T2
  c8::T2
  c9::T2
  c10::T2
  a0201::T
  a0301::T
  a0302::T
  a0401::T
  a0403::T
  a0501::T
  a0503::T
  a0504::T
  a0601::T
  a0604::T
  a0605::T
  a0701::T
  a0704::T
  a0705::T
  a0706::T
  a0801::T
  a0804::T
  a0805::T
  a0806::T
  a0807::T
  a0901::T
  a0904::T
  a0905::T
  a0906::T
  a0907::T
  a0908::T
  a1001::T
  a1004::T
  a1005::T
  a1006::T
  a1007::T
  a1008::T
  a1009::T
  a1101::T
  a1104::T
  a1105::T
  a1106::T
  a1107::T
  a1108::T
  a1109::T
  a1110::T
  a1201::T
  a1204::T
  a1205::T
  a1206::T
  a1207::T
  a1208::T
  a1209::T
  a1210::T
  a1211::T
  a1301::T
  a1304::T
  a1305::T
  a1306::T
  a1307::T
  a1308::T
  a1309::T
  a1310::T
  b1::T
  b6::T
  b7::T
  b8::T
  b9::T
  b10::T
  b11::T
  b12::T
  btilde1::T
  btilde6::T
  btilde7::T
  btilde8::T
  btilde9::T
  btilde10::T
  btilde11::T
  btilde12::T
  btilde13::T
end

"""
Cheap Error Estimation for Runge-Kutta methods, by Ch. Tsitouras and S.N. Papakostas,
 Siam Journal on Scientific Computing, Vol. 20, Issue 6, Nov 1999.
"""
function TsitPap8ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}

  c1    =T2(0.06338028169014084)
  c2    =T2(0.1027879458763643)
  c3    =T2(0.15418191881454646)
  c4    =T2(0.3875968992248062)
  c5    =T2(0.4657534246575342)
  c6    =T2(0.1554054054054054)
  c7    =T2(1.0070921985815602)
  c8    =T2(0.876141078561489)
  c9    =T2(0.9120879120879121)
  c10   =T2(0.959731543624161)
  a0201 =T(0.06338028169014084)
  a0301 =T(0.019438980427336498)
  a0302 =T(0.08334896544902781)
  a0401 =T(0.038545479703636615)
  a0403 =T(0.11563643911090984)
  a0501 =T(0.39436557770112496)
  a0503 =T(-1.4818719321673373)
  a0504 =T(1.4751032536910185)
  a0601 =T(0.045994489107698204)
  a0604 =T(0.23235070626395474)
  a0605 =T(0.18740822928588133)
  a0701 =T(0.06005228953244051)
  a0704 =T(0.11220383194636775)
  a0705 =T(-0.03357232951906142)
  a0706 =T(0.016721613445658576)
  a0801 =T(-1.5733292732086857)
  a0804 =T(-1.3167087730223663)
  a0805 =T(-11.723515296181773)
  a0806 =T(9.107825028173872)
  a0807 =T(6.512820512820513)
  a0901 =T(-0.48107625624391254)
  a0904 =T(-6.6506103607463904)
  a0905 =T(-4.530206099782572)
  a0906 =T(3.894414525020157)
  a0907 =T(8.634217645525526)
  a0908 =T(0.009401624788681498)
  a1001 =T(-0.7754121446230569)
  a1004 =T(-7.996604718235832)
  a1005 =T(-6.726558607230182)
  a1006 =T(5.532184454327406)
  a1007 =T(10.89757332024991)
  a1008 =T(0.020091650280045396)
  a1009 =T(-0.039186042680376856)
  a1101 =T(-1.1896363245449992)
  a1104 =T(-7.128368483301214)
  a1105 =T(-9.53722789710108)
  a1106 =T(7.574470108980868)
  a1107 =T(11.267486382070919)
  a1108 =T(0.051009801223058315)
  a1109 =T(0.08019413469508256)
  a1110 =T(-0.15819617839847347)
  a1201 =T(-0.39200039047127266)
  a1204 =T(3.916659042493856)
  a1205 =T(-2.8017459289080557)
  a1206 =T(2.441204566481742)
  a1207 =T(-2.4183655778824718)
  a1208 =T(-0.33943326290032927)
  a1209 =T(0.19496450383103364)
  a1210 =T(-0.19437176762508154)
  a1211 =T(0.5930888149805791)
  a1301 =T(-1.4847063081291894)
  a1304 =T(-2.390723588981498)
  a1305 =T(-11.184306772840532)
  a1306 =T(8.720804667459817)
  a1307 =T(7.33673830753461)
  a1308 =T(0.01289874999394761)
  a1309 =T(0.042583289842657704)
  a1310 =T(-0.05328834487981156)
  b1    =T(0.04441161093250152)
  b6    =T(0.35395063113733116)
  b7    =T(0.2485219684184965)
  b8    =T(-0.3326913171720666)
  b9    =T(1.921248828652836)
  b10   =T(-2.7317783000882523)
  b11   =T(1.4012004409899175)
  b12   =T(0.0951361371292365)
  # bhat1 =T(0.044484201850329544)
  # bhat6 =T(0.35502352274458154)
  # bhat7 =T(0.24825530158391707)
  # bhat8 =T(-2.4242252962684616)
  # bhat9 =T(1.5999301534099692)
  # bhat10=T(-1.8107646286929684)
  # bhat13=T(2.9872967453726327)
  btilde1 =T(-7.259091782802626e-5)
  btilde6 =T(-0.0010728916072503584)
  btilde7 =T(0.0002666668345794398)
  btilde8 =T(2.091533979096395)
  btilde9 =T(0.3213186752428666)
  btilde10=T(-0.921013671395284)
  btilde11=T(1.4012004409899175)
  btilde12=T(0.0951361371292365)
  btilde13=T(-2.9872967453726327)

  TsitPap8ConstantCache(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1310,b1,b6,b7,b8,b9,b10,b11,b12,btilde1,btilde6,btilde7,btilde8,btilde9,btilde10,btilde11,btilde12,btilde13)
end

"""
Cheap Error Estimation for Runge-Kutta methods, by Ch. Tsitouras and S.N. Papakostas,
 Siam Journal on Scientific Computing, Vol. 20, Issue 6, Nov 1999.
"""
function TsitPap8ConstantCache(T::Type,T2::Type)

  c1    =T2(9//142)
  c2    =T2(24514//238491)
  c3    =T2(12257//79497)
  c4    =T2(50//129)
  c5    =T2(34//73)
  c6    =T2(23//148)
  c7    =T2(142//141)
  c8    =T2(29104120198//33218531707)
  c9    =T2(83//91)
  c10   =T2(143//149)
  a0201 =T(9//142)
  a0301 =T(9950845450//511901613729)
  a0302 =T(42666469916//511901613729)
  a0401 =T(12257//317988)
  a0403 =T(12257//105996)
  a0501 =T(14131686489425//35833975601529)
  a0503 =T(-17700454220625//11944658533843)
  a0504 =T(17619604667500//11944658533843)
  a0601 =T(4418011//96055225)
  a0604 =T(9757757988//41995818067)
  a0605 =T(19921512441//106300094275)
  a0701 =T(480285555889619//7997789253816320)
  a0704 =T(52162106556696460538643//464887033284472899365888)
  a0705 =T(-16450769949384489//490009784398315520)
  a0706 =T(864813381633887//51718297665732608)
  a0801 =T(-7001048352088587137143//4449830351030385148800)
  a0804 =T(-2522765548294599044197935933//1915963195493758447677901792)
  a0805 =T(-318556182235222634647116091//27172411532484037214054400)
  a0806 =T(1205563885850790193966807//132365727505912221222912)
  a0807 =T(254//39)
  a0901 =T(-parse(BigInt,"20629396399716689122801179264428539394462855874226604463554767070845753369")//parse(BigInt,"42881759662770155956657513470012114488076017981128105307375594963152353200"))
  a0904 =T(-parse(BigInt,"3315443074343659404779422149387397712986453181141168247590906370819301077749322753")//parse(BigInt,"498517112641608872807821838566847987729514312398278633788185835348997643218553476"))
  a0905 =T(-parse(BigInt,"273749409411654060286948141164828452109898379203526945684314474186724062841643")//parse(BigInt,"60427583951377552503967840653825589117167443021628367691773162913607984889600"))
  a0906 =T(parse(BigInt,"16656372518874512738268060504309924900437672263609245028809229865738327731797537")//parse(BigInt,"4276990138533930522782076385771016009097627930550149628282203007913538394549504"))
  a0907 =T(parse(BigInt,"42008080033354305590804322944084805264441066760038302359736803632")//parse(BigInt,"4865302423216534910074823287811605599629170295030631799935804001"))
  a0908 =T(parse(BigInt,"668459780930716338000066627236927417191947396177093524377824")//parse(BigInt,"71100452948884643779799087713322002499799983434685642170251833"))
  a1001 =T(-parse(BigInt,"1793603946322260900828212460706877142477132870159527")//parse(BigInt,"2313097568511990753781649719556084665131024900300800"))
  a1004 =T(-parse(BigInt,"14776874123722838192315406145167687512425345723701")//parse(BigInt,"1847893530366076701102014146927696206329050105856"))
  a1005 =T(-parse(BigInt,"19587020919884661714856757105130246995757906603")//parse(BigInt,"2911893296942532038566182868170393629789388800"))
  a1006 =T(parse(BigInt,"6364380863259071677112259236455506477417699780613300364807")//parse(BigInt,"1150428174584942133406579091549443438814988349188394057728"))
  a1007 =T(parse(BigInt,"27725164402569748756040320433848245155581006369")//parse(BigInt,"2544159473655547770881695354241256106302348256"))
  a1008 =T(parse(BigInt,"10744247163960019876833255044784609639")//parse(BigInt,"534761804739901348825491947768304503296"))
  a1009 =T(-parse(BigInt,"50977737930792808232204417497248979399878217280011103197862899")//parse(BigInt,"1300915694564913675613280314081837358644964393191337994183389184"))
  a1101 =T(-parse(BigInt,"3587625717068952487214493441966897048737050755812600710793")//parse(BigInt,"3015733164033229624772006086429467685046639983706974412800"))
  a1104 =T(-parse(BigInt,"5453011711267804731211501837262816944661201619903")//parse(BigInt,"764973320899716397072448644710733101329290196992"))
  a1105 =T(-parse(BigInt,"884348836774584715070440485633026464653487653")//parse(BigInt,"92725983515963799584249219499787659014963200"))
  a1106 =T(parse(BigInt,"26823469063654084387375587616552322383082061411417182757389742951")//parse(BigInt,"3541299744763681675473620647087123057228744296123642301654630400"))
  a1107 =T(parse(BigInt,"142363419491686507162007051071007722765323162710521029")//parse(BigInt,"12634887202368261565807449771335728082273587905775840"))
  a1108 =T(parse(BigInt,"64747617454909275289531520412519442831235890581")//parse(BigInt,"1269317188117975960996670628974453443777854830080"))
  a1109 =T(parse(BigInt,"112633808253272720979874303367503891597499261046700689572459050065039333987335667")//parse(BigInt,"1404514291245034532812181377119034501014039830518218465064579291611563374608650240"))
  a1110 =T(-parse(BigInt,"10612202518573994431153697720606405883")//parse(BigInt,"67082546658259846778754594976831647575"))
  a1201 =T(-parse(BigInt,"7534081165544982478296202335922049210803875045423")//parse(BigInt,"19219575665440848756074598051658416387002479820800"))
  a1204 =T(parse(BigInt,"237696087452786717802270375283034262859273455")//parse(BigInt,"60688480889936839261131818793519097177034752"))
  a1205 =T(-parse(BigInt,"20610578209826329263318986584876108069323")//parse(BigInt,"7356333776438834884293014575840932659200"))
  a1206 =T(parse(BigInt,"51260471529841028040709654458903254781320136131844164563")//parse(BigInt,"20998023776318546907382302106788168158765514131585630208"))
  a1207 =T(-parse(BigInt,"3077214437173472971196810795615384000211457151011")//parse(BigInt,"1272435592582280820597059432060116893200684680384"))
  a1208 =T(-parse(BigInt,"1539218116260541896259682954580256454049")//parse(BigInt,"4534670830750360983422999946409310393344"))
  a1209 =T(parse(BigInt,"241886539350268429372116296787271276553970618941104594460614948326132797451456131")//parse(BigInt,"1240669632662465892528916120041123051054026283188324780101555345343662316090597376"))
  a1210 =T(-80556486832245966191717452425924975//414445409518676597565032008051106461)
  a1211 =T(parse(BigInt,"2944781680874500347594142792814463350")//parse(BigInt,"4965161383073676983610218096030654529"))
  a1301 =T(-parse(BigInt,"7757739937862944832927743694336116203639371542761")//parse(BigInt,"5225100678421794325654850845451340473577260032000"))
  a1304 =T(-parse(BigInt,"433889546009521405913741133329446636837810749")//parse(BigInt,"181488796115643001621310691169322233346938880"))
  a1305 =T(-parse(BigInt,"246044720162308748108107126829066792329071")//parse(BigInt,"21999103311417807170100413255475150848000"))
  a1306 =T(parse(BigInt,"2140331235425829844389060818616719848637810765257179167")//parse(BigInt,"245428182036011897638028252815369258646052549878743040"))
  a1307 =T(parse(BigInt,"1573990926219809229258666534611598771063240529")//parse(BigInt,"214535514317495538101650508596881057760818080"))
  a1308 =T(parse(BigInt,"62408280667309375445301959066100433563")//parse(BigInt,"4838320046251983849512355559327451645440"))
  a1309 =T(parse(BigInt,"1145609822249493677618113725359506642998153205603226883141207089968379")//parse(BigInt,"26902802166822768476367427840835743139559294105398677975568538466119680"))
  a1310 =T(-408950356875874683139089678053832//7674292714443204070455739595109785)
  b1    =T(55038446513529253801//1239280570055853383520)
  b6    =T(2335496795323782464411846394611//6598368783294020109895379936256)
  b7    =T(7636073376527143565375240869888//30725949199261296642046754748645)
  b8    =T(-4237087214169934312729607487//12735791394214625116604076160)
  b9    =T(parse(BigInt,"408505291291133241760995514121984335914363927884426780078325258228227984174126699")//parse(BigInt,"212624874612193697466655159405635202123821166475348052734199905546455860773575680"))
  b10   =T(-1108225296327029096435947//405679075893979729103310)
  b11   =T(2460988291206213825688467985//1756342789520947764222671739)
  b12   =T(4808707937311//50545545388065)
  # bhat1 =T(715953338020208413//16094552857868225760)
  # bhat6 =T(62284335162928966987066121//175437206755843240272669696)
  # bhat7 =T(184309146777472302831695872//742417767522160213324368465)
  # bhat8 =T(-509771598215811385123057257//210282269969085698264101760)
  # bhat9 =T(parse(BigInt,"2701602489646143640362891402924962500766379885231830470264478480621389")//parse(BigInt,"1688575269294196295712420798364925687791337062814161130291144634516480"))
  # bhat10=T(-7218534073012286740367561//3986456306153224954098780)
  # bhat13=T(5752173075461//1925544586212)
  btilde1 =T(-1124506425334925//15491007125698167294)
  btilde6 =T(-34035261967014004512968665//31722926842759712066804711232)
  btilde7 =T(24580774837247048845194838016//92177847597783889926140264245935)
  btilde8 =T(7658235379858959628545472335//3661540025836704721023671896)
  btilde9 =T(parse(BigInt,"1510930663253486646384296195586886863633653147150681174690536256975353069486325")//parse(BigInt,"4702280880846591386281796794547701585430660412435581935467882526508158459415616"))
  btilde10=T(-237184116991804346989794715//257525077377498332034781188)
  btilde11=T(2460988291206213825688467985//1756342789520947764222671739)
  btilde12=T(4808707937311//50545545388065)
  btilde13=T(-5752173075461//1925544586212)

  TsitPap8ConstantCache(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1310,b1,b6,b7,b8,b9,b10,b11,b12,btilde1,btilde6,btilde7,btilde8,btilde9,btilde10,btilde11,btilde12,btilde13)
end

struct DP8ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  c7::T2
  c8::T2
  c9::T2
  c10::T2
  c11::T2
  c6::T2
  c5::T2
  c4::T2
  c3::T2
  c2::T2
  b1::T
  b6::T
  b7::T
  b8::T
  b9::T
  b10::T
  b11::T
  b12::T
  btilde1::T
  btilde6::T
  btilde7::T
  btilde8::T
  btilde9::T
  btilde10::T
  btilde11::T
  btilde12::T
  er1::T
  er6::T
  er7::T
  er8::T
  er9::T
  er10::T
  er11::T
  er12::T
  a0201::T
  a0301::T
  a0302::T
  a0401::T
  a0403::T
  a0501::T
  a0503::T
  a0504::T
  a0601::T
  a0604::T
  a0605::T
  a0701::T
  a0704::T
  a0705::T
  a0706::T
  a0801::T
  a0804::T
  a0805::T
  a0806::T
  a0807::T
  a0901::T
  a0904::T
  a0905::T
  a0906::T
  a0907::T
  a0908::T
  a1001::T
  a1004::T
  a1005::T
  a1006::T
  a1007::T
  a1008::T
  a1009::T
  a1101::T
  a1104::T
  a1105::T
  a1106::T
  a1107::T
  a1108::T
  a1109::T
  a1110::T
  a1201::T
  a1204::T
  a1205::T
  a1206::T
  a1207::T
  a1208::T
  a1209::T
  a1210::T
  a1211::T
  c14::T2
  c15::T2
  c16::T2
  a1401::T
  a1407::T
  a1408::T
  a1409::T
  a1410::T
  a1411::T
  a1412::T
  a1413::T
  a1501::T
  a1506::T
  a1507::T
  a1508::T
  a1511::T
  a1512::T
  a1513::T
  a1514::T
  a1601::T
  a1606::T
  a1607::T
  a1608::T
  a1609::T
  a1613::T
  a1614::T
  a1615::T
  d401::T
  d406::T
  d407::T
  d408::T
  d409::T
  d410::T
  d411::T
  d412::T
  d413::T
  d414::T
  d415::T
  d416::T
  d501::T
  d506::T
  d507::T
  d508::T
  d509::T
  d510::T
  d511::T
  d512::T
  d513::T
  d514::T
  d515::T
  d516::T
  d601::T
  d606::T
  d607::T
  d608::T
  d609::T
  d610::T
  d611::T
  d612::T
  d613::T
  d614::T
  d615::T
  d616::T
  d701::T
  d706::T
  d707::T
  d708::T
  d709::T
  d710::T
  d711::T
  d712::T
  d713::T
  d714::T
  d715::T
  d716::T
end

function DP8ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  c7     = T2(0.25)
  c8     = T2(0.3076923076923077)
  c9     = T2(0.6512820512820513)
  c10    = T2(0.6)
  c11    = T2(0.8571428571428571)
  c6     = T2(0.3333333333333333)
  c5     = T2(0.2816496580927726)
  c4     = T2(0.1183503419072274)
  c3     = T2(0.0789002279381516)
  c2     = T2(0.05260015195876773)
  b1     =  T(0.054293734116568765)
  b6     =  T(4.450312892752409)
  b7     =  T(1.8915178993145003)
  b8     =  T(-5.801203960010585)
  b9     =  T(0.3111643669578199)
  b10    = T(-0.1521609496625161)
  b11    = T(0.20136540080403034)
  b12    = T(0.04471061572777259)
  # bhh1   = T(0.2440944881889764)
  # bhh2   = T(0.7338466882816118)
  # bhh3   = T(0.022058823529411766)
  btilde1 = T(-0.18980075407240762)
  btilde6 = T(4.450312892752409)
  btilde7 = T(1.8915178993145003)
  btilde8 = T(-5.801203960010585)
  btilde9 = T(-0.42268232132379197)
  btilde10 = T(-0.1521609496625161)
  btilde11 = T(0.20136540080403034)
  btilde12 = T(0.022651792198360825)
  er1    = T(0.01312004499419488)
  er6    = T(-1.2251564463762044)
  er7    = T(-0.4957589496572502)
  er8    = T(1.6643771824549864)
  er9    = T(-0.35032884874997366)
  er10   = T(0.3341791187130175)
  er11   = T(0.08192320648511571)
  er12   = T(-0.022355307863886294)
  a0201  =   T(0.05260015195876773)
  a0301  =   T(0.0197250569845379)
  a0302  =   T(0.0591751709536137)
  a0401  =   T(0.02958758547680685)
  a0403  =   T(0.08876275643042054)
  a0501  =   T(0.2413651341592667)
  a0503  =   T(-0.8845494793282861)
  a0504  =   T(0.924834003261792)
  a0601  =   T(0.037037037037037035)
  a0604  =   T(0.17082860872947386)
  a0605  =   T(0.12546768756682242)
  a0701  =   T(0.037109375)
  a0704  =   T(0.17025221101954405)
  a0705  =   T(0.06021653898045596)
  a0706  =   T(-0.017578125)
  a0801  =   T(0.03709200011850479)
  a0804  =   T(0.17038392571223998)
  a0805  =   T(0.10726203044637328)
  a0806  =   T(-0.015319437748624402)
  a0807  =   T(0.008273789163814023)
  a0901  =   T(0.6241109587160757)
  a0904  =   T(-3.3608926294469414)
  a0905  =   T(-0.868219346841726)
  a0906  =   T(27.59209969944671)
  a0907  =   T(20.154067550477894)
  a0908  =   T(-43.48988418106996)
  a1001  =  T(0.47766253643826434)
  a1004  =  T(-2.4881146199716677)
  a1005  =  T(-0.590290826836843)
  a1006  =  T(21.230051448181193)
  a1007  =  T(15.279233632882423)
  a1008  =  T(-33.28821096898486)
  a1009  =  T(-0.020331201708508627)
  a1101  =  T(-0.9371424300859873)
  a1104  =  T(5.186372428844064)
  a1105  =  T(1.0914373489967295)
  a1106  =  T(-8.149787010746927)
  a1107  =  T(-18.52006565999696)
  a1108  =  T(22.739487099350505)
  a1109  =  T(2.4936055526796523)
  a1110  = T(-3.0467644718982196)
  a1201  = T(2.273310147516538)
  a1204  = T(-10.53449546673725)
  a1205  = T(-2.0008720582248625)
  a1206  = T(-17.9589318631188)
  a1207  = T(27.94888452941996)
  a1208  = T(-2.8589982771350235)
  a1209  = T(-8.87285693353063)
  a1210  = T(12.360567175794303)
  a1211  = T(0.6433927460157636)
  c14,c15,c16,a1401,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1613,a1614,a1615 = DP8Interp(T,T2)
  d401,d406,d407,d408,d409,d410,d411,d412,d413,d414,d415,d416,d501,d506,d507,d508,d509,d510,d511,d512,d513,d514,d515,d516,d601,d606,d607,d608,d609,d610,d611,d612,d613,d614,d615,d616,d701,d706,d707,d708,d709,d710,d711,d712,d713,d714,d715,d716 = DP8Interp_polyweights(T)
  DP8ConstantCache(c7,c8,c9,c10,c11,c6,c5,c4,c3,c2,b1,b6,b7,b8,b9,b10,b11,b12,btilde1,btilde6,btilde7,btilde8,btilde9,btilde10,btilde11,btilde12,er1,er6,er7,er8,er9,er10,er11,er12,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,c14,c15,c16,a1401,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1613,a1614,a1615,d401,d406,d407,d408,d409,d410,d411,d412,d413,d414,d415,d416,d501,d506,d507,d508,d509,d510,d511,d512,d513,d514,d515,d516,d601,d606,d607,d608,d609,d610,d611,d612,d613,d614,d615,d616,d701,d706,d707,d708,d709,d710,d711,d712,d713,d714,d715,d716)
end

function DP8ConstantCache(T::Type,T2::Type)
  c7     = T2(1//4)
  c8     = T2(4//13)
  c9     = T2(127//195)
  c10    = T2(3//5)
  c11    = T2(6//7)
  c6     = T2(4//3 * c7)
  c5     = T2((6+sqrt(6))/10 * c6)
  c4     = T2((6-sqrt(6))/10 * c6)
  c3     = T2(2//3 * c4)
  c2     = T2(2//3 * c3)
  b1     =  T(parse(BigFloat," 5.42937341165687622380535766363e-2"))
  b6     =  T(parse(BigFloat," 4.45031289275240888144113950566"))
  b7     =  T(parse(BigFloat," 1.89151789931450038304281599044"))
  b8     =  T(parse(BigFloat,"-5.8012039600105847814672114227"))
  b9     =  T(parse(BigFloat," 3.1116436695781989440891606237e-1"))
  b10    = T(parse(BigFloat,"-1.52160949662516078556178806805e-1"))
  b11    = T(parse(BigFloat," 2.01365400804030348374776537501e-1"))
  b12    = T(parse(BigFloat," 4.47106157277725905176885569043e-2"))
  # bhh1   = T(parse(BigFloat,"0.244094488188976377952755905512"))
  # bhh2   = T(parse(BigFloat,"0.733846688281611857341361741547"))
  # bhh3   = T(parse(BigFloat,"0.220588235294117647058823529412e-01"))
  btilde1  = T(parse(BigFloat,"-1.898007540724076157147023288757e-1"))
  btilde6  = T(parse(BigFloat," 4.45031289275240888144113950566"))
  btilde7  = T(parse(BigFloat," 1.89151789931450038304281599044"))
  btilde8  = T(parse(BigFloat,"-5.8012039600105847814672114227"))
  btilde9  = T(parse(BigFloat,"-4.22682321323791962932445679177e-1"))
  btilde10 = T(parse(BigFloat,"-1.52160949662516078556178806805e-1"))
  btilde11 = T(parse(BigFloat," 2.01365400804030348374776537501e-1"))
  btilde12 = T(parse(BigFloat,"2.26517921983608258118062039631e-2"))
  er1    = T(parse(BigFloat," 0.1312004499419488073250102996e-01"))
  er6    = T(parse(BigFloat,"-0.1225156446376204440720569753e+01"))
  er7    = T(parse(BigFloat,"-0.4957589496572501915214079952"))
  er8    = T(parse(BigFloat," 0.1664377182454986536961530415e+01"))
  er9    = T(parse(BigFloat,"-0.3503288487499736816886487290"))
  er10   = T(parse(BigFloat," 0.3341791187130174790297318841"))
  er11   = T(parse(BigFloat," 0.8192320648511571246570742613e-01"))
  er12   = T(parse(BigFloat,"-0.2235530786388629525884427845e-01"))
  a0201  =   T(parse(BigFloat," 5.26001519587677318785587544488e-2"))
  a0301  =   T(parse(BigFloat," 1.97250569845378994544595329183e-2"))
  a0302  =   T(parse(BigFloat," 5.91751709536136983633785987549e-2"))
  a0401  =   T(parse(BigFloat," 2.95875854768068491816892993775e-2"))
  a0403  =   T(parse(BigFloat," 8.87627564304205475450678981324e-2"))
  a0501  =   T(parse(BigFloat," 2.41365134159266685502369798665e-1"))
  a0503  =   T(parse(BigFloat,"-8.84549479328286085344864962717e-1"))
  a0504  =   T(parse(BigFloat," 9.24834003261792003115737966543e-1"))
  a0601  =   T(parse(BigFloat," 3.7037037037037037037037037037e-2"))
  a0604  =   T(parse(BigFloat," 1.70828608729473871279604482173e-1"))
  a0605  =   T(parse(BigFloat," 1.25467687566822425016691814123e-1"))
  a0701  =   T(parse(BigFloat," 3.7109375e-2"))
  a0704  =   T(parse(BigFloat," 1.70252211019544039314978060272e-1"))
  a0705  =   T(parse(BigFloat," 6.02165389804559606850219397283e-2"))
  a0706  =   T(parse(BigFloat,"-1.7578125e-2"))
  a0801  =   T(parse(BigFloat," 3.70920001185047927108779319836e-2"))
  a0804  =   T(parse(BigFloat," 1.70383925712239993810214054705e-1"))
  a0805  =   T(parse(BigFloat," 1.07262030446373284651809199168e-1"))
  a0806  =   T(parse(BigFloat,"-1.53194377486244017527936158236e-2"))
  a0807  =   T(parse(BigFloat," 8.27378916381402288758473766002e-3"))
  a0901  =   T(parse(BigFloat," 6.24110958716075717114429577812e-1"))
  a0904  =   T(parse(BigFloat,"-3.36089262944694129406857109825"))
  a0905  =   T(parse(BigFloat,"-8.68219346841726006818189891453e-1"))
  a0906  =   T(parse(BigFloat," 2.75920996994467083049415600797e1"))
  a0907  =   T(parse(BigFloat," 2.01540675504778934086186788979e1"))
  a0908  =   T(parse(BigFloat,"-4.34898841810699588477366255144e1"))
  a1001  =  T(parse(BigFloat," 4.77662536438264365890433908527e-1"))
  a1004  =  T(parse(BigFloat,"-2.48811461997166764192642586468e0"))
  a1005  =  T(parse(BigFloat,"-5.90290826836842996371446475743e-1"))
  a1006  =  T(parse(BigFloat," 2.12300514481811942347288949897e1"))
  a1007  =  T(parse(BigFloat," 1.52792336328824235832596922938e1"))
  a1008  =  T(parse(BigFloat,"-3.32882109689848629194453265587e1"))
  a1009  =  T(parse(BigFloat,"-2.03312017085086261358222928593e-2"))
  a1101  =  T(parse(BigFloat,"-9.3714243008598732571704021658e-1"))
  a1104  =  T(parse(BigFloat," 5.18637242884406370830023853209e0"))
  a1105  =  T(parse(BigFloat," 1.09143734899672957818500254654e0"))
  a1106  =  T(parse(BigFloat,"-8.14978701074692612513997267357e0"))
  a1107  =  T(parse(BigFloat,"-1.85200656599969598641566180701e1"))
  a1108  =  T(parse(BigFloat," 2.27394870993505042818970056734e1"))
  a1109  =  T(parse(BigFloat," 2.49360555267965238987089396762e0"))
  a1110  = T(parse(BigFloat,"-3.0467644718982195003823669022e0"))
  a1201  = T(parse(BigFloat,"  2.27331014751653820792359768449e0"))
  a1204  = T(parse(BigFloat," -1.05344954667372501984066689879e1"))
  a1205  = T(parse(BigFloat," -2.00087205822486249909675718444e0"))
  a1206  = T(parse(BigFloat," -1.79589318631187989172765950534e1"))
  a1207  = T(parse(BigFloat,"  2.79488845294199600508499808837e1"))
  a1208  = T(parse(BigFloat," -2.85899827713502369474065508674e0"))
  a1209  = T(parse(BigFloat," -8.87285693353062954433549289258e0"))
  a1210  = T(parse(BigFloat," 1.23605671757943030647266201528e1"))
  a1211  = T(parse(BigFloat," 6.43392746015763530355970484046e-1"))
  c14,c15,c16,a1401,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1613,a1614,a1615 = DP8Interp(T,T2)
  d401,d406,d407,d408,d409,d410,d411,d412,d413,d414,d415,d416,d501,d506,d507,d508,d509,d510,d511,d512,d513,d514,d515,d516,d601,d606,d607,d608,d609,d610,d611,d612,d613,d614,d615,d616,d701,d706,d707,d708,d709,d710,d711,d712,d713,d714,d715,d716 = DP8Interp_polyweights(T)
  DP8ConstantCache(c7,c8,c9,c10,c11,c6,c5,c4,c3,c2,b1,b6,b7,b8,b9,b10,b11,b12,btilde1,btilde6,btilde7,btilde8,btilde9,btilde10,btilde11,btilde12,er1,er6,er7,er8,er9,er10,er11,er12,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,c14,c15,c16,a1401,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1613,a1614,a1615,d401,d406,d407,d408,d409,d410,d411,d412,d413,d414,d415,d416,d501,d506,d507,d508,d509,d510,d511,d512,d513,d514,d515,d516,d601,d606,d607,d608,d609,d610,d611,d612,d613,d614,d615,d616,d701,d706,d707,d708,d709,d710,d711,d712,d713,d714,d715,d716)
end

function DP8Interp(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  c14    = T2(0.1)
  c15    = T2(0.2)
  c16    = T2(0.7777777777777778)
  a1401  = T(0.056167502283047954)
  a1407  = T(0.25350021021662483)
  a1408  = T(-0.2462390374708025)
  a1409  = T(-0.12419142326381637)
  a1410  = T(0.15329179827876568)
  a1411  = T(0.00820105229563469)
  a1412  = T(0.007567897660545699)
  a1413  = T(-0.008298)
  a1501  = T(0.03183464816350214)
  a1506  = T(0.028300909672366776)
  a1507  = T(0.053541988307438566)
  a1508  = T(-0.05492374857139099)
  a1511  = T(-0.00010834732869724932)
  a1512  = T(0.0003825710908356584)
  a1513  = T(-0.00034046500868740456)
  a1514  = T(0.1413124436746325)
  a1601  = T(-0.42889630158379194)
  a1606  = T(-4.697621415361164)
  a1607  = T(7.683421196062599)
  a1608  = T(4.06898981839711)
  a1609  = T(0.3567271874552811)
  a1613  = T(-0.0013990241651590145)
  a1614  = T(2.9475147891527724)
  a1615  = T(-9.15095847217987)

  return c14,c15,c16,a1401,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1613,a1614,a1615
end

function DP8Interp(T::Type,T2::Type)
  c14    = T2(1//10)
  c15    = T2(2//10)
  c16    = T2(7//9)
  a1401  = T(parse(BigFloat," 5.61675022830479523392909219681e-2"))
  a1407  = T(parse(BigFloat," 2.53500210216624811088794765333e-1"))
  a1408  = T(parse(BigFloat,"-2.46239037470802489917441475441e-1"))
  a1409  = T(parse(BigFloat,"-1.24191423263816360469010140626e-1"))
  a1410  = T(parse(BigFloat," 1.5329179827876569731206322685e-1"))
  a1411  = T(parse(BigFloat," 8.20105229563468988491666602057e-3"))
  a1412  = T(parse(BigFloat," 7.56789766054569976138603589584e-3"))
  a1413  = T(parse(BigFloat,"-8.298e-3"))
  a1501  = T(parse(BigFloat," 3.18346481635021405060768473261e-2"))
  a1506  = T(parse(BigFloat," 2.83009096723667755288322961402e-2"))
  a1507  = T(parse(BigFloat," 5.35419883074385676223797384372e-2"))
  a1508  = T(parse(BigFloat,"-5.49237485713909884646569340306e-2"))
  a1511  = T(parse(BigFloat,"-1.08347328697249322858509316994e-4"))
  a1512  = T(parse(BigFloat," 3.82571090835658412954920192323e-4"))
  a1513  = T(parse(BigFloat,"-3.40465008687404560802977114492e-4"))
  a1514  = T(parse(BigFloat," 1.41312443674632500278074618366e-1"))
  a1601  = T(parse(BigFloat,"-4.28896301583791923408573538692e-1"))
  a1606  = T(parse(BigFloat,"-4.69762141536116384314449447206e0"))
  a1607  = T(parse(BigFloat," 7.68342119606259904184240953878e0"))
  a1608  = T(parse(BigFloat," 4.06898981839711007970213554331e0"))
  a1609  = T(parse(BigFloat," 3.56727187455281109270669543021e-1"))
  a1613  = T(parse(BigFloat,"-1.39902416515901462129418009734e-3"))
  a1614  = T(parse(BigFloat," 2.9475147891527723389556272149e0"))
  a1615  = T(parse(BigFloat,"-9.15095847217987001081870187138e0"))

  return c14,c15,c16,a1401,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1613,a1614,a1615
end

function DP8Interp_polyweights(::Type{T}) where T<:CompiledFloats
  d401 = T(-8.428938276109013)
  d406 = T(0.5667149535193777)
  d407 = T(-3.0689499459498917)
  d408 = T(2.38466765651207)
  d409 = T(2.117034582445028)
  d410 = T(-0.871391583777973)
  d411 = T(2.2404374302607883)
  d412 = T(0.6315787787694688)
  d413 = T(-0.08899033645133331)
  d414 = T(18.148505520854727)
  d415 = T(-9.194632392478356)
  d416 = T(-4.436036387594894)
  d501 = T(10.427508642579134)
  d506 = T(242.28349177525817)
  d507 = T(165.20045171727028)
  d508 = T(-374.5467547226902)
  d509 = T(-22.113666853125306)
  d510 = T(7.733432668472264)
  d511 = T(-30.674084731089398)
  d512 = T(-9.332130526430229)
  d513 = T(15.697238121770845)
  d514 = T(-31.139403219565178)
  d515 = T(-9.35292435884448)
  d516 = T(35.81684148639408)
  d601 = T(19.985053242002433)
  d606 = T(-387.0373087493518)
  d607 = T(-189.17813819516758)
  d608 = T(527.8081592054236)
  d609 = T(-11.57390253995963)
  d610 = T(6.8812326946963)
  d611 = T(-1.0006050966910838)
  d612 = T(0.7777137798053443)
  d613 = T(-2.778205752353508)
  d614 = T(-60.19669523126412)
  d615 = T(84.32040550667716)
  d616 = T(11.99229113618279)
  d701 = T(-25.69393346270375)
  d706 = T(-154.18974869023643)
  d707 = T(-231.5293791760455)
  d708 = T(357.6391179106141)
  d709 = T(93.40532418362432)
  d710 = T(-37.45832313645163)
  d711 = T(104.0996495089623)
  d712 = T(29.8402934266605)
  d713 = T(-43.53345659001114)
  d714 = T(96.32455395918828)
  d715 = T(-39.17726167561544)
  d716 = T(-149.72683625798564)

  return d401,d406,d407,d408,d409,d410,d411,d412,d413,d414,d415,d416,d501,d506,d507,d508,d509,d510,d511,d512,d513,d514,d515,d516,d601,d606,d607,d608,d609,d610,d611,d612,d613,d614,d615,d616,d701,d706,d707,d708,d709,d710,d711,d712,d713,d714,d715,d716
end

function DP8Interp_polyweights(T::Type)
  d401 = T(parse(BigFloat,"-0.84289382761090128651353491142e+01"))
  d406 = T(parse(BigFloat," 0.56671495351937776962531783590e+00"))
  d407 = T(parse(BigFloat,"-0.30689499459498916912797304727e+01"))
  d408 = T(parse(BigFloat," 0.23846676565120698287728149680e+01"))
  d409 = T(parse(BigFloat," 0.21170345824450282767155149946e+01"))
  d410 = T(parse(BigFloat,"-0.87139158377797299206789907490e+00"))
  d411 = T(parse(BigFloat," 0.22404374302607882758541771650e+01"))
  d412 = T(parse(BigFloat," 0.63157877876946881815570249290e+00"))
  d413 = T(parse(BigFloat,"-0.88990336451333310820698117400e-01"))
  d414 = T(parse(BigFloat," 0.18148505520854727256656404962e+02"))
  d415 = T(parse(BigFloat,"-0.91946323924783554000451984436e+01"))
  d416 = T(parse(BigFloat,"-0.44360363875948939664310572000e+01"))
  d501 = T(parse(BigFloat," 0.10427508642579134603413151009e+02"))
  d506 = T(parse(BigFloat," 0.24228349177525818288430175319e+03"))
  d507 = T(parse(BigFloat," 0.16520045171727028198505394887e+03"))
  d508 = T(parse(BigFloat,"-0.37454675472269020279518312152e+03"))
  d509 = T(parse(BigFloat,"-0.22113666853125306036270938578e+02"))
  d510 = T(parse(BigFloat," 0.77334326684722638389603898808e+01"))
  d511 = T(parse(BigFloat,"-0.30674084731089398182061213626e+02"))
  d512 = T(parse(BigFloat,"-0.93321305264302278729567221706e+01"))
  d513 = T(parse(BigFloat," 0.15697238121770843886131091075e+02"))
  d514 = T(parse(BigFloat,"-0.31139403219565177677282850411e+02"))
  d515 = T(parse(BigFloat,"-0.93529243588444783865713862664e+01"))
  d516 = T(parse(BigFloat," 0.35816841486394083752465898540e+02"))
  d601 = T(parse(BigFloat," 0.19985053242002433820987653617e+02"))
  d606 = T(parse(BigFloat,"-0.38703730874935176555105901742e+03"))
  d607 = T(parse(BigFloat,"-0.18917813819516756882830838328e+03"))
  d608 = T(parse(BigFloat," 0.52780815920542364900561016686e+03"))
  d609 = T(parse(BigFloat,"-0.11573902539959630126141871134e+02"))
  d610 = T(parse(BigFloat," 0.68812326946963000169666922661e+01"))
  d611 = T(parse(BigFloat,"-0.10006050966910838403183860980e+01"))
  d612 = T(parse(BigFloat," 0.77771377980534432092869265740e+00"))
  d613 = T(parse(BigFloat,"-0.27782057523535084065932004339e+01"))
  d614 = T(parse(BigFloat,"-0.60196695231264120758267380846e+02"))
  d615 = T(parse(BigFloat," 0.84320405506677161018159903784e+02"))
  d616 = T(parse(BigFloat," 0.11992291136182789328035130030e+02"))
  d701 = T(parse(BigFloat,"-0.25693933462703749003312586129e+02"))
  d706 = T(parse(BigFloat,"-0.15418974869023643374053993627e+03"))
  d707 = T(parse(BigFloat,"-0.23152937917604549567536039109e+03"))
  d708 = T(parse(BigFloat," 0.35763911791061412378285349910e+03"))
  d709 = T(parse(BigFloat," 0.93405324183624310003907691704e+02"))
  d710 = T(parse(BigFloat,"-0.37458323136451633156875139351e+02"))
  d711 = T(parse(BigFloat," 0.10409964950896230045147246184e+03"))
  d712 = T(parse(BigFloat," 0.29840293426660503123344363579e+02"))
  d713 = T(parse(BigFloat,"-0.43533456590011143754432175058e+02"))
  d714 = T(parse(BigFloat," 0.96324553959188282948394950600e+02"))
  d715 = T(parse(BigFloat,"-0.39177261675615439165231486172e+02"))
  d716 = T(parse(BigFloat,"-0.14972683625798562581422125276e+03"))

  return d401,d406,d407,d408,d409,d410,d411,d412,d413,d414,d415,d416,d501,d506,d507,d508,d509,d510,d511,d512,d513,d514,d515,d516,d601,d606,d607,d608,d609,d610,d611,d612,d613,d614,d615,d616,d701,d706,d707,d708,d709,d710,d711,d712,d713,d714,d715,d716
end


text:
using JLD
using DifferentialEquations
using Plots
using Printf
include("../../../../src/estimate/ct_filters/ct_kalman_simple.jl")
#include("../../../../src/estimate/kalman.jl")
#include("../../../../src/estimate/ct_filters/ct_kalman_filter_mdn.jl")

EX =2
# set up transition/meas eqs matrixes
if EX == 1
    #  AR(1) model
    n_vars = 1
    n_states = 1
    n_shocks = 1
    T = Array{Float64}(undef,n_states,n_states)
    T .= -10.0
    R = Array{Float64}(undef,n_states,n_shocks)
    R = 10.0
    Z = Array{Float64}(undef,n_vars,n_states)
    Z .= 1.0

elseif EX == 2
    #  independent AR(1) models -- load only on the first one
    n_vars = 1
    n_states = 2
    n_shocks = 1
    T = Array{Float64}(undef,n_states,n_states)
    T[1,1] = -10.0
    T[2,2] = -10.0
    R = Array{Float64}(undef,n_states,n_shocks)
    R .= 1.0
    @show R
    Z = Array{Float64}(undef,n_vars,n_states)
    Z[1,1] = 1.0
    Z[1,2] = 0.0
end

# generate data (solve using Euler-maruyama)

n_quarters = 200.
global s_t = zeros(Float64, n_states) #initial values

# #use Differential Equations
# # Set up SDE
# f(u,p,t) = T*u
# g(u,p,t) = R
# dt = 1.0/90.0 # daily frequency when one period is one quarter
# tspan = (0.0,n_quarters) # 50 years of data
# W = WienerProcess(0.0, 0., 0.)
# prob = SDEProblem(f, g, s_t, tspan, W)
#      sol = solve(prob, EM(), dt = dt)
# s_1T = vcat(sol.u...)
# data_d = s_1T*Z'

#    gr()
#     default(show = true)
# p = plot(sol.t, data_d)
# sleep(2)

data_d = Array{Float64}(undef,Int(round(n_quarters/dt)),n_vars)
time_d = Array{Float64}(undef,Int(round(n_quarters/dt)))
n_step = 10.0
for i_t = 1:Int(round(n_quarters/dt))
    for i = 1:Int(n_step)
        s_t = s_t + T*s_t*dt/n_step + sqrt(dt/n_step)*R*randn(n_shocks, 1);
    end
    data_d[i_t,:] = Z*s_t
    time_d[i_t] = i_t*dt
end

gr()
default(show = true)
p = plot(time_d, data_d, color = :red)

# created 50 years worth of daily data -> transform data sampled at quarterly frequency for discrete observations
Delta_t = 1.0#1.0/90.0##1.0
del_t = dt./Delta_t

@show n_y = Int(round.(length(data_d)*dt/Delta_t;digits = 0))
data = Array{Float64}(undef,n_y,n_vars)
for i = 1:n_y
    data[i,:] .= data_d[1+(i-1)*Int(floor(Delta_t/dt)),:]
end

E = zeros(n_vars, n_vars) # meas error
Q = Matrix{Float64}(I,n_shocks,n_shocks) # variance of shocks
mean_0 = zeros(n_states,1) # intial state mean
var_0 = 0.1*Matrix{Float64}(I,n_states,n_states) # initial state variance
for sig = 0.1:0.1:2.0
    prob_out = ct_kalman_simple(T, Z, R*(sig*Q)*R', E, mean_0, var_0, data, [Delta_t])
    @show [sig prob_out]
end

# not sure the Delta_t abd dt are correct
#E = zeros(n_vars, n_vars) * Delta_t # Shock is Browian motion -> variance is Δ_t b/c this is quarterly observation data
#Q = Matrix{Float64}(I,n_shocks,n_shocks) * dt    # Shock is Brownian motion -> variance is dt b/c this is state data
#out = ct_kalman_filter(data, T, R, C, Q, Z, D, E, Delta_t)
#@show true_lik = sum(out[1])

 #   #, label = "High vol",xlims=(glimpdf[1],glimpdf[2]))

#=
SeHyoun's function

% Test Kalman Filer <kalman_ct>
n_simul = 1000;
n_step = 50;
make_plots = false;
F = [-1, 0; 0, -1];
F_sim = F;

Q = 0.02*speye(2);
% Q = [0.1, 0.02; 0.02, 0.1];
R = 0.01*speye(2);
% R = [0.1, 0.02; 0.02, 0.1];
H = speye(2);

mean_init = [2; -2];
var_init = 0.01*speye(2);

mean_now = mean_init;
data_x = mean_now + sqrt(0.1)*randn(2,1);
var_now = Var_init;

dt = 0.25;

F_val = 0.1:0.1:1.9;
data_y_iter = zeros(2, n_simul);
for iter_time = 1:n_simul
  for i = 1:n_step
    data_x = data_x + F*data_x*dt/n_step + sqrt(dt/n_step)*sqrtm(Q)*randn(2, 1);
  end
  data_y = H*data_x + sqrtm(R)*randn(2, 1);
  data_y_iter(:, iter_time) = data_y;
end
figure(101);
% plot(0:dt:dt*(n_simul-1), data_y_iter);
plot(data_y_iter(1,:), data_y_iter(2, :), '.-');

time_step = dt*ones(n_simul, 1);
for iter_outer = 1:length(F_val)

  F_sim = F_val(iter_outer)*F;

  [mean_now, var_now, log_like] = kalman_ct(F_sim, H, Q, R, mean_init, var_init, data_y_iter, time_step);

  log_like_iter(iter_outer) = log_like;

  figure(4);
  clf;
  plot(F_val(1:iter_outer), log_like_iter(1:iter_outer));
  xlim([0.1, 1.9]);
  drawnow
end

=#

ext: jl


text:
##############################################
### Interactions with other point types    ###
##############################################

# convert to a geodesy type - overload this with user methods
"""
Convert to geodesy representations of a custom type

If the input is a Type, a Geodesy type should be returned
If the input is an object, a Geodesy object should be returned
"""
geodify{T}(::Type{T}) = error("Function to retrieve the appropriate geodesy type  for a type $(T) is not defined")
geodify(X) = error("No conversion is defined to convert from type $(typeof(X)) to a geodesy type")

# and set a default transformation to the desired output type
geotransform_ext{T}(::Type{T}, X) = T(geotransform(geodify(T), geodify(X)), X)

# The below was for when I was allowing extra arguments to pass to the output type constructor
# geotransform_ext{T}(::Type{T}, X...) = T(geotransform(geodify(T), geodify(X[1])), X...)

# define geodify for all known types
for geo_type in Geodesy.GeodesyTypes
    q = quote
        geodify(::Type{$(geo_type)})       = $(geo_type)
        geodify{T}(::Type{$(geo_type){T}}) = $(geo_type){T}
        geodify(X::$(geo_type))            = X
    end
    eval(q)     # TODO: eval? figure out the right way to do this
end


text:
module GNSSBenchmarks
    using
        CUDA,
        CSV,
        StructArrays,
        BenchmarkTools,
        LoopVectorization,
        Tracking,
        TrackingLoopFilters,
        GNSSSignals,
        DataFrames,
        PGFPlotsX,
        LinearAlgebra,
        StaticArrays

    using Unitful: upreferred, Hz, dBHz, ms, kHz, MHz

    import
        LinearAlgebra.dot,
        Base.length,
        Tracking.TrackingState,
        Tracking.NumAnts,
        Tracking.MomentsCN0Estimator,
        Tracking.AbstractCN0Estimator,
        Tracking.AbstractCorrelator,
        Tracking.EarlyPromptLateCorrelator,
        Tracking.SecondaryCodeOrBitDetector,
        Tracking.GainControlledSignal,
        Tracking.found,
        Tracking.TrackingResults,
        Tracking.BitBuffer,
        Tracking.get_default_post_corr_filter,
        Tracking.get_num_samples_left_to_integrate,
        Tracking.get_num_samples,
        Tracking.get_num_ants,
        Tracking.get_current_carrier_frequency,
        Tracking.get_current_code_frequency,
        TrackingLoopFilters.AbstractLoopFilter,
        GNSSSignals.AbstractGNSSSystem

    export
        main,
        plotall,
        benchmark_downconvert,
        benchmark_carrier_replica,
        benchmark_code_replica,
        benchmark_correlate,
        benchmark_tracking_loop,
        plot_carrier_replica,
        gpu_correlate,
        plot_correlate_all,
        sgpu_correlate,
        plot_code_replica_all

    include("main.jl")
    include("plot.jl")
    include("gpu_downconvert.jl")
    include("gpu_carrier_replica.jl")
    include("gpu_code_replica.jl")
    include("gpu_correlate.jl")
    include("gpu_tracking_state.jl")
    include("gpu_tracking_results.jl")
    include("gpu_tracking_loop.jl")
    include("benchmark_carrier_replica.jl")
    include("benchmark_code_replica.jl")
    include("benchmark_downconvert.jl")
    include("benchmark_correlate.jl")
    include("benchmark_gpu_tracking_loop.jl")

    const MAX_NUM_SAMPLES = 50000
    const SAMPLES = StepRange(2500,2500,MAX_NUM_SAMPLES)
    const ANTENNA = [1, 4, 16]
end


text:
abstract type AbstractARGNode end

"""
	mutable struct ARGNode

Node of an Ancestral Reassortment Graph representing two genes / segments.

## Fields:
```
	anc :: Vector{Union{Nothing, AbstractARGNode}} = [nothing, nothing] # Always of length 2
	color :: Vector{Bool} = [false, false]
	tau :: Vector{Union{Missing, Float64}} = [missing, missing]
	children :: Vector{ARGNode} = ARGNode[]
	label :: String = make_random_label()
	hybrid :: Bool = false # Hybrids can have the same ancestor for some reassortments
```
"""
@with_kw_noshow mutable struct ARGNode <: AbstractARGNode
	anc :: Vector{Union{Nothing, AbstractARGNode}} = [nothing, nothing] # Always of length 2
	color :: Vector{Bool} = [false, false]
	tau :: Vector{Union{Missing, Float64}} = [missing, missing]
	children :: Vector{ARGNode} = ARGNode[]
	label :: String = make_random_label()
	hybrid :: Bool = false # Hybrids can have the same ancestor for some reassortments
end

struct RootNode <: AbstractARGNode
	# children :: Vector{ARGNode}
end

# Equality
isequal(a::ARGNode, b::ARGNode) = (a.label==b.label)
isequal(a::ARGNode, b::RootNode) = false
isequal(a::RootNode, b::ARGNode) = false
isequal(a::ARGNode, b::Nothing) = false
isequal(a::Nothing, b::ARGNode) = false
Base.:(==)(x::ARGNode, y::ARGNode) = isequal(x,y)
Base.hash(x::ARGNode, h::UInt) = hash(x.label, h)

# Access
ancestors(n::ARGNode) = n.anc
function ancestor(n::ARGNode, c::Int)
	@assert hascolor(n, c)
	return n.anc[c]
end

children(n::ARGNode) = n.children
children(n::RootNode) = n.children
"""
	children(n::ARGNode, clr)

Children of `n` for color `clr`. Child `c` must meet the conditions:
- `hascolor(c, clr)`
- `ancestor(c, clr) == n`
"""
function children(n::ARGNode, clr)
	Iterators.filter(children(n)) do c
		hascolor(c, clr) && ancestor(c, clr) == n
	end
end

color(n::ARGNode) = n.color
color(::RootNode) = nothing

label(n::ARGNode) = n.label
label(::RootNode) = "_RootNode_"
label(::Nothing) = nothing

branch_length(n::ARGNode) = n.tau
branch_length(n::RootNode) = 0.
tau(n) = branch_length(n)
function branch_length(n::ARGNode, clr)
	@assert hascolor(n, clr)
	return n.tau[clr]
end
function branch_length(n::RootNode, clr)
	@assert hascolor(n, clr)
	return 0.
end

## isroot / isleaf
"""
	isroot(n)

Is `n` the root for one color?
"""
isroot(n::ARGNode) = in(RootNode(), n.anc)
"""
	isroot(n, c)

Is `n` the root for color `c`?
"""
isroot(n::ARGNode, c::Int) = (n.anc[c] == RootNode())
function is_global_root(n::ARGNode)
	if isroot(n)
		for c in 1:length(n.anc)
			if !isroot(n, c)
				return false
			end
		end
		return true
	end
	return false
end
function is_partial_root(n::ARGNode)
	if !isroot(n)
		return false
	else
		return !is_global_root(n)
	end
end

isleaf(n::ARGNode) = isempty(children(n))

ishybrid(n::ARGNode) = n.hybrid
ishybrid(::RootNode) = false

degree(n) = sum(color(n))
isshared(n) = (degree(n) == length(color(n)))

# Tests
hasancestor(n::ARGNode, c::Int) = !isnothing(ancestor(n,c)) && !isroot(n,c)
function hasancestor(n::ARGNode, a::ARGNode)
	for x in ancestors(n)
		if x == a
			return true
		end
	end
	return false
end

function haschild(a::ARGNode, n::ARGNode)
	for c in children(a)
		if c == n
			return true
		end
	end
	return false
end

hascolor(n::ARGNode, c::Int) = n.color[c]
hascolor(::RootNode, c::Int) = true
hascolor(::Nothing, c::Int) = true
function share_color(n1::ARGNode, n2::ARGNode)
	for (c1,c2) in zip(color(n1), color(n2))
		if c1 && c2
			return true
		end
	end
	return false
end
share_color(n::ARGNode, r::RootNode) = true
share_color(r::RootNode, n::ARGNode) = true

# Set / add children / ancestors
function add_child!(n::ARGNode, c::ARGNode)
	@assert !haschild(n, c) && share_color(n, c)
	push!(children(n), c)
end

function set_ancestor!(n, a, c::Int)
	@assert hascolor(n, c) && hascolor(a, c) "Nodes do not share color $c: \n $(label(n)): $(color(n)) \n $(label(a)): $(color(a))"
	n.anc[c] = a
end
function set_ancestor!(n, a)
	for (i,c) in color(n)
		if c
			set_ancestor!(n, a, i)
		end
	end
end

function add_color!(n, c::Int)
	@assert !hascolor(n, c) "Node $(n.label) already has color $c"
	n.color[c] = true
end

function set_branch_length!(n, τ::Real, c, clrs...)
	@assert hascolor(n, c)
	@assert hasancestor(n, c)
	n.tau[c] = τ + eps()
	for clr in clrs
		set_branch_length!(n, τ, clr)
	end
end
function set_branch_length!(n, ::Missing, c, clrs...)
	@assert hascolor(n, c)
	n.tau[c] = missing
	for clr in clrs
		set_branch_length!(n, missing, clr)
	end
end

set_label!(n::ARGNode, label) = (n.label = label)

# Utils.
function othercolor(c)
	@assert (c==1 || c==2) "Invalid color $c"
	return (c==1) ? 2 : 1
end
"""
	edgecolor(a, n)

Color of the edge going from `a` to `n`
"""
function edgecolor(a, n)
	clr = [false, false]
	for (i,c) in enumerate(color(n))
		if c && a == ancestor(n, i)
			clr[i] = true
		end
	end
	return clr
end

mutable struct ARG
	roots :: Dict{Int, ARGNode}
	nodes :: Dict{String, ARGNode}
	hybrids :: Dict{String, ARGNode}
	leaves :: Dict{String, ARGNode}
end
roots(arg::ARG) = arg.roots
nodes(arg::ARG) = values(arg.nodes)
hybrids(arg::ARG) = values(arg.hybrids)
leaves(arg::ARG) = values(arg.leaves)

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/RecombTools

max_forks_repo_name: PierreBarrat/RecombTools


text:
## Exercise 5-7
## https://benlauwens.github.io/ThinkJulia.jl/latest/images/fig52.svg
## Figure 8. A Koch curve

## The Koch curve is a fractal that looks something like A Koch curve. To draw a Koch curve with length x, all you have to do is

## 1. Draw a Koch curve with length x/3.
## 2. Turn left 60 degrees.
## 3. Draw a Koch curve with length x/3.
## 4. Turn right 120 degrees.
## 5. Draw a Koch curve with length x/3.
## 6. Turn left 60 degrees.
## 7. Draw a Koch curve with length x/3.

## The exception is if x is less than 3: in that case, you can just draw a straight line with length x.

using ThinkJulia

## 1. Write a function called koch that takes a turtle and a length as parameters, and that uses the turtle to draw a Koch curve with the given length.
println("Ans 1: ")

function koch(turtle::Turtle, length)
    if length < 3
        forward(turtle, length)
        return
    end

    koch(turtle, length/3)
    turn(turtle, 60)
    koch(turtle, length/3)
    turn(turtle, -120)
    koch(turtle, length/3)
    turn(turtle, 60)
    koch(turtle, length/3)
end

turtle = Turtle()
# @svg begin
#     koch(turtle, 100)
# end

## 2. Write a function called snowflake that draws three Koch curves to make the outline of a snowflake.
println("Ans 2: ")

function snowflake(turtle, length)
    koch(turtle, length)
    turn(turtle, -60)
    koch(turtle, length)
    turn(turtle, -60)
    koch(turtle, length)
    turn(turtle, -60)
    koch(turtle, length)
    turn(turtle, -60)
    koch(turtle, length)
    turn(turtle, -60)
    koch(turtle, length)
    turn(turtle, -60)
end

# @svg begin
#     snowflake(turtle, 100)
# end

## 3. The Koch curve can be generalized in several ways. See https://en.wikipedia.org/wiki/Koch_snowflake for examples and implement your favorite.
println("Ans 3: ")

function koch_square(turtle::Turtle, length)
    if length < 3
        forward(turtle, length)
        return
    end

    koch_square(turtle, length/3)
    turn(turtle, 90)
    koch_square(turtle, length/3)
    turn(turtle, -90)
    koch_square(turtle, length/3)
    turn(turtle, -90)
    koch_square(turtle, length/3)
    turn(turtle, 90)
    koch_square(turtle, length/3)
    turn(turtle, 90)
end

@svg begin
    koch_square(turtle, 60)
end

println("End.")


text:
import Base: convert, print, show, push!, length, getindex, sort!, sort, +,-,*,.*,==,/, setindex!,replace,start,next,done,zero

abstract Component
function ==(c1::Component, c2::Component)
	if isa(c1,typeof(c2))
		for n in fieldnames(c1)
			if getfield(c1,n)!=getfield(c2,n)
				return false
			end
		end
		return true
	end
	return false
end
abstract SingleArg <: Component
==(sa1::SingleArg,sa2::SingleArg)=isa(sa1,typeof(sa2))&&sa1.x==sa2.x
abstract NonAbelian <: Component
abstract Operator <: Component
function getargs(c::Component)
	n=fieldnames(c)
	ret=Any[]
	for nam in n
		push!(ret,getfield(c,nam))
	end
	return ret
end
function getarg(c::Component,n::Integer=1)
	getfield(c,fieldnames(c)[n])
end
function setarg!(c::Component,newarg,argn::Integer=1)
	setfield!(c,fieldnames(c)[argn],newarg)
	return c
end
setarg(c::Component,newarg,argn::Integer=1)=setarg!(deepcopy(c),newarg,argn)
type Components <: Component
	components
	coef
end #maybe this type should be deprecated and sumsym rewritten
function ==(cs1::Components,cs2::Components)
	nc=length(cs1.components)
	if nc != length(cs2.components)
		return false
	else
		for c in 1:nc
			if length(indsin(cs1.components,cs1.components[c]))!=length(indsin(cs2.components,cs1.components[c]))
				return false
			end
		end
	end
	return true
end

type Expression
	terms::Array{Array{Union{Number,Symbol,Component,Expression},1},1}
end
length(ex::Expression)=length(ex.terms)
getindex(ex::Expression,i::Integer)=getindex(ex.terms,i)
function getindex(ex::Expression,t::Array)
	p=ex[t[1]]
	if length(t)==1
		return p
	end
	for i in 2:length(t)-1
		p=p[t[i]]
	end
	return p[t[end]]
end
function setindex!(ex::Expression,a,t::Array)
	p=ex[t[1]]
	for i in 2:length(t)-1
		p=p[t[i]]
	end
	p[t[end]]=a
end
function getindex(c::Component,t::Array)
	p=getarg(c,t[1])
	for i in 2:length(t)-1
		p=p[t[i]]
	end
	return p[t[end]]
end
function setindex!(c::Component,a,t::Array)
	p=getarg(c,t[1])
	for i in 2:length(t)-1
		p=p[t[i]]
	end
	p[t[end]]=a
end
function print(io::IO,c::Union{Number,Component})
	if isa(c,Number)&&(isa(c,Complex)||c<0)
		print(io,'(')
		show(io,c)
		print(io,')')
	elseif isa(c,SingleArg)
		print(io,typeof(c),'(')
		print(io,getarg(c))
		print(io,')')
	elseif isa(c,Component)
		print(io,typeof(c),'(')
		ca=getargs(c)
		for a in 1:length(ca)-1
			print(io,ca[a])
			print(io,',')
		end
		print(io,ca[end])
		print(io,')')
	else
		show(io,c)
	end
end
function preprint(io::IO,fac)
	if isa(fac,Expression)
		print(io,'(')
		print(io,fac)
		print(io,')')
	else
		print(io,fac)
	end
end
function print(io::IO,ex::Expression)
	if isempty(ex.terms)
		print(io,ex.terms)
	else
		for term in 1:length(ex.terms)-1
			for fac in 1:length(ex.terms[term])-1
				preprint(io,ex.terms[term][fac])
				print(io,' ')
			end
			preprint(io,ex.terms[term][end])
			print(io," + ")
		end
		for fac in 1:length(ex.terms[end])-1
			preprint(io,ex.terms[end][fac])
			print(io,' ')
		end
		preprint(io,ex.terms[end][end])
	end
end
N=Union{Number,Symbol}
X=Union{Number,Symbol,Component}
Ex=Union{Symbol,Component,Expression}
EX=Union{Number,Symbol,Component,Expression}
typealias Factor EX
show(io::IO,x::Type{Factor})=print(io, "Factor")
typealias Term Array{Factor,1}
convert(::Type{Array{Array{Factor,1},1}},a::Array{Any,1})=Term[a]
zero(f::Factor)=0

macro delegate(source, targets) # by JMW
    typename = esc(source.args[1])
    fieldname = esc(Expr(:quote, source.args[2].args[1]))
    funcnames = targets.args
    n = length(funcnames)
    fdefs = Array(Any, n)
    for i in 1:n
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typename), args...) =
                       ($funcname)(a.($fieldname), args...)
                   end
    end
    return Expr(:block, fdefs...)
end

complexity(n::N)=1
function complexity(c::Component)
	tot=0
	for n in fieldnames(c)
		tot+=complexity(getfield(c,n))
	end
	return tot
end
function complexity(term::Term)
	tot=0
	for factor in term
		tot+=complexity(factor)
	end
	return tot
end
function complexity(ex::Expression)
	tot=0
	for term in ex
		tot+=complexity(term)
	end
	return tot
end
function complexity(a::Array)
	tot=0
	for item in a
		tot+=complexity(item)
	end
	return tot
end
expression(a::Term)=Expression(Term[a])
expression(a::Array{Term})=Expression(a)
function expression(cs::Array{Components})
	if isempty(cs)
		return 0
	end
	ex=Expression(Term[])
	for cc in cs
		if cc.coef!=0
			nterm=Factor[]
			if cc.coef!=1||isempty(cc.components)
				push!(nterm,cc.coef)
			end
			for c in cc.components
				push!(nterm,c)
			end
			push!(ex.terms,nterm)
		end
	end
	if length(ex)==0
		return 0
	else
		return ex
	end
end
expression(x::X)=expression(Factor[x])
==(ex1::Expression,ex2::Expression)=ex1.terms==ex2.terms
push!(ex::Expression,a)=push!(ex.terms,a)
push!(x::X,a)=expression(Factor[x,a])
+(ex1::Expression,ex2::Expression)=begin;ex=deepcopy(ex1);push!(ex.terms,[ex2]);ex;end
+(ex::Expression,a::X)=begin;ex=deepcopy(ex);push!(ex.terms,[a]);ex;end
-(ex::Expression,a::X)=begin;ex=deepcopy(ex);push!(ex.terms,[-1,a]);ex;end
-(a::X,ex::Expression)=a+(-1*ex)
-(ex1::Expression,ex2::Expression)=begin;ex=deepcopy(ex1);push!(ex,Factor[-1,ex2]);ex;end
+(a::X,ex::Expression)=begin;ex=deepcopy(ex);insert!(ex.terms,1,Factor[a]);ex;end
*(ex1::Expression,ex2::Expression)=expression(Factor[deepcopy(ex1),deepcopy(ex2)])
#.*(a::Array,ex::Ex)=ex.*a
function .*(ex::Ex,a::Array)
	na=EX[]
	for i in 1:length(a)
		push!(na,ex*a[i])
	end
	return na
end
.*(n::Number,ex::EX)=*(n,ex)
/(ex::Ex,n::Number)=*(1/n,ex)
function *(a::X,ex::Expression)
	ex=deepcopy(ex)
	for ti in 1:length(ex)
		insert!(ex[ti],1,a)
	end
	return ex
end

*(a::Number,c::Component)=Expression([a,c])
*(c1::Union{Component,Symbol},c2::Union{Component,Symbol})=Expression(Term[Factor[c1,c2]])
function *(ex::Expression,x::EX)
	ap=terms(deepcopy(ex))
	for t in ap
		push!(t,x)
	end
	return expression(ap)
end
-(c1::Component,c2::Component)=+(c1,-1*c2)
+(c1::Component,c2::Component)=Expression(Term[Factor[c1],Factor[c2]])
+(c::Component,ex::Expression)=Expression(Term[Factor[c],Factor[ex]])
+(ex::Expression,c::Component)=Expression(Term[Factor[ex],Factor[c]])

+(x1::X,x2::X)=Expression(Term[[x1],[x2]])
-(x1::X,x2::X)=Expression(Term[[x1],[-1,x2]])
-(x::X)=Expression([-1,x])
-(ex::Expression)=-1*ex
*(x1::X,x2::X)=Expression([x1,x2])

abstract Operation <: Component

function indin(array,item)
	ind=0
	for it in array
		ind+=1
		if it==item
			return ind
		end
	end
	return 0
end
function indin(array,typ::Type)
	for it in 1:length(array)
		if isa(array[it],typ)
			return it
		end
	end
	return 0
end
function expandindices(inds::Tuple,ninds::Array{Array{Integer}}=Array{Integer}[],nind::Array{Integer}=Integer[])
	for l in 1:length(inds)
		if isa(inds[l],Integer)
			push!(nind,inds[l])
		elseif isa(inds[l],Array)
			for a in inds[l]
				if isa(a,Tuple)
					tia=expandindices(a)
					for ti in deepcopy(tia)
						nnind=deepcopy(nind)
						pushall!(nnind,ti)
						push!(ninds,nnind)
					end
				else
					nnind=deepcopy(nind)
					push!(nnind,a)
					push!(ninds,nnind)
				end
			end
		else
			println(inds)
		end
	end
	return ninds
end
function expandindices(a::Array)
	na=Array{Integer}[]
	for i in 1:length(a)
		pushall!(na,expandindices(a[i]))
	end
	return na
end
function indsin(array::Array,item)
	ind=Int64[]
	for it in 1:length(array)
		if array[it]==item
			push!(ind,it)
		end
	end
	return ind
end
function indsin(ex::Expression,item)
	#if simplify(ex)<ex
	#	warn("Please call simplify on $ex")
	#end
	inds=Tuple[]
	for ti in 1:length(ex)
		ind=indsin(ex[ti],item)
		if !isempty(ind)
			push!(inds,(ti,ind))
		end
	end
	return inds
end
function indsin(c::Component,item)
	inds=Any[]
	args=getargs(c)
	for ai in 1:length(args)
		if args[ai]==item||isa(args[ai],item)
			push!(inds,ai)
			continue
		elseif isa(args[ai],N)
			continue
		end
		ind=indsin(args[ai],item)
		if !isempty(ind)
			push!(inds,(ai,ind))
		end
	end
	return inds
end
function indsin(array::Array,typ::Type)
	ind=Int64[]
	for it in 1:length(array)
		if isa(array[it],typ)
			push!(ind,it)
		end
	end
	return ind
end
indsin(n::Number,a)=[]
terms(ex::Expression)=ex.terms
terms(x::X)=x #Term[Factor[x]] #hmm
dcterms(ex::Expression)=terms(deepcopy(ex))
dcterms(x::X)=x
function has(a::Array,t::Type)
	for it in a
		if isa(it,t)
			return true
		end
	end
	return false
end
function has(a::Array,t::EX)
	for it in a
		if it==t
			return true
		end
	end
	return false
end
function has(term1::Term,term2::Term)
	if length(term1)<length(term2)
		return false
	end
	for p1 in permutations(term1)
		for p2 in permutations(term2)
			for l in 1:length(term2)
				if p1[1:l]==p2[1:l]
					return true
				end
			end
		end
	end
	return false
end
function has(ex::Expression,x::EX)
	for term in ex
		for c in term
			if c==x || (isa(c,Expression)&&has(c,x)) || (isa(c,Component)&&has(c,x))
				return true
			end
		end
	end
	return false
end
function has(ex::Expression,t::Type)
	for term in ex
		for c in term
			if isa(c,t) || (isa(c,Expression)&&has(c,t)) || (isa(c,Component)&&has(c,t))
				return true
			end
		end
	end
	return false
end
function has(c::Component,x::Symbol)
	for a in getargs(c)
		if has(a,x)
			return true
		end
	end
	return false
end
function has(c::Component,x::Type)
	if isa(c,x)
		return true
	end
	for a in getargs(c)
		if has(a,x)
			return true
		end
	end
	return false
end
has(n::N,x::Symbol)=n==x
has(::N, ::Type)=false
function maketype(c::Component,fun)
	args=getargs(c)
	for argi in 1:length(args)
		args[argi]=fun(args[argi])
	end
	return typeof(c)(args...)
end
function unnest(ex::Expression)
	nt=Term[]
	for term in ex.terms
		nf=Factor[]
		for fac in term
			if isa(term,Array)
				for nested in term
					push!(nf,nested)
				end
			else
				push!(nf,fac)
			end
		end
		push!(nt,nf)
	end
	return Expression(nt)
end
function componify(ex::Expression,raw=false)
	ap=dcterms(ex)
	lap=length(ap)
	for term in 1:lap
		exs=Expression[]
		xs=X[]
		for fac in ap[term]
			if isa(fac,Array)
				#warn("How did the array $fac end up in $ex?") #through replace!
				push!(exs,componify(Expression(fac)))
			elseif isa(fac,Expression)
				push!(exs,componify(fac))
			elseif isa(fac,Component)
				fac=maketype(fac,componify)
				push!(xs,fac)
			else
				push!(xs,fac)
			end
		end
		if isempty(exs)
			continue
		else
			tap=exs[1].terms
			for x in xs
				for tterm in tap
					unshift!(tterm,x)
				end
			end
			exs[1]=expression(tap)
			while length(exs)>1
				ap1=terms(exs[1])
				ap2=terms(exs[2])
				lap1,lap2=length(ap1),length(ap2)
				multed=Term[]
				for l in 1:lap1*lap2
					push!(multed,Factor[])
				end
				for l1 in 0:lap1-1
					for l2 in 1:lap2
						for fac in ap1[l1+1]
							push!(multed[l2+l1*lap2],fac)
						end
						for fac in ap2[l2]
							push!(multed[l2+l1*lap2],fac)
						end
					end
				end
				exs[2]=expression(multed)
				deleteat!(exs,1)
			end
			ap[term]=terms(exs[1])[1]
			for tterm in 2:length(exs[1])
				push!(ap,exs[1][tterm])
			end
		end
	end
	if raw
		return ap #this was needed for previous componify, maybe remove?
	else
		return expression(ap)
	end
end
function componify(a::Array)
	if length(a)==1
		return componify(extract(a))
	end
	a=deepcopy(a)
	for i in 1:length(a)
		a[i]=componify(a[i])
	end
	return a
end
componify(a::Array{Array})=componify(expression(a))
componify(c::Component)=maketype(c,componify)
componify(x::N)=x
function extract(ex::Expression)
	if length(ex.terms)==1&&length(ex.terms[1])==1
		return ex[1][1]
	end
	return ex
end
function extract(a::Array)
	if length(a)==1
		return a[1]
	end
	return a
end
import Base.isless
isless(ex::Expression,x::X)=false
isless(x::X,ex::Expression)=true
isless(c::Component,s::N)=false
isless(s::N,c::Component)=true
isless(s::Symbol,n::Number)=false
isless(n::Number,s::Symbol)=true
isless(s::Complex,n::Complex)=real(s)<real(n)?true:imag(s)<imag(n)?true:false
isless(n::Real,s::Complex)=true
isless(s::Complex,n::Real)=false
isless(na::NonAbelian,na2::NonAbelian)=false
isless(a::Symbol,op::Operator)=false
isless(op::Operator,a::Symbol)=false
function isless(c1::Component,c2::Component)
	if isa(c1,Operator)||isa(c2,Operator)
		return false
	end
	xi1=complexity(c1)
	xi2=complexity(c2)
	xi1==xi2?isless(string(c1),string(c2)):xi1<xi2
end
function isless(t1::Term,t2::Term)
	xi1=complexity(t1)
	xi2=complexity(t2)
	xi1==xi2?isless(string(t1),string(t2)):xi1<xi2
end
function isless(ex1::Expression,ex2::Expression)
	xi1=complexity(ex1)
	xi2=complexity(ex2)
	xi1==xi2?isless(string(ex1),string(ex2)):xi1<xi2
end
sort(n::N)=n
sort(c::Component)=maketype(c,sort)
function sort!(ex::Expression)
	ap=terms(ex)
	for term in ap
		sort!(term)
	end
	return expression(sort!(ap))
end
sort(ex::Expression)=sort!(deepcopy(ex))
include("tensors.jl")
function simplify(ex::Expression)
	tex=0
	nit=0
	while tex!=ex
		tex=ex
		ex=sumsym(sumnum(componify(ex)))
		ap=terms(ex)
		if isa(ap,X)
			return ap
		end
		for term in 1:length(ap)
			ap[term]=divify!(ap[term])
			ap[term]=divbine!(ap[term])
			ap[term]=divbinedify!(ap[term])
			for fac in 1:length(ap[term])
				ap[term][fac]=simplify(ap[term][fac])
			end
			unsqrt!(ap[term])
			#sort!(ap[term])
		end
		ex=extract(expression(ap)) #better to check if res::N before calling expression instead of extracting?
		nit+=1
		if has(ex,Ten)
			ex=simplify(ex,Ten)
		end
		if nit>90
			warn("Stuck in simplify! Iteration $nit: $ex")
			break
		end
	end
	return ex#sort!(ex)
end
function simplify(c::Component)
	args=deepcopy(getargs(c))
	for arg in 1:length(args)
		args[arg]=simplify(args[arg])
	end
	return typeof(c)(args...)
end
simplify(x::N)=x
simplify(x::N,a)=x
simplify!(x::N)=x
function simplify!(a::Array)
	if length(a)==1
		return simplify(a[1])
	else
		for i in 1:length(a)
			a[i]=simplify(a[i])
		end
	end
	return a
end
simplify(a::Array)=simplify!(deepcopy(a))
function simplify(d::Dict)
	nd=Dict()
	for k in keys(d)
		nk=simplify(k)
		nval=simplify(d[k])
		nd[nk]=nval
	end
	return nd
end
function sumnum(ex::Expression)
	terms=dcterms(ex)
	nterms=Term[]
	numsum=0
	for term in terms
		prod=1
		allnum=true
		nterm=Factor[]
		for t in term
			if typeof(t)<:Number
				prod*=t
			else
				if allnum
					allnum=false
				end
				push!(nterm,t)
			end
		end
		if prod==0
			continue
		elseif allnum
			numsum+=prod
			continue
		else
			if prod!=1
				unshift!(nterm,prod)
			end
			push!(nterms,nterm)
		end
	end
	if numsum!=0
		push!(nterms,[numsum])
	end
	if length(nterms)==1&&length(nterms[1])==1
		return nterms[1][1]
	elseif isempty(nterms)
		return 0
	else
		return expression(nterms)
	end
end
sumnum(c::Component)=typeof(c)(sumnum(getarg(c)))
sumnum(x::N)=x
function sumsym(ex::Expression)
	ap=terms(deepcopy(ex))
	nap=length(ap)
	cs=Array(Components,0)
	for add in 1:nap
		tcs=Array(Ex,0)
		coef=1
		for term in ap[add]
			if isa(term,Ex)
				push!(tcs,term)
			elseif isa(term,Number)
				coef*=term
			else
				error("Don't know how to handle $term")
			end
		end
		com=Components(tcs,coef)
		ind=indin(cs,com)
		if ind==0
			push!(cs,com)
		else
			cs[ind].coef+=com.coef
		end
	end
	ret=expression(cs)
	if isa(ret,Expression)&&length(ret.terms)==1&&length(ret.terms[1])==1
		return ret[1][1]
	else
		return ret
	end
end
sumsym(c::Component)=maketype(c,sumsym)
sumsym(x::N)=x
function sumsym!(a::Array)
	for i in 1:length(a)
		a[i]=sumsym(a[i])
	end
	return a
end
sumsym(a::Array)=sumsym!(deepcopy(a))
function findsyms(term::Array)
	syms=Dict()
	for fac in 1:length(term)
		if isa(term[fac],Symbol)
			if haskey(syms,term[fac])
				push!(syms[term[fac]],fac)
			else
				syms[term[fac]]=Any[fac]
			end
		end
	end
	return syms
end
function findsyms(ex::Expression)
	syms=Set{Symbol}()
	for term in ex
		for fac in term
			if isa(fac,Symbol)
				push!(syms,fac)
			elseif isa(fac,Expression)||isa(fac,Component)
				syms=union(syms,findsyms(fac))
			end
		end
	end
	return syms
end
findsyms(c::Component)=findsyms(getarg(c))
findsyms(s::Symbol)=Set([s])
findsyms(n::Number)=Set()
function findsyms(ex::Expression,symdic::Dict)
	syminds=Dict()
	for k in keys(symdic)
		inds=Tuple[]
		for ti in 1:length(ex)
			for ci in 1:length(ex[ti])
				if ex[ti][ci]==k
					push!(inds,(ti,ci))
				end
			end
		end
		syminds[k]=inds
	end
	return syminds
end
findsyms(C::Component,symdic::Dict)=findsyms(getarg(c),symdic)
function findcoms(term::Array)
	coms=Dict()
	for fac in 1:length(term)
		if isa(term[fac],Component)
			if haskey(coms,term[fac])
				push!(coms[term[fac]],fac)
			else
				coms[term[fac]]=Any[fac]
			end
		end
	end
	return coms
end
function hasex(symdic::Dict)
	for v in values(symdic)
		if isa(v,Ex)
			return true
		end
	end
	return false
end
function delexs(symdic::Dict)
	sd=deepcopy(symdic)
	for k in keys(symdic)
		if isa(symdic[k],Ex)
			delete!(sd,k)
		end
	end
	return sd
end
function replace!(ex::Expression,symdic::Dict)
	for ti in 1:length(ex)
		for c in 1:length(ex[ti])
			if isa(ex[ti][c],Ex)
				rep=replace(ex[ti][c],symdic)
				deleteat!(ex[ti],c)
				if isa(rep,Array)
					for r in rep
						insert!(ex[ti],c,r)
					end
				else
					insert!(ex[ti],c,rep)
				end
			end
		end
	end
	syminds=findsyms(ex,symdic)
	for tup in symdic
		sym,val=tup
		for i in syminds[sym]
			if isa(val,Array)
				deleteat!(ex[i[1]],i[2])
				for r in val
					insert!(ex[i[1]],i[2],r)
				end
			else
				ex[i[1]][i[2]]=val
			end
		end
	end
	if isa(ex,Array)
		ex=extract(expression(ex))
	end
	return componify(ex)
end
replace(ex::Expression,symdic::Dict)=replace!(deepcopy(ex),symdic)
function replace(c::Component,symdic::Dict)
	args=convert(Array{Any},getargs(c))
	for arg in 1:length(args)
		args[arg]=replace(args[arg],symdic)
	end
	return typeof(c)(args...)
end
function replace!(a::Array,symdic::Dict)
	for i in 1:length(a)
		a[i]=replace(a[i],symdic)
	end
	return a
end
replace(a::Array,symdic::Dict)=replace!(deepcopy(a),symdic)
function replace(s::Symbol,symdic::Dict)
	for tup in symdic
		sym,val=tup
		if s==sym
			return simplify(val)
		end
	end
	return s
end
replace(n::Number,symdic::Dict)=n
function evaluate(ex::Ex,symdic::Dict)
	ex=simplify(replace(simplify(ex),symdic))
	if isa(ex,Expression)&&hasex(symdic)
		symdic=delexs(symdic)
		ex=evaluate(ex,symdic)
	end
	return simplify(ex)
end
evaluate(x::Number,symdic::Dict)=x
function randeval(ex::Ex,seed=1)
	srand(seed)
	syms=findsyms(ex)
	d=Dict()
	for s in syms
		d[s]=rand()
	end
	evaluate(ex,d)
end
start(ex::Expression)=(1,terms(ex))
function next(ex::Expression,state)
	return (state[2][state[1]],(state[1]+1,state[2]))
end
done(ex::Expression,state)=state[1]>length(state[2])
function matches(ex::Expression,pat::Component)
	if length(ex)==1
		return matches(ex[1],pat)
	end
	termmds=Array{Dict}[]
	for term in ex
		push!(termmds,matches(extract(term),pat))
	end
	nterms=length(termmds)
	ndics=Integer[]
	for n in 1:nterms
		push!(ndics,length(termmds[n]))
	end
	ninc=reduce(*,ndics)-1
	validated=Dict[]
	for inc in 0:ninc
		indices=ones(Integer,nterms)
		tl=0
		for ti in 1:nterms
			tinc=floor(inc/reduce(*,ndics[1:ti-1]))%ndics[ti]
			tl+=ndics[ti]-1
			indices[ti]+=tinc
		end
		if clash(termmds[1][indices[1]],termmds[2][indices[2]])
			continue
		end
		comb=combine(termmds[1][indices[1]],termmds[2][indices[2]])
		clashed=false
		for tm in 3:length(termmds)
			if clash(comb,termmds[tm][indices[tm]])
				clashed=true
				break
			end
			comb=combine(comb,termmds[tm][indices[tm]])
		end
		if clashed
			continue
		end
		push!(validated,comb)
	end
	return validated
end
function pushall!(a::Array,b::Array)
	for c in b
		push!(a,c)
	end
	a
end
pushall!(a::Array,b)=push!(a,b)


text:
using Base.Test
using QuantumOptics

N = 500
xmin = -62.5
xmax = 70.1

basis_position = PositionBasis(xmin, xmax, N)
basis_momentum = MomentumBasis(basis_position)

# Test Gaussian wave-packet in both bases
x0 = 5.1
p0 = -3.2
sigma = 1.
sigma_x = sigma/sqrt(2)
sigma_p = 1./(sigma*sqrt(2))

psix0 = gaussianstate(basis_position, x0, p0, sigma)
psip0 = gaussianstate(basis_momentum, x0, p0, sigma)

@test_approx_eq 1.0 norm(psix0)
@test_approx_eq 1.0 norm(psip0)

opx_p = momentumoperator(basis_position)
opx_x = positionoperator(basis_position)

opp_p = momentumoperator(basis_momentum)
opp_x = positionoperator(basis_momentum)

@test_approx_eq x0 expect(opx_x, psix0)
@test_approx_eq_eps p0 expect(opx_p, psix0) 0.1
@test_approx_eq p0 expect(opp_p, psip0)
@test_approx_eq_eps x0 expect(opp_x, psip0) 0.1

# Test position and momentum operators
@test_approx_eq_eps expect(-opx_p^2, psix0) expect(laplace_x(basis_position), psix0) 0.3
@test_approx_eq expect(opp_p^2, psip0) expect(laplace_x(basis_momentum), psip0)
@test_approx_eq expect(opx_x^2, psix0) expect(laplace_p(basis_position), psix0)
@test_approx_eq_eps expect(-opp_x^2, psip0) expect(laplace_p(basis_momentum), psip0) 0.3

# Test FFT transformation
function transformation(b1::MomentumBasis, b2::PositionBasis, psi::Ket)
    Lp = (b1.pmax - b1.pmin)
    dx = particle.spacing(b2)
    if b1.N != b2.N || abs(2*pi/dx - Lp)/Lp > 1e-12
        throw(IncompatibleBases())
    end
    N = b1.N
    psi_shifted = exp(1im*b2.xmin*(particle.samplepoints(b1)-b1.pmin)).*psi.data
    psi_fft = exp(1im*b1.pmin*particle.samplepoints(b2)).*ifft(psi_shifted)*sqrt(N)
    return Ket(b2, psi_fft)
end

function transformation(b1::PositionBasis, b2::MomentumBasis, psi::Ket)
    Lx = (b1.xmax - b1.xmin)
    dp = particle.spacing(b2)
    if b1.N != b2.N || abs(2*pi/dp - Lx)/Lx > 1e-12
        throw(IncompatibleBases())
    end
    N = b1.N
    psi_shifted = exp(-1im*b2.pmin*(particle.samplepoints(b1)-b1.xmin)).*psi.data
    psi_fft = exp(-1im*b1.xmin*particle.samplepoints(b2)).*fft(psi_shifted)/sqrt(N)
    return Ket(b2, psi_fft)
end

psix0_fft = transformation(basis_position, basis_momentum, psix0)
psip0_fft = transformation(basis_momentum, basis_position, psip0)

@test_approx_eq_eps x0 expect(opp_x, psix0_fft) 0.1
@test_approx_eq p0 expect(opp_p, psix0_fft)
@test_approx_eq_eps p0 expect(opx_p, psip0_fft) 0.1
@test_approx_eq x0 expect(opx_x, psip0_fft)

@test_approx_eq_eps 0. norm(psix0_fft - psip0) 1e-12
@test_approx_eq_eps 0. norm(psip0_fft - psix0) 1e-12

Tpx = particle.FFTOperator(basis_momentum, basis_position)
Txp = particle.FFTOperator(basis_position, basis_momentum)
psix0_fft = Tpx*psix0
psip0_fft = Txp*psip0

@test_approx_eq_eps x0 expect(opp_x, psix0_fft) 0.1
@test_approx_eq p0 expect(opp_p, psix0_fft)
@test_approx_eq_eps p0 expect(opx_p, psip0_fft) 0.1
@test_approx_eq x0 expect(opx_x, psip0_fft)

@test_approx_eq_eps 0. norm(psix0_fft - psip0) 1e-12
@test_approx_eq_eps 0. norm(psip0_fft - psix0) 1e-12

@test_approx_eq_eps 0. norm(dagger(Tpx*psix0) - dagger(psix0)*dagger(Tpx)) 1e-12
@test_approx_eq_eps 0. norm(dagger(Txp*psip0) - dagger(psip0)*dagger(Txp)) 1e-12

psi_ = deepcopy(psip0)
operators.gemv!(Complex(1.), Tpx, psix0, Complex(0.), psi_)
@test_approx_eq_eps 0. norm(psi_ - psip0) 1e-12
operators.gemv!(Complex(1.), Tpx, psix0, Complex(1.), psi_)
@test_approx_eq_eps 0. norm(psi_ - 2*psip0) 1e-12

psi_ = deepcopy(psix0)
operators.gemv!(Complex(1.), Txp, psip0, Complex(0.), psi_)
@test_approx_eq_eps 0. norm(psi_ - psix0) 1e-12
operators.gemv!(Complex(1.), Txp, psip0, Complex(1.), psi_)
@test_approx_eq_eps 0. norm(psi_ - 2*psix0) 1e-12

rhox0x0 = tensor(psix0, dagger(psix0))
rhox0p0 = tensor(psix0, dagger(psip0))
rhop0x0 = tensor(psip0, dagger(psix0))
rhop0p0 = tensor(psip0, dagger(psip0))

rho_ = DenseOperator(basis_momentum, basis_position)
operators.gemm!(Complex(1.), Tpx, rhox0x0, Complex(0.), rho_)
@test_approx_eq_eps 0. tracedistance(rho_, rhop0x0) 1e-5
@test_approx_eq_eps 0. tracedistance(Tpx*rhox0x0, rhop0x0) 1e-5

rho_ = DenseOperator(basis_position, basis_momentum)
operators.gemm!(Complex(1.), rhox0x0, Txp, Complex(0.), rho_)
@test_approx_eq_eps 0. tracedistance(rho_, rhox0p0) 1e-5
@test_approx_eq_eps 0. tracedistance(rhox0x0*Txp, rhox0p0) 1e-5

rho_ = DenseOperator(basis_momentum, basis_momentum)
operators.gemm!(Complex(1.), Tpx, rhox0p0, Complex(0.), rho_)
@test_approx_eq_eps 0. tracedistance(rho_, rhop0p0) 1e-5
@test_approx_eq_eps 0. tracedistance(Tpx*rhox0x0*Txp, rhop0p0) 1e-5

rho_ = DenseOperator(basis_momentum, basis_momentum)
operators.gemm!(Complex(1.), rhop0x0, Txp, Complex(0.), rho_)
@test_approx_eq_eps 0. tracedistance(rho_, rhop0p0) 1e-5
@test_approx_eq_eps 0. tracedistance(Txp*rhop0p0*Tpx, rhox0x0) 1e-5

# Test FFT with lazy product
psi_ = deepcopy(psix0)
operators.gemv!(Complex(1.), LazyProduct(Txp, Tpx), psix0, Complex(0.), psi_)
@test_approx_eq_eps 0. norm(psi_ - psix0) 1e-12
@test_approx_eq_eps 0. norm(Txp*(Tpx*psix0) - psix0) 1e-12

psi_ = deepcopy(psix0)
I = dense_identity(basis_momentum)
operators.gemv!(Complex(1.), LazyProduct(Txp, I, Tpx), psix0, Complex(0.), psi_)
@test_approx_eq_eps 0. norm(psi_ - psix0) 1e-12
@test_approx_eq_eps 0. norm(Txp*I*(Tpx*psix0) - psix0) 1e-12

# Test dense FFT operator
Txp_dense = DenseOperator(Txp)
Tpx_dense = DenseOperator(Tpx)
@test typeof(Txp_dense) == DenseOperator
@test typeof(Tpx_dense) == DenseOperator
@test_approx_eq_eps 0. tracedistance(Txp_dense*rhop0p0*Tpx_dense, rhox0x0) 1e-5


text:
using Nuklear
using Nuklear.LibNuklear
using Nuklear.GLFWBackend
using GLFW, ModernGL

const WINDOW_WIDTH = 1200
const WINDOW_HEIGHT = 800

const MAX_VERTEX_BUFFER = 512 * 1024
const MAX_ELEMENT_BUFFER = 128 * 1024

@static if Sys.isapple()
    const VERSION_MAJOR = 3
    const VERSION_MINOR = 3
end

@static if Sys.isapple()
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, VERSION_MAJOR)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, VERSION_MINOR)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)
else
    GLFW.DefaultWindowHints()
end

# set up GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

# create window
win = GLFW.CreateWindow(WINDOW_WIDTH, WINDOW_HEIGHT, "Demo")
@assert win != C_NULL
GLFW.MakeContextCurrent(win)
glViewport(0, 0, WINDOW_WIDTH, WINDOW_HEIGHT)

# init context
ctx = nk_glfw3_init(win, NK_GLFW3_INSTALL_CALLBACKS, MAX_VERTEX_BUFFER, MAX_ELEMENT_BUFFER)

# init font
prefix = joinpath(@__DIR__, "extra_font")
nk_glfw3_font_stash_begin()
roboto = nk_font_atlas_add_from_file(glfw.atlas, joinpath(prefix, "Roboto-Bold.ttf"), 18, C_NULL)
nk_glfw3_font_stash_end()
roboto_handle = nk_get_font_handle(roboto)
nk_style_set_font(ctx, roboto_handle)

# states
const EASY = 0
const HARD = 1

op = EASY
property = Ref{Cint}(20)
bg = nk_colorf(0.10, 0.18, 0.24, 1.0)
while !GLFW.WindowShouldClose(win)
    global op
    global bg
    global ctx

    GLFW.PollEvents()
    nk_glfw3_new_frame()

    if Bool(nk_begin(ctx, "Demo", nk_rect(50, 50, 230, 250),
                     NK_WINDOW_BORDER|NK_WINDOW_MOVABLE|NK_WINDOW_SCALABLE|
                     NK_WINDOW_MINIMIZABLE|NK_WINDOW_TITLE))

        nk_layout_row_static(ctx, 30, 80, 1)
        nk_button_label(ctx, "button") == 1 && @show("button pressed")

        nk_layout_row_dynamic(ctx, 30, 2)
        Bool(nk_option_label(ctx, "easy", op == EASY)) && (op = EASY;)
        Bool(nk_option_label(ctx, "hard", op == HARD)) && (op = HARD;)

        nk_layout_row_dynamic(ctx, 25, 1)
        nk_property_int(ctx, "Compression:", 0, property, 100, 10, 1)

        nk_layout_row_dynamic(ctx, 20, 1)
        nk_label(ctx, "background:", NK_TEXT_LEFT)
        nk_layout_row_dynamic(ctx, 25, 1)
        if Bool(nk_combo_begin_color(ctx, nk_rgb_cf(bg), nk_vec2(nk_widget_width(ctx),400)))
            nk_layout_row_dynamic(ctx, 120, 1)
            bg = nk_color_picker(ctx, bg, NK_RGBA)
            nk_layout_row_dynamic(ctx, 25, 1)
            r = nk_propertyf(ctx, "#R:", 0, bg.r, 1.0, 0.01, 0.005)
            g = nk_propertyf(ctx, "#G:", 0, bg.g, 1.0, 0.01, 0.005)
            b = nk_propertyf(ctx, "#B:", 0, bg.b, 1.0, 0.01, 0.005)
            a = nk_propertyf(ctx, "#A:", 0, bg.a, 1.0, 0.01, 0.005)
            bg = nk_colorf(r, g, b, a)
            nk_combo_end(ctx)
        end
    end
    nk_end(ctx)

    # draw
    width, height = GLFW.GetWindowSize(win)
    glViewport(0, 0, width, height)
    glClear(GL_COLOR_BUFFER_BIT)
    glClearColor(bg.r, bg.g, bg.b, bg.a)
    nk_glfw3_render(NK_ANTI_ALIASING_ON, MAX_VERTEX_BUFFER, MAX_ELEMENT_BUFFER)
    GLFW.SwapBuffers(win)
end

nk_glfw3_shutdown()
GLFW.DestroyWindow(win)