diff --git a/.gitignore b/.gitignore
index 6c6fb97..3ff322a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,4 +8,5 @@
 /docs/src/examples/*
 /docs/build/*
 
-.vscode
\ No newline at end of file
+.vscode
+LocalPreferences.toml
diff --git a/Project.toml b/Project.toml
index 92c67de..c9daf2d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,23 +1,8 @@
+authors = ["James Fairbanks", "Sophie Libkind"]
 name = "AlgebraicDynamics"
 uuid = "5fd6ff03-a254-427e-8840-ba658f502e32"
-authors = ["James Fairbanks", "Sophie Libkind"]
 version = "0.1.5"
 
-[deps]
-Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
-Compose = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
-DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
-LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
-RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
-RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-
 [compat]
 Catlab = "0.11, 0.12, 0.13"
 Compose = "^0.9.1"
@@ -33,7 +18,25 @@ RecursiveArrayTools = "^2.10.0"
 StaticArrays = "0.12, 1.0"
 julia = "1.0"
 
+[deps]
+Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
+Compose = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
+DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634"
+LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
 [extras]
+CPUSummary = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/src/AlgebraicDynamics.jl b/src/AlgebraicDynamics.jl
index fc5a028..a1a9c92 100644
--- a/src/AlgebraicDynamics.jl
+++ b/src/AlgebraicDynamics.jl
@@ -8,5 +8,7 @@ using Catlab.Programs
 include("uwd_dynam.jl")
 include("dwd_dynam.jl")
 include("cpg_dynam.jl")
+include("dwd_dils.jl")
+
 
 end # module
diff --git a/src/dwd_dils.jl b/src/dwd_dils.jl
new file mode 100644
index 0000000..8669890
--- /dev/null
+++ b/src/dwd_dils.jl
@@ -0,0 +1,284 @@
+module DWDDILS
+export Monomial, tensor, positions, directions, DILS, input_polys, input_poly,
+       output_polys, output_poly, nstates, forward_readout, backward_readout,
+       eval_dynamics, isa_state, sequential_compose, layerize,
+       generalized_kroenecker, oapply, unit
+
+
+using Test
+using LinearAlgebra
+
+using Catlab.WiringDiagrams.DirectedWiringDiagrams
+using Catlab.Graphs, Catlab.CategoricalAlgebra
+
+
+"""
+
+A monomial is a pair of integers  A => B representing the poly ``\\mathbb{R}^By^{\\mathbb{R}^A}``
+"""
+
+const Monomial = Pair{Int,Int}
+tensor(ms::Vector{Monomial}) = sum(first.(ms)) => sum(last.(ms))
+
+positions(m::Monomial) = m.second
+directions(m::Monomial) = m.first
+positions(ms::Vector{Monomial}) = sum(first.(ms))
+
+
+struct DILS{T}
+  input_polys::Vector{Monomial}
+  output_polys::Vector{Monomial}
+  nstates::Int # dimension of the state space
+  forward_readout::Function
+  backward_readout::Function
+  eval_dynamics::Function
+end
+
+DILS{T}(input_poly::Monomial, output_poly::Monomial, nstates::Int,
+        forward_readout::Function, backward_readout::Function,
+        eval_dynamics::Function) where T =
+  DILS{T}([input_poly], [output_poly], nstates, forward_readout,
+          backward_readout, eval_dynamics)
+
+
+
+input_polys(d::DILS) = d.input_polys
+output_polys(d::DILS) = d.output_polys
+input_poly(d::DILS) = tensor(input_polys(d))
+output_poly(d::DILS) = tensor(output_polys(d))
+input_poly(d::DILS, i::Int) = input_polys(d)[i]
+output_poly(d::DILS, i::Int) = output_polys(d)[i]
+
+nstates(d::DILS) = d.nstates
+
+""" forward_readout(d::DILS, u::AbstractVector, x::AbstractVector)
+
+This is half of the map on positions for the DILS `d` given by
+  ``Sy^S \\to [p,q]``. It takes a state ``u \\in \\mathbb{R}^S``
+and a position in the input polynomial ``x \\in p(1)`` and returns a position of
+the output polynomial.
+"""
+
+forward_readout(d::DILS, u::AbstractVector, x::AbstractVector) =
+  d.forward_readout(u, x)
+
+""" backward_readout(d::DILS, u::AbstractVector, x::AbstractVector, y::AbstractVector)
+
+This is the other half of the map on positions. It takes a state, a position on
+the input polynomial and a position on the output polynomial and returns a
+direction on the input polynomial.
+
+``\\sum_{s \\in S} \\sum_{i\\in p(1)} q[f(s,i)] \\rightarrow p[i]``
+"""
+
+backward_readout(d::DILS, u::AbstractVector, x::AbstractVector, y::AbstractVector) =
+  d.backward_readout(u, x, y)
+
+""" eval_dynamics(d::DILS, u::AbstractVector, x::AbstractVector, y::AbstractVector)
+
+This is the map on directions of the DILS `d` given by
+``{\\mathbb{R}^S}y^{\\mathbb{R}^S} \\to [p,q]`. It takes a state
+``u \\in \\mathbb{R}^S``, a position ``x \\in p(1)``, a direction
+``y \\in q[f(u,x)]`` (where ``f`` is the `forward_readout`)
+and returns an updated state in ``\\mathbb{R}^S``
+"""
+
+eval_dynamics(d::DILS, u::AbstractVector, x::AbstractVector,
+              y::AbstractVector) = d.eval_dynamics(u, x, y)
+
+isa_state(d::DILS{T}, u::AbstractVector{T}) where T = length(u) == nstates(d)
+
+
+# A DILS based on a function
+function DILS{T}(f::Function, ninput::Int, noutput::Int; monomial=true) where T
+  DILS{T}(
+    monomial ? 0=>ninput : [0=>1 for _ in 1:ninput],  # input poly is ninput y^0
+    monomial ? 0=>noutput : [0=>1 for _ in 1:noutput], # out poly is noutput y^0
+    0,          # trivial state
+    (_, x) -> f(x),  #forward map
+    (_, _, _) -> T[], # backward map
+    (_, _, _) -> T[]  # eval dynamics
+  )
+end
+
+# A DILS that outputs constant values
+DILS{T}(c::AbstractVector; monomial=true) where T =
+  DILS{T}(_ -> c, 0, length(c); monomial=monomial)
+
+"""
+Basically a (dependent) lens
+"""
+function DILS{T}(forward::Function, backward::Function, ninput::Int,
+                 noutput::Int) where T
+  DILS{T}(
+    [ninput=>ninput],
+    [noutput=>noutput],
+    0,                     # trivial state
+    (_, x) -> forward(x),  #forward map
+    (_, x, y) -> backward(x,y), # backward map
+    (_, _, _) -> T[]  # eval dynamics
+    )
+end
+
+
+"""
+Compose two DILS in sequences
+
+TODO: Change braiding to be any wiring pattern from d1 to d2.
+If the output poly and input poly are of the form Ay^A (i.e. same inputs and outputs) then
+you can do this by constructing an integer-valued matrix based on the wiring pattern encoding copying, merging, delete, and create.
+For the forward pass multiply the positions of the output of d1 by the matrix.
+For the backward pass, multiply the directions of the input of d2 by the transpose of the matrix.
+"""
+
+function sequential_compose(d1::DILS{T}, d2::DILS{T},
+    weights::Union{Nothing,AbstractMatrix}=nothing)::DILS{T} where T
+
+  if isnothing(weights)
+    op, ip = output_polys(d1), input_polys(d2)
+    if op == ip
+      weights = I(length(op))
+    else
+      error("Cannot infer sequential compose weights between $op and $ip")
+    end
+  end
+
+  # check that the monomials match
+  for (i,k) in pairs(weights)
+    if k > 0
+      input_poly(d2, i[1]) == output_poly(d1, i[2])
+    end
+  end
+  d1_state = u -> u[1:d1.nstates]
+  d2_state = u -> u[d1.nstates+1:end]
+  d1_output = (u,x) -> forward_readout(d1, d1_state(u), x)
+  d2_backward = (u,x,y) -> backward_readout(d2, d2_state(u), d1_output(u,x), y)
+
+  # This produces the matrices that copy and add the things flowing along the wires correctly.
+  forward_matrix = generalized_kroenecker(weights,
+                                          map(positions, input_polys(d2)),
+                                          map(positions, output_polys(d1)))
+  backward_matrix = generalized_kroenecker(weights',
+                                           map(directions, output_polys(d1)),
+                                           map(directions, input_polys(d2)))
+  return DILS{T}(
+    input_polys(d1), # input monomials
+    output_polys(d2), # output monomials
+    nstates(d1)+nstates(d2), # number of states/dimensions
+    (u,x) -> forward_readout(d2, d2_state(u), forward_matrix*d1_output(u,x)),
+    (u,x,y) -> backward_readout(d1, d1_state(u), x, backward_matrix*d2_backward(u,x,y)),
+    (u,x,y) -> vcat(eval_dynamics(d1, d1_state(u), x, backward_matrix'*d2_backward(u, x, y)), # updated state for d1
+                    eval_dynamics(d2, d2_state(u), forward_matrix*d1_output(u,x), y))       # updated state for d2
+  )
+end
+
+"""TODO TEST"""
+function unit(::Type{DILS{T}}) where T
+  function fun2(u, x)
+    all(isempty.([u,x])) || error("Bad input")
+    T[]
+  end
+  function fun3(u, x, y)
+    all(isempty.([u,x,y])) || error("Bad input")
+    T[]
+  end
+
+  DILS{T}(Monomial[],Monomial[], 0, fun2, fun3, fun3)
+end
+
+"""TODO TEST BACKWARD READOUT AND DYNAMICS"""
+function tensor(d1::DILS{T}, d2::DILS{T}) where T
+  d1_state = u -> u[1:d1.nstates]
+  d2_state = u -> u[d1.nstates+1:end]
+  d1_in = x -> x[1:positions(input_poly(d1))]
+  d2_in = x -> x[positions(input_poly(d1))+1:end]
+  d1_out = y-> y[1:directions(output_poly(d1))]
+  d2_out = y-> y[directions(output_poly(d1))+1:end]
+  DILS{T}(vcat(input_polys(d1), input_polys(d2)),
+          vcat(output_polys(d1), output_polys(d2)),
+          d1.nstates + d2.nstates,
+          (u,x) -> vcat(forward_readout(d1, d1_state(u), d1_in(x)),
+                        forward_readout(d2, d2_state(u), d2_in(x))),
+          (u,x,y) -> vcat(backward_readout(d1, d1_state(u), d1_in(x), d1_out(y)),
+                          backward_readout(d2, d2_state(u), d2_in(x), d2_out(y))),
+          (u,x,y) -> vcat(eval_dynamics(d1, d1_state(u), d1_in(x),d1_out(y)),
+                          eval_dynamics(d2, d2_state(u), d2_in(x),d2_out(y)))
+  )
+end
+
+"""TODO TEST"""
+function tensor(ds::AbstractVector{<:DILS{T}}) where T
+  res = unit(DILS{T})
+  for d in ds
+    res = tensor(res, d)
+  end
+  res
+end # reduce(tensor, ds, unit(DILS{T}))
+
+function oapply(d::WiringDiagram, dils::Vector{<:DILS})::DILS
+  # check that dils[i] "fills" box i of d
+  @assert length(dils) == nboxes(d)
+  for (dil, box) in zip(dils, 1:nboxes(d))
+    @assert length(input_ports(d, box)) == length(input_polys(dil))
+    @assert length(output_ports(d, box)) == length(output_polys(dil))
+  end
+
+  layers = layerize(internal_graph(d))
+
+  t_layers = [tensor([dils[i] for i in layer]) for layer in layers]
+  res = t_layers[1]
+  for (i, (layer, dils_layer)) in enumerate(zip(layers[2:end], t_layers[2:end]))
+    weights = hcat(hcat(map(layer) do br # the right layer
+      vcat(map(incident(d.diagram, br, :in_port_box)) do pr
+        w_pr = incident(d.diagram, pr, :tgt)
+        map(layers[i]) do bl
+          map(incident(d.diagram, bl, :out_port_box)) do pl
+            length(w_pr ∩ incident(d.diagram, pl, :src))
+          end
+        end
+      end...)
+    end...)...)
+    res = sequential_compose(res, dils_layer, weights)
+  end
+  return res
+end
+
+
+"""
+Partition vertices of a graph by dependency. Get ordered layers.
+Can be used for composing morphisms in a SMC.
+"""
+function layerize(G::AbstractGraph)::Vector{Vector{Int}}
+  seen = Set{Int}()
+  layers = Vector{Int}[]
+  while length(seen) < nv(G)
+    next_layer = Int[]
+    for v in filter(x->x∉ seen, parts(G, :V))
+      if G[incident(G, v, :tgt) , :src] ⊆ seen
+        push!(next_layer, v)
+      end
+    end
+    push!(layers, next_layer)
+    union!(seen, next_layer)
+  end
+  return layers
+end
+
+
+"""
+Multiply each partition of a i×j partitioned matrix by a scalar, taken from a
+i×j scalar matrix.
+"""
+function generalized_kroenecker(weights::AbstractMatrix, rowdims::Vector{Int}, coldims::Vector{Int})
+  blockrow = map(enumerate(coldims)) do (col, coldim)
+    block = map(enumerate(rowdims)) do (row, rowdim)
+        weights[row,col] * (rowdim == coldim ? I(rowdim) : zeros(rowdim, coldim))
+    end
+    vcat(block...)
+  end
+  return hcat(blockrow...)
+end
+
+
+
+end # module
\ No newline at end of file
diff --git a/test/dwd_dils.jl b/test/dwd_dils.jl
new file mode 100644
index 0000000..8082041
--- /dev/null
+++ b/test/dwd_dils.jl
@@ -0,0 +1,164 @@
+using Revise
+using AlgebraicDynamics.DWDDILS
+
+using Test
+using LinearAlgebra
+using Plots: plot
+
+using Catlab.WiringDiagrams.DirectedWiringDiagrams
+
+
+
+# Functional DILS
+#################
+
+minus = DILS{Float64}(x->[x[1] - x[2]], 2, 1; monomial=false)
+constant = DILS{Float64}([1,2]; monomial=false)
+
+@test length(input_polys(minus)) == 2
+
+# Sequential composition
+########################
+one_minus_two = sequential_compose(constant, minus, [1 0; 0 1])
+two_minus_one = sequential_compose(constant, minus, [0 1; 1 0])
+two_minus_two = sequential_compose(constant, minus, [0 1; 2 0])
+
+@test forward_readout(one_minus_two, [], []) == [-1]
+@test forward_readout(two_minus_one, [], []) == [1]
+@test forward_readout(two_minus_two, [], []) == [0]
+@test backward_readout(one_minus_two, [], [], []) == []
+
+# Unit
+######
+U = unit(DILS{Float64})
+@test isempty(forward_readout(U, [], []))
+@test isempty(backward_readout(U, [], [], []))
+
+t = tensor(unit(typeof(one_minus_two)), one_minus_two);
+@test nstates(t) == nstates(one_minus_two)
+@test input_poly(t) == input_poly(one_minus_two)
+@test output_poly(t) == output_poly(one_minus_two)
+@test forward_readout(t, [], []) == [-1]
+@test backward_readout(t, [], [], []) == []
+
+# Oapply
+########
+
+d = WiringDiagram([], [:X])
+b1 = add_box!(d, Box(:const, [], [:one, :two]))
+b2 = add_box!(d, Box(:minus, [:x, :y], [:x_minus_y]))
+add_wires!(d, Pair[
+  (b1, 1) => (b2, 2),
+  (b1, 2) => (b2, 1),
+  (b1, 1) => (b2, 1),
+  (b2, 1) => (output_id(d), 1)
+])
+
+dils = [constant, minus]
+res = oapply(d, dils);
+@test forward_readout(res, [], []) == [2.]
+
+# NN
+####
+
+# A DILS with 1 -> [Ry^R,Ry^R] representing "relu" with gradient descent
+relu = DILS{Float64}(
+    (x) -> [max(0, x[1])],          # relu
+    (x,y)-> x[1] < 0 ? [0] : y,     # gradient of relu
+    1, 1
+)
+
+@test forward_readout(relu, [], [-1]) == [0]
+@test forward_readout(relu, [], [1]) == [1]
+@test backward_readout(relu, [], [10], [1]) == [1]
+@test backward_readout(relu, [], [-10], [1]) == [0]
+
+"""Gradient descent dynamics
+
+An unactivated_neuron with n inputs is a DILS:
+ R^{n+1} y^{R^{n+1}} -> [R^n y^{R^n}, Ry^R]
+"""
+function unactivated_neuron(n_inputs::Int)
+  DILS{Float64}(
+    [n_inputs => n_inputs],           # input
+    [1 => 1],                         # outputs
+    n_inputs+1,                       # state (bias is last element)
+    (u, x) -> [dot(u, vcat(x,[1]))],  #forward readout
+    (u, x, Δy) -> Δy .* u[1:end-1],          # backward readout - returns Δx
+    (u, x, Δy) -> Δy .* vcat(x, [1])         # eval_dynamics - returns Δu
+  )
+end
+
+"""Note, this tests the default sequential_compose (no matrix given)"""
+activated_neuron(ninputs::Int) =
+    sequential_compose(unactivated_neuron(ninputs), relu)
+
+my_fun(x) = [max(0, x[1] + 2*x[2] + 5)]
+my_fun(x,y) = my_fun([x,y])
+
+function training_data(f::Function, n_steps::Int=1000, low::Float64=-10., high::Float64=10.)
+  return map(1:n_steps) do _
+    x = [rand() * (high-low) + low for _ in 1:2]
+    x => f(x)
+  end
+end
+
+
+mutable struct SingleNeuron
+    weights::Vector
+    bias::Vector
+    d::DILS{Float64}
+end
+
+
+function initialize!(neuron::SingleNeuron)
+  neuron.weights = rand(2)
+  neuron.bias = rand(1)
+end
+
+function train!(neuron::SingleNeuron, training_data,
+                α::Float64=0.01)
+  ys = Float64[]
+  for (x,y) in training_data
+    u = vcat(neuron.weights, neuron.bias)
+    y′ =forward_readout(neuron.d, u, x)
+    Δy = y′.- y
+    push!(ys, Δy[1])
+    new_u = u - eval_dynamics(neuron.d, u, x, α*Δy)
+    neuron.weights = new_u[1:2]
+    neuron.bias = [new_u[3]]
+  end
+  return ys
+end
+
+(s::SingleNeuron)(x,y) = forward_readout(s.d, vcat(s.weights, s.bias), [x,y])
+
+
+weights = rand(2)
+bias = rand(1)
+nn = SingleNeuron(weights, bias, activated_neuron(2))
+nn(1.,1.)
+
+
+data = training_data(my_fun, 10000)
+initialize!(nn)
+errs = train!(nn, data)
+
+nonzero_errs = collect(filter(x->x!=(0),errs))
+plot(1:length(nonzero_errs), nonzero_errs)
+
+# plot(first.(data), last.(data)) # only if data is 1-D
+
+
+# MISC
+######
+@test generalized_kroenecker([1 0; 1 0; 0 1], [3,3,2],[3, 2]) == [
+    1 0 0 0 0;
+    0 1 0 0 0;
+    0 0 1 0 0;
+    1 0 0 0 0;
+    0 1 0 0 0;
+    0 0 1 0 0;
+    0 0 0 1 0;
+    0 0 0 0 1
+]