text:
module Test_ode_bouncingBall

#=
Try EventHandler of ModiaMath for DifferentialEquations.jl
to figure out what needs to be adapted.

Evaluate Measurements.jl for events (does not work yet in a good way; result is useless)

=#

using DifferentialEquations, Measurements, StaticArrays, PyPlot

MeanValue(v::Number) = v
MeanValue(v::Measurement{Float64}) = Measurements.value(v)

mutable struct EventHandler{M<:Real, F<:Real}   # M: Real type of variables where Measurement can be used.
                                                # F: Floating-point type for variables that are no measurements

   # Input values for the event functions
    time::M                    # Current simulation time
    initial::Bool              # = true, if model is called at initialization
                               #         (if initial, event=true)
    terminal::Bool             # = true, if model is called for termination (close files, streams, visualization, ...)
    event::Bool                # = true, if model is called at an event
    afterSimulationStart::Bool # = true, if model is called after simulation start
    crossing::Bool             # = true, if model is called to compute crossing function

    # Computed by the event functions
    # For time events:
    maxTime::F              # Integrate at most up to maxTime
                            # (if integrateToEvent == true, maxTime = nextEventTime,
                            #  otherwise maxTime > nextEventTime)
    nextEventTime::M        # Next time event instant (= typemax(Time) if no event)
    integrateToEvent::Bool  # = true, if integrator shall integrate to nextEventTime
                            # = false, if integrator can integrate beyond nextEventTime
    newEventIteration::Bool   # = true, if another event iteration; = false, if no event iteration anymore
    firstEventIteration::Bool # = true, if first iteration at an event.

   # For state events:
    nz::Int                 # Number of event indicators
    z::Vector{M}            # Vector of event indicators (zero crossings). If one of z[i] passes
                            # zero, that is beforeEvent(z[i])*z[i] < 0, an event is triggered
                            # (allocated during instanciation according to nz).
    zPositive::Vector{Bool} # = true if z > 0 at the last event instant, otherwise false
    zDir::Vector{Cint}      # zDir[i] =  0: Root is reported for both crossing directions
                            #         =  1: Root is reported when crossing from negative to positive direction
                            #         = -1: Root is reported when crossing from positive to negative direction

    function EventHandler{M,F}(;nz::Int=0) where {M<:Real, F<:Real}
        @assert(nz >= 0)
        new(0, false, false, false, false, false, M(floatmax(F)), M(floatmax(F)),
            true, false, false, nz, ones(M,nz), fill(false, nz), fill(0, nz))
    end
end

const zEps = 1e-12

function edge!(h::EventHandler, nr::Int, crossing)::Bool
    returnValue = false

    if h.initial
        h.zPositive[nr] = crossing > 0.0
        h.zDir[nr] = 1
        returnValue = false

    elseif h.event
        new_zPositive = crossing > 0.0
        edge = !h.zPositive[nr] && new_zPositive
        h.zPositive[nr] = new_zPositive

        if edge
            h.newEventIteration = true
            returnValue = true
        end

    elseif h.crossing && h.zPositive[nr] && crossing < 0.0
        # Eventually need to improve later
        # (need to check, whether h.get_z is good enough to switch h.zPositive)
        h.zPositive[nr] = false
    end
    h.z[nr] = crossing + (h.zPositive[nr] ? zEps : -zEps)
    return returnValue
end

mutable struct Model{M<:Real, F<:Real}   # M: Real type of variables where Measurement can be used.
                                         # F: Floating-point type for variables that are no measurements
    # parameters
    m::M   # mass
    e::M   # coefficient of restitution
    g::F   # gravity (not used as Measurment)

	# variables
	h::M  # height
	v::M  # velocity
	a::M  # acceleration
	flying::Bool   # = true, if flying

    # simulation state
	x_start::Vector{M}
	eventHandler::EventHandler

    function Model{M,F}(m::M, e::M, g::F, h0::M) where {M<:Real, F<:Real}
	    @assert(m >= 0)
	    @assert(e >= 0)
		@assert(h0 > 0)
		x_start = M[h0, 0]
        new(m, e, g, 0, 0, 0, true, x_start, EventHandler{M,F}(nz=1))
    end
end

Model(;m::M=1.0, e::M=0.7, g::F=9.81, h0::M=1.0) where {M<:Real, F<:Real} = Model{M,F}(m,e,g,h0)

"Compute m::Model at actual time instant"
function equations!(m::Model{M,F}, t, x)::Nothing where {M,F}
    m.h = x[1]
    m.v = x[2]
	if edge!(m.eventHandler, 1, -m.h)
		m.v = -m.e*m.v
		if m.v < 0.01
			m.flying = false
			m.v = 0.0
		end
		x[2] = m.v   # re-initialize state vector x
	end
	m.a = m.flying ? -m.g : 0.0
    return nothing
end

"derivatives!: Called by integrator"
function derivatives!(der_x, x, m::Model, t)::Nothing
    equations!(m, t, x)
    der_x[1] = m.v
    der_x[2] = m.a
    return nothing
end

"outputs!: Called at every communication point"
function outputs!(x::Vector{M}, t, integrator)::Model where {M<:Real}
    m::Model = integrator.p
    equations!(m, t, x)
    return deepcopy(m)
end

function outputs2!(x, t, integrator)::Nothing
    # println("outputs!2 called at time = ", t)
    return nothing
end

"condition_neg!: Called to compute the zero crossings from positive to negative"
function conditions_neg!(z, x, t, integrator)::Nothing
    m::Model = integrator.p
    h = m.eventHandler
    h.crossing = true
    equations!(m, t, x)
    h.crossing = false
    z[1] = -h.z[1]
    #println("conditions! called at time = ", t, ", z = ", z[1])
    return nothing
end

"affect!: Called at an event instant (either as affect! or here as affect_neg!)"
function affect!(integrator, event_index)::Nothing
    println("affect! called at time = ", integrator.t, ", h = ", integrator.u[1])
    m::Model = integrator.p
    h = m.eventHandler
    h.event = true
    equations!(m, integrator.t, integrator.u)
    h.event = false
	set_proposed_dt!(integrator, 1e-5)
    return nothing
end

function x₀(m::Model{M,F}, t0::M)::Vector{M} where {M<:Real, F<:Real}
    # Initialize model
    return deepcopy(m.x_start)
end

t_start = 0.0
t_inc   = 0.0001
t_end   = 5.0
tspan   = (t_start, t_end)
tspan2  = t_start:t_inc:t_end

# Problem 1 with Float64
model1 = Model()
println("typeof(model1) = ", typeof(model1))

# x1₀ = deepcopy(model1.x_start)

saved_values1=SavedValues(typeof(model1.x_start[1]), Model)
cb1 = SavingCallback(outputs!, saved_values1, saveat=tspan2)
cb2 = VectorContinuousCallback(conditions_neg!, nothing, 1; affect_neg! = affect!, interp_points=10)
cb3 = FunctionCallingCallback(outputs2!, funcat=tspan2)

prob1 = ODEProblem(derivatives!, x₀, tspan, model1)
sol1  = solve(prob1, Tsit5(), reltol = 1e-6, saveat=t_inc, callback=CallbackSet(cb1,cb2,cb3))

println("Simulation successful")

t = saved_values1.t
m = saved_values1.saveval
flying = getfield.(m, :flying)

figure(1)
clf()
plot(sol1.t, getindex.(sol1.u, 1), label="h")
plot(t, flying, label="flying")
grid(true)
legend()
title("Bouncing ball with Float64")

println("length(sol1.t) = ", length(sol1.t))
println("length(t) = ", length(t))

# Problem 2 with Measurement
Float64(v::Measurement{Float64}) = Measurements.value(v)

model2 = Model(m=1.0 ± 0.1, e=0.7 ± 0.1, h0=1.0± 0.1, g=9.81)
#println("model2 = $model2")
#println("typeof(model2) = ", typeof(model2))
x2₀ = deepcopy(model2.x_start)

t_start = 0.0 ± 0.0
t_inc   = 0.0001 ± 0.0
t_end   = 5.0 ± 0.0
tspan   = (t_start, t_end)
tspan2  = t_start:t_inc:t_end

saved_values2=SavedValues(typeof(model2.x_start[1]), Model)
cb21 = SavingCallback(outputs!, saved_values2, saveat=tspan2)
cb22 = VectorContinuousCallback(conditions_neg!, nothing, 1; affect_neg! = affect!, interp_points=10)
prob2 = ODEProblem(derivatives!, x2₀, tspan, model2)
sol2  = solve(prob2, Tsit5(), reltol = 1e-6, saveat=t_inc, callback=CallbackSet(cb21,cb22))

h = getindex.(sol2.u, 1)
hm = Measurements.value.(h)
hu = Measurements.uncertainty.(h)
hmax = hm + hu
hmin = hm - hu

figure(2)
clf()
plot(sol2.t, hmax, label="\$h_{max}\$")
plot(sol2.t, hmin, label="\$h_{min}\$")
plot(sol2.t, h, label="\$h_{mean}\$")
grid(true)
legend()
title("Bouncing ball with Measurement{Float64}")
#@show sol2.t
end


text:
module Viral

using ..MultiFidelityABC
using Random
using StatsBase
using LinearAlgebra

# Viral model specification contains:
# - Model in viral_model_generator - what to simulate
# - mf-abc specification (what is low-fidelity, high-fidelity, etc) in viral_mfabc_problem
# - Prior parameter distribution in viral_prior

const nu = Float64.([0 1 0 -1 0 0 ; 1 -1 0 0 0 -1 ; 0 0 1 0 -1 -1; 0 0 0 0 0 1])

function propensity!(v::Array{Float64,1}, x::Array{Float64,1}, k::NTuple{6,Float64})
    v[:] = k.*[x[1], x[2], x[1], x[1], x[3], x[2]*x[3]]
    return nothing
end

const x0 = [1.0, 0.0, 0.0, 0.0]
const k_nominal = (1.0, 0.025, 100.0, 0.25, 1.9985, 7.5e-5)
const T = 200.0

gm = GillespieModel(nu, propensity!, x0, k_nominal, T)

############ Add in requirements for a hybrid model

const stochastic_reactions = [true,true,false,true,false,true]
function reduced_propensity!(v::Array{Float64,1}, x::Array{Float64,1}, k::NTuple{6,Float64})
    propensity!(v, x, k)
    v[:] .*= stochastic_reactions
    return nothing
end

function deterministic_step(t,x,k)
    err = (k[3]/k[5])*x[1] - x[3]
    if abs(1/err)<1
        tau = (-1/k[5]) * log(1 - abs(1/err))
    else
        tau = T-t
        err = 0
    end
    return tau, [0, 0, sign(err), 0]
end

hm = HybridModel(nu, propensity!, x0, k_nominal, T, reduced_propensity!, deterministic_step)
pop_size = 10

############ Set up a multifidelity problem

Random.seed!(123)
t_pop, x_pop = simulate(gm, pop_size) # Simulate the nominal model with a fixed seed
Random.seed!()

infected_threshold = 3

viral_count(x) = x[4,end]
infected(x) = (viral_count(x)>infected_threshold)
infection_time(x,t) = (infected(x) && t[searchsortedfirst(x[4,:], infected_threshold)]/T)

function summary_statistics(t_pop::Array{Times,1}, x_pop::Array{States,1})
    inf_flag = infected.(x_pop)
    if mean(inf_flag)>0
        return mean(inf_flag), log(2,mean(viral_count.(x_pop)[inf_flag])), mean(infection_time.(x_pop,t_pop)[inf_flag])
    else
        return 0.0, Inf64, Inf64
    end
end
yo = summary_statistics(t_pop, x_pop)

distance(y1, y2) = norm(y1.-y2)
distance(y) = distance(y, yo)

function lofi(k::Parameters)
    t_pop, x_pop, pp_pop = simulate(hm, k, pop_size)
    return distance(summary_statistics(t_pop, x_pop)), pp_pop
end

# The following hifi *couples* low fidelity and high fidelity simulations:
function hifi(k::Parameters, coupling_input::Array{PP,1})
    t_pop, x_pop = complete(hm, k, coupling_input)
    return distance(summary_statistics(t_pop, x_pop))
end

# # This hifi version would produce independent (uncoupled) high fidelity simulations:
# function hifi(k::Parameters, coupling_input::Array{PP,1})
#     t_pop, x_pop = simulate(gm, k, pop_size)
#     return distance(summary_statistics(t_pop,x_pop))
# end

using Distributions
function prior()::Parameters
    return map(ki->ki*(1.5^rand(Uniform(-1,1))), k_nominal) # Scale nominal parameters by factor between 2/3 and 3/2
end

mf_prob = MFABC(prior, lofi, hifi)

end

ext: jl

max_stars_repo_name: tpprescott/mf-abc

max_issues_repo_name: tpprescott/mf-abc

max_forks_repo_name: tpprescott/mf-abc


text:
# see utils/types/BlockTemplate
const MD_BLOCKS = Dict{Symbol,BlockTemplate}(e.opening => e for e in [
   BlockTemplate(:COMMENT,         :COMMENT_OPEN,   :COMMENT_CLOSE ),
   BlockTemplate(:RAW_HTML,        :RAW_HTML,       :RAW_HTML      ),
   BlockTemplate(:RAW_LATEX,       :RAW_LATEX,      :RAW_LATEX     ),
   BlockTemplate(:MD_DEF_BLOCK,    :MD_DEF_BLOCK,   :MD_DEF_BLOCK  ),
   #
   BlockTemplate(:EMPH_EM,        :EM_OPEN,        (:EM_CLOSE,        :EM_MX        ), nesting=true),
   BlockTemplate(:EMPH_EM,        :EM_MX,          (:EM_CLOSE,        :EM_MX        ), nesting=true),
   BlockTemplate(:EMPH_STRONG,    :STRONG_OPEN,    (:STRONG_CLOSE,    :STRONG_MX    ), nesting=true),
   BlockTemplate(:EMPH_STRONG,    :STRONG_MX,      (:STRONG_CLOSE,    :STRONG_MX    ), nesting=true),
   BlockTemplate(:EMPH_EM_STRONG, :EM_STRONG_OPEN, (:EM_STRONG_CLOSE, :EM_STRONG_MX ), nesting=true),
   BlockTemplate(:EMPH_EM_STRONG, :EM_STRONG_MX,   (:EM_STRONG_CLOSE, :EM_STRONG_MX ), nesting=true),
   #
   BlockTemplate(:AUTOLINK,        :AUTOLINK_OPEN,  :AUTOLINK_CLOSE ),
   # these blocks are disabled in find_blocks if they're not attached in
   # a link/img/... context
   BlockTemplate(:BRACKETS,    :BRACKET_OPEN,    :BRACKET_CLOSE,    nesting=true),
   BlockTemplate(:SQ_BRACKETS, :SQ_BRACKET_OPEN, :SQ_BRACKET_CLOSE, nesting=true),
   # code
   BlockTemplate(:CODE_BLOCK,      :CODE_PENTA,   :CODE_PENTA  ),
   BlockTemplate(:CODE_BLOCK,      :CODE_QUAD,    :CODE_QUAD   ),
   BlockTemplate(:CODE_BLOCK,      :CODE_TRIPLE,  :CODE_TRIPLE ),
   BlockTemplate(:CODE_INLINE,     :CODE_DOUBLE,  :CODE_DOUBLE ),
   BlockTemplate(:CODE_INLINE,     :CODE_SINGLE,  :CODE_SINGLE ),
   # maths
   BlockTemplate(:MATH_INLINE,  :MATH_INLINE,       :MATH_INLINE        ),
   BlockTemplate(:MATH_DISPL_A, :MATH_DISPL_A,      :MATH_DISPL_A       ),
   BlockTemplate(:MATH_DISPL_B, :MATH_DISPL_B_OPEN, :MATH_DISPL_B_CLOSE ),
   # md def one line
   BlockTemplate(:MD_DEF, :MD_DEF_OPEN, END_OF_LINE ),
   # div and braces
   BlockTemplate(:DIV,         :DIV_OPEN, :DIV_CLOSE,               nesting=true),
   BlockTemplate(:CU_BRACKETS, :CU_BRACKET_OPEN, :CU_BRACKET_CLOSE, nesting=true),
   # headers
   BlockTemplate(:H1, :H1_OPEN, END_OF_LINE ),
   BlockTemplate(:H2, :H2_OPEN, END_OF_LINE ),
   BlockTemplate(:H3, :H3_OPEN, END_OF_LINE ),
   BlockTemplate(:H4, :H4_OPEN, END_OF_LINE ),
   BlockTemplate(:H5, :H5_OPEN, END_OF_LINE ),
   BlockTemplate(:H6, :H6_OPEN, END_OF_LINE ),
   # Direct blocks
   SingleTokenBlockTemplate(:LINEBREAK ),
   SingleTokenBlockTemplate(:HRULE     ),
   # Direct blocks -- latex objects
   SingleTokenBlockTemplate(:LX_NEWENVIRONMENT ),
   SingleTokenBlockTemplate(:LX_NEWCOMMAND     ),
   SingleTokenBlockTemplate(:LX_COMMAND        ),
   SingleTokenBlockTemplate(:LX_BEGIN          ),
   SingleTokenBlockTemplate(:LX_END            )
   ])

#
# NOTE: {{...}} blocks (DBB_BLOCK) are processed separately, see find_blocks
#

const CAN_BE_LEFT_OPEN = (
    :EM_OPEN, :EM_MX,
    :STRONG_OPEN, :STRONG_MX,
    :EM_STRONG_OPEN, :EM_STRONG_MX,
    :BRACKET_OPEN,
    :AUTOLINK_OPEN
)

const MD_PASS0 = Dict{Symbol,BlockTemplate}(
   :RAW => BlockTemplate(:RAW, :RAW, :RAW),
)

# First pass: container blocks etc
const MD_PASS1_TEMPLATES = Dict{Symbol,BlockTemplate}(
   o => bt for (o, bt) in MD_BLOCKS
   if bt.name in (
         :COMMENT,
         :RAW_HTML, :RAW_LATEX,
         :MD_DEF_BLOCK, :MD_DEF,
         :CODE_BLOCK_LANG, :CODE_BLOCK!, :CODE_BLOCK, :CODE_INLINE,
         :MATH_INLINE, :MATH_DISPL_A, :MATH_DISPL_B,
         :DIV,
         :AUTOLINK,
         :CU_BRACKETS,
         :H1, :H2, :H3, :H4, :H5, :H6,
         :LX_BEGIN, :LX_END,
      )
)

# Second pass: links etc
const MD_PASS2_TEMPLATES = Dict{Symbol,BlockTemplate}(
   o => bt for (o, bt) in MD_BLOCKS
   if bt.name in (
         :BRACKETS,
         :SQ_BRACKETS
      )
)

# Last pass: the rest
const MD_PASS3_TEMPLATES = Dict{Symbol,BlockTemplate}(
   o => bt for (o, bt) in MD_BLOCKS
   if o ∉ union(keys(MD_PASS1_TEMPLATES), keys(MD_PASS2_TEMPLATES))
)


text:
# # Rational numbers
#
# In julia rational numbers can be constructed with the `//` operator.
# Lets define two rational numbers, `x` and `y`:

## Define variable x and y
x = 1//3
y = 2//5

# When adding `x` and `y` together we obtain a new rational number:

z = x + y

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
"""
    @sync ex

Run expression `ex` and synchronize the GPU afterwards.

See also: [`synchronize`](@ref).
"""
macro sync(code)
    quote
        local ret = $(esc(code))
        synchronize()
        ret
    end
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using BenchmarkTools
using Plots
gr(fmt=:png)
using CUDA

# %%
plotjulia(j) = heatmap(j; c=:gist_earth,
    size=(540, 540), colorbar=false, ticks=false, frame=false)

plotmandelbrot(m) = heatmap(m; c=reverse(cgrad(:jet1)),
    size=(540, 540), colorbar=false, ticks=false, frame=false)

function julia(z, c, maxiters=2^10, threshold_abs2=Inf)
    for i in Base.OneTo(maxiters)
        z = z * z + c
        abs2(z) ≥ threshold_abs2 && return i
    end
    maxiters + oneunit(maxiters)
end

mandelbrot(c, maxiters=2^10, threshold_abs2=Inf) =
    julia(zero(c), c, maxiters, threshold_abs2)

# %%
n = 2^9
x = range(-1.5, 1.5; length=n)
y = range(-1.5, 1.5; length=n)
z = complex.(x', y)
c = -0.380 + 0.605im

@show typeof(z)
@time j = julia.(z, c, 2^8)
@time j = julia.(z, c, 2^8)
@show typeof(j)
plotjulia(j)

# %%
n = 2^9
x = range(-1.5, 1.5; length=n)
y = range(-1.5, 1.5; length=n)
z = complex.(x', y)
c = -0.380 + 0.605im

c32 = ComplexF32(c)
z_cuda32 = cu(z)

@show typeof(z_cuda32)
@time j_cuda32 = julia.(z_cuda32, c32, Int32(2^8), Inf32)
@time j_cuda32 = julia.(z_cuda32, c32, Int32(2^8), Inf32)
@show typeof(j_cuda32)
@time cj_cuda32 = collect(j_cuda32)
plotjulia(cj_cuda32)

# %%
n = 2^9
x = range(-1.5, 1.5; length=n)
y = range(-1.5, 1.5; length=n)
z = complex.(x', y)
c = -0.380 + 0.605im

z_cuda = CuMatrix{ComplexF64}(z)

@show typeof(z_cuda)
@time j_cuda = julia.(z_cuda, c, 2^8)
@time j_cuda = julia.(z_cuda, c, 2^8)
@show typeof(j_cuda)
@time cj_cuda = collect(j_cuda)
plotjulia(cj_cuda)

# %%
n = 2^9
x = range(-1.6, 0.5; length=n)
y = range(-1.1, 1.1; length=n)
c = complex.(x', y)

@show typeof(c)
@time m = mandelbrot.(c, 2^6)
@time m = mandelbrot.(c, 2^6)
@show typeof(m)
plotmandelbrot(m)

# %%
n = 2^9
x = range(-1.6, 0.5; length=n)
y = range(-1.1, 1.1; length=n)
c = complex.(x', y)
c_cuda = CuMatrix{ComplexF64}(c)

@show typeof(c_cuda)
@time m_cuda = mandelbrot.(c_cuda, 2^6)
@time m_cuda = mandelbrot.(c_cuda, 2^6)
@show typeof(m_cuda)
@time cm_cuda = collect(m_cuda)
@show typeof(cm_cuda)
plotmandelbrot(cm_cuda)

# %%
n = 2^9
x = range(-0.714689, -0.714679; length=n)
y = range( 0.299872,  0.299882; length=n)
c = complex.(x', y)

@show typeof(c)
@time m = mandelbrot.(c, 2^12)
@time m = mandelbrot.(c, 2^12)
@show typeof(m)
plotmandelbrot(m)

# %%
n = 2^9
x = range(-0.714689, -0.714679; length=n)
y = range( 0.299872,  0.299882; length=n)
c = complex.(x', y)
c_cuda = CuMatrix{ComplexF64}(c)

@show typeof(c_cuda)
@time m_cuda = mandelbrot.(c_cuda, 2^12)
@time m_cuda = mandelbrot.(c_cuda, 2^12)
@show typeof(m_cuda)
@time cm_cuda = collect(m_cuda)
@show typeof(cm_cuda)
plotmandelbrot(cm_cuda)

# %%

max_stars_repo_path: 0023/Julia and Mandelbrot sets.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0023/Julia and Mandelbrot sets.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0023/Julia and Mandelbrot sets.jl

max_forks_repo_name: genkuroki/public


text:
export Clearance, AbstractCollisionModule, AbstractCollisionModule3D
export checkCoords, checkCoordTable

abstract type AbstractCollisionModule <: VirtualDevice end
abstract type AbstractCollisionModule3D <: AbstractCollisionModule end

abstract type CollisionModuleType end
struct PositionCollisionType <: CollisionModuleType end
struct PathCollisionType <: CollisionModuleType end

@mustimplement collisionModuleType(cm::AbstractCollisionModule)

# Robot Constants
const minClearance = 0.5Unitful.mm;

struct Clearance
  distance::typeof(1.0Unitful.mm)
  Clearance(distance) = distance < minClearance ? error("Clearance below minimum") : new(distance)
end

convert(::Type{Clearance}, x::Unitful.Length) = Clearance(x)

@mustimplement _checkCoordinate(cm::AbstractCollisionModule, pos::AbstractVector{<:Unitful.Length})

function init(cm::AbstractCollisionModule)
  @debug "Initializing CollisionModule with ID `$(cm.deviceID)`."
  cm.present = true
end

checkDependencies(cm::AbstractCollisionModule) = true

Base.close(cm::AbstractCollisionModule) = nothing

include("CollisionGeometries.jl")
include("SimpleBoreCollisionModule.jl")

function checkCoords(cm::AbstractCollisionModule, pos::AbstractVector{<:Unitful.Length}; returnVerbose=false)
  result_status, result_dist = _checkCoordinate(cm, pos)
  if !returnVerbose
    return all(result_status)
  end
  return result_status, result_dist
end

"if the function is called on a vector of collision modules return the logic AND of all decisions"
checkCoords(cms::Vector{<:AbstractCollisionModule}, args...) = .&([checkCoords(cm, args...; returnVerbose=false) for cm in cms]...)

checkCoords(cm::AbstractCollisionModule, pos::RobotCoords; kwargs...) = error("The CollisionModule is defined on scanner coordinates! Please make sure you are using the correct coordinate system")
checkCoords(cm::AbstractCollisionModule, pos::ScannerCoords; kwargs...) = checkCoords(cm, pos.data; kwargs...)

checkCoords(cm::AbstractCollisionModule3D, positions::Positions) = [checkCoords(cm, pos; returnVerbose=false) for pos in positions]

"""Convenience function to check an array of nx3 position vectors at once"""
function checkCoords(cm::AbstractCollisionModule3D, positions::AbstractMatrix{<:Unitful.Length})
  numPos, dim = size(positions);

  if dim != 3
    error("Only 3-dimensional coordinates accepted!")
  end

  table = Array{Bool}(undef, numPos)
  for (i, pos) in enumerate(eachrow(positions))
    table[i] = checkCoords(cm, pos, returnVerbose=false)
  end
  return table
end

### Functions to check multiple cooridnates in a more verbose manner

checkCoordTable(cm::AbstractCollisionModule3D, positions::Positions; kwargs...) = checkCoordTable(cm, permutedims(hcat(collect(positions)...)); kwargs...)

function checkCoordTable(cm::AbstractCollisionModule3D, positions::AbstractMatrix{<:Unitful.Length}; plotresults=false)
  numPos, dim = size(positions);

  if dim != 3
    error("Only 3-dimensional coordinates accepted!")
  end

  table = Array{Any}(undef, numPos, 6)
  for (i, pos) in enumerate(eachrow(positions))
    tmp = checkCoords(cm, pos, returnVerbose=true)
    table[i,:] .= vcat(tmp...)
  end

  table = hcat(collect(1:numPos), positions, table)
  # table headlines
  headline = ["#" "x" "y" "z" "Status x" "Status y" "Status z" "delta_x" "delta_y" "delta_z"]
  # create final table
  coordTable = vcat(headline, table)

  errorIndices = unique(getindex.(findall(x -> x == false, table[:,5:7]), 1)) # get row index of all coordinates where at least one dimension in invalid

  if isempty(errorIndices)
    display("All coordinates are safe!");
  else
    display("Used geometry: $(cm.params.objGeometry.name)")
    display("Used scanner diameter: $(cm.params.scannerDiameter)")
    display("Following coordinates are dangerous and NOT valid!");
    @info errorIndices
    # ERROR: BoundsError
    errorTable = coordTable[[1,(errorIndices .+ 1)...], :];
    display(errorTable)

    plotresults ? plotSafetyErrors(cm, errorIndices, positions) : display("Plotting not chosen...")
  end
  return coordTable
end


text:
incx = abs(incx)
@inbounds for i = ix:incx:ix-1+n*incx
    x[i] = $op( a, x[i] )
end


text:
"""
    Cat(dims)

Operator that performs a concatenation over the dimensions `dims`.

# Examples

```jldoctest
julia> c1 = MPIMapReduce.Cat(1);

julia> c1([1], [2])
2-element Vector{Int64}:
 1
 2

julia> c2 = MPIMapReduce.Cat(2);

julia> c2([1], [2])
1×2 Matrix{Int64}:
 1  2
```
"""
struct Cat{D} <: Function
    dims :: D
end
Cat(c::Cat) = c

dims(c::Cat) = c.dims
dims(::typeof(vcat)) = 1
dims(::typeof(hcat)) = 2

(c::Cat)(x...) = cat(x...; dims = dims(c))

_cat_shape(catdims, sizes, lengths) = Base.cat_shape(catdims, (), sizes...)
_cat_shape(catdims, ::Nothing, lengths) = Base.cat_shape(catdims, (), map(Tuple, lengths)...)

function _similar_cat(c::Cat, x, sizes, lengths)
    catdims = Base.dims2cat(dims(c))
    shape = _cat_shape(catdims, sizes, lengths)
    Base.cat_similar(x, eltype(x), shape)
end
_similar_cat(c::Cat, x::Number, sizes, lengths) = _similar_cat(c, [x], sizes, lengths)


text:
using TensorOperations
using KrylovKit
using LinearAlgebra

# MPS A-matrix is a 3-index tensor, A[i,s,j]
#    s
#    |
# i -A- j
#
# [s] acts on the local Hilbert space
# [i,j] act on the virtual bonds

# MPO W-matrix is a 4-index tensor, W[i,j,s,t]
#     s
#     |
#  i -W- j
#     |
#     t
#
# [s,t] act on the local Hilbert space,
# [i,j] act on the virtual bonds

function contract_from_left_MPS(A,E,B)
    @tensor temp_1[:] := E[-1,1] * A[1,-2,-3]
    @tensor temp[:] := temp_1[1,2,-1] * conj( B[1,2,-2] )
    return temp
end

## tensor contraction from the left hand side
## +-    +--A-
## |     |  |
## L' =  L--R-
## |     |  |
## +-    +--B-

function contract_from_left(E_L,X,W)
    @tensor temp_1[:] := E_L[1,-3,-4]  * X[1,-2,-1]
    @tensor temp_2[:] := temp_1[-1, -2, -5, 1 ]* conj( X[1, -4,-3] )
    @tensor temp[:] := temp_2[ -1, 2, -3, 3, 1 ] * W[1,-2, 2, 3]
    return temp
end

## tensor contraction from the right hand side
##  -+     -A--+
##   |      |  |
##  -R' =  -W--R
##   |      |  |
##  -+     -B--+
function contract_from_right(E_R,Y,W)
    @tensor temp_1[:] :=  Y[-1,-2,1]  * E_R[ 1, -3,-4]
    @tensor temp_2[:] := temp_1[ -1, -2, -5,  1] * conj( Y[-3,-4,1] )
    @tensor temp[:] := temp_2[ -1, 1, -3, 2, 6] * W[ -2 6 1 2 ]
    return temp
end

max_stars_repo_name: aanniikkeettr/Tensor-Networks

max_issues_repo_name: aanniikkeettr/Tensor-Networks

max_forks_repo_name: aanniikkeettr/Tensor-Networks


text:
function rk4(f::Function, x₀::Float64, y₀::Float64, x₁::Float64, n)
    vx = Vector{Float64}(n + 1)
    vy = Vector{Float64}(n + 1)
    vx[1] = x = x₀
    vy[1] = y = y₀
    h = (x₁ - x₀) / n
    for i in 1:n
        k₁ = h * f(x, y)
        k₂ = h * f(x + 0.5h, y + 0.5k₁)
        k₃ = h * f(x + 0.5h, y + 0.5k₂)
        k₄ = h * f(x + h, y + k₃)
        vx[i + 1] = x = x₀ + i * h
        vy[i + 1] = y = y + (k₁ + 2k₂ + 2k₃ + k₄) / 6
    end
    return vx, vy
end

vx, vy = rk4(f, 0.0, 1.0, 10.0, 100)
for (x, y) in Iterators.take(zip(vx, vy), 10)
    @printf("%4.1f %10.5f %+12.4e\n", x, y, y - theoric(x))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
Code for writing for CHGCAR file containing 3D grid data
=#

module Chgcar

using Printf
using LinearAlgebra

"""Enumerate fortran float where a leading zero is written"""
function ffloat(f)
    if f == 0.
        return @sprintf "0.000000000"
    end
    string = @sprintf "%.10E" f
    if f > 0
        return @sprintf "0.%s%sE%+03d" string[1:1] string[3:12] Int(ceil(log10(f)))
    end
    return @sprintf "-0.%s%sE%+03d" string[2:2] string[4:13]  Int(ceil(log10(-f)))
end

"""
Write out a CHGCAR file with
"""
function write_chgcar(io::IO, datasets; lattice_col=Diagonal([10, 10, 10]),
                      species=(:Si, ), coords=[0., 0., 0.][:, :], comment="VolumetricData")
    write(io, "$(comment)\n")
    write(io, "   1.00000000000000\n")
    write(io, @sprintf(" %12.6f%12.6f%12.6f\n", lattice_col[:, 1]...))
    write(io, @sprintf(" %12.6f%12.6f%12.6f\n", lattice_col[:, 2]...))
    write(io, @sprintf(" %12.6f%12.6f%12.6f\n", lattice_col[:, 3]...))
    # Symbols and counts
    usymbols = unique(species)
    counts = zeros(Int, length(usymbols))
    for (i, symbol) in enumerate(unique(species))
        write(io, @sprintf "%5s" symbol)
        counts[i] = count(x -> x == symbol, species)
    end
    write(io, "\n")
    for c in counts
        write(io, @sprintf "%6d" c)
    end
    write(io, "\n")
    # Write the coordinates
    write(io, "Direct\n")
    for coord in eachcol(coords)
        @sprintf "%10.6f%10.6f%10.6f\n" coord[1] coord[2] coord[3]
    end
    write(io, " \n")
    # Now starts the volumetric data set
    for data in datasets
        write_data(io, data)
    end
end

"""
Write a CHGCAR file
"""

function write_chgcar(fname::AbstractString, datasets; lattice_col=Diagonal([10, 10, 10]),
                      species=(:Si, ), coords=[0., 0., 0.][:, :], comment="VolumetricData")
    open(fname, "w") do io
        write_chgcar(io, datasets; lattice_col, species, coords, comment)
    end
end

"""
Write the volumetric data block of the CHGCAR file
"""
function write_data(io::IO, data)

    ngx, ngy, ngz = size(data)
    write(io, "   $ngx   $ngy   $ngz\n")
    count = 1
    for k=1:ngz, j=1:ngy, i=1:ngx
        # Start of a new line
        count % 5 == 1 && write(io, " ")
        # Write data
        if count % 5 == 0
            write(io, ffloat(data[i, j, k]), "\n")
        else
            write(io, ffloat(data[i, j, k]), " ")
        end
        count += 1
    end
    # Write a final line break
    if count % 5 != 1
        write(io, "\n")
    end

end

export write_chgcar
end # module chgcar

import .Chgcar: write_chgcar

"""
"""
function write_chgcar(fname::AbstractString, density::ChargeDensity;sum_and_spin=true)
    if size(density.density, 4) == 2
        if sum_and_spin
        total = (density.density[:, :, :, 1] .+ density.density[:, :, :, 2]) ./ 2
        diff = (density.density[:, :, :, 1] .- density.density[:, :, :, 2]) ./ 2
        datasets = (total, diff)
        else
            datasets = [density.density[:, :, :, 1], density.density[:, :, :, 2]]
        end
    else
        datasets = [density.density[:, :, :, 1],]
    end

    write_chgcar(fname, datasets;lattice_col=density.basis, species=(:Si, ), coords=[0., 0., 0.][:, :], )
end


text:
# Outline:
# * Abstract types
#   - AbstractAnatomy
#   - PositionedAnatomy
#     - positiontype
#   - PositionedStructure
#     - structuretype
#   - WeightedPositionedStructure
#     - weighttype
# * Structure types
#   - AnatomicalStructure
#   - BilateralStructure
#   - PointStructure
#   - WeightedPointStructure
# * Link types
#   - SimpleLink
#   - StructureLink
#   - WeightedStructureLink
#   - StructureLinkPath
#   - WeightedStructureLinkPath
# * Hierarchical types
#
# - TODO should weights be the field for colors or should that be a separate option.
#

###
### Abstract types
###

"""
    AbstractAnatomy

Supertype for representing anatomical information.
"""
abstract type AbstractAnatomy end

"""
    PositionedAnatomy

An `AbstractAnatomy` subtype that contains positional information.
"""
abstract type PositionedAnatomy{P} <: AbstractAnatomy end

"""
    positioneltype(::PositionedAnatomy)

Returns the underlying type of the position in an instance of
`PositionedAnatomy`.
"""
positiontype(::PA) where {PA<:PositionedAnatomy} = positiontype(PA)
positiontype(::Type{<:PositionedAnatomy{P}}) where {P} = P

"""
    PositionedStructure{S,P}

An `PositionedAnatomy` subtype that contains positional and structural
information.
"""
abstract type PositionedStructure{S,P} <: PositionedAnatomy{S,P} end

"""
    structuretype(::PositionedStructure)

Returns the underlying structure type in a `PositionedStructure`.
"""
structuretype(::P) where {P<:PositionedAnatomy} = structuretype(P)
structuretype(::Type{<:PositionedAnatomy{S}}) where {S} = S

"""
    WeightedPositionedStructure{W<:Number,S,P}

A `PositionedStructure` subtype that contains a weight in addition to
positional and structural information.
"""
abstract type WeightedPositionedStructure{W,S,P} <: PositionedAnatomy{S,P} end

weighttype(::WPS) where {WPS<:WeightedPositionedStructure} = weighttype(WPS)
weighttype(::Type{<:WeightedPositionedStructure{W}}) where {W} = W

###
### Point/node/vertex structs
###

"""
    AnatomicalStructure

Represents an  anatomical structure.
"""
struct AnatomicalStructure <: AbstractAnatomy
    name::Symbol
end

AnatomicalStructure(n::String) = AnatomicalStructure(Symbol(n))

Base.names(s::AnatomicalStructure) = s.name

"""
    BilateralStructure
"""

struct BilateralStructure{S<:AnatomicalStructure} <: AbstractAnatomy
    structure::S
end

structure(s::BilateralStructure) = getproperty(s, :structure)

Base.names(s::) = names(structure(s))

"""
    PointStructure
"""
struct PointStructure{S<:AbstractAnatomy,T,N} <: PositionedStructure{S,Point{N,T}}
    structure::S
    position::Point{N,T}
end

position(s::PositionedStructure) = getproperty(s, :position)
structure(s::PositionedStructure) = getproperty(s, :structure)

"""
    WeightedPointStructure
"""
struct WeightedPointStructure{S<:AbstractAnatomy,T,N} <: WeightedPositionedStructure{W,S,Point{N,T}}
    structure::S
    position::Point{N,T}
    weight::W
end

position(s::WeightedPointStructure) = getproperty(s, :position)
structure(s::WeightedPointStructure) = getproperty(s, :structure)
weight(s::WeightedPointStructure) = getproperty(s, :weight)

###
### Link/edge structs
###

"""
    SimpleLink

A link containing only positional information that chains together anatomical
structures.
"""
struct SimpleLink{T,N} <: PositionedAnatomy{Point{N,T}}
    src::Point{N,T}
    dst::Point{N,T}
end

src(s::SimpleLink) = getproperty(s, :src)
dst(s::SimpleLink) = getproperty(s, :dst)

"""
    StructureLink

A link associating two distinct `PositionedStructures`.
"""
struct StructureLink{S,T,N,Src<:PositionedStructure{S,Point{N,T}},Dst<:PositionedStructure{S,Point{N,T}}} <: PositionedStructure{S,Point{N,T}}
    src::Src
    dst::Dst
end

src(s::StructureLink) = getproperty(s, :src)
dst(s::StructureLink) = getproperty(s, :dst)

"""
    WeightedStructureLink{W,S,T}
"""
struct WeightedStructureLink{W,T,N,Src<:PositionedStructure{S,Point{N,T}},Dst<:PositionedStructure{S,Point{N,T}}} <: WeightedPositionedStructure{S,Point{N,T}}
    src::Src
    dst::Dst
    weight::W
end

src(s::WeightedStructureLink) = getproperty(s, :src)
dst(s::WeightedStructureLink) = getproperty(s, :dst)
weight(s::WeightedStructureLink) = getproperty(s, :weight)

"""
    StructureLinkPath

A link associating two distinct `PositionedStructures` through series of
`SimpleLink` types.
"""
struct StructureLinkPath{S,T,Src<:PositionedStructure{S,Point{N,T}},Dst<:PositionedStructure{S,Point{N,T}}} <: PositionedStructure{S,Point{N,T}}
    src::Src
    dst::Dst
    pathway::Vector{SimpleLink{T}}
end

src(s::StructureLink) = getproperty(s, :src)
dst(s::StructureLink) = getproperty(s, :dst)

"""
    WeightedStructureLinkPath{W,S,T}

A link associating two distinct `PositionedStructures` through series of
`SimpleLink` types with accompanying weights.
"""
struct WeightedStructureLinkPath{W,T,N,Src<:PositionedStructure{S,Point{N,T}},Dst<:PositionedStructure{S,Point{N,T}}} <: WeightedPositionedStructure{S,Point{N,T}}
    src::Src
    dst::Dst
    pathway::Vector{SimpleLink{T}}
    weight::Vector{W}
end

src(s::WeightedStructureLink) = getproperty(s, :src)
dst(s::WeightedStructureLink) = getproperty(s, :dst)
weight(s::WeightedStructureLink) = getproperty(s, :weight)

const LinkeType = Union{StructureLink,StructureLinkPath,WeightedStructureLink,WeightedStructureLinkPath}

###
### Hierarchical links
###

"""
    HierarchicalLink
"""
Base.@kwdef struct HierarchicalLink{S<:AbstractAnatomy,P,C} <: AbstractAnatomy
    structure::S
    parent::P=nothing
    children::C=HierarchicalLink[]
end

const TopHierarchicalLink{S,C<:Vector} = HierarchicalLink{S,Nothing,C}
const BottomHierarchicalLink{S,P<:AbstractAnatomy} = HierarchicalLink{S,P,Nothing}

Base.parent(s::HierarchicalLink) = getproperty(s, :parent)

AbstractTrees.children(s::HierarchicalLink) = getproperty(s, :children)

Base.getindex(v::HierarchicalLink, i::Int) = children(v)[i]

# TODO find function for searching tree structure
function addchild!(tree:HierarchicalLink, parent::AbstractAnatomy, child::AbstractAnatomy)
end

function addchild!(parent::HierarchicalLink, child::AbstractAnatomy)
    push!(parent.children, HierarchicalLink(child, parent, HierarchicalLink[]))
end

function addchild!(parent::HierarchicalLink, child::Vector)
    for childe_i in child
        addchild!(parent, child)
    end
end

"""
    nchildren(::HierarchicalLink)

Returns the number of children that belong to a node.
"""
nchildren(s::HierarchicalLink) = length(children(s))

function AbstractTrees.printnode(io::IO, v::AbstractAnatomicalStructure)
    print(io, "$(typeof(v).name)($(names(v)))")
end

function Base.show(io::IO, v::AbstractAnatomicalStructure)
    print_tree(v)
end

#=
"""
    MetaStructure{C,M}
"""
struct MetaStructure{V<:AbstractAnatomicalStructure,M} <: AbstractAnatomicalStructure
    node::V
    meta::M
end

"""
    PointStructure
"""
struct PointStructure{V<:AbstractAnatomicalStructure,T<:Real,N} <: AbstractAnatomicalStructure
    vertex::V
    point::NTuple{N,T}
end

struct AnatomicalEdge
    parent::P
    child::C
end

struct SiblingEdge{W,S,D} <: AbstractEdge{W}
    src::S
    dst::D
    weight::W
end

struct ParentChildEdge{

weight(e::SiblingEdge) = e.weight
src(e::SiblingEdge) = e.src
dst(e::SiblingEdge) = e.dst

"""
    AnatomicalSiblings
"""
struct AnatomicalSiblings{names,T}
    siblings::NamedTuple{names,T}
end

siblings(x::AnatomicalSiblings) = getproperty(x, :siblings)

function Base.getproperty(x::AnatomicalSiblings, p::Symbol)
    getproperty(siblings(x), p)
end

Base.getindex(x::AnatomicalSiblings, i::Int) = getindex(siblings(x), i)
Base.getindex(x::AnatomicalSiblings, i::Symbol) = getproperty(x, i)

"""
    substructuretype(s::AbstractAnatomicalStructure) -> Type

Returns the type of the structure immediately wrapped by a given anaatomical
structure type. If the given type doesn't wrap anything it returns the type of
the original arugment.

# Examples
```
hippocampus = BilateralStructure(:hippocampus)

substructuretype(hippocampus)
```
"""
substructuretype(::S) where {S<:AbstractAnatomicalStructure} = S

=#


text:
# Autogenerated wrapper script for Notcurses_jll for x86_64-unknown-freebsd
export libnotcurses, libnotcurses_core, libnotcurses_ffi

using Ncurses_jll
using libunistring_jll
using libdeflate_jll
JLLWrappers.@generate_wrapper_header("Notcurses")
JLLWrappers.@declare_library_product(libnotcurses, "libnotcurses.so.3")
JLLWrappers.@declare_library_product(libnotcurses_core, "libnotcurses-core.so.3")
JLLWrappers.@declare_library_product(libnotcurses_ffi, "libnotcurses-ffi.so.3")
function __init__()
    JLLWrappers.@generate_init_header(Ncurses_jll, libunistring_jll, libdeflate_jll)
    JLLWrappers.@init_library_product(
        libnotcurses,
        "lib/libnotcurses.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnotcurses_core,
        "lib/libnotcurses-core.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnotcurses_ffi,
        "lib/libnotcurses-ffi.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Extend Zygote to work with static vectors and custom types on the GPU
# Here be dragons

using ForwardDiff: Chunk, Dual, dualize, partials, value
using Zygote: unbroadcast

Zygote.accum(x::AbstractArray{<:SizedVector}, ys::AbstractArray{<:SVector}...) = Zygote.accum.(convert(typeof(ys[1]), x), ys...)
Zygote.accum(x::AbstractArray{<:SVector}, ys::AbstractArray{<:SizedVector}...) = Zygote.accum.(x, convert.(typeof(x), ys)...)

Base.:+(x::Real, y::SizedVector) = x .+ y
Base.:+(x::SizedVector, y::Real) = x .+ y

Base.:+(x::Real, y::Zygote.OneElement) = x .+ y
Base.:+(x::Zygote.OneElement, y::Real) = x .+ y

function Base.:+(x::Atom{T, T, T, T}, y::Atom{T, T, T, T}) where T
    Atom{T, T, T, T}(0, x.charge + y.charge, x.mass + y.mass, x.σ + y.σ, x.ϵ + y.ϵ)
end

function Base.:-(x::Atom{T, T, T, T}, y::Atom{T, T, T, T}) where T
    Atom{T, T, T, T}(0, x.charge - y.charge, x.mass - y.mass, x.σ - y.σ, x.ϵ - y.ϵ)
end

function Zygote.accum(x::LennardJones{S, C, W, F, E}, y::LennardJones{S, C, W, F, E}) where {S, C, W, F, E}
    LennardJones{S, C, W, F, E}(x.cutoff, x.nl_only, x.lorentz_mixing, x.weight_14 + y.weight_14, x.force_units, x.energy_units)
end

function Zygote.accum(x::CoulombReactionField{D, S, W, T, F, E}, y::CoulombReactionField{D, S, W, T, F, E}) where {D, S, W, T, F, E}
    CoulombReactionField{D, S, W, T, F, E}(x.dist_cutoff + y.dist_cutoff, x.solvent_dielectric + y.solvent_dielectric, x.nl_only,
                x.weight_14 + y.weight_14, x.coulomb_const + y.coulomb_const, x.force_units, x.energy_units)
end

function Zygote.accum_sum(xs::AbstractArray{Tuple{LennardJones{S, C, W, F, E}, CoulombReactionField{D, SO, W, T, F, E}}}; dims=:) where {S, C, W, F, E, D, SO, T}
    reduce(Zygote.accum, xs, dims=dims; init=(
        LennardJones{S, C, W, F, E}(nothing, false, false, zero(W), NoUnits, NoUnits),
        CoulombReactionField{D, SO, W, T, F, E}(zero(D), zero(S), false, zero(W), zero(T), NoUnits, NoUnits),
    ))
end

function Zygote.accum(x::Tuple{NTuple{N, Int}, NTuple{N, T}, NTuple{N, E}},
                        y::Tuple{NTuple{N, Int}, NTuple{N, T}, NTuple{N, E}}) where {N, T, E}
    ntuple(n -> 0, N), x[2] .+ y[2], x[3] .+ y[3]
end

Base.zero(::Type{Atom{T, T, T, T}}) where {T} = Atom(0, zero(T), zero(T), zero(T), zero(T))
atom_or_empty(at::Atom, T) = at
atom_or_empty(at::Nothing, T) = zero(Atom{T, T, T, T})

Zygote.z2d(dx::AbstractArray{Union{Nothing, Atom{T, T, T, T}}}, primal::AbstractArray{Atom{T, T, T, T}}) where {T} = atom_or_empty.(dx, T)
Zygote.z2d(dx::SVector{3, T}, primal::T) where {T} = sum(dx)

Zygote.unbroadcast(x::Tuple{Any}, x̄::Nothing) = nothing

function Zygote.unbroadcast(x::AbstractArray{<:Real}, x̄::AbstractArray{<:StaticVector})
    if length(x) == length(x̄)
        Zygote._project(x, sum.(x̄))
    else
        dims = ntuple(d -> size(x, d) == 1 ? d : ndims(x̄) + 1, ndims(x̄))
        Zygote._project(x, accum_sum(x̄; dims=dims))
    end
end

Zygote._zero(xs::AbstractArray{<:StaticVector}, T) = fill!(similar(xs, T), zero(T))

function Zygote._zero(xs::AbstractArray{Atom{T, T, T, T}}, ::Type{Atom{T, T, T, T}}) where {T}
    fill!(similar(xs), Atom{T, T, T, T}(0, zero(T), zero(T), zero(T), zero(T)))
end

function Base.zero(::Type{Union{Nothing, SizedVector{D, T, Vector{T}}}}) where {D, T}
    zero(SizedVector{D, T, Vector{T}})
end

# Slower version than in Zygote but doesn't give wrong gradients on the GPU for repeated indices
# Here we just move it to the CPU then move it back
# See https://github.com/FluxML/Zygote.jl/pull/1131
Zygote.∇getindex(x::CuArray, inds::Tuple{AbstractArray{<:Integer}}) = dy -> begin
    inds1_cpu = Array(inds[1])
    dx = Zygote._zero(Array(x), eltype(dy))
    dxv = view(dx, inds1_cpu)
    dxv .= Zygote.accum.(dxv, Zygote._droplike(Array(dy), dxv))
    return Zygote._project(x, cu(dx)), nothing
end

# Extend to add extra empty partials before (B) and after (A) the SVector partials
@generated function ForwardDiff.dualize(::Type{T}, x::StaticArray, ::Val{B}, ::Val{A}) where {T, B, A}
    N = length(x)
    dx = Expr(:tuple, [:(Dual{T}(x[$i], chunk, Val{$i + $B}())) for i in 1:N]...)
    V = StaticArrays.similar_type(x, Dual{T, eltype(x), N + B + A})
    return quote
        chunk = Chunk{$N + $B + $A}()
        $(Expr(:meta, :inline))
        return $V($(dx))
    end
end

# Dualize a value with extra partials
macro dualize(x, n_partials::Integer, active_partial::Integer)
    ps = [i == active_partial for i in 1:n_partials]
    return :(ForwardDiff.Dual($(esc(x)), $(ps...)))
end

# Space for 4 duals given to interactions though only one used in this case
# No gradient for cutoff type
function dualize_fb(inter::LennardJones{S, C, W, F, E}) where {S, C, W, F, E}
    w14 = inter.weight_14
    dual_weight_14 = @dualize(w14, 21, 1)
    return LennardJones{S, C, typeof(dual_weight_14), F, E}(inter.cutoff, inter.nl_only,
                inter.lorentz_mixing, dual_weight_14, inter.force_units, inter.energy_units)
end

function dualize_fb(inter::CoulombReactionField{D, S, W, T, F, E}) where {D, S, W, T, F, E}
    dc, sd, w14, cc = inter.dist_cutoff, inter.solvent_dielectric, inter.weight_14, inter.coulomb_const
    dual_dist_cutoff        = @dualize(dc , 21, 1)
    dual_solvent_dielectric = @dualize(sd , 21, 2)
    dual_weight_14          = @dualize(w14, 21, 3)
    dual_coulomb_const      = @dualize(cc , 21, 4)
    return CoulombReactionField{typeof(dual_dist_cutoff), typeof(dual_solvent_dielectric), typeof(dual_weight_14), typeof(dual_coulomb_const), F, E}(
                                dual_dist_cutoff, dual_solvent_dielectric, inter.nl_only, dual_weight_14,
                                dual_coulomb_const, inter.force_units, inter.energy_units)
end

function dualize_fb(inter::HarmonicBond{D, K}) where {D, K}
    b0, kb = inter.b0, inter.kb
    dual_b0 = @dualize(b0, 11, 1)
    dual_kb = @dualize(kb, 11, 2)
    return HarmonicBond{typeof(dual_b0), typeof(dual_kb)}(dual_b0, dual_kb)
end

function dualize_fb(inter::HarmonicAngle{D, K}) where {D, K}
    th0, cth = inter.th0, inter.cth
    dual_th0 = @dualize(th0, 14, 1)
    dual_cth = @dualize(cth, 14, 2)
    return HarmonicAngle{typeof(dual_th0), typeof(dual_cth)}(dual_th0, dual_cth)
end

function dualize_fb(inter::PeriodicTorsion{6, T, E}) where {T, E}
    p1, p2, p3, p4, p5, p6 = inter.phases
    k1, k2, k3, k4, k5, k6 = inter.ks
    dual_phases = (
        @dualize(p1, 27,  1), @dualize(p2, 27,  2), @dualize(p3, 27,  3),
        @dualize(p4, 27,  4), @dualize(p5, 27,  5), @dualize(p6, 27,  6),
    )
    dual_ks = (
        @dualize(k1, 27,  7), @dualize(k2, 27,  8), @dualize(k3, 27,  9),
        @dualize(k4, 27, 10), @dualize(k5, 27, 11), @dualize(k6, 27, 12),
    )
    return PeriodicTorsion{6, eltype(dual_phases), eltype(dual_ks)}(inter.periodicities,
                            dual_phases, dual_ks)
end

function dualize_atom_fb1(at::Atom)
    charge, mass, σ, ϵ = at.charge, at.mass, at.σ, at.ϵ
    dual_charge = @dualize(charge, 21, 11)
    dual_mass = @dualize(mass, 21, 12)
    dual_σ = @dualize(σ, 21, 13)
    dual_ϵ = @dualize(ϵ, 21, 14)
    return Atom{typeof(dual_charge), typeof(dual_mass), typeof(dual_σ), typeof(dual_ϵ)}(
                at.index, dual_charge, dual_mass, dual_σ, dual_ϵ)
end

function dualize_atom_fb2(at::Atom)
    charge, mass, σ, ϵ = at.charge, at.mass, at.σ, at.ϵ
    dual_charge = @dualize(charge, 21, 15)
    dual_mass = @dualize(mass, 21, 16)
    dual_σ = @dualize(σ, 21, 17)
    dual_ϵ = @dualize(ϵ, 21, 18)
    return Atom{typeof(dual_charge), typeof(dual_mass), typeof(dual_σ), typeof(dual_ϵ)}(
                at.index, dual_charge, dual_mass, dual_σ, dual_ϵ)
end

function dual_function_svec(f::F) where F
    function (arg1)
        ds1 = dualize(Nothing, arg1, Val(0), Val(0))
        return f(ds1)
    end
end

function dual_function_svec_real(f::F) where F
    function (arg1::SVector{D, T}, arg2) where {D, T}
        ds1 = dualize(Nothing, arg1, Val(0), Val(1))
        # Leaving the integer type in here results in Float32 -> Float64 conversion
        ds2 = Zygote.dual(isa(arg2, Int) ? T(arg2) : arg2, (false, false, false, true))
        return f(ds1, ds2)
    end
end

function dual_function_svec_svec(f::F) where F
    function (arg1, arg2)
        ds1 = dualize(Nothing, arg1, Val(0), Val(3))
        ds2 = dualize(Nothing, arg2, Val(3), Val(0))
        return f(ds1, ds2)
    end
end

function dual_function_atom(f::F) where F
    function (arg1)
        charge, mass, σ, ϵ = arg1.charge, arg1.mass, arg1.σ, arg1.ϵ
        ds1 = Atom(arg1.index, @dualize(charge, 4, 1), @dualize(mass, 4, 2),
                    @dualize(σ, 4, 3), @dualize(ϵ, 4, 4))
        return f(ds1)
    end
end

function dual_function_force_broadcast(f::F) where F
    function (arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
        ds1 = (dualize_fb(arg1[1]), dualize_fb(arg1[2])) # Only works for this case
        ds2 = dualize(Nothing, arg2, Val(4), Val(14))
        ds3 = dualize(Nothing, arg3, Val(7), Val(11))
        ds4 = dualize_atom_fb1(arg4)
        ds5 = dualize_atom_fb2(arg5)
        ds6 = dualize(Nothing, arg6, Val(18), Val(0))
        ds7 = arg7
        ds8 = arg8
        return f(ds1, ds2, ds3, ds4, ds5, ds6, ds7, ds8)
    end
end

function dual_function_specific_2_atoms(f::F) where F
    function (arg1, arg2, arg3, arg4)
        ds1 = dualize_fb(arg1)
        ds2 = dualize(Nothing, arg2, Val(2), Val(6))
        ds3 = dualize(Nothing, arg3, Val(5), Val(3))
        ds4 = dualize(Nothing, arg4, Val(8), Val(0))
        return f(ds1, ds2, ds3, ds4)
    end
end

function dual_function_specific_3_atoms(f::F) where F
    function (arg1, arg2, arg3, arg4, arg5)
        ds1 = dualize_fb(arg1)
        ds2 = dualize(Nothing, arg2, Val( 2), Val(9))
        ds3 = dualize(Nothing, arg3, Val( 5), Val(6))
        ds4 = dualize(Nothing, arg4, Val( 8), Val(3))
        ds5 = dualize(Nothing, arg5, Val(11), Val(0))
        return f(ds1, ds2, ds3, ds4, ds5)
    end
end

function dual_function_specific_4_atoms(f::F) where F
    function (arg1, arg2, arg3, arg4, arg5, arg6)
        ds1 = dualize_fb(arg1)
        ds2 = dualize(Nothing, arg2, Val(12), Val(12))
        ds3 = dualize(Nothing, arg3, Val(15), Val( 9))
        ds4 = dualize(Nothing, arg4, Val(18), Val( 6))
        ds5 = dualize(Nothing, arg5, Val(21), Val( 3))
        ds6 = dualize(Nothing, arg6, Val(24), Val( 0))
        return f(ds1, ds2, ds3, ds4, ds5, ds6)
    end
end

@inline function sum_partials(sv::SVector{3, Dual{Nothing, T, P}}, y1, i::Integer) where {T, P}
    partials(sv[1], i) * y1[1] + partials(sv[2], i) * y1[2] + partials(sv[3], i) * y1[3]
end

@inline function Zygote.broadcast_forward(f, arg1::AbstractArray{SVector{D, T}}) where {D, T}
    out = dual_function_svec(f).(arg1)
    y = map(x -> value.(x), out)
    function bc_fwd_back(ȳ_in)
        ȳ = arg1 isa CuArray ? cu(ȳ_in) : ȳ_in
        barg1 = broadcast(ȳ, out) do y1, o1
            if length(y1) == 1
                y1 .* SVector{D, T}(partials(o1))
            else
                SVector{D, T}(sum_partials(o1, y1, 1), sum_partials(o1, y1, 2), sum_partials(o1, y1, 3))
            end
        end
        darg1 = unbroadcast(arg1, barg1)
        (nothing, nothing, darg1)
    end
    return y, bc_fwd_back
end

@inline function Zygote.broadcast_forward(f, arg1::AbstractArray{SVector{D, T}}, arg2) where {D, T}
    out = dual_function_svec_real(f).(arg1, arg2)
    y = map(x -> value.(x), out)
    function bc_fwd_back(ȳ_in)
        ȳ = arg1 isa CuArray ? cu(ȳ_in) : ȳ_in
        barg1 = broadcast(ȳ, out) do y1, o1
            if length(y1) == 1
                y1 .* SVector{D, T}(partials.((o1,), (1, 2, 3)))
            else
                SVector{D, T}(sum_partials(o1, y1, 1), sum_partials(o1, y1, 2), sum_partials(o1, y1, 3))
            end
        end
        darg1 = unbroadcast(arg1, barg1)
        darg2 = unbroadcast(arg2, broadcast((y1, o1) -> y1 .* partials.(o1, 4), ȳ, out))
        (nothing, nothing, darg1, darg2)
    end
    return y, bc_fwd_back
end

@inline function Zygote.broadcast_forward(f, arg1::AbstractArray{SVector{D, T}}, arg2::AbstractArray{SVector{D, T}}) where {D, T}
    out = dual_function_svec_svec(f).(arg1, arg2)
    y = map(x -> value.(x), out)
    function bc_fwd_back(ȳ_in)
        ȳ = arg1 isa CuArray ? cu(ȳ_in) : ȳ_in
        barg1 = broadcast(ȳ, out) do y1, o1
            if length(y1) == 1
                y1 .* SVector{D, T}(partials.((o1,), (1, 2, 3)))
            else
                SVector{D, T}(sum_partials(o1, y1, 1), sum_partials(o1, y1, 2), sum_partials(o1, y1, 3))
            end
        end
        darg1 = unbroadcast(arg1, barg1)
        barg2 = broadcast(ȳ, out) do y1, o1
            if length(y1) == 1
                y1 .* SVector{D, T}(partials.((o1,), (4, 5, 6)))
            else
                SVector{D, T}(sum_partials(o1, y1, 4), sum_partials(o1, y1, 5), sum_partials(o1, y1, 6))
            end
        end
        darg2 = unbroadcast(arg2, barg2)
        (nothing, nothing, darg1, darg2)
    end
    return y, bc_fwd_back
end

@inline function Zygote.broadcast_forward(f, arg1::AbstractArray{<:Atom})
    out = dual_function_atom(f).(arg1)
    y = map(x -> value.(x), out)
    function bc_fwd_back(ȳ_in)
        ȳ = arg1 isa CuArray ? cu(ȳ_in) : ȳ_in
        barg1 = broadcast(ȳ, out) do y1, o1
            ps = partials(o1)
            Atom(0, y1 * ps[1], y1 * ps[2], y1 * ps[3], y1 * ps[4])
        end
        darg1 = unbroadcast(arg1, barg1)
        (nothing, nothing, darg1)
    end
    return y, bc_fwd_back
end

function combine_dual_GeneralInteraction(y1::SVector{3, T}, o1::SVector{3, Dual{Nothing, T, P}}, i::Integer) where {T, P}
    (
        LennardJones{false, Nothing, T, typeof(NoUnits), typeof(NoUnits)}(
            nothing, false, false,
            y1[1] * partials(o1[1], i) + y1[2] * partials(o1[2], i) + y1[3] * partials(o1[3], i),
            NoUnits, NoUnits,
        ),
        CoulombReactionField{T, T, T, T, typeof(NoUnits), typeof(NoUnits)}(
            y1[1] * partials(o1[1], i    ) + y1[2] * partials(o1[2], i    ) + y1[3] * partials(o1[3], i    ),
            y1[1] * partials(o1[1], i + 1) + y1[2] * partials(o1[2], i + 1) + y1[3] * partials(o1[3], i + 1),
            false,
            y1[1] * partials(o1[1], i + 2) + y1[2] * partials(o1[2], i + 2) + y1[3] * partials(o1[3], i + 2),
            y1[1] * partials(o1[1], i + 3) + y1[2] * partials(o1[2], i + 3) + y1[3] * partials(o1[3], i + 3),
            NoUnits, NoUnits,
        ),
    )
end

function combine_dual_SpecificInteraction(inter::HarmonicBond, y1, o1, i::Integer)
    (y1.f1[1] * partials(o1.f1[1], i    ) + y1.f1[2] * partials(o1.f1[2], i    ) + y1.f1[3] * partials(o1.f1[3], i    ) + y1.f2[1] * partials(o1.f2[1], i    ) + y1.f2[2] * partials(o1.f2[2], i    ) + y1.f2[3] * partials(o1.f2[3], i    ),
     y1.f1[1] * partials(o1.f1[1], i + 1) + y1.f1[2] * partials(o1.f1[2], i + 1) + y1.f1[3] * partials(o1.f1[3], i + 1) + y1.f2[1] * partials(o1.f2[1], i + 1) + y1.f2[2] * partials(o1.f2[2], i + 1) + y1.f2[3] * partials(o1.f2[3], i + 1))
end

function combine_dual_SpecificInteraction(inter::HarmonicAngle, y1, o1, i::Integer)
    (y1.f1[1] * partials(o1.f1[1], i    ) + y1.f1[2] * partials(o1.f1[2], i    ) + y1.f1[3] * partials(o1.f1[3], i    ) + y1.f2[1] * partials(o1.f2[1], i    ) + y1.f2[2] * partials(o1.f2[2], i    ) + y1.f2[3] * partials(o1.f2[3], i    ) + y1.f3[1] * partials(o1.f3[1], i    ) + y1.f3[2] * partials(o1.f3[2], i    ) + y1.f3[3] * partials(o1.f3[3], i    ),
     y1.f1[1] * partials(o1.f1[1], i + 1) + y1.f1[2] * partials(o1.f1[2], i + 1) + y1.f1[3] * partials(o1.f1[3], i + 1) + y1.f2[1] * partials(o1.f2[1], i + 1) + y1.f2[2] * partials(o1.f2[2], i + 1) + y1.f2[3] * partials(o1.f2[3], i + 1) + y1.f3[1] * partials(o1.f3[1], i + 1) + y1.f3[2] * partials(o1.f3[2], i + 1) + y1.f3[3] * partials(o1.f3[3], i + 1))
end

function combine_dual_SpecificInteraction(inter::PeriodicTorsion{6}, y1, o1, i::Integer)
    (
        (0, 0, 0, 0, 0, 0),
        (
            y1.f1[1] * partials(o1.f1[1], i     ) + y1.f1[2] * partials(o1.f1[2], i     ) + y1.f1[3] * partials(o1.f1[3], i     ) + y1.f2[1] * partials(o1.f2[1], i     ) + y1.f2[2] * partials(o1.f2[2], i     ) + y1.f2[3] * partials(o1.f2[3], i     ) + y1.f3[1] * partials(o1.f3[1], i     ) + y1.f3[2] * partials(o1.f3[2], i     ) + y1.f3[3] * partials(o1.f3[3], i     ) + y1.f4[1] * partials(o1.f4[1], i     ) + y1.f4[2] * partials(o1.f4[2], i     ) + y1.f4[3] * partials(o1.f4[3], i     ),
            y1.f1[1] * partials(o1.f1[1], i +  1) + y1.f1[2] * partials(o1.f1[2], i +  1) + y1.f1[3] * partials(o1.f1[3], i +  1) + y1.f2[1] * partials(o1.f2[1], i +  1) + y1.f2[2] * partials(o1.f2[2], i +  1) + y1.f2[3] * partials(o1.f2[3], i +  1) + y1.f3[1] * partials(o1.f3[1], i +  1) + y1.f3[2] * partials(o1.f3[2], i +  1) + y1.f3[3] * partials(o1.f3[3], i +  1) + y1.f4[1] * partials(o1.f4[1], i +  1) + y1.f4[2] * partials(o1.f4[2], i +  1) + y1.f4[3] * partials(o1.f4[3], i +  1),
            y1.f1[1] * partials(o1.f1[1], i +  2) + y1.f1[2] * partials(o1.f1[2], i +  2) + y1.f1[3] * partials(o1.f1[3], i +  2) + y1.f2[1] * partials(o1.f2[1], i +  2) + y1.f2[2] * partials(o1.f2[2], i +  2) + y1.f2[3] * partials(o1.f2[3], i +  2) + y1.f3[1] * partials(o1.f3[1], i +  2) + y1.f3[2] * partials(o1.f3[2], i +  2) + y1.f3[3] * partials(o1.f3[3], i +  2) + y1.f4[1] * partials(o1.f4[1], i +  2) + y1.f4[2] * partials(o1.f4[2], i +  2) + y1.f4[3] * partials(o1.f4[3], i +  2),
            y1.f1[1] * partials(o1.f1[1], i +  3) + y1.f1[2] * partials(o1.f1[2], i +  3) + y1.f1[3] * partials(o1.f1[3], i +  3) + y1.f2[1] * partials(o1.f2[1], i +  3) + y1.f2[2] * partials(o1.f2[2], i +  3) + y1.f2[3] * partials(o1.f2[3], i +  3) + y1.f3[1] * partials(o1.f3[1], i +  3) + y1.f3[2] * partials(o1.f3[2], i +  3) + y1.f3[3] * partials(o1.f3[3], i +  3) + y1.f4[1] * partials(o1.f4[1], i +  3) + y1.f4[2] * partials(o1.f4[2], i +  3) + y1.f4[3] * partials(o1.f4[3], i +  3),
            y1.f1[1] * partials(o1.f1[1], i +  4) + y1.f1[2] * partials(o1.f1[2], i +  4) + y1.f1[3] * partials(o1.f1[3], i +  4) + y1.f2[1] * partials(o1.f2[1], i +  4) + y1.f2[2] * partials(o1.f2[2], i +  4) + y1.f2[3] * partials(o1.f2[3], i +  4) + y1.f3[1] * partials(o1.f3[1], i +  4) + y1.f3[2] * partials(o1.f3[2], i +  4) + y1.f3[3] * partials(o1.f3[3], i +  4) + y1.f4[1] * partials(o1.f4[1], i +  4) + y1.f4[2] * partials(o1.f4[2], i +  4) + y1.f4[3] * partials(o1.f4[3], i +  4),
            y1.f1[1] * partials(o1.f1[1], i +  5) + y1.f1[2] * partials(o1.f1[2], i +  5) + y1.f1[3] * partials(o1.f1[3], i +  5) + y1.f2[1] * partials(o1.f2[1], i +  5) + y1.f2[2] * partials(o1.f2[2], i +  5) + y1.f2[3] * partials(o1.f2[3], i +  5) + y1.f3[1] * partials(o1.f3[1], i +  5) + y1.f3[2] * partials(o1.f3[2], i +  5) + y1.f3[3] * partials(o1.f3[3], i +  5) + y1.f4[1] * partials(o1.f4[1], i +  5) + y1.f4[2] * partials(o1.f4[2], i +  5) + y1.f4[3] * partials(o1.f4[3], i +  5),
        ),
        (
            y1.f1[1] * partials(o1.f1[1], i +  6) + y1.f1[2] * partials(o1.f1[2], i +  6) + y1.f1[3] * partials(o1.f1[3], i +  6) + y1.f2[1] * partials(o1.f2[1], i +  6) + y1.f2[2] * partials(o1.f2[2], i +  6) + y1.f2[3] * partials(o1.f2[3], i +  6) + y1.f3[1] * partials(o1.f3[1], i +  6) + y1.f3[2] * partials(o1.f3[2], i +  6) + y1.f3[3] * partials(o1.f3[3], i +  6) + y1.f4[1] * partials(o1.f4[1], i +  6) + y1.f4[2] * partials(o1.f4[2], i +  6) + y1.f4[3] * partials(o1.f4[3], i +  6),
            y1.f1[1] * partials(o1.f1[1], i +  7) + y1.f1[2] * partials(o1.f1[2], i +  7) + y1.f1[3] * partials(o1.f1[3], i +  7) + y1.f2[1] * partials(o1.f2[1], i +  7) + y1.f2[2] * partials(o1.f2[2], i +  7) + y1.f2[3] * partials(o1.f2[3], i +  7) + y1.f3[1] * partials(o1.f3[1], i +  7) + y1.f3[2] * partials(o1.f3[2], i +  7) + y1.f3[3] * partials(o1.f3[3], i +  7) + y1.f4[1] * partials(o1.f4[1], i +  7) + y1.f4[2] * partials(o1.f4[2], i +  7) + y1.f4[3] * partials(o1.f4[3], i +  7),
            y1.f1[1] * partials(o1.f1[1], i +  8) + y1.f1[2] * partials(o1.f1[2], i +  8) + y1.f1[3] * partials(o1.f1[3], i +  8) + y1.f2[1] * partials(o1.f2[1], i +  8) + y1.f2[2] * partials(o1.f2[2], i +  8) + y1.f2[3] * partials(o1.f2[3], i +  8) + y1.f3[1] * partials(o1.f3[1], i +  8) + y1.f3[2] * partials(o1.f3[2], i +  8) + y1.f3[3] * partials(o1.f3[3], i +  8) + y1.f4[1] * partials(o1.f4[1], i +  8) + y1.f4[2] * partials(o1.f4[2], i +  8) + y1.f4[3] * partials(o1.f4[3], i +  8),
            y1.f1[1] * partials(o1.f1[1], i +  9) + y1.f1[2] * partials(o1.f1[2], i +  9) + y1.f1[3] * partials(o1.f1[3], i +  9) + y1.f2[1] * partials(o1.f2[1], i +  9) + y1.f2[2] * partials(o1.f2[2], i +  9) + y1.f2[3] * partials(o1.f2[3], i +  9) + y1.f3[1] * partials(o1.f3[1], i +  9) + y1.f3[2] * partials(o1.f3[2], i +  9) + y1.f3[3] * partials(o1.f3[3], i +  9) + y1.f4[1] * partials(o1.f4[1], i +  9) + y1.f4[2] * partials(o1.f4[2], i +  9) + y1.f4[3] * partials(o1.f4[3], i +  9),
            y1.f1[1] * partials(o1.f1[1], i + 10) + y1.f1[2] * partials(o1.f1[2], i + 10) + y1.f1[3] * partials(o1.f1[3], i + 10) + y1.f2[1] * partials(o1.f2[1], i + 10) + y1.f2[2] * partials(o1.f2[2], i + 10) + y1.f2[3] * partials(o1.f2[3], i + 10) + y1.f3[1] * partials(o1.f3[1], i + 10) + y1.f3[2] * partials(o1.f3[2], i + 10) + y1.f3[3] * partials(o1.f3[3], i + 10) + y1.f4[1] * partials(o1.f4[1], i + 10) + y1.f4[2] * partials(o1.f4[2], i + 10) + y1.f4[3] * partials(o1.f4[3], i + 10),
            y1.f1[1] * partials(o1.f1[1], i + 11) + y1.f1[2] * partials(o1.f1[2], i + 11) + y1.f1[3] * partials(o1.f1[3], i + 11) + y1.f2[1] * partials(o1.f2[1], i + 11) + y1.f2[2] * partials(o1.f2[2], i + 11) + y1.f2[3] * partials(o1.f2[3], i + 11) + y1.f3[1] * partials(o1.f3[1], i + 11) + y1.f3[2] * partials(o1.f3[2], i + 11) + y1.f3[3] * partials(o1.f3[3], i + 11) + y1.f4[1] * partials(o1.f4[1], i + 11) + y1.f4[2] * partials(o1.f4[2], i + 11) + y1.f4[3] * partials(o1.f4[3], i + 11),
        ),
    )
end

function combine_dual_Atom(y1::SVector{3, T}, o1::SVector{3, Dual{Nothing, T, P}}, i::Integer, j::Integer, k::Integer, l::Integer) where {T, P}
    ps1, ps2, ps3 = partials(o1[1]), partials(o1[2]), partials(o1[3])
    Atom(
        0,
        y1[1] * ps1[i] + y1[2] * ps2[i] + y1[3] * ps3[i],
        y1[1] * ps1[j] + y1[2] * ps2[j] + y1[3] * ps3[j],
        y1[1] * ps1[k] + y1[2] * ps2[k] + y1[3] * ps3[k],
        y1[1] * ps1[l] + y1[2] * ps2[l] + y1[3] * ps3[l],
    )
end

@inline function Zygote.broadcast_forward(f,
                                            arg1,
                                            arg2::AbstractArray{SVector{D, T}},
                                            arg3::AbstractArray{SVector{D, T}},
                                            arg4::AbstractArray{<:Atom},
                                            arg5::AbstractArray{<:Atom},
                                            arg6::Tuple{SVector{D, T}},
                                            arg7::Base.RefValue{<:Unitful.FreeUnits},
                                            arg8) where {D, T}
    out = dual_function_force_broadcast(f).(arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
    y = map(x -> value.(x), out)
    function bc_fwd_back(ȳ_in)
        ȳ = arg2 isa CuArray ? cu(ȳ_in) : ȳ_in
        darg1 = unbroadcast(arg1, broadcast(combine_dual_GeneralInteraction, ȳ, out, 1))
        darg2 = unbroadcast(arg2, broadcast((y1, o1) -> SVector{D, T}(sum_partials(o1, y1,  5),
                                sum_partials(o1, y1,  6), sum_partials(o1, y1,  7)), ȳ, out))
        darg3 = unbroadcast(arg3, broadcast((y1, o1) -> SVector{D, T}(sum_partials(o1, y1,  8),
                                sum_partials(o1, y1,  9), sum_partials(o1, y1, 10)), ȳ, out))
        darg4 = unbroadcast(arg4, broadcast(combine_dual_Atom, ȳ, out, 11, 12, 13, 14))
        darg5 = unbroadcast(arg5, broadcast(combine_dual_Atom, ȳ, out, 15, 16, 17, 18))
        darg6 = unbroadcast(arg6, broadcast((y1, o1) -> SVector{D, T}(sum_partials(o1, y1, 19),
                                sum_partials(o1, y1, 20), sum_partials(o1, y1, 21)), ȳ, out))
        darg7 = nothing
        darg8 = nothing
        return (nothing, nothing, darg1, darg2, darg3, darg4, darg5, darg6, darg7, darg8)
    end
    return y, bc_fwd_back
end

@inline function Zygote.broadcast_forward(f,
                                            arg1::AbstractArray{<:SpecificInteraction},
                                            arg2::AbstractArray{SVector{D, T}},
                                            arg3::AbstractArray{SVector{D, T}},
                                            arg4::Tuple{SVector{D, T}}) where {D, T}
    out = dual_function_specific_2_atoms(f).(arg1, arg2, arg3, arg4)
    y = broadcast(o1 -> SpecificForce2Atoms{D, T}(value.(o1.f1), value.(o1.f2)), out)
    function bc_fwd_back(ȳ_in)
        ȳ = arg1 isa CuArray ? cu(ȳ_in) : ȳ_in
        darg1 = unbroadcast(arg1, broadcast(combine_dual_SpecificInteraction, arg1, ȳ, out, 1))
        darg2 = unbroadcast(arg2, broadcast((y1, o1) -> SVector{D, T}(
                    sum_partials(o1.f1, y1.f1, 3) + sum_partials(o1.f2, y1.f2, 3),
                    sum_partials(o1.f1, y1.f1, 4) + sum_partials(o1.f2, y1.f2, 4),
                    sum_partials(o1.f1, y1.f1, 5) + sum_partials(o1.f2, y1.f2, 5)),
                    ȳ, out))
        darg3 = unbroadcast(arg3, broadcast((y1, o1) -> SVector{D, T}(
                    sum_partials(o1.f1, y1.f1, 6) + sum_partials(o1.f2, y1.f2, 6),
                    sum_partials(o1.f1, y1.f1, 7) + sum_partials(o1.f2, y1.f2, 7),
                    sum_partials(o1.f1, y1.f1, 8) + sum_partials(o1.f2, y1.f2, 8)),
                    ȳ, out))
        darg4 = unbroadcast(arg4, broadcast((y1, o1) -> SVector{D, T}(
                    sum_partials(o1.f1, y1.f1,  9) + sum_partials(o1.f2, y1.f2,  9),
                    sum_partials(o1.f1, y1.f1, 10) + sum_partials(o1.f2, y1.f2, 10),
                    sum_partials(o1.f1, y1.f1, 11) + sum_partials(o1.f2, y1.f2, 11)),
                    ȳ, out))
        return (nothing, nothing, darg1, darg2, darg3, darg4)
    end
    return y, bc_fwd_back
end

@inline function Zygote.broadcast_forward(f,
                                            arg1::AbstractArray{<:SpecificInteraction},
                                            arg2::AbstractArray{SVector{D, T}},
                                            arg3::AbstractArray{SVector{D, T}},
                                            arg4::AbstractArray{SVector{D, T}},
                                            arg5::Tuple{SVector{D, T}}) where {D, T}
    out = dual_function_specific_3_atoms(f).(arg1, arg2, arg3, arg4, arg5)
    y = broadcast(o1 -> SpecificForce3Atoms{D, T}(value.(o1.f1), value.(o1.f2), value.(o1.f3)), out)
    function bc_fwd_back(ȳ_in)
        ȳ = arg1 isa CuArray ? cu(ȳ_in) : ȳ_in
        darg1 = unbroadcast(arg1, broadcast(combine_dual_SpecificInteraction, arg1, ȳ, out, 1))
        darg2 = unbroadcast(arg2, broadcast((y1, o1) -> SVector{D, T}(
                    sum_partials(o1.f1, y1.f1, 3) + sum_partials(o1.f2, y1.f2, 3) + sum_partials(o1.f3, y1.f3, 3),
                    sum_partials(o1.f1, y1.f1, 4) + sum_partials(o1.f2, y1.f2, 4) + sum_partials(o1.f3, y1.f3, 4),
                    sum_partials(o1.f1, y1.f1, 5) + sum_partials(o1.f2, y1.f2, 5) + sum_partials(o1.f3, y1.f3, 5)),
                    ȳ, out))
        darg3 = unbroadcast(arg3, broadcast((y1, o1) -> SVector{D, T}(
                    sum_partials(o1.f1, y1.f1, 6) + sum_partials(o1.f2, y1.f2, 6) + sum_partials(o1.f3, y1.f3, 6),
                    sum_partials(o1.f1, y1.f1, 7) + sum_partials(o1.f2, y1.f2, 7) + sum_partials(o1.f3, y1.f3, 7),
                    sum_partials(o1.f1, y1.f1, 8) + sum_partials(o1.f2, y1.f2, 8) + sum_partials(o1.f3, y1.f3, 8)),
                    ȳ, out))
        darg4 = unbroadcast(arg4, broadcast((y1, o1) -> SVector{D, T}(
                    sum_partials(o1.f1, y1.f1,  9) + sum_partials(o1.f2, y1.f2,  9) + sum_partials(o1.f3, y1.f3,  9),
                    sum_partials(o1.f1, y1.f1, 10) + sum_partials(o1.f2, y1.f2, 10) + sum_partials(o1.f3, y1.f3, 10),
                    sum_partials(o1.f1, y1.f1, 11) + sum_partials(o1.f2, y1.f2, 11) + sum_partials(o1.f3, y1.f3, 11)),
                    ȳ, out))
        darg5 = unbroadcast(arg5, broadcast((y1, o1) -> SVector{D, T}(
                    sum_partials(o1.f1, y1.f1, 12) + sum_partials(o1.f2, y1.f2, 12) + sum_partials(o1.f3, y1.f3, 12),
                    sum_partials(o1.f1, y1.f1, 13) + sum_partials(o1.f2, y1.f2, 13) + sum_partials(o1.f3, y1.f3, 13),
                    sum_partials(o1.f1, y1.f1, 14) + sum_partials(o1.f2, y1.f2, 14) + sum_partials(o1.f3, y1.f3, 14)),
                    ȳ, out))
        return (nothing, nothing, darg1, darg2, darg3, darg4, darg5)
    end
    return y, bc_fwd_back
end

@inline function Zygote.broadcast_forward(f,
                                            arg1::AbstractArray{<:SpecificInteraction},
                                            arg2::AbstractArray{SVector{D, T}},
                                            arg3::AbstractArray{SVector{D, T}},
                                            arg4::AbstractArray{SVector{D, T}},
                                            arg5::AbstractArray{SVector{D, T}},
                                            arg6::Tuple{SVector{D, T}}) where {D, T}
    out = dual_function_specific_4_atoms(f).(arg1, arg2, arg3, arg4, arg5, arg6)
    y = broadcast(o1 -> SpecificForce4Atoms{D, T}(value.(o1.f1), value.(o1.f2), value.(o1.f3), value.(o1.f4)), out)
    function bc_fwd_back(ȳ_in)
        ȳ = arg1 isa CuArray ? cu(ȳ_in) : ȳ_in
        darg1 = unbroadcast(arg1, broadcast(combine_dual_SpecificInteraction, arg1, ȳ, out, 1))
        darg2 = unbroadcast(arg2, broadcast((y1, o1) -> SVector{D, T}(
                    sum_partials(o1.f1, y1.f1, 13) + sum_partials(o1.f2, y1.f2, 13) + sum_partials(o1.f3, y1.f3, 13) + sum_partials(o1.f4, y1.f4, 13),
                    sum_partials(o1.f1, y1.f1, 14) + sum_partials(o1.f2, y1.f2, 14) + sum_partials(o1.f3, y1.f3, 14) + sum_partials(o1.f4, y1.f4, 14),
                    sum_partials(o1.f1, y1.f1, 15) + sum_partials(o1.f2, y1.f2, 15) + sum_partials(o1.f3, y1.f3, 15) + sum_partials(o1.f4, y1.f4, 15)),
                    ȳ, out))
        darg3 = unbroadcast(arg3, broadcast((y1, o1) -> SVector{D, T}(
                    sum_partials(o1.f1, y1.f1, 16) + sum_partials(o1.f2, y1.f2, 16) + sum_partials(o1.f3, y1.f3, 16) + sum_partials(o1.f4, y1.f4, 16),
                    sum_partials(o1.f1, y1.f1, 17) + sum_partials(o1.f2, y1.f2, 17) + sum_partials(o1.f3, y1.f3, 17) + sum_partials(o1.f4, y1.f4, 17),
                    sum_partials(o1.f1, y1.f1, 18) + sum_partials(o1.f2, y1.f2, 18) + sum_partials(o1.f3, y1.f3, 18) + sum_partials(o1.f4, y1.f4, 18)),
                    ȳ, out))
        darg4 = unbroadcast(arg4, broadcast((y1, o1) -> SVector{D, T}(
                    sum_partials(o1.f1, y1.f1, 19) + sum_partials(o1.f2, y1.f2, 19) + sum_partials(o1.f3, y1.f3, 19) + sum_partials(o1.f4, y1.f4, 19),
                    sum_partials(o1.f1, y1.f1, 20) + sum_partials(o1.f2, y1.f2, 20) + sum_partials(o1.f3, y1.f3, 20) + sum_partials(o1.f4, y1.f4, 20),
                    sum_partials(o1.f1, y1.f1, 21) + sum_partials(o1.f2, y1.f2, 21) + sum_partials(o1.f3, y1.f3, 21) + sum_partials(o1.f4, y1.f4, 21)),
                    ȳ, out))
        darg5 = unbroadcast(arg5, broadcast((y1, o1) -> SVector{D, T}(
                    sum_partials(o1.f1, y1.f1, 22) + sum_partials(o1.f2, y1.f2, 22) + sum_partials(o1.f3, y1.f3, 22) + sum_partials(o1.f4, y1.f4, 22),
                    sum_partials(o1.f1, y1.f1, 23) + sum_partials(o1.f2, y1.f2, 23) + sum_partials(o1.f3, y1.f3, 23) + sum_partials(o1.f4, y1.f4, 23),
                    sum_partials(o1.f1, y1.f1, 24) + sum_partials(o1.f2, y1.f2, 24) + sum_partials(o1.f3, y1.f3, 24) + sum_partials(o1.f4, y1.f4, 24)),
                    ȳ, out))
        darg6 = unbroadcast(arg6, broadcast((y1, o1) -> SVector{D, T}(
                    sum_partials(o1.f1, y1.f1, 25) + sum_partials(o1.f2, y1.f2, 25) + sum_partials(o1.f3, y1.f3, 25) + sum_partials(o1.f4, y1.f4, 25),
                    sum_partials(o1.f1, y1.f1, 26) + sum_partials(o1.f2, y1.f2, 26) + sum_partials(o1.f3, y1.f3, 26) + sum_partials(o1.f4, y1.f4, 26),
                    sum_partials(o1.f1, y1.f1, 27) + sum_partials(o1.f2, y1.f2, 27) + sum_partials(o1.f3, y1.f3, 27) + sum_partials(o1.f4, y1.f4, 27)),
                    ȳ, out))
        return (nothing, nothing, darg1, darg2, darg3, darg4, darg5, darg6)
    end
    return y, bc_fwd_back
end

@inline function Zygote.broadcast_forward(f::typeof(getf1),
                                            arg1::AbstractArray{<:SpecificForce2Atoms}) where {D, T}
    return f.(arg1), ȳ -> (nothing, nothing, unbroadcast(arg1, broadcast(y1 -> (f1=y1, f2=zero(y1)), ȳ)))
end

@inline function Zygote.broadcast_forward(f::typeof(getf1),
                                            arg1::AbstractArray{<:SpecificForce3Atoms}) where {D, T}
    return f.(arg1), ȳ -> (nothing, nothing, unbroadcast(arg1, broadcast(y1 -> (f1=y1, f2=zero(y1), f3=zero(y1)), ȳ)))
end

@inline function Zygote.broadcast_forward(f::typeof(getf1),
                                            arg1::AbstractArray{<:SpecificForce4Atoms}) where {D, T}
    return f.(arg1), ȳ -> (nothing, nothing, unbroadcast(arg1, broadcast(y1 -> (f1=y1, f2=zero(y1), f3=zero(y1), f4=zero(y1)), ȳ)))
end

@inline function Zygote.broadcast_forward(f::typeof(getf2),
                                            arg1::AbstractArray{<:SpecificForce2Atoms}) where {D, T}
    return f.(arg1), ȳ -> (nothing, nothing, unbroadcast(arg1, broadcast(y1 -> (f1=zero(y1), f2=y1), ȳ)))
end

@inline function Zygote.broadcast_forward(f::typeof(getf2),
                                            arg1::AbstractArray{<:SpecificForce3Atoms}) where {D, T}
    return f.(arg1), ȳ -> (nothing, nothing, unbroadcast(arg1, broadcast(y1 -> (f1=zero(y1), f2=y1, f3=zero(y1)), ȳ)))
end

@inline function Zygote.broadcast_forward(f::typeof(getf2),
                                            arg1::AbstractArray{<:SpecificForce4Atoms}) where {D, T}
    return f.(arg1), ȳ -> (nothing, nothing, unbroadcast(arg1, broadcast(y1 -> (f1=zero(y1), f2=y1, f3=zero(y1), f4=zero(y1)), ȳ)))
end

@inline function Zygote.broadcast_forward(f::typeof(getf3),
                                            arg1::AbstractArray{<:SpecificForce3Atoms}) where {D, T}
    return f.(arg1), ȳ -> (nothing, nothing, unbroadcast(arg1, broadcast(y1 -> (f1=zero(y1), f2=zero(y1), f3=y1), ȳ)))
end

@inline function Zygote.broadcast_forward(f::typeof(getf3),
                                            arg1::AbstractArray{<:SpecificForce4Atoms}) where {D, T}
    return f.(arg1), ȳ -> (nothing, nothing, unbroadcast(arg1, broadcast(y1 -> (f1=zero(y1), f2=zero(y1), f3=y1, f4=zero(y1)), ȳ)))
end

@inline function Zygote.broadcast_forward(f::typeof(getf4),
                                            arg1::AbstractArray{<:SpecificForce4Atoms}) where {D, T}
    return f.(arg1), ȳ -> (nothing, nothing, unbroadcast(arg1, broadcast(y1 -> (f1=zero(y1), f2=zero(y1), f3=zero(y1), f4=y1), ȳ)))
end

# Use fast broadcast path on CPU
for op in (:+, :-, :*, :/, :force, :force_nounit, :mass, :remove_molar, :ustrip,
            :ustrip_vec, :wrap_coords_vec, :getf1, :getf2, :getf3, :getf4)
    @eval Zygote.@adjoint Broadcast.broadcasted(::Broadcast.AbstractArrayStyle, f::typeof($op), args...) = Zygote.broadcast_forward(f, args...)
    # Avoid ambiguous dispatch
    @eval Zygote.@adjoint Broadcast.broadcasted(::CUDA.AbstractGPUArrayStyle  , f::typeof($op), args...) = Zygote.broadcast_forward(f, args...)
end


text:
module TestLogMgrs
  using Base.Test
  using HelloDB.LogMgrs: LogMgr, append, flush, iterator
  using TestHelpers: set_test_db, drop_test_db
  using HelloDB: DB
  using HelloDB.Data: INT_SIZE, BLOCK_SIZE
  using HelloDB.LogIterators: LogIterator

  function test_init()
    # with no previous logfile
    set_test_db()
    filemgr = LogMgr()
    logfile = DB.filemgr.dbdir * "/log"
    @test isfile(logfile)
    @test filemgr.currentpos == INT_SIZE

    # with previous existing logfile
    file = open(logfile, "w")
    write(file, "foo bar baz")
    close(file)
    logmgr = LogMgr()
    @test logmgr.currentpos > INT_SIZE
    drop_test_db()
  end

  function test_append()
    set_test_db()
    logmgr = LogMgr()
    lsn = append(logmgr, ["FOO", "BAR"])
    @test lsn == 0
    lsn = append(logmgr, [ASCIIString(fill('a', BLOCK_SIZE))])
    @test lsn == 1
    drop_test_db()
  end

  function test_flush()
    set_test_db()
    logmgr = LogMgr()
    append(logmgr, ["FOO", "BAR"])
    flush(logmgr, 1)
    logmgr = LogMgr()
    @test logmgr.currentpos > INT_SIZE
    drop_test_db()
  end

  function test_iterator()
    set_test_db()
    logmgr = LogMgr()
    append(logmgr, ["FOO", "BAR"])
    it = iterator(logmgr)
    @test typeof(it) ==  LogIterator

    # test that prev block flushed
    logmgr = LogMgr()
    @test logmgr.currentpos > INT_SIZE
    drop_test_db()
  end

  test_init()
  test_append()
  test_flush()
  test_iterator()
end


text:
export set_ωpml!, set_boundft!, set_Npml!, set_kbloch!  # basic setter functions
export create_e⁻ⁱᵏᴸ
export clear_objs!  # plural because clearing both electric and magnetic quantities
export complete_fields

# Do not export Model; quality it with the package name MaxwellFDFD, because I would have
# similar types in other packages such as MaxwellSALT and MaxwellGuide.
Base.@kwdef mutable struct Model{K,Kₑ,Kₘ,
                                 Kₑₜ,Kₘₜ,Kₑₗ,Kₘₗ,
                                 K²,
                                 K₊₂,AK₊₂<:AbsArrComplexF{K₊₂}}
    # Frequency
    ωpml::Number = 0.0  # can be complex

    # Field type of equation
    ft_eq::FieldType

    # Grid
    grid::Grid{K}
    cmpₛ::SInt{K}  # Cartesian components of shape dimension
    cmpₑₜ::SInt{Kₑₜ}  # Cartesian components of transverse E-field
    cmpₘₜ::SInt{Kₘₜ}  # Cartesian components of transverse H-field
    cmpₑₗ::SInt{Kₑₗ}  # Cartesian components of longitudinal E-field
    cmpₘₗ::SInt{Kₘₗ}  # Cartesian components of longitudinal H-field

    # Boundary properties (in transverse directions)
    boundft::SVec{K,FieldType} = SVec(ntuple(k->EE, Val(K)))
    Npml::Tuple2{SInt{K}} = (SVec(ntuple(k->0, Val(K))), SVec(ntuple(k->0, Val(K))))
    kbloch::SFloat{K} = SVec(ntuple(k->0.0, Val(K)))  # [kx_bloch, ky_bloch, kz_bloch]

    # Material parameter arrays
    εₜarr::AK₊₂ = create_param_array(grid.N, ncmp=Kₑₜ)  # filled with zeros
    μₜarr::AK₊₂ = create_param_array(grid.N, ncmp=Kₘₜ)  # filled with zeros

    εₗarr::AK₊₂ = create_param_array(grid.N, ncmp=Kₑ-Kₑₜ)  # filled with zeros
    μₗarr::AK₊₂ = create_param_array(grid.N, ncmp=Kₘ-Kₘₜ)  # filled with zeros

    # Temporary storages for assignment and smoothing of material parameters
    oind2shp::Vector{Shape{K,K²}} = Shape{K,K²}[]
    oind2εind::Vector{ParamInd} = ParamInd[]
    oind2μind::Vector{ParamInd} = ParamInd[]
    εind2ε::Vector{S²ComplexF3} = S²ComplexF3[]
    μind2μ::Vector{S²ComplexF3} = S²ComplexF3[]

    # Boolean flags indicating if the transverse E- and H-field dimensions are orthogonal to
    # the shape dimensions.  Used for material parameter smoothing in calc_matparams!().
    iseₜ˔shp::Bool
    ishₜ˔shp::Bool

    # Indexing scheme for DOFs
    order_cmpfirst::Bool = true
end

function Base.size(mdl::Model{K,Kₑ,Kₘ,Kₑₜ,Kₘₜ}, ft::FieldType) where {K,Kₑ,Kₘ,Kₑₜ,Kₘₜ}
    Kf = ft==EE ? Kₑₜ : Kₘₜ
    sz_grid = size(mdl.grid)
    sz = mdl.order_cmpfirst ? (Kf, sz_grid...) : (sz_grid..., Kf)

    return sz
end

Base.length(mdl::Model, ft::FieldType) = prod(size(mdl, ft))

# Basic setters
set_ωpml!(mdl::Model, ωpml::Number) = (mdl.ωpml = ωpml; nothing)
set_boundft!(mdl::Model{K}, boundft::AbsVec{FieldType}) where {K} = (mdl.boundft = SVec{K}(boundft); nothing)
set_Npml!(mdl::Model{K}, Npml::Tuple2{AbsVecInteger}) where {K} = (mdl.Npml = SVec{K}.(Npml); nothing)
set_kbloch!(mdl::Model{K}, kbloch::AbsVecReal) where {K} = (mdl.kbloch = SVec{K}(kbloch); nothing)

create_e⁻ⁱᵏᴸ(mdl::Model) = exp.(-im .* mdl.kbloch .* mdl.grid.L)

# Main functions
function clear_objs!(mdl::Model)
    mdl.εₜarr .= 0
    mdl.μₜarr .= 0

    mdl.εₗarr .= 0
    mdl.μₗarr .= 0

    empty!(mdl.oind2shp)
    empty!(mdl.oind2εind)
    empty!(mdl.oind2μind)
    empty!(mdl.εind2ε)
    empty!(mdl.μind2μ)

    return nothing
end

MaxwellBase.add_obj!(mdl::Model{K}, matname::String, shapes::AbsVec{<:Shape{K}};
                     ε::MatParam=1.0, μ::MatParam=1.0) where {K} =
    add_obj!(mdl, matname, ε=ε, μ=μ, tuple(shapes...))

function MaxwellBase.add_obj!(mdl::Model{K}, matname::String, shapes::Shape{K}...;
                              ε::MatParam=1.0, μ::MatParam=1.0) where {K}
    mat = Material(matname, ε=ε, μ=μ)
    isbidirectional(EE, mat, mdl) || @error "ε = $ε should not couple transverse and longitudinal fields."
    isbidirectional(HH, mat, mdl) || @error "μ = $μ should not couple transverse and longitudinal fields."

    for shp = shapes  # shapes is tuple
        obj = Object(shp, mat)
        add_obj!(mdl.oind2shp, (mdl.oind2εind,mdl.oind2μind), (mdl.εind2ε,mdl.μind2μ), obj)
    end

    return nothing
end

# Returns true if the material parameter does not couple the transverse and longitudinal
# components of the fields; false otherwise.
function isbidirectional(ft::FieldType, mat::Material, mdl::Model)
    nft = Int(ft)
    cmpₜ = (mdl.cmpₑₜ, mdl.cmpₘₜ)[nft]
    cmpₗ = (mdl.cmpₑₗ, mdl.cmpₘₗ)[nft]

    param = mat.param[nft]

    return iszero(param[cmpₜ,cmpₗ]) && iszero(param[cmpₗ,cmpₜ])
end

# Implement the following functions for each specific Model.
function calc_matparams! end
function complete_fields end

function field2vec(E::NTuple{Kₑ,AbsVecNumber{K}}, H::NTuple{Kₘ,AbsVecNumber{K}},
                   mdl::Model{K,Kₑ,Kₘ}
                   ) where {K,Kₑ,Kₘ}
    U = mdl.ft_eq==EE ? E : H
    cmpᵤₜ = mdl.ft_eq==EE ? mdl.cmpₑₜ : mdl.cmpₘₜ
    Uₜ = cat(U[cmpᵤₜ]..., dims=K+1)

    order_cmpfirst = mdl.order_cmpfirst
    if order_cmpfirst
        Uₜ = permutedims(Uₜ, circshift(1:K, 1))
    end

    return field_arr2vec(Uₜ; order_cmpfirst)
end

include("poynting.jl")
include("calc_mode.jl")

# Specific models
include("full.jl")
include("te.jl")
include("tm.jl")


text:
MeanScorePlot=zeros(1,Iterations);
MinScorePlot=zeros(1,Iterations);
MaxScorePlot=zeros(1,Iterations);
for g=1:Iterations
    MeanScorePlot(1,g)=mean(ScoreLog{g,1}{1,1});
    MinScorePlot(1,g)=min(ScoreLog{g,1}{1,1});
    MaxScorePlot(1,g)=max(ScoreLog{g,1}{1,1});
end
figure
subplot(1,3,1)
plot((1:Iterations),MeanScorePlot)
title('Mean Score')
xlabel('Generation')
ylabel('Score')
axis([1,Iterations 0,100])

subplot(1,3,2)
plot((1:Iterations),MinScorePlot)
title('Min Score')
xlabel('Generation')
ylabel('Score')
axis([1,Iterations 0,100])

subplot(1,3,3)
plot((1:Iterations),MaxScorePlot)
title('Max Score')
xlabel('Generation')
ylabel('Score')
axis([1,Iterations 0,100])

max_stars_repo_name: DesignEngrLab/ColorGraph

max_issues_repo_name: DesignEngrLab/ColorGraph

max_forks_repo_name: DesignEngrLab/ColorGraph


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    SubgraphView, DiSubgraphView,
    subgraphview, supergraph, nodesubgraph, edgesubgraph

struct SubgraphView{T<:UndirectedGraph} <: UndirectedGraph
    graph::T
    nodes::Set{Int}
    edges::Set{Int}
end
SubgraphView(graph::UndirectedGraph, nodes, edges) = SubgraphView(graph, Set{Int}(nodes), Set{Int}(edges))

Base.:(==)(g1::SubgraphView, g2::SubgraphView) = g1.graph == g2.graph && g1.nodes == g2.nodes && g1.edges == g2.edges

struct DiSubgraphView{T<:DirectedGraph} <: DirectedGraph
    graph::T
    nodes::Set{Int}
    edges::Set{Int}
end
DiSubgraphView(graph::UndirectedGraph, nodes, edges) = DiSubgraphView(graph, Set{Int}(nodes), Set{Int}(edges))

Base.:(==)(g1::DiSubgraphView, g2::DiSubgraphView) = g1.graph == g2.graph && g1.nodes == g2.nodes && g1.edges == g2.edges

"""
    subgraphview(graph::UndirectedGraph, nodes, edges) -> SubgraphView
    subgraphview(graph::DirectedGraph, nodes, edges) -> DiSubgraphView

Generate subgraph view with the given nodes and edges.
"""
subgraphview(graph::UndirectedGraph, nodes, edges) = SubgraphView(graph, nodes, edges)
subgraphview(graph::DirectedGraph, nodes, edges) = DiSubgraphView(graph, nodes, edges)

"""
    supergraph(view::SubgraphView) -> UndirectedGraph
    supergraph(view::DiSubgraphView) -> DirectedGraph

Returns the graph on which the view is based.
"""
supergraph(view::SubgraphView) = view.graph
supergraph(view::DiSubgraphView) = view.graph

function nodesubgraphedges(graph::AbstractGraph, nodes)
    edges = Set{Int}()
    for n in nodes
        for (ninc, nadj) in neighbors(graph, n)
            nadj in nodes && push!(edges, ninc)
        end
    end
    return edges
end

"""
    nodesubgraph(graph::UndirectedGraph, nodes) -> SubgraphView
    nodesubgraph(graph::DirectedGraph, nodes) -> DiSubgraphView

Generate node-induced subgraph view.
"""
function nodesubgraph(graph, nodes)
    edges = nodesubgraphedges(graph, nodes)
    return subgraphview(graph, Set(nodes), edges)
end

function edgesubgraphnodes(graph::AbstractGraph, edges)
    nodes = Set{Int}()
    for e in edges
        union!(nodes, getedge(graph, e))
    end
    return nodes
end

"""
    edgesubgraph(graph::UndirectedGraph, edges) -> SubgraphView
    edgesubgraph(graph::DirectedGraph, edges) -> DiSubgraphView

Generate edge-induced subgraph view.
"""
function edgesubgraph(graph, edges)
    nodes = edgesubgraphnodes(graph, edges)
    return subgraphview(graph, nodes, Set(edges))
end

function neighbors(view::Union{SubgraphView,DiSubgraphView}, idx)
    (idx in view.nodes) || throw(KeyError(idx))
    return Dict(
        n for n in neighbors(view.graph, idx)
        if n.first in view.edges && n.second in view.nodes
    )
end

function outneighbors(view::DiSubgraphView, idx)
    (idx in view.nodes) || throw(KeyError(idx))
    return Dict(
        n for n in outneighbors(view.graph, idx)
        if n.first in view.edges && n.second in view.nodes
    )
end

function inneighbors(view::DiSubgraphView, idx)
    (idx in view.nodes) || throw(KeyError(idx))
    return Dict(
        n for n in inneighbors(view.graph, idx)
        if n.first in view.edges && n.second in view.nodes
    )
end

function getedge(view::Union{SubgraphView,DiSubgraphView}, idx)
    idx in view.edges || throw(KeyError(idx))
    return getedge(view.graph, idx)
end

function findedgekey(view::Union{SubgraphView,DiSubgraphView}, u, v)
    u in view.nodes || return
    v in view.nodes || return
    return findedgekey(view.graph, u, v)
end

function nodeattr(view::Union{SubgraphView,DiSubgraphView}, idx)
    idx in view.nodes || throw(KeyError(idx))
    return nodeattr(view.graph, idx)
end

function edgeattr(view::Union{SubgraphView,DiSubgraphView}, idx)
    idx in view.edges || throw(KeyError(idx))
    return edgeattr(view.graph, idx)
end

nodeset(view::Union{SubgraphView,DiSubgraphView}) = copy(view.nodes)
edgeset(view::Union{SubgraphView,DiSubgraphView}) = copy(view.edges)

nodecount(view::Union{SubgraphView,DiSubgraphView}) = length(view.nodes)
edgecount(view::Union{SubgraphView,DiSubgraphView}) = length(view.edges)

nodeattrs(view::Union{SubgraphView,DiSubgraphView}) = nodeattrs(view.graph)
edgeattrs(view::Union{SubgraphView,DiSubgraphView}) = edgeattrs(view.graph)

nodeattrtype(view::Union{SubgraphView,DiSubgraphView}) = nodeattrtype(view.graph)
edgeattrtype(view::Union{SubgraphView,DiSubgraphView}) = edgeattrtype(view.graph)

hascachestore(view::Union{SubgraphView,DiSubgraphView}) = hascachestore(view.graph)
hascache(view::Union{SubgraphView,DiSubgraphView}, key) = hascache(view.graph, key)
getcache(view::Union{SubgraphView,DiSubgraphView}, key) = getcache(view.graph, key)
setcache!(view::Union{SubgraphView,DiSubgraphView}, key; kwargs...) = setcache!(view.graph, key; kwargs...)
clearcache!(view::Union{SubgraphView,DiSubgraphView}) = clearcache!(view.graph)


text:
mutable struct MyState
  w
  x
  y
  z
end

@testset "Gibbs" begin
  init = MyState(0.0, 0.0, 0.0, 0.0)

  function callback_fn(state, iter, pbar)
    out = nothing
    iter % 20 == 0 && begin
      ll = round(normlogpdf(state.x), digits=3)
      set_description(pbar, "loglike: $(ll)")
      out = Dict(:ll => ll)
    end
    return out
  end

  function update!(state)
    state.w += 1
    state.x = randn()
  end

  chain, laststate, summary_stats = CDE.MCMC.gibbs(init, update!, thin=10,
                                                   monitors=[[:x, :y], [:z]],
                                                   callback_fn=callback_fn,
                                                   nsamps=[500, 10], nburn=500)
  @test length(chain) == 2
  @test length(chain[1]) == 500
  @test length(chain[2]) == 10
  @test length(summary_stats) == 275
  println("Chain start: $(chain[1][1])")
  println("Chain end: $(chain[1][end])")
  println("Last state: $(laststate)")
  println("Head summary stats: $(summary_stats[1])")
  println("Last summary stats: $(summary_stats[end])")
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation