From e0ac22a134136efc9752b7747a1d796ae6014008 Mon Sep 17 00:00:00 2001 From: Sophie Libkind <4992461+slibkind@users.noreply.github.com> Date: Mon, 24 Jan 2022 15:09:27 -0800 Subject: [PATCH 1/5] ehn: first draft of composing DILS --- src/dwd_dils.jl | 166 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 166 insertions(+) create mode 100644 src/dwd_dils.jl diff --git a/src/dwd_dils.jl b/src/dwd_dils.jl new file mode 100644 index 0000000..2d452e7 --- /dev/null +++ b/src/dwd_dils.jl @@ -0,0 +1,166 @@ +using Catlab.WiringDiagrams.DirectedWiringDiagrams +using LinearAlgebra +using Test +using Plots: plot + + +const Monomial = Pair{Int,Int} + + +struct DILS{T} + input_poly::Monomial + output_poly::Monomial + nstates::Int # dimension of the state space + forward_readout::Function + backward_readout::Function + eval_dynamics::Function +end + + +input_poly(d::DILS) = d.input_poly +output_poly(d::DILS) = d.output_poly +nstates(d::DILS) = d.nstates + +forward_readout(d::DILS, u::AbstractVector, x::AbstractVector) = + d.forward_readout(u, x) +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) = + d.eval_dynamics(u, x, y) + +isa_state(d::DILS{T}, u::AbstractVector{T}) where T = + length(u) == nstates(d) + + +function DILS{T}(f::Function, ninput::Int, noutput::Int) where T + DILS{T}( + 0=>ninput, + 0=>noutput, + 0, + (_, x) -> f(x)), #forward map + (_, _, _) -> T[], # backward map + (_, _, _) -> T[] # eval dynamics +end + + +""" +forward(x):: +backward(x,y) +""" +function DILS{T}(forward::Function, backward::Function, ninput::Int, noutput::Int) where T + DILS{T}( + ninput=>ninput, + noutput=>noutput, + 0, + (_, x) -> forward(x), #forward map + (_, x, y) -> backward(x,y), # backward map + (_, _, _) -> T[] # eval dynamics + ) +end + +relu = DILS{Float64}( + (x) -> [max(0, x[1])], + (x,y)-> x[1] < 0 ? [0] : y, + 1, 1 +) + +@test forward_readout(relu, [], [-1]) == [0] + +"""Gradient descent dynamics""" +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 + + +""" +Compose two DILS in sequences +""" +function sequential_compose(d1::DILS{T}, d2::DILS{T})::DILS{T} where T + 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) + return DILS{T}( + input_poly(d1), # input monomial + output_poly(d2), # output monomial + nstates(d1)+nstates(d2), # number of states/dimensions + (u,x) -> forward_readout(d2, d2_state(u), d1_output(u,x)), + (u,x,y) -> backward_readout(d1, d1_state(u), x, d2_backward(u,x,y)), + (u,x,y) -> vcat(eval_dynamics(d1, d1_state(u), x, d2_backward(u, x, y)), # updated state for d1 + eval_dynamics(d2, d2_state(u), d1_output(u,x), y)) # updated state for d2 + ) +end + +activated_neuron(ninputs::Int) = sequential_compose(unactivated_neuron(ninputs), relu) + + + + +mutable struct SingleNeuron + weights::Vector + bias::Vector + d::DILS{Float64} +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.) + + +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 + +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 + + +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)) + + +function oapply(d::WiringDiagram, dils::Vector{DILS})::DILS +end + + From 03fb03af684f814c7b2922b1f0ec574e2dfcd022 Mon Sep 17 00:00:00 2001 From: Sophie Libkind <4992461+slibkind@users.noreply.github.com> Date: Mon, 14 Mar 2022 15:31:19 -0700 Subject: [PATCH 2/5] Added doc strings and braiding --- src/dwd_dils.jl | 170 ++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 143 insertions(+), 27 deletions(-) diff --git a/src/dwd_dils.jl b/src/dwd_dils.jl index 2d452e7..c8414db 100644 --- a/src/dwd_dils.jl +++ b/src/dwd_dils.jl @@ -3,28 +3,66 @@ using LinearAlgebra using Test using Plots: plot +""" +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})::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_poly::Monomial - output_poly::Monomial + 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) = d.input_poly -output_poly(d::DILS) = d.output_poly 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) @@ -32,46 +70,56 @@ 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) where T DILS{T}( - 0=>ninput, - 0=>noutput, - 0, - (_, x) -> f(x)), #forward map + 0=>ninput, # input poly is ninput y^0 + 0=>noutput, # input poly is nouput 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) where T = DILS{T}(_ -> c, 0, length(c)) """ -forward(x):: -backward(x,y) +Basically a (dependent) lens """ function DILS{T}(forward::Function, backward::Function, ninput::Int, noutput::Int) where T DILS{T}( - ninput=>ninput, - noutput=>noutput, - 0, + [ninput=>ninput], + [noutput=>noutput], + 0, # trivial state (_, x) -> forward(x), #forward map (_, x, y) -> backward(x,y), # backward map (_, _, _) -> T[] # eval dynamics ) end +# A DILS with 1 -> [Ry^R,Ry^R] representing "relu" with gradient descent relu = DILS{Float64}( - (x) -> [max(0, x[1])], - (x,y)-> x[1] < 0 ? [0] : y, + (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""" +"""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) + [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 @@ -79,27 +127,46 @@ function unactivated_neuron(n_inputs::Int) 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})::DILS{T} where T +function sequential_compose(d1::DILS{T}, d2::DILS{T}, braiding::Union{Nothing,AbstractVector} = nothing)::DILS{T} where T + # @assert Base.isperm(braiding) + #assert positions(output_poly(d1)) == positions(input_poly(d2)) # TODO: check directions match + braiding = isnothing(braiding) ? (1:positions(output_poly(d1))) : braiding 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) return DILS{T}( - input_poly(d1), # input monomial - output_poly(d2), # output monomial + 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), d1_output(u,x)), - (u,x,y) -> backward_readout(d1, d1_state(u), x, d2_backward(u,x,y)), - (u,x,y) -> vcat(eval_dynamics(d1, d1_state(u), x, d2_backward(u, x, y)), # updated state for d1 - eval_dynamics(d2, d2_state(u), d1_output(u,x), y)) # updated state for d2 + (u,x) -> forward_readout(d2, d2_state(u), d1_output(u,x)[braiding]), + (u,x,y) -> backward_readout(d1, d1_state(u), x, d2_backward(u,x,y)[invperm(braiding)]), + (u,x,y) -> vcat(eval_dynamics(d1, d1_state(u), x, d2_backward(u, x, y)[invperm(braiding)]), # updated state for d1 + eval_dynamics(d2, d2_state(u), d1_output(u,x)[braiding], y)) # updated state for d2 ) end -activated_neuron(ninputs::Int) = sequential_compose(unactivated_neuron(ninputs), relu) +minus = DILS{Float64}(x->[x[1] - x[2]], 2, 1) + +constant = DILS{Float64}([1,2]) + +one_minus_two = sequential_compose(constant, minus) +two_minus_one = sequential_compose(constant, minus, [2,1]) +forward_readout(one_minus_two, [], []) == [-1] +forward_readout(two_minus_one, [], []) == [1] + +activated_neuron(ninputs::Int) = sequential_compose(unactivated_neuron(ninputs), relu) @@ -159,8 +226,57 @@ plot(1:length(nonzero_errs), nonzero_errs) plot(first.(data), last.(data)) +"""TODO TEST""" +function unit(::DILS{T}) + DILS{T}(Monomial[],Monomial[], 0, id, id, id) +end + +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] + DILS{T}(tensor(input_polys(d1), input_polys(d2)), + tensor(output_polys(d1), output_polys(d2)), + d1.nstates + d2.nstates, + (u,x) -> vcat(forward_readout(d1, u[1:d1.nstates], x[1:first(input_poly(d1))], forward_readout(d2, ))) + ) +end + +"""TODO TEST""" +tensor(ds::Vector(DILS{T})) where {T} = 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, boxes(d)) + @assert length(input_ports(d, box)) == length(input_polys(dil)) + @assert length(input_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] + return reduce(compose, t_layers) +end + + +""" +Partition vertices of a graph by dependency. Get ordered layers. +Can be used for composing morphisms in a SMC. +""" +function layerize(G::Graph)::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 From 88c33fd4843d64235ec59c178008e9d651264d49 Mon Sep 17 00:00:00 2001 From: Sophie Libkind <4992461+slibkind@users.noreply.github.com> Date: Mon, 21 Mar 2022 20:42:44 -0700 Subject: [PATCH 3/5] sequential composition --- src/dwd_dils.jl | 111 +++++++++++++++++++++++++++++++++++++----------- 1 file changed, 87 insertions(+), 24 deletions(-) diff --git a/src/dwd_dils.jl b/src/dwd_dils.jl index c8414db..52d5434 100644 --- a/src/dwd_dils.jl +++ b/src/dwd_dils.jl @@ -1,4 +1,6 @@ using Catlab.WiringDiagrams.DirectedWiringDiagrams +using Catlab.Graphs, Catlab.CategoricalAlgebra + using LinearAlgebra using Test using Plots: plot @@ -7,8 +9,9 @@ using Plots: plot 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})::Monomial = sum(first.(ms)) => sum(last.(ms)) +tensor(ms::Vector{Monomial}, ns::Vector{Monomial}) = sum(first.(ms)) => sum(last.(ms)) positions(m::Monomial) = m.second directions(m::Monomial) = m.first @@ -71,10 +74,10 @@ isa_state(d::DILS{T}, u::AbstractVector{T}) where T = # A DILS based on a function -function DILS{T}(f::Function, ninput::Int, noutput::Int) where T +function DILS{T}(f::Function, ninput::Int, noutput::Int; monomial=true) where T DILS{T}( - 0=>ninput, # input poly is ninput y^0 - 0=>noutput, # input poly is nouput y^0 + monomial ? 0=>ninput : [0=>1 for _ in 1:ninput], # input poly is ninput y^0 + monomial ? 0=>noutput : [0=>1 for _ in 1:noutput], # input poly is nouput y^0 0, # trivial state (_, x) -> f(x), #forward map (_, _, _) -> T[], # backward map @@ -83,7 +86,7 @@ function DILS{T}(f::Function, ninput::Int, noutput::Int) where T end # A DILS that outputs constant values -DILS{T}(c::AbstractVector) where T = DILS{T}(_ -> c, 0, length(c)) +DILS{T}(c::AbstractVector; monomial=true) where T = DILS{T}(_ -> c, 0, length(c); monomial=monomial) """ Basically a (dependent) lens @@ -137,34 +140,68 @@ you can do this by constructing an integer-valued matrix based on the wiring pat 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}, braiding::Union{Nothing,AbstractVector} = nothing)::DILS{T} where T - # @assert Base.isperm(braiding) - #assert positions(output_poly(d1)) == positions(input_poly(d2)) # TODO: check directions match - braiding = isnothing(braiding) ? (1:positions(output_poly(d1))) : braiding + +# function sequential_compose(d1::DILS{T}, d2::DILS{T}, braiding::Nothing)::DILS{T} where T +# sequential_compose(d1::DILS{T}) +# end + +input_poly(d::DILS, i::Int) = input_polys(d)[i] +output_poly(d::DILS, i::Int) = output_polys(d)[i] + +function sequential_compose(d1::DILS{T}, d2::DILS{T}, weights::AbstractMatrix)::DILS{T} where T + # 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), d1_output(u,x)[braiding]), - (u,x,y) -> backward_readout(d1, d1_state(u), x, d2_backward(u,x,y)[invperm(braiding)]), - (u,x,y) -> vcat(eval_dynamics(d1, d1_state(u), x, d2_backward(u, x, y)[invperm(braiding)]), # updated state for d1 - eval_dynamics(d2, d2_state(u), d1_output(u,x)[braiding], y)) # updated state for d2 + (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 -minus = DILS{Float64}(x->[x[1] - x[2]], 2, 1) +function generalized_kroenecker(weights::AbstractMatrix, rowdims::Vector{Int}, coldims::Vector{Int}) + + block = map(enumerate(coldims)) do (col, coldim) + vcat(map(enumerate(rowdims)) do (row, rowdim) + weights[row,col] * (rowdim == coldim ? I(rowdim) : zeros(rowdim, coldim)) + end...) + end + return hcat(block...) +end + +generalized_kroenecker([1 0; 1 0; 0 1], [3,3,2],[3, 2]) + +minus = DILS{Float64}(x->[x[1] - x[2]], 2, 1; monomial=false) +@assert length(input_polys(minus)) == 2 + +constant = DILS{Float64}([1,2]; monomial=false) -constant = DILS{Float64}([1,2]) +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]) -one_minus_two = sequential_compose(constant, minus) -two_minus_one = sequential_compose(constant, minus, [2,1]) forward_readout(one_minus_two, [], []) == [-1] forward_readout(two_minus_one, [], []) == [1] +forward_readout(two_minus_two, [], []) == [0] + +backward_readout(one_minus_two, [], [], []) == [] activated_neuron(ninputs::Int) = sequential_compose(unactivated_neuron(ninputs), relu) @@ -227,22 +264,22 @@ plot(1:length(nonzero_errs), nonzero_errs) plot(first.(data), last.(data)) """TODO TEST""" -function unit(::DILS{T}) +function unit(::DILS{T}) where T DILS{T}(Monomial[],Monomial[], 0, id, id, id) end 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] - DILS{T}(tensor(input_polys(d1), input_polys(d2)), - tensor(output_polys(d1), output_polys(d2)), + 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, u[1:d1.nstates], x[1:first(input_poly(d1))], forward_readout(d2, ))) + (u,x) -> vcat(forward_readout(d1, d1_state(u), x[1:positions(input_poly(d1))], forward_readout(d2, d2_state(u), x[positions(input_poly(d1))+1,end]))) ) end """TODO TEST""" -tensor(ds::Vector(DILS{T})) where {T} = reduce(tensor, ds, unit(DILS{T}})) +tensor(ds::Vector{D}) where {T, D<:DILS{T}} = reduce(tensor, ds, unit(DILS{T})) function oapply(d::WiringDiagram, dils::Vector{DILS})::DILS # check that dils[i] "fills" box i of d @@ -255,15 +292,41 @@ function oapply(d::WiringDiagram, dils::Vector{DILS})::DILS layers = layerize(internal_graph(d)) t_layers = [tensor(dils[i] for i in layer) for layer in layers] - return reduce(compose, t_layers) + res = t_layers[1] + for (i, layer) in t_layers[2:end] + weights = map(in_ports(d, )) + res = compose(res, layer, weights) + end + return res + +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) +]) + +g = internal_graph(d) +layers = layerize(internal_graph(d)) + +dils = [constant, minus] + +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 +oapply(d, dils) + """ Partition vertices of a graph by dependency. Get ordered layers. Can be used for composing morphisms in a SMC. """ -function layerize(G::Graph)::Vector{Vector{Int}} +function layerize(G::AbstractGraph)::Vector{Vector{Int}} seen = Set{Int}() layers = Vector{Int}[] while length(seen) < nv(G) From c594dd16ed89b6972f01c5352fb81ee74d203c6e Mon Sep 17 00:00:00 2001 From: Kris Date: Wed, 23 Mar 2022 08:56:04 -0400 Subject: [PATCH 4/5] finish oapply, move tests to separate file, default seq compose, Fix tensor monomial Fix LaTeX Fix unit DILS (can't use 'id' for the do-nothing functions b/c wrong arity) --- .gitignore | 3 +- Project.toml | 36 ++--- src/AlgebraicDynamics.jl | 1 + src/dwd_dils.jl | 290 +++++++++++++++------------------------ test/dwd_dils.jl | 164 ++++++++++++++++++++++ 5 files changed, 296 insertions(+), 198 deletions(-) create mode 100644 test/dwd_dils.jl 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 b86feb0..5fdfe10 100644 --- a/Project.toml +++ b/Project.toml @@ -1,24 +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" -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" - [compat] Catlab = "0.11, 0.12, 0.13" Compose = "^0.9.1" @@ -35,7 +19,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 cc97226..7dc078e 100644 --- a/src/AlgebraicDynamics.jl +++ b/src/AlgebraicDynamics.jl @@ -9,5 +9,6 @@ include("uwd_dynam.jl") include("dwd_dynam.jl") include("cpg_dynam.jl") include("trajectories.jl") +include("dwd_dils.jl") end # module diff --git a/src/dwd_dils.jl b/src/dwd_dils.jl index 52d5434..204b42e 100644 --- a/src/dwd_dils.jl +++ b/src/dwd_dils.jl @@ -1,17 +1,24 @@ +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 -using LinearAlgebra -using Test -using Plots: plot -""" +""" -A monomial is a pair of integers A => B representing the poly ``\mathbb{R}^By^{\mathbb{R}^A}`` +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}, ns::Vector{Monomial}) = sum(first.(ms)) => sum(last.(ms)) +tensor(ms::Vector{Monomial}) = sum(first.(ms)) => sum(last.(ms)) positions(m::Monomial) = m.second directions(m::Monomial) = m.first @@ -27,8 +34,8 @@ struct DILS{T} 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::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) @@ -37,23 +44,25 @@ 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. +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 +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]`` +``\\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) = @@ -61,9 +70,9 @@ backward_readout(d::DILS, u::AbstractVector, x::AbstractVector, y::AbstractVecto """ 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`` +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) = @@ -102,56 +111,32 @@ function DILS{T}(forward::Function, backward::Function, ninput::Int, noutput::In ) end -# 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 - - """ 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. +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}, braiding::Nothing)::DILS{T} where T -# sequential_compose(d1::DILS{T}) -# end +function sequential_compose(d1::DILS{T}, d2::DILS{T}, + weights::Union{Nothing,AbstractMatrix}=nothing)::DILS{T} where T -input_poly(d::DILS, i::Int) = input_polys(d)[i] -output_poly(d::DILS, i::Int) = output_polys(d)[i] + 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 -function sequential_compose(d1::DILS{T}, d2::DILS{T}, weights::AbstractMatrix)::DILS{T} where T # check that the monomials match - for (i,k) in pairs(weights) - if k > 0 + for (i,k) in pairs(weights) + if k > 0 input_poly(d2, i[1]) == output_poly(d1, i[2]) end end @@ -161,8 +146,12 @@ function sequential_compose(d1::DILS{T}, d2::DILS{T}, weights::AbstractMatrix):: 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))) + 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 @@ -174,154 +163,78 @@ function sequential_compose(d1::DILS{T}, d2::DILS{T}, weights::AbstractMatrix):: ) end - -function generalized_kroenecker(weights::AbstractMatrix, rowdims::Vector{Int}, coldims::Vector{Int}) - - block = map(enumerate(coldims)) do (col, coldim) - vcat(map(enumerate(rowdims)) do (row, rowdim) - weights[row,col] * (rowdim == coldim ? I(rowdim) : zeros(rowdim, coldim)) - end...) - end - return hcat(block...) -end - -generalized_kroenecker([1 0; 1 0; 0 1], [3,3,2],[3, 2]) - -minus = DILS{Float64}(x->[x[1] - x[2]], 2, 1; monomial=false) -@assert length(input_polys(minus)) == 2 - -constant = DILS{Float64}([1,2]; monomial=false) - -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]) - -forward_readout(one_minus_two, [], []) == [-1] -forward_readout(two_minus_one, [], []) == [1] -forward_readout(two_minus_two, [], []) == [0] - -backward_readout(one_minus_two, [], [], []) == [] - -activated_neuron(ninputs::Int) = sequential_compose(unactivated_neuron(ninputs), relu) - - - -mutable struct SingleNeuron - weights::Vector - bias::Vector - d::DILS{Float64} -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.) - - -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) +"""TODO TEST""" +function unit(::Type{DILS{T}}) where T + function fun2(u, x) + all(isempty.([u,x])) || error("Bad input") + T[] end -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]] + function fun3(u, x, y) + all(isempty.([u,x,y])) || error("Bad input") + T[] end - return ys -end - -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)) - -"""TODO TEST""" -function unit(::DILS{T}) where T - DILS{T}(Monomial[],Monomial[], 0, id, id, id) + 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] - DILS{T}(vcat(input_polys(d1), input_polys(d2)), + 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), x[1:positions(input_poly(d1))], forward_readout(d2, d2_state(u), x[positions(input_poly(d1))+1,end]))) + 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""" -tensor(ds::Vector{D}) where {T, D<:DILS{T}} = reduce(tensor, ds, unit(DILS{T})) +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 +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, boxes(d)) + for (dil, box) in zip(dils, 1:nboxes(d)) @assert length(input_ports(d, box)) == length(input_polys(dil)) - @assert length(input_ports(d, box)) == length(output_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] + t_layers = [tensor([dils[i] for i in layer]) for layer in layers] res = t_layers[1] - for (i, layer) in t_layers[2:end] - weights = map(in_ports(d, )) - res = compose(res, layer, weights) + 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 - -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) -]) - -g = internal_graph(d) -layers = layerize(internal_graph(d)) - -dils = [constant, minus] - -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)) + return res end -oapply(d, dils) - """ Partition vertices of a graph by dependency. Get ordered layers. Can be used for composing morphisms in a SMC. @@ -343,3 +256,20 @@ function layerize(G::AbstractGraph)::Vector{Vector{Int}} 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 +] From 678b3e148d795cf2c5bd0ab13a10965701221ffb Mon Sep 17 00:00:00 2001 From: Kris Date: Wed, 23 Mar 2022 13:47:37 -0400 Subject: [PATCH 5/5] Respecting 80 char line limit --- src/dwd_dils.jl | 39 ++++++++++++++++++++++++--------------- 1 file changed, 24 insertions(+), 15 deletions(-) diff --git a/src/dwd_dils.jl b/src/dwd_dils.jl index 204b42e..8669890 100644 --- a/src/dwd_dils.jl +++ b/src/dwd_dils.jl @@ -35,8 +35,10 @@ struct DILS{T} 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) + forward_readout::Function, backward_readout::Function, + eval_dynamics::Function) where T = + DILS{T}([input_poly], [output_poly], nstates, forward_readout, + backward_readout, eval_dynamics) @@ -51,8 +53,10 @@ 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. +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) = @@ -60,8 +64,10 @@ forward_readout(d::DILS, u::AbstractVector, x::AbstractVector) = """ 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. +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]`` """ @@ -70,23 +76,24 @@ backward_readout(d::DILS, u::AbstractVector, x::AbstractVector, y::AbstractVecto """ 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`) +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) +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) +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], # input poly is nouput 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 @@ -95,12 +102,14 @@ function DILS{T}(f::Function, ninput::Int, noutput::Int; monomial=true) where T end # A DILS that outputs constant values -DILS{T}(c::AbstractVector; monomial=true) where T = DILS{T}(_ -> c, 0, length(c); monomial=monomial) +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 +function DILS{T}(forward::Function, backward::Function, ninput::Int, + noutput::Int) where T DILS{T}( [ninput=>ninput], [noutput=>noutput],