text:
#=
@muladd function perform_step!(integrator,cache::SRICache,f=integrator.f)
  @unpack c₀,c₁,A₀,A₁,B₀,B₁,α,β₁,β₂,β₃,β₄,stages,error_terms = cache.tab
  @unpack H0,H1,A0temp,A1temp,B0temp,B1temp,A0temp2,A1temp2,B0temp2,B1temp2,atemp,btemp,E₁,E₂,E₁temp,ftemp,gtemp,chi1,chi2,chi3,tmp = cache
  @unpack t,dt,uprev,u,W,p = integrator
  @.. chi1 = .5*(W.dW.^2 - dt)/integrator.sqdt #I_(1,1)/sqrt(h)
  @.. chi2 = .5*(W.dW + W.dZ/sqrt3) #I_(1,0)/h
  @.. chi3 = 1/6 * (W.dW.^3 - 3*W.dW*dt)/dt #I_(1,1,1)/h
  for i=1:stages
    fill!(H0[i],zero(eltype(integrator.u)))
    fill!(H1[i],zero(eltype(integrator.u)))
  end
  for i = 1:stages
    fill!(A0temp,zero(eltype(integrator.u)))
    fill!(B0temp,zero(eltype(integrator.u)))
    fill!(A1temp,zero(eltype(integrator.u)))
    fill!(B1temp,zero(eltype(integrator.u)))
    for j = 1:i-1
      integrator.f((t + c₀[j]*dt),H0[j],ftemp)
      integrator.g((t + c₁[j]*dt),H1[j],gtemp)
      @.. A0temp = A0temp + A₀[j,i] * ftemp
      @.. B0temp = B0temp + B₀[j,i] * gtemp
      @.. A1temp = A1temp + A₁[j,i] * ftemp
      @.. B1temp = B1temp + B₁[j,i] * gtemp
    end
    @.. H0[i] = uprev + A0temp*dt + B0temp*chi2
    @.. H1[i] = uprev + A1temp*dt + B1temp*integrator.sqdt
  end
  fill!(atemp,zero(eltype(integrator.u)))
  fill!(btemp,zero(eltype(integrator.u)))
  fill!(E₂,zero(eltype(integrator.u)))
  fill!(E₁temp,zero(eltype(integrator.u)))
  for i = 1:stages
    integrator.f((t+c₀[i]*dt),H0[i],ftemp)
    integrator.g((t+c₁[i]*dt),H1[i],gtemp)
    @.. atemp = atemp + α[i] * ftemp
    @.. btemp = btemp + (β₁[i] * W.dW + β₂[i] * chi1) * gtemp
    @.. E₂    = E₂    + (β₃[i] * chi2 + β₄[i] * chi3) * gtemp
    if i <= error_terms
      @.. E₁temp += ftemp
    end
  end

  @.. u = uprev + (dt*atemp + btemp) + E₂

  if integrator.opts.adaptive
    @.. E₁ = dt * E₁temp

    calculate_residuals!(tmp, E₁, E₂, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)
  end
  integrator.u = u
end
=#

@muladd function perform_step!(integrator,cache::SRICache,f=integrator.f)
  @unpack c₀,c₁,A₀,A₁,B₀,B₁,α,β₁,β₂,β₃,β₄,stages,error_terms = cache.tab
  @unpack H0,H1,A0temp,A1temp,B0temp,B1temp,A0temp2,A1temp2,B0temp2,B1temp2,atemp,btemp,E₁,E₂,E₁temp,ftemp,gtemp,chi1,chi2,chi3,tmp = cache
  @unpack t,dt,uprev,u,W,p = integrator

  sqrt3 = sqrt(3one(eltype(W.dW)))
  if typeof(W.dW) <: Union{SArray,Number}
    chi1 = (W.dW.^2 - dt)/2integrator.sqdt #I_(1,1)/sqrt(h)
    chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
    chi3 = (W.dW.^3 - 3W.dW*dt)/6dt #I_(1,1,1)/h
  else
    @.. chi1 = (W.dW.^2 - dt)/2integrator.sqdt #I_(1,1)/sqrt(h)
    @.. chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
    @.. chi3 = (W.dW.^3 - 3W.dW*dt)/6dt #I_(1,1,1)/h
  end

  for i=1:stages
    fill!(H0[i],zero(eltype(integrator.u)))
    fill!(H1[i],zero(eltype(integrator.u)))
  end
  for i = 1:stages
    fill!(A0temp,zero(eltype(integrator.u)))
    fill!(B0temp,zero(eltype(integrator.u)))
    fill!(A1temp,zero(eltype(integrator.u)))
    fill!(B1temp,zero(eltype(integrator.u)))
    for j = 1:i-1
      integrator.f(ftemp,H0[j],p,t + c₀[j]*dt)
      integrator.g(gtemp,H1[j],p,t + c₁[j]*dt)
      @.. A0temp = A0temp + A₀[j,i]*ftemp
      @.. B0temp = B0temp + B₀[j,i]*gtemp
      @.. A1temp = A1temp + A₁[j,i]*ftemp
      @.. B1temp = B1temp + B₁[j,i]*gtemp
    end
    @.. H0[i] = uprev + A0temp*dt + B0temp*chi2
    @.. H1[i] = uprev + A1temp*dt + B1temp*integrator.sqdt
  end
  fill!(atemp,zero(eltype(integrator.u)))
  fill!(btemp,zero(eltype(integrator.u)))
  fill!(E₂,zero(eltype(integrator.u)))
  fill!(E₁temp,zero(eltype(integrator.u)))
  for i = 1:stages
    integrator.f(ftemp,H0[i],p,t+c₀[i]*dt)
    integrator.g(gtemp,H1[i],p,t+c₁[i]*dt)
    @.. atemp = atemp + α[i]*ftemp
    @.. btemp = btemp + (β₁[i]*W.dW + β₂[i]*chi1)*gtemp
    @.. E₂    = E₂    + (β₃[i]*chi2 + β₄[i]*chi3)*gtemp
    if i <= error_terms
      @.. E₁temp += ftemp
    end
  end

  @.. u = uprev + (dt*atemp + btemp) + E₂

  if integrator.opts.adaptive
    @.. E₁ = dt * E₁temp

    calculate_residuals!(tmp, E₁, E₂, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)
  end
end

#=
@muladd function perform_step!(integrator,cache::SRIW1Cache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack chi1,chi2,chi3,fH01o4,g₁o2,H0,H11,H12,H13,g₂o3,Fg₂o3,g₃o3,Tg₃o3,mg₁,E₁,E₂,fH01,fH02,g₁,g₂,g₃,g₄,tmp = cache
  @.. chi1 = (W.dW.^2 - dt)/2integrator.sqdt #I_(1,1)/sqrt(h)
  @.. chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
  @.. chi3 = (W.dW.^3 - 3W.dW*dt)/6dt #I_(1,1,1)/h
  integrator.f(fH01,uprev,p,t)
  @.. fH01 = dt*fH01
  integrator.g(t,uprev,g₁)
  dto4 = dt/4
  @.. fH01o4 = fH01/4
  @.. g₁o2 = g₁/2
  @.. H0 =  uprev + 3 * (fH01o4 + chi2 * g₁o2)
  @.. H11 = uprev + fH01o4 + integrator.sqdt * g₁o2
  @.. H12 = uprev + fH01 - integrator.sqdt * g₁
  integrator.g(t+dto4,H11,g₂)
  integrator.g(t+dt,H12,g₃)
  @.. H13 = uprev + fH01o4 + integrator.sqdt * (-5 * g₁ + 3 * g₂ + g₃ / 2)

  integrator.g(t+dto4,H13,g₄)
  integrator.f(fH02,H0,p,t+3dto4)

  @.. fH02 = fH02*dt
  @.. g₂o3 = g₂/3
  @.. Fg₂o3 = 4g₂o3
  @.. g₃o3 = g₃/3
  @.. Tg₃o3 = 2g₃o3
  @.. mg₁ = -g₁
  @.. E₂ = chi2 * (2 * g₁ - Fg₂o3 - Tg₃o3) + chi3 * (2 * mg₁ + 5 * g₂o3 - Tg₃o3 + g₄)

  @.. u = uprev +  (fH01 + 2 * fH02) / 3 + W.dW * (mg₁ + Fg₂o3 + Tg₃o3) + chi1 * (mg₁ + Fg₂o3 - g₃o3) + E₂

  if integrator.opts.adaptive
    @.. E₁ = fH01 + fH02

    calculate_residuals!(tmp, E₁, E₂, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)
  end
  integrator.u = u
end
=#

@muladd function perform_step!(integrator,cache::SRIW1Cache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack chi1,chi2,chi3,fH01o4,g₁o2,H0,H11,H12,H13,g₂o3,Fg₂o3,g₃o3,Tg₃o3,mg₁,E₁,E₂,fH01,fH02,g₁,g₂,g₃,g₄,tmp = cache

  sqrt3 = sqrt(3one(eltype(W.dW)))
  if typeof(W.dW) <: Union{SArray,Number}
    chi1 = (W.dW.^2 - dt)/2integrator.sqdt #I_(1,1)/sqrt(h)
    chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
    chi3 = (W.dW.^3 - 3W.dW*dt)/6dt #I_(1,1,1)/h
  else
    @.. chi1 = (W.dW.^2 - dt)/2integrator.sqdt #I_(1,1)/sqrt(h)
    @.. chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
    @.. chi3 = (W.dW.^3 - 3W.dW*dt)/6dt #I_(1,1,1)/h
  end

  integrator.f(fH01,uprev,p,t)
  @.. fH01 = dt*fH01
  integrator.g(g₁,uprev,p,t)
  dto4 = dt/4

  @.. fH01o4 = fH01/4
  @.. g₁o2 = g₁/2
  @.. H0 =  uprev + 3*(fH01o4  + chi2*g₁o2)
  @.. H11 = uprev + fH01o4   + integrator.sqdt*g₁o2
  @.. H12 = uprev + fH01     - integrator.sqdt*g₁

  integrator.g(g₂,H11,p,t+dto4)
  integrator.g(g₃,H12,p,t+dt)

  @.. H13 = uprev + fH01o4 + integrator.sqdt*(-5g₁ + 3g₂ + g₃/2)

  integrator.g(g₄,H13,p,t+dto4)
  integrator.f(fH02,H0,p,t+3dto4)

  @.. fH02 = fH02*dt
  @.. g₂o3 = g₂/3
  @.. Fg₂o3 = 4g₂o3
  @.. g₃o3 = g₃/3
  @.. Tg₃o3 = 2g₃o3
  @.. mg₁ = -g₁

  @.. E₂ = chi2 * (2 * g₁ - Fg₂o3 - Tg₃o3) + chi3 * (2 * mg₁ + 5 * g₂o3 - Tg₃o3 + g₄)

  @.. u = uprev + (fH01 + 2 * fH02) / 3 + W.dW * (mg₁ + Fg₂o3 + Tg₃o3) + chi1 * (mg₁ + Fg₂o3 - g₃o3) + E₂

  if integrator.opts.adaptive
    @.. E₁ = fH01 + fH02

    calculate_residuals!(tmp, E₁, E₂, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)
  end
end

@muladd function perform_step!(integrator,cache::SRIW1ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  sqrt3 = sqrt(3one(eltype(W.dW)))
  chi1 = @.. (W.dW.^2 - dt)/2integrator.sqdt #I_(1,1)/sqrt(h)
  chi2 = @.. (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
  chi3 = @.. (W.dW.^3 - 3W.dW*dt)/6dt #I_(1,1,1)/h
  fH01 = dt*integrator.f(uprev,p,t)

  g₁ = integrator.g(uprev,p,t)
  fH01o4 = fH01/4
  dto4 = dt/4
  g₁o2 = g₁/2
  H0 =  @.. uprev + 3 * (fH01o4 + chi2 * g₁o2)
  H11 = @.. uprev + fH01o4 + integrator.sqdt * g₁o2
  H12 = @.. uprev + fH01 - integrator.sqdt * g₁
  g₂ = integrator.g(H11,p,t+dto4)
  g₃ = integrator.g(H12,p,t+dt)
  H13 = @.. uprev + fH01o4 + integrator.sqdt * (-5 * g₁ + 3 * g₂ + g₃ / 2)

  g₄ = integrator.g(H13,p,t+dto4)
  fH02 = dt*integrator.f(H0,p,t+3dto4)

  g₂o3 = g₂/3
  Fg₂o3 = 4g₂o3
  g₃o3 = g₃/3
  Tg₃o3 = 2g₃o3
  mg₁ = -g₁
  E₂ = @.. chi2 * (2 * g₁ - Fg₂o3 - Tg₃o3) + chi3 * (2 * mg₁ + 5 * g₂o3 - Tg₃o3 + g₄)

  u = uprev + (fH01 + 2fH02)/3 + W.dW.*(mg₁ + Fg₂o3 + Tg₃o3) + chi1.*(mg₁ + Fg₂o3 - g₃o3) + E₂
  if integrator.opts.adaptive
    E₁ = fH01 .+ fH02

    resids = calculate_residuals(E₁, E₂, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::SRIConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack c₀,c₁,A₀,A₁,B₀,B₁,α,β₁,β₂,β₃,β₄,stages,H0,H1,error_terms = cache
  sqrt3 = sqrt(3one(eltype(W.dW)))
  chi1 = .5*(W.dW.^2 - dt)/integrator.sqdt #I_(1,1)/sqrt(h)
  chi2 = .5*(W.dW + W.dZ/sqrt3) #I_(1,0)/h
  chi3 = 1/6 * (W.dW.^3 - 3*W.dW*dt)/dt #I_(1,1,1)/h

  fill!(H0,zero(typeof(u)))
  fill!(H1,zero(typeof(u)))
  @inbounds for i in 1:stages
    A0temp = zero(u)
    B0temp = zero(u)
    A1temp = zero(u)
    B1temp = zero(u)
    for j in 1:i-1
      ftmp = integrator.f(H0[j],p,t + c₀[j]*dt)
      A0temp = @.. A0temp + A₀[j,i] * ftmp
      A1temp = @.. A1temp + A₁[j,i] * ftmp

      gtmp = integrator.g(H1[j],p,t + c₁[j]*dt)
      B0temp = @.. B0temp + B₀[j,i] * gtmp
      B1temp = @.. B1temp + B₁[j,i] * gtmp
    end
    H0[i] = @.. uprev + dt * A0temp + chi2 * B0temp
    H1[i] = @.. uprev + dt * A1temp + integrator.sqdt * B1temp
  end
  atemp = zero(u)
  btemp = zero(u)
  E₂    = zero(u)
  E₁temp= zero(u)
  @inbounds for i in 1:stages
    ftmp = integrator.f(H0[i],p,t+c₀[i]*dt)
    atemp = @.. atemp + α[i] * ftmp

    gtmp = integrator.g(H1[i],p,t+c₁[i]*dt)
    btemp = @.. btemp + (β₁[i] * W.dW + β₂[i] * chi1) * gtmp
    E₂ = @.. E₂ + (β₃[i] * chi2 + β₄[i] * chi3) * gtmp
    if i <= error_terms #1 or 2
      E₁temp = E₁temp .+ ftmp
    end
  end

  u = (uprev + dt*atemp) + btemp + E₂

  if integrator.opts.adaptive
    E₁ = dt .* E₁temp

    resids = calculate_residuals(E₁, E₂, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::FourStageSRIConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack a021,a031,a032,a041,a042,a043,a121,a131,a132,a141,a142,a143,b021,b031,b032,b041,b042,b043,b121,b131,b132,b141,b142,b143,c02,c03,c04,c11,c12,c13,c14,α1,α2,α3,α4,beta11,beta12,beta13,beta14,beta21,beta22,beta23,beta24,beta31,beta32,beta33,beta34,beta41,beta42,beta43,beta44 = cache

  sqrt3 = sqrt(3one(eltype(W.dW)))
  chi1 = (W.dW.^2 .- dt)/(2integrator.sqdt) #I_(1,1)/sqrt(h)
  chi2 = (W.dW .+ W.dZ./sqrt3)./2 #I_(1,0)/h
  chi3 = (W.dW.^3 .- 3*W.dW*dt)/(6dt) #I_(1,1,1)/h

  k1 = integrator.f(uprev,p,t)
  g1 = integrator.g(uprev,p,t+c11*dt)

  H01 = uprev + dt*a021*k1 + b021*chi2.*g1
  H11 = uprev + dt*a121*k1 + integrator.sqdt*b121*g1

  k2 = integrator.f(H01,p,t+c02*dt)
  g2 = integrator.g(H11,p,t+c12*dt)

  H02 = uprev + dt*(a031*k1 + a032*k2) + chi2.*(b031*g1 + b032*g2)
  H12 = uprev + dt*(a131*k1 + a132*k2) + integrator.sqdt*(b131*g1 + b132*g2)

  k3 = integrator.f(H02,p,t+c03*dt)
  g3 = integrator.g(H12,p,t+c13*dt)

  H03 = uprev + dt*(a041*k1 + a042*k2 + a043*k3) + chi2.*(b041*g1 + b042*g2 + b043*g3)
  H13 = uprev + dt*(a141*k1 + a142*k2 + a143*k3) + integrator.sqdt*(b141*g1 + b142*g2 + b143*g3)

  k4 = integrator.f(H03,p,t+c04*dt)
  g4 = integrator.g(H13,p,t+c14*dt)

  E₂ = chi2.*(beta31*g1 + beta32*g2 + beta33*g3 + beta34*g4) + chi3.*(beta41*g1 + beta42*g2 + beta43*g3 + beta44*g4)

  u = uprev + dt*(α1*k1 + α2*k2 + α3*k3 + α4*k4) + E₂ + W.dW.*(beta11*g1 + beta12*g2 + beta13*g3 + beta14*g4) + chi1.*(beta21*g1 + beta22*g2 + beta23*g3 + beta24*g4)

  E₁ = dt*(k1 + k2 + k3 + k4)

  if typeof(integrator.alg) <: StochasticCompositeAlgorithm && typeof(integrator.alg.algs[1]) <: SOSRI2
    ϱu = integrator.opts.internalnorm(k4 - k3, t)
    ϱd = integrator.opts.internalnorm(H03 - H02, t)
    integrator.eigen_est = ϱu/ϱd
  end

  if integrator.opts.adaptive
    E₁ = @.. dt * (k1 + k2 + k3 + k4)

    resids = calculate_residuals(E₁, E₂, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::FourStageSRICache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack chi1,chi2,chi3,tab,g1,g2,g3,g4,k1,k2,k3,k4,E₁,E₂,tmp,H02,H03 = cache
  @unpack a021,a031,a032,a041,a042,a043,a121,a131,a132,a141,a142,a143,b021,b031,b032,b041,b042,b043,b121,b131,b132,b141,b142,b143,c02,c03,c04,c11,c12,c13,c14,α1,α2,α3,α4,beta11,beta12,beta13,beta14,beta21,beta22,beta23,beta24,beta31,beta32,beta33,beta34,beta41,beta42,beta43,beta44 = cache.tab

  sqdt = integrator.sqdt
  sqrt3 = sqrt(3one(eltype(W.dW)))
  if typeof(W.dW) <: Union{SArray,Number}
    chi1 = (W.dW.^2 - dt)/2integrator.sqdt #I_(1,1)/sqrt(h)
    chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
    chi3 = (W.dW.^3 - 3W.dW*dt)/6dt #I_(1,1,1)/h
  else
    @.. chi1 = (W.dW.^2 - dt)/2integrator.sqdt #I_(1,1)/sqrt(h)
    @.. chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
    @.. chi3 = (W.dW.^3 - 3W.dW*dt)/6dt #I_(1,1,1)/h
  end

  integrator.f(k1,uprev,p,t)
  integrator.g(g1,uprev,p,t+c11*dt)

  @.. tmp = uprev + dt*a021*k1 + chi2*b021*g1
  integrator.f(k2,tmp,p,t+c02*dt)

  @.. tmp = uprev + dt*a121*k1 + sqdt*b121*g1
  integrator.g(g2,tmp,p,t+c12*dt)

  @.. H02 = uprev + dt*(a031*k1 + a032*k2) + chi2*(b031*g1 + b032*g2)
  integrator.f(k3,H02,p,t+c03*dt)
  @.. tmp = uprev + dt*(a131*k1 + a132*k2) + sqdt*(b131*g1 + b132*g2)
  integrator.g(g3,tmp,p,t+c13*dt)

  @.. H03 = uprev + dt*(a041*k1 + a042*k2 + a043*k3) + chi2*(b041*g1 + b042*g2 + b043*g3)
  integrator.f(k4,H03,p,t+c04*dt)

  @.. tmp = uprev + dt*(a141*k1 + a142*k2 + a143*k3) + sqdt*(b141*g1 + b142*g2 + b143*g3)
  integrator.g(g4,tmp,p,t+c14*dt)

  if typeof(integrator.alg) <: StochasticCompositeAlgorithm && typeof(integrator.alg.algs[1]) <: SOSRI2
    @.. tmp = k4 - k3
    ϱu = integrator.opts.internalnorm(tmp, t)
    @.. tmp = H03 - H02
    ϱd = integrator.opts.internalnorm(tmp, t)
    integrator.eigen_est = ϱu/ϱd
  end

  @.. E₂ = chi2*(beta31*g1 + beta32*g2 + beta33*g3 + beta34*g4) + chi3*(beta41*g1 + beta42*g2 + beta43*g3 + beta44*g4)
  @.. u = uprev + dt*(α1*k1 + α2*k2 + α3*k3 + α4*k4) + E₂ + W.dW*(beta11*g1 + beta12*g2 + beta13*g3 + beta14*g4) + chi1*(beta21*g1 + beta22*g2 + beta23*g3 + beta24*g4)

  if integrator.opts.adaptive
    @.. E₁ = dt * (k1 + k2 + k3 + k4)

    calculate_residuals!(tmp, E₁, E₂, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)
  end
end


text:
@doc """ Compares the lengths and maximum fitnesses of random, greedy, and reluctant adaptive walks.
Outputs a table of means and standard deviations for these walks.
The parameters shown in the table are set in compare_walks_driver.jl.
See compare_walks_driver.jl for information on how to run on multiple cores.
"""
@everywhere immutable Job
  trial::Int64
  n::Int64
  k::Int64
  q::Int64  # q  is the number of integer random values for NKq landscapes (F in the paper)
            # q==0 means NK landscape, q > 0 means NKq landscape
  b::Int64
end

function run_jobs(jobs)
  results = pmap(run,jobs)
  all_results=accumulate_results(results)
  print_results(all_results)
end

@everywhere immutable JobResult
  job::Job
  rand_len::Int64
  rand_fit::Float64
  greedy_len::Int64
  greedy_fit::Float64
  reluct_len::Int64
  reluct_fit::Float64
  fiteq_len::Int64
  fiteq_fit::Float64
end

@everywhere immutable AllResults
  n::Int64
  k::Int64
  q::Int64
  rand_len::Array{Int64,1}
  rand_fit::Array{Float64,1}
  greedy_len::Array{Int64,1}
  greedy_fit::Array{Float64,1}
  reluct_len::Array{Int64,1}
  reluct_fit::Array{Float64,1}
  fiteq_len::Array{Int64,1}
  fiteq_fit::Array{Float64,1}
end

@everywhere function run(job::Job)
  if job.q == 0
    l = NKLandscape(job.n, job.k, near=job.b == 1)
  else
    l = NKqLandscape(job.n, job.k, job.q, near=job.b==1)
  end
  g = rand(Genotype, l)
  wrand = random_adaptive_walk(g, l)
  wgreedy = greedy_adaptive_walk(g, l)
  wreluct = reluctant_adaptive_walk(g, l)
  wfiteq = fitter_then_neutral_walk(g, l)
  res = JobResult(
    job,
    wrand.length,
    wrand.fitnesses[end],
    wgreedy.length,
    wgreedy.fitnesses[end],
    wreluct.length,
    wreluct.fitnesses[end],
    wfiteq.length,
    wfiteq.fitnesses[end])

  return res
end

@everywhere function print_header()
  println("Length of Paths to Optima and Fitness of Optima for random, greedy, reluctant, and fitter_or_equal walks")
  println(BOUNDARIES[B_VALUE], " Interactions")
  println(RUNS, " runs per line of output")
  println("rnd:  random")
  println("grd:  greedy")
  println("rlt:  reluctant")
  println("ftq:  fitter or equal")
  println("ln:   length")
  println("ft:   fitness")
  println("N\tK\tq\tμ rndln\tσ rndln\tμ rndft\tσ rndft\tμ grdln\tσ grdln\tμ grdft\tσ grdft\tμ rltln\tσ rltln\tμ rltft\tσ rltft\tμ ftqln\tσ ftqln\tμ ftqft\tσ ftqft")
end

function accumulate_results(result_list::Array{Any,1})
  len = length(result_list)
  all_results = AllResults(
    result_list[1].job.n,
    result_list[1].job.k,
    result_list[1].job.q,
    zeros(Int64,len),
    zeros(Float64,len),
    zeros(Int64,len),
    zeros(Float64,len),
    zeros(Int64,len),
    zeros(Float64,len),
    zeros(Int64,len),
    zeros(Float64,len)
  )
  i = 1
  for res in result_list
    all_results.rand_len[i] = res.rand_len
    all_results.rand_fit[i] = res.rand_fit
    all_results.greedy_len[i] = res.greedy_len
    all_results.greedy_fit[i] = res.greedy_fit
    all_results.reluct_len[i] = res.reluct_len
    all_results.reluct_fit[i] = res.reluct_fit
    all_results.fiteq_len[i] = res.fiteq_len
    all_results.fiteq_fit[i] = res.fiteq_fit
    i += 1
  end
  return all_results
end

@everywhere function print_results(a::AllResults)
    @printf("%d\t%d\t%d\t%.3f\t%.3f\t%.3f\t%.3f", a.n, a.k, a.q, mean(a.rand_len), std(a.rand_len), mean(a.rand_fit), std(a.rand_fit))
    @printf("\t%.3f\t%.3f\t%.3f\t%.3f", mean(a.greedy_len), std(a.greedy_len), mean(a.greedy_fit), std(a.greedy_fit))
    @printf("\t%.3f\t%.3f\t%.3f\t%.3f", mean(a.reluct_len), std(a.reluct_len), mean(a.reluct_fit), std(a.reluct_fit))
    @printf("\t%.3f\t%.3f\t%.3f\t%.3f\n", mean(a.fiteq_len), std(a.fiteq_len), mean(a.fiteq_fit), std(a.fiteq_fit))
end

ext: jl


text:
using Turing: @VarName, randcat
using Distributions: Normal
using Test, StatsFuns

i = 1
@test @VarName(s)[1:end-1] == (:s,())
@test @VarName(x[1,2][1+5][45][3][i])[1:end-1] == (:x,([1,2],[6],[45],[3],[1]))
@test StatsFuns.logistic(1.1) == 1.0 / (exp.(-1.1) + 1.0)
@test StatsFuns.logit(0.3) ≈ -0.8472978603872036 atol=1e-9
@test isnan.(StatsFuns.logit(1.0)) == false
@test isinf.(StatsFuns.logit(1.0)) == true
@test isnan.(StatsFuns.logit(0.)) == false
@test isinf.(StatsFuns.logit(0.)) == true
randcat([0.1, 0.9])

@test Turing.data(Dict(:y => 1, :x => 2, :z => 1), [:x, :y]) == (x = 2, y = 1)


text:
location="Belgium"
df = DataFrame(XLSX.readtable("/Users/shardy/Documents/GitHub/hv_offshore_topology_optimization/v1.0/layout/input_data/"*location*"/data.xlsx", "wind_data")...)
wnd="Norther"
cbl_data=get_220kV_cables()
cable_tables=Array{Array{cable,1},1}()
km=72.5
#ps=[50.0,100.0,150.0,200.0,250.0,300.0,350.0,400.0,450.0,500.0,550.0,600.0,650.0,700.0,750.0,800.0,850.0,900.0,950.0,1000.0,1050.0,1100.0,1150.0,1200.0,1250.0,1300.0,1350.0,1400.0,1450.0,1500.0]
#ps=[50.0,100.0,150.0,200.0,250.0,300.0,350.0,400.0,450.0,500.0,550.0,600.0,650.0,700.0]
#
ps=[1200.0]

osss=Vector{}()
cbls=Vector{}()
osss_costs=Vector{}()
cbls_costs=Vector{}()
ks=get_Cost_Data()
xfo_data=get_Xfo_Data()
xfo_on=transformer()
for p in ps
    oss=bus()
    xfo_off=transformer()
    xfo_off.relia=get_offshore_xfo_failure_data(xfo_off)
    xfo_off.wnd=wnd
    xfo_off.mva=p
    xfo_on=deepcopy(xfo_off)
    xfo_off=xfo_oss(xfo_off,ks,xfo_data)
    push!(oss.xfmrs,xfo_off)
    xfo_pcc(xfo_on,ks,xfo_data)

    plat=platform()
    plat.mva=p
    plat.wnd=wnd
    plat=cost_ac_platform(plat,ks)
    plat=adjust_base_ac_platform(plat,ks)
    push!(oss.plat,plat)
    push!(osss,oss)

cable()

    push!(cable_tables,solve_cable_interval((1,km,p),cbl_data,wnd,ks))
    cable_dictionary=make_dictionaries_cables([cable_tables])
    sg_qc=cable_dictionary["220.0"][string(p)][length(cable_dictionary["220.0"][string(p)])].costs.sg+cable_dictionary["220.0"][string(p)][length(cable_dictionary["220.0"][string(p)])].costs.qc
    cbl_cost=cable_dictionary["220.0"][string(p)][length(cable_dictionary["220.0"][string(p)])].costs.ttl-sg_qc-cable_dictionary["220.0"][string(p)][length(cable_dictionary["220.0"][string(p)])].costs.rlc-cable_dictionary["220.0"][string(p)][length(cable_dictionary["220.0"][string(p)])].costs.eens-cable_dictionary["220.0"][string(p)][length(cable_dictionary["220.0"][string(p)])].costs.cm
    push!(cbls,cable_dictionary["220.0"][string(p)][length(cable_dictionary["220.0"][string(p)])])
    println(sg_qc)
    println(sg_qc)
    println(sg_qc)
    push!(cbls_costs,cbl_cost/km)
    push!(osss_costs,(plat.costs.ttl-plat.costs.cm+xfo_off.costs.ttl-xfo_off.costs.cm-xfo_off.costs.eens-xfo_off.costs.tlc+sg_qc/2))
    push!(osss_costs,(plat.costs.ttl-plat.costs.cm+sg_qc/2))
    #push!(osss_costs,xfo_on.costs.ttl+sg_qc/2)
end

cbls_costs[1]*km

println(osss_costs)

plotly()
p=plot()
plot!(p,[ps],[cbls_costs],xticks = 0:200:800,xlims=(0,800),ylims=(0,3),color = :black,label="",xaxis = ("MW", font(20, "Courier")),yaxis = ("ME", font(20, "Courier")))
gui()=#

plotly()
p=plot()
plot!(p,[ps],[osss_costs],xticks = 0:100:800,xlims=(0,700),ylims=(0,250),color = :black,label="",xaxis = ("MW", font(20, "Courier")),yaxis = ("ME", font(20, "Courier")))
gui()

osss_costs_orig=deepcopy(osss_costs)

GMIN=minimum([wf.mva for wf in owpps])
A=range(1,stop=30)
base_b=zeros(Int8,1,length(owpps))
B=Array{Tuple{Array{Int8,2},Int64},1}()
combos=collect(Iterators.flatten([combinations(A, k) for k = 1:9]))


text:
abstract type AbstractDisjunctiveFormulation end

struct NaiveBigM <: AbstractDisjunctiveFormulation
    activity_method::AbstractActivityMethod
end

struct NaiveBigMState
    z_vis::Vector{VI}
    sum_ci::CI{SAF,ET}
    lt_cis::Matrix{Union{Nothing,CI{SAF,LT}}}
    gt_cis::Matrix{Union{Nothing,CI{SAF,GT}}}
    proven_infeasible::Bool
end

function formulate!(
    model::MOI.ModelLike,
    method::NaiveBigM,
    disjunction::Disjunction,
)
    z_v, z_c = MOI.add_constrained_variables(
        model,
        [ZO() for i in 1:num_alternatives(disjunction)],
    )
    return formulate!(model, method, disjunction, z_v)
end

function formulate!(
    model::MOI.ModelLike,
    method::NaiveBigM,
    disjunction::Disjunction,
    z_vis::Vector{VI},
)
    # We need to copy the model (do we?) because we are intermixing adding
    # constraints and optimizing over the model.
    # TODO: Cache this activity_model inside the method for resuse.
    activity_model = MOIU.Model{Float64}()
    MOI.copy_to(activity_model, model, copy_names = false)
    m = MOI.output_dimension(disjunction.f)
    @assert m ==
            MOI.dimension(disjunction.s) ==
            size(disjunction.s.lbs, 1) ==
            size(disjunction.s.ubs, 1)
    d = length(z_vis)
    @assert d ==
            num_alternatives(disjunction) ==
            size(disjunction.s.lbs, 2) ==
            size(disjunction.s.ubs, 2)

    sum_ci = MOI.add_constraint(
        model,
        SAF([SAT(1.0, vi) for vi in z_vis], 0.0),
        MOI.EqualTo(1.0),
    )

    f_scalar = MOIU.scalarize(disjunction.f)

    lt_cis = similar(disjunction.s.lbs, Union{Nothing,CI{SAF,LT}})
    gt_cis = similar(disjunction.s.ubs, Union{Nothing,CI{SAF,GT}})
    proven_infeasible = false
    for i in 1:m
        f = f_scalar[i]
        for j in 1:d
            let lb = disjunction.s.lbs[i, j]
                # TODO: Can potentially run activity method to tighten lb
                gt_cis[i, j] = (
                    if lb == Inf
                        throw(
                            ErrorException(
                                "The disjunction is trivially infeasible.",
                            ),
                        )
                    elseif lb == -Inf
                        # Do nothing
                        nothing
                    else
                        m_val = minimum_activity(
                            method.activity_method,
                            activity_model,
                            f,
                        )
                        if m_val == -Inf
                            throw(
                                ErrorException(
                                    "Encountered infinite big-M coefficient, cannot formulate.",
                                ),
                            )
                        elseif m_val == Inf
                            proven_infeasible = true
                            nothing
                        else
                            MOIU.normalize_and_add_constraint(
                                model,
                                f + (lb - m_val) * (1.0 - SV(z_vis[j])),
                                GT(lb),
                            )
                        end
                    end
                )
            end

            let ub = disjunction.s.ubs[i, j]
                # TODO: Can potentially run activity method to tighten ub
                lt_cis[i, j] = (
                    if ub == -Inf
                        throw(
                            ErrorException(
                                "The disjunction is trivially infeasible.",
                            ),
                        )
                    elseif ub == Inf
                        # Do nothing
                        nothing
                    else
                        m_val = maximum_activity(
                            method.activity_method,
                            activity_model,
                            f,
                        )
                        if m_val == Inf
                            throw(
                                ErrorException(
                                    "Encountered infinite big-M coefficient, cannot formulate.",
                                ),
                            )
                        elseif m_val == -Inf
                            proven_infeasible = true
                            nothing
                        else
                            MOIU.normalize_and_add_constraint(
                                model,
                                f + (ub - m_val) * (1.0 - SV(z_vis[j])),
                                LT(ub),
                            )
                        end
                    end
                )
            end
        end
    end
    return NaiveBigMState(z_vis, sum_ci, lt_cis, gt_cis, proven_infeasible)
end


text:
#############################################
# BBBackpropTraceState (for backprop_trace) #
#############################################

struct SelectedLeafNode
    addr::Symbol
    value_ref::Expr
    gradient_var::Symbol
end

struct InternalNode
    addr::Symbol
    values_var::Symbol
    gradients_var::Symbol
end

struct BBBackpropTraceState
    generator_type::Type
    ir::BasicBlockIR
    trace::Symbol
    stmts::Vector{Expr}
    schema::Union{StaticAddressSchema,EmptyAddressSchema}
    value_refs::Dict{ValueNode,Expr}
    grad_vars::Dict{ValueNode,Symbol}
    leaf_nodes::Set{SelectedLeafNode}
    internal_nodes::Set{InternalNode}
end

function BBBackpropTraceState(generator_type, trace, stmts, schema, value_refs, grad_vars)
    leaf_nodes = Set{SelectedLeafNode}()
    internal_nodes = Set{InternalNode}()
    ir = get_ir(generator_type)
    BBBackpropTraceState(generator_type, ir, trace, stmts, schema, value_refs, grad_vars,
        leaf_nodes, internal_nodes)
end

###############################################
# BBBackpropParamsState (for backprop_params) #
###############################################

struct BBBackpropParamsState
    generator_type::Type
    ir::BasicBlockIR
    trace::Symbol
    stmts::Vector{Expr}
    value_refs::Dict{ValueNode,Expr}
    grad_vars::Dict{ValueNode,Symbol}
end

function BBBackpropParamsState(generator_type, trace, stmts, value_refs, grad_vars)
    ir = get_ir(generator_type)
    BBBackpropParamsState(generator_type, ir, trace, stmts, value_refs, grad_vars)
end

#####################################################
# code shared by backprop_trace and backprop_params #
#####################################################

function process!(state::Union{BBBackpropParamsState,BBBackpropTraceState}, node::JuliaNode)

    # if the output node is not floating point, don't do anything
    if !haskey(state.grad_vars, node.output)
        return # no statements
    end

    # call gradient function
    input_value_refs = [state.value_refs[in_node] for in_node in node.input_nodes]
    output_value_ref = state.value_refs[node.output]
    input_grad_increments = [gensym("incr") for _ in node.input_nodes]
    output_grad_var = state.grad_vars[node.output]
    gradient_fn = get_grad_fn(state.generator_type, node)
    push!(state.stmts, quote
        ($(input_grad_increments...),) = $gradient_fn(
            $output_grad_var, $output_value_ref, $(input_value_refs...))
    end)

    # increment input gradients
    inputs_do_ad = map((in_node) -> haskey(state.grad_vars, in_node), node.input_nodes)
    for (in_node, incr, do_ad) in zip(node.input_nodes, input_grad_increments, inputs_do_ad)
        if do_ad
            grad_var = state.grad_vars[in_node]
            push!(state.stmts, quote
                $grad_var += $incr
            end)
        end
    end
end

function process!(state::Union{BBBackpropParamsState,BBBackpropTraceState}, node::ArgsChangeNode)
    # skip
end

function process!(state::Union{BBBackpropParamsState,BBBackpropTraceState}, node::AddrChangeNode)
    # skip
end

function initialize_backprop!(ir::BasicBlockIR, stmts::Vector{Expr})

    # create a gradient variable for each value node, initialize them to zero.
    value_refs = Dict{ValueNode,Expr}()
    grad_vars = Dict{ValueNode,Symbol}()
    for (name, node) in ir.value_nodes
        value_refs[node] = value_trace_ref(:trace, node)
        grad_var = gensym("grad_$name")
        if node !== ir.output_node || !ir.output_ad
            if is_differentiable(get_type(node))
                grad_vars[node] = grad_var
                push!(stmts, quote
                    $grad_var = zero($(value_refs[node]))
                end)
            end
        end
    end

    # initialize the gradient variable for the output node
    if ir.output_node !== nothing && ir.output_ad
        grad_var = grad_vars[something(ir.output_node)]
        push!(stmts, quote
            $grad_var = retval_grad
        end)
    end

    (value_refs, grad_vars)
end

function input_gradients(ir::BasicBlockIR, grad_vars)
    input_grads_var = gensym("input_grads")
    input_grads = []
    for (node, has_grad) in zip(ir.arg_nodes, ir.args_ad)
        if has_grad
            push!(input_grads, grad_vars[node])
        else
            push!(input_grads, QuoteNode(nothing))
        end
    end
    Expr(:tuple, input_grads...)
end

function increment_input_gradients!(stmts, node, dist_or_gen, grad_vars, increments)
    inputs_do_ad = has_argument_grads(dist_or_gen)
    for (in_node, do_ad, incr) in zip(node.input_nodes, inputs_do_ad, increments)
        if do_ad
            if !haskey(grad_vars, in_node)
                error("$(dist_or_gen) has AD for an input that is not floating point, node: $node")
            end
            grad_var = grad_vars[in_node]
            push!(stmts, quote
                $grad_var += $incr
            end)
        end
    end
end

function increment_output_gradient!(stmts, node::AddrDistNode, grad_vars, increment)
    # NOTE: if the output has a gradient, then it must be a float...
    # but it may be a float and the may not hae a gradient, currently this is
    # silent; could warn?
    if has_output_grad(node.dist)
        if !haskey(grad_vars, node.output)
            error("Distribution $(node.dist) has AD but the return value is not floating point, node: $node")
        end
        output_grad_var = grad_vars[node.output]
        push!(stmts, quote
            $output_grad_var += $increment
        end)
    end
end

###################
# backprop_params #
###################

function process!(state::BBBackpropParamsState, node::AddrDistNode)

    # get gradient of log density with respect to output and inputs
    input_value_refs = [state.value_refs[in_node] for in_node in node.input_nodes]
    output_value_ref = Expr(:(.), :trace, QuoteNode(node.address))# state.value_refs[node.output]
    output_grad_incr = gensym("incr")
    input_grad_incrs = [gensym("incr") for _ in node.input_nodes]
    push!(state.stmts, quote
        ($output_grad_incr, $(input_grad_incrs...),) = Gen.logpdf_grad(
            $(QuoteNode(node.dist)), $output_value_ref, $(input_value_refs...))
    end)

    increment_output_gradient!(state.stmts, node, state.grad_vars, output_grad_incr)
    increment_input_gradients!(state.stmts, node, node.dist, state.grad_vars, input_grad_incrs)
end

function process!(state::BBBackpropParamsState, node::AddrGeneratorNode)

    # get gradients from generator
    output_do_ad = accepts_output_grad(node.gen)
    output_grad = output_do_ad ? state.grad_vars[node.output] : QuoteNode(nothing)
    addr = node.address
    input_grad_incrs = [gensym("incr") for _ in node.input_nodes]
    subtrace = Expr(:(.), :trace, QuoteNode(addr))
    push!(state.stmts, quote
        ($(input_grad_incrs...),) = backprop_params($(QuoteNode(node.gen)), $subtrace, $output_grad)
    end)

    increment_input_gradients!(state.stmts, node, node.gen, state.grad_vars, input_grad_incrs)
end

function codegen_backprop_params(gen::Type{T}, trace, retval_grad) where {T <: BasicGenFunction}
    ir = get_ir(gen)
    stmts = Expr[]

    # create a gradient variable for each value node, initialize them to zero.
    # also get trace references for each value node
    (value_refs, grad_vars) = initialize_backprop!(ir, stmts)

    # visit statements in reverse topological order
    state = BBBackpropParamsState(gen, :trace, stmts, value_refs, grad_vars)
    for node in reverse(ir.expr_nodes_sorted)
        process!(state, node)
    end

    # increment gradient accumulators for parameters
    for param in ir.params
        value_node = ir.value_nodes[param.name]
        grad_var = grad_vars[value_node]
        push!(stmts, quote
            gen.params_grad[$(QuoteNode(param.name))] += $grad_var
        end)
    end

    # return statement
    push!(stmts, quote
        return $(input_gradients(ir, grad_vars))
    end)
    Expr(:block, stmts...)
end

push!(Gen.generated_functions, quote
@generated function Gen.backprop_params(gen::Gen.BasicGenFunction, trace, retval_grad)
    Gen.codegen_backprop_params(gen, trace, retval_grad)
end
end)

##################
# backprop_trace #
##################

function process!(state::BBBackpropTraceState, node::AddrDistNode)

    # get gradient of log density with respect to output and inputs
    input_value_refs = [state.value_refs[in_node] for in_node in node.input_nodes]
    output_value_ref = state.value_refs[node.output]
    output_grad_incr = gensym("incr")
    input_grad_incrs = [gensym("incr") for _ in node.input_nodes]
    push!(state.stmts, quote
        ($output_grad_incr, $(input_grad_incrs...),) = Gen.logpdf_grad(
            $(QuoteNode(node.dist)), $output_value_ref, $(input_value_refs...))
    end)

    increment_output_gradient!(state.stmts, node, state.grad_vars, output_grad_incr)
    increment_input_gradients!(state.stmts, node, node.dist, state.grad_vars, input_grad_incrs)

    # handle selected address
    addr = node.address
    if isa(state.schema, StaticAddressSchema) && addr in leaf_node_keys(state.schema)
        if !haskey(state.grad_vars, node.output)
            error("Selected a random choice that is not floating point for gradient: $node")
        end
        output_grad_var = state.grad_vars[node.output]
        push!(state.leaf_nodes, SelectedLeafNode(addr, Expr(:(.), :trace, QuoteNode(addr)), output_grad_var))
    end
end

function process!(state::BBBackpropTraceState, node::AddrGeneratorNode)

    # get gradients from generator and handle selection
    output_do_ad = accepts_output_grad(node.gen)
    output_grad = output_do_ad ? state.grad_vars[node.output] : QuoteNode(nothing)
    addr = node.address
    input_grad_incrs = [gensym("incr") for _ in node.input_nodes]
    subtrace = Expr(:(.), :trace, QuoteNode(addr))
    has_selection = isa(state.schema, StaticAddressSchema) && addr in internal_node_keys(state.schema)
    if has_selection
        selection = Expr(:call, :static_get_internal_node, :selection, QuoteNode(Val(addr)))
        value_trie = gensym("value_trie_$addr")
        gradient_trie = gensym("gradient_trie_$addr")
        push!(state.internal_nodes, InternalNode(addr, value_trie, gradient_trie))
    else
        selection = quote Gen.EmptyAddressSet() end
        value_trie = :_
        gradient_trie = :_
    end
    push!(state.stmts, quote
        (($(input_grad_incrs...),), $value_trie, $gradient_trie) = backprop_trace(
            $(QuoteNode(node.gen)), $subtrace, $selection, $output_grad)
    end)

    increment_input_gradients!(state.stmts, node, node.gen, state.grad_vars, input_grad_incrs)
end

const backprop_values_trie = gensym("values")
const backprop_gradients_trie = gensym("gradients")

function choice_trie_construction(leaf_nodes_set, internal_nodes_set)
    leaf_nodes = collect(leaf_nodes_set)
    quoted_leaf_keys = map((node) -> QuoteNode(node.addr), leaf_nodes)
    leaf_values = map((node) -> node.value_ref, leaf_nodes)
    leaf_gradients = map((node) -> node.gradient_var, leaf_nodes)
    internal_nodes = collect(internal_nodes_set)
    quoted_internal_keys = map((node) -> QuoteNode(node.addr), internal_nodes)
    internal_values = map((node) -> node.values_var, internal_nodes)
    internal_gradients = map((node) -> node.gradients_var, internal_nodes)
    quote
        $backprop_values_trie = StaticAssignment(
            NamedTuple{($(quoted_leaf_keys...),)}(($(leaf_values...),)),
            NamedTuple{($(quoted_internal_keys...),)}(($(internal_values...),)))
        $backprop_gradients_trie = StaticAssignment(
            NamedTuple{($(quoted_leaf_keys...),)}(($(leaf_gradients...),)),
            NamedTuple{($(quoted_internal_keys...),)}(($(internal_gradients...),)))
    end
end

function codegen_backprop_trace(gen::Type{T}, trace, selection, retval_grad) where {T <: BasicGenFunction}
    schema = get_address_schema(selection)
    ir = get_ir(gen)
    stmts = Expr[]

    # create a gradient variable for each value node, initialize them to zero.
    # also get trace references for each value node
    (value_refs, grad_vars) = initialize_backprop!(ir, stmts)

    # visit statements in reverse topological order
    state = BBBackpropTraceState(gen, :trace, stmts, schema, value_refs, grad_vars)
    for node in reverse(ir.expr_nodes_sorted)
        process!(state, node)
    end

    # construct values and gradients static choice tries
    push!(stmts, choice_trie_construction(state.leaf_nodes, state.internal_nodes))

    # return statement
    push!(stmts, quote
        return ($(input_gradients(ir, grad_vars)), $backprop_values_trie, $backprop_gradients_trie)
    end)
    Expr(:block, stmts...)
end

push!(Gen.generated_functions, quote
@generated function Gen.backprop_trace(gen::Gen.BasicGenFunction{T,U}, trace::U, selection, retval_grad) where {T,U}
    schema = get_address_schema(selection)
    if !(isa(schema, StaticAddressSchema) || isa(schema, EmptyAddressSchema))
        # try to convert it to a static address set
        return quote backprop_trace(gen, trace, StaticAddressSet(selection), retval_grad) end
    end
    Gen.codegen_backprop_trace(gen, trace, selection, retval_grad)
end
end)

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
module MaxHelpingHandHintsFlagController

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/HintsFlagController" HintsFlagController(x::Integer, y::Integer, outputFlag::String="hints", not::Bool=false)

const placements = Ahorn.PlacementDict(
    "Hints Flag Controller (max480's Helping Hand)" => Ahorn.EntityPlacement(
        HintsFlagController
    )
)

function Ahorn.selection(entity::HintsFlagController)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

function Ahorn.editingOptions(effect::HintsFlagController)
    return Dict{String, Any}(
    )
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::HintsFlagController, room::Maple.Room)
    if entity.not
        Ahorn.drawSprite(ctx, "ahorn/MaxHelpingHand/hints_flag_controller_inv", 0, 0)
    else
        Ahorn.drawSprite(ctx, "ahorn/MaxHelpingHand/hints_flag_controller", 0, 0)
    end
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
@system Root(Organ) begin
end


text:
"""
    Defines the parameters used in the smoothing process.
"""
mutable struct StructuralParameters
    interest_yield::Float64
    intercept::Bool
    intercept_varying_maturity::Bool
    seasonality_trading_date::Bool
    seasonality_delivery_date::Bool
    run_lasso::Bool
    minimum_maturity::Int64
    λ::Float64
end

"""
    ContractInfo

Defines a negotiated forward contract.
"""
mutable struct ContractInfo
    trading_day::DateTime
    time_to_maturity::Int64 # H
    duration::Int64 # ΔT
    price::Float64 # F
end

"""
    ContractsByDay

Defines a vector of contracts traded in a trading day.
"""
mutable struct ContractsByDay
    trading_day::DateTime
    contracts::Vector{ContractInfo}
end

function number_of_contracts(vec_contracts_by_day::Vector{ContractsByDay})
    return number_of_contracts(vec_contracts_by_day[1])
end

function number_of_contracts(contracts_by_day::ContractsByDay)
    return length(contracts_by_day.contracts)
end

function date_of_trades(vec_contracts_by_day::Vector{ContractsByDay})
    dates = Vector{DateTime}(undef, length(vec_contracts_by_day))

    for (t, contracts_by_day) in enumerate(vec_contracts_by_day)
        dates[t] = contracts_by_day.trading_day
    end

    return dates
end

"""
    OptimizationInputs

parameters = Dict(
    "interest_yield" = 0.0
)
"""
mutable struct StructuralModelInputs
    parameters::Dict{String,Any}
    vec_contracts_by_day::Vector{ContractsByDay}
    X::Array{Float64,3}

    function StructuralModelInputs(parameters::Dict{String,Any},
                                   vec_contracts_by_day::Vector{ContractsByDay},
                                   X::Array{Float64,3})

        dates = date_of_trades(vec_contracts_by_day)

        T, J, L = size(X)
        if length(vec_contracts_by_day) != T
            error("X and vec_contracts_by_day must have same days")
        end

        return new(parameters, vec_contracts_by_day, X)
    end
end

function StructuralModelInputs(parameters::Dict{String,Any},
                                vec_contracts_by_day::Vector{ContractsByDay},
                                X::Nothing)

    max_maturity = maximum_maturity_from_vec_contracts_by_day(vec_contracts_by_day)
    max_maturity = max(max_maturity, parameters["minimum_maturity"])
    T = length(vec_contracts_by_day)

    return return StructuralModelInputs(parameters, vec_contracts_by_day, zeros(T, max_maturity, 0))
end

mutable struct ThirdLevelSolveTime
    maturity_smooth::Float64
    time_smooth::Float64
    estimate_structural_error::Float64
end

mutable struct StructuralSolveTime
    ground_level::Float64
    first_level::Float64
    second_level::Float64
    third_level::ThirdLevelSolveTime
end

mutable struct StructuralModelOutputs
    glm_output::Any
    args_structural_expression::Dict
    daily_minimum_arbitrage::Vector{Float64}
    maturity_weigth::Float64
    time_weigth::Float64
    structural_error_obj_function::Float64
    structural_error::Matrix{Float64}
    structural_dinamic::Matrix{Float64}
    reconstructed_forward_prices::Vector{ContractsByDay}
    solve_time::StructuralSolveTime
    adalasso_model
    lasso_selected_variables_idx::Vector{Int64}
    minimum_arbitrage_by_contract::Vector
end

"""
    read_ContractsByDay(df::DataFrame)

Converts the input data frame to the adequate format used in the input struct.
"""
function read_ContractsByDay(df::DataFrame)
    vec_contracts_by_day = ContractsByDay[]
    for row in eachrow(df)
        contracs_by_day = get_contracts_by_day(row)
        push!(vec_contracts_by_day, contracs_by_day)
    end
    return vec_contracts_by_day
end

function get_contracts_by_day(row::DataFrameRow)
    trading_day = row[1]
    num_contracts = Int(length(row[2:end]) / 3)
    vec_contract_info = ContractInfo[]

    for contract = 1:num_contracts
        price                   = row[2 + 3 * (contract - 1)]
        start_delivering_period = row[3 + 3 * (contract - 1)]
        end_delivering_period   = row[4 + 3 * (contract - 1)]
        contract_info = get_contract_info(trading_day, start_delivering_period, end_delivering_period, price)
        push!(vec_contract_info, contract_info)
    end

    return ContractsByDay(trading_day, vec_contract_info)
end

contract_price(price) = ismissing(price) ? NaN : price
duration_between_dates_by_day(dt1::DateTime, dt2::DateTime) = Dates.days(convert(Dates.Day, (dt2 - dt1)))

function get_contract_info(trading_day::DateTime,
                        start_delivering_period::DateTime,
                        end_delivering_period::DateTime,
                        price)
    time_to_maturity = duration_between_dates_by_day(trading_day, start_delivering_period)
    duration         = duration_between_dates_by_day(start_delivering_period, end_delivering_period)
    price            = contract_price(price)

    return ContractInfo(trading_day,
                    time_to_maturity,
                    duration,
                    price)
end

max_stars_repo_name: marinadietze/ForwardCurveSmoother

max_issues_repo_name: marinadietze/ForwardCurveSmoother

max_forks_repo_name: marinadietze/ForwardCurveSmoother


text:
abstract type AbstractSQLHandle end

#===================================================================================================
    <SQLEnvironement>
===================================================================================================#
struct SQLEnvironement <: AbstractSQLHandle
    ptr::Ptr{Void}
end

handletype(::Type{SQLEnvironement}) = SQL_HANDLE_ENV
#===================================================================================================
    </SQLEnvironement>
===================================================================================================#

#===================================================================================================
    <SQLConnection>
===================================================================================================#
struct SQLConnection <: AbstractSQLHandle
    ptr::Ptr{Void}
end

handletype(::Type{SQLConnection}) = SQL_HANDLE_DBC
#===================================================================================================
    <SQLConnection>
===================================================================================================#

#===================================================================================================
    <SQLStatement>
===================================================================================================#
struct SQLStatement <: AbstractSQLHandle
    ptr::Ptr{Void}
end

handletype(::Type{SQLStatement}) = SQL_HANDLE_STMT

# TODO not sure whether this will stay, this is a higher level wrapper of SQLStatement
type Statement
    dsn::DSN
    stmt::Ptr{Void}
    query::String
    task::Task
end
#===================================================================================================
    </SQLStatement>
===================================================================================================#

#===================================================================================================
    <SQLDescriptor>
===================================================================================================#
struct SQLDescriptor <: AbstractSQLHandle
    ptr::Ptr{Void}
end

handletype(::Type{SQLDescriptor}) = SQL_HANDLE_DESC
#===================================================================================================
    </SQLDescriptor>
===================================================================================================#

# TODO can't find any information on this...
#===================================================================================================
    <SQLWindow>
===================================================================================================#
struct SQLWindow <: AbstractSQLHandle
    ptr::Ptr{Void}
end
#===================================================================================================
    </SQLWindow>
===================================================================================================#

const HANDLETYPE_DICT = Dict(SQL_HANDLE_ENV=>SQLEnvironement,
                             SQL_HANDLE_DBC=>SQLConnection,
                             SQL_HANDLE_STMT=>SQLStatement,
                             SQL_HANDLE_DESC=>SQLDescriptor)

# this is a low level function for constructing handles
makehandle(htype::Integer, ptr::Ptr{Void}) = HANDLETYPE_DICT[Int16(htype)](ptr)
makehandle{T<:AbstractSQLHandle}(::Type{T}, ptr::Ptr{Void}) = T(ptr)

# "Allocate ODBC handles for interacting with the ODBC Driver Manager"
function ODBCAllocHandle(handletype, parenthandle)
    handle = Ref{Ptr{Void}}()
    ODBC.API.SQLAllocHandle(handletype, parenthandle, handle)
    handle = handle[]
    if handletype == ODBC.API.SQL_HANDLE_ENV
        ODBC.API.SQLSetEnvAttr(handle, ODBC.API.SQL_ATTR_ODBC_VERSION, ODBC.API.SQL_OV_ODBC3)
    end
    return handle
end

#===================================================================================================
    <block>
===================================================================================================#
"just a block of memory; T is the element type, `len` is total # of **bytes** pointed to, and `elsize` is size of each element"
type Block{T}
    ptr::Ptr{T}    # pointer to a block of memory
    len::Int       # total # of bytes in block
    elsize::Int    # size between elements in bytes
end

"""
Block allocator:
    -Takes an element type, and number of elements to allocate in a linear block
    -Optionally specify an extra dimension of elements that make up each element (i.e. container types)
"""
function Block{T}(::Type{T}, elements::Int, extradim::Integer=1)
    len = sizeof(T) * elements * extradim
    block = Block{T}(convert(Ptr{T}, Libc.malloc(len)), len, sizeof(T) * extradim)
    finalizer(block, x->Libc.free(x.ptr))
    return block
end

# used for getting messages back from ODBC driver manager; SQLDrivers, SQLError, etc.
Base.string(block::Block, len::Integer) = String(transcode(UInt8, unsafe_wrap(Array, block.ptr, len, false)))
#===================================================================================================
    </block>
===================================================================================================#

#===================================================================================================
    <ODBCError>
===================================================================================================#
type ODBCError <: Exception
    msg::String
end

const BUFLEN = 1024

function ODBCError(handle::Ptr{Void}, handletype::Int16)
    i = Int16(1)
    state = ODBC.Block(ODBC.API.SQLWCHAR, 6)
    native = Ref{ODBC.API.SQLINTEGER}()
    error_msg = ODBC.Block(ODBC.API.SQLWCHAR, BUFLEN)
    msg_length = Ref{ODBC.API.SQLSMALLINT}()
    while ODBC.API.SQLGetDiagRec(handletype, handle, i, state.ptr, native, error_msg.ptr, BUFLEN, msg_length) == ODBC.API.SQL_SUCCESS
        st  = string(state, 5)
        msg = string(error_msg, msg_length[])
        println("[ODBC] $st: $msg")
        i += 1
    end
    return true
end
ODBCError(hdl::SQLHandle) = ODBCError(hdl.ptr, handletype(hdl))
#===================================================================================================
    <ODBCError>
===================================================================================================#

"List ODBC drivers that have been installed and registered"
function listdrivers()
    descriptions = String[]
    attributes   = String[]
    driver_desc = Block(ODBC.API.SQLWCHAR, BUFLEN)
    desc_length = Ref{ODBC.API.SQLSMALLINT}()
    driver_attr = Block(ODBC.API.SQLWCHAR, BUFLEN)
    attr_length = Ref{ODBC.API.SQLSMALLINT}()
    dir = ODBC.API.SQL_FETCH_FIRST
    while ODBC.API.SQLDrivers(ENV, dir, driver_desc.ptr, BUFLEN, desc_length, driver_attr.ptr, BUFLEN, attr_length) == ODBC.API.SQL_SUCCESS
        push!(descriptions, string(driver_desc, desc_length[]))
        push!(attributes,   string(driver_attr, attr_length[]))
        dir = ODBC.API.SQL_FETCH_NEXT
    end
    return [descriptions attributes]
end

"List ODBC DSNs, both user and system, that have been previously defined"
function listdsns()
    descriptions = String[]
    attributes   = String[]
    dsn_desc    = Block(ODBC.API.SQLWCHAR, BUFLEN)
    desc_length = Ref{ODBC.API.SQLSMALLINT}()
    dsn_attr    = Block(ODBC.API.SQLWCHAR, BUFLEN)
    attr_length = Ref{ODBC.API.SQLSMALLINT}()
    dir = ODBC.API.SQL_FETCH_FIRST
    while ODBC.API.SQLDataSources(ENV, dir, dsn_desc.ptr, BUFLEN, desc_length, dsn_attr.ptr, BUFLEN, attr_length) == ODBC.API.SQL_SUCCESS
        push!(descriptions, string(dsn_desc, desc_length[]))
        push!(attributes,   string(dsn_attr, attr_length[]))
        dir = ODBC.API.SQL_FETCH_NEXT
    end
    return [descriptions attributes]
end

"""
A DSN represents an established ODBC connection.
It is passed to most other ODBC methods as a first argument
"""
type DSN
    dsn::String
    dbc_ptr::Ptr{Void}
    stmt_ptr::Ptr{Void}
    stmt_ptr2::Ptr{Void}
end

Base.show(io::IO,conn::DSN) = print(io, "ODBC.DSN($(conn.dsn))")

"""
Construct a `DSN` type by connecting to a valid ODBC DSN or by specifying a valid connection string.
Takes optional 2nd and 3rd arguments for `username` and `password`, respectively.
1st argument `dsn` can be either the name of a pre-defined ODBC DSN or a valid connection string.
A great resource for building valid connection strings is [http://www.connectionstrings.com/](http://www.connectionstrings.com/).
"""
function DSN(dsn::AbstractString, username::AbstractString=String(""), password::AbstractString=String(""); prompt::Bool=true)
    dbc = ODBC.ODBCAllocHandle(ODBC.API.SQL_HANDLE_DBC, ODBC.ENV)
    dsns = ODBC.listdsns()
    found = false
    for d in dsns[:,1]
        dsn == d && (found = true)
    end
    if found
        @CHECK dbc ODBC.API.SQL_HANDLE_DBC ODBC.API.SQLConnect(dbc, dsn, username, password)
    else
        dsn = ODBCDriverConnect!(dbc, dsn, prompt)
    end
    stmt = ODBCAllocHandle(ODBC.API.SQL_HANDLE_STMT, dbc)
    stmt2 = ODBCAllocHandle(ODBC.API.SQL_HANDLE_STMT, dbc)
    conn = DSN(dsn, dbc, stmt, stmt2)
    return conn
end

"disconnect a connected `DSN`"
function disconnect!(conn::DSN)
    ODBCFreeStmt!(conn.stmt_ptr)
    ODBCFreeStmt!(conn.stmt_ptr2)
    ODBC.API.SQLDisconnect(conn.dbc_ptr)
    return nothing
end

"An `ODBC.Source` type executes a `query` string upon construction and prepares data for streaming to an appropriate `Data.Sink`"
type Source <: Data.Source
    schema::Data.Schema
    dsn::DSN
    query::String
    columns::Vector{Any}
    status::Int
    rowsfetched::Ref{ODBC.API.SQLLEN}
    rowoffset::Int
    boundcols::Vector{Any}
    indcols::Vector{Vector{ODBC.API.SQLLEN}}
    sizes::Vector{ODBC.API.SQLULEN}
    ctypes::Vector{ODBC.API.SQLSMALLINT}
    jltypes::Vector{DataType}
end

Base.show(io::IO, source::Source) = print(io, "ODBC.Source:\n\tDSN: $(source.dsn)\n\tstatus: $(source.status)\n\tschema: $(source.schema)")

function __init__()
    global const SQLEnv = ODBCAllocHandle(SQL_HANDLE_ENV, SQL_NULL_HANDLE)
end

# used to 'clear' a statement of bound columns, resultsets,
# and other bound parameters in preparation for a subsequent query
function ODBCFreeStmt!(stmt)
    ODBC.API.SQLFreeStmt(stmt,ODBC.API.SQL_CLOSE)
    ODBC.API.SQLFreeStmt(stmt,ODBC.API.SQL_UNBIND)
    ODBC.API.SQLFreeStmt(stmt,ODBC.API.SQL_RESET_PARAMS)
end


text:
import Base: push!
using JLD2, FileIO

"""
these function save variables into the output folder.
"""
function save_variables(m::ContainerTypes, vars::Vector{String})
	for v in vars
		save_variables(m,v)
	end
end

"""
takes a variable, checks if we have a file for that variable name, and saves data in that file
"""
function save_variables(m::ContainerTypes, varname::String)
	newvar = get_variable(m,varname)
	fname = "C:\\Users\\James\\.julia\\dev\\MushroomBody\\src\\output\\Variables\\" * varname * ".jld2"
	if isfile(fname)
		oldvar = load(fname, varname)
		oldvar = !isa(oldvar,Array) ? [oldvar] : oldvar
		push!(oldvar, newvar)
	else
		oldvar = newvar
	end

	save(fname, varname, oldvar)

end

get_variable(m::ContainerTypes, var::String) = getproperty(m,var)

#=Base.push!(layers::Union{NeuronLayers,SynapseLayers},newvar) = layers=[layers]; push!([layers],newvar)
=#

@inline function return_variable(m::ContainerTypes, varname::String)
	copy(get_variable(m,varname))
end

@inline function initialise_return_variable(numsteps, m::ContainerTypes, varname::String)

	Array{typeof(get_variable(m,varname)),1}(undef,numsteps)

end


text:
function make_sky(light_pos::A, power::Integer, ambient::T, diffusePortion::T) where {T<:AbstractFloat, A<:AbstractArray{T,1}}
  skymap::AbstractArray{T,2} = zeros(3600, 900)
  azi::T = 0
  zen::T = 0
  for i=1:3600
      for j=1:900
          azi = float(i)/1800.0*pi
          zen = float(j)/1800.0*pi
          xi, yi, zi = spherical2carthesian(azi, zen)
          skymap[i,j] = background_radiance([xi, yi, zi], light_pos, power, ambient, diffusePortion)
      end
  end
  return skymap
end

function background_radiance(hemi_sample::Vector{T}, light_pos::Vector{T}, power::Int, ambient::T, diffusePortion) where {T<:AbstractFloat}
  # expand this function with sun position, atmospheric model, etc.
  direct::T = max(0.0, dot(light_pos, hemi_sample))^power
  diffuse::T = max(0.0, dot(light_pos, hemi_sample))^1
  directPortion::T = 1. - diffusePortion
  dirrad::T =  directPortion * direct
  diffrad::T = diffusePortion * diffuse
  radiance::T = (1.0-ambient) * (dirrad + diffrad) + ambient
  return radiance
end

function sample_skymap(skymap::AbstractArray{T, 2}, dx::T, dy::T, dz::T) where {T<:AbstractFloat}
    phi::T, theta::T = carthesian2spherical(dx, dy, dz)
    v::T = phi*1800.0/pi # convert to an index...
    i::Int = round(Int, v)
    w::T = theta*1800.0/pi
    j::Int = round(Int, w)
    if i === 0
        i = 1
    elseif i > 3600
        i = 3600
    end
    if j === 0
        j = 1
    elseif j > 900
        j = 900
    end
    radiance::T = skymap[i, j]
    return radiance
end

max_stars_repo_name: martinvanleeuwen/RenderJay

max_issues_repo_name: martinvanleeuwen/RenderJay

max_forks_repo_name: martinvanleeuwen/RenderJay


text:
module TestMockOptimizer

using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

# function test_optimizer_solve_no_result()
#     optimizer = MOIU.MockOptimizer(MOIU.Model{Float64}())

#     v1 = MOI.add_variable(optimizer)

#     # Load fake solution
#     MOI.set(optimizer, MOI.TerminationStatus(), MOI.INFEASIBLE)
#     MOI.optimize!(optimizer)
#     @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.INFEASIBLE

#     MOI.set(optimizer, MOI.ConflictStatus(), MOI.CONFLICT_FOUND)
#     MOI.compute_conflict!(optimizer)
#     @test MOI.get(optimizer, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
# end

function test_optimizer_solve_with_result()
    optimizer = MOIU.MockOptimizer(
        MOIU.Model{Float64}(),
        eval_objective_value = false,
        eval_variable_constraint_dual = false,
    )

    v = MOI.add_variables(optimizer, 2)
    c1 = MOI.add_constraint(optimizer, v[1], MOI.GreaterThan(1.0))
    soc = MOI.add_constraint(
        optimizer,
        MOI.VectorOfVariables(v),
        MOI.SecondOrderCone(2),
    )
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.VariableIndex}(), v[1])
    MOI.set(optimizer, MOI.ResultCount(), 1)
    @test_throws(
        ErrorException("No mock primal is set for variable `$(v[1])`."),
        MOI.get(optimizer, MOI.VariablePrimal(), v[1])
    )
    @test_throws(
        MOI.InvalidIndex(MOI.VariableIndex(-1)),
        MOI.get(optimizer, MOI.VariablePrimal(), MOI.VariableIndex(-1))
    )
    MOI.set(optimizer, MOI.ResultCount(), 0)
    err = MOI.ResultIndexBoundsError(MOI.VariablePrimal(1), 0)
    @test_throws err MOI.get(optimizer, MOI.VariablePrimal(), v[1])
    err = MOI.ResultIndexBoundsError(MOI.ConstraintDual(1), 0)
    @test_throws err MOI.get(optimizer, MOI.ConstraintDual(), c1)

    # Load fake solution
    # TODO: Provide a more compact API for this.
    MOI.set(optimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(optimizer, MOI.ObjectiveValue(), 1.0)
    MOI.set(optimizer, MOI.DualObjectiveValue(2), 2.0)
    MOI.set(optimizer, MOI.ResultCount(), 2)
    MOI.set(optimizer, MOI.PrimalStatus(), MOI.FEASIBLE_POINT)
    MOI.set(optimizer, MOI.DualStatus(2), MOI.FEASIBLE_POINT)
    MOI.set(optimizer, MOI.VariablePrimal(), v, [1.0, 2.0])
    MOI.set(optimizer, MOI.VariablePrimal(), v[1], 3.0)
    MOI.set(optimizer, MOI.ConstraintDual(2), c1, 5.9)
    MOI.set(optimizer, MOI.ConstraintDual(2), soc, [1.0, 2.0])

    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.ResultCount()) == 2
    @test MOI.get(optimizer, MOI.ObjectiveValue()) == 1.0
    @test isnan(MOI.get(optimizer, MOI.ObjectiveValue(2)))
    optimizer.eval_objective_value = true
    @test MOI.get(optimizer, MOI.ObjectiveValue()) == 3.0
    @test_throws(
        ErrorException(
            "No mock primal is set for variable `$(v[1])` at result index `2`.",
        ),
        MOI.get(optimizer, MOI.ObjectiveValue(2))
    )
    @test_throws(
        ErrorException(
            "No mock dual is set for constraint `$c1` at result index `1`.",
        ),
        MOI.get(optimizer, MOI.DualObjectiveValue())
    )
    @test MOI.get(optimizer, MOI.DualObjectiveValue(2)) == 5.9
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(optimizer, MOI.DualStatus(2)) == MOI.FEASIBLE_POINT
    @test MOI.get(optimizer, MOI.VariablePrimal(), v) == [3.0, 2.0]
    @test MOI.get(optimizer, MOI.VariablePrimal(), v[1]) == 3.0
    @test_throws(
        ErrorException(
            "No mock primal is set for variable `$(v[1])` at result index `2`.",
        ),
        MOI.get(optimizer, MOI.VariablePrimal(2), v[1])
    )
    @test MOI.get(optimizer, MOI.ConstraintPrimal(), c1) == 3.0
    @test MOI.get(optimizer, MOI.ConstraintPrimal(), soc) == [3.0, 2.0]
    @test_throws(
        ErrorException(
            "No mock primal is set for variable `$(v[1])` at result index `2`.",
        ),
        @show MOI.get(optimizer, MOI.ConstraintPrimal(2), c1)
    )
    @test_throws(
        ErrorException(
            "No mock primal is set for variable `$(v[1])` at result index `2`.",
        ),
        MOI.get(optimizer, MOI.ConstraintPrimal(2), soc)
    )
    @test MOI.get(optimizer, MOI.ConstraintDual(2), c1) == 5.9
    @test MOI.get(optimizer, MOI.ConstraintDual(2), soc) == [1.0, 2.0]
    @test_throws(
        ErrorException(
            "No mock dual is set for constraint `$c1` at result index `1`.",
        ),
        MOI.get(optimizer, MOI.ConstraintDual(1), c1)
    )
end

function test_CanonicalConstraintFunction()
    mock = MOIU.MockOptimizer(MOIU.Model{Int}())
    fx, fy = MOI.add_variables(mock, 2)
    cx = MOI.add_constraint(mock, fx, MOI.LessThan(0))
    c = MOI.add_constraint(mock, 1fx + fy, MOI.LessThan(1))
    @test MOIU.is_canonical(MOI.get(mock, MOI.ConstraintFunction(), cx))
    @test MOIU.is_canonical(
        MOI.get(mock, MOI.CanonicalConstraintFunction(), cx),
    )
    @test !MOIU.is_canonical(MOI.get(mock, MOI.ConstraintFunction(), c))
    func = MOI.get(mock, MOI.CanonicalConstraintFunction(), c)
    @test MOIU.is_canonical(func)
    # Check that indices have been xored
    for term in func.terms
        @test MOI.is_valid(mock, term.variable)
    end
end

function test_conflict_access()
    mock = MOIU.MockOptimizer(MOIU.Model{Int}())
    fx, fy = MOI.add_variables(mock, 2)
    cx = MOI.add_constraint(mock, fx, MOI.LessThan(0))
    c = MOI.add_constraint(mock, 1fx + fy, MOI.LessThan(1))
    MOI.set(mock, MOI.ConstraintConflictStatus(), cx, MOI.NOT_IN_CONFLICT)
    MOI.set(mock, MOI.ConstraintConflictStatus(), c, MOI.IN_CONFLICT)
    MOI.compute_conflict!(mock)

    @test MOI.get(mock, MOI.ConstraintConflictStatus(), cx) ==
          MOI.NOT_IN_CONFLICT
    @test MOI.get(mock, MOI.ConstraintConflictStatus(), c) == MOI.IN_CONFLICT
end

function test_MockVariableAttribute()
    mock = MOIU.MockOptimizer(MOIU.Model{Int}())
    x = MOI.add_variable(mock)
    MOI.set(mock, MOI.Utilities.MockVariableAttribute(), x, 1)
    @test MOI.get(mock, MOI.Utilities.MockVariableAttribute(), x) == 1
    return
end

function test_MockConstraintAttribute()
    mock = MOIU.MockOptimizer(MOIU.Model{Int}())
    x = MOI.add_variable(mock)
    c = MOI.add_constraint(mock, x, MOI.LessThan(0))
    MOI.set(mock, MOI.Utilities.MockConstraintAttribute(), c, 1)
    @test MOI.get(mock, MOI.Utilities.MockConstraintAttribute(), c) == 1
    return
end

function test_DualObjectiveValue()
    mock =
        MOIU.MockOptimizer(MOIU.Model{Int}(); eval_dual_objective_value = false)
    @test isnan(MOI.get(mock, MOI.DualObjectiveValue()))
    return
end

end  # module

TestMockOptimizer.runtests()


text:
# Standard plugin template

function pluginCallback(vis::MeshCat.Visualizer,
                        params::Dict{Symbol, Any},
                        rose_fgl::Union{FactorGraph, Tuple{<:AbstractString, <:AbstractString}} )
    #
    @show length(params[:cachevars])

    @error "Implement the plugin functionality here"

    nothing
end


text:
module IntervalLinearAlgebra

using StaticArrays, Requires, Reexport
using LinearAlgebra: checksquare

import Base: +, -, *, /, \, ==,
            show, convert, promote_rule, zero, one,
            getindex, IndexStyle, setindex!, size
import CommonSolve: solve
import IntervalArithmetic: mid

@reexport using LinearAlgebra, IntervalArithmetic

const IA = IntervalArithmetic

export
    set_multiplication_mode, get_multiplication_mode,
    LinearKrawczyk, Jacobi, GaussSeidel, GaussianElimination, HansenBliekRohn, NonLinearOettliPrager, LinearOettliPrager,
    NoPrecondition, InverseMidpoint, InverseDiagonalMidpoint,
    solve, enclose, epsilon_inflation,
    comparison_matrix, interval_norm, interval_isapprox, Orthants,
    is_H_matrix, is_strongly_regular, is_strictly_diagonally_dominant, is_Z_matrix, is_M_matrix,
    rref,
    eigenbox, Rohn, Hertz, verify_eigen, bound_perron_frobenius_eigenvalue,
    AffineExpression, @linvars,
    AffineParametricArray, AffineParametricMatrix, AffineParametricVector

include("linear_systems/enclosures.jl")
include("linear_systems/precondition.jl")
include("linear_systems/solve.jl")
include("linear_systems/verify.jl")
include("linear_systems/oettli.jl")
include("multiplication.jl")
include("utils.jl")
include("classify.jl")
include("rref.jl")
include("pils/linexpr.jl")
include("pils/affine_parametric_array.jl")

include("eigenvalues/interval_eigenvalues.jl")
include("eigenvalues/verify_eigs.jl")

function  __init__()
    @require IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153" include("linear_systems/oettli_nonlinear.jl")
    @require LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" include("linear_systems/oettli_linear.jl")
end

set_multiplication_mode(config[:multiplication])

end

max_stars_repo_name: lucaferranti/IntervalLinearAlgebra

max_issues_repo_name: lucaferranti/IntervalLinearAlgebra

max_forks_repo_name: lucaferranti/IntervalLinearAlgebra


text:
module Wpm

import GreyDecision.GreyNumbers: GreyNumber
import GreyDecision.Utility: unitize

function wpm(decisionMat::Array{GreyNumber{T}, 2}, weights::Array{Float64,1},
    fns::Array{Function,1}) where {T <: Real}

    row, col = size(decisionMat)
    normalizedDecisionMat = similar(decisionMat)
    w = unitize(weights)
    colminmax = zeros(GreyNumber{T}, col)

    @inbounds for i in 1:col
        colminmax[i] = decisionMat[:, i] |> fns[i]
        if fns[i] == maximum
            normalizedDecisionMat[:, i] = map(x -> x / colminmax[i], decisionMat[:, i])
        elseif fns[i] == minimum
            normalizedDecisionMat[:, i] = map(x -> colminmax[i] / x, decisionMat[:, i])
        end
    end

    scoreMat = similar(normalizedDecisionMat)
    for i in 1:col
        scoreMat[:, i] = normalizedDecisionMat[:, i].^w[i]
    end

     scores = zeros(GreyNumber{Float64}, row)
    for i in 1:row
        scores[i] = prod(scoreMat[i, :])
    end
    rankings = sortperm(scores)

    bestIndex = rankings |> last

    result = Dict(
        :decisionMat => decisionMat,
        :normalizedDecisionMat => normalizedDecisionMat,
        :w => w,
        :scores => scores,
        :ranking => rankings,
        :bestIndex => bestIndex
    )

    return result
end

end # End of module


text:
include("PerFlagFuncs.jl")
include("FeaturedFuncs.jl")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Qt6Svg"
version = v"6.3.0"

# Collection of sources required to build qt6
sources = [
    ArchiveSource("https://download.qt.io/official_releases/qt/$(version.major).$(version.minor)/$version/submodules/qtsvg-everywhere-src-$version.tar.xz",
                  "3164504d7e3f640439308235739b112605ab5fc9cc517ca0b28f9fb93a8db0e3"),
    ArchiveSource("https://github.com/phracker/MacOSX-SDKs/releases/download/11.0-11.1/MacOSX11.1.sdk.tar.xz",
                  "9b86eab03176c56bb526de30daa50fa819937c54b280364784ce431885341bf6"),
]

script = raw"""
cd $WORKSPACE/srcdir

mkdir build
cd build/
qtsrcdir=`ls -d ../qtsvg-*`

case "$target" in

    x86_64-linux-musl*)
        cmake -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_FIND_ROOT_PATH=$prefix -DCMAKE_BUILD_TYPE=Release $qtsrcdir
    ;;

    *apple-darwin*)
        apple_sdk_root=$WORKSPACE/srcdir/MacOSX11.1.sdk
        sed -i "s!/opt/x86_64-apple-darwin14/x86_64-apple-darwin14/sys-root!$apple_sdk_root!" $CMAKE_TARGET_TOOLCHAIN
        cmake -DQT_HOST_PATH=$host_prefix \
            -DPython_ROOT_DIR=/usr \
            -DCMAKE_INSTALL_PREFIX=${prefix} \
            -DCMAKE_PREFIX_PATH=$host_prefix \
            -DCMAKE_FIND_ROOT_PATH=$prefix \
            -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
            -DCMAKE_SYSROOT=$apple_sdk_root -DCMAKE_FRAMEWORK_PATH=$apple_sdk_root/System/Library/Frameworks -DCMAKE_OSX_DEPLOYMENT_TARGET=10.15 \
            -DCMAKE_BUILD_TYPE=Release \
        $qtsrcdir
    ;;

    *)
        cmake -DQT_HOST_PATH=$host_prefix -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_FIND_ROOT_PATH=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release $qtsrcdir
    ;;

esac

cmake --build . --parallel ${nproc}
cmake --install .
install_license $WORKSPACE/srcdir/qt*-src-*/LICENSE.LGPL3
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(filter(!Sys.isapple, supported_platforms()))
filter!(p -> arch(p) != "armv6l", platforms) # No OpenGL on armv6
platforms_macos = [ Platform("x86_64", "macos"), Platform("aarch64", "macos") ]

# The products that we will ensure are always built
products = [
    LibraryProduct(["Qt6Svg", "libQt6Svg", "QtSvg"], :libqt6svg),
    LibraryProduct(["Qt6SvgWidgets", "libQt6SvgWidgets", "QtSvgWidgets"], :libqt6svgwidgets),
]

products_macos = [
    FrameworkProduct("QtSvg", :libqt6svg),
    FrameworkProduct("QtSvgWidgets", :libqt6svgwidgets),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    HostBuildDependency("Qt6Base_jll"),
    Dependency("Qt6Base_jll"),
]

include("../../fancy_toys.jl")

if any(should_build_platform.(triplet.(platforms_macos)))
    build_tarballs(ARGS, name, version, sources, script, platforms_macos, products_macos, dependencies; preferred_gcc_version = v"9", julia_compat="1.6")
end
if any(should_build_platform.(triplet.(platforms)))
    build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"9", julia_compat="1.6")
end


text:
if VERSION >= v"0.5-"
    using BaseTestDeprecated
    const Test = BaseTestDeprecated
else
    using Base.Test
end
include(joinpath(dirname(@__FILE__), "kabsch.jl"))
include(joinpath(dirname(@__FILE__), "plot.jl"))
include(joinpath(dirname(@__FILE__), "pdb.jl"))
include(joinpath(dirname(@__FILE__), "hit.jl"))
# mafft tests do not run reliably (especially with julia v0.5)
if VERSION < v"0.5-"
	#include(joinpath(dirname(@__FILE__), "mafft.jl"))
end
include(joinpath(dirname(@__FILE__), "webblast.jl"))


text:
# RwaJur3

export VERTICES_RWAJUR3,
		RwaJur3Exposure,
		RwaJur3Parameters,
		getRwaJur3Parameters,
		getLongExposureArray,
		getShortExposureArray,
		getRwaJur3NetExposure,
		getRwaJur3VerticalGap,
		getRwaJur3Zones,
		getRwaJur3HorizontalGapWithinZones,
		getRwaJur3HorizontalGapBetweenZones,
		getPJur3,
		getRwaJur3

"""
Standardized vertices for mapping price index coupon rate exposures (RWAJUR3)
"""
const VERTICES_RWAJUR3 = VERTICES_RWAJUR2

"""
    RwaJur3Exposure

Holds exposures in price index coupon rates.

# Fields
* `date::Date`: date of evaluation
* `VXX::Float64`: exposure to vertex XX,
given by the sum of present value in BRL
of all cashflows that matures in XX business days (BD).

# Constructors
    RwaJur3Exposure(date::Date)

Create an instance with all exposures equal to zero.

	RwaJur2Exposure(date::Date, currency::AbstractString, longExposures::Array{Float64, 1}, shortExposures::Array{Float64, 1})

Create an instance with exposures given by two vectors of length 10 containing long and short exposures.
Exposures are assigned in sequence from vertex 1 BD to vertex 2520 BD.

"""
type RwaJur3Exposure
	# Date of evaluation
	date::Date
	# Currency
	currency::AbstractString
	# Exposures in standard vertices: VXX, where XX is a maturity calculated in business days
	# LONG EXPOSURES
	V1Long::Float64
	V21Long::Float64
	V42Long::Float64
	V63Long::Float64
	V126Long::Float64
	V252Long::Float64
	V504Long::Float64
	V756Long::Float64
	V1008Long::Float64
	V1260Long::Float64
	V2520Long::Float64
	# SHORT EXPOSURES
	V1Short::Float64
	V21Short::Float64
	V42Short::Float64
	V63Short::Float64
	V126Short::Float64
	V252Short::Float64
	V504Short::Float64
	V756Short::Float64
	V1008Short::Float64
	V1260Short::Float64
	V2520Short::Float64

	# Constructors
	RwaJur3Exposure(date::Date, currency::AbstractString) = new(date, currency, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

	function RwaJur3Exposure(date::Date, currency::AbstractString, longExposures::Array{Float64, 1}, shortExposures::Array{Float64, 1})
		if size(longExposures, 1) != 11
			error("long exposures must have 11 values")
		end #if
		if any(longExposures .< 0)
			error("all long exposures must be positive")
		end #if

		if size(shortExposures, 1) != 11
			error("short exposures must have 11 values")
		end #if
		if any(shortExposures .< 0)
			error("all short exposures must be positive")
		end #if

		new(date,
			currency,
		    longExposures[1], longExposures[2], longExposures[3], longExposures[4],
			longExposures[5], longExposures[6], longExposures[7], longExposures[8],
			longExposures[9], longExposures[10], longExposures[11],
			shortExposures[1], shortExposures[2], shortExposures[3], shortExposures[4],
			shortExposures[5], shortExposures[6], shortExposures[7], shortExposures[8],
			shortExposures[9], shortExposures[10], shortExposures[11]
			)
	end #function

end #type

"""
	RWAJUR3Parameters

Holds BCB parameters to calculate price index coupon stantardized capital (RWA JUR3)

"""
type RwaJur3Parameters
	date::Date
	M3::Float64
end #type

"""
    getRwaJur3Parameters(date::Date)

Return RWA JUR3 parameters supplied by BCB for `date::Date`.
"""
function getRwaJur3Parameters(date::Date)
	if date <= Date(2012, 4, 29)
		M3 = 1.93
	elseif date <= Date(2012, 8, 30)
		M3 = 2.19
	elseif date <= Date(2012, 12, 30)
		M3 = 2.45
	else
		M3 = 2.7
	end #if
	return(RwaJur3Parameters(date, M3))
end #function

"""
    getLongExposureArray(rExp::RwaJur3Exposure)

Return a vector with all long exposures in sequence,
from vertex 1 BD to vertex 2520 BD.

"""
getLongExposureArray(rExp::RwaJur3Exposure) = [rExp.V1Long, rExp.V21Long, rExp.V42Long, rExp.V63Long, rExp.V126Long, rExp.V252Long, rExp.V504Long, rExp.V756Long, rExp.V1008Long, rExp.V1260Long, rExp.V2520Long]

"""
    getShortExposureArray(rExp::RwaJur3Exposure)

Return a vector with all short exposures in sequence,
from vertex 1 BD to vertex 2520 BD.

"""
getShortExposureArray(rExp::RwaJur3Exposure) = [rExp.V1Short, rExp.V21Short, rExp.V42Short, rExp.V63Short, rExp.V126Short, rExp.V252Short, rExp.V504Short, rExp.V756Short, rExp.V1008Short, rExp.V1260Short, rExp.V2520Short]

# Casts
"""
    convert(::Type{RwaJur2Exposure}, x::RwaJur3Exposure)

Casts an RwaJur3Exposure into an RwaJur2Exposure.
"""
convert(::Type{RwaJur2Exposure}, x::RwaJur3Exposure) = RwaJur2Exposure(x.date, x.currency, getLongExposureArray(x), getShortExposureArray(x))

"""
    convert(::Type{RwaJur2Parameters}, x::RwaJur3Parameters)

Casts an RwaJur3Parameters into an RwaJur2Parameters.
"""
convert(::Type{RwaJur2Parameters}, x::RwaJur3Parameters) = RwaJur2Parameters(x.date, x.M3)

"""
    getRwaJur3NetExposure(rExp::RwaJur3Exposure)

Return a vector containing RWA JUR3 weighted net exposures.

"""
getRwaJur3NetExposure(rExp::RwaJur3Exposure) = getRwaJur2NetExposure(convert(RwaJur2Exposure, rExp))

"""
    getRwaJur3VerticalGap(rExp::RwaJur3Exposure)

Return a vector containing RWA JUR3 vertical gaps.

"""
getRwaJur3VerticalGap(rExp::RwaJur3Exposure) = getRwaJur2VerticalGap(convert(RwaJur2Exposure, rExp))

"""
    getRwaJur3Zones(rExp::RwaJur3Exposure)

Return a matrix containing RWA JUR3 zone exposures.
Column 1 contains long exposures, column2 contains short exposures.
Zones are represented in each of the 3 rows.

"""
getRwaJur3Zones(rExp::RwaJur3Exposure) = getRwaJur2Zones(convert(RwaJur2Exposure, rExp))

"""
    getRwaJur3HorizontalGapWithinZones(rExp::RwaJur3Exposure)

Return a vector containing RWA JUR3 horizontal gap within zones.

"""
getRwaJur3HorizontalGapWithinZones(rExp::RwaJur3Exposure) = getRwaJur2HorizontalGapWithinZones(convert(RwaJur2Exposure, rExp))

"""
    getRwaJur3HorizontalGapBetweenZones(rExp::RwaJur3Exposure)

Return a vector containing RWA JUR3 horizontal gap between zones.

"""
getRwaJur3HorizontalGapBetweenZones(rExp::RwaJur3Exposure) = getRwaJur2HorizontalGapBetweenZones(convert(RwaJur2Exposure, rExp))

"""
	getPJur3(rExp::RwaJur3Exposure, par::RwaJur3Parameters=getRwaJur3Parameters(rExp.date))

Calculate PJUR3 using maturity ladder methodology.

"""
getPJur3(rExp::RwaJur3Exposure, par=getRwaJur3Parameters(rExp.date)) = getPJur2(convert(RwaJur2Exposure, rExp), convert(RwaJur2Parameters, par))

"""
	getRwaJur3(rExp::RwaJur3Exposure, par::RwaJur3Parameters=getRwaJur3Parameters(rExp.date))

Calculate risk-weighted assets for price index coupon exposures (RWA JUR3),
using maturity ladder methodology.

"""
getRwaJur3(rExp::RwaJur3Exposure, par=getRwaJur3Parameters(rExp.date)) = getRwaJur2(convert(RwaJur2Exposure, rExp), convert(RwaJur2Parameters, par))


text:
# LSP indexing solution
# https://github.com/julia-vscode/julia-vscode/issues/800#issuecomment-650085983
if isdefined(@__MODULE__, :LanguageServer)
    include("../src/SpinDoctor.jl")
    using .SpinDoctor
end

using SpinDoctor
using LinearAlgebra
using GLMakie
using OrdinaryDiffEq: Rodas4

## Chose a plotting theme
set_theme!(theme_light())
set_theme!(theme_dark())
set_theme!(theme_black())

## Create model from setup recipe
# include("setups/axon.jl")
# include("setups/sphere.jl")
# include("setups/plates.jl")
# include("setups/cylinders.jl")
# include("setups/spheres.jl")
include("setups/neuron.jl")

mesh, = @time create_geometry(setup; recreate = true);
model = Model(; mesh, coeffs...);
volumes = get_cmpt_volumes(model.mesh)
D_avg = 1 / 3 * tr.(model.D)' * volumes / sum(volumes)
@info "Number of nodes per compartment:" length.(model.mesh.points)

## Plot mesh
plot_mesh(model.mesh)

## Assemble finite element matrices
matrices = @time assemble_matrices(model);

## Magnetic field gradient
dir = [1.0, 0.0, 0.0]
profile = PGSE(2000.0, 6000.0)
# profile = CosOGSE(5000.0, 5000.0, 2)
b = 1000
g = √(b / int_F²(profile)) / coeffs.γ
gradient = ScalarGradient(dir, profile, g)

## Solve BTPDE

# Callbacks for time stepping (plot solution, save time series)
printer = Printer(; nupdate = 1, verbosity = 2)
writer = VTKWriter(; nupdate = 5)
plotter = Plotter{T}(; nupdate = 5)
# callbacks = [printer, plotter]
callbacks = [printer, plotter, writer]

# General BTPDE for all gradients
btpde = GeneralBTPDE(;
    model,
    matrices,
    reltol = 1e-4,
    abstol = 1e-6,
    odesolver = Rodas4(autodiff = false),
)

# BTPDE specialized for `ScalarGradient`s with constant interval profiles (PGSE, DoublePGSE)
btpde = IntervalConstantBTPDE{T}(; model, matrices, θ = 0.5, timestep = 5)

# Solve BTPDE
ξ = @time solve(btpde, gradient; callbacks)

## Plot magnetization
plot_field(model.mesh, ξ)

## Compute signal
compute_signal(matrices.M, ξ)
compute_signal.(matrices.M_cmpts, split_field(model.mesh, ξ))

## Save magnetization
savefield(model.mesh, ξ, "output/magnetization")

## Matrix Formalism

# Perform Laplace eigendecomposition
laplace = Laplace{T}(; model, matrices, neig_max = 400)
lap_eig = @time solve(laplace)
length_scales = eig2length.(lap_eig.values, D_avg)

# Truncate basis at minimum length scale
length_scale = 3
λ_max = length2eig(length_scale, D_avg)
lap_eig = limit_lengthscale(lap_eig, λ_max)

# Compute magnetization using the matrix formalism reduced order model
mf = MatrixFormalism(; model, matrices, lap_eig, ninterval = 500)
ξ = @time solve(mf, gradient)

## Solve HADC
hadc = HADC(; model, matrices, reltol = 1e-4, abstol = 1e-6)
adc_cmpts = @time solve(hadc, gradient)

## Solve Karger model

# Compute HADC and fit difftensors
directions = unitsphere(50)
gradients = [
    ScalarGradient(collect(d), gradient.profile, gradient.amplitude) for
    d ∈ eachcol(directions)
]
hadc = HADC(; model, matrices, reltol = 1e-4, abstol = 1e-6)
adcs, = @time solve_multigrad(hadc, gradients)
difftensors = fit_tensors(directions, adcs)

# Solve Karger
karger = Karger(; model, difftensors, timestep = 5.0)
signal = @time solve(karger, gradient)

## Solve analytical model
# Compute analytical Laplace eigenfunctions
length_scale = 0.3
eigstep = 1e-8
eiglim = length2eig(length_scale, D_avg)
analytical_coeffs = analytical_coefficients(setup, coeffs)
analytical_laplace = AnalyticalLaplace(; analytical_coeffs..., eiglim, eigstep)
lap_mat = @time solve(analytical_laplace)

# Compute analytical matrix formalism signal truncation
analytical_mf = AnalyticalMatrixFormalism(; analytical_laplace, lap_mat, volumes)
signal = solve(analytical_mf, gradient)


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

mutable struct VirtualMachineScaleSetStorageProfile <: SwaggerModel
    imageReference::Any # spec type: Union{ Nothing, ImageReference } # spec name: imageReference
    osDisk::Any # spec type: Union{ Nothing, VirtualMachineScaleSetOSDisk } # spec name: osDisk
    dataDisks::Any # spec type: Union{ Nothing, Vector{VirtualMachineScaleSetDataDisk} } # spec name: dataDisks

    function VirtualMachineScaleSetStorageProfile(;imageReference=nothing, osDisk=nothing, dataDisks=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetStorageProfile, Symbol("imageReference"), imageReference)
        setfield!(o, Symbol("imageReference"), imageReference)
        validate_property(VirtualMachineScaleSetStorageProfile, Symbol("osDisk"), osDisk)
        setfield!(o, Symbol("osDisk"), osDisk)
        validate_property(VirtualMachineScaleSetStorageProfile, Symbol("dataDisks"), dataDisks)
        setfield!(o, Symbol("dataDisks"), dataDisks)
        o
    end
end # type VirtualMachineScaleSetStorageProfile

const _property_map_VirtualMachineScaleSetStorageProfile = Dict{Symbol,Symbol}(Symbol("imageReference")=>Symbol("imageReference"), Symbol("osDisk")=>Symbol("osDisk"), Symbol("dataDisks")=>Symbol("dataDisks"))
const _property_types_VirtualMachineScaleSetStorageProfile = Dict{Symbol,String}(Symbol("imageReference")=>"ImageReference", Symbol("osDisk")=>"VirtualMachineScaleSetOSDisk", Symbol("dataDisks")=>"Vector{VirtualMachineScaleSetDataDisk}")
Base.propertynames(::Type{ VirtualMachineScaleSetStorageProfile }) = collect(keys(_property_map_VirtualMachineScaleSetStorageProfile))
Swagger.property_type(::Type{ VirtualMachineScaleSetStorageProfile }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineScaleSetStorageProfile[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetStorageProfile }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetStorageProfile[property_name]

function check_required(o::VirtualMachineScaleSetStorageProfile)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetStorageProfile }, name::Symbol, val)
end


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

# Factorials

const _fact_table64 = Vector{Int64}(undef, 20)
_fact_table64[1] = 1
for n in 2:20
    _fact_table64[n] = _fact_table64[n-1] * n
end

const _fact_table128 = Vector{UInt128}(undef, 34)
_fact_table128[1] = 1
for n in 2:34
    _fact_table128[n] = _fact_table128[n-1] * n
end

function factorial_lookup(n::Integer, table, lim)
    n < 0 && throw(DomainError(n, "`n` must not be negative."))
    n > lim && throw(OverflowError(string(n, " is too large to look up in the table; consider using `factorial(big(", n, "))` instead")))
    n == 0 && return one(n)
    @inbounds f = table[n]
    return oftype(n, f)
end

factorial(n::Int128) = factorial_lookup(n, _fact_table128, 33)
factorial(n::UInt128) = factorial_lookup(n, _fact_table128, 34)
factorial(n::Union{Int64,UInt64}) = factorial_lookup(n, _fact_table64, 20)

if Int === Int32
    factorial(n::Union{Int8,UInt8,Int16,UInt16}) = factorial(Int32(n))
    factorial(n::Union{Int32,UInt32}) = factorial_lookup(n, _fact_table64, 12)
else
    factorial(n::Union{Int8,UInt8,Int16,UInt16,Int32,UInt32}) = factorial(Int64(n))
end

# Basic functions for working with permutations

@inline function _foldoneto(op, acc, ::Val{N}) where N
    @assert N::Integer > 0
    if @generated
        quote
            acc_0 = acc
            Base.Cartesian.@nexprs $N i -> acc_{i} = op(acc_{i-1}, i)
            return $(Symbol(:acc_, N))
        end
    else
        for i in 1:N
            acc = op(acc, i)
        end
        return acc
    end
end

"""
    isperm(v) -> Bool

Return `true` if `v` is a valid permutation.

# Examples
```jldoctest
julia> isperm([1; 2])
true

julia> isperm([1; 3])
false
```
"""
isperm(A) = _isperm(A)

function _isperm(A)
    n = length(A)
    used = falses(n)
    for a in A
        (0 < a <= n) && (used[a] ⊻= true) || return false
    end
    true
end

isperm(p::Tuple{}) = true
isperm(p::Tuple{Int}) = p[1] == 1
isperm(p::Tuple{Int,Int}) = ((p[1] == 1) & (p[2] == 2)) | ((p[1] == 2) & (p[2] == 1))

function isperm(P::Tuple)
    valn = Val(length(P))
    _foldoneto(true, valn) do b,i
        s = _foldoneto(false, valn) do s, j
            s || P[j]==i
        end
        b&s
    end
end

isperm(P::Any16) = _isperm(P)

# swap columns i and j of a, in-place
function swapcols!(a::AbstractMatrix, i, j)
    i == j && return
    cols = axes(a,2)
    @boundscheck i in cols || throw(BoundsError(a, (:,i)))
    @boundscheck j in cols || throw(BoundsError(a, (:,j)))
    for k in axes(a,1)
        @inbounds a[k,i],a[k,j] = a[k,j],a[k,i]
    end
end
# like permute!! applied to each row of a, in-place in a (overwriting p).
function permutecols!!(a::AbstractMatrix, p::AbstractVector{<:Integer})
    require_one_based_indexing(a, p)
    count = 0
    start = 0
    while count < length(p)
        ptr = start = findnext(!iszero, p, start+1)::Int
        next = p[start]
        count += 1
        while next != start
            swapcols!(a, ptr, next)
            p[ptr] = 0
            ptr = next
            next = p[next]
            count += 1
        end
        p[ptr] = 0
    end
    a
end

function permute!!(a, p::AbstractVector{<:Integer})
    require_one_based_indexing(a, p)
    count = 0
    start = 0
    while count < length(a)
        ptr = start = findnext(!iszero, p, start+1)::Int
        temp = a[start]
        next = p[start]
        count += 1
        while next != start
            a[ptr] = a[next]
            p[ptr] = 0
            ptr = next
            next = p[next]
            count += 1
        end
        a[ptr] = temp
        p[ptr] = 0
    end
    a
end

"""
    permute!(v, p)

Permute vector `v` in-place, according to permutation `p`. No checking is done
to verify that `p` is a permutation.

To return a new permutation, use `v[p]`. Note that this is generally faster than
`permute!(v,p)` for large vectors.

See also [`invpermute!`](@ref).

# Examples
```jldoctest
julia> A = [1, 1, 3, 4];

julia> perm = [2, 4, 3, 1];

julia> permute!(A, perm);

julia> A
4-element Array{Int64,1}:
 1
 4
 3
 1
```
"""
permute!(a, p::AbstractVector) = permute!!(a, copymutable(p))

function invpermute!!(a, p::AbstractVector{<:Integer})
    require_one_based_indexing(a, p)
    count = 0
    start = 0
    while count < length(a)
        start = findnext(!iszero, p, start+1)::Int
        temp = a[start]
        next = p[start]
        count += 1
        while next != start
            temp_next = a[next]
            a[next] = temp
            temp = temp_next
            ptr = p[next]
            p[next] = 0
            next = ptr
            count += 1
        end
        a[next] = temp
        p[next] = 0
    end
    a
end

"""
    invpermute!(v, p)

Like [`permute!`](@ref), but the inverse of the given permutation is applied.

# Examples
```jldoctest
julia> A = [1, 1, 3, 4];

julia> perm = [2, 4, 3, 1];

julia> invpermute!(A, perm);

julia> A
4-element Array{Int64,1}:
 4
 1
 3
 1
```
"""
invpermute!(a, p::AbstractVector) = invpermute!!(a, copymutable(p))

"""
    invperm(v)

Return the inverse permutation of `v`.
If `B = A[v]`, then `A == B[invperm(v)]`.

# Examples
```jldoctest
julia> v = [2; 4; 3; 1];

julia> invperm(v)
4-element Array{Int64,1}:
 4
 1
 3
 2

julia> A = ['a','b','c','d'];

julia> B = A[v]
4-element Array{Char,1}:
 'b': ASCII/Unicode U+0062 (category Ll: Letter, lowercase)
 'd': ASCII/Unicode U+0064 (category Ll: Letter, lowercase)
 'c': ASCII/Unicode U+0063 (category Ll: Letter, lowercase)
 'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)

julia> B[invperm(v)]
4-element Array{Char,1}:
 'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
 'b': ASCII/Unicode U+0062 (category Ll: Letter, lowercase)
 'c': ASCII/Unicode U+0063 (category Ll: Letter, lowercase)
 'd': ASCII/Unicode U+0064 (category Ll: Letter, lowercase)
```
"""
function invperm(a::AbstractVector)
    require_one_based_indexing(a)
    b = zero(a) # similar vector of zeros
    n = length(a)
    @inbounds for (i, j) in enumerate(a)
        ((1 <= j <= n) && b[j] == 0) ||
            throw(ArgumentError("argument is not a permutation"))
        b[j] = i
    end
    b
end

function invperm(p::Union{Tuple{},Tuple{Int},Tuple{Int,Int}})
    isperm(p) || throw(ArgumentError("argument is not a permutation"))
    p  # in dimensions 0-2, every permutation is its own inverse
end

function invperm(P::Tuple)
    valn = Val(length(P))
    ntuple(valn) do i
        s = _foldoneto(nothing, valn) do s, j
            s !== nothing && return s
            P[j]==i && return j
            nothing
        end
        s === nothing && throw(ArgumentError("argument is not a permutation"))
        s
    end
end

invperm(P::Any16) = Tuple(invperm(collect(P)))

#XXX This function should be moved to Combinatorics.jl but is currently used by Base.DSP.
"""
    nextprod(factors::Union{Tuple,AbstractVector}, n)

Next integer greater than or equal to `n` that can be written as ``\\prod k_i^{p_i}`` for integers
``p_1``, ``p_2``, etcetera, for factors ``k_i`` in `factors`.

# Examples
```jldoctest
julia> nextprod((2, 3), 105)
108

julia> 2^2 * 3^3
108
```

!!! compat "Julia 1.6"
    The method that accepts a tuple requires Julia 1.6 or later.
"""
function nextprod(a::Union{Tuple{Vararg{<:Integer}},AbstractVector{<:Integer}}, x::Real)
    if x > typemax(Int)
        throw(ArgumentError("unsafe for x > typemax(Int), got $x"))
    end
    k = length(a)
    v = fill(1, k)                    # current value of each counter
    mx = map(a -> nextpow(a,x), a)   # maximum value of each counter
    v[1] = mx[1]                      # start at first case that is >= x
    p::widen(Int) = mx[1]             # initial value of product in this case
    best = p
    icarry = 1

    while v[end] < mx[end]
        if p >= x
            best = p < best ? p : best  # keep the best found yet
            carrytest = true
            while carrytest
                p = div(p, v[icarry])
                v[icarry] = 1
                icarry += 1
                p *= a[icarry]
                v[icarry] *= a[icarry]
                carrytest = v[icarry] > mx[icarry] && icarry < k
            end
            if p < x
                icarry = 1
            end
        else
            while p < x
                p *= a[1]
                v[1] *= a[1]
            end
        end
    end
    # might overflow, but want predictable return type
    return mx[end] < best ? Int(mx[end]) : Int(best)
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
module FileRules

using FilePathsBase: Path, hasparent, iswritable

export Touch, Write

struct Touch{T}
    path::T
    filemode::Integer
    dirmode::Unsigned
end
Touch(path; filemode = 0o777, dirmode = 0o777) = Touch(path, filemode, dirmode)
function (t::Touch)()
    p = expanduser(t.path) |> Path
    if hasparent(p)
        mkpath(joinpath(splitpath(p)[1:end-1]...); mode = t.dirmode)
    end
    return chmod(touch(p), t.filemode)
end

struct Write{T}
    path::T
    mode::AbstractString
    function Write{T}(path, mode) where {T}
        @assert string(lowercase(mode)) ∈ ("r", "w", "a", "r+", "w+", "a+")
        return new(path, lowercase(mode))
    end
end
Write(path::T, mode) where {T} = Write{T}(path, mode)
function (w::Write)(x)
    p = expanduser(w.path) |> Path
    if !isfile(p)
        Touch(p)()
    end
    if !iswritable(p)
        error("file $p is not writable!")
    end
    if w.mode ∈ ("w", "w+")
        @info "try to overwrite file $p, a backup file will be created!"
        Touch(p * ".bkp")()
    end
    open(p, w.mode) do io
        write(io, x)
    end
end

end


text:
#! /home/kunal/julia/julia

using RobotOS
@rosimport geometry_msgs.msg: Point, Pose, PoseStamped
rostypegen()
using .geometry_msgs.msg

function poseCB(msg::PoseStamped, pose)
    # extract the pose
    pose = msg.pose
end

function loop(quadPoses)
    loop_rate = Rate(5.0)
    while ! is_shutdown()
        for pose in quadPoses
            print("x:", pose.position.x, " y:", pose.position.y, " z:", pose.position.z, "\n")
        end
        rossleep(loop_rate)
    end
end

function main()
    init_node("getPose0_jl")
    quadNames = ["rexquad0"]
    quadPoses = [Pose()]
    poseTopic = "/mavros/local_position/pose"
    quadSubs = [Subscriber{PoseStamped}(name*poseTopic, poseCB, (pose,), queue_size=10) for (name, pose) in zip(quadNames, quadPoses)]
    loop(quadPoses)
end

main()

max_stars_repo_name: StanfordMSL/mslquad

max_forks_repo_name: StanfordMSL/mslquad