text:
test01_eigvals = [
    -0.15883100381194412 - 13.264252860693972im
    -0.15883100381194412 + 13.264252860693972im
]

test01_eigvals_psat = [
    -0.15883 - 13.2643im
    -0.15883 + 13.2643im
]

test02_eigvals = [
    -999.9999922288165 + 0.0im
    -999.9999824863403 + 0.0im
    -5.667091663572526 + 0.0im
    -4.939539293199701 - 7.861631594468724im
    -4.939539293199701 + 7.861631594468724im
    -4.653008913694498 + 0.0im
    -4.144476689106577 - 7.59715123723563im
    -4.144476689106577 + 7.59715123723563im
    -1.3798285681206308 - 16.03637615711796im
    -1.3798285681206308 + 16.03637615711796im
    -0.9405263132225796 - 13.949752815351614im
    -0.9405263132225796 + 13.949752815351614im
    -0.4805073944481576 - 0.5379518986751338im
    -0.4805073944481576 + 0.5379518986751338im
    -0.27347592115017016 - 0.4565372544427486im
    -0.27347592115017016 + 0.4565372544427486im
]

test02_eigvals_psat = [
    -999.999990000000 + 0.00000000000000im
    -999.999980000000 + 0.00000000000000im
    -5.67379000000000 + 0.00000000000000im
    -4.93393000000000 - 7.85614000000000im
    -4.93393000000000 + 7.85614000000000im
    -4.64127000000000 + 0.00000000000000im
    -4.14297000000000 - 7.59557000000000im
    -4.14297000000000 + 7.59557000000000im
    -1.36440000000000 - 14.5583400000000im
    -1.36440000000000 + 14.5583400000000im
    -0.967800000000000 - 12.6558000000000im
    -0.967800000000000 + 12.6558000000000im
    -0.471260000000000 - 0.579790000000000im
    -0.471260000000000 + 0.579790000000000im
    -0.283050000000000 - 0.460360000000000im
    -0.283050000000000 + 0.460360000000000im
]

test03_eigvals = [
    -1000.0000000405296 + 0.0im
    -1000.0000000190737 + 0.0im
    -56.28815836405992 + 0.0im
    -53.12282741877609 + 0.0im
    -5.957028010206276 + 0.0im
    -5.02141964268554 - 7.898682332162664im
    -5.02141964268554 + 7.898682332162664im
    -4.8668814493140715 + 0.0im
    -4.218933906222015 - 7.617379500840438im
    -4.218933906222015 + 7.617379500840438im
    -2.117324030780644 + 0.0im
    -2.029229441492322 + 0.0im
    -1.7025637072900843 - 16.974819568195652im
    -1.7025637072900843 + 16.974819568195652im
    -1.1756081480320293 - 14.621543797706247im
    -1.1756081480320293 + 14.621543797706247im
    -0.3806667287614366 - 0.5866699857779792im
    -0.3806667287614366 + 0.5866699857779792im
    -0.223947575088223 - 0.4971805954842863im
    -0.223947575088223 + 0.4971805954842863im
]

test03_eigvals_psat = [
    -1000.00000000000 + 0.00000000000000im
    -1000.00000000000 + 0.00000000000000im
    -56.4167900000000 + 0.00000000000000im
    -53.3586000000000 + 0.00000000000000im
    -5.95708000000000 + 0.00000000000000im
    -5.02325000000000 - 7.89808000000000im
    -5.02325000000000 + 7.89808000000000im
    -4.82461000000000 + 0.00000000000000im
    -4.21943000000000 - 7.61745000000000im
    -4.21943000000000 + 7.61745000000000im
    -2.12296000000000 + 0.00000000000000im
    -2.03127000000000 + 0.00000000000000im
    -1.63531000000000 - 15.4006100000000im
    -1.63531000000000 + 15.4006100000000im
    -1.19151000000000 - 13.1976800000000im
    -1.19151000000000 + 13.1976800000000im
    -0.351320000000000 - 0.628360000000000im
    -0.351320000000000 + 0.628360000000000im
    -0.237770000000000 - 0.501610000000000im
    -0.237770000000000 + 0.501610000000000im
]

test04_eigvals = [
    -1000.0000000605402 + 0.0im
    -1000.0000000226693 + 0.0im
    -57.04885681143213 - 637.6983559917089im
    -57.04885681143213 + 637.6983559917089im
    -56.354194606909395 + 0.0im
    -53.18155023445317 + 0.0im
    -13.892570544392896 - 469.6478362124919im
    -13.892570544392896 + 469.6478362124919im
    -5.958255480381696 + 0.0im
    -5.021533592090426 - 7.898654827966522im
    -5.021533592090426 + 7.898654827966522im
    -4.868492180624505 + 0.0im
    -4.21902585809358 - 7.617326092491787im
    -4.21902585809358 + 7.617326092491787im
    -2.117386441281952 + 0.0im
    -2.029241696160179 + 0.0im
    -1.6959817848903833 - 16.973813636829668im
    -1.6959817848903833 + 16.973813636829668im
    -1.1682442396700994 - 14.616878713067868im
    -1.1682442396700994 + 14.616878713067868im
    -0.3805794727021142 - 0.586700331651268im
    -0.3805794727021142 + 0.586700331651268im
    -0.22394101115665865 - 0.49721792171561435im
    -0.22394101115665865 + 0.49721792171561435im
]

test05_eigvals = [
    -1000.0000037690621 + 0.0im
    -1000.000001656506 + 0.0im
    -126.28458815493319 + 0.0im
    -120.85493880185086 + 0.0im
    -46.48594631374737 + 0.0im
    -44.373301051079025 + 0.0im
    -10.0 + 0.0im
    -10.0 + 0.0im
    -3.9151560397679868 - 7.4523722134285535im
    -3.9151560397679868 + 7.4523722134285535im
    -3.836126674242972 - 7.360635983756108im
    -3.836126674242972 + 7.360635983756108im
    -0.3447395139088272 - 12.01547546088681im
    -0.3447395139088272 + 12.01547546088681im
    -0.20916948755870068 - 10.607466623279306im
    -0.20916948755870068 + 10.607466623279306im
    -0.16185759374101516 - 1.055686298197309im
    -0.16185759374101516 + 1.055686298197309im
    -0.0354521777174975 - 0.684452609868359im
    -0.0354521777174975 + 0.684452609868359im
]

test06_eigvals = [
    -1000.0000053062265 + 0.0im
    -1000.0000018990927 + 0.0im
    -126.75458543123136 + 0.0im
    -121.57403285596196 + 0.0im
    -46.52052665063294 + 0.0im
    -44.44016429276653 + 0.0im
    -34.69415568906223 - 514.6812519452002im
    -34.69415568906223 + 514.6812519452002im
    -10.0 + 0.0im
    -10.0 + 0.0im
    -9.032779692453872 - 426.00371528493883im
    -9.032779692453872 + 426.00371528493883im
    -3.915125641859613 - 7.452252518206685im
    -3.915125641859613 + 7.452252518206685im
    -3.8358938132681137 - 7.35974852776193im
    -3.8358938132681137 + 7.35974852776193im
    -0.3416921618276846 - 12.015383170233138im
    -0.3416921618276846 + 12.015383170233138im
    -0.2070154019752105 - 10.607728100520474im
    -0.2070154019752105 + 10.607728100520474im
    -0.16160669596957933 - 1.0557462673615183im
    -0.16160669596957933 + 1.0557462673615183im
    -0.03542058981778998 - 0.6844527516049774im
    -0.03542058981778998 + 0.6844527516049774im
]

test07_eigvals = [
    -999.9999914519317 + 0.0im
    -999.9999811538516 + 0.0im
    -5.027556777415359 + 0.0im
    -4.718819382076426 - 7.793357006494638im
    -4.718819382076426 + 7.793357006494638im
    -4.102660291475314 + 0.0im
    -4.052932897118745 - 7.552138242722714im
    -4.052932897118745 + 7.552138242722714im
    -1.786149751424446 - 14.793097612110634im
    -1.786149751424446 + 14.793097612110634im
    -1.2196573076881576 - 10.886078373690756im
    -1.2196573076881576 + 10.886078373690756im
    -0.8706541640279326 - 305.1344520297189im
    -0.8706541640279326 + 305.1344520297189im
    -0.5028032785077488 - 0.5630005288071646im
    -0.5028032785077488 + 0.5630005288071646im
    -0.29876013022174563 - 0.4533947142923075im
    -0.29876013022174563 + 0.4533947142923075im
    -0.263305765310502 - 177.9037730410636im
    -0.263305765310502 + 177.9037730410636im
    -0.1738090829698704 - 62.759642648170946im
    -0.1738090829698704 + 62.759642648170946im
    -0.15473139663326663 - 136.7874068069445im
    -0.15473139663326663 + 136.7874068069445im
]

test08_eigvals = [
    -2285.9589650663725 - 6825.3997886670795im
    -2285.9589650663725 + 6825.3997886670795im
    -2095.818014258244 - 6533.140681149019im
    -2095.818014258244 + 6533.140681149019im
    -1595.4562252195476 - 233.8877642497275im
    -1595.4562252195476 + 233.8877642497275im
    -986.7365548177571 + 0.0im
    -500.00000000000034 + 0.0im
    -471.8999089077167 + 0.0im
    -220.89110522830256 + 0.0im
    -50.32121120551376 + 0.0im
    -50.245608339393456 + 0.0im
    -34.244848184557476 - 260.53071319880445im
    -34.244848184557476 + 260.53071319880445im
    -11.303322380605753 + 0.0im
    -11.19036178279552 + 0.0im
    -7.691789694217273 - 28.802262759948757im
    -7.691789694217273 + 28.802262759948757im
    -4.513698394541721 + 0.0im
]

test09_eigvals = [
    -2266.7564379695177 - 6842.581227177472im
    -2266.7564379695177 + 6842.581227177472im
    -2102.099837452439 - 6516.907948869909im
    -2102.099837452439 + 6516.907948869909im
    -1603.381787508888 - 299.73078962017183im
    -1603.381787508888 + 299.73078962017183im
    -999.9999853200115 + 0.0im
    -981.5890587115106 + 0.0im
    -499.9999999999994 + 0.0im
    -472.94184918788915 + 0.0im
    -222.99856689050083 + 0.0im
    -58.916888833164734 - 297.3751281905745im
    -58.916888833164734 + 297.3751281905745im
    -50.391739774950395 + 0.0im
    -50.28496360614687 + 0.0im
    -11.312374832548043 + 0.0im
    -11.192736801973947 + 0.0im
    -6.304389558457761 - 25.99314750112348im
    -6.304389558457761 + 25.99314750112348im
    -4.936441589468676 - 7.860345016545073im
    -4.936441589468676 + 7.860345016545073im
    -4.874648728527242 + 0.0im
    -4.628324873366437 + 0.0im
    -1.3057837545018711 - 14.985883650327404im
    -1.3057837545018711 + 14.985883650327404im
    -0.522548405283088 - 0.5356190189108261im
    -0.522548405283088 + 0.5356190189108261im
]

test10_eigvals = [
    -2272.361404997375 - 6851.776097141858im
    -2272.361404997375 + 6851.776097141858im
    -2124.253409194327 - 6511.932840117861im
    -2124.253409194327 + 6511.932840117861im
    -1575.608572024949 - 348.98204359994224im
    -1575.608572024949 + 348.98204359994224im
    -999.9999861307571 + 0.0im
    -960.833224373504 + 0.0im
    -499.9999999999994 + 0.0im
    -477.8792003131099 + 0.0im
    -218.19978493862791 + 0.0im
    -75.81020867999811 - 305.1528459259956im
    -75.81020867999811 + 305.1528459259956im
    -50.41186110201852 + 0.0im
    -50.352254884885475 + 0.0im
    -11.32379102275944 + 0.0im
    -11.22342662568646 + 0.0im
    -7.581779910900623 - 27.826786728047516im
    -7.581779910900623 + 27.826786728047516im
    -5.484429063458374 + 0.0im
    -3.9465322340765994 - 7.508941380094648im
    -3.9465322340765994 + 7.508941380094648im
    -2.907232875387719 + 0.0im
    -0.7994143393863338 - 14.039239192684892im
    -0.7994143393863338 + 14.039239192684892im
    -0.13279249045498032 - 0.4749686774987308im
    -0.13279249045498032 + 0.4749686774987308im
]

test11_eigvals = [
    -2280.2076336006976 - 6841.438512443791im
    -2280.2076336006976 + 6841.438512443791im
    -1997.9475701466922 - 6457.618671478096im
    -1997.9475701466922 + 6457.618671478096im
    -1538.5966307902208 - 345.7913984886787im
    -1538.5966307902208 + 345.7913984886787im
    -1250.2936918022115 - 9608.999772834375im
    -1250.2936918022115 + 9608.999772834375im
    -999.9999874645969 + 0.0im
    -956.6099185810597 + 0.0im
    -863.7185028125181 - 6165.596120417795im
    -863.7185028125181 + 6165.596120417795im
    -500.0 + 0.0im
    -478.235862371146 + 0.0im
    -218.28706882452354 + 0.0im
    -58.31526552633564 - 299.61738004327435im
    -58.31526552633564 + 299.61738004327435im
    -50.414620572174904 + 0.0im
    -50.35016783858917 + 0.0im
    -42.437976630212574 - 467.0560264577345im
    -42.437976630212574 + 467.0560264577345im
    -11.323801014242834 + 0.0im
    -11.223449045909994 + 0.0im
    -7.577419561243316 - 27.833300571016597im
    -7.577419561243316 + 27.833300571016597im
    -5.483728552089268 + 0.0im
    -3.9464440004154135 - 7.508758038269464im
    -3.9464440004154135 + 7.508758038269464im
    -2.9069124155527093 + 0.0im
    -0.8024184353347599 - 14.03943882277388im
    -0.8024184353347599 + 14.03943882277388im
    -0.13278021568161175 - 0.4749742367984249im
    -0.13278021568161175 + 0.4749742367984249im
]

test12_eigvals = [
    -0.8055565836807297 - 0.5525321844563806im,
    -0.8055565836807297 + 0.5525321844563806im,
    -0.5562604430028794 - 13.482773238020515im,
    -0.5562604430028794 + 13.482773238020515im,
    0.0 + 0.0im,
]

test13_eigvals = [
    -10459.358851638417 + 0.0im
    -9518.115670592897 + 0.0im
    -187.23828129745846 + 0.0im
    -10.106065550045505 + 0.0im
    -5.224596052187488 + 0.0im
    -4.501324080186281 + 0.0im
    -2.2096561921652595 + 0.0im
    -1.364291810141959 - 16.041013087225934im
    -1.364291810141959 + 16.041013087225934im
    -0.9689898367154597 - 0.5406190107493515im
    -0.9689898367154597 + 0.5406190107493515im
    -0.912054975024942 - 14.019102730055716im
    -0.912054975024942 + 14.019102730055716im
    -0.7804036299927885 - 3.4337139420134424im
    -0.7804036299927885 + 3.4337139420134424im
    -0.01997649074097721 + 0.0im
]

test14_eigvals = [
    -2211.5592709440352 - 6867.397622238308im,
    -2211.5592709440352 + 6867.397622238308im,
    -2153.3431924567208 - 6491.972141691999im,
    -2153.3431924567208 + 6491.972141691999im,
    -1462.9335729586219 - 695.7972901292562im,
    -1462.9335729586219 + 695.7972901292562im,
    -970.0160993678469 + 0.0im,
    -500.00000000000045 + 0.0im,
    -490.8250468494144 + 0.0im,
    -202.2656322427033 - 420.17747462821575im,
    -202.2656322427033 + 420.17747462821575im,
    -189.55369688202828 + 0.0im,
    -50.54441602439222 + 0.0im,
    -50.45006132952855 + 0.0im,
    -17.63867329949132 - 46.366884762058845im,
    -17.63867329949132 + 46.366884762058845im,
    -11.355851160484914 + 0.0im,
    -11.335928224977632 + 0.0im,
    -1.4563561997747976 - 8.058315342439121im,
    -1.4563561997747976 + 8.058315342439121im,
    -0.5056731355961069 + 0.0im,
    0.0 + 0.0im,
]

test15_eigvals = [
    -39.22517135042766 - 0.23707459763811717im
    -39.22517135042766 + 0.23707459763811717im
    -8.162038320970622 + 0.0im
    -0.6894464243695174 - 8.596453673443301im
    -0.6894464243695174 + 8.596453673443301im
    -0.20829262548711688 + 0.0im
]

test15_eigvals_no_sat = [
    -41.28060502979358 + 0.0im
    -38.77124317334057 + 0.0im
    -6.080781526447125 + 0.0im
    -0.6716198831075708 - 8.713016722797422im
    -0.6716198831075708 + 8.713016722797422im
    -0.12135746559823009 + 0.0im
]

test15_eigvals_high_sat = [
    -38.90919918657826 - 0.5990898792027552im
    -38.90919918657826 + 0.5990898792027552im
    -8.716992122605445 + 0.0im
    -0.7802024926568868 - 8.344701265008863im
    -0.7802024926568868 + 8.344701265008863im
    -0.34108463674219724 + 0.0im
]

test16_eigvals = [
    -40.86489059167709 + 0.0im
    -38.79830220459659 + 0.0im
    -6.574706325237592 + 0.0im
    -0.6823066948004968 - 8.675007074248816im
    -0.6823066948004968 + 8.675007074248816im
    -0.14426163605639922 + 0.0im
]

test16_eigvals_high_sat = [
    -40.47662880667318 + 0.0im
    -38.83801972051611 + 0.0im
    -7.051748409012963 + 0.0im
    -0.6774433789575647 - 8.665077736802273im
    -0.6774433789575647 + 8.665077736802273im
    -0.15891752736421258 + 0.0im
]

test17_eigvals = [
    -999.9999972744532 + 0.0im
    -41.280911913534204 + 0.0im
    -38.76959077602596 + 0.0im
    -6.091730762853858 + 0.0im
    -4.397466628592957 - 7.701876240796052im
    -4.397466628592957 + 7.701876240796052im
    -0.667225952071195 - 8.709471358908846im
    -0.667225952071195 + 8.709471358908846im
    -0.19039981492553382 - 0.3625458335906406im
    -0.19039981492553382 + 0.3625458335906406im
]

test18_eigvals = [
    -28.693211883382347 + 0.0im
    -9.84238486492247 + 0.0im
    -1.7031770286738355 - 8.629174022229318im
    -1.7031770286738355 + 8.629174022229318im
    -0.7640548534736492 + 0.0im
]

test19_eigvals = [
    -28.76860514709647 + 0.0im
    -11.070928271739678 + 0.0im
    -1.665595867723869 - 8.572069631777202im
    -1.665595867723869 + 8.572069631777202im
    -0.4928568982924553 + 0.0im
]

test20_eigvals = [
    -100.00021842764602 + 0.0im
    -39.23707786473363 - 0.30958964511252446im
    -39.23707786473363 + 0.30958964511252446im
    -22.649793103911897 + 0.0im
    -8.543255967128875 + 0.0im
    -4.268315666497509 - 5.271829751631857im
    -4.268315666497509 + 5.271829751631857im
    -0.6848864366573055 - 8.58151201322762im
    -0.6848864366573055 + 8.58151201322762im
    -0.3129633046035981 - 0.5454321262659743im
    -0.3129633046035981 + 0.5454321262659743im
]

test20_eigvals_sat = [
    -100.00021924248418 + 0.0im
    -39.237422024473375 - 0.3111853351055383im
    -39.237422024473375 + 0.3111853351055383im
    -23.53165375119335 + 0.0im
    -8.410026826125367 + 0.0im
    -4.083256130308611 - 6.7290588685710615im
    -4.083256130308611 + 6.7290588685710615im
    -0.6809735490286486 - 8.581164994672914im
    -0.6809735490286486 + 8.581164994672914im
    -0.2887994589525554 - 0.4637034811911777im
    -0.2887994589525554 + 0.4637034811911777im
]

test21_eigvals = [
    -100.00020053262892 + 0.0im
    -41.279576521715754 + 0.0im
    -38.7876974005338 + 0.0im
    -22.69484347984886 + 0.0im
    -7.165287493953102 + 0.0im
    -6.28758745738531 + 0.0im
    -4.319362829032215 - 5.205827008438781im
    -4.319362829032215 + 5.205827008438781im
    -1.6673562530234518 - 1.8236637052912814im
    -1.6673562530234518 + 1.8236637052912814im
    -0.7102160509747345 - 8.694108632428092im
    -0.7102160509747345 + 8.694108632428092im
    -0.2846615017660007 - 0.607222973124181im
    -0.2846615017660007 + 0.607222973124181im
]

test22_eigvals = [
    -100.00020053242672 + 0.0im
    -41.2798088140291 + 0.0im
    -38.788275874605674 + 0.0im
    -22.69485159080186 + 0.0im
    -6.302325814485422 + 0.0im
    -4.931267158171579 + 0.0im
    -4.31725705837116 - 5.200834867430221im
    -4.31725705837116 + 5.200834867430221im
    -1.2329330244784797 + 0.0im
    -0.7463112762985753 - 8.85884300731678im
    -0.7463112762985753 + 8.85884300731678im
    -0.28579333866010825 - 0.6025484098882521im
    -0.28579333866010825 + 0.6025484098882521im
]

test23_eigvals = [
    -2285.984972487398 - 6825.409970531314im
    -2285.984972487398 + 6825.409970531314im
    -2095.8292998741636 - 6533.139948007082im
    -2095.8292998741636 + 6533.139948007082im
    -1596.153193009681 - 233.6009166562695im
    -1596.153193009681 + 233.6009166562695im
    -986.970453826522 + 0.0im
    -50.61566828819503 + 0.0im
    -50.24422063249634 + 0.0im
    -33.8405217631388 - 255.4659930547156im
    -33.8405217631388 + 255.4659930547156im
    -16.35252074296158 - 35.19268456048173im
    -16.35252074296158 + 35.19268456048173im
    -11.402718862758675 + 0.0im
    -11.30330507373553 + 0.0im
]

test24_eigvals = [
    -16167.452158969038 - 8.043416230174378im
    -16167.452158969038 + 8.043416230174378im
    -306.82201246319016 - 5867.887293375189im
    -306.82201246319016 + 5867.887293375189im
    -258.57819319671626 - 4795.264357598769im
    -258.57819319671626 + 4795.264357598769im
    -202.1765718936583 - 520.5322326644553im
    -202.1765718936583 + 520.5322326644553im
    -26.89823743203937 - 74.8261131632425im
    -26.89823743203937 + 74.8261131632425im
    -14.146571409094378 - 5.185851435997753im
    -14.146571409094378 + 5.185851435997753im
    -10.131535087046432 + 0.0im
    -1.8914933786938954 - 0.005547477442692941im
    -1.8914933786938954 + 0.005547477442692941im
]

test26_eigvals = [
    -39.22608968461495 - 0.24435922038292715im
    -39.22608968461495 + 0.24435922038292715im
    -8.153259227989246 + 0.0im
    -0.9387139673515449 + 0.0im
    -0.6886612978044744 - 8.596257367969528im
    -0.6886612978044744 + 8.596257367969528im
    -0.2390456681797083 - 0.17819776543201055im
    -0.2390456681797083 + 0.17819776543201055im
]

test26_eigvals_noTE = [
    -39.47393476296245 + 0.0im
    -38.90941405619594 + 0.0im
    -8.225656744581503 + 0.0im
    -0.6876080511263236 - 8.603146074166037im
    -0.6876080511263236 + 8.603146074166037im
    -0.24592182992188777 - 0.14919699127432942im
    -0.24592182992188777 + 0.14919699127432942im
]

test29_eigvals = [
    -59.7320535837085 + 0.0im
    -50.00000000000001 + 0.0im
    -50.0 + 0.0im
    -50.0 + 0.0im
    -45.133973208145804 - 8.692497819063023im
    -45.133973208145804 + 8.692497819063023im
    -20.0 + 0.0im
    -5.0 + 0.0im
    0.0 + 0.0im
]

test29_eigvals_fflag = [
    -59.772220372205 + 0.0im
    -51.0447940895546 - 7.0658041406756436im
    -51.0447940895546 + 7.0658041406756436im
    -50.0 + 0.0im
    -45.18265590251619 - 8.61499070758468im
    -45.18265590251619 + 8.61499070758468im
    -20.0 + 0.0im
    -5.873878451784012 - 6.347167941328094im
    -5.873878451784012 + 6.347167941328094im
    -5.0 + 0.0im
    -0.02512274008534303 + 0.0im
    0.0 + 0.0im
]

test29_eigvals_Rflag = [
    -76.91753318136166 - 32.2358556345301im
    -76.91753318136166 + 32.2358556345301im
    -50.0 + 0.0im
    -49.999999999999986 - 9.959757193250047e-8im
    -49.999999999999986 + 9.959757193250047e-8im
    -7.881077881909685 - 31.52570878189592im
    -7.881077881909685 + 31.52570878189592im
    -5.0 + 0.0im
    -0.40277787345725297 + 0.0im
    0.0 + 0.0im
]

test29_eigvals_Rflag_no_Vflag = [
    -76.91753318136163 - 32.235855634530125im
    -76.91753318136163 + 32.235855634530125im
    -50.0 + 0.0im
    -49.999999999999986 - 2.669704642995615e-7im
    -49.999999999999986 + 2.669704642995615e-7im
    -7.8810778819096985 - 31.525708781895904im
    -7.8810778819096985 + 31.525708781895904im
    -5.0 + 0.0im
    -0.4027778734572379 + 0.0im
]

test29_eigvals_Rflag_Qflag = [
    -50.00000000000023 + 0.0im
    -50.0 + 0.0im
    -49.99977208639873 + 0.0im
    -48.583188396862354 - 43.26816561002791im
    -48.583188396862354 + 43.26816561002791im
    -19.851972845705998 + 0.0im
    -5.0 + 0.0im
    -2.616189310673934 + 0.0im
    -0.1828444817481023 - 0.361230702596121im
    -0.1828444817481023 + 0.361230702596121im
]

test30_eigvals_with_filt = [
    -76.92307692307699 + 0.0im
    -49.99999999999998 + 0.0im
    -39.22608968461492 - 0.24435922038297828im
    -39.22608968461492 + 0.24435922038297828im
    -8.15325922798925 + 0.0im
    -1.0000000000000009 + 0.0im
    -0.938713967351545 + 0.0im
    -0.6886612978044739 - 8.596257367969525im
    -0.6886612978044739 + 8.596257367969525im
    -0.6060606060606061 + 0.0im
    -0.23904566817970954 - 0.17819776543200883im
    -0.23904566817970954 + 0.17819776543200883im
    0.0 + 0.0im
]

test30_eigvals_no_filt = [
    -50.0 + 0.0im
    -39.226089684614934 - 0.24435922038295738im
    -39.226089684614934 + 0.24435922038295738im
    -8.153259227989247 + 0.0im
    -0.9387139673515453 + 0.0im
    -0.688661297804473 - 8.59625736796952im
    -0.688661297804473 + 8.59625736796952im
    -0.6060606060606061 + 0.0im
    -0.2390456681797072 - 0.17819776543201066im
    -0.2390456681797072 + 0.17819776543201066im
    0.0 + 0.0im
    0.0 + 0.0im
]

test31_eigvals = [
    -41.28028496346849 + 0.0im
    -38.772141289217586 + 0.0im
    -6.0734966129815335 + 0.0im
    -3.2959209145486295 + 0.0im
    -0.9334113090863949 + 0.0im
    -0.714334242846485 - 8.712902325176582im
    -0.714334242846485 + 8.712902325176582im
    -0.196671279155743 - 0.19929469863615853im
    -0.196671279155743 + 0.19929469863615853im
    -0.1790842443333912 + 0.0im
    -0.0491682970453784 + 0.0im
    -0.03336518379541009 + 0.0im
]

test33_eigvals_constantP = [
    -39.45761150626904 + 0.0im
    -38.80018214770829 + 0.0im
    -7.781678192463298 + 0.0im
    -0.9355834221192385 - 8.472947434030154im
    -0.9355834221192385 + 8.472947434030154im
    -0.535017710389672 + 0.0im
]

test33_eigvals_constantI = [
    -39.154101609632264 - 0.30688836000778646im
    -39.154101609632264 + 0.30688836000778646im
    -7.780866011433047 + 0.0im
    -0.9320771142348544 - 8.476045052644412im
    -0.9320771142348544 + 8.476045052644412im
    -0.536696069017935 + 0.0im
]

test33_eigvals_constantZ = [
    -39.17600477596733 - 0.5179363938174002im
    -39.17600477596733 + 0.5179363938174002im
    -7.780182413917664 + 0.0im
    -0.9290443071358351 - 8.478692979297731im
    -0.9290443071358351 + 8.478692979297731im
    -0.5381468520266435 + 0.0im
]

test37_eigvals = [
    -29.354494961308443 - 928.1009201162458im
    -29.354494961308443 + 928.1009201162458im
    -21.358584482366716 + 0.0im
    -9.446872682532018 + 0.0im
    -5.415417142234835 + 0.0im
    -3.5158065492012964 - 15.750279074009407im
    -3.5158065492012964 + 15.750279074009407im
    -3.405242872115318 - 2.228869205181562im
    -3.405242872115318 + 2.228869205181562im
    -3.1776970108237883 + 0.0im
    -1.0 + 0.0im
    -0.39349870582862 - 8.82930733740199im
    -0.39349870582862 + 8.82930733740199im
    -0.07411941090418851 + 0.0im
    -0.001992211534244684 + 0.0im
    -0.001001001001001001 + 0.0im
]

test39_eigvals = [
    -100.45338295445788 + 0.0im
    -21.35498851941608 + 0.0im
    -10.0 + 0.0im
    -5.463171684799997 + 0.0im
    -4.91399966284345 - 4.485021668563468im
    -4.91399966284345 + 4.485021668563468im
    -3.149582560233519 + 0.0im
    -0.34584314540905314 - 9.434904079322248im
    -0.34584314540905314 + 9.434904079322248im
    -0.09719979123415864 + 0.0im
    -0.07412657316639956 + 0.0im
    -0.002014246353738706 + 0.0im
    -0.001001001001001001 + 0.0im
]

test40_eigvals = [
    -99.94136467100992 + 0.0im
    -21.355217306408555 + 0.0im
    -11.231106548123972 - 9.575595271086215im
    -11.231106548123972 + 9.575595271086215im
    -5.4404844263731915 + 0.0im
    -3.157024438013972 + 0.0im
    -1.4915618476418233 - 4.836401833717598im
    -1.4915618476418233 + 4.836401833717598im
    -1.0 + 0.0im
    -0.5592551382672033 - 9.485811647183036im
    -0.5592551382672033 + 9.485811647183036im
    -0.07411245925908164 + 0.0im
    -0.0020142681327209597 + 0.0im
    -0.001001001001001001 + 0.0im
]


text:
module M

export f

f(x) = x

g(x = 1, y = 2, z = 3; kwargs...) = x

h(x::Int, y::Int = 2, z::Int = 3; kwargs...) = x

const A{T} = Union{Vector{T}, Matrix{T}}

h_1(x::A) = x
h_2(x::A{Int}) = x
h_3(x::A{T}) where {T} = x

i_1(x; y = x) = x * y
i_2(x::Int; y = x) = x * y
i_3(x::T; y = x) where {T} = x * y
i_4(x; y::T = zero(T), z::U = zero(U)) where {T, U} = x + y + z

j_1(x, y) = x * y # two arguments, no keyword arguments
j_1(x; y = x) = x * y # one argument, one keyword argument

mutable struct T
    a
    b
    c
end

struct K
    K(; a = 1) = new()
end

abstract type AbstractType <: Integer end

struct CustomType{S, T <: Integer} <: Integer
end

primitive type BitType8 8 end

primitive type BitType32 <: Real 32 end

end

ext: jl


text:
export CS_model
"""
    CS_model(df, y, grouping, d, link)

Form the compound symmetric (CS) model for intercept only regression with the specified base distribution (d) and link function (link).

# Arguments
- `df`: A named `DataFrame`
- `y`: Ouctcome variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `grouping`: Grouping or Clustering variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `d`: Base `Distribution` of outcome from `Distributions.jl`.
- `link`: Canonical `Link` function of the base distribution specified in `d`, from `GLM.jl`.

# Optional Arguments
- `penalized`: Boolean to specify whether or not to add an L2 Ridge penalty on the variance parameter for the CS structured covariance.
    One can turn this option on by specifying `penalized = true` to add this penalty for numerical stability. (default `penalized = false`).
"""
function CS_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Normal, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GaussianCopulaCSObs{Float64}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        gcs[i] = GaussianCopulaCSObs(Y, X)
    end
    gcm = GaussianCopulaCSModel(gcs; penalized = penalized);
    return gcm
end

function CS_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Union{Poisson, Bernoulli}, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GLMCopulaCSObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        gcs[i] = GLMCopulaCSObs(Y, X, d, link)
    end
    gcm = GLMCopulaCSModel(gcs; penalized = penalized);
    return gcm
end

function CS_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:NegativeBinomial, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{NBCopulaCSObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        gcs[i] = NBCopulaCSObs(Y, X, d, link)
    end
    gcm = NBCopulaCSModel(gcs; penalized = penalized);
    return gcm
end

"""
    CS_model(df, y, grouping, covariates, d, link; penalized = penalized)

Form the compound symmetric (CS) model for regression with the specified base distribution (d) and link function (link).

# Arguments
- `df`: A named `DataFrame`
- `y`: Ouctcome variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `grouping`: Grouping or Clustering variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `covariates`: Covariate names of interest as a vector of `Symbol`s.
    Each variable name must be present in `df`.
- `d`: Base `Distribution` of outcome from `Distributions.jl`.
- `link`: Canonical `Link` function of the base distribution specified in `d`, from `GLM.jl`.

# Optional Arguments
- `penalized`: Boolean to specify whether or not to add an L2 Ridge penalty on the variance parameter for the CS structured covariance.
    One can turn this option on by specifying `penalized = true` to add this penalty for numerical stability. (default `penalized = false`).
"""
function CS_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    covariates::Vector{Symbol},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Normal, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GaussianCopulaCSObs{Float64}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        @inbounds for i in 1:length(covariates)
            U = Float64.(df[gidx, covariates[i]])
            X = hcat(X, U)
        end
        gcs[i] = GaussianCopulaCSObs(Y, X)
    end
    gcm = GaussianCopulaCSModel(gcs; penalized = penalized);
    return gcm
end

function CS_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    covariates::Vector{Symbol},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Union{Poisson, Bernoulli}, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GLMCopulaCSObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        @inbounds for i in 1:length(covariates)
            U = Float64.(df[gidx, covariates[i]])
            X = hcat(X, U)
        end
        gcs[i] = GLMCopulaCSObs(Y, X, d, link)
    end
    gcm = GLMCopulaCSModel(gcs; penalized = penalized);
    return gcm
end

function CS_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    covariates::Vector{Symbol},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:NegativeBinomial, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{NBCopulaCSObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        @inbounds for i in 1:length(covariates)
            U = Float64.(df[gidx, covariates[i]])
            X = hcat(X, U)
        end
        gcs[i] = NBCopulaCSObs(Y, X, d, link)
    end
    gcm = NBCopulaCSModel(gcs; penalized = penalized);
    return gcm
end


text:
include("../src/AsymptoticBoltzmann.jl")
using Plots

function mx_scan(mxs, mv, gvxx, eps)
    rds = Array{Float64,2}(undef, length(mxs), 3)

    Threads.@threads for i = 1:length(mxs)
        mx = mxs[i]
        model = AsymptoticBoltzmann.KineticMixing(mx, mv, gvxx, eps)

        try
            rds[i, 1] = AsymptoticBoltzmann.relic_density(
                model,
                AsymptoticBoltzmann.ODE(
                    xstart = 1,
                    xend = 5e4,
                    reltol = 1e-14,
                    abstol = 2e-19,
                ),
            )
        catch
            rds[i, 1] = NaN
        end

        try
            rds[i, 2] = AsymptoticBoltzmann.relic_density(
                model,
                AsymptoticBoltzmann.MPU(),
            )
        catch
            rds[i, 2] = NaN
        end

        try
            if 2mx < mv
                α = max((mv - 2mx) / mx, 0.0)
            else
                α = 2 * (mv - mx) / mx
            end
            rds[i, 3] = AsymptoticBoltzmann.relic_density(
                model,
                AsymptoticBoltzmann.MPU(α = max(α, 0.0)),
            )
        catch
            rds[i, 3] = NaN
        end
    end
    rds
end

rds = begin
    MV = 1e3
    EPS = 1e-3
    GVXX = 1.0
    rs = 10 .^ LinRange(log10(0.1), log10(1.2), 500)
    mxs = rs .* MV

    mx_scan(mxs, MV, GVXX, EPS)
end

xfs = begin
    MV = 1e3
    EPS = 1e-3
    GVXX = 1.0
    rs = 10 .^ LinRange(log10(0.1), log10(1.2), 500)
    mxs = rs .* MV

    _xfs = zeros(eltype(mxs), length(mxs))

    for (i, mx) in enumerate(mxs)
        model = AsymptoticBoltzmann.KineticMixing(mx, MV, GVXX, EPS)
        _xfs[i] = AsymptoticBoltzmann.compute_xf(model)
    end
    _xfs
end

tcs = begin

    MV = 1e3
    EPS = 1e-3
    GVXX = 1.0
    MX = MV / 2 * 1.1
    model = AsymptoticBoltzmann.KineticMixing(MX, MV, GVXX, EPS)

    xs = 10 .^ LinRange(-1.0, 2.0, 500)
    _tcs = zeros(eltype(mxs), length(mxs))
    for (i,x) in enumerate(xs)
        _tcs[i] = AsymptoticBoltzmann.dm_thermal_cross_section(x,model)
    end
    _tcs
end

cs = begin

    MV = 1e3
    EPS = 1e-3
    GVXX = 1.0
    MX = 100.0
    model = AsymptoticBoltzmann.KineticMixing(MX, MV, GVXX, EPS)

    Qs = 10 .^ LinRange(log10(2.001*MX), log10(50*MX), 500)
    _cs = zeros(eltype(mxs), length(mxs))
    for (i,Q) in enumerate(Qs)
        _cs[i] = AsymptoticBoltzmann.dm_annihilation_cross_section(Q,model)
    end
    _cs
end

log.(rds[:, 1])

plot(rs, abs.(rds[:, 1] - rds[:, 2]) ./ rds[:, 1], yaxis = :log)
plot!(rs, abs.(rds[:, 1] - rds[:, 3]) ./ rds[:, 1], yaxis = :log)

plot(rs, log.(rds[:, 1]))

@show log.(rds[:, 1])

@show cs

max_stars_repo_name: LoganAMorrison/AsymptoticBoltzmann

max_issues_repo_name: LoganAMorrison/AsymptoticBoltzmann

max_forks_repo_name: LoganAMorrison/AsymptoticBoltzmann


text:
function build_adiabatic_frame_evolution(hp, stage; polaron=false)
    ilvl = init_lvl(hp)
    gap_lvl = ilvl == 1 ? 1 : 2
    s_axis = range(0, 1, length=1001)
    H = build_H(hp, stage, polaron=polaron)

    # calculate the position of crossing
    sm, = min_gap(H, gap_lvl, gap_lvl + 1)

    dH = build_dH(hp, stage, polaron=polaron)
    coupling = build_couplings(polaron=polaron)
    if stage == 1
        inds = sm < 0.5 ? findlast((x) -> x < sm, s_axis) : nothing
        proj = project_to_lowlevel(H, s_axis, coupling, dH, lvl=5, direction=:backward)
        adiabatic_H = build_adiabatic_H(proj)
        couplings = interp_couplings(proj, inds, sm, polaron, stage)
    elseif stage == 2
        w, v = eigen_decomp(H, 0.0, lvl=5)
        adiabatic_H = AdiabaticFrameHamiltonian((s) -> w[1:4], nothing)
        couplings = rotate_couplings(coupling, v[:, 1:4], polaron)
    elseif stage == 3
        inds = sm > 0.5 ? findlast((x) -> x < sm, s_axis) : nothing
        proj = project_to_lowlevel(H, s_axis, coupling, dH, lvl=5)
        adiabatic_H = build_adiabatic_H(proj)
        couplings = interp_couplings(proj, inds, sm, polaron, stage)
    else
        throw(ArgumentError("Stage $stage not defined."))
    end
    cb = build_diabatic_callback(sm, gap_lvl, stage)
    adiabatic_H, couplings, cb
end

function build_diabatic_callback(sm, swap_lvl, stage)
    if stage == 1
        if sm < 0.5
            op = swap_operator(swap_lvl)
            cb = InstPulseCallback([2 * τ1 * sm], (c, i) -> c .= op * c * op')
        else
            cb = nothing
        end
    elseif stage == 2
        cb = nothing
    elseif stage == 3
        if sm > 0.5
            op = swap_operator(swap_lvl)
            cb = InstPulseCallback([2 * τ1 * sm], (c, i) -> c .= op * c * op')
        else
            cb = nothing
        end
    else
        throw(ArgumentError("Stage $stage not defined."))
    end
    cb
end

function build_adiabatic_H(proj)
    D = construct_interpolations(proj.s, hcat(proj.ev...)[1:4, :], order=2)
    Dfun(s) = [D(i, s) for i in 1:4]
    AdiabaticFrameHamiltonian(Dfun, nothing)
end

function min_gap(H, i=1, j=2)
    gap = function (s)
        w, = eigen_decomp(H, s, lvl=4)
        w[j] - w[i]
    end
    res1 = optimize(gap, 0, 0.5)
    res2 = optimize(gap, 0.5, 1)
    v, ind = findmin([Optim.minimum(res1), Optim.minimum(res2)])
    Optim.minimizer([res1, res2][ind]), v
end

function swap_operator(swap_lvl)
    a1 = zeros(4)
    a2 = zeros(4)
    a1[swap_lvl] = 1
    a2[swap_lvl + 1] = 1
    op = a1 * a2' + a2 * a1'
    resi = [i for i = 1:4 if !(i in [swap_lvl, swap_lvl + 1])]
    for i in resi
        op[i, i] = 1
    end
    op
end


text:
type RandomSelectionDefinition <: SelectionDefinition; end
type RandomSelection <: Selection; end

"""
Composes a random selection operator from its definition.
"""
compose!(d::RandomSelectionDefinition) = RandomSelection()

"""
Random selection operates by selecting individuals at random with
replacement.
"""
random() = RandomSelectionDefinition()

select_ids{F}(::RandomSelection, ::FitnessScheme, ::Vector{Tuple{Int, F}}, n::Int) =
  rand(1:length(ic.fitnesses), n)


text:
import CLBLAS
import OpenCL

    const clblas = CLBLAS
    clblas.setup()

    const cl = OpenCL

    device, ctx, queue = clblas.get_next_compute_context()

    N = unsigned(5)
    data = Array(Float32, 5)
    data[1] = 1.1
    data[2] = 2.22
    data[3] = 3.333
    data[4] = 4.4444
    data[5] = 5.55555

    bufX = cl.Buffer(Float32, ctx, (:r, :copy),hostbuf=data)
    bufAsum = cl.Buffer(Float32, ctx, :w, length(bufX))
    scratchBuff = cl.Buffer(Float32, ctx, :rw, length(bufX))

    offx = unsigned(0)
    incx = cl.cl_int(1)
    ncq = cl.cl_uint(1)
    clq = queue.id
    newl = cl.cl_uint(0)

    event = cl.UserEvent(ctx, retain=true)
    ptrEvent = [event.id]

    clblas.clblasSasum( N, bufAsum.id, 0, bufX.id, 0, incx, scratchBuff.id,
                                    ncq, [clq], newl, C_NULL, ptrEvent);

    cl.api.clWaitForEvents(cl.cl_uint(1), ptrEvent)

    result = Array(Float32, length(1))
    cl.enqueue_read_buffer(queue, bufAsum, result, unsigned(0), nothing, true)

    expected = data[1:end] |> sum
    println(result)
    println("------------")
    println(expected)

    try
        @assert(isapprox(norm(expected - result[1]), zero(Float32)))
        info("success!")
    finally
        clblas.teardown()
    end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.7.0
#     language: julia
#     name: julia-1.7
# ---

# %%
struct Foo{A, B}
    a::A
    b::B
end

# %%
methods(Foo)

# %%
methods(Foo{Int, Int})

# %%
foo = Foo(1.2, 34)

# %%
mult(x::Foo{<:Number, <:Number}) = x.a * x.b

# %%
methods(mult)

# %%
mult(Foo(2, 3 + 4im))

# %%
mult(x::Foo{<:Integer, <:AbstractString}) = x.b ^ x.a

# %%
methods(mult)

# %%
mult(Foo(10, "hoge! "))

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using Conda
#===============================================================================
NOTE: pyxid is a conda package I built from a PyPI package,
following the directions at https://conda.io/docs/build_tutorials/pkgs.html.

This process requried one small modification for the package to work on windows.
I had to change the file meta.yaml generated by `conda skeleton pypi pyxid`
to contain the following.

    preserve_egg_dir: True
===============================================================================#
Conda.add_channel("haberdashPI"); Conda.add("pyxid")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayFrontendIPConfigurationPropertiesFormat <: SwaggerModel
    privateIPAddress::Any # spec type: Union{ Nothing, String } # spec name: privateIPAddress
    privateIPAllocationMethod::Any # spec type: Union{ Nothing, String } # spec name: privateIPAllocationMethod
    subnet::Any # spec type: Union{ Nothing, SubResource } # spec name: subnet
    publicIPAddress::Any # spec type: Union{ Nothing, SubResource } # spec name: publicIPAddress
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function ApplicationGatewayFrontendIPConfigurationPropertiesFormat(;privateIPAddress=nothing, privateIPAllocationMethod=nothing, subnet=nothing, publicIPAddress=nothing, provisioningState=nothing)
        o = new()
        validate_property(ApplicationGatewayFrontendIPConfigurationPropertiesFormat, Symbol("privateIPAddress"), privateIPAddress)
        setfield!(o, Symbol("privateIPAddress"), privateIPAddress)
        validate_property(ApplicationGatewayFrontendIPConfigurationPropertiesFormat, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        setfield!(o, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        validate_property(ApplicationGatewayFrontendIPConfigurationPropertiesFormat, Symbol("subnet"), subnet)
        setfield!(o, Symbol("subnet"), subnet)
        validate_property(ApplicationGatewayFrontendIPConfigurationPropertiesFormat, Symbol("publicIPAddress"), publicIPAddress)
        setfield!(o, Symbol("publicIPAddress"), publicIPAddress)
        validate_property(ApplicationGatewayFrontendIPConfigurationPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ApplicationGatewayFrontendIPConfigurationPropertiesFormat

const _property_map_ApplicationGatewayFrontendIPConfigurationPropertiesFormat = Dict{Symbol,Symbol}(Symbol("privateIPAddress")=>Symbol("privateIPAddress"), Symbol("privateIPAllocationMethod")=>Symbol("privateIPAllocationMethod"), Symbol("subnet")=>Symbol("subnet"), Symbol("publicIPAddress")=>Symbol("publicIPAddress"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ApplicationGatewayFrontendIPConfigurationPropertiesFormat = Dict{Symbol,String}(Symbol("privateIPAddress")=>"String", Symbol("privateIPAllocationMethod")=>"String", Symbol("subnet")=>"SubResource", Symbol("publicIPAddress")=>"SubResource", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ ApplicationGatewayFrontendIPConfigurationPropertiesFormat }) = collect(keys(_property_map_ApplicationGatewayFrontendIPConfigurationPropertiesFormat))
Swagger.property_type(::Type{ ApplicationGatewayFrontendIPConfigurationPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ApplicationGatewayFrontendIPConfigurationPropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationGatewayFrontendIPConfigurationPropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationGatewayFrontendIPConfigurationPropertiesFormat[property_name]

const _allowed_ApplicationGatewayFrontendIPConfigurationPropertiesFormat_privateIPAllocationMethod = ["Static", "Dynamic"]

function check_required(o::ApplicationGatewayFrontendIPConfigurationPropertiesFormat)
    true
end

function validate_property(::Type{ ApplicationGatewayFrontendIPConfigurationPropertiesFormat }, name::Symbol, val)
    if name === Symbol("privateIPAllocationMethod")
        Swagger.validate_param(name, "ApplicationGatewayFrontendIPConfigurationPropertiesFormat", :enum, val, _allowed_ApplicationGatewayFrontendIPConfigurationPropertiesFormat_privateIPAllocationMethod)
    end
end


text:
using Documenter, Distance

makedocs(sitename="Distance.jl")


text:
module PowerModelsWildfire

import InfrastructureModels
import InfrastructureModels: nw_id_default
import PowerModels
import PowerModelsRestoration
import JuMP
import MathOptInterface
import Memento

const _MOI = MathOptInterface
const _PM = PowerModels
const _PMR = PowerModelsRestoration
const _IM = InfrastructureModels

include("core/variable.jl")
include("core/constraint_template.jl")
include("core/data.jl")

include("form/dcp.jl")

include("prob/ops.jl")
include("util/risk_heuristic.jl")

include("core/export.jl")

end


text:
using EditorUtils
using Test

@testset "All Tests" begin

    @testset "capitalize strings" begin
        @test capitalize("hello") == "Hello"
        @test capitalize("HELLO") == "Hello"
        @test capitalize("heLLo") == "Hello"
        @test capitalize("HELLO_WORLD") == "Hello_world"
        @test capitalize.(["hello", "world"]) == ["Hello", "World"]
    end

    @testset "camel case to snake case" begin
        @test snake_case("hello") == "hello"
        @test snake_case("HelloWorld") == "hello_world"
        @test snake_case("helloWorld") == "hello_world"
    end

    @testset "snake case to camel case" begin
        @test snake_case("hello") == "hello"
        @test snake_case("HelloWorld") == "hello_world"
    end

end


text:
export BrickletIndustrialQuadRelayV2Monoflop
struct BrickletIndustrialQuadRelayV2Monoflop
    value::Bool
    time::Integer
    time_remaining::Integer
end

export BrickletIndustrialQuadRelayV2SPITFPErrorCount
struct BrickletIndustrialQuadRelayV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletIndustrialQuadRelayV2Identity
struct BrickletIndustrialQuadRelayV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIndustrialQuadRelayV2
"""
4 galvanically isolated solid state relays
"""
mutable struct BrickletIndustrialQuadRelayV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIndustrialQuadRelayV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_industrial_quad_relay_v2")
        deviceInternal = package.BrickletIndustrialQuadRelayV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_value
"""
    $(SIGNATURES)

Sets the value of all four relays. A value of *true* closes the
relay and a value of *false* opens the relay.

Use :func:`Set Selected Value` to only change one relay.

All running monoflop timers will be aborted if this function is called.
"""
function set_value(device::BrickletIndustrialQuadRelayV2, value)
    device.deviceInternal.set_value(value)
end

export get_value
"""
    $(SIGNATURES)

Returns the values as set by :func:`Set Value`.
"""
function get_value(device::BrickletIndustrialQuadRelayV2)
    return device.deviceInternal.get_value()
end

export set_monoflop
"""
    $(SIGNATURES)

Configures a monoflop of the specified channel.

The second parameter is the desired value of the specified
channel. A *true* means relay closed and a *false* means relay open.

The third parameter indicates the time that the channels should hold
the value.

If this function is called with the parameters (0, 1, 1500) channel 0 will
close and in 1.5s channel 0 will open again

A monoflop can be used as a fail-safe mechanism. For example: Lets assume you
have a RS485 bus and a Industrial Quad Relay Bricklet 2.0 connected to one of
the slave stacks. You can now call this function every second, with a time
parameter of two seconds and channel 0 closed. Channel 0 will be closed all the
time. If now the RS485 connection is lost, then channel 0 will be opened in at
most two seconds.
"""
function set_monoflop(device::BrickletIndustrialQuadRelayV2, channel, value, time)
    device.deviceInternal.set_monoflop(channel, value, time)
end

export get_monoflop
"""
    $(SIGNATURES)

Returns (for the given channel) the current value and the time as set by
:func:`Set Monoflop` as well as the remaining time until the value flips.

If the timer is not running currently, the remaining time will be returned
as 0.
"""
function get_monoflop(device::BrickletIndustrialQuadRelayV2, channel)
    return device.deviceInternal.get_monoflop(channel)
end

export set_selected_value
"""
    $(SIGNATURES)

Sets the output value of the specified channel without affecting the other
channels.

A running monoflop timer for the specified channel will be aborted if this
function is called.
"""
function set_selected_value(device::BrickletIndustrialQuadRelayV2, channel, value)
    device.deviceInternal.set_selected_value(channel, value)
end

export set_channel_led_config
"""
    $(SIGNATURES)

Each channel has a corresponding LED. You can turn the LED off, on or show a
heartbeat. You can also set the LED to "Channel Status". In this mode the
LED is on if the channel is high and off otherwise.
"""
function set_channel_led_config(device::BrickletIndustrialQuadRelayV2, channel, config)
    device.deviceInternal.set_channel_led_config(channel, config)
end

export get_channel_led_config
"""
    $(SIGNATURES)

Returns the channel LED configuration as set by :func:`Set Channel LED Config`
"""
function get_channel_led_config(device::BrickletIndustrialQuadRelayV2, channel)
    return device.deviceInternal.get_channel_led_config(channel)
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletIndustrialQuadRelayV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletIndustrialQuadRelayV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletIndustrialQuadRelayV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletIndustrialQuadRelayV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletIndustrialQuadRelayV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletIndustrialQuadRelayV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletIndustrialQuadRelayV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletIndustrialQuadRelayV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletIndustrialQuadRelayV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletIndustrialQuadRelayV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletIndustrialQuadRelayV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIndustrialQuadRelayV2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletIndustrialQuadRelayV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
import VoronoiCells
import Deldir
using GeometryBasics
using Plots

# Number of generators: Deldir does not go all the way
Nsmall = [1000; 2000:2000:30_000]
Nbig = 40_000:10_000:100_000
#= Nsmall = [1000; 2000:2000:10_000] =#
#= Nbig = [10_000, 20_000] =#
N = vcat(Nsmall, Nbig)

VCtime = Vector{Float64}(undef, length(N))
Dtime = similar(VCtime)
fill!(Dtime, NaN)

# Compile methods
Deldir.voronoiarea(rand(8), rand(8))
VoronoiCells.voronoicells(rand(8), rand(8), VoronoiCells.Rectangle(Point2(0, 0), Point2(1, 1)))

# Simulation window
W = [0.0 ; 10.0 ; 0.0 ; 10.0]
rect = VoronoiCells.Rectangle(Point2(0, 0), Point2(10, 10))

for (idx, n) in enumerate(N)
	println(n)

	local x = W[1] .+ W[2]*rand(n)
	local y = W[3] .+ W[4]*rand(n)

    if n <= Nsmall[end]
		Dtime[idx] = @elapsed Deldir.voronoiarea(x, y, W)
	end

    local points = [Point2(x[k], y[k]) for k in 1:n]

    VCtime[idx] = @elapsed begin
        tess = VoronoiCells.voronoicells(points, rect)
        VoronoiCells.voronoiarea(tess)
    end

    GC.gc()
end

# ------------------------------------------------------------

x = div.(N, 1000)
maxy = 5 * ceil(maximum(filter(!isnan, Dtime)) / 5)

scatter(x, Dtime, label = "Deldir")
scatter!(x, VCtime, label = "VoronoiCells")
plot!(xlabel = "number of points in 1000s", ylabel = "time in seconds", xticks = 0:20:100, yticks = 0:5:maxy)
plot!(tickfont = font(13), legendfont = font(12))

savefig("comparison.png")


text:
module Mathy

export @$, @eduction

using Transducers
using Transducers: air

const _atmath_syntax = raw"""
    @$      op dot_expr
    @$ init op dot_expr
    @$      op { dot_expr | filter_expr }
    @$ init op { dot_expr | filter_expr }
"""

"""
$_atmath_syntax

Reduce dot call expression `dot_expr` [^dot_syntax] with binary
operator `op` without allocating any intermediate arrays.  The result
of the `dot_expr` may be filtered by `filter_expr` using set-builder
like notation `{ dot_expr | filter_expr }`.  Expression `dot_expr`
itself may contain the set-builder like notation.  Optionally, the
initial value `init` for operator `op` can be specified.

For example, the expression `@\$ + { f.(a) | p(_) }` is equivalent to
the mathematical notation

```math
\\sum \\{ y \\in f.(a) \\,|\\, p(y) \\}
=
\\sum_{x \\in a \\,|\\, p(f(x))} f(x)
```

where predicate `p` is a function that takes the output of `f` and
returns a Boolean.  In general `filter_expr` is evaluated by
"substituting" each output element of `dot_expr` to `_` in
`filter_expr`.  Likewise, `@\$ * { f.(a) | p(_) }` is equivalent to
``\\prod \\{ y \\in f.(a) \\,|\\, p(y) \\}``.

!!! note

    Although set-builder like notation is used, the order of
    application of `op` is guaranteed to be left-to-right (`foldl`).
    Thus, unlike the mathematical notation, operator `op` does not
    have to be commutative or associative.

Under the hood, `@\$ init op { dot_expr | filter_expr }` is converted
to an expression that is conceptually equivalent to

```
mapfoldl(Filter(_ -> filter_expr), op, dot_expr, init=init, simd=true)
```

where `mapfoldl` and `Filter` are implemented in
[Transducers.jl](https://github.com/tkf/Transducers.jl).  However, the
output of `dot_expr` is never "materialized"; each output element is
computed on-the-fly.

[^dot_syntax]:
    Dot Syntax for Vectorizing Functions:
    <https://docs.julialang.org/en/latest/manual/functions/#man-vectorized-1>

# Examples
```jldoctest
julia> using Mathy

julia> @\$ + (1:10) .^ 2
385

julia> ans == sum((1:10) .^ 2)
true

julia> @\$ 0.0 + (1:10) .^ 2
385.0

julia> @\$ 1000 + (1:10) .^ 2
1385

julia> @\$ + { (1:10) .^ 2 | isodd(_) }
165

julia> ans == sum(x for x in (1:10) .^ 2 if isodd(x))
true

julia> @\$ + { (1:10) .^ 2 | (_ + 2) % 3 == 0 }
259

julia> ans == sum(x for x in (1:10) .^ 2 if (x + 2) % 3 == 0)
true

julia> ∑ = +
       ∏ = *;

julia> @\$ ∑ { (1:10) .^ 2 | (_ + 2) % 3 == 0 }
259

julia> @\$ ∏ { (1:10) .^ 2 | (_ + 2) % 3 == 0 }
501760000
```

Curly braces can be nested (although it is not optimized yet as of
Transducers 0.2.1):

```jldoctest; setup = :(using Mathy)
julia> @\$ + { 1:10 | isodd(_) } .^ 2
165

julia> @\$ 0 + { 2 .* ({ 1:10 | isodd(_) } .^ 2 .+ 1) .- 1 | _ % 3 == 0 }
153

julia> ans == sum(filter(x -> x % 3 == 0, 2 .* (filter(isodd, 1:10) .^ 2 .+ 1) .- 1))
true
```
"""
macro ($)(args...)
    esc(atmath_impl(args...))
end

function atmath_impl(expr::Expr)
    @assert expr.head == :call
    if length(expr.args) == 3
        op, init, body = expr.args
        return atmath_impl(init, op, body)
    elseif length(expr.args) == 2
        op, body = expr.args
        return atmath_impl(op, body)
    else
        error("Unsupported expression ", expr, "\n",
              "Supported syntax:\n",
              _atmath_syntax)
    end
end

atmath_impl(op, body) =
    atmath_impl(Transducers.MissingInit(), op, body)

function compile_body(body)
    if body.head === :braces
        @assert length(body.args) == 1
        parsed = parse_dot_expr(body.args[1])
        if parsed isa Some
            dot_expr, filter_expr = something(parsed)
            xf_expr = compile_filter(filter_expr)
        else
            dot_expr = body.args[1]
            xf_expr = Map(identity)
        end
    else
        dot_expr = body
        xf_expr = Map(identity)
    end
    return brace_to_eduction(dot_expr), xf_expr
end

function atmath_impl(init, op, body)
    dot_expr, xf_expr = compile_body(body)
    return quote
        $mapfoldl($xf_expr, $op, $air.($dot_expr); init=$init, simd=true)
    end
end

parse_dot_expr(x) = x
function parse_dot_expr(expr::Expr)
    if expr.head === :call
        if expr.args[1] == :|
            @assert length(expr.args) == 3
            _, l, r = expr.args
            return Some((l, r))
        elseif length(expr.args) === 3  # binary operator
            op, x, y = expr.args
            a = parse_dot_expr(x)
            if a isa Some
                l, r = something(a)
                return Some((l, Expr(:call, op, r, y)))
            end
            b = parse_dot_expr(y)
            if b isa Some
                l, r = something(b)
                return Some((Expr(:call, op, x, l), r))
            end
            return expr
        else
            return expr
        end
    elseif expr.head === :if
        a1 = parse_dot_expr(expr.args[1])
        if a1 isa Some
            l, r = something(a1)
            return Some((l, Expr(expr.head, r, expr.args[2:end]...)))
        else
            return expr
        end
    else
        return expr
    end
end

function compile_filter(filter_expr)
    var = gensym("x")
    subs(x) = x
    subs(x::Expr) = Expr(x.head, subs.(x.args)...)
    subs(x::Symbol) = x == :_ ? var : x
    return :($Filter($var -> $(subs(filter_expr))))
end

"""
    isdotcall(x)

# Examples
```jldoctest
julia> using Mathy: isdotcall

julia> isdotcall(:(f.(x)))
true

julia> isdotcall(:(f(x)))
false

julia> isdotcall(:(x .+ y))  # handled separately
false
```
"""
isdotcall(::Any) = false
isdotcall(ex::Expr) =
    ex.head === :. && length(ex.args) == 2 && ex.args[2] isa Expr &&
    ex.args[2].head === :tuple

"""
    isdotbracecall(x)

# Examples
```jldoctest
julia> using Mathy: isdotbracecall

julia> isdotbracecall(:(f.{x}))
true

julia> isdotbracecall(:(f.(x)))
false

julia> isdotbracecall(:(f{x}))
false
```
"""
isdotbracecall(::Any) = false
isdotbracecall(ex::Expr) =
    ex.head === :. && length(ex.args) == 2 && ex.args[2] isa Expr &&
    ex.args[2].head === :quote && length(ex.args[2].args) == 1 &&
    ex.args[2].args[1] isa Expr &&
    ex.args[2].args[1].head == :braces && length(ex.args[2].args[1].args) == 1

brace_to_eduction(x) = x
function brace_to_eduction(ex::Expr)
    if isdotcall(ex)
        args = brace_to_eduction.(ex.args[2].args)
        return Expr(:., ex.args[1], Expr(:tuple, args...))
    #=
    elseif isdotbracecall(ex)
        return Expr(:., ex.args[1],
                    Expr(:tuple, brace_to_eduction(ex.args[2].args[1])))
    =#
    elseif ex.head === :braces
        @assert length(ex.args) == 1
        return :($(@__MODULE__).@eduction $ex)
    elseif ex.head === :call
        return Expr(:call, brace_to_eduction.(ex.args)...)
    else
        return ex
    end
end

"""
    @eduction dot_expr
    @eduction { dot_expr | filter_expr }

Like [`@\$`](@ref) but without `op`.

# Examples
```jldoctest
julia> using Mathy

julia> collect(@eduction { (1:10) .^ 2 | isodd(_) })
5-element Array{Int64,1}:
  1
  9
 25
 49
 81
```
"""
macro eduction(ex)
    esc(ateduction_impl(ex))
end

function ateduction_impl(ex)
    dot_expr, xf_expr = compile_body(ex)
    return :($eduction($xf_expr, $air.($dot_expr)))
end

end # module


text:
# canonical data version: 1.2.1

using Test

include("collatz-conjecture.jl")

# canonical data
@testset "Canonical data" begin
    @test collatz_steps(1) == 0
    @test collatz_steps(16) == 4
    @test collatz_steps(12) == 9
    @test collatz_steps(1000000) == 152
    @test_throws DomainError collatz_steps(0)
    @test_throws DomainError collatz_steps(-15)
end


text:
module ReusableResourceAllocation

using JuMP
using Gurobi

using LinearAlgebra
using MathOptInterface

export reusable_resource_allocation

function reusable_resource_allocation(
		initial_supply::Array{<:Real,1},
		supply::Array{<:Real,2},
		demand::Array{<:Real,2},
		adj_matrix::BitArray{2};
		obj_dir::Symbol=:shortage,
		send_new_only::Bool=false,
		sendrecieve_switch_time::Int=0,
		min_send_amt::Real=0,
		smoothness_penalty::Real=0,
		setup_cost::Real=0,
		sent_penalty::Real=0,
		verbose::Bool=false,
)
	N, T = size(supply)
	@assert(size(initial_supply, 1) == N)
	@assert(size(demand, 1) == N)
	@assert(size(demand, 2) == T)
	@assert(size(adj_matrix, 1) == N)
	@assert(size(adj_matrix, 2) == N)
	@assert(obj_dir in [:shortage, :overflow])

	model = Model(Gurobi.Optimizer)
	if !verbose set_silent(model) end

	@variable(model, sent[1:N,1:N,1:T])
	@variable(model, obj_dummy[1:N,1:T] >= 0)

	if min_send_amt <= 0
		@constraint(model, sent .>= 0)
	else
		@constraint(model, [i=1:N,j=1:N,t=1:T], sent[i,j,t] in MOI.Semicontinuous(Float64(min_send_amt), Inf))
	end

	objective = @expression(model, sum(obj_dummy))
	if sent_penalty > 0
		add_to_expression!(objective, sent_penalty*sum(sent))
	end
	if smoothness_penalty > 0
		@variable(model, smoothness_dummy[i=1:N,j=1:N,t=1:T-1] >= 0)
		@constraint(model, [t=1:T-1],  (sent[:,:,t] - sent[:,:,t+1]) .<= smoothness_dummy[:,:,t])
		@constraint(model, [t=1:T-1], -(sent[:,:,t] - sent[:,:,t+1]) .<= smoothness_dummy[:,:,t])

		add_to_expression!(objective, smoothness_penalty * sum(smoothness_dummy))
		add_to_expression!(objective, smoothness_penalty * sum(sent[:,:,1]))
	end
	if setup_cost > 0
		@variable(model, setup_dummy[i=1:N,j=i+1:N], Bin)
		@constraint(model, [i=1:N,j=i+1:N], [1-setup_dummy[i,j], sum(sent[i,j,:])+sum(sent[j,i,:])] in MOI.SOS1([1.0, 1.0]))
		add_to_expression!(objective, setup_cost*sum(setup_dummy))
	end
	@objective(model, Min, objective)

	if send_new_only
		@constraint(model, [t=1:T],
			sum(sent[:,:,t], dims=2) .<= max.(0, supply[:,t])
		)
	else
		@constraint(model, [i=1:N,t=1:T],
			sum(sent[i,:,t]) <=
				initial_supply[i]
				+ sum(supply[i,1:t])
				- sum(sent[i,:,1:t-1])
				+ sum(sent[:,i,1:t-1])
		)
	end

	for i = 1:N
		for j = 1:N
			if ~adj_matrix[i,j]
				@constraint(model, sum(sent[i,j,:]) .== 0)
			end
		end
	end

	if sendrecieve_switch_time > 0
		@constraint(model, [i=1:N,t=1:T-1],
			[sum(sent[:,i,t]), sum(sent[i,:,t:min(t+sendrecieve_switch_time,T)])] in MOI.SOS1([1.0, 1.0])
		)
		@constraint(model, [i=1:N,t=1:T-1],
			[sum(sent[:,i,t:min(t+sendrecieve_switch_time,T)]), sum(sent[i,:,t])] in MOI.SOS1([1.0, 1.0])
		)
	end

	flip_sign = (obj_dir == :shortage) ? 1 : -1
	z1, z2 = (obj_dir == :shortage) ? (0, -1) : (-1, 0)
	@constraint(model, [i=1:N,t=1:T],
		obj_dummy[i,t] >= flip_sign * (
			demand[i,t] - (
				initial_supply[i]
				+ sum(supply[i,1:t])
				- sum(sent[i,:,1:t+z1])
				+ sum(sent[:,i,1:t+z2])
			)
		)
	)

	optimize!(model)
	return model
end

end;

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


text:
# The Dual Arrhenius and Michealis-Menten (DAMM) model, Davidson et al. 2012

"""
    DAMM(x::VecOrMat{<: Real}, p::NTuple{7, Float64})

Calculate respiration as a function of soil temperature (Tₛ) and moisture (θ).

# Examples
```julia-repl
julia> df = DAMMfdata(100) # generates a fake dataset
100×3 DataFrame
 Row │ Tₛ       θ        Rₛ
     │ Float64  Float64  Float64
─────┼─────────────────────────────
   1 │    15.5      0.3   1.72216
   2 │    22.3      0.6   1.8213
  ⋮  │    ⋮        ⋮         ⋮
  99 │     9.5      0.2   0.223677
 100 │     6.6      0.6   0.730627
julia> fp # parameters: αₛₓ, Eaₛₓ, kMₛₓ, kMₒ₂, Sxₜₒₜ, Q10kM
(1.0e9, 64.0, 3.46e-8, 0.002, 0.7, 0.02, 1.0)
julia> DAMM(hcat(df.Tₛ, df.θ), fp) # μmolCO₂ m⁻² s⁻¹
100-element Vector{Float64}:
 6.023429035220588
 0.9298933641647085
 ⋮
 0.8444248717855868
 3.805243237387702
```
"""
function DAMM(x::VecOrMat{<: Real}, p::NTuple{7, Float64})
# Independent variables
  Tₛ = x[:, 1]°C # Soil temperature
  Tₛ = Tₛ .|> K # Tₛ in Kelvin
  θ = x[:, 2]m^3*m^-3 # Soil moisture, m³ m⁻³
# Parameters
  αₛₓ = p[1]mg*cm^-3*hr^-1 # Pre-exponential factor, mgC cm⁻³ h⁻¹
  Eaₛₓ = p[2]kJ*mol^-1 # Activation energy, kJ mol⁻¹
  kMₛₓ = p[3]g*cm^-3 # Michaelis constant, gC cm⁻³
  kMₒ₂ = p[4]L*L^-1 # Michaelis constant for O₂, L L⁻¹
  porosity = p[5] # 1 - soil buld density / soil particle density
  Sxₜₒₜ = p[6]g*cm^-3 # Total soil carbon, gC cm⁻³
  Q10Km = p[7]
# DAMM model
  Tref = 293.15K
  Vmax = @. αₛₓ * exp(-Eaₛₓ/(R * Tₛ)) # Maximum potential rate of respiration
  Sₓ = @. pₛₓ * Sxₜₒₜ * Dₗᵢ * θ^3 # All soluble substrate, gC cm⁻³
  MMₛₓ = @. Sₓ / (kMₛₓ * Q10Km^(ustrip(Tₛ - Tref)/10.0) + Sₓ) # Availability of substrate factor, 0-1
  porosityₐᵢᵣ = porosity .- θ # Air filled porosity
  porosityₐᵢᵣ[porosityₐᵢᵣ .< 0.0] .= NaN # Moisture cannot be greater than porosity
  O₂ = @. Dₒₐ * O₂ₐ * (porosityₐᵢᵣ^(4/3)) # Oxygen concentration
  MMₒ₂ = @. O₂ / (kMₒ₂ + O₂) # Oxygen limitation factor, 0-1
  Rₛₘ = @. Vmax * MMₛₓ * MMₒ₂ # Respiration, mg C cm⁻³ hr⁻¹
  Rₛₘₛ = @. Rₛₘ * Eₛ # Respiration, effective depth 10 cm, mg C cm⁻² hr⁻¹
  Rₛₛₘ = Rₛₘₛ .|> molC*m^-2*s^-1
  Rₛ = Rₛₛₘ .|> μmolCO₂*m^-2*s^-1
  Rₛ = Float64.(ustrip(Rₛ)) # no units
  return Rₛ
end


text:
function input(prompt::AbstractString)
    print(prompt)
    return readline()
end

n = input("Upper bound for dimension 1: ") |>
    x -> parse(Int, x)
m = input("Upper bound for dimension 2: ") |>
    x -> parse(Int, x)

x = rand(n, m)
display(x)
x[3, 3]         # overloads `getindex` generic function
x[3, 3] = 5.0   # overloads `setindex!` generic function
x::Matrix # `Matrix{T}` is an alias for `Array{T, 2}`
x = 0; gc() # Julia has no `del` command, rebind `x` and call the garbage collector

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module TestXGBoost

using MLJ
using Test
using Pkg
Pkg.clone("https://github.com/dmlc/XGBoost.jl")
import XGBoost

function readlibsvm(fname::String, shape)
    dmx = zeros(Float32, shape)
    label = Float32[]
    fi = open(fname, "r")
    cnt = 1
    for line in eachline(fi)
        line = split(line, " ")
        push!(label, parse(Float64, line[1]))
        line = line[2:end]
        for itm in line
            itm = split(itm, ":")
            dmx[cnt, parse(Int, itm[1]) + 1] = float(parse(Int, itm[2]))
        end
        cnt += 1
    end
    close(fi)
    return (dmx, label)
end

train_X, train_Y = readlibsvm("../data/agaricus.txt.train", (6513, 126))
test_X, test_Y = readlibsvm("../data/agaricus.txt.test", (1611, 126))

dtrain = XGBoost.DMatrix("../data/agaricus.txt.train")
dtest = XGBoost.DMatrix("../data/agaricus.txt.test")

#import XGBoost_
bareboost = XGBoostRegressor(num_round=6)
bfit,= MLJ.fit(bareboost,1,train_X,train_Y)

num_round=6
bst = XGBoost.xgboost(train_X, num_round, label = train_Y)
bst_pred = XGBoost.predict(bst,test_X)
bst_DMatrix, = MLJ.fit(bareboost,1,dtrain)

bpredict_fulltree = MLJ.predict(bareboost,bfit,test_X)
bpredict_onetree = MLJ.predict(bareboost,bfit,test_X,ntree_limit=1)
bpredict_DMatrix = MLJ.predict(bareboost,bst_DMatrix,dtest)

@test  bpredict_fulltree !=bpredict_onetree
@test  bpredict_fulltree ≈ bpredict_DMatrix atol=0.000001
@test_logs (:warn,"updater has been changed to shotgun, the default option for booster=\"gblinear\"") XGBoostRegressor(num_round=1,booster="gblinear",updater="grow_colmaker")
@test_logs (:warn,"\n num_class has been changed to 2") XGBoostClassifier(eval_metric="mlogloss",objective="multi:softprob")
@test bst_pred ≈ bpredict_fulltree atol=0.000001
end
true

max_stars_repo_name: ablaom/MLJModels

max_issues_repo_name: ablaom/MLJModels

max_forks_repo_name: ablaom/MLJModels


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    CubicVariogram(range=r, sill=s, nugget=n)
    CubicVariogram(ball; sill=s, nugget=n)

A cubic variogram with range `r`, sill `s` and nugget `n`.
Optionally, use a custom metric `ball`.
"""
struct CubicVariogram{V,B} <: Variogram
  sill::V
  nugget::V
  ball::B
end

CubicVariogram(ball; sill=1.0, nugget=zero(typeof(sill))) =
  CubicVariogram(sill, nugget, ball)

CubicVariogram(; range=1.0, sill=1.0, nugget=zero(typeof(sill))) =
  CubicVariogram(sill, nugget, MetricBall(range))

function (γ::CubicVariogram)(h::T) where {T}
  r = radius(γ.ball)
  s = γ.sill
  n = γ.nugget

  # constants
  c1 = T(35) / T(4)
  c2 = T(7)  / T(2)
  c3 = T(3)  / T(4)
  s1 = 7*(h/r)^2 - c1*(h/r)^3 + c2*(h/r)^5 - c3*(h/r)^7
  s2 = T(1)

  (h < r) * (s - n) * s1 +
  (h ≥ r) * (s - n) * s2 +
  (h > 0) * n
end

isstationary(::Type{<:CubicVariogram}) = true


text:
#
# Parallel thread spliiter
#
splitter(first,nbatches,n) = first:nbatches:n

"""

```
reduce(output, output_threaded)
```

Functions to reduce the output of common options (vectors of numbers
and vectors of vectors). This function can be replacted by custom
reduction methods. It always must both receive the `output` variable
as a parameter, and return it at the end.

"""
reduce(output::Number, output_threaded::Vector{<:Number}) = sum(output_threaded)
function reduce(output::AbstractVector, output_threaded::AbstractVector{<:AbstractVector})
    @. output = output_threaded[1]
    for ibatch in 2:length(output_threaded)
        @. output += output_threaded[ibatch]
    end
    return output
end

#
# Serial version for self-pairwise computations
#
function map_pairwise_serial!(
    f::F, output, box::Box, cl::CellList{N,T};
    show_progress::Bool=false
) where {F,N,T}
    @unpack n_cells_with_real_particles = cl
    show_progress && (p = Progress(n_cells_with_real_particles, dt=1))
    ibatch = 1
    for i in 1:n_cells_with_real_particles
        cellᵢ = cl.cells[cl.cell_indices_real[i]]
        output = inner_loop!(f, box, cellᵢ, cl, output, ibatch)
        show_progress && next!(p)
    end
    return output
end

#
# Parallel version for self-pairwise computations
#
function map_pairwise_parallel!(
    f::F1, output, box::Box, cl::CellList{N,T};
    output_threaded=output_threaded,
    reduce::F2=reduce,
    show_progress::Bool=false
) where {F1,F2,N,T}
    @unpack n_cells_with_real_particles = cl
    show_progress && (p = Progress(n_cells_with_real_particles, dt=1))
    nbatches = cl.nbatches.map_computation
    @sync for ibatch in 1:nbatches
        Threads.@spawn begin
            for i in splitter(ibatch, nbatches, n_cells_with_real_particles)
                cellᵢ = cl.cells[cl.cell_indices_real[i]]
                output_threaded[ibatch] =
                    inner_loop!(f, box, cellᵢ, cl, output_threaded[ibatch], ibatch)
                show_progress && next!(p)
            end
        end
    end
    output = reduce(output, output_threaded)
    return output
end

"""

```
partition!(x::AbstractVector,by)
```

Internal function or structure - interface may change.

# Extended help

Function that reorders `x` vector by putting in the first positions the
elements with values satisfying `by(el)`. Returns the number of elements
that satisfy the condition.

"""
function partition!(by, x::AbstractVector)
    iswap = 1
    @inbounds for i in eachindex(x)
        if by(x[i])
            if iswap != i
                x[iswap], x[i] = x[i], x[iswap]
            end
            iswap += 1
        end
    end
    return iswap - 1
end

#
# Serial version for cross-interaction computations
#
function map_pairwise_serial!(
    f::F, output, box::Box,
    cl::CellListPair{N,T};
    show_progress=show_progress
) where {F,N,T}
    show_progress && (p = Progress(length(cl.ref), dt=0))
    for i in eachindex(cl.ref)
        output = inner_loop!(f, output, i, box, cl)
        show_progress && next!(p)
    end
    return output
end

#
# Parallel version for cross-interaction computations
#
function map_pairwise_parallel!(
    f::F1, output, box::Box,
    cl::CellListPair{N,T};
    output_threaded=output_threaded,
    reduce::F2=reduce,
    show_progress=show_progress
) where {F1,F2,N,T}
    show_progress && (p = Progress(length(cl.ref), dt=1))
    nbatches = cl.target.nbatches.map_computation
    @sync for ibatch in 1:nbatches
        Threads.@spawn begin
            for i in splitter(ibatch, nbatches, length(cl.ref))
                output_threaded[ibatch] = inner_loop!(f, output_threaded[ibatch], i, box, cl)
                show_progress && next!(p)
            end
        end
    end
    output = reduce(output, output_threaded)
    return output
end

#
# Inner loop for Orthorhombic cells is faster because we can guarantee that
# there are not repeated computations even if running over half of the cells.
#
function inner_loop!(
    f,box::Box{TriclinicCell},cellᵢ,
    cl::CellList{N,T},
    output,
    ibatch
) where {N,T}
    @unpack cutoff, cutoff_sq, nc = box

    for neighbor_cell in neighbor_cells(box)
        jc_cartesian = cellᵢ.cartesian_index + neighbor_cell
        jc_linear = cell_linear_index(nc,jc_cartesian)
        # if cellⱼ is empty, cycle
        if cl.cell_indices[jc_linear] == 0
            continue
        end
        cellⱼ = cl.cells[cl.cell_indices[jc_linear]]

        # same cell
        if cellⱼ.linear_index == cellᵢ.linear_index
            for i in 1:cellᵢ.n_particles
                @inbounds pᵢ = cellᵢ.particles[i]
                (!pᵢ.real) && continue
                xpᵢ = pᵢ.coordinates
                for j in 1:cellᵢ.n_particles
                    @inbounds pⱼ = cellᵢ.particles[j]
                    if pᵢ.index < pⱼ.index
                        xpⱼ = pⱼ.coordinates
                        d2 = norm_sqr(xpᵢ - xpⱼ)
                        if d2 <= cutoff_sq
                            output = f(xpᵢ, xpⱼ, pᵢ.index, pⱼ.index, d2, output)
                        end
                    end
                end
            end
        # neighbor cells
        else
            # Vector connecting cell centers
            Δc = cellⱼ.center - cellᵢ.center
            Δc_norm = norm(Δc)
            Δc = Δc / Δc_norm
            pp = project_particles!(cl.projected_particles[ibatch],cellⱼ,cellᵢ,Δc,Δc_norm,box)
            for i in 1:cellᵢ.n_particles
                @inbounds pᵢ = cellᵢ.particles[i]
                (!pᵢ.real) && continue
                xpᵢ = pᵢ.coordinates

                # Partition pp array according to the current projections. This
                # is faster than it looks, because after the first partitioning, the
                # array will be almost correct, and essentially the algorithm
                # will only run over most elements. Avoiding this partitioning or
                # trying to sort the array before always resulted to be much more
                # expensive.
                xproj = dot(xpᵢ - cellᵢ.center, Δc)
                n = partition!(el -> abs(el.xproj - xproj) <= cutoff, pp)

                for j in 1:n
                    @inbounds pⱼ = pp[j]
                    if pᵢ.index < pⱼ.index
                        xpⱼ = pⱼ.coordinates
                        d2 = norm_sqr(xpᵢ - xpⱼ)
                        if d2 <= cutoff_sq
                            output = f(xpᵢ, xpⱼ, pᵢ.index, pⱼ.index, d2, output)
                        end
                    end
                end
            end
        end
    end

    return output
end

#
# Inner loop for Orthorhombic cells is faster because we can guarantee that
# there are not repeated computations even if running over half of the cells.
#
function inner_loop!(
    f,box::Box{OrthorhombicCell},cellᵢ,
    cl::CellList{N,T},
    output,
    ibatch
) where {N,T}
    @unpack cutoff_sq = box

    # loop over list of non-repeated particles of cell ic
    for i in 1:cellᵢ.n_particles - 1
        @inbounds pᵢ = cellᵢ.particles[i]
        xpᵢ = pᵢ.coordinates
        for j in i + 1:cellᵢ.n_particles
            @inbounds pⱼ = cellᵢ.particles[j]
            xpⱼ = pⱼ.coordinates
            d2 = norm_sqr(xpᵢ - xpⱼ)
            if d2 <= cutoff_sq
                output = f(xpᵢ, xpⱼ, pᵢ.index, pⱼ.index, d2, output)
            end
        end
    end

    for jcell in neighbor_cells_forward(box)
        jc_linear = cell_linear_index(box.nc,cellᵢ.cartesian_index + jcell)
        if cl.cell_indices[jc_linear] != 0
            cellⱼ = cl.cells[cl.cell_indices[jc_linear]]
            output = cell_output!(f, box, cellᵢ, cellⱼ, cl, output, ibatch)
        end
    end

    return output
end

function cell_output!(
    f,
    box::Box{OrthorhombicCell},
    cellᵢ,
    cellⱼ,
    cl::CellList{N,T},
    output,
    ibatch
) where {N,T}
    @unpack cutoff, cutoff_sq, nc = box

    # Vector connecting cell centers
    Δc = cellⱼ.center - cellᵢ.center
    Δc_norm = norm(Δc)
    Δc = Δc / Δc_norm
    pp = project_particles!(cl.projected_particles[ibatch],cellⱼ,cellᵢ,Δc,Δc_norm,box)
    if length(pp) == 0
        return output
    end

    # Loop over particles of cell icell
    for i in 1:cellᵢ.n_particles
        @inbounds pᵢ = cellᵢ.particles[i]
        xpᵢ = pᵢ.coordinates
        xproj = dot(xpᵢ - cellᵢ.center, Δc)

        # Partition pp array according to the current projections
        n = partition!(el -> abs(el.xproj - xproj) <= cutoff, pp)

        # Compute the interactions
        for j in 1:n
            @inbounds pⱼ = pp[j]
            xpⱼ = pⱼ.coordinates
            d2 = norm_sqr(xpᵢ - xpⱼ)
            if d2 <= cutoff_sq
                output = f(xpᵢ, xpⱼ, pᵢ.index, pⱼ.index, d2, output)
            end
        end
    end

    return output
end

"""

```
project_particles!(projected_particles,cellⱼ,cellᵢ,Δc,Δc_norm,box)
```

Internal function or structure - interface may change.

# Extended help

Projects all particles of the cell `cellⱼ` into unnitary vector `Δc` with direction
connecting the centers of `cellⱼ` and `cellᵢ`. Modifies `projected_particles`, and
returns a view of `projected particles, where only the particles for which
the projection on the direction of the cell centers still allows the particle
to be within the cutoff distance of any point of the other cell.

"""
function project_particles!(
    projected_particles,cellⱼ,cellᵢ,
    Δc,Δc_norm,box::Box{UnitCellType,N}
) where {UnitCellType,N}
    if box.lcell == 1
        margin = box.cutoff + Δc_norm/2 # half of the distance between centers
    else
        margin = box.cutoff*(1 + sqrt(N)/2) # half of the diagonal of the cutoff-box
    end
    iproj = 0
    @inbounds for j in 1:cellⱼ.n_particles
        pⱼ = cellⱼ.particles[j]
        xproj = dot(pⱼ.coordinates - cellᵢ.center, Δc)
        if abs(xproj) <= margin
            iproj += 1
            projected_particles[iproj] = ProjectedParticle(pⱼ.index, xproj, pⱼ.coordinates)
        end
    end
    pp = @view(projected_particles[1:iproj])
    return pp
end

#
# Inner loop of cross-interaction computations. If the two sets
# are large, this might(?) be improved by computing two separate
# cell lists and using projection and partitioning.
#
function inner_loop!(
    f,output,i,box,
    cl::CellListPair{N,T}
) where {N,T}
    @unpack nc, cutoff_sq = box
    xpᵢ = wrap_to_first(cl.ref[i], box)
    ic = particle_cell(xpᵢ, box)
    for neighbor_cell in neighbor_cells(box)
        jc_cartesian = neighbor_cell + ic
        jc_linear = cell_linear_index(nc,jc_cartesian)
        # If cellⱼ is empty, cycle
        if cl.target.cell_indices[jc_linear] == 0
            continue
        end
        cellⱼ = cl.target.cells[cl.target.cell_indices[jc_linear]]
        # loop over particles of cellⱼ
        for j in 1:cellⱼ.n_particles
            @inbounds pⱼ = cellⱼ.particles[j]
            xpⱼ = pⱼ.coordinates
            d2 = norm_sqr(xpᵢ - xpⱼ)
            if d2 <= cutoff_sq
                output = f(xpᵢ, xpⱼ, i, pⱼ.index, d2, output)
            end
        end
    end
    return output
end


text:
using Observers
using Plots

export AQObserver

struct AQObserver <: Observer
    snode::Int
    actions::Vector{Any}
    ns::Vector{Any}
    Xs::Vector{Any}
    ys::Vector{Any}
end
AQObserver(snode::Int) = AQObserver(snode, [], [], [], [])

function Observers.notify_observer!(o::AQObserver, b::ModularBandit; kwargs...)
    snode = kwargs_get(kwargs, :snode)
    o.snode == snode || return  #only monitor snode

    p = kwargs_get(kwargs, :planner)
    tree = get(p.tree)

    sanode = kwargs_get(kwargs, :sanode)
    push!(o.actions, tree.a_labels[sanode])

    children = tree.children[snode]
    actions = [tree.a_labels[c] for c in children]
    qs = [tree.q[c] for c in children]
    ns = [tree.n[c] for c in children]

    push!(o.Xs, actions)
    push!(o.ys, qs)
    push!(o.ns, ns)
end

@recipe function plot(o::AQObserver, i, ylim=nothing)
    @series begin
        seriestype := :scatter
        xlabel := "action"
        ylabel := "q"
        title := "n=$i"
        if ylim != nothing
            ylim := ylim
        end
        x = o.Xs[i]
        y = o.ys[i]
        x, y
    end
end

Base.length(o::AQObserver) = length(o.Xs)

function Plots.animate(o::AQObserver, fname="./aqobserver.gif"; fps=2, ylim=nothing)
    anim = @animate for i=1:length(o)
        Plots.plot(o, i, ylim=ylim)
    end
    gif(anim, fname; fps=fps)
end


text:
const GRPC_STATIC_HEADERS = Ref{Ptr{Nothing}}(C_NULL)

const StatusCode = (
    OK                  = (code=0,  message="Success"),
    CANCELLED           = (code=1,  message="The operation was cancelled"),
    UNKNOWN             = (code=2,  message="Unknown error"),
    INVALID_ARGUMENT    = (code=3,  message="Client specified an invalid argument"),
    DEADLINE_EXCEEDED   = (code=4,  message="Deadline expired before the operation could complete"),
    NOT_FOUND           = (code=5,  message="Requested entity was not found"),
    ALREADY_EXISTS      = (code=6,  message="Entity already exists"),
    PERMISSION_DENIED   = (code=7,  message="No permission to execute the specified operation"),
    RESOURCE_EXHAUSTED  = (code=8,  message="Resource exhausted"),
    FAILED_PRECONDITION = (code=9,  message="Operation was rejected because the system is not in a state required for the operation's execution"),
    ABORTED             = (code=10, message="Operation was aborted"),
    OUT_OF_RANGE        = (code=11, message="Operation was attempted past the valid range"),
    UNIMPLEMENTED       = (code=12, message="Operation is not implemented or is not supported/enabled in this service"),
    INTERNAL            = (code=13, message="Internal error"),
    UNAVAILABLE         = (code=14, message="The service is currently unavailable"),
    DATA_LOSS           = (code=15, message="Unrecoverable data loss or corruption"),
    UNAUTHENTICATED     = (code=16, message="The request does not have valid authentication credentials for the operation")
)

grpc_status_info(code) = StatusCode[code+1]
grpc_status_message(code) = (grpc_status_info(code)).message
grpc_status_code_str(code) = string(propertynames(StatusCode)[code+1])

#=
const SEND_BUFFER_SZ = 1024 * 1024
function buffer_send_data(input::Channel{T}) where T <: ProtoType
    data = nothing
    if isready(input)
        iob = IOBuffer()
        while isready(input) && (iob.size < SEND_BUFFER_SZ)
            write(iob, to_delimited_message_bytes(take!(input)))
            yield()
        end
        data = take!(iob)
    elseif isopen(input)
        data = UInt8[]
    end
    data
end
=#

function send_data(easy::Curl.Easy, input::Channel{T}, max_send_message_length::Int) where T <: ProtoType
    while true
        yield()
        data = isready(input) ? to_delimited_message_bytes(take!(input), max_send_message_length) : isopen(input) ? UInt8[] : nothing
        easy.input === nothing && break
        easy.input = data
        Curl.curl_easy_pause(easy.handle, Curl.CURLPAUSE_CONT)
        wait(easy.ready)
        easy.input === nothing && break
        easy.ready = Threads.Event()
    end
end

function grpc_timeout_header_val(timeout::Real)
    if round(Int, timeout) == timeout
        timeout_secs = round(Int64, timeout)
        return "$(timeout_secs)S"
    end
    timeout *= 1000
    if round(Int, timeout) == timeout
        timeout_millisecs = round(Int64, timeout)
        return "$(timeout_millisecs)m"
    end
    timeout *= 1000
    if round(Int, timeout) == timeout
        timeout_microsecs = round(Int64, timeout)
        return "$(timeout_microsecs)u"
    end
    timeout *= 1000
    timeout_nanosecs = round(Int64, timeout)
    return "$(timeout_nanosecs)n"
end

function grpc_headers(; timeout::Real=Inf)
    headers = C_NULL
    headers = LibCURL.curl_slist_append(headers, "User-Agent: $(Curl.USER_AGENT)")
    headers = LibCURL.curl_slist_append(headers, "Content-Type: application/grpc+proto")
    headers = LibCURL.curl_slist_append(headers, "Content-Length:")
    headers = LibCURL.curl_slist_append(headers, "te: trailers")
    if timeout !== Inf
        headers = LibCURL.curl_slist_append(headers, "grpc-timeout: $(grpc_timeout_header_val(timeout))")
    end
    headers
end

function grpc_request_header(request_timeout::Real)
    if request_timeout == Inf
        GRPC_STATIC_HEADERS[]
    else
        grpc_headers(; timeout=request_timeout)
    end
end

function easy_handle(maxage::Clong, keepalive::Clong, negotiation::Symbol, revocation::Bool, request_timeout::Real)
    easy = Curl.Easy()
    http_version = (negotiation === :http2) ? CURL_HTTP_VERSION_2_0 :
                   (negotiation === :http2_tls) ? CURL_HTTP_VERSION_2TLS :
                   (negotiation === :http2_prior_knowledge) ? CURL_HTTP_VERSION_2_PRIOR_KNOWLEDGE :
                   throw(ArgumentError("unsupported HTTP2 negotiation mode $negotiation"))
    Curl.setopt(easy, CURLOPT_HTTP_VERSION, http_version)
    Curl.setopt(easy, CURLOPT_PIPEWAIT, Clong(1))
    Curl.setopt(easy, CURLOPT_POST, Clong(1))
    Curl.setopt(easy, CURLOPT_HTTPHEADER, grpc_request_header(request_timeout))
    if !revocation
        Curl.setopt(easy, CURLOPT_SSL_OPTIONS, CURLSSLOPT_NO_REVOKE)
    end
    if maxage > 0
        Curl.setopt(easy, CURLOPT_MAXAGE_CONN, maxage)
    end
    if keepalive > 0
        Curl.setopt(easy, CURLOPT_TCP_KEEPALIVE, Clong(1))
        Curl.setopt(easy, CURLOPT_TCP_KEEPINTVL, keepalive);
        Curl.setopt(easy, CURLOPT_TCP_KEEPIDLE, keepalive);
    end
    easy
end

function recv_data(easy::Curl.Easy, output::Channel{T}, max_recv_message_length::Int) where T <: ProtoType
    iob = PipeBuffer()
    waiting_for_header = true
    msgsize = 0
    compressed = UInt8(0)
    datalen = UInt32(0)
    need_more = true
    for buf in easy.output
        write(iob, buf)
        need_more = false
        while !need_more
            if waiting_for_header
                if bytesavailable(iob) >= 5
                    compressed = read(iob, UInt8)       # compression
                    datalen = ntoh(read(iob, UInt32))   # message length

                    if datalen > max_recv_message_length
                        throw(gRPCMessageTooLargeException(max_recv_message_length, datalen))
                    end

                    waiting_for_header = false
                else
                    need_more = true
                end
            end

            if !waiting_for_header
                if bytesavailable(iob) >= datalen
                    msgbytes = IOBuffer(view(iob.data, iob.ptr:(iob.ptr+datalen-1)))
                    put!(output, readproto(msgbytes, T()))  # decode message bytes
                    iob.ptr += datalen
                    waiting_for_header = true
                else
                    need_more = true
                end
            end
        end
    end
    close(output)
end

function set_connect_timeout(easy::Curl.Easy, timeout::Real)
    timeout >= 0 ||
        throw(ArgumentError("timeout must be positive, got $timeout"))
    if timeout ≤ typemax(Clong) ÷ 1000
        timeout_ms = round(Clong, timeout * 1000)
        Curl.setopt(easy, CURLOPT_CONNECTTIMEOUT_MS, timeout_ms)
    else
        timeout = timeout ≤ typemax(Clong) ? round(Clong, timeout) : Clong(0)
        Curl.setopt(easy, CURLOPT_CONNECTTIMEOUT, timeout)
    end
end

function grpc_request(downloader::Downloader, url::String, input::Channel{T1}, output::Channel{T2};
        maxage::Clong = typemax(Clong),
        keepalive::Clong = 60,
        negotiation::Symbol = :http2_prior_knowledge,
        revocation::Bool = true,
        request_timeout::Real = Inf,
        connect_timeout::Real = 0,
        max_recv_message_length::Int = DEFAULT_MAX_RECV_MESSAGE_LENGTH,
        max_send_message_length::Int = DEFAULT_MAX_SEND_MESSAGE_LENGTH,
        verbose::Bool = false)::gRPCStatus where {T1 <: ProtoType, T2 <: ProtoType}
    Curl.with_handle(easy_handle(maxage, keepalive, negotiation, revocation, request_timeout)) do easy
        # setup the request
        Curl.set_url(easy, url)
        Curl.set_timeout(easy, request_timeout)
        set_connect_timeout(easy, connect_timeout)
        Curl.set_verbose(easy, verbose)
        Curl.add_upload_callbacks(easy)
        Downloads.set_ca_roots(downloader, easy)

        # do the request
        Curl.add_handle(downloader.multi, easy)

        function cleanup()
            Curl.remove_handle(downloader.multi, easy)
            # though remove_handle sets easy.handle to C_NULL, it does not close output and progress channels
            # we need to close them here to unblock anything waiting on them
            close(easy.output)
            close(easy.progress)
            close(output)
            close(input)
            nothing
        end

        # do send recv data
        if VERSION < v"1.5"
            cleaned_up = false
            exception = nothing
            cleanup_once = (ex)->begin
                if !cleaned_up
                    cleaned_up = true
                    exception = ex
                    cleanup()
                end
            end

            @sync begin
                @async try
                    recv_data(easy, output, max_recv_message_length)
                catch ex
                    cleanup_once(ex)
                end
                @async try
                    send_data(easy, input, max_send_message_length)
                catch ex
                    cleanup_once(ex)
                end
            end

            if exception !== nothing
                throw(exception)
            end
        else
            try
                Base.Experimental.@sync begin
                    @async recv_data(easy, output, max_recv_message_length)
                    @async send_data(easy, input, max_send_message_length)
                end
            finally # ensure handle is removed
                cleanup()
            end
        end

        @debug("response headers", easy.res_hdrs)

        # parse the grpc headers
        grpc_status = StatusCode.OK.code
        grpc_message = ""
        for hdr in easy.res_hdrs
            if startswith(hdr, "grpc-status")
                grpc_status = parse(Int, strip(last(split(hdr, ':'; limit=2))))
            elseif startswith(hdr, "grpc-message")
                grpc_message = string(strip(last(split(hdr, ':'; limit=2))))
            end
        end
        if (easy.code == CURLE_OPERATION_TIMEDOUT) && (grpc_status == StatusCode.OK.code)
            grpc_status = StatusCode.DEADLINE_EXCEEDED.code
        end
        if (grpc_status != StatusCode.OK.code) && isempty(grpc_message)
            grpc_message = grpc_status_message(grpc_status)
        end

        if ((easy.code == CURLE_OK) && (grpc_status == StatusCode.OK.code))
            gRPCStatus(true, grpc_status, "")
        else
            gRPCStatus(false, grpc_status, isempty(grpc_message) ? Curl.get_curl_errstr(easy) : grpc_message)
        end
    end
end


text:
module SparsityDetection

using SpecialFunctions
using Cassette, LinearAlgebra, SparseArrays
using Cassette: tag, untag, Tagged, metadata, hasmetadata, istagged, canrecurse
using Cassette: tagged_new_tuple, ContextTagged, BindingMeta, DisableHooks, nametype
using Core: SSAValue

export Sparsity, jacobian_sparsity, hessian_sparsity, hsparsity, sparsity!

include("util.jl")
include("controlflow.jl")
include("propagate_tags.jl")
include("linearity.jl")
include("jacobian.jl")
include("hessian.jl")
include("blas.jl")

sparsity!(args...; kwargs...) = jacobian_sparsity(args...; kwargs...)
hsparsity(args...; kwargs...) = hessian_sparsity(args...; kwargs...)

end