text:
@inline function initialize!(integrator,cache::TanYam7ConstantCache,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::TanYam7ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  @unpack c1,c2,c3,c4,c5,c6,c7,a21,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,a91,a93,a94,a95,a96,a97,a98,a101,a103,a104,a105,a106,a107,a108,b1,b4,b5,b6,b7,b8,b9,bhat1,bhat4,bhat5,bhat6,bhat7,bhat8,bhat10 = cache
  k1 = integrator.fsalfirst
  a = dt*a21
  k2 = f(@muladd(t + c1*dt),@muladd(uprev+a*k1))
  k3 = f(@muladd(t + c2*dt),@muladd(uprev+dt*(a31*k1+a32*k2)))
  k4 = f(@muladd(t + c3*dt),@muladd(uprev+dt*(a41*k1       +a43*k3)))
  k5 = f(@muladd(t + c4*dt),@muladd(uprev+dt*(a51*k1       +a53*k3+a54*k4)))
  k6 = f(@muladd(t + c5*dt),@muladd(uprev+dt*(a61*k1       +a63*k3+a64*k4+a65*k5)))
  k7 = f(@muladd(t + c6*dt),@muladd(uprev+dt*(a71*k1       +a73*k3+a74*k4+a75*k5+a76*k6)))
  k8 = f(@muladd(t + c7*dt),@muladd(uprev+dt*(a81*k1       +a83*k3+a84*k4+a85*k5+a86*k6+a87*k7)))
  k9 = f(t+dt,@muladd(uprev+dt*(a91*k1       +a93*k3+a94*k4+a95*k5+a96*k6+a97*k7+a98*k8)))
  k10= f(t+dt,@muladd(uprev+dt*(a101*k1      +a103*k3+a104*k4+a105*k5+a106*k6+a107*k7+a108*k8)))
  u = @muladd uprev + dt*(k1*b1+k4*b4+k5*b5+k6*b6+k7*b7+k8*b8+k9*b9)
  if integrator.opts.adaptive
    utilde = @muladd uprev + dt*(k1*bhat1+k4*bhat4+k5*bhat5+k6*bhat6+k7*bhat7+k8*bhat8+k10*bhat10)
    integrator.EEst = integrator.opts.internalnorm( ((utilde-u)./@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::TanYam7Cache,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::TanYam7Cache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack fsalfirst,k2,k3,k4,k5,k6,k7,k8,k9,k10,utilde,tmp,atmp,k = cache
  @unpack c1,c2,c3,c4,c5,c6,c7,a21,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,a91,a93,a94,a95,a96,a97,a98,a101,a103,a104,a105,a106,a107,a108,b1,b4,b5,b6,b7,b8,b9,bhat1,bhat4,bhat5,bhat6,bhat7,bhat8,bhat10 = cache.tab
  k1 = fsalfirst
  f(t,uprev,k1)
  a = dt*a21
  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*(a31*k1[i]+a32*k2[i])
  end
  f(@muladd(t + c2*dt),tmp,k3)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a41*k1[i]+a43*k3[i])
  end
  f(@muladd(t + c3*dt),tmp,k4)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a51*k1[i]+a53*k3[i]+a54*k4[i])
  end
  f(@muladd(t + c4*dt),tmp,k5)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a61*k1[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
  end
  f(@muladd(t + c5*dt),tmp,k6)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a71*k1[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
  end
  f(@muladd(t + c6*dt),tmp,k7)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a81*k1[i]+a83*k3[i]+a84*k4[i]+a85*k5[i]+a86*k6[i]+a87*k7[i])
  end
  f(@muladd(t + c7*dt),tmp,k8)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a91*k1[i]+a93*k3[i]+a94*k4[i]+a95*k5[i]+a96*k6[i]+a97*k7[i]+a98*k8[i])
  end
  f(t+dt,tmp,k9)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a101*k1[i]+a103*k3[i]+a104*k4[i]+a105*k5[i]+a106*k6[i]+a107*k7[i]+a108*k8[i])
  end
  f(t+dt,tmp,k10)
  for i in uidx
    u[i] = uprev[i] + dt*(k1[i]*b1+k4[i]*b4+k5[i]*b5+k6[i]*b6+k7[i]*b7+k8[i]*b8+k9[i]*b9)
  end
  if integrator.opts.adaptive
    for i in uidx
      utilde[i] = uprev[i] + dt*(k1[i]*bhat1+k4[i]*bhat4+k5[i]*bhat5+k6[i]*bhat6+k7[i]*bhat7+k8[i]*bhat8+k10[i]*bhat10)
      atmp[i] = ((utilde[i]-u[i])./(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::DP8ConstantCache,f=integrator.f)
  integrator.kshortsize = 7
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.fsalfirst = f(integrator.t,integrator.uprev) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::DP8ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  @unpack c7,c8,c9,c10,c11,c6,c5,c4,c3,c2,b1,b6,b7,b8,b9,b10,b11,b12,bhh1,bhh2,bhh3,er1,er6,er7,er8,er9,er10,er11,er12,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211 = cache
  k1 = integrator.fsalfirst
  a = dt*a0201
  k2 = f(@muladd(t + c2*dt),uprev+a*k1)
  k3 = f(@muladd(t + c3*dt),@muladd(uprev+dt*(a0301*k1+a0302*k2)))
  k4 = f(@muladd(t + c4*dt),@muladd(uprev+dt*(a0401*k1       +a0403*k3)))
  k5 = f(@muladd(t + c5*dt),@muladd(uprev+dt*(a0501*k1       +a0503*k3+a0504*k4)))
  k6 = f(@muladd(t + c6*dt),@muladd(uprev+dt*(a0601*k1                +a0604*k4+a0605*k5)))
  k7 = f(@muladd(t + c7*dt),@muladd(uprev+dt*(a0701*k1                +a0704*k4+a0705*k5+a0706*k6)))
  k8 = f(@muladd(t + c8*dt),@muladd(uprev+dt*(a0801*k1                +a0804*k4+a0805*k5+a0806*k6+a0807*k7)))
  k9 = f(@muladd(t + c9*dt),@muladd(uprev+dt*(a0901*k1                +a0904*k4+a0905*k5+a0906*k6+a0907*k7+a0908*k8)))
  k10 =f(@muladd(t + c10*dt),@muladd(uprev+dt*(a1001*k1                +a1004*k4+a1005*k5+a1006*k6+a1007*k7+a1008*k8+a1009*k9)))
  k11= f(@muladd(t + c11*dt),@muladd(uprev+dt*(a1101*k1                +a1104*k4+a1105*k5+a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10)))
  k12= f(t+dt,@muladd(uprev+dt*(a1201*k1                +a1204*k4+a1205*k5+a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11)))
  kupdate= @muladd b1*k1+b6*k6+b7*k7+b8*k8+b9*k9+b10*k10+b11*k11+b12*k12
  update = dt*kupdate
  u = uprev + update
  if integrator.opts.adaptive
    err5 = integrator.opts.internalnorm(dt*(@muladd(k1*er1 + k6*er6 + k7*er7 + k8*er8 + k9*er9 + k10*er10 + k11*er11 + k12*er12))./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)) # Order 5
    err3 = integrator.opts.internalnorm(@muladd(update - dt*(bhh1*k1 + bhh2*k9 + bhh3*k12))./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)) # Order 3
    err52 = err5*err5
    if err5 ≈ 0 && err3 ≈ 0
      integrator.EEst = zero(typeof(integrator.EEst))
    else
      integrator.EEst = err52./sqrt(err52 + 0.01*err3*err3)
    end
  end
  k13 = f(t+dt,u)
  integrator.fsallast = k13
  if integrator.opts.calck
    @unpack c14,c15,c16,a1401,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1613,a1614,a1615 = cache
    @unpack d401,d406,d407,d408,d409,d410,d411,d412,d413,d414,d415,d416,d501,d506,d507,d508,d509,d510,d511,d512,d513,d514,d515,d516,d601,d606,d607,d608,d609,d610,d611,d612,d613,d614,d615,d616,d701,d706,d707,d708,d709,d710,d711,d712,d713,d714,d715,d716 = cache
    k14 = f(@muladd(t + c14*dt),@muladd(uprev+dt*(a1401*k1         +a1407*k7+a1408*k8+a1409*k9+a1410*k10+a1411*k11+a1412*k12+a1413*k13)))
    k15 = f(@muladd(t + c15*dt),@muladd(uprev+dt*(a1501*k1+a1506*k6+a1507*k7+a1508*k8                   +a1511*k11+a1512*k12+a1513*k13+a1514*k14)))
    k16 = f(@muladd(t + c16*dt),@muladd(uprev+dt*(a1601*k1+a1606*k6+a1607*k7+a1608*k8+a1609*k9                              +a1613*k13+a1614*k14+a1615*k15)))
    udiff = kupdate
    integrator.k[1] = udiff
    bspl = k1 - udiff
    integrator.k[2] = bspl
    integrator.k[3] = udiff - k13 - bspl
    integrator.k[4] = @muladd(d401*k1+d406*k6+d407*k7+d408*k8+d409*k9+d410*k10+d411*k11+d412*k12+d413*k13+d414*k14+d415*k15+d416*k16)
    integrator.k[5] = @muladd(d501*k1+d506*k6+d507*k7+d508*k8+d509*k9+d510*k10+d511*k11+d512*k12+d513*k13+d514*k14+d515*k15+d516*k16)
    integrator.k[6] = @muladd(d601*k1+d606*k6+d607*k7+d608*k8+d609*k9+d610*k10+d611*k11+d612*k12+d613*k13+d614*k14+d615*k15+d616*k16)
    integrator.k[7] = @muladd(d701*k1+d706*k6+d707*k7+d708*k8+d709*k9+d710*k10+d711*k11+d712*k12+d713*k13+d714*k14+d715*k15+d716*k16)
  end
  @pack integrator = t,dt,u,k
end

@inline function initialize!(integrator,cache::DP8Cache,f=integrator.f)
  integrator.kshortsize = 7
  integrator.k = [cache.udiff,cache.bspl,cache.dense_tmp3,cache.dense_tmp4,cache.dense_tmp5,cache.dense_tmp6,cache.dense_tmp7]
  integrator.fsalfirst = cache.k1
  integrator.fsallast = cache.k13
  f(integrator.t,integrator.uprev,integrator.fsalfirst) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::DP8Cache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack c7,c8,c9,c10,c11,c6,c5,c4,c3,c2,b1,b6,b7,b8,b9,b10,b11,b12,bhh1,bhh2,bhh3,er1,er6,er7,er8,er9,er10,er11,er12,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211 = cache.tab
  @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,update,udiff,bspl,dense_tmp3,dense_tmp4,dense_tmp5,dense_tmp6,dense_tmp7,kupdate,utilde,tmp,atmp,atmp2 = cache
  f(t,uprev,k1)
  a = dt*a0201
  for i in uidx
    tmp[i] = @muladd uprev[i]+a*k1[i]
  end
  f(@muladd(t + c2*dt),tmp,k2)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a0301*k1[i]+a0302*k2[i])
  end
  f(@muladd(t + c3*dt),tmp,k3)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a0401*k1[i]+a0403*k3[i])
  end
  f(@muladd(t + c4*dt),tmp,k4)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a0501*k1[i]+a0503*k3[i]+a0504*k4[i])
  end
  f(@muladd(t + c5*dt),tmp,k5)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a0601*k1[i]+a0604*k4[i]+a0605*k5[i])
  end
  f(@muladd(t + c6*dt),tmp,k6)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a0701*k1[i]+a0704*k4[i]+a0705*k5[i]+a0706*k6[i])
  end
  f(@muladd(t + c7*dt),tmp,k7)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a0801*k1[i]+a0804*k4[i]+a0805*k5[i]+a0806*k6[i]+a0807*k7[i])
  end
  f(@muladd(t + c8*dt),tmp,k8)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a0901*k1[i]+a0904*k4[i]+a0905*k5[i]+a0906*k6[i]+a0907*k7[i]+a0908*k8[i])
  end
  f(@muladd(t + c9*dt),tmp,k9)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a1001*k1[i]+a1004*k4[i]+a1005*k5[i]+a1006*k6[i]+a1007*k7[i]+a1008*k8[i]+a1009*k9[i])
  end
  f(@muladd(t + c10*dt),tmp,k10)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a1101*k1[i]+a1104*k4[i]+a1105*k5[i]+a1106*k6[i]+a1107*k7[i]+a1108*k8[i]+a1109*k9[i]+a1110*k10[i])
  end
  f(@muladd(t + c11*dt),tmp,k11)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a1201*k1[i]+a1204*k4[i]+a1205*k5[i]+a1206*k6[i]+a1207*k7[i]+a1208*k8[i]+a1209*k9[i]+a1210*k10[i]+a1211*k11[i])
  end
  f(t+dt,tmp,k12)
  for i in uidx
    kupdate[i] = @muladd b1*k1[i]+b6*k6[i]+b7*k7[i]+b8*k8[i]+b9*k9[i]+b10*k10[i]+b11*k11[i]+b12*k12[i]
    update[i] = dt*kupdate[i]
    u[i] = uprev[i] + update[i]
  end
  if integrator.opts.adaptive
    for i in uidx
      atmp[i] = (dt*(@muladd(k1[i]*er1 + k6[i]*er6 + k7[i]*er7 + k8[i]*er8 + k9[i]*er9 + k10[i]*er10 + k11[i]*er11 + k12[i]*er12))./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
      atmp2[i]= (@muladd(update[i] - dt*(bhh1*k1[i] + bhh2*k9[i] + bhh3*k12[i]))./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
    end
    err5 = integrator.opts.internalnorm(atmp) # Order 5
    err3 = integrator.opts.internalnorm(atmp2) # Order 3
    err52 = err5*err5
    if err5 ≈ 0 && err3 ≈ 0
      integrator.EEst = zero(typeof(integrator.EEst))
    else
      integrator.EEst = err52./sqrt(err52 + 0.01*err3*err3)
    end
  end
  f(t+dt,u,k13)
  if integrator.opts.calck
    @unpack c14,c15,c16,a1401,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1613,a1614,a1615 = cache.tab
    @unpack d401,d406,d407,d408,d409,d410,d411,d412,d413,d414,d415,d416,d501,d506,d507,d508,d509,d510,d511,d512,d513,d514,d515,d516,d601,d606,d607,d608,d609,d610,d611,d612,d613,d614,d615,d616,d701,d706,d707,d708,d709,d710,d711,d712,d713,d714,d715,d716 = cache.tab
    for i in uidx
      tmp[i] = @muladd uprev[i]+dt*(a1401*k1[i]+a1407*k7[i]+a1408*k8[i]+a1409*k9[i]+a1410*k10[i]+a1411*k11[i]+a1412*k12[i]+a1413*k13[i])
    end
    f(@muladd(t + c14*dt),tmp,k14)
    for i in uidx
      tmp[i] = @muladd uprev[i]+dt*(a1501*k1[i]+a1506*k6[i]+a1507*k7[i]+a1508*k8[i]+a1511*k11[i]+a1512*k12[i]+a1513*k13[i]+a1514*k14[i])
    end
    f(@muladd(t + c15*dt),tmp,k15)
    for i in uidx
      tmp[i] = @muladd uprev[i]+dt*(a1601*k1[i]+a1606*k6[i]+a1607*k7[i]+a1608*k8[i]+a1609*k9[i]+a1613*k13[i]+a1614*k14[i]+a1615*k15[i])
    end
    f(@muladd(t + c16*dt),tmp,k16)
    for i in uidx
      udiff[i]= kupdate[i]
      bspl[i] = k1[i] - udiff[i]
      integrator.k[3][i] = udiff[i] - k13[i] - bspl[i]
      integrator.k[4][i] = @muladd(d401*k1[i]+d406*k6[i]+d407*k7[i]+d408*k8[i]+d409*k9[i]+d410*k10[i]+d411*k11[i]+d412*k12[i]+d413*k13[i]+d414*k14[i]+d415*k15[i]+d416*k16[i])
      integrator.k[5][i] = @muladd(d501*k1[i]+d506*k6[i]+d507*k7[i]+d508*k8[i]+d509*k9[i]+d510*k10[i]+d511*k11[i]+d512*k12[i]+d513*k13[i]+d514*k14[i]+d515*k15[i]+d516*k16[i])
      integrator.k[6][i] = @muladd(d601*k1[i]+d606*k6[i]+d607*k7[i]+d608*k8[i]+d609*k9[i]+d610*k10[i]+d611*k11[i]+d612*k12[i]+d613*k13[i]+d614*k14[i]+d615*k15[i]+d616*k16[i])
      integrator.k[7][i] = @muladd(d701*k1[i]+d706*k6[i]+d707*k7[i]+d708*k8[i]+d709*k9[i]+d710*k10[i]+d711*k11[i]+d712*k12[i]+d713*k13[i]+d714*k14[i]+d715*k15[i]+d716*k16[i])
    end
  end
  @pack integrator = t,dt,u,k
end

@inline function initialize!(integrator,cache::TsitPap8ConstantCache,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::TsitPap8ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  @unpack c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1310,b1,b6,b7,b8,b9,b10,b11,b12,bhat1,bhat6,bhat7,bhat8,bhat9,bhat10,bhat13 = cache
  k1 = integrator.fsalfirst
  a = dt*a0201
  k2 = f(@muladd(t + c1*dt),uprev+a*k1)
  k3 = f(@muladd(t + c2*dt),@muladd(uprev+dt*(a0301*k1+a0302*k2)))
  k4 = f(@muladd(t + c3*dt),@muladd(uprev+dt*(a0401*k1       +a0403*k3)))
  k5 = f(@muladd(t + c4*dt),@muladd(uprev+dt*(a0501*k1       +a0503*k3+a0504*k4)))
  k6 = f(@muladd(t + c5*dt),@muladd(uprev+dt*(a0601*k1                +a0604*k4+a0605*k5)))
  k7 = f(@muladd(t + c6*dt),@muladd(uprev+dt*(a0701*k1                +a0704*k4+a0705*k5+a0706*k6)))
  k8 = f(@muladd(t + c7*dt),@muladd(uprev+dt*(a0801*k1                +a0804*k4+a0805*k5+a0806*k6+a0807*k7)))
  k9 = f(@muladd(t + c8*dt),@muladd(uprev+dt*(a0901*k1                +a0904*k4+a0905*k5+a0906*k6+a0907*k7+a0908*k8)))
  k10 =f(@muladd(t + c9*dt),@muladd(uprev+dt*(a1001*k1                +a1004*k4+a1005*k5+a1006*k6+a1007*k7+a1008*k8+a1009*k9)))
  k11= f(@muladd(t + c10*dt),@muladd(uprev+dt*(a1101*k1                +a1104*k4+a1105*k5+a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10)))
  k12= f(t+dt,@muladd(uprev+dt*(a1201*k1                +a1204*k4+a1205*k5+a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11)))
  k13= f(t+dt,@muladd(uprev+dt*(a1301*k1                +a1304*k4+a1305*k5+a1306*k6+a1307*k7+a1308*k8+a1309*k9+a1310*k10)))
  update = dt*(@muladd(b1*k1+b6*k6+b7*k7+b8*k8+b9*k9+b10*k10+b11*k11+b12*k12))
  u = uprev + update
  if integrator.opts.adaptive
    integrator.EEst = integrator.opts.internalnorm(@muladd(update - dt*(k1*bhat1 + k6*bhat6 + k7*bhat7 + k8*bhat8 + k9*bhat9 + k10*bhat10 + k13*bhat13))./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
  end
  k = f(t+dt,u)
  integrator.fsallast = k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  @pack integrator = t,dt,u
end

@inline function initialize!(integrator,cache::TsitPap8Cache,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::TsitPap8Cache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1310,b1,b6,b7,b8,b9,b10,b11,b12,bhat1,bhat6,bhat7,bhat8,bhat9,bhat10,bhat13 = cache.tab
  @unpack k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,utilde,update,tmp,atmp,k = cache
  k1 = cache.fsalfirst
  f(t,uprev,k1)
  a = dt*a0201
  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*(a0301*k1[i]+a0302*k2[i])
  end
  f(@muladd(t + c2*dt),tmp,k3)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a0401*k1[i]+a0403*k3[i])
  end
  f(@muladd(t + c3*dt),tmp,k4)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a0501*k1[i]+a0503*k3[i]+a0504*k4[i])
  end
  f(@muladd(t + c4*dt),tmp,k5)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a0601*k1[i]+a0604*k4[i]+a0605*k5[i])
  end
  f(@muladd(t + c5*dt),tmp,k6)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a0701*k1[i]+a0704*k4[i]+a0705*k5[i]+a0706*k6[i])
  end
  f(@muladd(t + c6*dt),tmp,k7)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a0801*k1[i]+a0804*k4[i]+a0805*k5[i]+a0806*k6[i]+a0807*k7[i])
  end
  f(@muladd(t + c7*dt),tmp,k8)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a0901*k1[i]+a0904*k4[i]+a0905*k5[i]+a0906*k6[i]+a0907*k7[i]+a0908*k8[i])
  end
  f(@muladd(t + c8*dt),tmp,k9)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a1001*k1[i]+a1004*k4[i]+a1005*k5[i]+a1006*k6[i]+a1007*k7[i]+a1008*k8[i]+a1009*k9[i])
  end
  f(@muladd(t + c9*dt),tmp,k10)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a1101*k1[i]+a1104*k4[i]+a1105*k5[i]+a1106*k6[i]+a1107*k7[i]+a1108*k8[i]+a1109*k9[i]+a1110*k10[i])
  end
  f(@muladd(t + c10*dt),tmp,k11)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a1201*k1[i]+a1204*k4[i]+a1205*k5[i]+a1206*k6[i]+a1207*k7[i]+a1208*k8[i]+a1209*k9[i]+a1210*k10[i]+a1211*k11[i])
  end
  f(t+dt,tmp,k12)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a1301*k1[i]+a1304*k4[i]+a1305*k5[i]+a1306*k6[i]+a1307*k7[i]+a1308*k8[i]+a1309*k9[i]+a1310*k10[i])
  end
  f(t+dt,tmp,k13)
  for i in uidx
    update[i] = dt*@muladd(b1*k1[i]+b6*k6[i]+b7*k7[i]+b8*k8[i]+b9*k9[i]+b10*k10[i]+b11*k11[i]+b12*k12[i])
    u[i] = uprev[i] + update[i]
  end
  if integrator.opts.adaptive
    for i in uidx
      atmp[i] = (@muladd(update[i] - dt*(k1[i]*bhat1 + k6[i]*bhat6 + k7[i]*bhat7 + k8[i]*bhat8 + k9[i]*bhat9 + k10[i]*bhat10 + k13[i]*bhat13))./@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


text:
import ModelSanitizer
import Test

struct Foo end

Test.@test ModelSanitizer._get_property(Foo(), :a) isa Nothing
Test.@test !ModelSanitizer._is_iterable(typeof(1.0))
Test.@test !ModelSanitizer._is_iterable(typeof('a'))
Test.@test !ModelSanitizer._is_iterable(Foo)
Test.@test !ModelSanitizer._is_indexable(typeof(1.0))
Test.@test !ModelSanitizer._is_indexable(typeof('a'))
Test.@test !ModelSanitizer._is_indexable(Foo)

ModelSanitizer._sanitize_indexable_with_check_assigned!([1,2,3], ModelSanitizer.Data[], ModelSanitizer._DataElements(0))
ModelSanitizer._sanitize_indexable_with_check_assigned!(Foo(), ModelSanitizer.Data[], ModelSanitizer._DataElements(0))
ModelSanitizer._sanitize_indexable_without_check_assigned!([1,2,3], ModelSanitizer.Data[], ModelSanitizer._DataElements(0))
ModelSanitizer._sanitize_indexable_without_check_assigned!(Foo(), ModelSanitizer.Data[], ModelSanitizer._DataElements(0))

ModelSanitizer._is_iterable(::Type{<:Foo}) = true
ModelSanitizer._sanitize_iterable!(Foo(), ModelSanitizer.Data[], ModelSanitizer._DataElements(0))

ModelSanitizer._has_isassigned(::Type{<:Foo}) = true
ModelSanitizer._sanitize_indexable!(Foo(), ModelSanitizer.Data[], ModelSanitizer._DataElements(0))

ext: jl


text:
using Combinatorics
n = 4
m = 3
for i in combinations(0:n,m)
    println(i')
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
const engvall = let res_init=zeros(5), jac_init=zeros(5,3), x_init=[1.0, 2.0, 0.0]

    function res(x, r)
        r[1] = x[1]^2 + x[2]^2 + x[3]^2 - 1.0
        r[2] = x[1]^2 + x[2]^2 + (x[3] - 2.0)^2 - 1.0
        r[3] = x[1] + x[2] + x[3] - 1.0
        r[4] = x[1] + x[2] - x[3] + 1.0
        r[5] = x[1]^3 + 3.0 * x[2] ^ 2 + (5.0 * x[3] - x[1] + 1.0)^2 - 36.0
        return r
    end

    function jac(x, jac)
        jac[1, 1] = 2.0 * x[1]
        jac[1, 2] = 2.0 * x[2]
        jac[1, 3] = 2.0 * x[3]
        jac[2, 1] = jac[1, 1]
        jac[2, 2] = jac[1, 2]
        jac[2, 3] = 2.0 * (x[3] - 2.0)
        jac[3, 1] = 1.0
        jac[3, 2] = 1.0
        jac[3, 3] = 1.0
        jac[4, 1] = 1.0
        jac[4, 2] = 1.0
        jac[4, 3] = -1.0
        t = 2.0 * (5.0 * x[3] - x[1] + 1.0)
        jac[5, 1] = 3.0 * x[1] ^ 2 - t
        jac[5, 2] = 6.0 * x[2]
        jac[5, 3] = 5.0 * t
        return jac
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("engvall", res, jac,
                                                   res_init, jac_init, x_init;
                                                   scale=scale, verbose=verbose,
                                                   print_steps=print_steps)
end


text:
using Polylogarithms
using SpecialFunctions
using Test
using DataFrames, CSV
import Base.MathConstants: π, pi, ℯ, e, γ, eulergamma, catalan, φ, golden
include("test_defs.jl")

@testset "Bernoulli numbers" begin
    @testset "    throws errors" begin
        @test_throws DomainError bernoulli(-1)
        @test_throws DomainError bernoulli(36)
        @test_throws MethodError bernoulli(1.0)
    end

    @testset "    types" begin
        @test typeof(bernoulli(2)) == Rational{Int64}
    end

    @testset "    explicit cases" begin
        for i=2:34
            # zeta isn't quite as accurate as we would like
            @test Float64(bernoulli(i)) ≈ -i*zeta( 1 - i )
        end
    end

    @testset "    dataset 1" begin
        B = Symbol("B_n")
        data1 = CSV.read("../data/bernoulli_test_data_1.csv", DataFrame)
        m = size(data1,1)
        for i=1:m
            @test bernoulli( data1[i,:n] ) ≅ parse(Rational{Int64},data1[i,B])
        end
    end
end

@testset "Bernoulli polynomials" begin
    @testset "    throws errors" begin
        @test_throws DomainError bernoulli(-1, 0.0)
        @test_throws MethodError bernoulli(1.0, 0.0)
        @test_throws MethodError bernoulli(1.0, 1)
    end

    @testset "    types" begin
        @test typeof(bernoulli(2, 1.0)) == Float64
    end

    @testset "    consistency with Bernoulli numbers" begin
        for i=0:14
            @test Float64(bernoulli(i)) ≅ bernoulli(i, 0.0)
        end
        for i=15:34
            # somewhat reliant here on accuracy of SpecialFunctions, so can't make these better without rewriting that
            @test Float64(bernoulli(i)) ≈ bernoulli(i, 0.0)
        end
     end

    @testset "    explicit cases" begin
        X = collect(-2.0 : 0.1 : 3.0)
        for i=1:length(X)
            x = X[i]
            # println("   x = $x")
            @test bernoulli(0, x) ≅ 1.0
            @test bernoulli(1, x) ≅ x   - 0.5
            @test bernoulli(2, x) ≅ x^2 -     x   + 1.0/6.0
            @test bernoulli(3, x) ≅ x^3 - 1.5*x^2 + 0.5*x
            @test bernoulli(4, x) ≅ x^4 - 2.0*x^3 +     x^2 - 1/30.0
            @test bernoulli(5, x) ≅ x^5 - 2.5*x^4 + (5.0/3.0)*x^3 - x/6.0
            # @test bernoulli(6, x) ≈ x^6 - 3.0*x^5 + (5.0/2.0)*x^4 - 0.5*x^2 + 1/42.0
            # Hurwitz-zeta doesn't work well enough here
        end
    end

    @testset "    symmetry identities" begin
        X = collect(0.0: 0.1 : 1.0)
        for i=1:length(X)
            x = X[i]
            # for n=1:9
            for n=1:6
                # println("   n = $n, x = $x")
                # start to get errors of order 1.0e-14 for n=5, 1.0e-13 around n=10, ...
                @test bernoulli(n, 1-x)  ≈  (-1.0)^n * bernoulli(n, x)
                @test bernoulli(n, x+1)  ≈  bernoulli(n, x) + n*x^(n-1)
                # SpecialFunctions.zeta(-6, -1.0) = NaN (probably should be 2)
                @test (-1)^n * bernoulli(n, -x)  ≈   bernoulli(n, x) + n*x^(n-1)

                # Raabe (1851)
                m = 6.0
                k = 0:m-1
                @test sum( bernoulli.(n, x .+ k./m ) )/m  ≈  m^(-n)*bernoulli(n, m*x )
            end
        end
    end

    @testset " dataset 2" begin
        B = Symbol("B_n(x)")
        data2 = CSV.read("../data/bernoulli_test_data_2.csv", DataFrame; delim=",", type=String)
        data2[!,:n] = parse.(Float64, data2[!,:n] )
        data2[!,:x] = parse.(Float64, data2[!,:x] )
        data2[!,B] = parse.(Float64, data2[!,B] )
        m = size(data2,1)
        for i=1:m
            @test bernoulli( Int(data2[i,:n]),  data2[i,:x]) ≈ data2[i,B]
        end
    end

    @testset " dataset 3" begin
        B = Symbol("B_n(x)")
        data3 = CSV.read("../data/bernoulli_test_data_3.csv", DataFrame; delim=",")
        m = size(data3,1)
        for i=1:m
            @test bernoulli( Int(data3[i,:n]),  data3[i,:x]) ≈ data3[i,B]
        end
   end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

function Base.length(blob::GitBlob)
    ensure_initialized()
    return ccall((:git_blob_rawsize, :libgit2), Int64, (Ptr{Cvoid},), blob.ptr)
end

"""
    rawcontent(blob::GitBlob) -> Vector{UInt8}

Fetch the *raw* contents of the [`GitBlob`](@ref) `blob`. This is an
`Array` containing the contents of the blob, which may be binary or may be Unicode.
If you write to this `Array` the blob on disk will not be updated.
`rawcontent` will allow the user to load the raw binary data into
the output `Array` and will not check to ensure it is valid Unicode, so errors
may occur if the result is passed to functions which expect valid Unicode data.

See also [`content`](@ref), which *will* throw an error if the content of the `blob`
is binary and not valid Unicode.
"""
function rawcontent(blob::GitBlob)
    ensure_initialized()
    ptr = ccall((:git_blob_rawcontent, :libgit2), Ptr{UInt8}, (Ptr{Cvoid},), blob.ptr)
    copy(unsafe_wrap(Array, ptr, (length(blob),), own = false))
end

"""
    content(blob::GitBlob) -> String

Fetch the contents of the [`GitBlob`](@ref) `blob`. If the `blob` contains
binary data (which can be determined using [`isbinary`](@ref)),
throw an error. Otherwise, return a `String` containing the contents
of the `blob`.
"""
function content(blob::GitBlob)
    s = String(rawcontent(blob))
    isvalid(s) || error("Blob does not contain valid UTF-8 data")
    return s
end

"""
    isbinary(blob::GitBlob) -> Bool

Use a heuristic to guess if a file is binary: searching for NULL bytes and
looking for a reasonable ratio of printable to non-printable characters among
the first 8000 bytes.
"""
function isbinary(blob::GitBlob)
    ensure_initialized()
    bin_flag = ccall((:git_blob_is_binary, :libgit2), Cint, (Ptr{Cvoid},), blob.ptr)
    return bin_flag == 1
end

"""
    LibGit2.addblob!(repo::GitRepo, path::AbstractString)

Read the file at `path` and adds it to the object database of `repo` as a loose blob.
Return the [`GitHash`](@ref) of the resulting blob.

# Examples
```julia
hash_str = string(commit_oid)
blob_file = joinpath(repo_path, ".git", "objects", hash_str[1:2], hash_str[3:end])
id = LibGit2.addblob!(repo, blob_file)
```
"""
function addblob!(repo::GitRepo, path::AbstractString)
    ensure_initialized()
    id_ref = Ref{GitHash}()
    @check ccall((:git_blob_create_from_disk, :libgit2), Cint,
                 (Ptr{GitHash}, Ptr{Cvoid}, Cstring),
                 id_ref, repo.ptr, path)
    return id_ref[]
end

function Base.show(io::IO, blob::GitBlob)
    if !isbinary(blob)
        conts   = split(content(blob), "\n")
        showlen = max(length(conts), 3)
        println(io, "GitBlob:\nBlob id: ", GitHash(blob), "\nContents:")
        for i in 1:showlen
            println(io, conts[i])
        end
    else
        println(io, "GitBlob:\nBlob id: ", GitHash(blob), "\nContents are binary.")
    end
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using Colors, Highlights
using Highlights.Format
using Highlights.Tokens
using Highlights.Themes

css2color(str) = parse(RGBA{Float32}, string("#", str))
css2color(c::Themes.RGB) = RGBA{Float32}(c.r / 255, c.g / 255, c.b / 255, 1.0)

function style2color(style, default)
    if style.fg.active #has_fg(style)
        css2color(style.fg)
    else
        default
    end
end

function render_str(
    ctx::Format.Context, theme::Format.Theme
)
    defaultcolor = if theme.base.fg.active # #has_fg(theme.base)
        css2color(theme.base.fg)
    else
        RGBA(0.0f0, 0.0f0, 0.0f0, 1.0f0)
    end
    colormap = map(s -> style2color(s, defaultcolor), theme.styles)
    tocolor = Dict(zip(Tokens.__TOKENS__, colormap))
    colors = RGBA{Float32}[]
    io = IOBuffer()
    for token in ctx.tokens
        t = Tokens.__TOKENS__[token.value.value]
        str = SubString(ctx.source, token.first, token.last)
        print(io, str)
        append!(colors, fill(tocolor[t], length(str)))
        # push!(colors, tocolor[t])
    end
    String(take!(io)), colors
end

function highlight_text(src::AbstractString, theme = Themes.DefaultTheme)
    io = IOBuffer()
    render_str(
        Highlights.Compiler.lex(src, Lexers.JuliaLexer),
        Themes.theme(theme)
    )
end


text:
using Plots

include("CollisionMatrix.jl")
include("Thresholds.jl")

"""
    sort_decreasing(arr)

Returns the indices needed to sort input array in decreasing order.
"""
function sort_decreasing(arr)
    liste = copy(arr)
    sorted_indexes = Int.(zeros(length(arr)))
    for i in 1:length(liste)
        max_idx = findfirst(isequal(maximum(liste)),liste)
        sorted_indexes[i] = max_idx
        liste[max_idx] = 0
    end
    return sorted_indexes
end

"""
    A class whose instances contain all usefull informations about a given motif
found in a categorical time-series.
"""
struct pattern
    shape
    instances
    positions
end

"""
    detect_motifs(ts, w, d, t = w - d; iters = 1000, tolerance = 0.95)

Detects all motifs of length 'w' occuring more often than chance, being identical up to 'd' differences.
Input:

    w : length of motifs to look for
    d : allowed errors (differences) between motifs
    t : size of the masks to use for random projection in the detection
    iters : the numbers of iterations for the random projection process (defaults to 1000)
    tolerance : threshold of motif identification.

returns :
    motifs : a list of motif instances containing all usefull informations about the found motifs (positions, frequency, shapes ...)
"""
function detect_motifs(ts, w, d, t = w - d; iters = 800, tolerance = 0.7)
    CM = collision_matrix(ts, w, d, t; iters = iters)
    min_occurences = least_occurence_threshold(ts, w, d)
    # println("Expected similar matches by chance : $min_occurences")
    expected_matches = iters * masked_match(w, d, t) * tolerance
    occurences = sum(CM .> expected_matches, dims = 1) .+ 1
    idx = sort_decreasing(occurences[1,:]) #the indices needed to sort the final dict and list
    sorted_indices = [idx[i] for i in 1:length(idx)-1 if abs(idx[i+1] - idx[i]) > w] # removing trivial matches
    if isempty(occurences[occurences .>= min_occurences])
        println("No motif found more frequently than chance.")
    end
    list = []
    positions = []
    for j in sorted_indices
        tmp_list = []
        tmp_positions = []
        for i in 1:size(CM)[1]-div(w,2)
            if CM[i,j] > expected_matches
                push!(tmp_list, ts[i:i+w-1]); push!(tmp_positions, i)
                CM[i+1:i+div(w,2), j+1:j+div(w,2)] .= 0 # removing trivial matches
            end
        end
        if isempty(tmp_list) == false
            push!(tmp_list, ts[j:j+w-1]); push!(list, tmp_list)
            push!(tmp_positions, j); push!(positions, tmp_positions)
        end
    end
    motifs = [pattern(list[i][end], list[i], positions[i]) for i in 1:length(list)]
    return motifs
end

"""
    find_motifs(ts, w, d)

Given a motif of shape 'shape' (array{any,1}), look for all the repetitions of it which differ only up to 'd' differences.
Input:

    ts : time-series in which to look for motifs
    shape : shape (aray{any,1}) of the motif to look for.
    d : allowed errors (differences) between motifs

returns :
    motif : an instance of 'pattern' containing the found repetition of the input 'shape'.

"""
function find_motifs(ts, shape, d)
    w = length(shape)
    ts_copy = copy(ts)
    found_motifs = []
    positions = []
    for err in 0:d
        for i in 1:length(ts)-w
            if error_dist(ts_copy[i:i+w-1], shape) <= err
                push!(found_motifs, ts_copy[i:i+w-1]); push!(positions, i)
                ts_copy[i:i+w-1] .= 0
            end
        end
    end
    return pattern(shape, found_motifs, positions)
end

"""
    plot_motif(motif, ts; subplots = false)

Plots all instances of the given motif. If the corresponding time series 'ts'
is provided, plots them on top of it, preserving the time-orderding. Otherwise,
plots all instances of 'motif' on top of each other to facilitate their comparison.
if subplots is given false, will only plot the motifs on top of the time-series.
"""
function plot_motif(motif, ts; subplots = true)
    len = length(motif.shape)
    p2 = plot(1:length(ts), ts, color = "grey", xlabel = "Time", ylabel = "Value", label = "", title = "motif positions")
    for (idx, p) in enumerate(motif.positions)
        plot!(p2, p:p+len-1, ts[p:p+len-1], lw = 3, label = "")
    end
    if subplots
        p1 = plot_motif(motif)
        a = plot(p1, p2, layout = (2,1))
    else
        a = p2
    end
    display(a)
    return a
end

function plot_motif(motif)
    a = plot(motif.shape, lw = 3, yticks = minimum(motif.shape):maximum(motif.shape), label = "#1", title = "Occurences of input motif")
    for (idx, m) in enumerate(motif.instances[1:end-1])
        idx += 1
        plot!(a, m, lw = 3, label = "#$idx")
    end
    display(a)
    return a
end

"""
    map_to_integer(ts)
returns an integer mapping of input time-series 'ts'.
"""
function map_to_integer(ts)
    categories = unique(ts)
    arr = copy(ts)
    for (idx, elt) in enumerate(arr)
        for (pos, c) in enumerate(categories)
            if elt == c
                arr[idx] = pos
            end
        end
    end
    return arr
end

 export detect_motifs, map_to_integer, find_motifs, plot_motif, expected_matches, least_occurence_threshold


text:
using DrWatson
using DataFrames
using CSV
using BSON
using ValueHistories
# because of vae
using GenerativeAD
using GenerativeAD.Models
using ConditionalDists
using GenerativeModels
import GenerativeModels: VAE
using Distributions
using DistributionsAD

"""
	helper functions for two stage models
"""
function string2dict(params)
	function get_values(value)
		try
			return parse(Int64, value)
		catch e
			if occursin("f0", value)
				return parse(Float32, value[1:end-2])
			else
				return value[2:end-1]
			end
		end
	end

	params = map(x->split(x, " = "), split(params, ","))
	key = map(x->x[1][2:end], params)
	value = map(x->get_values(x[2]), params)
	return Dict(map(x->x[1]=>x[2], zip(key, value)))
end

function return_best_n(df, rev=false, n=10, dataset="MNIST")
	df = df[df.dataset .== dataset, :]
	top = first(sort(by(df, :params, :criterion => mean), :criterion_mean, rev=rev), n)
	df_top = []
	for i=1:n
		tmp = df[df.params .== top.params[i],:]
		ind = i.*ones(size(tmp,1))
		push!(df_top, hcat(tmp, DataFrame(ind=ind)))
	end
	return vcat(df_top...)
end

function load_encoding(model, data; dataset::String="iris", seed::Int=1, model_index::Int=1)
	# model = "vae_AUC_tabular"
	# load csv
	df = CSV.read(datadir("tables/$(model)_best_tab.csv"))
	df = df[df.dataset .== dataset, :]
	#df = return_best_n(df, 10, dataset)
	# checking if model configuration was trained on all classes and all seeds
	n_comb = 50
	check_comb = sum([sum(df.ind .== i) for i=1:10])
	if check_comb < n_comb
		@info "One of chosen models does not include all seeds!! $(check_comb) out of $(n_comb)"
	end
	# get correct model path
	encoding_path = df[(df.seed .==seed) .& (df.ind .==model_index),:][:path][1]
	encoding_info = df[(df.seed .==seed) .& (df.ind .==model_index),:][:info_path][1]
	#encoder_params = string2dict(df[(df.seed .==seed) .& (df.ind .==model_index),:][:params][1])
	encoder_params = BSON.load(encoding_info)[:parameters]
	# load model and encode data
	model = BSON.load(encoding_path)
	encodings = map(x->cpu(GenerativeAD.Models.encode_mean(model["model"], x)), (data[1][1], data[2][1], data[3][1]))
	data = (encodings[1], data[1][2]), (encodings[2], data[2][2]), (encodings[3], data[3][2])

	return data, split(encoding_path, "/")[end], encoder_params, model["fit_t"]
end

function load_encoding(model, data, anomaly_class; dataset::String="MNIST", seed::Int=1, model_index::Int=1)
	# model ="vae_AUC_images"
	# load csv
	df = CSV.read(datadir("tables/$(model)_best_tab.csv"))
	#df = return_best_n(df, 10, dataset)
	df = df[df.dataset .== dataset, :]
	# checking if model configuration was trained on all classes and all seeds
	n_comb = 100
	check_comb = sum([sum(df.ind .== i) for i=1:10])
	if check_comb < n_comb
		@info "One of chosen models does not include all anomaly classes and all seeds!! $(check_comb) out of $(n_comb)"
	end
	# get correct model path
	encoding_path = df[(df.ac .== anomaly_class) .& (df.seed .==seed) .& (df.ind .==model_index),:][:path][1]
	encoding_info = df[(df.ac .== anomaly_class) .& (df.seed .==seed) .& (df.ind .==model_index),:][:info_path][1]
	#encoder_params = string2dict(df[(df.ac .== anomaly_class) .& (df.seed .==seed) .& (df.ind .==model_index),:][:params][1])
	encoder_params = BSON.load(encoding_info)[:parameters]
	# load model and encodings
	model = BSON.load(encoding_path)
	data = (model["tr_encodings"], data[1][2]), (model["val_encodings"], data[2][2]), ( model["tst_encodings"], data[3][2])

	return data, split(encoding_path, "/")[end], encoder_params
end


text:
abstract type MSM_Timeseries{T<:AbstractFloat} end

struct MSM_Timeseries_Point{T<:AbstractFloat} <: MSM_Timeseries{T}
    x::AbstractArray{T,2}
    residuals::AbstractArray{T,3}
    timesteps::AbstractArray{T,1}

    MSM_Timeseries_Point{T}(x::AbstractArray{T,2},
        residuals::AbstractArray{T,3},
        timesteps::AbstractArray{T,1}) where T <: Real =
            new(x,residuals,timesteps)
end

struct MSM_Timeseries_Dist{T<:AbstractFloat} <: MSM_Timeseries{T}
    x::AbstractArray{T,3}
    residuals::AbstractArray{T,4}
    timesteps::AbstractArray{T,1}

    ts_points::Array{MSM_Timeseries_Point{T},1}

    MSM_Timeseries_Dist{T}(x::AbstractArray{T,3},
        residuals::AbstractArray{T,4},
        timesteps::AbstractArray{T,1}) where T <: Real =
        new(x,residuals,timesteps,
            [MSM_Timeseries_Point{T}(
                view(x,:,:,i),
                view(residuals,:,:,:,i),
                view(timesteps,:))
                for i in 1:size(x,3)
            ])
end

function MSM_Timeseries_Point{T}(num_obs::Integer, num_params::Integer,
    num_layers::Integer, timesteps::Array{T,1}) where T <: Real

    MSM_Timeseries_Point{T}(Array{T}(undef,num_obs, num_params),
        Array{T}(undef, num_obs, num_params, num_layers), timesteps)
end

function MSM_Timeseries_Point{T}(num_obs::Integer, num_params::Integer,
    num_layers::Integer, timeStep::T = one(T)) where T <: Real

    MSM_Timeseries_Point{T}(num_obs, num_params, num_layers,
        repeat([timeStep], num_obs))
end

MSM_Timeseries_Point{T}(timeseries::MSM_Timeseries_Dist{T},
    aggregate_fun::Function)  where T <: Real =
        MSM_Timeseries_Point{T}(
            mapslices(aggregate_fun, timeseries.x; dims=(3))[:,:,1],
            mapslices(aggregate_fun, timeseries.residuals; dims=(4))[:,:,:,1],
            timeseries.timesteps
        )

MSM_Timeseries_Point{T}(timeseries::Array{MSM_Timeseries_Dist{T},1},
    aggregate_fun::Function)  where T <: Real =
                [MSM_Timeseries_Point{T}(ts,aggregate_fun)
                    for ts in timeseries]

function MSM_Timeseries_Dist{T}(num_obs::Integer, num_params::Integer,
    num_layers::Integer, num_samples::Integer, timesteps::Array{T,1}) where T <: Real

    MSM_Timeseries_Dist{T}(Array{T}(undef,num_obs, num_params, num_samples),
        Array{T}(undef, num_obs, num_params, num_layers, num_samples), timesteps)
end

function MSM_Timeseries(timeseries::S,
    num_preds::Integer,
    timeStep::T = one(T),
    num_samples::Integer = 1) where S <: MSM_Timeseries{T} where T<:Real

    MSM_Timeseries(S, params(timeseries),
        layers(timeseries), repeat([timeStep], num_preds), num_samples)
end

function MSM_Timeseries(timeseries::S,
    timeSteps::Array{T,1},
    num_samples::Integer = 1) where S <: MSM_Timeseries{T} where T<:Real

    MSM_Timeseries(S, params(timeseries),
        layers(timeseries), timeSteps, num_samples)
end

MSM_Timeseries(::Type{MSM_Timeseries_Point{T}},
    num_params::Integer, num_layers::Integer, timeSteps::Array{T,1},
    num_samples::Integer = 1) where T <: Real =
    ifelse(num_samples == 1,
        MSM_Timeseries_Point{T}(
            length(timeSteps),
            num_params,
            num_layers,
            timeSteps),
        MSM_Timeseries_Dist{T}(length(timeSteps),
            num_params,
            num_layers,
            num_samples,
            timeSteps
            )
    )

    MSM_Timeseries(::Type{MSM_Timeseries_Dist{T}},
        num_params::Integer, num_layers::Integer, timeSteps::Array{T,1},
        num_samples::Integer = 1) where T <: Real =
            MSM_Timeseries_Dist{T}(length(timeSteps),
                num_params,
                num_layers,
                num_samples,
                timeSteps
                )

timesteps(ts::S) where {S<:MSM_Timeseries} = ts.timesteps
values(ts::S) where {S<:MSM_Timeseries} = ts.x
residuals(ts::S) where {S<:MSM_Timeseries}  = ts.residuals

Base.length(ts::S) where {S<:MSM_Timeseries} = size(values(ts), 1)

params(ts::S) where {S<:MSM_Timeseries} = size(values(ts), 2)
layers(ts::S) where {S<:MSM_Timeseries} = size(residuals(ts), 3)
samples(ts::MSM_Timeseries_Dist) = size(values(ts),3)

function Base.copy(ts::S) where {S<:MSM_Timeseries{T}} where {T <:Real}
    S(copy(ts.x), copy(ts.residuals), copy(ts.timesteps))
end

function DataFrame(ts::MSM_Timeseries_Point)
    df = DataFrames.DataFrame()

    for p in 1:params(ts)
        df[!,Symbol("x-" * string(p))] = ts.x[:,p]
    end

    for p in 1:params(ts)
        for l in 1:layers(ts)
            df[!,Symbol("r-" * string(p) * "-" * string(l))] =
                ts.residuals[:,p,l]
        end
    end

    df[!,Symbol("t")] = accumulate(+, ts.timesteps) .- ts.timesteps[1]

    df
end

function DataFrame(ts::MSM_Timeseries_Dist)
    DataFrame.(ts.ts_points)
end

write(output::String, ts::MSM_Timeseries_Point) =
    CSV.write(output * ".csv", DataFrame(ts))

function write(output::String, ts::MSM_Timeseries_Dist)
    dfs = DataFrame(ts)

    for i in 1:length(dfs)
        CSV.write(output * "-" * string(i) * ".csv", dfs[i])
    end
end


text:
djldist(::T, params...) where {T <: RandVar} = djltype(T)(params...)

djltype(::Type{<:Prim.Uniform}) = Djl.Uniform
djltype(::Type{<:Prim.Normal}) = Djl.Normal
djltype(::Type{<:Prim.Beta}) = Djl.Beta
djltype(::Type{<:Prim.Bernoulli}) = Djl.Bernoulli
djltype(::Type{<:ReplaceRandVar{Prim}}) where Prim = djltype(Prim)

mayberand(x::RandVar) = rand(x)
mayberand(c) = c

"Convert `rv` into a Distributions.jl `Distribution``"
function distributionsjl(rv::RandVar)
  θs = Prim.params(rv)
  θisconst = Prim.isconstant.(θs)
  @pre all(θisconst) "All params must be constant to convert to Distributions"
  θsc = mayberand.(θs)
  djldist(rv, θsc...)
end

# distribution(f::F, args) where F = djltype(F)(args...)
# djltype(::typeof(normal)) = Distributions.Normal
# djltype(::typeof(betarv)) = Distributions.Beta
# djltype(::typeof(uniform)) = Distributions.Uniform


text:
function forward(backend::GPUBackend, state::TiedInnerProductLayerState, inputs::Vector{Blob})
  recon_dim  = size(state.W, 1)
  hidden_dim = size(state.W, 2)
  dtype = eltype(state.W)
  for i = 1:length(inputs)
    input = inputs[i]
    N = get_num(input)   # batch size
    output = state.blobs[i]
    # output = (W^T)^T * X
    CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_N, CuBLAS.OP_N, recon_dim, N, hidden_dim, one(dtype),
                state.W.ptr, recon_dim, input.ptr, hidden_dim, zero(dtype), output.ptr, recon_dim)
    # output += bias
    CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_N, CuBLAS.OP_N, recon_dim, N, 1, one(dtype),
                state.b.ptr, recon_dim, state.bias_multipliers[i].ptr, 1, one(dtype), output.ptr, recon_dim)
  end
end

function backward(backend::GPUBackend, state::TiedInnerProductLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  recon_dim  = size(state.W, 1)
  hidden_dim = size(state.W, 2)
  data_type  = eltype(state.W)

  # used in BLAS, at first it is zero, indicating overwriting the data
  # then it becomes one, indicating adding to the data
  zero_and_then_one = zero(data_type)

  for i = 1:length(inputs)
    input = inputs[i]
    batch_size = get_num(input)
    ∂f_∂o = state.blobs_diff[i]

    if !state.frozen
      # ∂f/∂b = sum(∂f/∂o, 2)
      CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_N, CuBLAS.OP_N, recon_dim, 1, batch_size,
          convert(data_type, 1), ∂f_∂o.ptr, recon_dim, state.bias_multipliers[i].ptr, batch_size, zero_and_then_one, state.∇b.ptr, recon_dim)
    end

    zero_and_then_one = one(data_type)

    # if back propagate down
    if isa(diffs[i], CuTensorBlob)
      # ∂f/∂x = W^T * [∂f/∂o]
      CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_T, CuBLAS.OP_N, hidden_dim, batch_size, recon_dim,
          convert(data_type, 1), state.W.ptr, recon_dim, ∂f_∂o.ptr, recon_dim, convert(data_type, 0), diffs[i].ptr, hidden_dim)
    end
  end
end


text:
export Analytical, CG,
        Newton, NewtonCG,
        LBFGS,
        ProxGrad, FISTA, ISTA,
        IWLSCG
        # ADMM, FADMM # NOTE: these do not work currently

# =====
# TODO
# * all - pick linesearch
# * NewtonCG number of inner iter
# * FISTA field to enforce descent
# ====

abstract type Solver end

# ===================== analytical.jl

"""
$SIGNATURES

Analytical solver (Cholesky). If the `iterative` parameter is set to `true`
then a CG solver is used. The CG solver is matrix-free and should be preferred
in "large scale" cases (when the hat matrix `X'X` is "big").

## Parameters

* `iterative` (Bool): whether to use CG (iterative) or not
* `max_inner` (Int): in the iterative mode, how many inner iterations to do.
"""
@with_kw struct Analytical <: Solver
    iterative::Bool = false
    max_inner::Int  = 200
end
CG() = Analytical(; iterative=true)

# ===================== newton.jl

"""
$SIGNATURES

Newton solver. This is a full Hessian solver and should be avoided for
"large scale" cases.
"""
struct Newton <: Solver end

"""
$SIGNATURES

Newton CG solver. This is the same as the Newton solver except that instead
of solving systems of the form `H\\b` where `H` is the full Hessian, it uses
a matrix-free conjugate gradient approach to solving that system. This should
generally be preferred for larger scale cases.
"""
struct NewtonCG <: Solver end

"""
$SIGNATURES

LBFGS quasi-Newton solver. See [the wikipedia entry](https://en.wikipedia.org/wiki/Limited-memory_BFGS).
"""
struct LBFGS <: Solver end

# struct BFGS <: Solver end

# ===================== pgrad.jl

"""
$SIGNATURES

Proximal Gradient solver for non-smooth objective functions.

## Parameters

* `accel` (Bool): whether to use Nesterov-style acceleration
* `max_iter` (Int): number of overall iterations
* `tol` (Float64): tolerance for the relative change θ ie `norm(θ-θ_)/norm(θ)`
* `max_inner`: number of inner steps when searching for a stepsize in the
               backtracking step
* `beta`: rate of shrinkage in the backtracking step (between 0 and 1)
"""
@with_kw struct ProxGrad <: Solver
    accel::Bool    = false # use Nesterov style acceleration (see also FISTA)
    max_iter::Int  = 1000  # max number of overall iterations
    tol::Float64   = 1e-4  # tol relative change of θ i.e. norm(θ-θ_)/norm(θ)
    max_inner::Int = 100   # β^max_inner should be > 1e-10
    beta::Float64  = 0.8   # in (0, 1); shrinkage in the backtracking step
end

FISTA(; kwa...) = ProxGrad(;accel = true, kwa...)
ISTA(; kwa...)  = ProxGrad(;accel = false, kwa...)

# ===================== iwls.jl

"""
$SIGNATURES

Iteratively Reweighted Least Square with Conjugate Gradient. This is the
standard (expensive) IWLS but with more efficient solves to avoid full matrix
computations.

## Parameters

* `max_iter` (Int): number of max iterations (outer)
* `max_inner` (Int): number of iterations for the CG solves
* `tol` (Float64): tolerance for the relative change θ ie `norm(θ-θ_)/norm(θ)`
* `damping` (Float64): how much to trust iterates (1=full trust)
* `threshold` (Float64): threshold for the residuals
"""
@with_kw struct IWLSCG <: Solver
    max_iter::Int      = 100
    max_inner::Int     = 200
    tol::Float64       = 1e-4
    damping::Float64   = 1.0   # should be between 0 and 1, 1 = trust iterates
    threshold::Float64 = 1e-6  # thresh for residuals; used eg in quantile reg
end

# ===================== admm.jl

# @with_kw struct ADMM <: Solver
#     max_iter::Int  = 100
#     tol::Float64   = 1e-4
#     alpha::Float64 = 1.5  # over-relaxation (recommended between 1.5 and 1.8)
#     rho::Float64   = 1.0  # Lagrangian parameter (should be decreased if poor condition)
# end
#
# @with_kw struct FADMM <: Solver
#     max_iter::Int  = 100
#     tol::Float64   = 1e-4
#     eta::Float64   = 0.999 # restart parameter
#     rho::Float64   = 1.0   # Lagrangian parameter (should be decreased if poor condition)
#     tau::Float64   = 2.0   # Increase / decrease of ρ should be > 1
#     mu::Float64    = 10.0  #
# end


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLSVCProjFO

using Test, RLinearAlgebra, LinearAlgebra, Random

@testset "LSVC Projection Full Orthogonal -- Procedural" begin

    # Supertype
    @test supertype(LinSysVecColProjFO) == LinSysVecColProjection

    # Testing context
    Random.seed!(1010)

    A = rand(15,10)
    x = rand(10)
    b = A * x

    z = rand(10)

    rsub = LinSysVecColProjFO()

    # Verify that S (storage) is nothing
    @test isnothing(rsub.S)

    # Verify that S (storage) is appropriately initialized
    @test let
        e = zeros(10)
        e[1] = 1.0
        res = dot(A[:,1], A * z - b)
        RLinearAlgebra.rsubsolve!(rsub, z, (e, A, res), 1)
        typeof(rsub.S) <: Matrix
    end

    # Verify that solution is orthogonalized against all previous hyperplanes
    @test let
        for i = 1:13
            e = zeros(10)
            e[mod(i, 1:10)] = 1.0
            res = dot(A[:,mod(i, 1:10)], A * z - b)
            RLinearAlgebra.rsubsolve!(rsub, z, (e, A, res), i)
        end

        norm(A' * (A * z - b) ) < 1e-11
    end
end

end # End Module


text:
endpoint_tuple = [
            #no required arguments: works fine
            (:get_oauth, :get_help_configuration, "help/configuration.json", nothing),
            (:get_oauth, :get_help_languages, "help/languages.json", nothing),
            (:get_oauth, :get_help_privacy, "help/privacy.json", nothing),
            (:get_oauth, :get_help_tos, "help/tos.json", nothing),
            (:get_oauth, :get_application_rate_limit_status, "application/rate_limit_status.json", nothing),
            (:get_oauth, :get_profile_banner, "users/profile_banner.json", nothing),
            (:get_oauth, :get_blocks_ids, "blocks/ids.json", nothing),
            (:get_oauth, :get_account_settings, "account/settings.json", nothing),
            (:get_oauth, :get_oembed, "statuses/oembed.json", nothing),
            (:get_oauth, :get_trends_available, "trends/available.json", nothing),
            (:get_oauth, :get_user_suggestions, "users/suggestions.json", nothing),
            (:get_oauth, :get_saved_searches_list, "saved_searches/list.json", nothing),
            (:get_oauth, :get_lists_subscribers_show, "lists/subscribers/show.json", nothing),
            (:get_oauth, :get_friendships_no_retweets, "friendships/no_retweets/ids.json", nothing),
            (:get_oauth, :get_friends_ids, "friends/ids.json", nothing),
            (:get_oauth, :get_followers_ids, "followers/ids.json", nothing),
            (:get_oauth, :get_friendships_lookup, "friendships/lookup.json", nothing),
            (:get_oauth, :get_friendships_incoming, "friendships/incoming.json", nothing),
            (:get_oauth, :get_friendships_outgoing, "friendships/outgoing.json", nothing),
            (:get_oauth, :get_friendships_show, "friendships/show.json", nothing),
            (:get_oauth, :get_direct_messages, "direct_messages.json", :Tweets),
            (:get_oauth, :get_direct_messages_sent, "direct_messages/sent.json", :Tweets),
            (:get_oauth, :get_favorites_list, "favorites/list.json", :Tweets),
            (:get_oauth, :get_friends_list, "friends/list.json", nothing),
            (:get_oauth, :get_followers_list, "followers/list.json", nothing),
            (:get_oauth, :get_lists, "lists/list.json", :Lists),
            (:get_oauth, :get_lists_statuses, "lists/statuses.json", :Tweets),
            (:get_oauth, :get_lists_memberships, "lists/memberships.json", nothing),
            (:get_oauth, :get_lists_subscribers, "lists/subscribers.json", nothing),
            (:get_oauth, :get_lists_members_show, "members/show.json", :Users),
            (:get_oauth, :get_lists_members, "lists/members.json", :Users),
            (:get_oauth, :get_lists_show, "lists/show.json", :Lists),
            (:get_oauth, :get_lists_subscriptions, "lists/subscriptions.json", nothing),
            (:get_oauth, :get_lists_ownerships, "lists/ownerships.json", nothing),
            (:get_oauth, :get_mentions_timeline, "statuses/mentions_timeline.json", :Tweets),
            (:get_oauth, :get_user_timeline, "statuses/user_timeline.json", :Tweets),
            (:get_oauth, :get_home_timeline, "statuses/home_timeline.json", :Tweets),
            (:get_oauth, :get_retweets_of_me, "statuses/retweets_of_me.json", :Tweets),
            (:get_oauth, :get_verify_credentials, "account/verify_credentials.json", :Users),
            (:get_oauth, :get_blocks_list, "blocks/list.json", nothing),
            (:get_oauth, :get_users_lookup, "users/lookup.json", :Users),
            (:get_oauth, :get_users_show, "users/show.json", :Users),
            (:get_oauth, :get_users_contributees, "users/contributees.json", nothing),
            (:get_oauth, :get_users_contributors, "users/contributors.json", nothing),
            (:post_oauth, :post_friendships_create, "friendships/create.json", :Users),
            (:post_oauth, :post_friendships_destroy, "friendships/destroy.json", :Users),
            (:post_oauth, :post_friendships_update, "friendships/update.json", nothing),
            (:post_oauth, :post_lists_members_destroy, "lists/members/destroy.json", :Lists),
            (:post_oauth, :post_lists_subscribers_create, "subscribers/create.json", :Lists),
            (:post_oauth, :post_lists_subscribers_destroy, "subscribers/destroy.json", :Lists),
            (:post_oauth, :post_lists_members_createall, "members/create_all.json", :Lists),
            (:post_oauth, :post_lists_members_create, "members/create.json", :Lists),
            (:post_oauth, :post_lists_destroy, "lists/destroy.json", :Lists),
            (:post_oauth, :post_lists_update, "lists/update.json", :Lists),
            (:post_oauth, :post_lists_members_destroy_all, "members/destroy_all.json", :Lists),
            (:post_oauth, :post_users_report_spam, "users/report_spam.json", :Users),
            (:post_oauth, :post_account_settings, "account/settings.json", nothing),
            (:post_oauth, :post_update_profile, "account/update_profile.json", nothing),
            (:post_oauth, :post_update_profile_background, "account/update_profile_background_image.json", nothing),
            (:post_oauth, :post_blocks_create, "blocks/create.json", nothing),
            (:post_oauth, :post_blocks_destroy, "blocks/destroy.json", nothing),
            (:post_oauth, :post_account_remove_profile_banner, "account/remove_profile_banner.json", nothing),
            #have required arguments: no argument checking currently
            (:get_oauth, :get_direct_messages_show, "direct_messages/events/show.json", :Tweets),
            (:post_oauth, :post_direct_messages_destroy, "direct_messages/destroy.json", :Tweets),
            (:post_oauth, :post_direct_messages_send, "direct_messages/events/new.json", :Tweets),
            (:post_oauth, :post_favorites_destroy, "favorites/destroy.json", :Tweets),
            (:post_oauth, :post_favorites_create, "favorites/create.json", :Tweets),
            (:get_oauth, :get_geo_reverse_geocode, "geo/reverse_geocode.json", nothing),
            (:get_oauth, :get_geo_search, "geo/search.json", nothing),
            (:get_oauth, :get_geo_similar_places, "geo/similar_places.json", nothing),
            (:post_oauth, :post_lists_create, "lists/create.json", nothing),
            (:get_oauth, :get_search_tweets, "search/tweets.json", nothing),
            (:get_oauth, :get_trends_place, "trends/place.json", nothing),
            (:get_oauth, :get_trends_closest, "trends/closest.json", nothing),
            (:get_oauth, :get_single_tweet_id, "statuses/show.json", :Tweets),
            (:post_oauth, :post_status_update, "statuses/update.json", :Tweets),
            (:get_oauth, :get_retweeters_id, "statuses/retweeters/ids.json", nothing),
            (:get_oauth, :get_users_search, "users/search.json", :Users),
            (:post_oauth, :post_update_profile_image, "account/update_profile_image.json", nothing),
            (:post_oauth, :post_account_update_profile_banner, "account/update_profile_banner.json", nothing)
]

#dynamically build methods
for (verb, func, endp, t) in endpoint_tuple
    @eval begin
            function ($func)(;kwargs...)

                #Take kwargs array of tuples, create dict needed by OAuth process
                options = Dict{String, Any}()
                for arg in kwargs
                    options[string(arg[1])] = string(arg[2])
                end

                #Call endpoint
                r = ($verb)($"https://api.twitter.com/1.1/$endp", options)

                #If successful API call, return JSON as Julia data structure, otherwise return error
                if r.status == 200
                    success = JSON.parse(String(r.body))

                    #If type known, parse
                    if ($t) != nothing && ($t) <: TwitterType
                        return ($t)(success)
                    else
                        return success
                    end
                else
                    error("Twitter API returned $(r.status) status")
                end

            end
    end
end


text:
using LinearAlgebra

"""
This file contains a list of cost functions that are ready to be used in optimisations
I'm going to try and keep some standard inputs that make sense
"""

"""
Target-gate fidelity
"""
function C1(KT, KN, D)
    1 - abs2(tr(KT' * KN)/D)
end

"""
Target-state infidelity
"""
function C2(ψT, ψN)
    1 - abs2(ψT' * ψN)
end

"""
Control amplitudes
"""
function C3(u)
    sum(abs2.(u))
end

"""
Control variations
"""
function C4(u)
    (K, J) = size(u)
    sum(abs2.(diff(u, dims = 2)))
end

"""
Occupation of forbidden states
"""
function C5(ψF, ψj)
    # since ψj is the intermediate states we can assume its an array
    sum(abs2(ψF' * ψ) for ψ in ψj)
end

"""
Evolution time (target gate)
"""
function C6(KT, KJ, N, D)
    1 - 1/N * sum(abs2(tr(KT' * Kj)/D) for Kj in KJ)
end

"""
Evolution time (target state)
"""
function C7(ψT, ψJ, N)
    1 - 1/N * sum(abs2(ψT, ψj) for ψj in ψJ)
end


text:
save_results = false

@testset "Quadrotor" begin
# load expected results
res = load(joinpath(@__DIR__,"quadrotor.jld2"))

# Set up problem and solver
prob, opts = Problems.Quadrotor()
times = TO.get_times(prob)
u0 = zeros(prob.model)[2]
U = [u0 + 0.1*SA[sin(t), sin(t), -sin(t), -sin(t)] for t in times]
initial_controls!(prob, U)
rollout!(prob)
states(prob)[end]

solver = ALTROSolver(prob, opts, save_S=true)
n,m,N = size(solver)
ilqr = Altro.get_ilqr(solver)
Altro.initialize!(ilqr)
X = states(ilqr)
U = controls(ilqr)
@test X ≈ res["X"] atol=1e-6
@test U ≈ res["U"] atol=1e-6

# State diff Jacobian
TO.state_diff_jacobian!(ilqr.G, ilqr.model, ilqr.Z)
G = ilqr.G
@test G[1] != zeros(n,n)
@test G[N] != zeros(n,n)
@test G ≈ res["G"] atol=1e-6

# Dynamics Jacobians
TO.dynamics_expansion!(TO.integration(ilqr), ilqr.D, ilqr.model, ilqr.Z)
TO.error_expansion!(ilqr.D, ilqr.model, ilqr.G)
A_ = [Matrix(D.A_) for D in ilqr.D]
B_ = [Matrix(D.B_) for D in ilqr.D]
A = [Matrix(D.A) for D in ilqr.D]
B = [Matrix(D.B) for D in ilqr.D]
@test A_ ≈ res["A_"] atol=1e-6
@test B_ ≈ res["B_"] atol=1e-6
@test A ≈ res["A"] atol=1e-6
@test B ≈ res["B"] atol=1e-6

# Unconstrained Cost Expansion
TO.cost_expansion!(ilqr.quad_obj, ilqr.obj.obj, ilqr.Z, init=true, rezero=true)
TO.error_expansion!(ilqr.E, ilqr.quad_obj, ilqr.model, ilqr.Z, ilqr.G)
hess0 = [[E.Q E.H'; E.H E.R] for E in ilqr.E]
grad0 = [Vector([E.q; E.r]) for E in ilqr.E]
@test hess0 ≈ res["hess0"] atol=1e-6
@test grad0 ≈ res["grad0"] atol=1e-6

# AL Cost Expansion
TO.cost_expansion!(ilqr.quad_obj, ilqr.obj, ilqr.Z, init=true, rezero=true)
TO.error_expansion!(ilqr.E, ilqr.quad_obj, ilqr.model, ilqr.Z, ilqr.G)
hess = [[E.Q E.H'; E.H E.R] for E in ilqr.E]
grad = [Vector([E.q; E.r]) for E in ilqr.E]
@test hess ≈ res["hess"] atol=1e-6
@test grad ≈ res["grad"] atol=1e-6

# Backward Pass
ilqr.ρ[1] = ilqr.opts.bp_reg_min*0
ΔV = Altro.static_backwardpass!(ilqr)
K = Matrix.(ilqr.K)
d = Vector.(ilqr.d)
S = [Matrix(S.Q) for S in ilqr.S]
s = [Vector(S.q) for S in ilqr.S]
Qzz = [[E.Q E.H'; E.H E.R] for E in ilqr.Q]
Qz = [Vector([E.q; E.r]) for E in ilqr.Q]
@test norm(K - res["K"]) < 1e-6
@test norm(d - res["d"]) < 1e-5
@test S ≈ res["S"] atol=1e-6
@test s ≈ res["s"] atol=1e-6
@test Qzz ≈ res["Qzz"] atol=1e-6
@test Qz  ≈ res["Qz"] atol=1e-6
@test ΔV ≈ res["ΔV"] atol=1e-6

if save_results
    @save joinpath(@__DIR__,"quadrotor.jld2") X U G A B A_ B_ hess0 grad0 hess grad K d S s Qzz Qz ΔV
end

end


text:
using JuMP, GLPK

# wrap all of our code inside a function (for better performance)
function example_basic()

    # define model and optimizer
    model = Model(GLPK.Optimizer)

    # define variables
    @variable(model, 0 <= x <= 2)
    @variable(model, 0 <= y <= 30)

    # define objective function
    @objective(model, Max, 5x + 3y)

    # add additional constraints
    @constraint(model, 1x + 5y <= 3.0)

    # display the model
    print(model)

    # optimize the model
    JuMP.optimize!(model)

    # return and print objective function and optimal values of variables
    obj_value = JuMP.objective_value(model)
    x_value = JuMP.value(x)
    y_value = JuMP.value(y)
    println("Objective value: ", obj_value)
    println("x = ", x_value)
    println("y = ", y_value)
end

# call the function defined above
example_basic()

using JuMP, GLPK
function example_sudoku()

    # input the initial puzzle board (0s mean blanks)
    initial_grid = [
                    3 1 0 0 5 8 0 0 4;
                    0 0 9 3 2 0 0 0 0;
                    0 2 5 1 0 4 0 9 0;
                    0 0 0 0 0 0 3 8 9;
                    0 0 8 0 0 0 5 0 0;
                    5 4 6 0 0 0 0 0 0;
                    0 8 0 2 0 3 6 5 0;
                    0 0 0 0 7 1 4 0 0;
                    7 0 0 4 8 0 0 2 1
                    ]

    # use GLPK Optimizer
    model = Model(GLPK.Optimizer)

    # Set up the variables: each one can only take on binary values, so we add "Bin" to the end as a constraint
    @variable(model, x[1:9, 1:9, 1:9], Bin)

    # Add the constraints
    @constraints(model, begin
                     # Constraint 1 - Only one value appears in each cell
                     cell[i in 1:9, j in 1:9], sum(x[i, j, :]) == 1
                     # Constraint 2 - Each value appears in each row once only
                     row[i in 1:9, k in 1:9], sum(x[i, :, k]) == 1
                     # Constraint 3 - Each value appears in each column once only
                     col[j in 1:9, k in 1:9], sum(x[:, j, k]) == 1
                     # Constraint 4 - Each value appears in each 3x3 subgrid once only
                     subgrid[i=1:3:7, j=1:3:7, val=1:9], sum(x[i:i + 2, j:j + 2, val]) == 1
                 end)

    # Add additional constraints that reflect the starting point of the puzzle board
    # (i.e. don't attempt to update the numbers that were given as part of the puzzle)
    for row in 1:9, col in 1:9
        if initial_grid[row, col] != 0
            @constraint(model, x[row, col, initial_grid[row, col]] == 1)
        end
    end

    # Solve it
    JuMP.optimize!(model)

    term_status = JuMP.termination_status(model)
    primal_status = JuMP.primal_status(model)
    is_optimal = term_status == MOI.OPTIMAL

    # Check solution
    if is_optimal
        mip_solution = JuMP.value.(x)
        sol = zeros(Int, 9, 9)
        for row in 1:9, col in 1:9, val in 1:9
            if mip_solution[row, col, val] >= 0.9
                sol[row, col] = val
            end
        end
        return sol
    else
        error("The solver did not find an optimal solution.")
    end
end

function print_sudoku_solution(solution)
    println("Solution:")
    println("[-----------------------]")
    for row in 1:9
        print("[ ")
        for col in 1:9
            print(solution[row, col], " ")
            if col % 3 == 0 && col < 9
                print("| ")
            end
        end
        println("]")
        if row % 3 == 0
            println("[-----------------------]")
        end
    end
end

sol = example_sudoku()
print_sudoku_solution(sol)

using HTTP, JuliaDB, JuMP, GLPK

# first function: read in the player data from the class GitHub repository
function read_in_data(url)
    newtable  = csvread(IOBuffer(HTTP.get(url).body), skiplines_begin=0, header_exists=true)
    players   = newtable[1][2]
    salaries  = newtable[1][4]./1000000
    ppg       = newtable[1][11]
    return players,salaries,ppg
end

function SolveModel(players,salary,points)
    N = length(salary)

    m = Model(GLPK.Optimizer)

    # define the variables: they are 0 if the player did not make the team, 1 if the player did make the team
    @variable(m, picked[1:N], Bin)

    # categories:
    @objective(m, Max, sum( points[i] * picked[i] for i in 1:N))

    @constraints m begin
        # Constraint 1 - payroll <= 132.6m
        sum(salary[i] * picked[i] for i in 1:N) <= 132.6
        # Constraint 2 - must have exactly 15 players on roster
        sum(picked[i] for i in 1:N) == 15
    end

    # Solve it
    JuMP.optimize!(m);
    pck     = convert(BitArray,JuMP.value.(picked))
    lineup  = players[pck]
    points  = JuMP.objective_value(m)
    payroll = sum(salary[pck])
    return lineup,points,payroll
end

# call first function (to import data)
players,salaries,pts = read_in_data("https://raw.githubusercontent.com/tyleransom/DScourseS20/master/WebData/playerSalaryStats.csv")

# pass data into second function to get optimal lineup
lineup,total_points,payroll = SolveModel(players,salaries,pts)
println("team: ",lineup)
println("total points scored per game: ",total_points)
println("payroll: ",payroll)

# first function: read in the player data from the class GitHub repository
function read_in_data(url)
    newtable = csvread(IOBuffer(HTTP.get(url).body), skiplines_begin=0, header_exists=true)
    players  = newtable[1][2]
    salaries = newtable[1][4]./1000000
    mpg      = newtable[1][5]
    fgaG5    = 1.0.*((newtable[1][9]).>5)
    fga      = newtable[1][9]
    ppg      = newtable[1][11]
    return players,salaries,mpg,fgaG5,fga,ppg
end

function SolveModel(players,salary,minutes,field_goals_over5,field_goals,points)
    N = length(salary)

    m = Model(GLPK.Optimizer)

    # define the variables: they are 0 if the player did not make the team, 1 if the player did make the team
    @variable(m, picked[1:N], Bin)

    # categories:
    @objective(m, Max, sum( points[i] * picked[i] for i in 1:N))

    @constraints m begin
        # Constraint 1 - payroll <= 132.6m
        sum(salary[i] * picked[i] for i in 1:N) <= 132.6
        # Constraint 2 - must have exactly 15 players on roster
        sum(picked[i] for i in 1:N) == 15
        # Constraint 3 - total minutes must not exceed 240
        sum(minutes[i] * picked[i] for i in 1:N) <= 240
        # Constraint 4 - total shot attempts must be lower than 80
        sum(field_goals[i] * picked[i] for i in 1:N) <= 80
    end

    # Solve it
    JuMP.optimize!(m);
    pck      = convert(BitArray,JuMP.value.(picked))
    lineup   = players[pck]
    totmin   = sum(minutes[pck])
    payroll  = sum(salary[pck])
    totshots = sum(field_goals[pck])
    points   = JuMP.objective_value(m)
    return lineup,points,totmin,payroll,totshots
end

# call first function (to import data)
players,salaries,minutes,over5fg,fga,pts = read_in_data("https://raw.githubusercontent.com/tyleransom/DScourseS20/master/WebData/playerSalaryStats.csv")

# pass data into second function to get optimal lineup
lineup,total_points,total_minutes,payroll,totshots = SolveModel(players,salaries,minutes,over5fg,fga,pts)
println("team: ",lineup)
println("total points scored per game: ",total_points)
println("payroll: ",payroll)
println("total shots per game: ",totshots)
println("total minutes per game: ",total_minutes)

using HTTP, JuliaDB, JuMP, Ipopt
function import_auto(url)
    newtable  = csvread(IOBuffer(HTTP.get(url).body), skiplines_begin=0, header_exists=true)
    depvar    = log.(newtable[1][2]) # log price
    indepvars = cat(ones(size(depvar)),newtable[1][3],newtable[1][5],newtable[1][6]; dims=2) # constant, mpg, headroom, trunk
    return depvar,indepvars
end

Y,X = import_auto("https://tyleransom.github.io/teaching/MetricsLabs/auto.csv")

function jumpOLS(Y,X,startval=zeros(size(X,2),1))
    OLS = Model(Ipopt.Optimizer)

    # Declare the variables you are optimizing over
    @variable(OLS, b[i=1:size(X,2)], start = startval[i])

    # Write your objective function
    @NLobjective(OLS, Min, sum( (Y[i]-sum( X[i,k]*b[k] for k in 1:size(X,2) ))^2 for i in 1:size(X,1) ) )

    # Solve the objective function
    JuMP.optimize!(OLS)

    SSR = JuMP.objective_value(OLS)
    b_value = JuMP.value.(b)
    println("Objective value: ", SSR)
    println("beta hat = ", b_value)
    println("RMSE = ", sqrt(SSR/(size(X,1)-size(X,2))))
end

jumpOLS(Y,X)


text:
using Robotlib, Robotlib.Frames, MAT, Test, Plots

R = orthonormal(randn(3, 3))
t = randn(3)
f = Frame(R, t)
f = Frame([R t; 0 0 0 1])
F = f

path = joinpath(dirname(@__FILE__), "..", "applications", "frames/")
add_frame_name!("SEAM", "Weld seam frame")
add_frame_name!("TAB", "Table frame")

T_RB_Tm = MAT.matread(path * "T_RB_T.mat")["T_RB_T"]
T_TF_TCPm = MAT.matread(path * "T_TF_TCP.mat")["T_TF_TCP"]
T_T_TABm = MAT.matread(path * "T_T_Table.mat")["T_T_Table"]

T_RB_T = Frame(T_RB_Tm, "RB", "T")
T_S_D = Frame(T_TF_TCPm, "S", "D")
T_T_TAB = Frame(T_T_TABm, "T", "TAB")

cloud_seam = readcloud(path * "CloudSeam_edge.txt")
plane_seam = readplane(path * "PlaneSeam_edge.txt")
cloud_seam_projected = project(plane_seam, cloud_seam)
line_seam = fitline(cloud_seam_projected)

T_T_SEAM = framefromfeatures(
    ("z+", line_seam),
    ("y-", plane_seam),
    cloud_seam_projected[1],
    "SEAM",
)
T_RB_SEAM = T_RB_T * T_T_SEAM
T_RB_TAB = T_RB_T * T_T_TAB
T_TAB_SEAM = inv(T_T_TAB) * T_T_SEAM

cloud_seam_RB = T_RB_T * cloud_seam
cloud_seam_projected_RB = T_RB_T * cloud_seam_projected
plane_seam_RB = T_RB_T * plane_seam
line_seam_RB = T_RB_T * line_seam

plot(Frame(I4,"RB","U"),200)
plot!(cloud_seam_RB,label="Cloud seam_RB")
plot!(cloud_seam_projected_RB,label="Cloud seam_RB projected")
plot!(line_seam_RB,500,label="Line seam")
plot!(plane_seam_RB,200,label="Plane seam")
plot!(T_RB_SEAM,200)
plot!(T_RB_TAB,200)
# gui()

T = convert(Matrix, f)
@test convert(Array, f) == T
@test promote(T, f) == (T, T)
@test promote_type(Matrix{Float64}, Frame{Float64}) == Matrix{Float64}

max_forks_repo_name: baggepinnen/Robotlib


text:
using FileIO
include("julia2ims.jl")
field = ["s4-2" "s9-1" "s9-2" "s10-5" "s24-2" "s6-1" "s6-2" "s7-5" "s11-1" "s15-1" "s28-2" "s32-1" "s33-1"]
data_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF/20200318";
output_dir = "$(data_dir)_visualization"

try
    mkdir("$output_dir")
catch
end

function max_projection(cell)
    println("Processing $cell")
    #raw_img = load(File(format"TIFF", "$data_dir/HE7-11-1-80uw-const_3_$pos.ome.tiff"));
    raw_img = loadims2("$(data_dir)_result/output/d18$cell.ims");
    img_size = size(raw_img);
    img_max = zeros(N0f16, img_size[1], img_size[2], Int(img_size[3]/20));
    for i in 1:Int(img_size[3]/20)
        img_max[:, :, i] = maximum(raw_img[:, :, (i-1)*20+1:i*20], dims=3);
    end
    save("$output_dir/$cell.tiff", img_max);
end

for cell in field
    max_projection(cell)
end


text:
export elu, relu, selu, sigm, invx
using AutoGrad: AutoGrad, @primitive

"""
    elu(x)

Return `(x > 0 ? x : exp(x)-1)`.

Reference: Fast and Accurate Deep Network Learning by Exponential Linear Units (ELUs) (https://arxiv.org/abs/1511.07289).
"""
elu(x::T) where T = (x >= 0 ? x : exp(x)-T(1))
eluback(dy::T,y::T) where T = (y >= 0 ? dy : dy * (T(1)+y))
@primitive elu(x),dy,y eluback.(dy,y)
@primitive eluback(dy,y),ddx  (ddx.*((y.>=0).+(y.<0).*(y.+1)))  (ddx.*dy.*(y.<0))

"""
    relu(x)
Return `max(0,x)`.

References:
* [Nair and Hinton, 2010](https://icml.cc/Conferences/2010/abstracts.html#432). Rectified Linear Units Improve Restricted Boltzmann Machines. ICML.
* [Glorot, Bordes and Bengio, 2011](http://proceedings.mlr.press/v15/glorot11a). Deep Sparse Rectifier Neural Networks. AISTATS.
"""
relu(x::T) where T = max(x,T(0))
reluback(dy::T,y::T) where T = (y>0 ? dy : T(0))
@primitive relu(x),dy,y reluback.(dy,y)
@primitive reluback(dy,y),ddx  (ddx.*(y.>0))  nothing

"""
    selu(x)

Return `λ01 * (x > 0 ? x : α01 * (exp(x)-1))` where `λ01=1.0507009873554805` and `α01=1.6732632423543778`.

Reference: Self-Normalizing Neural Networks (https://arxiv.org/abs/1706.02515).
"""
selu(x::T) where T = (x >= 0 ? T(λ01)*x : T(λα01)*(exp(x)-T(1)))
seluback(dy::T,y::T) where T = (y >= 0 ? dy * T(λ01) : dy * (y + T(λα01)))
@primitive selu(x),dy,y seluback.(dy,y)
@primitive seluback(dy,y),ddx  (T=eltype(y); ddx.*((y.>=0).*T(λ01).+(y.<0).*(y.+T(λα01))))  (ddx.*dy.*(y.<0))
const λ01 = 1.0507009873554805  # (1-erfc(1/sqrt(2))*sqrt(exp(1)))*sqrt(2pi)*(2*erfc(sqrt(2))*exp(2)+pi*erfc(1/sqrt(2))^2*exp(1)-2*(2+pi)*erfc(1/sqrt(2))*sqrt(exp(1))+pi+2)^(-0.5)
const α01 = 1.6732632423543778  # -sqrt(2/pi)/(erfc(1/sqrt(2))*exp(1/2)-1)
const λα01 = 1.7580993408473773 # λ01 * α01

"""
    sigm(x)

Return `1/(1+exp(-x))`.

Reference: Numerically stable sigm implementation from http://timvieira.github.io/blog/post/2014/02/11/exp-normalize-trick.
"""
sigm(x::T) where T = (x >= 0 ? T(1)/(T(1)+exp(-x)) : (z=exp(x); z/(T(1)+z)))
sigmback(dy::T,y::T) where T = (dy*y*(T(1)-y))
@primitive sigm(x),dy,y  sigmback.(dy,y)
@primitive sigmback(dy,y),ddx  ddx.*y.*(1 .- y)  ddx.*dy.*(1 .- 2 .* y)

function invx(x)
    @warn "invx() is deprecated, please use 1/x instead" maxlog=1
    1/x
end


text:
using Giggle
using Test

@testset "Giggle.jl" begin
    # Write your tests here.
end


text:
using Distances
using LinearAlgebra
using PointPatternStatistics: iscovered

window = (x=(0.0, 100.0), y=(30.0, 200.0))
x = rand(SimpleSequentialInhibition(window, 10.0))
d = pairwise(Euclidean(), reshape(reinterpret(Float64, x.data), 2, :), dims=2)
d[diagind(d)] .= Inf
@test minimum(d) >= 10.0
@test all(x -> inside(window)(x) >= 0, x.data)

@test !iscovered(0.0, 1.0, 0.0, 1.0, [(0.5, 0.5)], sqrt(0.5) - 0.00001)
@test iscovered(0.0, 1.0, 0.0, 1.0, [(0.5, 0.5)], sqrt(0.5) + 0.00001)
@test iscovered(0.0, 1.0, 0.0, 1.0, [(1.0, 0.5), (0.0,0.5)], sqrt(0.5)+0.00001)
@test !iscovered(0.0, 1.0, 0.0, 1.0, [(1.0, 0.5), (0.0,0.5)], sqrt(0.5)-0.00001)
@test !iscovered(0.0, 1.0, 0.0, 1.0, [(1.5, 0.5), (-0.5,0.5),(0.5,1.5),(0.5,-0.5)], 1-0.00001)
@test iscovered(0.0, 1.0, 0.0, 1.0, [(1.5, 0.5), (-0.5,0.5),(0.5,1.5),(0.5,-0.5)], 1+0.00001)


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

JLLWrappers.@generate_main_file_header("SHOT")
JLLWrappers.@generate_main_file("SHOT", UUID("c1ab834c-c4a5-50f5-9156-8f0fe7758b0e"))
end  # module SHOT_jll


text:
# # Visualizing Acset Schemas with Graphviz
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/graphics/graphviz_wiring_diagrams.ipynb)
#
# It is convenient to visualize schemas in a non-textual way; often schemas can have many objects, homs and attributes, and it is hard to keep these all in your head.
#
# For this reason, we provide a function which takes a schema and produces a visual representation of it, using the Catlab Graphviz interface.
#
# This is as simple as the following code. These figures are by no means publication-quality, but they can be useful for interactive development.

using Catlab.Present, Catlab.Graphics, Catlab.Graphs

to_graphviz(Graphs.BasicGraphs.TheoryWeightedGraph)


text:
@testset "Test utils" begin

using HDiscontinuousGalerkin
#Test affine maps
A,b = get_affine_map([[-1.0,-1.0],[1.0,-1.0],[-1.0,1.0]],[[0.0,0.0],[1.0,0.0],[0.0,1.0]])
@test A*[0.0,0.0]+b ≈ [0.5,0.5]

end


text:
using CmdStan, Test

cnames = ["x", "y.1", "y.2", "z.1.1", "z.2.1", "z.3.1", "z.1.2", "z.2.2", "z.3.2", "k.1.1.1.1.1"]
draws = 100
vars = length(cnames)
chains = 2
chns = randn(draws, vars, chains)
nt = extract(chns, cnames)

@testset "extract" begin
    # Write your tests here.

    @test size(nt.x) == (draws, chains)
    @test size(nt.y) == (2, draws, chains)
    @test size(nt.z) == (3, 2, draws, chains)
    @test size(nt.k) == (1, 1, 1, 1, 1, draws, chains)

    key_to_idx = Dict(name => idx for (idx, name) in enumerate(cnames))

    @test nt.x[2,1] == chns[2, key_to_idx["x"], 1]
    @test nt.y[2,3,2] == chns[3, key_to_idx["y.2"], 2]
    @test nt.z[3, 1, 10, 1] == chns[10, key_to_idx["z.3.1"], 1]
    @test nt.k[1,1,1,1,1,draws,2] == chns[draws, key_to_idx["k.1.1.1.1.1"], 2]

    @test keys(nt) == (:y, :k, :z, :x)
    @test size(values(nt.z)) == (3, 2, 100, 2)
    @test size(nt.z) == (3, 2, 100, 2)

end


text:
function vtk_grid(
        dtype::VTKImageData, filename::AbstractString, Ns::Dims{N};
        origin::NTuple{N} = ntuple(d -> 0.0, Val(N)),
        spacing::NTuple{N} = ntuple(d -> 1.0, Val(N)),
        extent = Ns, whole_extent = extent, kwargs...,
    ) where {N}
    Npts = prod(Ns)
    Ncls = num_cells_structured(Ns)

    xvtk = XMLDocument()
    vtk = DatasetFile(dtype, xvtk, filename, Npts, Ncls; kwargs...)

    # VTKFile node
    xroot = vtk_xml_write_header(vtk)

    # ImageData node
    xGrid = new_child(xroot, vtk.grid_type)
    let ext = extent_attribute(whole_extent)
        set_attribute(xGrid, "WholeExtent", ext)
    end

    origin_str = _tuple_to_str3(origin, zero(eltype(origin)))
    spacing_str = _tuple_to_str3(spacing, one(eltype(origin)))

    set_attribute(xGrid, "Origin", origin_str)
    set_attribute(xGrid, "Spacing", spacing_str)

    # Piece node
    xPiece = new_child(xGrid, "Piece")
    let ext = extent_attribute(extent)
        set_attribute(xPiece, "Extent", ext)
    end

    vtk
end

_tuple_to_str3(ts::NTuple{3, T}, default::T) where {T} = join(ts, " ")

function _tuple_to_str3(ts::NTuple{N, T}, default::T) where {N, T}
    @assert N < 3
    M = 3 - N
    tt = (ts..., ntuple(d -> default, Val(M))...)
    _tuple_to_str3(tt, default)
end

function vtk_grid(filename::AbstractString, Ns::Vararg{Integer}; kwargs...)
    # We put the origin at (0.0, 0.0, 0.0)
    vtk_grid(filename, map(N -> range(0.0; length = N), Ns)...; kwargs...)
end

"""
    vtk_grid(filename, x::AbstractRange{T}, y::AbstractRange{T}, [z::AbstractRange{T}];
             kwargs...)

Create image data (`.vti`) file.

Along each direction, the grid is specified in terms of an AbstractRange object.

# Examples

```jldoctest
julia> vtk = vtk_grid("abc", 1:0.2:5, 2:1.:3, 4:1.:5)  # 3D dataset
VTK file 'abc.vti' (ImageData file, open)

julia> vtk = vtk_grid("abc", 1:5, 2:1.:3, range(4, 5; length = 3))  # different kinds of ranges
VTK file 'abc.vti' (ImageData file, open)

julia> vtk = vtk_grid("abc", 1:0.2:5, 2:1.:3)  # 2D dataset
VTK file 'abc.vti' (ImageData file, open)

julia> vtk = vtk_grid("def",
                      LinRange(0., 5., 10),
                      LinRange(0., 2π, 16),
                      LinRange(1., 10., 12))
VTK file 'def.vti' (ImageData file, open)

```

"""
function vtk_grid(filename::AbstractString, xyz::Vararg{AbstractRange}; kwargs...)
    Nxyz = promote(length.(xyz)...)
    spacing = promote(step.(xyz)...)
    origin = promote(first.(xyz)...)

    if (extent = get(kwargs, :extent, nothing)) !== nothing
        # Shift origin accordingly
        length.(extent) == Nxyz ||
            throw(DimensionMismatch("`extent` argument doesn't match grid dimensions"))
        origin_new = @. origin + (1 - first(extent)) * spacing
        origin = oftype(origin, origin_new)
    end

    vtk_grid(VTKImageData(), filename, Nxyz;
             origin=origin, spacing=spacing, kwargs...)
end


text:
"""
    pore_size(array, phase = 0; nbins = 10, periodic = false)

Calculate pore size correlation function for one-, two- or
three-dimensional multiphase system.

Pore size correlation function `P(x)` equals to probability of
inserting a ball with radius `R ∈ [x, x + δx]` into a system so that
it lies entirely in the phase `phase`.

This implementation divides the range of possible radii into `nbins`
subranges and returns a normalized histogram of radii. This is roughly
equal to integrating `P(x)` for each subrange.
"""
function pore_size(array    :: AbstractArray,
                   phase               = 0;
                   nbins    :: Integer = 10,
                   periodic :: Bool    = false)
    indicator = map(x -> x ≠ phase, array)
    distances = distance_transform(indicator, periodic ? Torus() : Plane())
    distances = filter(x -> x != 0, distances)
    hist = fit(Histogram, distances; nbins = nbins)
    return normalize(hist; mode = :probability)
end


text:
using Luna

a = 13e-6
gas = :Ar
pres = 7.5
L = 15e-2

τfwhm = 30e-15
λ0 = 800e-9
energy = 1e-6

coren, densityfun = Capillary.gradient(gas, L, pres, 0);

modes = (
    Capillary.MarcatiliMode(a, coren, n=1, m=1, kind=:HE, ϕ=0.0, loss=false),
    Capillary.MarcatiliMode(a, coren, n=1, m=2, kind=:HE, ϕ=0.0, loss=false)
)
nmodes = length(modes)

grid = Grid.EnvGrid(L, λ0, (160e-9, 3000e-9), 1e-12)

energyfun = Fields.energyfuncs(grid)[1]

responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),)

inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)

Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
                               modes, :y; full=false)

linop = LinearOps.make_linop(grid, modes, λ0)
statsfun = Stats.default(grid, Eω, modes, linop, transform; gas=gas, windows=((150e-9, 300e-9),))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)

Luna.run(Eω, grid, linop, transform, FT, output)

##
Plotting.pygui(true)
Plotting.stats(output)
Plotting.prop_2D(output)
Plotting.prop_2D(output, bandpass=(150e-9, 300e-9))
Plotting.time_1D(output, [5e-2, 9.8e-2])
Plotting.time_1D(output, [5e-2, 9e-2], modes=:sum, bandpass=(150e-9, 300e-9))
Plotting.spec_1D(output, [5e-2, 9.8e-2])
Plotting.spec_1D(output, [5e-2, 9.8e-2], modes=:sum)


text:
using JQL
using Base.Test

# Column name skins
include("DataView.jl")


text:
"""
    struct Domain

Represent a physical domain as a union of entities.
"""
struct Domain
    entities::Vector{AbstractEntity}
end
Domain(ω::AbstractEntity) = Domain([ω,])
Domain() = Domain(AbstractEntity[])

"""
    entities(Ω::Domain)

Return a vector of all elementary entities making up a domain.
"""
entities(Ω::Domain) = Ω.entities

"""
    geometric_dimension(Ω::Domain)

If all entities forming the domain have the same `geometric_dimension`, return
that value; otherwise throw an assertion error.
"""
function geometric_dimension(Ω::Domain)
    dims = [geometric_dimension(ent) for ent in entities(Ω)] |> unique!
    msg = "currently not able to handle a `Domain` with entities of different
           geometrical dimension. Got entities of dimensions ($dims)"
    @assert length(dims) == 1 msg
    return first(dims)
end

"""
    length(Ω:::Domain)

The length of a domain corresponds to the number of elementary entities that make it.
"""
Base.length(Ω::Domain) = length(entities(Ω))

"""
    skeleton(Ω::Domain)

Return all the boundaries of the domain, i.e. the domain's skeleton.
"""
skeleton(Ω::Domain) = union(Domain.(boundary.(entities(Ω)))...)

"""
    ===(Ω1::Domain,Ω2::Domain)

Two `Domain`s are equal if all their entities are equal (regardless of order).
"""
function Base.:(==)(Ω1::Domain, Ω2::Domain)
    return issetequal(entities(Ω1),entities(Ω2))
end

"""
    in(ω::ElementaryEntity,Ω::Domain)

Check whether an `ElementaryEntity` belongs to a `Domain`. Not that his only
does a *shallow* comparisson, meaning it only checks that `ω` is not one of the
entities in `entities(Ω)`; thus `ω` could still belong equal one of the
boundaries of an entity in `entities(Ω)`.
"""
Base.in(ω::ElementaryEntity, Ω::Domain) = in(ω, entities(Ω))

Base.getindex(Ω::Domain, i) = entities(Ω)[i]
Base.lastindex(Ω::Domain) = length(Ω)

"""
    iterate()

Iterating over domain means iterating over its entities.
"""
function Base.iterate(Ω::Domain, state=1)
    # Check if we are done iterating
    if state > length(Ω)
        return nothing
    end
    # Return (result, state)
    return (Ω[state], state+1)
end

"""
Check whether a Domain is empty.
"""
Base.isempty(Ω::Domain) = length(entities(Ω)) == 0

"""
Domain which is not in the intersection of domains Ω1 and Ω2.
"""
function Base.setdiff(Ω1::Domain, Ω2::Domain)
    Domain(setdiff(entities(Ω1), entities(Ω2)))
end

"""
    union(Ωs::Domain...)

Union of domains. This does not perform a union on the true geometric objects in
`entities(Ω)`, but simply on their identification through `(dim,tag)`.
"""
function Base.union(Ω1::Domain,Ωs::Domain...)
    Domain(Vector{AbstractEntity}(unique(vcat(entities(Ω1),entities.(Ωs)...))))
end
Base.union(Ω::Domain) = Domain(unique(entities(Ω)))

# NOTE: this version was replaced due to type piracy
# function Base.union(Ωs::Domain...)
#     Domain(Vector{ElementaryEntity}(unique(vcat(entities.(Ωs)...))))
# end

"""
    assertequaldim(Ω1::Domain,Ω2::Domain)

Check that two domains have same dimension.

If one of the domain (or both) are empty, the assertion is assumed to be true.
"""
function assertequaldim(Ω1::Domain, Ω2::Domain)
    if isempty(Ω1) || isempty(Ω2) return nothing end
    msg = "The dimension of the first domain should be equal to the dimension
    of the second domain."
    @assert geometric_dimension(Ω1) == geometric_dimension(Ω2) msg
end

"""
Intersection between domain Ω1 and domain Ω2.
"""
function Base.intersect(Ω1::Domain, Ω2::Domain)
    assertequaldim(Ω1, Ω2)
    Ωinter = Domain(intersect(entities(Ω1), entities(Ω2)))
    if isempty(Ωinter)
        if isempty(Ω1) || isempty(Ω2)
            return Domain()
        else
            return intersect(skeleton(Ω1), skeleton(Ω2))
        end
    else
        return Ωinter
    end
end

function Base.push!(Ω::Domain,ent::ElementaryEntity)
    push!(entities(Ω),ent)
end

"""
Determine whether every element of domain Ω1 is also in domain Ω2.
"""
function Base.issubset(Ω1::Domain, Ω2::Domain)
    assertequaldim(Ω1, Ω2)
    return issubset(entities(Ω1), entities(Ω2))
end

"""
Remove domain Ω1 from domain Ω2.
"""
function remove(Ω1::Domain, Ω2::Domain)
    assertequaldim(Ω1, Ω2)
    return Domain(setdiff(entities(Ω2), intersect(entities(Ω1), entities(Ω2))))
end

"""Return the internal boundaries inside a domain."""
function internal_boundary(Ω::Domain)
    Ω1 = Domain(Ω[1])
    γ  = Domain()
    for ω2 in Ω[2:end]
        Ω2 = Domain(ω2)
        γ1 = Domain(vcat(boundary.(entities(Ω1))...))
        γ2 = Domain(vcat(boundary.(entities(Ω2))...))
        γ = union(γ, intersect(γ1, γ2))
        Ω1 = union(Ω1, Ω2)
    end
    return γ
end

"""Return the external boundaries inside a domain."""
function external_boundary(Ω::Domain)
    return remove(internal_boundary(Ω), skeleton(Ω))
end

"""
    boundary(Ω)

Return a domain comprising the external boundary of Ω.

See also: [`external_boundary`](@ref)
"""
boundary(Ω::Domain) = external_boundary(Ω::Domain)

"""
Return all tags of the elementary entities in the domain `Ω` corresponding to the dimension `d`.
"""
function keys(Ω::Domain, d::Integer)
    if isempty(Ω)
        return Tuple{Int64,Int64}[]
    elseif d == geometric_dimension(Ω)
        return Vector{Tuple{Int64,Int64}}(vcat(key.(entities(Ω))))
    elseif d < geometric_dimension(Ω)
        return unique(keys(skeleton(Ω),d))
    else
        error("Asking for tags with dimension > dimension of domain")
    end
end
function keys(Ω::Domain)
    isempty(Ω) ? Tuple{Int64, Int64}[] : keys(Ω, geometric_dimension(Ω))
end

"""
Return all tags of the elementary entities in the domain `Ω` corresponding to the dimensions contained in `dims`.
"""
function keys(Ω::Domain, dims::Vector{T}) where T <: Integer
    tgs = Vector{Tuple{Int64, Int64}}(undef, 0)
    for d in dims push!(tgs, keys(Ω, d)...) end
    return tgs
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# This file is part of Kpax3. License is MIT.

include("joint/densities.jl")

include("likelihoods/log_likelihood.jl")
include("likelihoods/marginal_likelihood.jl")
include("likelihoods/merge.jl")
include("likelihoods/split.jl")
include("likelihoods/biased_random_walk.jl")

include("partitioncols/aminoacids/densities.jl")
include("partitioncols/aminoacids/simulate.jl")

include("partitionrows/ewenspitman/densities.jl")
include("partitionrows/ewenspitman/mcmc.jl")

include("loss_functions/binder.jl")


text:
# # Types
#
# Topics:
# 1. How types work in Julia
# 2. Defining types
# 3. Type parameters
# 4. Example: `OneHotVector`
# 5. Tuple types
# 6. Summary of kinds of types
# 7. The types of functions

#-

# ## How types work in Julia
#
# Every Julia value has a unique concrete type associated with it:

typeof(3)
typeof('a')

# Types are also first-class values:

typeof(Int)
typeof(DataType)

# The type of a value describes its memory layout as well as describing,
# semantically, "what kind of thing it is".

#-

# There are also *abstract types* that serve to group other types together.

supertype(Int)
supertype(Signed)
supertype(Integer)
supertype(Real)
supertype(Number)

# `isa` queries whether a value belongs to a type.
# `<:` queries whethe one type is a subtype (subset) of another

2 isa Int
2 isa Real
2 isa String
Int <: Integer

# `::` asserts that a value has a certain type, or else throws an error

2::Int
2::String

#-

# ## Defining types

struct Point
    x::Float64
    y::Float64
end

Point(1,2)

#-

# Types can be "called" and have methods added like functions.
# When used that way they are sometimes called *constructors*.

Point() = Point(0,0)

Point()

Int(2.0)

Int(2.1)

# ### Abstract types

abstract type PointLike end

struct Point2 <: PointLike
    x::Float64
    y::Float64
end

# ### Mutability
#
# Instances of `Point` cannot be modified after construction:

p = Point(1,2)
p.x = 0

# Use `mutable struct` to declare a mutable type:

mutable struct MutablePoint
    x  # note: default field type is `Any`
    y
end

# ### Primitive types
#
# Primitive types are "scalar"-like types whose representation is just a string
# of bits. They have no fields.
# These are not defined often in user code, but are needed in some cases.

primitive type MyInt32 <: Signed 32 end

# ## Type parameters
#
# Some types actually define a family of related types.
# The following is a concrete type, but the *parameters* in curly
# braces indicate that it is one of many possible related types:

Array{Float64,2}

# What is `Array` by itself?

typeof(Array)

# `Array` is an *iterated union* of dense, in-memory arrays over all
# element types and numbers of dimensions.

Array == Array{T,N} where N where T

# The curly braces substitute values for these variables:

Array{Int} == Array{Int,N} where N

Array{Int} <: Array

(Array{T,3} where T) <: Array

[1,2] isa Array{Int}

[1,2] isa Array{T,1} where T

[1,2] isa Array{T,2} where T

# `Vector` is an alias for 1-d Arrays:

Vector == Array{T,1} where T

#-

# Types with different parameters are just different, and have no subtype
# relationship. This is called *invariance*.

Array{Int} <: Array{Real}

# This is surprising at first, but makes sense if you think about memory
# representations.
# `Array{Real}` can hold any kind of real number so it needs to be an array of
# pointers (or some other structure containing type information for each
# element).
# `Array{Int}` only has `Int`s, so it can use an efficient "flat"
# representation.

# ### Defining types with parameters

struct GenericPoint{T<:Real}
    x::T
    y::T
end

GenericPoint(1,2)

GenericPoint(1,2.0)

GenericPoint(1.0,2.0)

GenericPoint{Float64}(1,2)

# Constructors can be defined for these too, using syntax that looks like
# how they are called:

function GenericPoint{T}() where T
    return GenericPoint{T}(0,0)
end

GenericPoint{Int}()

GenericPoint()

GenericPoint() = GenericPoint{Float64}()

GenericPoint()

# ## Example: `OneHotVector`
#
# Often used in machine learning, a "one hot" vector is a vector of all zeros,
# except for a single `1` entry.
# Representing it as a standard vector is memory-inefficient, so it cries out
# for a special implementation.

struct OneHotVector <: AbstractVector{Int}
    idx::Int
    len::Int
end

Base.size(v::OneHotVector) = (v.len,)

Base.getindex(v::OneHotVector, i::Integer) = Int(i == v.idx)

OneHotVector(3, 10)

A = rand(5,5)

A * OneHotVector(3, 5)

# ### Exercise: make it generic over element type

# ## Tuple types
#
# Tuples are very much like structs, but their types a bit special.
# 1. They have fields, but no names, only indices.
# 2. They are always immutable.
# 3. They can have an arbitrary number of parameters/fields.
# 4. They are *covariant*.

t = ("a", "b")
typeof(t)

# Field access is syntax for calls to `getfield` and `setfield!`

getfield(t, 2)

# For tuples `getindex` == `getfield`:

t[1]

getfield(p, :x)
getfield(p, 2)

Tuple{String} <: Tuple{Any}

# Remember: not true for `Vector{String} <: Vector{Any}`!

# ## Summary of kinds of types
#
# 1. Concrete types: struct, mutable struct, tuple, primitive
# 2. abstract types
# 3. UnionAll types, `Array{T} where T`
# 4. Union types, `Union{Int, Missing}`
# 5. The empty, or "bottom" type `Union{}`

# ## The types of functions
#
# How do functions fit into this scheme?
# Somewhat non-obviously, they are structs!
# In Julia, every "function" has its own struct type.

#-

f(x) = 2x

# Is equivalent to

struct typeof_f
end

(::typeof_f)(x) = 2x

const f = typeof_f()

f(3)

#-

# A simple, top-level function is a struct with 0 fields.

typeof(+)
typeof(cos)

# *Closures*, or inner functions, that capture surrounding
# variables, keep those variables as fields.

adder(x) = y->x+y

a = adder(100)

typeof(a)
dump(a)

a.x

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials