text:
@inline function initialize!(integrator,cache::Feagin10ConstantCache,f=integrator.f)
  integrator.fsalfirst = f(integrator.t,integrator.uprev) # Pre-start fsal
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
end

@inline function perform_step!(integrator,cache::Feagin10ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  @unpack adaptiveConst,a0100,a0200,a0201,a0300,a0302,a0400,a0402,a0403,a0500,a0503,a0504,a0600,a0603,a0604,a0605,a0700,a0704,a0705,a0706,a0800,a0805,a0806,a0807,a0900,a0905,a0906,a0907,a0908,a1000,a1005,a1006,a1007,a1008,a1009,a1100,a1105,a1106,a1107,a1108,a1109,a1110,a1200,a1203,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1300,a1302,a1303,a1305,a1306,a1307,a1308,a1309,a1310,a1311,a1312,a1400,a1401,a1404,a1406,a1412,a1413,a1500,a1502,a1514,a1600,a1601,a1602,a1604,a1605,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1613,a1614,a1615,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16 = cache
  k1  = integrator.fsalfirst
  a = dt*a0100
  k2  = f(@muladd(@muladd(t + c1*dt)),@muladd(uprev  + a*k1))
  k3  = f(@muladd(@muladd(t + c2*dt)),@muladd(uprev + dt*(a0200*k1 + a0201*k2)))
  k4  = f(@muladd(@muladd(t + c3*dt)),@muladd(uprev  + dt*(a0300*k1              + a0302*k3)))
  k5  = f(@muladd(@muladd(t + c4*dt)),@muladd(uprev  + dt*(a0400*k1              + a0402*k3 + a0403*k4)))
  k6  = f(@muladd(@muladd(t + c5*dt)),@muladd(uprev  + dt*(a0500*k1                           + a0503*k4 + a0504*k5)))
  k7  = f(@muladd(@muladd(t + c6*dt)),@muladd(uprev  + dt*(a0600*k1                           + a0603*k4 + a0604*k5 + a0605*k6)))
  k8  = f(@muladd(@muladd(t + c7*dt)),@muladd(uprev  + dt*(a0700*k1                                        + a0704*k5 + a0705*k6 + a0706*k7)))
  k9  = f(@muladd(@muladd(t + c8*dt)),@muladd(uprev  + dt*(a0800*k1                                                     + a0805*k6 + a0806*k7 + a0807*k8)))
  k10 = f(@muladd(@muladd(t + c9*dt)),@muladd(uprev  + dt*(a0900*k1                                                     + a0905*k6 + a0906*k7 + a0907*k8 + a0908*k9)))
  k11 = f(@muladd(@muladd(t + c10*dt)),@muladd(uprev + dt*(a1000*k1                                                     + a1005*k6 + a1006*k7 + a1007*k8 + a1008*k9 + a1009*k10)))
  k12 = f(@muladd(@muladd(t + c11*dt)),@muladd(uprev + dt*(a1100*k1                                                     + a1105*k6 + a1106*k7 + a1107*k8 + a1108*k9 + a1109*k10 + a1110*k11)))
  k13 = f(@muladd(@muladd(t + c12*dt)),@muladd(uprev + dt*(a1200*k1                           + a1203*k4 + a1204*k5 + a1205*k6 + a1206*k7 + a1207*k8 + a1208*k9 + a1209*k10 + a1210*k11 + a1211*k12)))
  k14 = f(@muladd(@muladd(t + c13*dt)),@muladd(uprev + dt*(a1300*k1              + a1302*k3 + a1303*k4              + a1305*k6 + a1306*k7 + a1307*k8 + a1308*k9 + a1309*k10 + a1310*k11 + a1311*k12 + a1312*k13)))
  k15 = f(@muladd(@muladd(t + c14*dt)),@muladd(uprev + dt*(a1400*k1 + a1401*k2                           + a1404*k5              + a1406*k7 +                                                                     a1412*k13 + a1413*k14)))
  k16 = f(@muladd(@muladd(t + c15*dt)),@muladd(uprev + dt*(a1500*k1              + a1502*k3                                                                                + a1514*k15)))
  k17 = f(@muladd(@muladd(t + c16*dt)),@muladd(uprev + dt*(a1600*k1 + a1601*k2 + a1602*k3              + a1604*k5 + a1605*k6 + a1606*k7 + a1607*k8 + a1608*k9 + a1609*k10 + a1610*k11 + a1611*k12 + a1612*k13 + a1613*k14 + a1614*k15 + a1615*k16)))
  u = @muladd uprev + dt*(b1*k1 + b2*k2 + b3*k3 + b5*k5 + b7*k7 + b9*k9 + b10*k10 + b11*k11 + b12*k12 + b13*k13 + b14*k14 + b15*k15 + b16*k16 + b17*k17)
  if integrator.opts.adaptive
    integrator.EEst = integrator.opts.internalnorm((dt*(k2 - k16) * adaptiveConst)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
  end
  k = f(t+dt,u) # For the interpolation, needs k at the updated point
  integrator.fsallast = k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  @pack integrator = t,dt,u
end

@inline function initialize!(integrator,cache::Feagin10Cache,f=integrator.f)
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  f(integrator.t,integrator.uprev,integrator.fsalfirst) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::Feagin10Cache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack adaptiveConst,a0100,a0200,a0201,a0300,a0302,a0400,a0402,a0403,a0500,a0503,a0504,a0600,a0603,a0604,a0605,a0700,a0704,a0705,a0706,a0800,a0805,a0806,a0807,a0900,a0905,a0906,a0907,a0908,a1000,a1005,a1006,a1007,a1008,a1009,a1100,a1105,a1106,a1107,a1108,a1109,a1110,a1200,a1203,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1300,a1302,a1303,a1305,a1306,a1307,a1308,a1309,a1310,a1311,a1312,a1400,a1401,a1404,a1406,a1412,a1413,a1500,a1502,a1514,a1600,a1601,a1602,a1604,a1605,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1613,a1614,a1615,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16 = cache.tab
  @unpack k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,tmp,atmp,uprev,k = cache
  k1 = cache.fsalfirst
  a =  dt*a0100
  for i in uidx
    tmp[i] = @muladd uprev[i] + a*k1[i]
  end
  f(@muladd(t + c1*dt),tmp,k2)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0200*k1[i] + a0201*k2[i])
  end
  f(@muladd(t + c2*dt) ,tmp,k3)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0300*k1[i] + a0302*k3[i])
  end
  f(@muladd(t + c3*dt),tmp,k4)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0400*k1[i] + a0402*k3[i] + a0403*k4[i])
  end
  f(@muladd(t + c4*dt),tmp,k5)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0500*k1[i] + a0503*k4[i] + a0504*k5[i])
  end
  f(@muladd(t + c5*dt),tmp,k6)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0600*k1[i] + a0603*k4[i] + a0604*k5[i] + a0605*k6[i])
  end
  f(@muladd(t + c6*dt),tmp,k7)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0700*k1[i] + a0704*k5[i] + a0705*k6[i] + a0706*k7[i])
  end
  f(@muladd(t + c7*dt),tmp,k8)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0800*k1[i] + a0805*k6[i] + a0806*k7[i] + a0807*k8[i])
  end
  f(@muladd(t + c8*dt),tmp,k9)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0900*k1[i] + a0905*k6[i] + a0906*k7[i] + a0907*k8[i] + a0908*k9[i])
  end
  f(@muladd(t + c9*dt),tmp,k10)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1000*k1[i] + a1005*k6[i] + a1006*k7[i] + a1007*k8[i] + a1008*k9[i] + a1009*k10[i])
  end
  f(@muladd(t + c10*dt),tmp,k11)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1100*k1[i] + a1105*k6[i] + a1106*k7[i] + a1107*k8[i] + a1108*k9[i] + a1109*k10[i] + a1110*k11[i])
  end
  f(@muladd(t + c11*dt),tmp,k12)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1200*k1[i] + a1203*k4[i] + a1204*k5[i] + a1205*k6[i] + a1206*k7[i] + a1207*k8[i] + a1208*k9[i] + a1209*k10[i] + a1210*k11[i] + a1211*k12[i])
  end
  f(@muladd(t + c12*dt),tmp,k13)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1300*k1[i] + a1302*k3[i] + a1303*k4[i] + a1305*k6[i] + a1306*k7[i] + a1307*k8[i] + a1308*k9[i] + a1309*k10[i] + a1310*k11[i] + a1311*k12[i] + a1312*k13[i])
  end
  f(@muladd(t + c13*dt),tmp,k14)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1400*k1[i] + a1401*k2[i] + a1404*k5[i] + a1406*k7[i] + a1412*k13[i] + a1413*k14[i])
  end
  f(@muladd(t + c14*dt),tmp,k15)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1500*k1[i] + a1502*k3[i] + a1514*k15[i])
  end
  f(@muladd(t + c15*dt),tmp,k16)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1600*k1[i] + a1601*k2[i] + a1602*k3[i] + a1604*k5[i] + a1605*k6[i] + a1606*k7[i] + a1607*k8[i] + a1608*k9[i] + a1609*k10[i] + a1610*k11[i] + a1611*k12[i] + a1612*k13[i] + a1613*k14[i] + a1614*k15[i] + a1615*k16[i])
  end
  f(@muladd(t + c16*dt),tmp,k17)
  for i in uidx
    u[i] = @muladd uprev[i] + dt*(b1*k1[i] + b2*k2[i] + b3*k3[i] + b5*k5[i] + b7*k7[i] + b9*k9[i] + b10*k10[i] + b11*k11[i] + b12*k12[i] + b13*k13[i] + b14*k14[i] + b15*k15[i] + b16*k16[i] + b17*k17[i])
  end
  if integrator.opts.adaptive
    for i in uidx
      atmp[i] = (dt*(k2[i] - k16[i]) * adaptiveConst)./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
    end
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  f(t+dt,u,integrator.fsallast) # For the interpolation, needs k at the updated point
  @pack integrator = t,dt,u
end

@inline function initialize!(integrator,cache::Feagin12ConstantCache,f=integrator.f)
  integrator.fsalfirst = f(integrator.t,integrator.uprev) # Pre-start fsal
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
end

@inline function perform_step!(integrator,cache::Feagin12ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  @unpack adaptiveConst,a0100,a0200,a0201,a0300,a0302,a0400,a0402,a0403,a0500,a0503,a0504,a0600,a0603,a0604,a0605,a0700,a0704,a0705,a0706,a0800,a0805,a0806,a0807,a0900,a0905,a0906,a0907,a0908,a1000,a1005,a1006,a1007,a1008,a1009,a1100,a1105,a1106,a1107,a1108,a1109,a1110,a1200,a1208,a1209,a1210,a1211,a1300,a1308,a1309,a1310,a1311,a1312,a1400,a1408,a1409,a1410,a1411,a1412,a1413,a1500,a1508,a1509,a1510,a1511,a1512,a1513,a1514,a1600,a1608,a1609,a1610,a1611,a1612,a1613,a1614,a1615,a1700,a1705,a1706,a1707,a1708,a1709,a1710,a1711,a1712,a1713,a1714,a1715,a1716,a1800,a1805,a1806,a1807,a1808,a1809,a1810,a1811,a1812,a1813,a1814,a1815,a1816,a1817,a1900,a1904,a1905,a1906,a1908,a1909,a1910,a1911,a1912,a1913,a1914,a1915,a1916,a1917,a1918,a2000,a2003,a2004,a2005,a2007,a2009,a2010,a2017,a2018,a2019,a2100,a2102,a2103,a2106,a2107,a2109,a2110,a2117,a2118,a2119,a2120,a2200,a2201,a2204,a2206,a2220,a2221,a2300,a2302,a2322,a2400,a2401,a2402,a2404,a2406,a2407,a2408,a2409,a2410,a2411,a2412,a2413,a2414,a2415,a2416,a2417,a2418,a2419,a2420,a2421,a2422,a2423,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19,c20,c21,c22,c23,c24,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18,b19,b20,b21,b22,b23,b24,b25 = cache
  k1   = integrator.fsalfirst
  a =  dt*a0100
  k2  = f(@muladd(t + c1*dt),@muladd(uprev  + a*k1))
  k3  = f(@muladd(t + c2*dt),@muladd(uprev + dt*(a0200*k1 + a0201*k2)))
  k4  = f(@muladd(t + c3*dt),@muladd(uprev  + dt*(a0300*k1              + a0302*k3)))
  k5  = f(@muladd(t + c4*dt),@muladd(uprev  + dt*(a0400*k1              + a0402*k3 + a0403*k4)))
  k6  = f(@muladd(t + c5*dt),@muladd(uprev  + dt*(a0500*k1                           + a0503*k4 + a0504*k5)))
  k7  = f(@muladd(t + c6*dt),@muladd(uprev  + dt*(a0600*k1                           + a0603*k4 + a0604*k5 + a0605*k6)))
  k8  = f(@muladd(t + c7*dt),@muladd(uprev  + dt*(a0700*k1                                        + a0704*k5 + a0705*k6 + a0706*k7)))
  k9  = f(@muladd(t + c8*dt),@muladd(uprev  + dt*(a0800*k1                                                     + a0805*k6 + a0806*k7 + a0807*k8)))
  k10 = f(@muladd(t + c9*dt),@muladd(uprev  + dt*(a0900*k1                                                     + a0905*k6 + a0906*k7 + a0907*k8 + a0908*k9)))
  k11 = f(@muladd(t + c10*dt),@muladd(uprev + dt*(a1000*k1                                                     + a1005*k6 + a1006*k7 + a1007*k8 + a1008*k9 + a1009*k10)))
  k12 = f(@muladd(t + c11*dt),@muladd(uprev + dt*(a1100*k1                                                     + a1105*k6 + a1106*k7 + a1107*k8 + a1108*k9 + a1109*k10 + a1110*k11)))
  k13 = f(@muladd(t + c12*dt),@muladd(uprev + dt*(a1200*k1                                                                                + a1208*k9 + a1209*k10 + a1210*k11 + a1211*k12)))
  k14 = f(@muladd(t + c13*dt),@muladd(uprev + dt*(a1300*k1                                                                                + a1308*k9 + a1309*k10 + a1310*k11 + a1311*k12 + a1312*k13)))
  k15 = f(@muladd(t + c14*dt),@muladd(uprev + dt*(a1400*k1                                                                                + a1408*k9 + a1409*k10 + a1410*k11 + a1411*k12 + a1412*k13 + a1413*k14)))
  k16 = f(@muladd(t + c15*dt),@muladd(uprev + dt*(a1500*k1                                                                                + a1508*k9 + a1509*k10 + a1510*k11 + a1511*k12 + a1512*k13 + a1513*k14 + a1514*k15)))
  k17 = f(@muladd(t + c16*dt),@muladd(uprev + dt*(a1600*k1                                                                                + a1608*k9 + a1609*k10 + a1610*k11 + a1611*k12 + a1612*k13 + a1613*k14 + a1614*k15 + a1615*k16)))
  k18 = f(@muladd(t + c17*dt),@muladd(uprev + dt*(a1700*k1                                                     + a1705*k6 + a1706*k7 + a1707*k8 + a1708*k9 + a1709*k10 + a1710*k11 + a1711*k12 + a1712*k13 + a1713*k14 + a1714*k15 + a1715*k16 + a1716*k17)))
  k19 = f(@muladd(t + c18*dt),@muladd(uprev + dt*(a1800*k1                                                     + a1805*k6 + a1806*k7 + a1807*k8 + a1808*k9 + a1809*k10 + a1810*k11 + a1811*k12 + a1812*k13 + a1813*k14 + a1814*k15 + a1815*k16 + a1816*k17 + a1817*k18)))
  k20 = f(@muladd(t + c19*dt),@muladd(uprev + dt*(a1900*k1                                        + a1904*k5 + a1905*k6 + a1906*k7              + a1908*k9 + a1909*k10 + a1910*k11 + a1911*k12 + a1912*k13 + a1913*k14 + a1914*k15 + a1915*k16 + a1916*k17 + a1917*k18 + a1918*k19)))
  k21 = f(@muladd(t + c20*dt),@muladd(uprev + dt*(a2000*k1                           + a2003*k4 + a2004*k5 + a2005*k6              + a2007*k8              + a2009*k10 + a2010*k11                                                                                + a2017*k18 + a2018*k19 + a2019*k20)))
  k22 = f(@muladd(t + c21*dt),@muladd(uprev + dt*(a2100*k1              + a2102*k3 + a2103*k4                           + a2106*k7 + a2107*k8              + a2109*k10 + a2110*k11                                                                                + a2117*k18 + a2118*k19 + a2119*k20 + a2120*k21)))
  k23 = f(@muladd(t + c22*dt),@muladd(uprev + dt*(a2200*k1 + a2201*k2                           + a2204*k5              + a2206*k7                                                                                + a2220*k21 + a2221*k22)))
  k24 = f(@muladd(t + c23*dt),@muladd(uprev + dt*(a2300*k1              + a2302*k3                                                                                + a2322*k23)))
  k25 = f(@muladd(t + c24*dt),@muladd(uprev + dt*(a2400*k1 + a2401*k2 + a2402*k3              + a2404*k5              + a2406*k7 + a2407*k8 + a2408*k9 + a2409*k10 + a2410*k11 + a2411*k12 + a2412*k13 + a2413*k14 + a2414*k15 + a2415*k16 + a2416*k17 + a2417*k18 + a2418*k19 + a2419*k20 + a2420*k21 + a2421*k22 + a2422*k23 + a2423*k24)))

  u = @muladd uprev + dt*(b1*k1 + b2*k2 + b3*k3 + b5*k5 + b7*k7 + b8*k8 + b10*k10 + b11*k11 + b13*k13 + b14*k14 + b15*k15 + b16*k16 + b17*k17 + b18*k18 + b19*k19 + b20*k20 + b21*k21 + b22*k22 + b23*k23 + b24*k24 + b25*k25)
  k = f(t+dt,u)
  integrator.fsallast = k
  if integrator.opts.adaptive
    integrator.EEst = integrator.opts.internalnorm((dt*(k2 - k24) * adaptiveConst)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  @pack integrator = t,dt,u
end

@inline function initialize!(integrator,cache::Feagin12Cache,f=integrator.f)
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  f(integrator.t,integrator.uprev,integrator.fsalfirst) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::Feagin12Cache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack adaptiveConst,a0100,a0200,a0201,a0300,a0302,a0400,a0402,a0403,a0500,a0503,a0504,a0600,a0603,a0604,a0605,a0700,a0704,a0705,a0706,a0800,a0805,a0806,a0807,a0900,a0905,a0906,a0907,a0908,a1000,a1005,a1006,a1007,a1008,a1009,a1100,a1105,a1106,a1107,a1108,a1109,a1110,a1200,a1208,a1209,a1210,a1211,a1300,a1308,a1309,a1310,a1311,a1312,a1400,a1408,a1409,a1410,a1411,a1412,a1413,a1500,a1508,a1509,a1510,a1511,a1512,a1513,a1514,a1600,a1608,a1609,a1610,a1611,a1612,a1613,a1614,a1615,a1700,a1705,a1706,a1707,a1708,a1709,a1710,a1711,a1712,a1713,a1714,a1715,a1716,a1800,a1805,a1806,a1807,a1808,a1809,a1810,a1811,a1812,a1813,a1814,a1815,a1816,a1817,a1900,a1904,a1905,a1906,a1908,a1909,a1910,a1911,a1912,a1913,a1914,a1915,a1916,a1917,a1918,a2000,a2003,a2004,a2005,a2007,a2009,a2010,a2017,a2018,a2019,a2100,a2102,a2103,a2106,a2107,a2109,a2110,a2117,a2118,a2119,a2120,a2200,a2201,a2204,a2206,a2220,a2221,a2300,a2302,a2322,a2400,a2401,a2402,a2404,a2406,a2407,a2408,a2409,a2410,a2411,a2412,a2413,a2414,a2415,a2416,a2417,a2418,a2419,a2420,a2421,a2422,a2423,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19,c20,c21,c22,c23,c24,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18,b19,b20,b21,b22,b23,b24,b25 = cache.tab
  @unpack k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,k18,k19,k20,k21,k22,k23,k24,k25,tmp,atmp,uprev,k = cache
  k1 = cache.fsalfirst
  a = dt*a0100
  for i in uidx
    tmp[i] = @muladd uprev[i] + a*k1[i]
  end
  f(@muladd(t + c1*dt),tmp,k2)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0200*k1[i] + a0201*k2[i])
  end
  f(@muladd(t + c2*dt) ,tmp,k3)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0300*k1[i] + a0302*k3[i])
  end
  f(@muladd(t + c3*dt),tmp,k4)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0400*k1[i] + a0402*k3[i] + a0403*k4[i])
  end
  f(@muladd(t + c4*dt),tmp,k5)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0500*k1[i] + a0503*k4[i] + a0504*k5[i])
  end
  f(@muladd(t + c5*dt),tmp,k6)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0600*k1[i] + a0603*k4[i] + a0604*k5[i] + a0605*k6[i])
  end
  f(@muladd(t + c6*dt),tmp,k7)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0700*k1[i] + a0704*k5[i] + a0705*k6[i] + a0706*k7[i])
  end
  f(@muladd(t + c7*dt),tmp,k8)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0800*k1[i] + a0805*k6[i] + a0806*k7[i] + a0807*k8[i])
  end
  f(@muladd(t + c8*dt),tmp,k9)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0900*k1[i] + a0905*k6[i] + a0906*k7[i] + a0907*k8[i] + a0908*k9[i])
  end
  f(@muladd(t + c9*dt),tmp,k10)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1000*k1[i] + a1005*k6[i] + a1006*k7[i] + a1007*k8[i] + a1008*k9[i] + a1009*k10[i])
  end
  f(@muladd(t + c10*dt),tmp,k11)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1100*k1[i] + a1105*k6[i] + a1106*k7[i] + a1107*k8[i] + a1108*k9[i] + a1109*k10[i] + a1110*k11[i])
  end
  f(@muladd(t + c11*dt),tmp,k12)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1200*k1[i] + a1208*k9[i] + a1209*k10[i] + a1210*k11[i] + a1211*k12[i])
  end
  f(@muladd(t + c12*dt),tmp,k13)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1300*k1[i] + a1308*k9[i] + a1309*k10[i] + a1310*k11[i] + a1311*k12[i] + a1312*k13[i])
  end
  f(@muladd(t + c13*dt),tmp,k14)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1400*k1[i] + a1408*k9[i] + a1409*k10[i] + a1410*k11[i] + a1411*k12[i] + a1412*k13[i] + a1413*k14[i])
  end
  f(@muladd(t + c14*dt),tmp,k15)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1500*k1[i] + a1508*k9[i] + a1509*k10[i] + a1510*k11[i] + a1511*k12[i] + a1512*k13[i] + a1513*k14[i] + a1514*k15[i])
  end
  f(@muladd(t + c15*dt),tmp,k16)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*((a1600*k1[i] + a1608*k9[i] + a1609*k10[i]) + (a1610*k11[i] + a1611*k12[i] + a1612*k13[i] + a1613*k14[i]) + (a1614*k15[i] + a1615*k16[i]))
  end
  f(@muladd(t + c16*dt),tmp,k17)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*((a1700*k1[i] + a1705*k6[i] + a1706*k7[i]) + (a1707*k8[i] + a1708*k9[i] + a1709*k10[i] + a1710*k11[i]) + (a1711*k12[i] + a1712*k13[i] + a1713*k14[i] + a1714*k15[i]) + (a1715*k16[i] + a1716*k17[i]))
  end
  f(@muladd(t + c17*dt),tmp,k18)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*((a1800*k1[i] + a1805*k6[i] + a1806*k7[i]) + (a1807*k8[i] + a1808*k9[i] + a1809*k10[i] + a1810*k11[i]) + (a1811*k12[i] + a1812*k13[i] + a1813*k14[i] + a1814*k15[i]) + (a1815*k16[i] + a1816*k17[i] + a1817*k18[i]))
  end
  f(@muladd(t + c18*dt),tmp,k19)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*((a1900*k1[i] + a1904*k5[i] + a1905*k6[i]) + (a1906*k7[i] + a1908*k9[i] + a1909*k10[i] + a1910*k11[i]) + (a1911*k12[i] + a1912*k13[i] + a1913*k14[i] + a1914*k15[i]) + (a1915*k16[i] + a1916*k17[i] + a1917*k18[i] + a1918*k19[i]))
  end
  f(@muladd(t + c19*dt),tmp,k20)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*((a2000*k1[i] + a2003*k4[i] + a2004*k5[i]) + (a2005*k6[i] + a2007*k8[i] + a2009*k10[i] + a2010*k11[i]) + (a2017*k18[i] + a2018*k19[i] + a2019*k20[i]))
  end
  f(@muladd(t + c20*dt),tmp,k21)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*((a2100*k1[i] + a2102*k3[i] + a2103*k4[i]) + (a2106*k7[i] + a2107*k8[i] + a2109*k10[i] + a2110*k11[i]) + (a2117*k18[i] + a2118*k19[i] + a2119*k20[i] + a2120*k21[i]))
  end
  f(@muladd(t + c21*dt),tmp,k22)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*((a2200*k1[i] + a2201*k2[i] + a2204*k5[i]) + (a2206*k7[i] + a2220*k21[i] + a2221*k22[i]))
  end
  f(@muladd(t + c22*dt),tmp,k23)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a2300*k1[i] + a2302*k3[i] + a2322*k23[i])
  end
  f(@muladd(t + c23*dt),tmp,k24)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*((a2400*k1[i] + a2401*k2[i] + a2402*k3[i]) + (a2404*k5[i] + a2406*k7[i] + a2407*k8[i] + a2408*k9[i]) + (a2409*k10[i] + a2410*k11[i] + a2411*k12[i] + a2412*k13[i]) + (a2413*k14[i] + a2414*k15[i] + a2415*k16[i] + a2416*k17[i]) + (a2417*k18[i] + a2418*k19[i] + a2419*k20[i] + a2420*k21[i]) + (a2421*k22[i] + a2422*k23[i] + a2423*k24[i]))
  end
  f(@muladd(t + c24*dt),tmp,k25)
  for i in uidx
    u[i] = @muladd uprev[i] + dt*(b1*k1[i] + b2*k2[i] + b3*k3[i] + b5*k5[i] + b7*k7[i] + b8*k8[i] + b10*k10[i] + b11*k11[i] + b13*k13[i] + b14*k14[i] + b15*k15[i] + b16*k16[i] + b17*k17[i] + b18*k18[i] + b19*k19[i] + b20*k20[i] + b21*k21[i] + b22*k22[i] + b23*k23[i] + b24*k24[i] + b25*k25[i])
  end
  if integrator.opts.adaptive
    for i in uidx
      atmp[i] = (dt*(k2[i] - k24[i]) * adaptiveConst)./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
    end
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  f(t+dt,u,k)
  @pack integrator = t,dt,u
end

@inline function initialize!(integrator,cache::Feagin14ConstantCache,f=integrator.f)
  integrator.fsalfirst = f(integrator.t,integrator.uprev) # Pre-start fsal
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
end

@inline function perform_step!(integrator,cache::Feagin14ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  @unpack adaptiveConst,a0100,a0200,a0201,a0300,a0302,a0400,a0402,a0403,a0500,a0503,a0504,a0600,a0603,a0604,a0605,a0700,a0704,a0705,a0706,a0800,a0805,a0806,a0807,a0900,a0905,a0906,a0907,a0908,a1000,a1005,a1006,a1007,a1008,a1009,a1100,a1105,a1106,a1107,a1108,a1109,a1110,a1200,a1208,a1209,a1210,a1211,a1300,a1308,a1309,a1310,a1311,a1312,a1400,a1408,a1409,a1410,a1411,a1412,a1413,a1500,a1508,a1509,a1510,a1511,a1512,a1513,a1514,a1600,a1608,a1609,a1610,a1611,a1612,a1613,a1614,a1615,a1700,a1712,a1713,a1714,a1715,a1716,a1800,a1812,a1813,a1814,a1815,a1816,a1817,a1900,a1912,a1913,a1914,a1915,a1916,a1917,a1918,a2000,a2012,a2013,a2014,a2015,a2016,a2017,a2018,a2019,a2100,a2112,a2113,a2114,a2115,a2116,a2117,a2118,a2119,a2120,a2200,a2212,a2213,a2214,a2215,a2216,a2217,a2218,a2219,a2220,a2221,a2300,a2308,a2309,a2310,a2311,a2312,a2313,a2314,a2315,a2316,a2317,a2318,a2319,a2320,a2321,a2322,a2400,a2408,a2409,a2410,a2411,a2412,a2413,a2414,a2415,a2416,a2417,a2418,a2419,a2420,a2421,a2422,a2423,a2500,a2508,a2509,a2510,a2511,a2512,a2513,a2514,a2515,a2516,a2517,a2518,a2519,a2520,a2521,a2522,a2523,a2524,a2600,a2605,a2606,a2607,a2608,a2609,a2610,a2612,a2613,a2614,a2615,a2616,a2617,a2618,a2619,a2620,a2621,a2622,a2623,a2624,a2625,a2700,a2705,a2706,a2707,a2708,a2709,a2711,a2712,a2713,a2714,a2715,a2716,a2717,a2718,a2719,a2720,a2721,a2722,a2723,a2724,a2725,a2726,a2800,a2805,a2806,a2807,a2808,a2810,a2811,a2813,a2814,a2815,a2823,a2824,a2825,a2826,a2827,a2900,a2904,a2905,a2906,a2909,a2910,a2911,a2913,a2914,a2915,a2923,a2924,a2925,a2926,a2927,a2928,a3000,a3003,a3004,a3005,a3007,a3009,a3010,a3013,a3014,a3015,a3023,a3024,a3025,a3027,a3028,a3029,a3100,a3102,a3103,a3106,a3107,a3109,a3110,a3113,a3114,a3115,a3123,a3124,a3125,a3127,a3128,a3129,a3130,a3200,a3201,a3204,a3206,a3230,a3231,a3300,a3302,a3332,a3400,a3401,a3402,a3404,a3406,a3407,a3409,a3410,a3411,a3412,a3413,a3414,a3415,a3416,a3417,a3418,a3419,a3420,a3421,a3422,a3423,a3424,a3425,a3426,a3427,a3428,a3429,a3430,a3431,a3432,a3433,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19,c20,c21,c22,c23,c24,c25,c26,c27,c28,c29,c30,c31,c32,c33,c34,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18,b19,b20,b21,b22,b23,b24,b25,b26,b27,b28,b29,b30,b31,b32,b33,b34,b35 = cache
  k1  = integrator.fsalfirst
  a = dt*a0100
  k2  = f(@muladd(t + c1*dt),@muladd(uprev  + a*k1))
  k3  = f(@muladd(t + c2*dt),@muladd(uprev + dt*(a0200*k1 + a0201*k2)))
  k4  = f(@muladd(t + c3*dt),@muladd(uprev  + dt*(a0300*k1              + a0302*k3)))
  k5  = f(@muladd(t + c4*dt),@muladd(uprev  + dt*(a0400*k1              + a0402*k3 + a0403*k4)))
  k6  = f(@muladd(t + c5*dt),@muladd(uprev  + dt*(a0500*k1                           + a0503*k4 + a0504*k5)))
  k7  = f(@muladd(t + c6*dt),@muladd(uprev  + dt*(a0600*k1                           + a0603*k4 + a0604*k5 + a0605*k6)))
  k8  = f(@muladd(t + c7*dt),@muladd(uprev  + dt*(a0700*k1                                        + a0704*k5 + a0705*k6 + a0706*k7)))
  k9  = f(@muladd(t + c8*dt),@muladd(uprev  + dt*(a0800*k1                                                     + a0805*k6 + a0806*k7 + a0807*k8)))
  k10 = f(@muladd(t + c9*dt),@muladd(uprev  + dt*(a0900*k1                                                     + a0905*k6 + a0906*k7 + a0907*k8 + a0908*k9)))
  k11 = f(@muladd(t + c10*dt),@muladd(uprev + dt*(a1000*k1                                                     + a1005*k6 + a1006*k7 + a1007*k8 + a1008*k9 + a1009*k10)))
  k12 = f(@muladd(t + c11*dt),@muladd(uprev + dt*(a1100*k1                                                     + a1105*k6 + a1106*k7 + a1107*k8 + a1108*k9 + a1109*k10 + a1110*k11)))
  k13 = f(@muladd(t + c12*dt),@muladd(uprev + dt*(a1200*k1                                                                                + a1208*k9 + a1209*k10 + a1210*k11 + a1211*k12)))
  k14 = f(@muladd(t + c13*dt),@muladd(uprev + dt*(a1300*k1                                                                                + a1308*k9 + a1309*k10 + a1310*k11 + a1311*k12 + a1312*k13)))
  k15 = f(@muladd(t + c14*dt),@muladd(uprev + dt*(a1400*k1                                                                                + a1408*k9 + a1409*k10 + a1410*k11 + a1411*k12 + a1412*k13 + a1413*k14)))
  k16 = f(@muladd(t + c15*dt),@muladd(uprev + dt*(a1500*k1                                                                                + a1508*k9 + a1509*k10 + a1510*k11 + a1511*k12 + a1512*k13 + a1513*k14 + a1514*k15)))
  k17 = f(@muladd(t + c16*dt),@muladd(uprev + dt*(a1600*k1                                                                                + a1608*k9 + a1609*k10 + a1610*k11 + a1611*k12 + a1612*k13 + a1613*k14 + a1614*k15 + a1615*k16)))
  k18 = f(@muladd(t + c17*dt),@muladd(uprev + dt*(a1700*k1                                                                                + a1712*k13 + a1713*k14 + a1714*k15 + a1715*k16 + a1716*k17)))
  k19 = f(@muladd(t + c18*dt),@muladd(uprev + dt*(a1800*k1                                                                                + a1812*k13 + a1813*k14 + a1814*k15 + a1815*k16 + a1816*k17 + a1817*k18)))
  k20 = f(@muladd(t + c19*dt),@muladd(uprev + dt*(a1900*k1                                                                                + a1912*k13 + a1913*k14 + a1914*k15 + a1915*k16 + a1916*k17 + a1917*k18 + a1918*k19)))
  k21 = f(@muladd(t + c20*dt),@muladd(uprev + dt*(a2000*k1                                                                                + a2012*k13 + a2013*k14 + a2014*k15 + a2015*k16 + a2016*k17 + a2017*k18 + a2018*k19 + a2019*k20)))
  k22 = f(@muladd(t + c21*dt),@muladd(uprev + dt*(a2100*k1                                                                                + a2112*k13 + a2113*k14 + a2114*k15 + a2115*k16 + a2116*k17 + a2117*k18 + a2118*k19 + a2119*k20 + a2120*k21)))
  k23 = f(@muladd(t + c22*dt),@muladd(uprev + dt*(a2200*k1                                                                                + a2212*k13 + a2213*k14 + a2214*k15 + a2215*k16 + a2216*k17 + a2217*k18 + a2218*k19 + a2219*k20 + a2220*k21 + a2221*k22)))
  k24 = f(@muladd(t + c23*dt),@muladd(uprev + dt*(a2300*k1                                                                                + a2308*k9 + a2309*k10 + a2310*k11 + a2311*k12 + a2312*k13 + a2313*k14 + a2314*k15 + a2315*k16 + a2316*k17 + a2317*k18 + a2318*k19 + a2319*k20 + a2320*k21 + a2321*k22 + a2322*k23)))
  k25 = f(@muladd(t + c24*dt),@muladd(uprev + dt*(a2400*k1                                                                                + a2408*k9 + a2409*k10 + a2410*k11 + a2411*k12 + a2412*k13 + a2413*k14 + a2414*k15 + a2415*k16 + a2416*k17 + a2417*k18 + a2418*k19 + a2419*k20 + a2420*k21 + a2421*k22 + a2422*k23 + a2423*k24)))
  k26 = f(@muladd(t + c25*dt),@muladd(uprev + dt*(a2500*k1                                                                                + a2508*k9 + a2509*k10 + a2510*k11 + a2511*k12 + a2512*k13 + a2513*k14 + a2514*k15 + a2515*k16 + a2516*k17 + a2517*k18 + a2518*k19 + a2519*k20 + a2520*k21 + a2521*k22 + a2522*k23 + a2523*k24 + a2524*k25)))
  k27 = f(@muladd(t + c26*dt),@muladd(uprev + dt*(a2600*k1                                                     + a2605*k6 + a2606*k7 + a2607*k8 + a2608*k9 + a2609*k10 + a2610*k11               + a2612*k13 + a2613*k14 + a2614*k15 + a2615*k16 + a2616*k17 + a2617*k18 + a2618*k19 + a2619*k20 + a2620*k21 + a2621*k22 + a2622*k23 + a2623*k24 + a2624*k25 + a2625*k26)))
  k28 = f(@muladd(t + c27*dt),@muladd(uprev + dt*(a2700*k1                                                     + a2705*k6 + a2706*k7 + a2707*k8 + a2708*k9 + a2709*k10               + a2711*k12 + a2712*k13 + a2713*k14 + a2714*k15 + a2715*k16 + a2716*k17 + a2717*k18 + a2718*k19 + a2719*k20 + a2720*k21 + a2721*k22 + a2722*k23 + a2723*k24 + a2724*k25 + a2725*k26 + a2726*k27)))
  k29 = f(@muladd(t + c28*dt),@muladd(uprev + dt*(a2800*k1                                                     + a2805*k6 + a2806*k7 + a2807*k8 + a2808*k9               + a2810*k11 + a2811*k12               + a2813*k14 + a2814*k15 + a2815*k16                                                                                + a2823*k24 + a2824*k25 + a2825*k26 + a2826*k27 + a2827*k28)))
  k30 = f(@muladd(t + c29*dt),@muladd(uprev + dt*(a2900*k1                                        + a2904*k5 + a2905*k6 + a2906*k7                           + a2909*k10 + a2910*k11 + a2911*k12               + a2913*k14 + a2914*k15 + a2915*k16                                                                                + a2923*k24 + a2924*k25 + a2925*k26 + a2926*k27 + a2927*k28 + a2928*k29)))
  k31 = f(@muladd(t + c30*dt),@muladd(uprev + dt*(a3000*k1                           + a3003*k4 + a3004*k5 + a3005*k6              + a3007*k8              + a3009*k10 + a3010*k11                             + a3013*k14 + a3014*k15 + a3015*k16                                                                                + a3023*k24 + a3024*k25 + a3025*k26               + a3027*k28 + a3028*k29 + a3029*k30)))
  k32 = f(@muladd(t + c31*dt),@muladd(uprev + dt*(a3100*k1              + a3102*k3 + a3103*k4                           + a3106*k7 + a3107*k8              + a3109*k10 + a3110*k11                             + a3113*k14 + a3114*k15 + a3115*k16                                                                                + a3123*k24 + a3124*k25 + a3125*k26               + a3127*k28 + a3128*k29 + a3129*k30 + a3130*k31)))
  k33 = f(@muladd(t + c32*dt),@muladd(uprev + dt*(a3200*k1 + a3201*k2                           + a3204*k5              + a3206*k7                                                                                + a3230*k31 + a3231*k32)))
  k34 = f(@muladd(t + c33*dt),@muladd(uprev + dt*(a3300*k1              + a3302*k3                                                                                + a3332*k33)))
  k35 = f(@muladd(t + c34*dt),@muladd(uprev + dt*(a3400*k1 + a3401*k2 + a3402*k3              + a3404*k5              + a3406*k7 + a3407*k8              + a3409*k10 + a3410*k11 + a3411*k12 + a3412*k13 + a3413*k14 + a3414*k15 + a3415*k16 + a3416*k17 + a3417*k18 + a3418*k19 + a3419*k20 + a3420*k21 + a3421*k22 + a3422*k23 + a3423*k24 + a3424*k25 + a3425*k26 + a3426*k27 + a3427*k28 + a3428*k29 + a3429*k30 + a3430*k31 + a3431*k32 + a3432*k33 + a3433*k34)))
  u = @muladd uprev + dt*(b1*k1 + b2*k2 + b3*k3 + b5*k5 + b7*k7 + b8*k8 + b10*k10 + b11*k11 + b12*k12 + b14*k14 + b15*k15 + b16*k16 + b18*k18 + b19*k19 + b20*k20 + b21*k21 + b22*k22 + b23*k23 + b24*k24 + b25*k25 + b26*k26 + b27*k27 + b28*k28 + b29*k29 + b30*k30 + b31*k31 + b32*k32 + b33*k33 + b34*k34 + b35*k35)
  if integrator.opts.adaptive
    integrator.EEst = integrator.opts.internalnorm((dt*(k2 - k34) * adaptiveConst)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
  end
  k = f(t+dt,u) # For the interpolation, needs k at the updated point
  integrator.fsallast = k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  @pack integrator = t,dt,u
end

@inline function initialize!(integrator,cache::Feagin14Cache,f=integrator.f)
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  f(integrator.t,integrator.uprev,integrator.fsalfirst) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::Feagin14Cache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack adaptiveConst,a0100,a0200,a0201,a0300,a0302,a0400,a0402,a0403,a0500,a0503,a0504,a0600,a0603,a0604,a0605,a0700,a0704,a0705,a0706,a0800,a0805,a0806,a0807,a0900,a0905,a0906,a0907,a0908,a1000,a1005,a1006,a1007,a1008,a1009,a1100,a1105,a1106,a1107,a1108,a1109,a1110,a1200,a1208,a1209,a1210,a1211,a1300,a1308,a1309,a1310,a1311,a1312,a1400,a1408,a1409,a1410,a1411,a1412,a1413,a1500,a1508,a1509,a1510,a1511,a1512,a1513,a1514,a1600,a1608,a1609,a1610,a1611,a1612,a1613,a1614,a1615,a1700,a1712,a1713,a1714,a1715,a1716,a1800,a1812,a1813,a1814,a1815,a1816,a1817,a1900,a1912,a1913,a1914,a1915,a1916,a1917,a1918,a2000,a2012,a2013,a2014,a2015,a2016,a2017,a2018,a2019,a2100,a2112,a2113,a2114,a2115,a2116,a2117,a2118,a2119,a2120,a2200,a2212,a2213,a2214,a2215,a2216,a2217,a2218,a2219,a2220,a2221,a2300,a2308,a2309,a2310,a2311,a2312,a2313,a2314,a2315,a2316,a2317,a2318,a2319,a2320,a2321,a2322,a2400,a2408,a2409,a2410,a2411,a2412,a2413,a2414,a2415,a2416,a2417,a2418,a2419,a2420,a2421,a2422,a2423,a2500,a2508,a2509,a2510,a2511,a2512,a2513,a2514,a2515,a2516,a2517,a2518,a2519,a2520,a2521,a2522,a2523,a2524,a2600,a2605,a2606,a2607,a2608,a2609,a2610,a2612,a2613,a2614,a2615,a2616,a2617,a2618,a2619,a2620,a2621,a2622,a2623,a2624,a2625,a2700,a2705,a2706,a2707,a2708,a2709,a2711,a2712,a2713,a2714,a2715,a2716,a2717,a2718,a2719,a2720,a2721,a2722,a2723,a2724,a2725,a2726,a2800,a2805,a2806,a2807,a2808,a2810,a2811,a2813,a2814,a2815,a2823,a2824,a2825,a2826,a2827,a2900,a2904,a2905,a2906,a2909,a2910,a2911,a2913,a2914,a2915,a2923,a2924,a2925,a2926,a2927,a2928,a3000,a3003,a3004,a3005,a3007,a3009,a3010,a3013,a3014,a3015,a3023,a3024,a3025,a3027,a3028,a3029,a3100,a3102,a3103,a3106,a3107,a3109,a3110,a3113,a3114,a3115,a3123,a3124,a3125,a3127,a3128,a3129,a3130,a3200,a3201,a3204,a3206,a3230,a3231,a3300,a3302,a3332,a3400,a3401,a3402,a3404,a3406,a3407,a3409,a3410,a3411,a3412,a3413,a3414,a3415,a3416,a3417,a3418,a3419,a3420,a3421,a3422,a3423,a3424,a3425,a3426,a3427,a3428,a3429,a3430,a3431,a3432,a3433,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19,c20,c21,c22,c23,c24,c25,c26,c27,c28,c29,c30,c31,c32,c33,c34,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18,b19,b20,b21,b22,b23,b24,b25,b26,b27,b28,b29,b30,b31,b32,b33,b34,b35 = cache.tab
  @unpack k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,k18,k19,k20,k21,k22,k23,k24,k25,k26,k27,k28,k29,k30,k31,k32,k33,k34,k35,tmp,atmp,uprev,k = cache
  k1 = cache.fsalfirst
  f(t,uprev,k1)
  a = dt*a0100
  for i in uidx
    tmp[i] = @muladd uprev[i] + a*k1[i]
  end
  f(@muladd(t + c1*dt),tmp,k2)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0200*k1[i] + a0201*k2[i])
  end
  f(@muladd(t + c2*dt) ,tmp,k3)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0300*k1[i] + a0302*k3[i])
  end
  f(@muladd(t + c3*dt),tmp,k4)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0400*k1[i] + a0402*k3[i] + a0403*k4[i])
  end
  f(@muladd(t + c4*dt),tmp,k5)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0500*k1[i] + a0503*k4[i] + a0504*k5[i])
  end
  f(@muladd(t + c5*dt),tmp,k6)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0600*k1[i] + a0603*k4[i] + a0604*k5[i] + a0605*k6[i])
  end
  f(@muladd(t + c6*dt),tmp,k7)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0700*k1[i] + a0704*k5[i] + a0705*k6[i] + a0706*k7[i])
  end
  f(@muladd(t + c7*dt),tmp,k8)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0800*k1[i] + a0805*k6[i] + a0806*k7[i] + a0807*k8[i])
  end
  f(@muladd(t + c8*dt),tmp,k9)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a0900*k1[i] + a0905*k6[i] + a0906*k7[i] + a0907*k8[i] + a0908*k9[i])
  end
  f(@muladd(t + c9*dt),tmp,k10)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1000*k1[i] + a1005*k6[i] + a1006*k7[i] + a1007*k8[i] + a1008*k9[i] + a1009*k10[i])
  end
  f(@muladd(t + c10*dt),tmp,k11)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1100*k1[i] + a1105*k6[i] + a1106*k7[i] + a1107*k8[i] + a1108*k9[i] + a1109*k10[i] + a1110*k11[i])
  end
  f(@muladd(t + c11*dt),tmp,k12)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1200*k1[i] + a1208*k9[i] + a1209*k10[i] + a1210*k11[i] + a1211*k12[i])
  end
  f(@muladd(t + c12*dt),tmp,k13)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1300*k1[i] + a1308*k9[i] + a1309*k10[i] + a1310*k11[i] + a1311*k12[i] + a1312*k13[i])
  end
  f(@muladd(t + c13*dt),tmp,k14)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1400*k1[i] + a1408*k9[i] + a1409*k10[i] + a1410*k11[i] + a1411*k12[i] + a1412*k13[i] + a1413*k14[i])
  end
  f(@muladd(t + c14*dt),tmp,k15)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1500*k1[i] + a1508*k9[i] + a1509*k10[i] + a1510*k11[i] + a1511*k12[i] + a1512*k13[i] + a1513*k14[i] + a1514*k15[i])
  end
  f(@muladd(t + c15*dt),tmp,k16)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1600*k1[i] + a1608*k9[i] + a1609*k10[i] + a1610*k11[i] + a1611*k12[i] + a1612*k13[i] + a1613*k14[i] + a1614*k15[i] + a1615*k16[i])
  end
  f(@muladd(t + c16*dt),tmp,k17)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1700*k1[i] + a1712*k13[i] + a1713*k14[i] + a1714*k15[i] + a1715*k16[i] + a1716*k17[i])
  end
  f(@muladd(t + c17*dt),tmp,k18)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1800*k1[i] + a1812*k13[i] + a1813*k14[i] + a1814*k15[i] + a1815*k16[i] + a1816*k17[i] + a1817*k18[i])
  end
  f(@muladd(t + c18*dt),tmp,k19)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a1900*k1[i] + a1912*k13[i] + a1913*k14[i] + a1914*k15[i] + a1915*k16[i] + a1916*k17[i] + a1917*k18[i] + a1918*k19[i])
  end
  f(@muladd(t + c19*dt),tmp,k20)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a2000*k1[i] + a2012*k13[i] + a2013*k14[i] + a2014*k15[i] + a2015*k16[i] + a2016*k17[i] + a2017*k18[i] + a2018*k19[i] + a2019*k20[i])
  end
  f(@muladd(t + c20*dt),tmp,k21)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a2100*k1[i] + a2112*k13[i] + a2113*k14[i] + a2114*k15[i] + a2115*k16[i] + a2116*k17[i] + a2117*k18[i] + a2118*k19[i] + a2119*k20[i] + a2120*k21[i])
  end
  f(@muladd(t + c21*dt),tmp,k22)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a2200*k1[i] + a2212*k13[i] + a2213*k14[i] + a2214*k15[i] + a2215*k16[i] + a2216*k17[i] + a2217*k18[i] + a2218*k19[i] + a2219*k20[i] + a2220*k21[i] + a2221*k22[i])
  end
  f(@muladd(t + c22*dt),tmp,k23)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a2300*k1[i] + a2308*k9[i] + a2309*k10[i] + a2310*k11[i] + a2311*k12[i] + a2312*k13[i] + a2313*k14[i] + a2314*k15[i] + a2315*k16[i] + a2316*k17[i] + a2317*k18[i] + a2318*k19[i] + a2319*k20[i] + a2320*k21[i] + a2321*k22[i] + a2322*k23[i])
  end
  f(@muladd(t + c23*dt),tmp,k24)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a2400*k1[i] + a2408*k9[i] + a2409*k10[i] + a2410*k11[i] + a2411*k12[i] + a2412*k13[i] + a2413*k14[i] + a2414*k15[i] + a2415*k16[i] + a2416*k17[i] + a2417*k18[i] + a2418*k19[i] + a2419*k20[i] + a2420*k21[i] + a2421*k22[i] + a2422*k23[i] + a2423*k24[i])
  end
  f(@muladd(t + c24*dt),tmp,k25)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a2500*k1[i] + a2508*k9[i] + a2509*k10[i] + a2510*k11[i] + a2511*k12[i] + a2512*k13[i] + a2513*k14[i] + a2514*k15[i] + a2515*k16[i] + a2516*k17[i] + a2517*k18[i] + a2518*k19[i] + a2519*k20[i] + a2520*k21[i] + a2521*k22[i] + a2522*k23[i] + a2523*k24[i] + a2524*k25[i])
  end
  f(@muladd(t + c25*dt),tmp,k26)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a2600*k1[i] + a2605*k6[i] + a2606*k7[i] + a2607*k8[i] + a2608*k9[i] + a2609*k10[i] + a2610*k11[i] + a2612*k13[i] + a2613*k14[i] + a2614*k15[i] + a2615*k16[i] + a2616*k17[i] + a2617*k18[i] + a2618*k19[i] + a2619*k20[i] + a2620*k21[i] + a2621*k22[i] + a2622*k23[i] + a2623*k24[i] + a2624*k25[i] + a2625*k26[i])
  end
  f(@muladd(t + c26*dt),tmp,k27)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a2700*k1[i] + a2705*k6[i] + a2706*k7[i] + a2707*k8[i] + a2708*k9[i] + a2709*k10[i] + a2711*k12[i] + a2712*k13[i] + a2713*k14[i] + a2714*k15[i] + a2715*k16[i] + a2716*k17[i] + a2717*k18[i] + a2718*k19[i] + a2719*k20[i] + a2720*k21[i] + a2721*k22[i] + a2722*k23[i] + a2723*k24[i] + a2724*k25[i] + a2725*k26[i] + a2726*k27[i])
  end
  f(@muladd(t + c27*dt),tmp,k28)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a2800*k1[i] + a2805*k6[i] + a2806*k7[i] + a2807*k8[i] + a2808*k9[i] + a2810*k11[i] + a2811*k12[i] + a2813*k14[i] + a2814*k15[i] + a2815*k16[i] + a2823*k24[i] + a2824*k25[i] + a2825*k26[i] + a2826*k27[i] + a2827*k28[i])
  end
  f(@muladd(t + c28*dt),tmp,k29)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a2900*k1[i] + a2904*k5[i] + a2905*k6[i] + a2906*k7[i] + a2909*k10[i] + a2910*k11[i] + a2911*k12[i] + a2913*k14[i] + a2914*k15[i] + a2915*k16[i] + a2923*k24[i] + a2924*k25[i] + a2925*k26[i] + a2926*k27[i] + a2927*k28[i] + a2928*k29[i])
  end
  f(@muladd(t + c29*dt),tmp,k30)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a3000*k1[i] + a3003*k4[i] + a3004*k5[i] + a3005*k6[i] + a3007*k8[i] + a3009*k10[i] + a3010*k11[i] + a3013*k14[i] + a3014*k15[i] + a3015*k16[i] + a3023*k24[i] + a3024*k25[i] + a3025*k26[i] + a3027*k28[i] + a3028*k29[i] + a3029*k30[i])
  end
  f(@muladd(t + c30*dt),tmp,k31)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a3100*k1[i] + a3102*k3[i] + a3103*k4[i] + a3106*k7[i] + a3107*k8[i] + a3109*k10[i] + a3110*k11[i] + a3113*k14[i] + a3114*k15[i] + a3115*k16[i] + a3123*k24[i] + a3124*k25[i] + a3125*k26[i] + a3127*k28[i] + a3128*k29[i] + a3129*k30[i] + a3130*k31[i])
  end
  f(@muladd(t + c31*dt),tmp,k32)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a3200*k1[i] + a3201*k2[i] + a3204*k5[i] + a3206*k7[i] + a3230*k31[i] + a3231*k32[i])
  end
  f(@muladd(t + c32*dt),tmp,k33)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a3300*k1[i] + a3302*k3[i] + a3332*k33[i])
  end
  f(@muladd(t + c33*dt),tmp,k34)
  for i in uidx
    tmp[i] = @muladd uprev[i] + dt*(a3400*k1[i] + a3401*k2[i] + a3402*k3[i] + a3404*k5[i] + a3406*k7[i] + a3407*k8[i] + a3409*k10[i] + a3410*k11[i] + a3411*k12[i] + a3412*k13[i] + a3413*k14[i] + a3414*k15[i] + a3415*k16[i] + a3416*k17[i] + a3417*k18[i] + a3418*k19[i] + a3419*k20[i] + a3420*k21[i] + a3421*k22[i] + a3422*k23[i] + a3423*k24[i] + a3424*k25[i] + a3425*k26[i] + a3426*k27[i] + a3427*k28[i] + a3428*k29[i] + a3429*k30[i] + a3430*k31[i] + a3431*k32[i] + a3432*k33[i] + a3433*k34[i])
  end
  f(@muladd(t + c34*dt),tmp,k35)
  for i in uidx
    u[i] = @muladd uprev[i] + dt*(b1*k1[i] + b2*k2[i] + b3*k3[i] + b5*k5[i] + b7*k7[i] + b8*k8[i] + b10*k10[i] + b11*k11[i] + b12*k12[i] + b14*k14[i] + b15*k15[i] + b16*k16[i] + b18*k18[i] + b19*k19[i] + b20*k20[i] + b21*k21[i] + b22*k22[i] + b23*k23[i] + b24*k24[i] + b25*k25[i] + b26*k26[i] + b27*k27[i] + b28*k28[i] + b29*k29[i] + b30*k30[i] + b31*k31[i] + b32*k32[i] + b33*k33[i] + b34*k34[i] + b35*k35[i])
  end
  if integrator.opts.adaptive
    for i in uidx
      atmp[i] = (dt*(k2[i] - k34[i]) * adaptiveConst)./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
    end
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  f(t+dt,u,integrator.fsallast) # For the interpolation, needs k at the updated point
  @pack integrator = t,dt,u
end


text:
"""
Modia module with rotational component models (inspired from Modelica Standard Library).

* Developer: Hilding Elmqvist, Mogram AB, Martin Otter, DLR
* Copyright (c) 2016-2021: Hilding Elmqvist, Martin Otter
* License: MIT (expat)

"""
#module Rotational

#export Flange, Inertia, Spring, SpringDamper, EMF, IdealGear, Torque, CurrentSensor, Fixed, Damper, IdealGear_withSupport, SpeedSensor

using TinyModia
using Unitful

# Connector for 1D rotational systems
Flange = Model(
    phi = potential,
    tau = flow
)

# Flange fixed in housing at a given angle
#=
Fixed = Model(
    flange = Flange,
    phi0   = 0.0u"rad",
    equations = :[
        flange.phi = 0.0u"rad"]
)
=#
Fixed = Model(
    flange = Flange,
    equations = :[
        flange.phi = 0.0]
)

# 1D-rotational component with inertia
Inertia = Model(
    flange_a = Flange, # (info = "Left flange of shaft")
    flange_b = Flange, # (info = "Right flange of shaft")
    J = 1.0u"kg*m^2", # (0, min=0, info = "Moment of inertia") #, T = u"kg*m^2")
    phi = Var(init=0.0u"rad"),
    w = Var(init=0.0u"rad/s"),
    equations = :[
        phi = flange_a.phi
        phi = flange_b.phi
        w = der(phi)
        a = der(w)
        J * a = flange_a.tau + flange_b.tau ]
)

# Partial model for the compliant connection of two rotational 1-dim. shaft flanges
PartialCompliant = Model(
    flange_a = Flange,
    flange_b = Flange,
    partialEquations = :[
        phi_rel = flange_b.phi - flange_a.phi
        flange_b.tau = tau
        flange_a.tau = -tau ]
)

# Linear 1D rotational spring
Spring = PartialCompliant | Model(
    c = 1.0u"N*m/rad",   # (min = 0, info = "Spring constant")
    phi_rel0 = 0u"rad", # (info = "Unstretched spring angle")
    equations = :[
        tau = c * (phi_rel - phi_rel0) ]
)

# Linear 1D rotational spring with damper
SpringDamper = PartialCompliant | Model(
    c = 1.0*u"N*m/rad", # (min = 0, info = "Spring constant")
    d = 0.0u"N*m*s/rad", # (info = "Damping constant")
    phi_rel0 = 0u"rad", # (info = "Unstretched spring angle")
    equations = :[
        tau = c * (phi_rel - phi_rel0) + d * der(phi_rel) ]
)

# Electromotoric force (electric/mechanic) transformer
EMF = Model(
    k = 1.0u"N*m/A", # (info = "Transformation coefficient")
    p = Pin,
    n = Pin,
    flange = Flange,

    equations = :[
        v = p.v - n.v
        0 = p.i + n.i
        i = p.i

        phi = flange.phi
        w = der(phi)
        k * w = v
        flange.tau = -k * i ]
)

# Ideal gear
IdealGear = Model(
    flange_a = Flange, # "Left flange of shaft"
    flange_b = Flange, # "Right flange of shaft"
    ratio    = 1.0, # "Transmission ratio (flange_a.phi/flange_b.phi)"
    equations = :[
        phi_a = flange_a.phi
        phi_b = flange_b.phi
        phi_a = ratio * phi_b
        0 = ratio * flange_a.tau + flange_b.tau ]
)

# Ideal gear with support
IdealGear_withSupport = Model(
    flange_a = Flange, # "Left flange of shaft"
    flange_b = Flange, # "Right flange of shaft"
    support  = Flange, # "Support flange"
    ratio    = 1.0,  # "Transmission ratio"
    equations = :[
        phi_a = flange_a.phi - support.phi
        phi_b = flange_b.phi - support.phi
        phi_a = ratio * phi_b
        0 = ratio * flange_a.tau + flange_b.tau
        0 = flange_a.tau + flange_b.tau + support.tau ]
)

# Partial input signal acting as external torque on a flange
PartialTorque = Model(
    tau = input,
    flange = Flange
)

# Input signal acting as external torque on a flange
Torque         = PartialTorque | Model(equations = :[flange.tau = -tau])
UnitlessTorque = PartialTorque | Model(equations = :[flange.tau = -tau*u"N*m"])

#=
# Partial model for the compliant connection of two rotational 1-dim. shaft flanges where the relative angle and speed are used as preferred states
PartialCompliantWithRelativeStates = Model(
    flange_a = Flange,
    flange_b = Flange,
    init     = Map(phi_rel=0.0u"rad", w_rel=0.0u"rad/s"),
    equations = :[
        phi_rel      = flange_b.phi - flange_a.phi
        w_rel        = der(phi_rel)
        a_rel        = der(w_rel)
        flange_b.tau = tau
        flange_a.tau = -tau ]
)

# Linear 1D rotational damper
Damper = PartialCompliantWithRelativeStates | Model(
    d = 1.0u"N*m*s/rad", # (info = "Damping constant"),
    equation = :[
        tau = d * w_rel ]
)
=#
# Linear 1D rotational damper
Damper = Model(
    flange_a = Flange,
    flange_b = Flange,
    phi_rel  = Var(start=0.0u"rad"),
    d        = 1.0u"N*m*s/rad", # (info = "Damping constant"),
    equation = :[
        phi_rel      = flange_b.phi - flange_a.phi
        w_rel        = der(phi_rel)
        flange_b.tau = tau
        flange_a.tau = -tau
        tau          = d * w_rel ]
)

# Partial model to measure a single absolute flange variable
PartialAbsoluteSensor = Model(
    flange   = Flange,
    equation = :[flange.tau = 0]
)

# Ideal sensor to measure the absolute flange angular velocity
SpeedSensor         = PartialAbsoluteSensor | Model(w = output, equations = :[w = der(flange.phi)])
UnitlessSpeedSensor = PartialAbsoluteSensor | Model(w = output, equations = :[w = der(flange.phi)*u"s/rad"])

#end

ext: jl


text:
#=
cons(a, b) constructs a pair, and car(pair) and cdr(pair) returns the first and last element of that pair. For example, car(cons(3, 4)) returns 3, and cdr(cons(3, 4)) returns 4.

Given this implementation of cons:

def cons(a, b):
    def pair(f):
        return f(a, b)
    return pair

Implement car and cdr.
=#

function cons(a,b)
    function pair(f)
        return f(a,b)
    end
    return pair
end

function car(pair)
    return getfield(pair, 1)
end

function cdr(pair)
    return getfield(pair,2)
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using DiffEqBase, DiffEqOperators, LinearAlgebra, Zygote, Test
const A = rand(Float32, 300, 300)
f(du, u, p, t) = mul!(du, A, u)
f(u, p, t) = A * u
x = rand(Float32, 300)
v = rand(Float32, 300)
du = similar(x)

J = VecJacOperator(f, x)
actual_vjp = Zygote.jacobian(x -> f(x, nothing, nothing), x)[1]' * v
@test J * v ≈ actual_vjp
J = VecJacOperator(f, x; autodiff = false)
@test J * v ≈ actual_vjp


text:
struct Element{T}<:AbstractElement{T}
    𝓒::Vector{Node}
    𝓖::Vector{Node}
end
Element{T}(𝓒::Vector{Node}) where T = Element{T}(𝓒,Node[])
function Element{T}(data::Dict{Symbol,Vector{Float64}},index::Int...) where T
    𝓒 = [Node(i,data) for i in index]
    𝓖 = Node[]
    return Element{T}(𝓒,𝓖)
end

## convert
Element{T}(a::S) where {T,S<:AbstractElement} = Element{T}(a.𝓒)
function Element{T}(as::Vector{S};renumbering::Bool=false) where {T,S<:AbstractElement}
    aps = Element{T}[]
    if renumbering
        index, data = renumber(as)
        for a in as
            𝓒 = [Node(index[x.id],data) for x in a.𝓒]
            𝓖 = Node[]
            push!(aps,Element{T}(𝓒,𝓖))
        end
    else
        for a in as
            push!(aps,Element{T}(a))
        end
    end
    return aps
end

function Element{T}(a::AbstractElement,b::AbstractElement) where T
    𝓒 = a.𝓒
    𝓖 = get𝓖(a,b)
    𝓖 ≠ nothing ? Element{T}(𝓒,𝓖) : nothing
end

function Element{T}(as::Vector{A},bs::Vector{B}) where {T,A<:AbstractElement,B<:AbstractElement}
    aps = Element{T}[]
    for a in as
        for b in bs
            ap = Element{T}(a,b)
            ap ≠ nothing ? push!(aps,ap) : nothing
        end
    end
    return aps
end

function renumber(aps::Vector{T}) where T<:AbstractElement
    index = Dict{Int,Int}()
    n = 0
    for ap in aps
        for x in ap.𝓒
            I = x.id
            if ~haskey(index,I)
                n += 1
                index[I] = n
            end
        end
    end
    data_ = aps[1].𝓒[1].data
    data = Dict(:x=>zeros(n),:y=>zeros(n),:z=>zeros(n))
    for (j,i) in index
        data[:x][i] = data_[:x][j]
        data[:y][i] = data_[:y][j]
        data[:z][i] = data_[:z][j]
    end
    return index, data
end
## get𝒙
@inline get𝒙(ap::T,::Any) where T<:AbstractElement{:Poi1} = (ap.𝓒[1].x,ap.𝓒[1].y,ap.𝓒[1].z)
@inline get𝒙(ap::T,ξ::𝝃) where {T<:AbstractElement{:Seg2},𝝃<:AbstractNode} = get𝒙(ap,ξ.ξ)
@inline get𝒙(ap::T,ξ::𝝃) where {T<:AbstractElement{:Tri3},𝝃<:AbstractNode} = get𝒙(ap,ξ.ξ,ξ.η)
@inline get𝒙(ap::T,ξ::𝝃) where {T<:AbstractElement{:Quad},𝝃<:AbstractNode} = get𝒙(ap,ξ.ξ,ξ.η)

function get𝒙(ap::T,ξ::Float64) where T<:AbstractElement{:Seg2}
    x₁ = ap.𝓒[1].x
    y₁ = ap.𝓒[1].y
    z₁ = ap.𝓒[1].z
    x₂ = ap.𝓒[2].x
    y₂ = ap.𝓒[2].y
    z₂ = ap.𝓒[2].z
    N₁ = 0.5*(1-ξ)
    N₂ = 0.5*(1+ξ)
    return (x₁*N₁+x₂*N₂,y₁*N₁+y₂*N₂,z₁*N₁+z₂*N₂)
end
function get𝒙(ap::T,ξ::Float64,η::Float64) where T<:AbstractElement{:Tri3}
    x₁ = ap.𝓒[1].x
    y₁ = ap.𝓒[1].y
    z₁ = ap.𝓒[1].z
    x₂ = ap.𝓒[2].x
    y₂ = ap.𝓒[2].y
    z₂ = ap.𝓒[2].z
    x₃ = ap.𝓒[3].x
    y₃ = ap.𝓒[3].y
    z₃ = ap.𝓒[3].z
    N₁ = ξ
    N₂ = η
    N₃ = 1.0-ξ-η
    return (x₁*N₁+x₂*N₂+x₃*N₃,y₁*N₁+y₂*N₂+y₃*N₃,z₁*N₁+z₂*N₂+z₃*N₃)
end

function get𝒙(ap::T,ξ::Float64,η::Float64) where T<:AbstractElement{:Quad}
    x₁ = ap.𝓒[1].x
    y₁ = ap.𝓒[1].y
    z₁ = ap.𝓒[1].z
    x₂ = ap.𝓒[2].x
    y₂ = ap.𝓒[2].y
    z₂ = ap.𝓒[2].z
    x₃ = ap.𝓒[3].x
    y₃ = ap.𝓒[3].y
    z₃ = ap.𝓒[3].z
    x₄ = ap.𝓒[4].x
    y₄ = ap.𝓒[4].y
    z₄ = ap.𝓒[4].z
    N₁,N₂,N₃,N₄ = get𝝭(ap,ξ,η)
    return (x₁*N₁+x₂*N₂+x₃*N₃+x₄*N₄,y₁*N₁+y₂*N₂+y₃*N₃+y₄*N₄,z₁*N₁+z₂*N₂+z₃*N₃+z₄*N₄)
end
## get∇𝒙
function get𝑱(ap::T,ξ::𝝃) where {T<:AbstractElement{:Quad},𝝃<:AbstractNode}
    x₁ = ap.𝓒[1].x
    x₂ = ap.𝓒[2].x
    x₃ = ap.𝓒[3].x
    x₄ = ap.𝓒[4].x
    y₁ = ap.𝓒[1].y
    y₂ = ap.𝓒[2].y
    y₃ = ap.𝓒[3].y
    y₄ = ap.𝓒[4].y
    ∂N₁∂ξ,∂N₂∂ξ,∂N₃∂ξ,∂N₄∂ξ = get∂𝝭∂ξ(ap,ξ)
    ∂N₁∂η,∂N₂∂η,∂N₃∂η,∂N₄∂η = get∂𝝭∂η(ap,ξ)
    J₁₁ = ∂N₁∂ξ*x₁ + ∂N₂∂ξ*x₂ + ∂N₃∂ξ*x₃ + ∂N₄∂ξ*x₄
    J₁₂ = ∂N₁∂η*x₁ + ∂N₂∂η*x₂ + ∂N₃∂η*x₃ + ∂N₄∂η*x₄
    J₂₁ = ∂N₁∂ξ*y₁ + ∂N₂∂ξ*y₂ + ∂N₃∂ξ*y₃ + ∂N₄∂ξ*y₄
    J₂₂ = ∂N₁∂η*y₁ + ∂N₂∂η*y₂ + ∂N₃∂η*y₃ + ∂N₄∂η*y₄
    return J₁₁,J₂₁,J₁₂,J₂₂
end
## get𝐽
@inline get𝐽(  ::T,::Any) where T<:AbstractElement{:Poi1} = 1.0
@inline get𝐽(ap::T,::Any) where T<:AbstractElement{:Seg2} = 0.5*get𝐿(ap)
@inline get𝐽(ap::T,::Any) where T<:AbstractElement{:Tri3} = 2.0*get𝐴(ap)
@inline function get𝐽(ap::T,ξ::𝝃) where {T<:AbstractElement{:Quad},𝝃<:AbstractNode}
    J₁₁,J₂₁,J₁₂,J₂₂ = get𝑱(ap,ξ)
    return J₁₁*J₂₂-J₂₁*J₁₂
end
## get𝑤
@inline get𝑤(  ::T,::Any) where T<:AbstractElement{:Poi1} = 1.0
@inline get𝑤(ap::T,ξ::𝝃) where {T<:AbstractElement{:Seg2},𝝃<:AbstractNode} = 0.5*get𝐿(ap)*ξ.w
@inline get𝑤(ap::T,ξ::𝝃) where {T<:AbstractElement{:Tri3},𝝃<:AbstractNode} = get𝐴(ap)*ξ.w
@inline get𝑤(ap::T,ξ::𝝃) where {T<:AbstractElement{:Quad},𝝃<:AbstractNode} = get𝐽(ap,ξ)*ξ.w

## get𝐿 get𝐴 get𝑉
@inline function get𝐿(ap::T) where T<:AbstractElement{:Seg2}
    x₁ = ap.𝓒[1].x
    y₁ = ap.𝓒[1].y
    z₁ = ap.𝓒[1].z
    x₂ = ap.𝓒[2].x
    y₂ = ap.𝓒[2].y
    z₂ = ap.𝓒[2].z
    return ((x₂-x₁)^2+(y₂-y₁)^2+(z₂-z₁)^2)^0.5
end
function get𝐴(ap::T) where T<:AbstractElement{:Tri3}
    x₁ = ap.𝓒[1].x
    y₁ = ap.𝓒[1].y
    z₁ = ap.𝓒[1].z
    x₂ = ap.𝓒[2].x
    y₂ = ap.𝓒[2].y
    z₂ = ap.𝓒[2].z
    x₃ = ap.𝓒[3].x
    y₃ = ap.𝓒[3].y
    z₃ = ap.𝓒[3].z
    𝐴₁ = 0.5*(y₁*z₂+y₂*z₃+y₃*z₁-y₂*z₁-y₃*z₂-y₁*z₃)
    𝐴₂ = 0.5*(z₁*x₂+z₂*x₃+z₃*x₁-z₂*x₁-z₃*x₂-z₁*x₃)
    𝐴₃ = 0.5*(x₁*y₂+x₂*y₃+x₃*y₁-x₂*y₁-x₃*y₂-x₁*y₃)
    return (𝐴₁^2 + 𝐴₂^2 + 𝐴₃^2)^0.5
end
## get𝒏
@inline get𝒏(ap::T) where T<:AbstractElement{:Poi1} = 1.0
@inline function get𝒏(ap::T) where T<:AbstractElement{:Seg2}
    x₁ = ap.𝓒[1].x
    y₁ = ap.𝓒[1].y
    x₂ = ap.𝓒[2].x
    y₂ = ap.𝓒[2].y
    𝐿 = get𝐿(ap)
    return (y₂-y₁)/𝐿,(x₁-x₂)/𝐿,(x₂-x₁)/𝐿,(y₂-y₁)/𝐿
end

@inline get𝒏(ap::T,ξ::𝝃) where {T<:AbstractElement{:Seg2},𝝃<:AbstractNode} = get𝒏(ap,ξ.ξ)
@inline get𝒏(ap::T,ξ::𝝃) where {T<:AbstractElement{:Tri3},𝝃<:AbstractNode} = get𝒏(ap,ξ.ξ,ξ.η)

function get𝒏(ap::T,ξ::Float64) where T<:AbstractElement{:Seg2}
    n₁ = 0.0
    n₁ += ξ == -1.0 ? -1.0 : 0.0
    n₁ += ξ ==  1.0 ?  1.0 : 0.0
    return n₁
end
function get𝒏(ap::T,ξ::Float64,η::Float64) where T<:AbstractElement{:Tri3}
    n₁ = 0
    n₂ = 0
    x₁ = ap.𝓒[1].x
    y₁ = ap.𝓒[1].y
    x₂ = ap.𝓒[2].x
    y₂ = ap.𝓒[2].y
    x₃ = ap.𝓒[3].x
    y₃ = ap.𝓒[3].y
    γ = 1.0-ξ-η
    n₁ += ξ == 0.0 ? y₃-y₂ : 0.0
    n₁ += η == 0.0 ? y₁-y₃ : 0.0
    n₁ += γ == 0.0 ? y₂-y₁ : 0.0
    n₂ += ξ == 0.0 ? x₂-x₃ : 0.0
    n₂ += η == 0.0 ? x₃-x₁ : 0.0
    n₂ += γ == 0.0 ? x₁-x₂ : 0.0
    return n₁,n₂
end
## set𝒏!
function set𝒏!(aps::Vector{T}) where T<:AbstractElement
    for ap in aps
        set𝒏!(ap)
    end
end

function set𝒏!(ap::T) where T<:AbstractElement{:Seg2}
    𝓖 = ap.𝓖
    for ξ in 𝓖
        ξ.n₁ = get𝒏(ap,ξ)
    end
end

function set𝒏!(ap::T) where T<:AbstractElement{:Tri3}
    𝓖 = ap.𝓖
    for ξ in 𝓖
        ξ.n₁, ξ.n₂ = get𝒏(ap,ξ)
    end
end

## shape functions
# ------------- Poi1 ---------------
@inline get𝝭(::Element{:Poi1},::Any) = 1.0
# ------------- Seg2 ---------------
@inline get𝝭(ap::Element{:Seg2},ξ::𝝃) where 𝝃<:AbstractNode = get𝝭(ap,ξ.ξ)
@inline get𝝭(ap::Element{:Seg2},ξ::Float64) = (0.5*(1-ξ),0.5*(1+ξ))
@inline function get∂𝝭∂x(ap::Element{:Seg2},::Any)
    𝐿 = get𝐿(ap)
    return (-1.0/𝐿,1.0/𝐿)
end
@inline get∂𝝭∂y(ap::Element{:Seg2},::Any) = (0.0,0.0)
@inline get∂𝝭∂z(ap::Element{:Seg2},::Any) = (0.0,0.0)
@inline get∇𝝭(ap::Element{:Seg2},ξ::𝝃) where 𝝃<:AbstractNode = (get𝝭(ap,ξ),get∂𝝭∂x(ap,ξ),(0.0,0.0),(0.0,0.0))
@inline function get∂𝝭∂𝑛(ap::Element{:Seg2},ξ::𝝃) where 𝝃<:AbstractNode
    n₁ = get𝒏(ap,ξ)
    𝐿 = get𝐿(ap)
    return (-n₁/𝐿,n₁/𝐿)
end
@inline get∇𝑛𝝭(ap::Element{:Seg2},ξ::𝝃) where 𝝃<:AbstractNode = (get𝝭(ap,ξ),get∂𝝭∂𝑛(ap,ξ))
# ------------- Tri3 ---------------
@inline get𝝭(ap::Element{:Tri3},ξ::𝝃) where 𝝃<:AbstractNode = (ξ.ξ,ξ.η,1.0-ξ.ξ-ξ.η)
@inline function get∂𝝭∂x(ap::Element{:Tri3},ξ::𝝃) where 𝝃<:AbstractNode
    y₁ = ap.𝓒[1].y
    y₂ = ap.𝓒[2].y
    y₃ = ap.𝓒[3].y
    𝐴 = get𝐴(ap)
    return (y₂-y₃)/2.0/𝐴,(y₃-y₁)/2.0/𝐴,(y₁-y₂)/2.0/𝐴
end

@inline function get∂𝝭∂y(ap::Element{:Tri3},ξ::𝝃) where 𝝃<:AbstractNode
    x₁ = ap.𝓒[1].x
    x₂ = ap.𝓒[2].x
    x₃ = ap.𝓒[3].x
    𝐴 = get𝐴(ap)
    return (x₃-x₂)/2.0/𝐴,(x₁-x₃)/2.0/𝐴,(x₂-x₁)/2.0/𝐴
end
@inline get∂𝝭∂z(ap::Element{:Tri3},::Any) = (0.0,0.0,0.0)
@inline get∇𝝭(ap::Element{:Tri3},ξ::𝝃) where 𝝃<:AbstractNode = (get𝝭(ap,ξ),get∂𝝭∂x(ap,ξ),get∂𝝭∂y(ap,ξ),(0.0,0.0,0.0))

# ------------- Quad ---------------
@inline get𝝭(ap::Element{:Quad},ξ::𝝃) where 𝝃<:AbstractNode = get𝝭(ap,ξ.ξ,ξ.η)
@inline get∂𝝭∂ξ(ap::Element{:Quad},ξ::𝝃) where 𝝃<:AbstractNode = get∂𝝭∂ξ(ap,ξ.η)
@inline get∂𝝭∂η(ap::Element{:Quad},ξ::𝝃) where 𝝃<:AbstractNode = get∂𝝭∂η(ap,ξ.ξ)

function get𝝭(ap::Element{:Quad},ξ::Float64,η::Float64)
    N₁ = 0.25*(1.0-ξ)*(1.0-η)
    N₂ = 0.25*(1.0+ξ)*(1.0-η)
    N₃ = 0.25*(1.0+ξ)*(1.0+η)
    N₄ = 0.25*(1.0-ξ)*(1.0+η)
    return (N₁,N₂,N₃,N₄)
end
function get∂𝝭∂ξ(ap::Element{:Quad},η::Float64)
    ∂N₁∂ξ = - 0.25*(1-η)
    ∂N₂∂ξ =   0.25*(1-η)
    ∂N₃∂ξ =   0.25*(1+η)
    ∂N₄∂ξ = - 0.25*(1+η)
    return (∂N₁∂ξ,∂N₂∂ξ,∂N₃∂ξ,∂N₄∂ξ)
end
function get∂𝝭∂η(ap::Element{:Quad},ξ::Float64)
    ∂N₁∂η = - 0.25*(1-ξ)
    ∂N₂∂η = - 0.25*(1+ξ)
    ∂N₃∂η =   0.25*(1+ξ)
    ∂N₄∂η =   0.25*(1-ξ)
    return (∂N₁∂η,∂N₂∂η,∂N₃∂η,∂N₄∂η)
end
function get∂𝝭∂x∂𝝭∂y(ap::Element{:Quad},ξ::𝝃) where 𝝃<:AbstractNode
    x₁ = ap.𝓒[1].x
    x₂ = ap.𝓒[2].x
    x₃ = ap.𝓒[3].x
    x₄ = ap.𝓒[4].x
    y₁ = ap.𝓒[1].y
    y₂ = ap.𝓒[2].y
    y₃ = ap.𝓒[3].y
    y₄ = ap.𝓒[4].y
    ∂N₁∂ξ,∂N₂∂ξ,∂N₃∂ξ,∂N₄∂ξ = get∂𝝭∂ξ(ap,ξ)
    ∂N₁∂η,∂N₂∂η,∂N₃∂η,∂N₄∂η = get∂𝝭∂η(ap,ξ)
    ∂x∂ξ = ∂N₁∂ξ*x₁ + ∂N₂∂ξ*x₂ + ∂N₃∂ξ*x₃ + ∂N₄∂ξ*x₄
    ∂x∂η = ∂N₁∂η*x₁ + ∂N₂∂η*x₂ + ∂N₃∂η*x₃ + ∂N₄∂η*x₄
    ∂y∂ξ = ∂N₁∂ξ*y₁ + ∂N₂∂ξ*y₂ + ∂N₃∂ξ*y₃ + ∂N₄∂ξ*y₄
    ∂y∂η = ∂N₁∂η*y₁ + ∂N₂∂η*y₂ + ∂N₃∂η*y₃ + ∂N₄∂η*y₄
    detJ = ∂x∂ξ*∂y∂η - ∂x∂η*∂y∂ξ
    ∂ξ∂x =   ∂y∂η/detJ
    ∂η∂x = - ∂y∂ξ/detJ
    ∂ξ∂y = - ∂x∂η/detJ
    ∂η∂y =   ∂x∂ξ/detJ
    ∂N₁∂x = ∂N₁∂ξ*∂ξ∂x + ∂N₁∂η*∂η∂x
    ∂N₂∂x = ∂N₂∂ξ*∂ξ∂x + ∂N₂∂η*∂η∂x
    ∂N₃∂x = ∂N₃∂ξ*∂ξ∂x + ∂N₃∂η*∂η∂x
    ∂N₄∂x = ∂N₄∂ξ*∂ξ∂x + ∂N₄∂η*∂η∂x
    ∂N₁∂y = ∂N₁∂ξ*∂ξ∂y + ∂N₁∂η*∂η∂y
    ∂N₂∂y = ∂N₂∂ξ*∂ξ∂y + ∂N₂∂η*∂η∂y
    ∂N₃∂y = ∂N₃∂ξ*∂ξ∂y + ∂N₃∂η*∂η∂y
    ∂N₄∂y = ∂N₄∂ξ*∂ξ∂y + ∂N₄∂η*∂η∂y
    return (∂N₁∂x,∂N₂∂x,∂N₃∂x,∂N₄∂x),(∂N₁∂y,∂N₂∂y,∂N₃∂y,∂N₄∂y)
end
@inline get∇𝝭(ap::Element{:Quad},ξ::𝝃) where 𝝃<:AbstractNode = get𝝭(ap,ξ),get∂𝝭∂x∂𝝭∂y(ap,ξ)...,(0.0,0.0,0.0,0.0)

## ⊆,∩
function issubset(a::T,b::S) where {T<:AbstractElement{:Poi1},S<:AbstractElement{:Seg2}}
    i = findfirst(x->x==a.𝓒[1],b.𝓒)
    return i ≠ nothing && i ≤ 2
end

@inline intersect(a::T,b::T) where T<:AbstractElement = a.𝓒 == b.𝓒 ? a : nothing
@inline function intersect(a::T,b::S) where {T<:AbstractElement{:Seg2},S<:AbstractElement{:Poi1}}
    i = findfirst(x->x==b.𝓒[1],a.𝓒)
    return i ≠ nothing && i ≤ 2 ? a : nothing
end
@inline function intersect(a::T,b::S) where {T<:AbstractElement{:Tri3},S<:AbstractElement{:Seg2}}
    i = findfirst(x->x==b.𝓒[1],a.𝓒)
    j = findfirst(x->x==b.𝓒[2],a.𝓒)
    return i ≠ nothing && j ≠ nothing && i ≤ 3 && j ≤ 3 ? a : nothing
end
function intersect(as::Vector{T},bs::Vector{S}) where {T<:AbstractElement,S<:AbstractElement}
    aps = T[]
    for b in bs
        for a in as
            ap = a∩b
            ap ≠ nothing ? push!(aps,ap) : nothing
        end
    end
    return aps
end

## getnₚ,getnₑ
getnₚ(a::T) where T<:AbstractElement = length(a.𝓒[1].data[:x])
getnₚ(as::Vector{T}) where T<:AbstractElement = getnₚ(as[1])


text:
function FieldODEState(arg0::RealFieldElement, arg1::Vector{RealFieldElement})
    return FieldODEState((RealFieldElement, Vector{RealFieldElement}), arg0, arg1)
end

function FieldODEState(arg0::RealFieldElement, arg1::Vector{RealFieldElement}, arg2::Vector{Vector{RealFieldElement}})
    return FieldODEState((RealFieldElement, Vector{RealFieldElement}, Vector{Vector{RealFieldElement}}), arg0, arg1, arg2)
end

function get_complete_state(obj::FieldODEState)
    return jcall(obj, "getCompleteState", Vector{RealFieldElement}, ())
end

function get_complete_state_dimension(obj::FieldODEState)
    return jcall(obj, "getCompleteStateDimension", jint, ())
end

function get_number_of_secondary_states(obj::FieldODEState)
    return jcall(obj, "getNumberOfSecondaryStates", jint, ())
end

function get_primary_state(obj::FieldODEState)
    return jcall(obj, "getPrimaryState", Vector{RealFieldElement}, ())
end

function get_primary_state_dimension(obj::FieldODEState)
    return jcall(obj, "getPrimaryStateDimension", jint, ())
end

function get_secondary_state(obj::FieldODEState, arg0::jint)
    return jcall(obj, "getSecondaryState", Vector{RealFieldElement}, (jint,), arg0)
end

function get_secondary_state_dimension(obj::FieldODEState, arg0::jint)
    return jcall(obj, "getSecondaryStateDimension", jint, (jint,), arg0)
end

function get_time(obj::FieldODEState)
    return jcall(obj, "getTime", RealFieldElement, ())
end


text:
#__precompile__(true)

## License follows that of https://github.com/eiscor/eiscor
## This codebase is derived in large part from software at https://people.cs.kuleuven.be/~raf.vandebril/
## which (presumably) has found its way into the eiscor code base and released under the MIT license

module AMVW
using LinearAlgebra
import Base: adjoint

#using IntervalArithmetic
#Base.eps{T}(::Type{Interval{T}}) = eps(T)

# package code goes here

include("types.jl")
include("utils.jl")
include("transformations.jl")
include("factorization.jl")
include("diagonal-block.jl")
include("bulge.jl")
include("deflation.jl")
include("AMVW_algorithm.jl")

include("diagnostics.jl")

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Gingko_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Gingko")
JLLWrappers.@generate_main_file("Gingko", UUID("ebc12d66-1ae6-5e16-bc1b-c04a92e425de"))
end  # module Gingko_jll


text:
# sample from a double well in 1D using HMC
using Printf
using Random
using LinearAlgebra
using ForwardDiff
using QuadGK
using BasicMD
using Statistics

include("potentials.jl")

β = 5.0;
x₀ = [-1.0];
seed = 100;
Δt = 1e-1;
n_iters = 10^4; # number of samples

V = x->DoubleWell(x);

# second moment and energy obsrevables
f₁ = x-> x[1]^2;
f₂ = x-> V(x[1]);
observables = (f₁, f₂);

sampler = RWM(V, β, Δt);

Random.seed!(100);
Random.seed!(100);
Xvals,_ = sample_trajectory(x₀, sampler, options=MDOptions(n_iters=n_iters,n_save_iters=n_save_iters));
f₁_estimate = mean(f₁.(Xvals));
f₂_estimate = mean(f₂.(Xvals));
@printf("f₁ estimate with %d samples: %g\n",length(Xvals), f₁_estimate);
@printf("f₂ estimate with %d samples: %g\n", length(Xvals), f₂_estimate);
Random.seed!(100);
observable_samples = sample_observables(x₀, sampler, observables, options=MDOptions(n_iters=n_iters,n_save_iters=n_save_iters));
@printf("f₁ estimate with %d samples: %g\n",length(observable_samples), mean(observable_samples[1,:]));
@printf("f₂ estimate with %d samples: %g\n", length(observable_samples), mean(observable_samples[2,:]));

max_stars_repo_name: gideonsimpson/JuBasicMD


text:
function F(x)
    # 1 + sum(x.^2)
    y = one(eltype(x))
    for _x in x
        y += _x^2
    end
    return y
end

function ∇F!(g, x)
    g .= 2 .* x
end

function H!(g, x)
    g .= 0
    for i in eachindex(x)
        g[i,i] = 2
    end
end


text:
using DifferentialEquations
using LSODA
using JLD2
using ArgParse

parser = ArgParseSettings(allow_ambiguous_opts=false)
@add_arg_table! parser begin
    "output"
        help = "output JLD2"
        required = true
    "pset"
        help = "parameter set: snic, sn, hopf, sl"
        required = true
    "form"
        help = "noise form: add or mul"
        required = true
    "initstate"
        help = "initial state: hi or lo"
        required = true
    "signals"
        help = "signal values"
        required = true
    "--cells"
        help = "number of cells to simulate"
        default = 500
        arg_type = Int
    "--noises"
        help = "noise values"
        default = "0.5:0.5:1.5"
    "--stochtime"
        help = "stochastic simulation endpoint"
        arg_type = Float64
        default = 200.0
end
args = parse_args(parser)

function model_mmi2_add(du, u, p, t)
    sR, a1, a2, b1, b2, kR, kr, K1, K2, koff, sr, sR0, σr, σR, τ = p

    sig = sR

    r, R, c1A, c1B, c2 = u[1], u[2], u[3], u[4], u[5]

    du[1] = dr = sr - kr * r + koff * (c1A+c1B) - K1 * koff * 2R * r +
    koff * 2c2 - K2 * koff * (c1A+c1B) * r +
    (c1A+c1B) * kR * a1 + 2c2 * kR * a2 + u[6]

    du[2] = dR = sig + sR0 - kR * R + koff * (c1A+c1B) - K1 * koff * 2R * r +
    (c1A+c1B) * kr * b1 + u[7]

    du[3] = dc1A = K1 * koff * R * r - koff * c1A +
    koff * c2 - K2 * koff * c1A * r +
    c2 * 1kr * b2 - c1A * kR * a1 - c1A * kr * b1

    du[4] = dc1B = K1 * koff * R * r - koff * c1B +
    koff * c2 - K2 * koff * c1B * r +
    c2 * 1kr * b2 - c1B * kR * a1 - c1B * kr * b1

    du[5] = dc2 = K2 * koff * (c1A+c1B) * r - koff * 2c2 -
    c2 * 2kr * b2 - c2 * kR * a2

    du[6] = - u[6]/τ
    du[7] = - u[7]/τ
end

function stoch_add(du,u,p,t)
    sR, a1, a2, b1, b2, kR, kr, K1, K2, koff, sr, sR0, σr, σR, τ = p
    du[1] = 0
    du[2] = 0
    du[3] = 0
    du[4] = 0
    du[5] = 0
    du[6] = sqrt((2.0*σr^2.0/τ))
    du[7] = sqrt((2.0*σR^2.0/τ))
end

function model_mmi2_mul(du, u, p, t)
    sR, a1, a2, b1, b2, kR, kr, K1, K2, koff, sr, sR0, σr, σR, τ = p

    sig = sR

    r, R, c1A, c1B, c2 = u[1], u[2], u[3], u[4], u[5]

    du[1] = dr = sr - kr * r + koff * (c1A+c1B) - K1 * koff * 2R * r +
    koff * 2c2 - K2 * koff * (c1A+c1B) * r +
    (c1A+c1B) * kR * a1 + 2c2 * kR * a2 + u[6]

    du[2] = dR = sig + sR0 - kR * R + koff * (c1A+c1B) - K1 * koff * 2R * r +
    (c1A+c1B) * kr * b1 + u[7]

    du[3] = dc1A = K1 * koff * R * r - koff * c1A +
    koff * c2 - K2 * koff * c1A * r +
    c2 * 1kr * b2 - c1A * kR * a1 - c1A * kr * b1 + u[8]

    du[4] = dc1B = K1 * koff * R * r - koff * c1B +
    koff * c2 - K2 * koff * c1B * r +
    c2 * 1kr * b2 - c1B * kR * a1 - c1B * kr * b1 + u[9]

    du[5] = dc2 = K2 * koff * (c1A+c1B) * r - koff * 2c2 -
    c2 * 2kr * b2 - c2 * kR * a2 + u[10]

    du[6] = - u[6]/τ
    du[7] = - u[7]/τ
    du[8] = - u[8]/τ
    du[9] = - u[9]/τ
    du[10] = - u[10]/τ
end

function stoch_mul(du,u,p,t)
    sR, a1, a2, b1, b2, kR, kr, K1, K2, koff, sr, sR0, σr, σR, τ = p
    du[1] = 0
    du[2] = 0
    du[3] = 0
    du[4] = 0
    du[5] = 0
    du[6] = sqrt((2.0*σR^2.0/τ))*u[1]
    du[7] = sqrt((2.0*σR^2.0/τ))*u[2]
    du[8] = sqrt((2.0*σR^2.0/τ))*u[3]
    du[9] = sqrt((2.0*σR^2.0/τ))*u[4]
    du[10] = sqrt((2.0*σR^2.0/τ))*u[5]
end

p_snic = [ # SNIC, like Figure 5B
    0.3, #1 sR
    1.0, #2 a1
    12.0, #3 a2
    1.0, #4 b1
    4.0, #5 b2
    1.0, #6 kR
    0.25, #7 kr
    1000.0, #8 K1
    1000.0, #9 K2
    100.0, #10 koff
    1.0, #11 sr
    5.7, #12 sR0    total 6
    0.00, #13 σr    replaced
    0.4, #14 σR     replaced
    0.1 #15 τ
   ];

p_hopf = [ # Diverging Hopf, like Figure 3B
   0.5, #1 sR
   1.0, #2 a1
   12.0, #3 a2
   1.0, #4 b1
   7.0, #5 b2
   1.0, #6 kR
   0.25, #7 kr
   1000.0, #8 K1
   1000.0, #9 K2
   100.0, #10 koff
   1.0, #11 sr
   3.1, #12 sR0    total 3.6
   0.00, #13 σr    replaced
   0.4, #14 σR     replaced
   0.1 #15 τ
  ];

p_sn = [ # SN, like Figure 3F
    3.0, #1 sR
    1.0, #2 a1
    4.0, #3 a2
    0.5, #4 b1
    0.1, #5 b2
    1.0, #6 kR
    2.0, #7 kr
    1000.0, #8 K1
    1000.0, #9 K2
    100.0, #10 koff
    1.0, #11 sr
    0.0, #12 sR0    total 3
    0.00, #13 σr    replaced
    0.4, #14 σR     replaced
    0.1 #15 τ
   ];

p_sl = [ # Saddle-loop, like Figure 5G
   1.0, #1 sR
   1.0, #2 a1
   12.0, #3 a2
   1.0, #4 b1
   3.0, #5 b2
   1.0, #6 kR
   0.25, #7 kr
   1000.0, #8 K1
   1000.0, #9 K2
   100.0, #10 koff
   1.0, #11 sr
   6.9, #12 sR0    total 7.9
   0.00, #13 σr    replaced
   0.4, #14 σR     replaced
   0.1 #15 τ
  ];

function filter_sols(sols)
    nvars, ntps, nsols = size(sols)
    good_sols = []
    for i in 1:nsols
        if isassigned(sols[i][1,:], ntps) == true
            push!(good_sols, i)
        end
    end
    println("Obtained $(length(good_sols))/$nsols converged solutions.")
    sols[:, :, good_sols]
end

function sim_pop_ic(model, model_stoch, u0, pars; popsize=1000, selectR::String="lo", saveat=0.1)
    selectR == "lo" || selectR == "hi" || error("invalid selectR")
    prob_ode = ODEProblem(model, u0, (0., 200.), pars);
    prob_sde = SDEProblem(model, model_stoch, u0, tspan, pars);
    r0s = range(0.001, 10.0, length=20)
    R0s = range(0.001, 10.0, length=20)
    function prob_func_cells_init(prob,i,repeat)
        m, n = i ÷ length(r0s), i % length(r0s)
        if n != 0
            r0, R0 = r0s[m+1],  R0s[n]
        else
            r0, R0 = r0s[m], R0s[end]
        end
        remake(prob, u0=[r0; R0; u0[3:end]])
    end
    prob_cells_init = EnsembleProblem(prob_ode, prob_func=prob_func_cells_init)
    @time sim_init = solve(prob_cells_init, lsoda(), saveat=saveat, dt=0.01,
        reltol=1e-6, abstol=1e-6, trajectories=length(r0s)*length(R0s));
    sol_init = sim_init[:, 400:end, :] # The selection of state is sensitive to the start index
    sol_reduced = dropdims(sum(sol_init[2:5, :, :]; dims=1), dims=1)
    if selectR == "lo"
        sol_reduced[isnan.(sol_reduced)] .= Inf
        idsel = argmin(sol_reduced)
    else
        sol_reduced[isnan.(sol_reduced)] .= -Inf
        idsel = argmax(sol_reduced)
    end
    u0_s = sol_init[:, idsel[1], idsel[2]]
    println("Starting stochastic simulation with init cond ", u0_s)
    function prob_func_cells(prob,i,repeat)
        remake(prob,u0=u0_s)
    end
    prob_cells = EnsembleProblem(prob_sde, prob_func=prob_func_cells)
    @time sim = solve(prob_cells, STrapezoid(), saveat=saveat, maxiter=1e12, dt=0.0001, trajectories=popsize);
    filter_sols(sim)
end

u0s = Dict("add" => [0.8, 0.1, 0., 0., 0., 0., 0.], "mul" => [0.8, 0.1, 0., 0., 0., 0., 0., 0., 0., 0.])
models = Dict("add" => model_mmi2_add, "mul" => model_mmi2_mul)
stochs = Dict("add" => stoch_add, "mul" => stoch_mul)
psets = Dict("snic" => p_snic, "sn" => p_sn, "hopf" => p_hopf, "sl" => p_sl)

model = models[args["form"]]
stoch = stochs[args["form"]]
u0 = u0s[args["form"]]
p = psets[args["pset"]]
tspan = (0., args["stochtime"])
noises = eval(Meta.parse(args["noises"]))
signals = eval(Meta.parse(args["signals"]))

all_results = map(noises) do noise
    results = map(signals) do signal
        signal => sim_pop_ic(model, stoch, u0, [signal; p[2:end-3]; noise*0.25; noise*1.0; p[end]]; popsize=args["cells"], selectR=args["initstate"])
    end
    noise => Dict(results)
end
save(args["output"], "results", Dict(all_results))


text:
# run locally
using TextSearch

include(joinpath(ENV["JULIA_PROJECT"], "src/io.jl"))

const VALID = Set(['😄', '😍', '😘', '😊', '😌', '🥺', '😡', '😒', '😭', '😢',  '❤', '💔', '👌', '👏', '🤔'])

function create_emoji_dataset(filename)
    trainfile = replace(filename, ".txt" => "") * "-emojis-train.txt"
    testfile = replace(filename, ".txt" => "") * "-emojis-test.txt"
    @info filename => [trainfile, testfile]
    #isfile(outfile) && return
    s = 0
    train = open(trainfile, "w")
    test = open(testfile, "w")
    F = [train, test]
    for line in eachline(filename)
        emolist = emojis(line)
        e = intersect(VALID, emolist)
        if length(e) == 1
            e = first(e)
            line = replace(line, e => " _emo ")
            println(rand(F), "__label__", e, ' ', line)
            s += 1
        end
    end

    close(train)
    close(test)
end

for filename in ARGS
    create_emoji_dataset(filename)
end

max_stars_repo_name: sadit/RegionalSpanish

max_issues_repo_name: sadit/RegionalSpanish

max_forks_repo_name: sadit/RegionalSpanish


text:
using MPI
using ClimateMachine
using Logging
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using LinearAlgebra
using Printf
using Test

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

const output = parse(Bool, lowercase(get(ENV, "JULIA_CLIMA_OUTPUT", "false")))

include("advection_diffusion_model.jl")

struct Pseudo1D{n1, n2, α, β, μ, δ} <: AdvectionDiffusionProblem end

function init_velocity_diffusion!(
    ::Pseudo1D{n1, n2, α, β},
    aux::Vars,
    geom::LocalGeometry,
) where {n1, n2, α, β}
    # Direction of flow is n1 (resp n2) with magnitude α
    aux.advection.u = hcat(α * n1, α * n2)

    # diffusion of strength β in the n1 and n2 directions
    aux.diffusion.D = hcat(β * n1 * n1', β * n2 * n2')
end

function initial_condition!(
    ::Pseudo1D{n1, n2, α, β, μ, δ},
    state,
    aux,
    localgeo,
    t,
) where {n1, n2, α, β, μ, δ}
    ξn1 = dot(n1, localgeo.coord)
    ξn2 = dot(n2, localgeo.coord)
    ρ1 = exp(-(ξn1 - μ - α * t)^2 / (4 * β * (δ + t))) / sqrt(1 + t / δ)
    ρ2 = exp(-(ξn2 - μ - α * t)^2 / (4 * β * (δ + t))) / sqrt(1 + t / δ)
    state.ρ = (ρ1, ρ2)
end

Dirichlet_data!(P::Pseudo1D, x...) = initial_condition!(P, x...)

function Neumann_data!(
    ::Pseudo1D{n1, n2, α, β, μ, δ},
    ∇state,
    aux,
    x,
    t,
) where {n1, n2, α, β, μ, δ}
    ξn1 = dot(n1, x)
    ξn2 = dot(n2, x)
    ∇ρ1 =
        -(
            2n1 * (ξn1 - μ - α * t) / (4 * β * (δ + t)) *
            exp(-(ξn1 - μ - α * t)^2 / (4 * β * (δ + t))) / sqrt(1 + t / δ)
        )
    ∇ρ2 =
        -(
            2n2 * (ξn2 - μ - α * t) / (4 * β * (δ + t)) *
            exp(-(ξn2 - μ - α * t)^2 / (4 * β * (δ + t))) / sqrt(1 + t / δ)
        )
    ∇state.ρ = hcat(∇ρ1, ∇ρ2)
end

function test_run(mpicomm, dim, polynomialorders, level, ArrayType, FT)

    n_hd =
        dim == 2 ? SVector{3, FT}(1, 0, 0) :
        SVector{3, FT}(1 / sqrt(2), 1 / sqrt(2), 0)

    n_vd = dim == 2 ? SVector{3, FT}(0, 1, 0) : SVector{3, FT}(0, 0, 1)

    α = FT(1)
    β = FT(1 // 100)
    μ = FT(-1 // 2)
    δ = FT(1 // 10)

    # Grid/topology information
    base_num_elem = 4
    Ne = 2^(level - 1) * base_num_elem
    brickrange = ntuple(j -> range(FT(-1); length = Ne + 1, stop = 1), dim)
    periodicity = ntuple(j -> false, dim)
    bc = ntuple(j -> (1, 2), dim)

    topl = StackedBrickTopology(
        mpicomm,
        brickrange;
        periodicity = periodicity,
        boundary = bc,
    )

    dt = (α / 4) / (Ne * maximum(polynomialorders)^2)
    timeend = 1
    @info "time step" dt

    @info @sprintf """Test parameters:
    ArrayType                   = %s
    FloatType                   = %s
    Dimension                   = %s
    Horizontal polynomial order = %s
    Vertical polynomial order   = %s
      """ ArrayType FT dim polynomialorders[1] polynomialorders[end]

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = polynomialorders,
    )

    # Model being tested
    model = AdvectionDiffusion{dim}(
        Pseudo1D{n_hd, n_vd, α, β, μ, δ}(),
        num_equations = 2,
    )

    # Main DG discretization
    dg = DGModel(
        model,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
        direction = EveryDirection(),
    )

    # Initialize all relevant state arrays and create solvers
    Q = init_ode_state(dg, FT(0))

    eng0 = norm(Q, dims = (1, 3))
    @info @sprintf """Starting
    norm(Q₀) = %.16e""" eng0[1]

    solver = LSRK54CarpenterKennedy(dg, Q; dt = dt, t0 = 0)
    solve!(Q, solver; timeend = timeend)

    # Reference solution
    engf = norm(Q, dims = (1, 3))
    Q_ref = init_ode_state(dg, FT(timeend))

    engfe = norm(Q_ref, dims = (1, 3))
    errf = norm(Q_ref .- Q, dims = (1, 3))

    metrics = @. (engf, engf / eng0, engf - eng0, errf, errf / engfe)

    @info @sprintf """Finished
    Horizontal field:
      norm(Q)                 = %.16e
      norm(Q) / norm(Q₀)      = %.16e
      norm(Q) - norm(Q₀)      = %.16e
      norm(Q - Qe)            = %.16e
      norm(Q - Qe) / norm(Qe) = %.16e
    Vertical field:
      norm(Q)                 = %.16e
      norm(Q) / norm(Q₀)      = %.16e
      norm(Q) - norm(Q₀)      = %.16e
      norm(Q - Qe)            = %.16e
      norm(Q - Qe) / norm(Qe) = %.16e
      """ first.(metrics)... last.(metrics)...

    return errf
end

"""
    main()

Run this test problem
"""
function main()

    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()
    mpicomm = MPI.COMM_WORLD

    # Dictionary keys: dim, level, polynomial order, FT, and direction
    expected_result = Dict()

    # Dim 2, degree 4 in the horizontal, Float64
    expected_result[2, 1, 4, Float64, HorizontalDirection] = 0.0467837436192571
    expected_result[2, 2, 4, Float64, HorizontalDirection] =
        0.004066556782723549
    expected_result[2, 3, 4, Float64, HorizontalDirection] =
        5.3144336694234015e-5
    expected_result[2, 4, 4, Float64, HorizontalDirection] =
        3.978000110046181e-7

    # Dim 2, degree 2 in the vertical, Float64
    expected_result[2, 1, 2, Float64, VerticalDirection] = 0.15362016594121006
    expected_result[2, 2, 2, Float64, VerticalDirection] = 0.04935353328794371
    expected_result[2, 3, 2, Float64, VerticalDirection] = 0.015530511948609192
    expected_result[2, 4, 2, Float64, VerticalDirection] = 0.0006275095484456197

    # Dim 2, degree 2 in the horizontal, Float64
    expected_result[2, 1, 2, Float64, HorizontalDirection] = 0.15362016594121003
    expected_result[2, 2, 2, Float64, HorizontalDirection] = 0.04935353328794369
    expected_result[2, 3, 2, Float64, HorizontalDirection] =
        0.015530511948609204
    expected_result[2, 4, 2, Float64, HorizontalDirection] =
        0.0006275095484455967

    # Dim 2, degree 4 in the vertical, Float64
    expected_result[2, 1, 4, Float64, VerticalDirection] = 0.04678374361925714
    expected_result[2, 2, 4, Float64, VerticalDirection] = 0.0040665567827235
    expected_result[2, 3, 4, Float64, VerticalDirection] = 5.3144336694109365e-5
    expected_result[2, 4, 4, Float64, VerticalDirection] = 3.978000109805811e-7

    # Dim 3, degree 4 in the horizontal, Float64
    expected_result[3, 1, 4, Float64, HorizontalDirection] =
        0.017475667486259432
    expected_result[3, 2, 4, Float64, HorizontalDirection] =
        0.0012502148161420109
    expected_result[3, 3, 4, Float64, HorizontalDirection] =
        6.999081063570052e-5
    expected_result[3, 4, 4, Float64, HorizontalDirection] =
        2.8724182090419642e-6

    # Dim 3, degree 2 in the vertical, Float64
    expected_result[3, 1, 2, Float64, VerticalDirection] = 0.2172517221280645
    expected_result[3, 2, 2, Float64, VerticalDirection] = 0.06979643612684193
    expected_result[3, 3, 2, Float64, VerticalDirection] = 0.02196346062832051
    expected_result[3, 4, 2, Float64, VerticalDirection] = 0.0008874325139302493

    # Dim 3, degree 2 in the horizontal, Float64
    expected_result[3, 1, 2, Float64, HorizontalDirection] = 0.10343354980172516
    expected_result[3, 2, 2, Float64, HorizontalDirection] = 0.03415137756593495
    expected_result[3, 3, 2, Float64, HorizontalDirection] =
        0.0035959803480493553
    expected_result[3, 4, 2, Float64, HorizontalDirection] =
        0.0002714157844893719

    # Dim 3, degree 4 in the vertical, Float64
    expected_result[3, 1, 4, Float64, VerticalDirection] = 0.06616220472493903
    expected_result[3, 2, 4, Float64, VerticalDirection] = 0.005750979754288175
    expected_result[3, 3, 4, Float64, VerticalDirection] = 7.515744171591452e-5
    expected_result[3, 4, 4, Float64, VerticalDirection] = 5.625741705890895e-7

    # Dim 2, degree 4 in the horizontal, Float32
    expected_result[2, 1, 4, Float32, HorizontalDirection] = 0.046783954f0
    expected_result[2, 2, 4, Float32, HorizontalDirection] = 0.004066328f0
    expected_result[2, 3, 4, Float32, HorizontalDirection] = 5.327546f-5

    # Dim 2, degree 2 in the vertical, Float32
    expected_result[2, 1, 2, Float32, VerticalDirection] = 0.15362015f0
    expected_result[2, 2, 2, Float32, VerticalDirection] = 0.04935346f0
    expected_result[2, 3, 2, Float32, VerticalDirection] = 0.015530386f0

    # Dim 2, degree 2 in the horizontal, Float32
    expected_result[2, 1, 2, Float32, HorizontalDirection] = 0.1536202f0
    expected_result[2, 2, 2, Float32, HorizontalDirection] = 0.04935346f0
    expected_result[2, 3, 2, Float32, HorizontalDirection] = 0.015530357f0

    # Dim 2, degree 4 in the vertical, Float32
    expected_result[2, 1, 4, Float32, VerticalDirection] = 0.04678398f0
    expected_result[2, 2, 4, Float32, VerticalDirection] = 0.0040662177f0
    expected_result[2, 3, 4, Float32, VerticalDirection] = 5.3401447f-5

    # Dim 3, degree 4 in the horizontal, Float32
    expected_result[3, 1, 4, Float32, HorizontalDirection] = 0.01747554f0
    expected_result[3, 2, 4, Float32, HorizontalDirection] = 0.0012502924f0
    expected_result[3, 3, 4, Float32, HorizontalDirection] = 7.00218f-5

    # Dim 3, degree 2 in the vertical, Float32
    expected_result[3, 1, 2, Float32, VerticalDirection] = 0.21725166f0
    expected_result[3, 2, 2, Float32, VerticalDirection] = 0.06979626f0
    expected_result[3, 3, 2, Float32, VerticalDirection] = 0.021963252f0

    # Dim 3, degree 2 in the horizontal, Float32
    expected_result[3, 1, 2, Float32, HorizontalDirection] = 0.10343349f0
    expected_result[3, 2, 2, Float32, HorizontalDirection] = 0.034151305f0
    expected_result[3, 3, 2, Float32, HorizontalDirection] = 0.0035958516f0

    # Dim 3, degree 4 in the vertical, Float32
    expected_result[3, 1, 4, Float32, VerticalDirection] = 0.06616244f0
    expected_result[3, 2, 4, Float32, VerticalDirection] = 0.005750495f0
    expected_result[3, 3, 4, Float32, VerticalDirection] = 7.538217f-5

    @testset "Variable degree DG: advection diffusion model" begin
        for FT in (Float32, Float64)
            numlevels =
                integration_testing ||
                ClimateMachine.Settings.integration_testing ?
                (FT == Float64 ? 4 : 3) : 1
            for dim in 2:3
                for polynomialorders in ((4, 2), (2, 4))
                    result = Dict()
                    for level in 1:numlevels
                        result[level] = test_run(
                            mpicomm,
                            dim,
                            polynomialorders,
                            level,
                            ArrayType,
                            FT,
                        )
                        horiz_poly = polynomialorders[1]
                        vert_poly = polynomialorders[2]
                        @test result[level][1] ≈ FT(expected_result[
                            dim,
                            level,
                            horiz_poly,
                            FT,
                            HorizontalDirection,
                        ])
                        @test result[level][2] ≈ FT(expected_result[
                            dim,
                            level,
                            vert_poly,
                            FT,
                            VerticalDirection,
                        ])
                    end
                    @info begin
                        msg = ""
                        for l in 1:(numlevels - 1)
                            rate = @. log2(result[l]) - log2(result[l + 1])
                            msg *= @sprintf(
                                "\n  rates for level %d Horizontal = %e",
                                l,
                                rate[1]
                            )
                            msg *= @sprintf(", Vertical = %e\n", rate[2])
                        end
                        msg
                    end
                end
            end
        end
    end
end

main()

max_forks_repo_name: jm-c/CLIMA


text:
using Pkg
@info "Activating project environment"
Pkg.activate(ENV["JULIA_PROJECT"])
@info "installing dependencies"
@time Pkg.add([
    "PackageCompiler",
    "MLJ",
    "MLJScikitLearnInterface",
    "MLJDecisionTreeInterface",
    "MLJClusteringInterface",
    "MLJGLMInterface",
    "MLJLIBSVMInterface",
    "MLJMultivariateStatsInterface",
    "MLJXGBoostInterface",
    "Pandas",
    "DataFrames",
    "CSV",
    "PyCall",
])

using PackageCompiler,
    MLJ,
    MLJScikitLearnInterface,
    MLJDecisionTreeInterface,
    MLJClusteringInterface,
    MLJGLMInterface,
    MLJLIBSVMInterface,
    MLJMultivariateStatsInterface,
    MLJScikitLearnInterface,
    MLJXGBoostInterface,
    Pandas,
    DataFrames,
    CSV,
    PyCall
@info "creating system image at $(ENV["JULIA_SYS_IMAGE"])"
@time create_sysimage(
    [
        :MLJ,
        :MLJScikitLearnInterface,
        :MLJDecisionTreeInterface,
        :MLJClusteringInterface,
        :MLJGLMInterface,
        :MLJLIBSVMInterface,
        :MLJMultivariateStatsInterface,
        :MLJXGBoostInterface,
        :Pandas,
        :DataFrames,
        :CSV,
        :PyCall,
    ];
    sysimage_path = ENV["JULIA_SYS_IMAGE"],
)
@info "Julia environment setup complete! \U0001F37E"

max_stars_repo_name: andreakropp/datarobot-user-models

max_issues_repo_name: andreakropp/datarobot-user-models

max_forks_repo_name: andreakropp/datarobot-user-models


text:
module RAISDKResultsWrapper

import RAI
import JSON3
using DataStructures: OrderedSet

export output, tuples, relations
export RelationIterator, TuplesIterator, PhysicalRelation

include("tuples.jl")
include("physical-relations.jl")
include("trie-iterator.jl")
include("relation-iterator.jl")

end


text:
struct VehicleState
    posG::VecSE2{Float64} # global
    posF::Frenet # lane-relative frame
    v::Float64
end

VehicleState() = VehicleState(VecSE2(), NULL_FRENET, NaN)
VehicleState(posG::VecSE2{Float64}, v::Float64) = VehicleState(posG, NULL_FRENET, v)
VehicleState(posG::VecSE2{Float64}, roadway::Roadway, v::Float64) = VehicleState(posG, Frenet(posG, roadway), v)
VehicleState(posG::VecSE2{Float64}, lane::Lane, roadway::Roadway, v::Float64) = VehicleState(posG, Frenet(posG, lane, roadway), v)
VehicleState(posF::Frenet, roadway::Roadway, v::Float64) = VehicleState(get_posG(posF, roadway), posF, v)

Base.show(io::IO, s::VehicleState) = print(io, "VehicleState(", s.posG, ", ", s.posF, ", ", @sprintf("%.3f", s.v), ")")

function Base.write(io::IO, ::MIME"text/plain", s::VehicleState)
    @printf(io, "%.16e %.16e %.16e", s.posG.x, s.posG.y, s.posG.θ)
    @printf(io, " %d %.16e %d %d", s.posF.roadind.ind.i, s.posF.roadind.ind.t, s.posF.roadind.tag.segment, s.posF.roadind.tag.lane)
    @printf(io, " %.16e %.16e %.16e", s.posF.s, s.posF.t, s.posF.ϕ)
    @printf(io, " %.16e", s.v)
end
function Base.read(io::IO, ::MIME"text/plain", ::Type{VehicleState})
    tokens = split(strip(readline(io)), ' ')
    i = 0
    posG = VecSE2(parse(Float64, tokens[i+=1]), parse(Float64, tokens[i+=1]), parse(Float64, tokens[i+=1]))
    roadind = RoadIndex(CurveIndex(parse(Int, tokens[i+=1]), parse(Float64, tokens[i+=1])),
                        LaneTag(parse(Int, tokens[i+=1]), parse(Int, tokens[i+=1])))
    posF = Frenet(roadind, parse(Float64, tokens[i+=1]), parse(Float64, tokens[i+=1]), parse(Float64, tokens[i+=1]))
    v = parse(Float64, tokens[i+=1])
    return VehicleState(posG, posF, v)
end

function Vec.lerp(a::VehicleState, b::VehicleState, t::Float64, roadway::Roadway)
    posG = lerp(a.posG, b.posG, t)
    v = lerp(a.v, b.v, t)
    VehicleState(posG, roadway, v)
end

get_vel_s(s::VehicleState) = s.v * cos(s.posF.ϕ) # velocity along the lane
get_vel_t(s::VehicleState) = s.v * sin(s.posF.ϕ) # velocity ⟂ to lane

function move_along(vehstate::VehicleState, roadway::Roadway, Δs::Float64;
    ϕ₂::Float64=vehstate.posF.ϕ, t₂::Float64=vehstate.posF.t, v₂::Float64=vehstate.v
    )

    roadind = move_along(vehstate.posF.roadind, roadway, Δs)
    try
        footpoint = roadway[roadind]
    catch
        println(roadind)
    end
    footpoint = roadway[roadind]
    posG = convert(VecE2, footpoint.pos) + polar(t₂, footpoint.pos.θ + π/2)
    posG = VecSE2(posG.x, posG.y, footpoint.pos.θ + ϕ₂)
    VehicleState(posG, roadway, v₂)
end

const Vehicle = Entity{VehicleState,VehicleDef,Int}

Base.show(io::IO, v::Vehicle) = print(io, "Vehicle(", v.id, ", ", v.state, ", ", v.def, ")")

get_center(veh::Vehicle) = veh.state.posG
get_footpoint(veh::Vehicle) = veh.state.posG + polar(veh.state.posF.t, veh.state.posG.θ-veh.state.posF.ϕ-π/2)
get_front(veh::Vehicle) = veh.state.posG + polar(veh.def.length/2, veh.state.posG.θ)
get_rear(veh::Vehicle) = veh.state.posG - polar(veh.def.length/2, veh.state.posG.θ)

function get_lane_width(veh::Vehicle, roadway::Roadway)
    # println("veh.state.posF.roadind.tag")
    # println(veh.state.posF.roadind.tag)
    lane = roadway[veh.state.posF.roadind.tag]

    if n_lanes_left(lane, roadway) > 0
        footpoint = get_footpoint(veh)
        lane_left = roadway[LaneTag(lane.tag.segment, lane.tag.lane + 1)]
        return -proj(footpoint, lane_left, roadway).curveproj.t
    else
        return lane.width
    end
end

function get_markerdist_left(veh::Vehicle, roadway::Roadway)
    t = veh.state.posF.t
    lane_width = get_lane_width(veh, roadway)
    return lane_width/2 - t
end
function get_markerdist_right(veh::Vehicle, roadway::Roadway)
    t = veh.state.posF.t
    lane_width = get_lane_width(veh, roadway)
    return lane_width/2 + t
end

function propagate(veh::Vehicle, action::LaneFollowingAccel, roadway::Roadway, ΔT::Float64)

    a_lon = action.a

    ds = veh.state.v

    ΔT² = ΔT*ΔT
    Δs = ds*ΔT + 0.5*a_lon*ΔT²

    v₂ = ds + a_lon*ΔT

    roadind = move_along(veh.state.posF.roadind, roadway, Δs)
    posG = roadway[roadind].pos
    posF = Frenet(roadind, roadway, t = veh.state.posF.t, ϕ = veh.state.posF.ϕ)
    VehicleState(posG, posF, v₂)
end

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
using Adapt

using Oceananigans.BoundaryConditions
using Oceananigans.BoundaryConditions: AbstractBoundaryConditionClassification

import Base: show
import Oceananigans.BoundaryConditions: bc_str, print_condition

struct CubedSphereExchange <: AbstractBoundaryConditionClassification end

const CubedSphereExchangeBC = BoundaryCondition{<:CubedSphereExchange}

bc_str(::CubedSphereExchangeBC) ="CubedSphereExchange"

CubedSphereExchangeBoundaryCondition(val; kwargs...) = BoundaryCondition(CubedSphereExchange, val; kwargs...)

struct CubedSphereExchangeInformation{F, S}
    from_face :: F
      to_face :: F
    from_side :: S
      to_side :: S
end

CubedSphereExchangeInformation(; from_face, to_face, from_side, to_side) =
    CubedSphereExchangeInformation(from_face, to_face, from_side, to_side)

Base.show(io::IO, ex::CubedSphereExchangeInformation) =
    print(io, "CubedSphereExchangeInformation: (from: face $(ex.from_face) $(ex.from_side) side, to: face $(ex.to_face) $(ex.to_side) side)")

print_condition(info::CubedSphereExchangeInformation) =
    "(from: face $(info.from_face) $(info.from_side) side, to: face $(info.to_face) $(info.to_side) side)"

function inject_cubed_sphere_exchange_boundary_conditions(field_bcs, face_number, face_connectivity)

    west_exchange_info = CubedSphereExchangeInformation(
        from_face = face_number,
        from_side = :west,
          to_face = face_connectivity[face_number].west.face,
          to_side = face_connectivity[face_number].west.side
    )

    east_exchange_info = CubedSphereExchangeInformation(
        from_face = face_number,
        from_side = :east,
          to_face = face_connectivity[face_number].east.face,
          to_side = face_connectivity[face_number].east.side
    )

    south_exchange_info = CubedSphereExchangeInformation(
        from_face = face_number,
        from_side = :south,
          to_face = face_connectivity[face_number].south.face,
          to_side = face_connectivity[face_number].south.side
    )

    north_exchange_info = CubedSphereExchangeInformation(
        from_face = face_number,
        from_side = :north,
          to_face = face_connectivity[face_number].north.face,
          to_side = face_connectivity[face_number].north.side
    )

    west_exchange_bc = CubedSphereExchangeBoundaryCondition(west_exchange_info)
    east_exchange_bc = CubedSphereExchangeBoundaryCondition(east_exchange_info)
    south_exchange_bc = CubedSphereExchangeBoundaryCondition(south_exchange_info)
    north_exchange_bc = CubedSphereExchangeBoundaryCondition(north_exchange_info)

    return FieldBoundaryConditions(west_exchange_bc,
                                   east_exchange_bc,
                                   south_exchange_bc,
                                   north_exchange_bc,
                                   field_bcs.bottom,
                                   field_bcs.top,
                                   field_bcs.immersed)
end

Adapt.adapt_structure(to, ::CubedSphereExchangeInformation) = nothing
Adapt.adapt_structure(to, ::CubedSphereExchangeBC) = nothing