Skip to content

Instantly share code, notes, and snippets.

@slwu89

slwu89/ode-spline-help.jl

Last active December 1, 2025 00:12
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save slwu89/e6d51cef69a694ec78969e202fecf9f3 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save slwu89/e6d51cef69a694ec78969e202fecf9f3 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
ode-spline-help.jl
import DifferentialEquations as DE
using DataFrames
using Distributions
using LinearAlgebra, Statistics
using SciMLSensitivity
using Optimization, OptimizationLBFGSB
using Mooncake
# ------------------------------------------------------------
# make our constrained basis functions
k = 10
xk = [0.0, 0.1111111111111111, 0.22222222222222224, 0.3333333333333333, 0.4444444444444444, 0.5555555555555556, 0.6666666666666666, 0.7777777777777778, 0.8888888888888888, 1.0]
D = [9.0 -18.0 8.999999999999998 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 8.999999999999998 -18.0 9.000000000000004 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 9.000000000000004 -18.000000000000004 9.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 9.0 -17.999999999999996 8.999999999999996 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 8.999999999999996 -18.0 9.000000000000005 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 9.000000000000005 -18.0 8.999999999999996 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 8.999999999999996 -18.0 9.000000000000005 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 9.000000000000005 -18.0 8.999999999999996]
B = [0.07407407407407408 0.01851851851851852 0.0 0.0 0.0 0.0 0.0 0.0; 0.01851851851851852 0.07407407407407407 0.018518518518518514 0.0 0.0 0.0 0.0 0.0; 0.0 0.018518518518518514 0.07407407407407406 0.018518518518518517 0.0 0.0 0.0 0.0; 0.0 0.0 0.018518518518518517 0.07407407407407408 0.018518518518518528 0.0 0.0 0.0; 0.0 0.0 0.0 0.018518518518518528 0.07407407407407407 0.018518518518518507 0.0 0.0; 0.0 0.0 0.0 0.0 0.018518518518518507 0.07407407407407407 0.018518518518518528 0.0; 0.0 0.0 0.0 0.0 0.0 0.018518518518518528 0.07407407407407407 0.018518518518518507; 0.0 0.0 0.0 0.0 0.0 0.0 0.018518518518518507 0.07407407407407407]
F = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 130.223307436182 -295.3398446170921 209.35937846836848 -56.09766925638182 15.031298557158713 -4.027524972253052 1.0788013318534968 -0.2876803551609323 0.07192008879023315 -0.011986681465038842; -34.89322974472808 209.35937846836848 -351.43751387347396 224.39067702552728 -60.125194228634854 16.110099889012208 -4.315205327413987 1.1507214206437293 -0.2876803551609326 0.047946725860155366; 9.349611542730301 -56.0976692563818 224.3906770255273 -355.4650388457271 225.4694783573807 -60.412874583795784 16.182019977802454 -4.315205327413985 1.078801331853497 -0.1798002219755826; -2.505216426193119 15.031298557158712 -60.12519422863487 225.4694783573807 -355.75271920088784 225.54139844617086 -60.4128745837958 16.1100998890122 -4.027524972253054 0.6712541620421748; 0.6712541620421758 -4.027524972253055 16.110099889012222 -60.41287458379582 225.5413984461709 -355.75271920088795 225.46947835738075 -60.12519422863483 15.031298557158719 -2.505216426193117; -0.1798002219755827 1.0788013318534961 -4.315205327413985 16.18201997780244 -60.41287458379575 225.46947835738075 -355.465038845727 224.3906770255272 -56.097669256381835 9.349611542730294; 0.0479467258601554 -0.2876803551609324 1.15072142064373 -4.315205327413986 16.110099889012204 -60.12519422863489 224.3906770255272 -351.43751387347396 209.35937846836853 -34.893229744728046; -0.011986681465038843 0.07192008879023305 -0.2876803551609323 1.0788013318534957 -4.027524972253048 15.031298557158712 -56.09766925638176 209.3593784683685 -295.3398446170921 130.22330743618195; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
X = [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.8860799076581576 0.14424955405105438 -0.03845621620421753 0.010304310765815764 -0.002761026859045505 0.0007397966703662596 -0.00019815982241953395 5.2842619311875695e-5 -1.3210654827968936e-5 2.2017758046614863e-6; 0.7733318250832408 0.285841049500555 -0.07502819800221976 0.020103742508324092 -0.005386772031076582 0.001443345615982242 -0.0003866104328523865 0.00010309611542730301 -2.5774028856825777e-5 4.295671476137624e-6; 0.6629277620421754 0.4221164277469478 -0.10783171098779133 0.028893416204217544 -0.007741953829078801 0.002074399112097669 -0.0005556426193118759 0.00014817136514983352 -3.7042841287458414e-5 6.17380688124306e-6; 0.5560397283018867 0.5504176301886793 -0.13498252075471698 0.03616845283018869 -0.009691290566037737 0.002596709433962264 -0.0006955471698113212 0.0001854792452830189 -4.636981132075477e-5 7.72830188679245e-6; 0.45383973362930075 0.6680865982241954 -0.15459639289678132 0.041423973362930085 -0.011099500554938956 0.0029740288568257485 -0.0007966148723640403 0.00021243063263041065 -5.310765815760271e-5 8.851276359600439e-6; 0.35749978779134295 0.7724652732519423 -0.16478909300776912 0.0441550987791343 -0.011831302108768033 0.003170109655937846 -0.0008491365149833522 0.00022643640399556045 -5.6609100998890166e-5 9.434850166481681e-6; 0.26819190055493886 0.8608955966703664 -0.16367638668146506 0.04385695005549392 -0.011751413540510547 0.0031487041065482797 -0.0008434028856825755 0.00022490743618202002 -5.622685904550506e-5 9.37114317425083e-6; 0.18708808168701443 0.9307195098779135 -0.14937403951165368 0.040024648168701445 -0.01072455316315205 0.002873564483906769 -0.000769704772475028 0.0002052546059933407 -5.131365149833522e-5 8.552275249722526e-6; 0.115360340954495 0.97927895427303 -0.11999781709211983 0.032153314095449505 -0.008615439289678133 0.00230844306326304 -0.0006183329633740291 0.00016488879023307431 -4.122219755826862e-5 6.870366259711427e-6; 0.054180688124306285 1.003915871254162 -0.07366348501664807 0.019738068812430617 -0.005288790233074356 0.0014170921198668127 -0.00037957824639289653 0.00010122086570477235 -2.5305216426193115e-5 4.2175360710321796e-6; 0.0047211329633739755 1.001972202219756 -0.008486808879023297 0.0022740332963374008 -0.0006093243063263034 0.00016326392896781332 -4.373140954495003e-5 1.1661709211986668e-5 -2.9154273029966696e-6 4.859045504994443e-7; -0.03215745467998522 0.9716327280799113 0.07660508768035516 -0.02038907880133187 0.005463227524972256 -0.0014638312985571595 0.00039209766925638224 -0.00010455937846836854 2.6139844617092156e-5 -4.356640769515354e-6; -0.057335251572327085 0.9157985094339622 0.17814496226415102 -0.04641735849056608 0.012437471698113219 -0.0033325283018867946 0.0008926415094339634 -0.00023803773584905682 5.950943396226426e-5 -9.918238993710697e-6; -0.07229769885312617 0.8390101931187568 0.29108722752497246 -0.07328710321864604 0.019637185349611563 -0.0052616381798002276 0.0014093673695893473 -0.0003758312985571591 9.395782463928986e-5 -1.5659637439881624e-5; -0.07853084535701076 0.7458100721420644 0.4103847114317425 -0.09847391786903448 0.026385960044395126 -0.007069922308546062 0.0018937291897891252 -0.000504994450610433 0.00012624861265260836 -2.1041435442101366e-5; -0.07752073991860897 0.6407404395116538 0.5309902419533852 -0.11945340732519433 0.03200738734739181 -0.008576142064372922 0.0022971809100998915 -0.0006125815760266374 0.0001531453940066595 -2.552423233444321e-5; -0.07075343137254905 0.528343588235294 0.6478566470588236 -0.13370117647058832 0.03582505882352942 -0.009599058823529414 0.0025711764705882378 -0.0006856470588235297 0.00017141176470588258 -2.856862745098039e-5; -0.05971496855345916 0.4131618113207549 0.7559367547169811 -0.13869283018867934 0.037162566037735864 -0.009957433962264154 0.0026671698113207576 -0.0007112452830188682 0.00017781132075471722 -2.963522012578616e-5; -0.04589140029596747 0.29973740177580477 0.8501833928967814 -0.13190397336293017 0.03534350055493897 -0.009470028856825751 0.002536614872364042 -0.0006764306326304108 0.0001691076581576029 -2.8184609692933774e-5; -0.030768775434702193 0.19261265260821311 0.9255493895671477 -0.11081021087680365 0.029691453940066607 -0.00795560488346282 0.0021309655937846858 -0.0005682574916759159 0.0001420643729189791 -2.3677395486496488e-5; -0.015833142804291555 0.09632985682574932 0.9769875726970033 -0.07288714761376267 0.01953001775804666 -0.005232923418423985 0.0014016759156492825 -0.0003737802441731418 9.344506104328554e-5 -1.5574176840547567e-5; -0.002570551239363695 0.015431307436182169 0.9994507702552718 -0.015610388457269751 0.004182783573806894 -0.001120745837957828 0.0003001997780244184 -8.005327413984486e-5 2.0013318534961235e-5 -3.3355530891602013e-6; 0.0076807081761006216 -0.04608424905660372 0.9886939962264151 0.06293726415094339 -0.016772052830188663 0.004493947169811318 -0.0012037358490566037 0.00032099622641509423 -8.024905660377356e-5 1.3374842767295579e-5; 0.014755790751017375 -0.08853474450610424 0.9468429780244174 0.15985083240843484 -0.04173430765815752 0.011182398224195318 -0.002995285238623748 0.0007987427302996656 -0.0001996856825749165 3.3280947095819366e-5; 0.019049435812060662 -0.11429661487236395 0.8790614594894561 0.26992577691453934 -0.0681395671476137 0.018257491675915633 -0.004890399556048835 0.0013041065482796885 -0.0003260266370699223 5.433777284498697e-5; 0.0209598289308176 -0.1257589735849056 0.7905318943396225 0.38794339622641516 -0.09341747924528297 0.02503052075471697 -0.006704603773584909 0.0017878943396226414 -0.00044697358490566045 7.449559748427665e-5; 0.020885155678875313 -0.1253109340732519 0.6864367362930075 0.508684988901221 -0.11499769189789122 0.030812778690344055 -0.008253422863485021 0.0022009127635960043 -0.0005502281908990014 9.17046984831668e-5; 0.01922360162782092 -0.1153416097669255 0.5719584390677022 0.6269318534961157 -0.13030985305216422 0.03491555871254161 -0.009352381798002223 0.0024939684794672585 -0.000623492119866815 0.00010391535331113568; 0.01637335234924158 -0.09824011409544949 0.45227945638179806 0.7374652885682573 -0.13678361065482791 0.036650154051054366 -0.00981700554938957 0.0026178681465038843 -0.0006544670366259714 0.00010907783943766173; 0.012732593414724372 -0.07639556048834623 0.33258224195338504 0.8350665926748058 -0.13184861265260817 0.035327857935627066 -0.009462819089900114 0.0025234184239733624 -0.0006308546059933409 0.00010514243433222335; 0.008699510395856447 -0.052197062375138686 0.2180492495005548 0.914517064372919 -0.11293450699223075 0.030259963596004413 -0.008105347391786902 0.0021614259711431726 -0.0005403564927857935 9.005941546429879e-5; 0.004672288864224926 -0.028033733185349557 0.11386293274139825 0.9705980022197559 -0.07747094162042158 0.020757764261931138 -0.005560115427302989 0.001482697447280796 -0.0003706743618201992 6.177906030336646e-5; 0.0010491143914169337 -0.006294686348501601 0.025205745394006412 0.9980907047724751 -0.02288756448390652 0.0061325531631519865 -0.0016426481687014261 0.0004380395116537134 -0.00010950987791342841 1.825164631890471e-5; -0.0017967597484276848 0.010780558490566108 -0.04312223396226445 0.9922923773584905 0.053000724528302214 -0.014143275471698195 0.0037883773584905882 -0.001010233962264156 0.00025255849056603935 -4.2093081761006506e-5; -0.00378191823899371 0.022691509433962254 -0.09076603773584907 0.954497641509434 0.14790047169811316 -0.038724528301886774 0.010372641509433959 -0.0027660377358490534 0.0006915094339622641 -0.0001152515723270439; -0.0050099379948205704 0.030059627968923423 -0.12023851187569373 0.8898704195338513 0.2566288337402886 -0.06505775449500553 0.017426184239733633 -0.004646982463928965 0.0011617456159822423 -0.00019362426933037353; -0.005587512467628562 0.03352507480577137 -0.13410029922308553 0.8036391220865704 0.3739548108768037 -0.09056936559378466 0.024259651498335193 -0.006469240399556046 0.0016173100998890131 -0.0002695516833148352; -0.0056213351091379945 0.033728010654827965 -0.13491204261931192 0.7010321598224194 0.494647403329634 -0.11268577314095446 0.03018368923418424 -0.00804898379578246 0.0020122459489456167 -0.0003353743248242691; -0.005218099371069182 0.03130859622641509 -0.1252343849056604 0.5872779433962261 0.6134756113207551 -0.12883338867924524 0.03450894339622642 -0.00920238490566037 0.0023005962264150952 -0.00038343270440251544; -0.004484498705142432 0.026906992230854583 -0.10762796892341839 0.46760488346281864 0.7252084350721425 -0.1364386237513873 0.03654605993340733 -0.009745615982241945 0.002436403995560489 -0.00040606733259341437; -0.0035272265630780612 0.021163359378468367 -0.0846534375138735 0.3472413906770256 0.8246148748057713 -0.13292788990011092 0.03560568479467259 -0.00949484927857935 0.0023737123196448397 -0.00039561871994080617; -0.0024529763965963737 0.01471785837957824 -0.05887143351831299 0.23141587569367364 0.9064639307436182 -0.11572759866814644 0.030998463928967813 -0.008266257047724743 0.002066564261931188 -0.00034442737698853085; -0.0013684416574176821 0.008210649944506094 -0.032842599778024395 0.12535674916759143 0.9655246031076582 -0.0822641615982241 0.022035043285238608 -0.005876011542730291 0.0014690028856825742 -0.00024483381428042875; -0.00038031579726229887 0.002281894783573793 -0.00912757913429518 0.0342924217536069 0.9965658921198668 -0.029963990233074174 0.008026068812430587 -0.0021402850166481547 0.0005350712541620392 -8.917854236033975e-5; 0.0004085738068812458 -0.0024514428412874744 0.009805771365149901 -0.036771642619312116 0.9946557991120977 0.04352344617092147 -0.011624583795782547 0.0030998890122086783 -0.00077497225305217 0.00012916204217536156; 0.0009648035516093255 -0.0057888213096559515 0.023155285238623816 -0.08683231964483928 0.9602299933407323 0.13654434628190942 -0.03583937846836862 0.009557167591564962 -0.0023892918978912414 0.0003982153163152065; 0.0013149822419533858 -0.007889893451720316 0.03155957380688127 -0.11834840177580473 0.8983670332963376 0.24398126859045494 -0.062093107658157634 0.01655816204217536 -0.004139540510543843 0.0006899234184239732; 0.0014876981132075483 -0.008926188679245288 0.03570475471698117 -0.13389283018867934 0.814298566037736 0.3606025660377358 -0.08781283018867933 0.023416754716981143 -0.00585418867924529 0.0009756981132075473; 0.0015115394006659274 -0.009069236403995564 0.03627694561598228 -0.13603854605993349 0.7132562386237515 0.481176591564928 -0.11042560488346295 0.029446827968923436 -0.007361706992230863 0.0012269511653718096; 0.0014150943396226427 -0.008490566037735854 0.03396226415094343 -0.12735849056603785 0.6004716981132077 0.6004716981132078 -0.12735849056603787 0.03396226415094342 -0.008490566037735861 0.001415094339622642; 0.00122695116537181 -0.007361706992230859 0.029446827968923447 -0.1104256048834629 0.48117659156492787 0.7132562386237516 -0.13603854605993354 0.03627694561598226 -0.009069236403995571 0.0015115394006659268; 0.0009756981132075477 -0.005854188679245286 0.023416754716981146 -0.08781283018867928 0.36060256603773577 0.8142985660377361 -0.13389283018867937 0.035704754716981146 -0.008926188679245294 0.0014876981132075473; 0.0006899234184239734 -0.004139540510543839 0.01655816204217536 -0.06209310765815759 0.24398126859045488 0.8983670332963377 -0.11834840177580477 0.03155957380688126 -0.00788989345172032 0.0013149822419533854; 0.00039821531631520503 -0.0023892918978912297 0.009557167591564924 -0.03583937846836845 0.1365443462819088 0.9602299933407328 -0.08683231964483909 0.023155285238623743 -0.00578882130965594 0.0009648035516093223; 0.00012916204217536026 -0.0007749722530521614 0.0030998890122086475 -0.011624583795782425 0.04352344617092102 0.9946557991120977 -0.03677164261931176 0.009805771365149799 -0.0024514428412874515 0.0004085738068812415; -8.917854236034078e-5 0.0005350712541620447 -0.002140285016648179 0.008026068812430669 -0.029963990233074476 0.9965658921198667 0.0342924217536074 -0.0091275791342953 0.0022818947835738265 -0.00038031579726230386; -0.00024483381428042826 0.0014690028856825695 -0.0058760115427302795 0.022035043285238542 -0.08226416159822382 0.9655246031076583 0.12535674916759093 -0.03284259977802422 0.00821064994450606 -0.001368441657417675; -0.00034442737698853074 0.002066564261931184 -0.008266257047724738 0.030998463928967768 -0.11572759866814622 0.9064639307436184 0.23141587569367308 -0.05887143351831279 0.01471785837957821 -0.002452976396596364; -0.00039561871994080617 0.002373712319644837 -0.009494849278579352 0.035605684794672555 -0.13292788990011076 0.8246148748057714 0.3472413906770251 -0.08465343751387333 0.021163359378468347 -0.0035272265630780526; -0.0004060673325934146 0.0024364039955604875 -0.009745615982241954 0.036546059933407316 -0.1364386237513872 0.7252084350721422 0.4676048834628187 -0.10762796892341828 0.026906992230854583 -0.0044844987051424244; -0.00038343270440251566 0.0023005962264150935 -0.009202384905660378 0.034508943396226406 -0.12883338867924515 0.6134756113207547 0.5872779433962262 -0.12523438490566025 0.031308596226415075 -0.005218099371069172; -0.0003353743248242692 0.0020122459489456154 -0.008048983795782463 0.030183689234184234 -0.11268577314095438 0.4946474033296336 0.7010321598224195 -0.13491204261931175 0.03372801065482796 -0.005621335109137986; -0.0002695516833148352 0.0016173100998890114 -0.0064692403995560465 0.02425965149833517 -0.09056936559378455 0.37395481087680327 0.8036391220865705 -0.13410029922308533 0.033525074805771354 -0.0055875124676285515; -0.0001936242693303734 0.00116174561598224 -0.004646982463928962 0.017426184239733605 -0.0650577544950054 0.2566288337402881 0.8898704195338514 -0.12023851187569348 0.030059627968923384 -0.005009937994820558; -0.0001152515723270436 0.0006915094339622616 -0.002766037735849047 0.010372641509433924 -0.03872452830188662 0.14790047169811266 0.9544976415094342 -0.0907660377358488 0.022691509433962215 -0.0037819182389936977; -4.2093081761005835e-5 0.000252558490566035 -0.00101023396226414 0.0037883773584905253 -0.014143275471697947 0.053000724528301305 0.9922923773584907 -0.043122233962263765 0.010780558490565948 -0.0017967597484276556; 1.8251646318905332e-5 -0.000109509877913432 0.00043803951165372805 -0.0016426481687014797 0.006132553163152189 -0.022887564483907292 0.998090704772475 0.025205745394007245 -0.006294686348501815 0.001049114391416968; 6.177906030336693e-5 -0.00037067436182020155 0.0014826974472808064 -0.005560115427303023 0.020757764261931273 -0.07747094162042215 0.9705980022197555 0.11386293274139928 -0.028033733185349834 0.004672288864224967; 9.005941546429895e-5 -0.0005403564927857936 0.002161425971143175 -0.008105347391786904 0.030259963596004427 -0.1129345069922309 0.9145170643729194 0.21804924950055468 -0.052197062375138714 0.008699510395856444; 0.00010514243433222355 -0.0006308546059933411 0.0025234184239733655 -0.009462819089900116 0.03532785793562709 -0.13184861265260833 0.8350665926748062 0.3325822419533849 -0.07639556048834628 0.012732593414724365; 0.00010907783943766192 -0.0006544670366259715 0.002617868146503887 -0.009817005549389572 0.03665015405105439 -0.13678361065482808 0.7374652885682578 0.4522794563817979 -0.09824011409544955 0.016373352349241574; 0.00010391535331113588 -0.0006234921198668153 0.0024939684794672615 -0.009352381798002228 0.034915558712541636 -0.13030985305216441 0.6269318534961157 0.5719584390677025 -0.11534160976692573 0.019223601627820934; 9.170469848316696e-5 -0.0005502281908990017 0.002200912763596007 -0.008253422863485025 0.03081277869034408 -0.11499769189789139 0.5086849889012212 0.6864367362930078 -0.12531093407325206 0.020885155678875313; 7.449559748427679e-5 -0.0004469735849056607 0.0017878943396226429 -0.00670460377358491 0.025030520754716987 -0.0934174792452831 0.38794339622641516 0.7905318943396227 -0.12575897358490581 0.02095982893081761; 5.4337772844987095e-5 -0.0003260266370699226 0.0013041065482796905 -0.004890399556048837 0.01825749167591565 -0.06813956714761382 0.26992577691453945 0.8790614594894562 -0.11429661487236414 0.01904943581206067; 3.328094709581948e-5 -0.00019968568257491687 0.0007987427302996675 -0.0029952852386237517 0.011182398224195336 -0.04173430765815764 0.159850832408435 0.9468429780244175 -0.08853474450610442 0.014755790751017385; 1.3374842767295593e-5 -8.024905660377356e-5 0.00032099622641509423 -0.001203735849056603 0.004493947169811317 -0.016772052830188677 0.06293726415094333 0.9886939962264152 -0.04608424905660377 0.00768070817610062; -3.3355530891602136e-6 2.001331853496128e-5 -8.005327413984513e-5 0.0003001997780244192 -0.001120745837957831 0.004182783573806909 -0.015610388457269794 0.999450770255272 0.015431307436182127 -0.002570551239363685; -1.557417684054756e-5 9.344506104328536e-5 -0.00037378024417314155 0.0014016759156492806 -0.0052329234184239765 0.019530017758046646 -0.07288714761376255 0.9769875726970031 0.09632985682574942 -0.01583314280429155; -2.3677395486496498e-5 0.00014206437291897897 -0.0005682574916759161 0.002130965593784685 -0.007955604883462818 0.029691453940066618 -0.11081021087680357 0.9255493895671472 0.19261265260821353 -0.03076877543470221; -2.818460969293378e-5 0.00016910765815760266 -0.000676430632630411 0.0025366148723640403 -0.009470028856825744 0.03534350055493897 -0.13190397336293005 0.8501833928967808 0.29973740177580527 -0.04589140029596747; -2.963522012578614e-5 0.0001778113207547168 -0.0007112452830188676 0.002667169811320753 -0.009957433962264136 0.03716256603773583 -0.1386928301886791 0.7559367547169813 0.41316181132075425 -0.059714968553458994; -2.8568627450980375e-5 0.00017141176470588225 -0.0006856470588235293 0.0025711764705882343 -0.009599058823529402 0.035825058823529404 -0.13370117647058813 0.6478566470588237 0.5283435882352938 -0.0707534313725489; -2.5524232334443203e-5 0.00015314539400665923 -0.0006125815760266372 0.002297180910099889 -0.008576142064372912 0.03200738734739179 -0.11945340732519416 0.5309902419533852 0.6407404395116535 -0.07752073991860883; -2.1041435442101356e-5 0.00012624861265260812 -0.0005049944506104328 0.0018937291897891226 -0.007069922308546051 0.026385960044395106 -0.09847391786903432 0.4103847114317424 0.7458100721420643 -0.0785308453570106; -1.5659637439881597e-5 9.395782463928956e-5 -0.00037583129855715835 0.0014093673695893438 -0.005261638179800212 0.01963718534961153 -0.07328710321864584 0.291087227524972 0.839010193118757 -0.07229769885312604; -9.918238993710668e-6 5.950943396226401e-5 -0.00023803773584905608 0.0008926415094339604 -0.0033325283018867825 0.01243747169811318 -0.04641735849056591 0.17814496226415055 0.9157985094339623 -0.05733525157232689; -4.356640769515327e-6 2.6139844617091963e-5 -0.00010455937846836788 0.00039209766925637953 -0.0014638312985571491 0.005463227524972222 -0.020389078801331725 0.07660508768035472 0.9716327280799114 -0.032157454679984994; 4.859045504994717e-7 -2.91542730299683e-6 1.1661709211987325e-5 -4.373140954495245e-5 0.00016326392896782237 -0.0006093243063263376 0.0022740332963375265 -0.008486808879023776 1.001972202219756 0.0047211329633743; 4.217536071032212e-6 -2.530521642619327e-5 0.00010122086570477312 -0.000379578246392899 0.0014170921198668222 -0.005288790233074395 0.01973806881243074 -0.07366348501664864 1.003915871254162 0.05418068812430667; 6.87036625971145e-6 -4.12221975582687e-5 0.00016488879023307486 -0.0006183329633740304 0.0023084430632630457 -0.00861543928967816 0.03215331409544957 -0.11999781709212023 0.9792789542730298 0.11536034095449546; 8.552275249722546e-6 -5.131365149833527e-5 0.00020525460599334116 -0.000769704772475029 0.0028735644839067735 -0.010724553163152075 0.04002464816870149 -0.14937403951165404 0.9307195098779132 0.18708808168701496; 9.371143174250844e-6 -5.622685904550505e-5 0.0002249074361820203 -0.0008434028856825757 0.003148704106548281 -0.011751413540510559 0.04385695005549392 -0.16367638668146528 0.8608955966703657 0.2681919005549396; 9.434850166481693e-6 -5.660910099889015e-5 0.00022643640399556072 -0.0008491365149833522 0.0031701096559378465 -0.011831302108768045 0.0441550987791343 -0.1647890930077693 0.7724652732519428 0.35749978779134267; 8.85127635960045e-6 -5.31076581576027e-5 0.0002124306326304109 -0.0007966148723640404 0.0029740288568257493 -0.011099500554938966 0.04142397336293009 -0.15459639289678154 0.6680865982241957 0.4538397336293005; 7.72830188679246e-6 -4.6369811320754756e-5 0.0001854792452830191 -0.0006955471698113213 0.002596709433962265 -0.009691290566037746 0.0361684528301887 -0.13498252075471717 0.5504176301886796 0.5560397283018865; 6.1738068812430705e-6 -3.704284128745842e-5 0.00014817136514983374 -0.0005556426193118763 0.0020743991120976703 -0.007741953829078811 0.028893416204217554 -0.1078317109877915 0.4221164277469482 0.6629277620421752; 4.295671476137629e-6 -2.5774028856825774e-5 0.00010309611542730314 -0.0003866104328523866 0.0014433456159822424 -0.005386772031076588 0.020103742508324095 -0.07502819800221985 0.28584104950055517 0.7733318250832407; 2.2017758046614897e-6 -1.3210654827968936e-5 5.284261931187576e-5 -0.00019815982241953403 0.0007397966703662599 -0.002761026859045508 0.010304310765815766 -0.03845621620421759 0.1442495540510545 0.8860799076581575; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0]
S = [0.08057422395181214 -0.18273824594520433 0.12953878824948 -0.03470980928704722 0.009300448898708805 -0.0024919863077880046 0.0006674963324432159 -0.00017799902198485748 4.44997554962144e-5 -7.416625916035724e-6; -0.18273824594520433 0.49501528013988866 -0.47652563173121154 0.2082588557222833 -0.055802693392252814 0.014951917846728027 -0.004004977994659296 0.0010679941319091447 -0.0002669985329772864 4.449975549621434e-5; 0.12953878824948 -0.47652563173121154 0.703274135862172 -0.5323283251234647 0.22321077356901134 -0.05980767138691213 0.01601991197863719 -0.0042719765276365805 0.0010679941319091458 -0.00017799902198485743; -0.03470980928704722 0.2082588557222833 -0.5323283251234647 0.7182260537089002 -0.5363333031181238 0.22427876770092042 -0.06007466991988943 0.016019911978637173 -0.004004977994659296 0.0006674963324432152; 0.009300448898708805 -0.055802693392252814 0.22321077356901134 -0.5363333031181238 0.7192940478408089 -0.5366003016511008 0.22427876770092042 -0.05980767138691207 0.014951917846728028 -0.002491986307788002; -0.0024919863077880046 0.014951917846728027 -0.05980767138691213 0.22427876770092042 -0.5366003016511008 0.7192940478408092 -0.5363333031181239 0.22321077356901134 -0.05580269339225286 0.0093004488987088; 0.0006674963324432159 -0.004004977994659296 0.01601991197863719 -0.06007466991988943 0.22427876770092042 -0.5363333031181239 0.7182260537089 -0.5323283251234643 0.2082588557222833 -0.03470980928704717; -0.00017799902198485748 0.0010679941319091447 -0.0042719765276365805 0.016019911978637173 -0.05980767138691207 0.22321077356901134 -0.5323283251234643 0.7032741358621719 -0.4765256317312117 0.12953878824948; 4.44997554962144e-5 -0.0002669985329772864 0.0010679941319091458 -0.004004977994659296 0.014951917846728028 -0.05580269339225286 0.2082588557222833 -0.4765256317312117 0.49501528013988877 -0.18273824594520421; -7.416625916035724e-6 4.449975549621434e-5 -0.00017799902198485743 0.0006674963324432152 -0.002491986307788002 0.0093004488987088 -0.03470980928704717 0.12953878824948 -0.18273824594520421 0.08057422395181209]
A = [-15.588457269700335 12.530743618202001 3.8770255271920075 -1.0388457269700335 0.27835738068812405 -0.07458379578246387 0.019977802441731415 -0.005327413984461705 0.0013318534961154276 -0.00022197558268590426; 4.176914539400668 -25.061487236404 19.245948945615982 2.077691453940067 -0.5567147613762481 0.14916759156492773 -0.03995560488346283 0.01065482796892341 -0.0026637069922308553 0.0004439511653718085; -1.1192008879023312 6.715205327413981 -26.86082130965594 19.728079911209772 1.9485016648168694 -0.5220865704772472 0.13984461709211993 -0.03729189789123195 0.009322974472807992 -0.00155382907880133; 0.29988901220865716 -1.7993340732519418 7.197336293007772 -26.99001109877914 19.762708102108775 1.939178690344061 -0.5194228634850168 0.13851276359600437 -0.03462819089900111 0.005771365149833511; -0.08035516093229747 0.4821309655937846 -1.9285238623751395 7.231964483906773 -26.999334073251937 19.76537180910099 1.9378468368479478 -0.5167591564927856 0.1291897891231965 -0.021531631520532724; 0.021531631520532737 -0.12918978912319637 0.5167591564927858 -1.937846836847947 7.234628190898997 -27.000665926748063 19.76803551609324 1.9285238623751384 -0.48213096559378493 0.08035516093229739; -0.00577136514983352 0.034628190899001106 -0.1385127635960045 0.5194228634850169 -1.9391786903440615 7.237291897891236 -27.009988901220865 19.802663706992224 1.7993340732519425 -0.2998890122086566; 0.0015538290788013316 -0.009322974472807983 0.037291897891231954 -0.13984461709211984 0.522086570477247 -1.94850166481687 7.271920088790229 -27.139178690344067 20.284794672586024 1.1192008879023299; -0.0004439511653718092 0.0026637069922308535 -0.01065482796892342 0.03995560488346282 -0.14916759156492776 0.5567147613762488 -2.0776914539400657 7.754051054384021 -28.938512763596 22.82308546059933; -22.82308546059933 28.938512763595998 -7.754051054384015 2.077691453940067 -0.5567147613762481 0.14916759156492773 -0.03995560488346283 0.01065482796892341 -0.0026637069922308553 0.0004439511653718085; -1.1192008879023312 -20.284794672586017 27.139178690344064 -7.271920088790236 1.9485016648168694 -0.5220865704772472 0.13984461709211993 -0.03729189789123195 0.009322974472807992 -0.00155382907880133; 0.29988901220865716 -1.7993340732519418 -19.802663706992234 27.00998890122087 -7.237291897891227 1.939178690344061 -0.5194228634850168 0.13851276359600437 -0.03462819089900111 0.005771365149833511; -0.08035516093229747 0.4821309655937846 -1.9285238623751395 -19.768035516093228 27.000665926748052 -7.234628190898997 1.9378468368479478 -0.5167591564927856 0.1291897891231965 -0.021531631520532724; 0.021531631520532737 -0.12918978912319637 0.5167591564927858 -1.937846836847947 -19.76537180910099 26.99933407325194 -7.231964483906776 1.9285238623751384 -0.48213096559378493 0.08035516093229739; -0.00577136514983352 0.034628190899001106 -0.1385127635960045 0.5194228634850169 -1.9391786903440615 -19.76270810210878 26.99001109877914 -7.197336293007766 1.7993340732519425 -0.2998890122086566; 0.0015538290788013316 -0.009322974472807983 0.037291897891231954 -0.13984461709211984 0.522086570477247 -1.94850166481687 -19.72807991120976 26.860821309655936 -6.715205327413989 1.1192008879023299; -0.0004439511653718092 0.0026637069922308535 -0.01065482796892342 0.03995560488346282 -0.14916759156492776 0.5567147613762488 -2.0776914539400657 -19.245948945615993 25.061487236404 -4.17691453940066; 0.00022197558268590461 -0.0013318534961154268 0.0053274139844617105 -0.019977802441731415 0.0745837957824639 -0.2783573806881244 1.038845726970033 -3.877025527192011 -12.530743618201987 15.588457269700324; -1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 -1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 -1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 -1.0 1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 -1.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 -1.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 1.0; -11.411542730299665 14.469256381797999 -3.8770255271920075 1.0388457269700335 -0.27835738068812405 0.07458379578246387 -0.019977802441731415 0.005327413984461705 -0.0013318534961154276 0.00022197558268590426; -4.176914539400668 -1.9385127635959987 7.754051054384015 -2.077691453940067 0.5567147613762481 -0.14916759156492773 0.03995560488346283 -0.01065482796892341 0.0026637069922308553 -0.0004439511653718085; 1.1192008879023312 -6.715205327413981 -0.1391786903440682 7.271920088790236 -1.9485016648168694 0.5220865704772472 -0.13984461709211993 0.03729189789123195 -0.009322974472807992 0.00155382907880133; -0.29988901220865716 1.7993340732519418 -7.197336293007772 -0.009988901220859612 7.237291897891227 -1.939178690344061 0.5194228634850168 -0.13851276359600437 0.03462819089900111 -0.005771365149833511; 0.08035516093229747 -0.4821309655937846 1.9285238623751395 -7.231964483906773 -0.0006659267480504719 7.234628190898997 -1.9378468368479478 0.5167591564927856 -0.1291897891231965 0.021531631520532724; -0.021531631520532737 0.12918978912319637 -0.5167591564927858 1.937846836847947 -7.234628190898997 0.0006659267480471557 7.231964483906776 -1.9285238623751384 0.48213096559378493 -0.08035516093229739; 0.00577136514983352 -0.034628190899001106 0.1385127635960045 -0.5194228634850169 1.9391786903440615 -7.237291897891236 0.009988901220875031 7.197336293007766 -1.7993340732519425 0.2998890122086566; -0.0015538290788013316 0.009322974472807983 -0.037291897891231954 0.13984461709211984 -0.522086570477247 1.94850166481687 -7.271920088790229 0.13917869034405173 6.715205327413989 -1.1192008879023299; 0.0004439511653718092 -0.0026637069922308535 0.01065482796892342 -0.03995560488346282 0.14916759156492776 -0.5567147613762488 2.0776914539400657 -7.754051054384021 1.9385127635960142 4.17691453940066; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0]
b = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0]
C = [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
# needed to build prediction matrix https://rdrr.io/cran/mgcv/man/Predict.matrix.html
function model_matrix_cr(xk, F, x = nothing)
F = transpose(F)
if isnothing(x)
x = xk
end
k = length(xk)
n = length(x)
# model matrix
X = zeros(Float64, n, k)
for i in 1:n
# find interval containing x[i] or extrapolate
if (x[i] < xk[1] || x[i] > xk[end])
# extrapolate
if x[i] < xk[1]
# below knot range
hj = xk[2] - xk[1]
xik = x[i] - xk[1]
cjm = -xik*hj/3
cjp = -xik*hj/6
# set model matrix elements
X[i,:] = (cjm * F[:,1]) + (cjp * F[:,2])
X[i,1] += 1 - xik/hj
X[i,2] += xik/hj
else
# above knot range
hj = xk[end] - xk[end-1]
xik = x[i] - xk[end]
cjm = xik*hj/6
cjp = xik*hj/3
X[i,:] = (cjm * F[:,end-1]) + (cjp * F[:,end])
X[i,end-1] += -xik/hj
X[i,end] += 1 + xik/hj
end
else
# inside knot sequence
j = searchsortedlast(xk, x[i])
j = min(j,k-1) # j is the start of the interval, so cannot be greater than k-1
# if j == k
# j -= 1
# end
hj = xk[j+1] - xk[j] # interval width
ajm = (xk[j+1] - x[i]) # basis fn's; table 5.1 in GAM book
ajp = (x[i] - xk[j])
cjm = (ajm*(ajm*ajm/hj - hj))/6
cjp = (ajp*(ajp*ajp/hj - hj))/6
ajm /= hj
ajp /= hj
# set i-th row of X
X[i,:] = (cjm * F[:,j]) + (cjp * F[:,j+1])
X[i,j] += ajm
X[i,j+1] += ajp
end
end
return X
end
# ------------------------------------------------------------
# SIR model with nonlinear force of infection
# we know the true nonlinearity and get some training data
function sir!(du, u, p, t)
du[1] = -p.β*u[1]*(u[2]^p.α)
du[2] = p.β*u[1]*(u[2]^p.α) - p.γ*u[2]
du[3] = p.γ*u[2]
end
N = 1000
sir_pars = (β=0.7, γ=0.25, α=0.7)
u0 = [0.99, 0.01, 0]
tspan = (0.0, 40.0)
sir_prob = DE.ODEProblem(sir!, u0, tspan, sir_pars)
sir_sol = DE.solve(sir_prob)
# training data
training_data = let
time=tspan[1]:1:30
cdata = diff(sir_sol.(time, idxs=3))
noisy_data = rand.(Poisson.(N .* cdata))
DataFrame(time=time[2:end], cdata=cdata, noisy_data=noisy_data)
end
# ------------------------------------------------------------
# SIR model; we want to estimate force of infection using constrained
# basis fn approach
"""
I: the place to evaluate the spline
p: spline coefficients
F: matrix from `basis_cr`
xk: knot points
"""
function foi_eval(I, p, F, xk)
only(model_matrix_cr(xk, F, [I]) * p)
end
function sir_psm!(du, u, p, t, pars)
foi = foi_eval(u[2], p, pars.F, pars.xk)
du[1] = -foi*u[1]
du[2] = foi*u[2] - pars.γ*u[2]
du[3] = pars.γ*u[2]
end
function get_feasible_theta(A, b, C)
k = size(A,2)
θ = [0; cumsum(rand(Uniform(0,0.1), k-1))]
while any(A * θ .< b) || any(C * θ .!= zeros(size(C,1)))
θ = [0; cumsum(rand(Uniform(0,0.1), k-1))]
end
return θ
end
θ0 = get_feasible_theta(A, b, C)
static_pars = (sir_pars..., F=F, xk=xk, A=A, b=b, C=C) # wont change over entire fitting procedure
sir_psm_prob = DE.ODEProblem((du, u, p, t) -> sir_psm!(du, u, p, t, static_pars), u0, tspan, θ0)
function sir_pnll(θ, Sb)
prob = DE.remake(sir_psm_prob, p = θ)
soln = DE.solve(prob, saveat=1, save_idxs=3)
pred = N .* diff(soln.u[1:31])
ll = sum([logpdf(Poisson(pred[i]), training_data.noisy_data[i]) for i in eachindex(pred)])
ll -= only(transpose(θ) * Sb * θ)
return -ll
end
# pars: static pars
function cons(res, u, _, pars)
icons = pars.A*u - pars.b
econs = pars.C*u
res[1:size(A,1)] .= icons
res[size(A,1)+1:end] .= econs
end
Sb = S * 0.5
sir_optfn = OptimizationFunction((u, _) -> sir_pnll(u, Sb), AutoMooncake(), cons = (res, u, p) -> cons(res, u, p, static_pars))
sir_optprob = OptimizationProblem(
sir_optfn, θ0,
lcons=zeros(size(A,1) + size(C,1)), ucons=[fill(Inf, size(A,1)); zeros(size(C, 1))]
)
# error here
sol = solve(sir_optprob, OptimizationLBFGSB.LBFGSB())
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment