diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml deleted file mode 100644 index b6df2f6c..00000000 --- a/.github/workflows/ci.yml +++ /dev/null @@ -1,46 +0,0 @@ -name: CI - -on: - push: - paths-ignore: - - 'LICENSE' - - 'README.md' - branches: - - master - -jobs: - test: - name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} - runs-on: ${{ matrix.os }} - if: "!contains(github.event.head_commit.message, 'skip ci')" - strategy: - fail-fast: false - matrix: - version: - - '1' - os: - - ubuntu-latest - arch: - - x64 - env: - JULIA_PKG_SERVER: '' - steps: - - name: Check out repository - uses: actions/checkout@v4 - - name: Set up Julia - uses: julia-actions/setup-julia@v1 - with: - version: ${{ matrix.version }} - arch: ${{ matrix.arch }} - - name: Cache artifacts - uses: julia-actions/cache@v1 - - name: Build package - uses: julia-actions/julia-buildpkg@v1 - - name: Run tests - uses: julia-actions/julia-runtest@v1 - - name: Process coverage - uses: julia-actions/julia-processcoverage@v1 - - name: Upload coverage - uses: codecov/codecov-action@v3 - with: - files: lcov.info diff --git a/.github/workflows/ci_PR.yml b/.github/workflows/ci_PR.yml deleted file mode 100644 index f0e842af..00000000 --- a/.github/workflows/ci_PR.yml +++ /dev/null @@ -1,49 +0,0 @@ -name: CI-PR - -on: - pull_request: - paths-ignore: - - 'LICENSE' - - 'README.md' - -jobs: - test: - name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} - runs-on: ${{ matrix.os }} - if: "!contains(github.event.head_commit.message, 'skip ci')" - strategy: - fail-fast: false - matrix: - version: - - '1' - os: - - ubuntu-latest - - windows-latest - arch: - - x64 - include: - - version: '1.6' # test on oldest supported version - arch: x64 - os: ubuntu-latest - env: - JULIA_PKG_SERVER: '' - steps: - - name: Check out repository - uses: actions/checkout@v4 - - name: Set up Julia - uses: julia-actions/setup-julia@v1 - with: - version: ${{ matrix.version }} - arch: ${{ matrix.arch }} - - name: Cache artifacts - uses: julia-actions/cache@v1 - - name: Build package - uses: julia-actions/julia-buildpkg@v1 - - name: Run tests - uses: julia-actions/julia-runtest@v1 - - name: Process coverage - uses: julia-actions/julia-processcoverage@v1 - - name: Upload coverage - uses: codecov/codecov-action@v3 - with: - files: lcov.info diff --git a/Project.toml b/Project.toml index 625af647..703aba2c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,11 +1,12 @@ name = "ClosedLoopReachability" uuid = "73c0b437-a350-4e9b-97ac-9adb151c271b" -version = "0.3.0" +version = "0.4.0" [deps] CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" ControllerFormats = "02ac4b2c-022a-44aa-84a5-ea45a5754bcc" LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" +NeuralNetworkReachability = "5de1e908-2f08-45bf-a571-ac88a54f7e7f" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" ReachabilityAnalysis = "1e97bd63-91d1-579d-8e8d-501d2b57c93f" ReachabilityBase = "379f33d0-9447-4353-bd03-d664070e549f" @@ -15,10 +16,11 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df" [compat] CommonSolve = "0.1 - 0.2" ControllerFormats = "0.1 - 0.2" -LazySets = "2.7.6" +LazySets = "2.11.1" +NeuralNetworkReachability = "0.1" Parameters = "0.12" ReachabilityAnalysis = "0.14.5 - 0.22" -ReachabilityBase = "0.1.1 - 0.2" +ReachabilityBase = "0.2.1" Reexport = "0.2, 1" Requires = "0.5, 1" julia = "1.6" diff --git a/models/VertCAS/VertCAS.jl b/models/VertCAS/VertCAS.jl index 5b53a50d..936639ec 100644 --- a/models/VertCAS/VertCAS.jl +++ b/models/VertCAS/VertCAS.jl @@ -181,7 +181,7 @@ end function forward_adv(X::Singleton, τ, adv; alg=nothing) v = vcat(element(X), τ) v = normalize(v) - u = forward(CONTROLLERS[adv], v) + u = forward(v, CONTROLLERS[adv]) imax = argmax(u) return CTRL_IDX[imax] end @@ -191,7 +191,7 @@ function forward_adv(X::AbstractZonotope, τ, adv; alg=DeepZ()) Y = cartesian_product(X, Singleton([τ])) Y = normalize(Y) - out = forward(alg, CONTROLLERS[adv], Y) + out = forward(Y, CONTROLLERS[adv], alg) imax = argmax(high(out)) return CTRL_IDX[imax] diff --git a/src/ClosedLoopReachability.jl b/src/ClosedLoopReachability.jl index 4947c369..7810cd1f 100644 --- a/src/ClosedLoopReachability.jl +++ b/src/ClosedLoopReachability.jl @@ -2,7 +2,6 @@ module ClosedLoopReachability include("init.jl") include("problem.jl") -include("nnops.jl") include("setops.jl") include("split.jl") include("utils.jl") @@ -14,19 +13,16 @@ export ControlledPlant, BlackBoxController # splitters -export BoxSplitter, ZonotopeSplitter, +export BoxSplitter, + ZonotopeSplitter, IndexedSplitter, SignSplitter # solvers -export solve, forward, simulate, - DeepZ, SampledApprox, VertexSolver, BoxSolver, SplitSolver, BlackBoxSolver +export solve, simulate # utility functions export @modelpath, print_timed -# parsers -@reexport using ControllerFormats.FileFormats - end diff --git a/src/init.jl b/src/init.jl index bb7450e2..7d64d251 100644 --- a/src/init.jl +++ b/src/init.jl @@ -3,6 +3,11 @@ using Requires, Reexport using ControllerFormats # namespace conflict using ControllerFormats: Id +# parsers +@reexport using ControllerFormats.FileFormats + +@reexport using NeuralNetworkReachability.ForwardAlgorithms +using NeuralNetworkReachability.ForwardAlgorithms: ForwardAlgorithm @reexport using ReachabilityAnalysis # namespace conflict @@ -15,6 +20,8 @@ using ReachabilityAnalysis: _check_dim, _get_tspan, _get_cpost, _default_cpost, using ReachabilityBase.Timing, ReachabilityBase.Require +using ReachabilityBase.Arrays: SingleEntryVector +using ReachabilityBase.Comparison: isapproxzero using Parameters: @with_kw @@ -22,9 +29,6 @@ using Parameters: @with_kw const RA = ReachabilityAnalysis const IA = IntervalArithmetic -using LazySets: _leq, _geq, isapproxzero, _isapprox, array, - remove_zero_generators, remove_zero_columns, subtypes, SingleEntryVector - import CommonSolve: solve # optional dependencies diff --git a/src/nnops.jl b/src/nnops.jl deleted file mode 100644 index 9e72e1e6..00000000 --- a/src/nnops.jl +++ /dev/null @@ -1,399 +0,0 @@ -abstract type Solver end - -# ================================================ -# Internal forward network functions -# ================================================ - -# output of neural network for a single input -function forward(nnet::FeedforwardNetwork, x0::Vector{<:Number}) - x = x0 - @inbounds for layer in nnet.layers - W = layer.weights - b = layer.bias - x = layer.activation(W * x + b) - end - return x -end - -function forward(solver::Solver, nnet::FeedforwardNetwork, X0) - X = X0 - for layer in nnet.layers - _, X = forward_layer(solver, layer, X) - end - return X -end - -function forward_layer(solver, layer, reach) - Yl = forward_linear(solver, layer, reach) - Y = forward_act(solver, layer, Yl) - return Yl, Y -end - -# ================================================ -# Composite methods to compute the network output -# ================================================ - -@with_kw struct SplitSolver{S<:Solver,FS,FM} <: Solver - solver::S - split_fun::FS - merge_fun::FM -end - -function SplitSolver(solver) - # default: box approximation and split in two sets per dimension - split_fun = X -> split(box_approximation(X), 2 * ones(Int, dim(X))) - # default: box approximation of the union - merge_fun = X -> box_approximation(X) - return SplitSolver(solver, split_fun, merge_fun) -end - -function forward(solver::SplitSolver, nnet::FeedforwardNetwork, X0) - X0_split = solver.split_fun(X0) - Y_union = UnionSetArray() - for X in X0_split - Y = forward(solver.solver, nnet, X) - push!(array(Y_union), Y) - end - Y_merged = solver.merge_fun(Y_union) - return Y_merged -end - -# ============================================================================ -# Methods to approximate the network output (without mathematical guarantees) -# ============================================================================ - -# solver using the CH of the sampled outputs as an inner approx of the real output -@with_kw struct SampledApprox <: Solver - nsamples::Int = 10000 - include_vertices::Bool = true - directions = OctDirections -end - -function forward(solver::SampledApprox, nnet::FeedforwardNetwork, input) - samples = sample(input, solver.nsamples; - include_vertices=solver.include_vertices) - - m = dim_out(nnet) - if m == 1 - MIN = Inf - MAX = -Inf - for sample in samples - output = first(forward(nnet, sample)) - MIN = min(MIN, output) - MAX = max(MAX, output) - end - return LazySets.Interval(MIN, MAX) - else - vlist = Vector{Vector{eltype(samples[1])}}(undef, length(samples)) - @inbounds for (i, sample) in enumerate(samples) - vlist[i] = forward(nnet, sample) - end - convex_hull!(vlist) - P = VPolytope(vlist) - Z = overapproximate(P, Zonotope, solver.directions) - return Z - end -end - -# ============================================================================== -# Method to handle networks with ReLU, sigmoid, and Tanh activation functions -# from [FAS18] -# -# [FAS18]: Singh, Gagandeep, et al. "Fast and Effective Robustness -# Certification." NeurIPS 2018. -# ============================================================================== - -struct DeepZ <: Solver end - -function forward_linear(solver::DeepZ, L::DenseLayerOp, Z::AbstractZonotope) - return affine_map(L.weights, Z, L.bias) -end - -function forward_act(solver::DeepZ, L::DenseLayerOp{Id}, Z::AbstractZonotope) - return Z -end - -function forward_act(solver::DeepZ, L::DenseLayerOp{ReLU}, Z::AbstractZonotope) - return overapproximate(Rectification(Z), Zonotope) # implemented in LazySets -end - -function sigmoid(x::Number) - ex = exp(x) - return ex / (1 + ex) -end - -function sigmoid2(x::Number) - ex = exp(x) - return ex / (1 + ex)^2 -end - -function _overapproximate_zonotope(Z::AbstractZonotope{N}, act, act′) where {N} - c = copy(center(Z)) - G = copy(genmat(Z)) - n, m = size(G) - row_idx = Vector{Int}() - μ_idx = Vector{N}() - - @inbounds for i in 1:n - lx, ux = low(Z, i), high(Z, i) - ly, uy = act(lx), act(ux) - - if _isapprox(lx, ux) - c[i] = uy - for j in 1:m - G[i, j] = zero(N) - end - else - λ = min(act′(lx), act′(ux)) - μ₁ = (uy + ly - λ * (ux + lx)) / 2 - # Note: there is a typo in the paper (missing parentheses) - μ₂ = (uy - ly - λ * (ux - lx)) / 2 - c[i] = c[i] * λ + μ₁ - for j in 1:m - G[i, j] = G[i, j] * λ - end - push!(row_idx, i) - push!(μ_idx, μ₂) - end - end - - q = length(row_idx) - if q >= 1 - Gnew = zeros(N, n, q) - j = 1 - @inbounds for i in row_idx - Gnew[i, j] = μ_idx[j] - j += 1 - end - Gout = hcat(G, Gnew) - else - Gout = G - end - - return Zonotope(c, remove_zero_columns(Gout)) -end - -function forward_act(solver::DeepZ, L::DenseLayerOp{Sigmoid}, Z::AbstractZonotope) - act(x) = sigmoid(x) - act′(x) = sigmoid2(x) - return _overapproximate_zonotope(Z, act, act′) -end - -function forward_act(solver::DeepZ, L::DenseLayerOp{Tanh}, Z::AbstractZonotope) - act(x) = tanh(x) - act′(x) = 1 - tanh(x)^2 - return _overapproximate_zonotope(Z, act, act′) -end - -# ========================================================== -# Methods to handle networks with ReLU activation functions -# ========================================================== - -# solver that computes the box approximation in each layer -# it exploits that box(relu(X)) == relu(box(X)) -struct BoxSolver <: Solver end - -function forward(solver::BoxSolver, nnet::FeedforwardNetwork, X0) - X = X0 - for layer in nnet.layers - # affine map and box approximation - W = layer.weights - b = layer.bias - X_am = AffineMap(W, X, b) - - # activation function - if layer.activation isa Id - X = X_am - continue - end - - @assert layer.activation isa ReLU "unsupported activation function" - X = rectify(box_approximation(X_am)) - end - return X -end - -@with_kw struct ConcreteReLU <: Solver - concrete_intersection::Bool = false - convexify::Bool = false -end - -function forward(solver::ConcreteReLU, nnet::FeedforwardNetwork, X0) - X = [X0] - for layer in nnet.layers - if typeof(X[1]) <: UnionSetArray - X = [x.array for x in X] - X = reduce(vcat, X) - end - X = affine_map.(Ref(layer.weights), X, Ref(layer.bias)) - - # activation function - if layer.activation isa Id - continue - end - @assert layer.activation isa ReLU "unsupported activation function" - X = rectify.(X, solver.concrete_intersection) - end - return solver.convexify ? ConvexHullArray(X) : X -end - -# solver that propagates the vertices, computes their convex hull, and applies -# some postprocessing to the result -@with_kw struct VertexSolver{T} <: Solver - postprocessing::T = x -> x # default: identity (= no postprocessing) - apply_convex_hull::Bool = false -end - -function forward(solver::VertexSolver, nnet::FeedforwardNetwork, X0) - N = eltype(X0) - P = X0 - - for layer in nnet.layers - # apply affine map - W = layer.weights - b = layer.bias - P = convert(VPolytope, P) - Q_am = affine_map(W, P, b) - - # activation function - if layer.activation isa Id - P = Q_am - continue - end - - @assert layer.activation isa ReLU "unsupported activation function" - - # compute Q_chull = convex_hull(Q_am, rectify(Q_am)) - vlist = Vector{Vector{N}}() - vlist_rect = Vector{Vector{N}}() - for v in vertices(Q_am) - push!(vlist, v) - v_rect = rectify(v) - push!(vlist_rect, v_rect) - if v != v_rect - push!(vlist, v_rect) - end - end - if solver.apply_convex_hull || true - convex_hull!(vlist) - end - Q_chull = VPolytope(vlist) - - # filter out negative part - # Q_pos = box_approximation(VPolytope(vlist_rect)) # alternative - n = dim(Q_am) - Q_pos = HPolyhedron([HalfSpace(SingleEntryVector(i, n, -one(N)), zero(N)) for i in 1:n]) - P = intersection(Q_chull, Q_pos) - end - Q = solver.postprocessing(P) - return Q -end - -# ============================================================================== -# Methods to handle networks with sigmoid activation functions from [VER19] -# -# [VER19]: Ivanov, Radoslav, et al. "Verisig: verifying safety properties of -# hybrid systems with neural network controllers." Proceedings of the 22nd ACM -# International Conference on Hybrid Systems: Computation and Control. 2019. -# ============================================================================== - -# ref. Eq (6) in [VER19] -# d(σ(x))/dx = σ(x)*(1-σ(x)) -# g(t, x) = σ(tx) = 1 / (1 + exp(-tx)) -# dg(t, x)/dt = g'(t, x) = x * g(t, x) * (1 - g(t, x)) -@taylorize function sigmoid!(dx, x, p, t) - xᴶ, xᴾ = x - dx[1] = zero(xᴶ) - return dx[2] = xᴶ * (xᴾ - xᴾ^2) -end - -# footnote (3) in [VER19] -# d(tanh(x))/dx = 1 - tanh(x)^2 -# g(t, x) = tanh(tx) -# dg(t, x)/dt = g'(t, x) = x * (1 - g(t, x)^2) -@taylorize function tanh!(dx, x, p, t) - xᴶ, xᴾ = x - dx[1] = zero(xᴶ) - return dx[2] = xᴶ * (1 - xᴾ^2) -end - -const HALFINT = IA.Interval(0.5, 0.5) -const ZEROINT = IA.Interval(0.0, 0.0) -const ACTFUN = Dict(Tanh() => (tanh!, ZEROINT), - Sigmoid() => (sigmoid!, HALFINT)) - -# Method: Cartesian decomposition (intervals for each one-dimensional subspace) -# Only Tanh, Sigmoid and Id functions are supported -function forward(nnet::FeedforwardNetwork, X0::LazySet; - alg=TMJets(; abstol=1e-14, orderQ=2, orderT=6)) - - # initial states - xᴾ₀ = _decompose_1D(X0) - xᴾ₀ = array(xᴾ₀) # see https://github.com/JuliaReach/ReachabilityAnalysis.jl/issues/254 - xᴾ₀ = [x.dat for x in xᴾ₀] # use concrete interval matrix-vector operations - - for layer in nnet.layers # loop over layers - W = layer.weights - m, n = size(W) - b = layer.bias - act = layer.activation - - xᴶ′ = W * xᴾ₀ + b # (scalar matrix) * (interval vector) + (scalar vector) - - if act == Id() - xᴾ₀ = copy(xᴶ′) - continue - end - - activation!, ival = ACTFUN[act] - xᴾ′ = fill(ival, m) - - for i in 1:m # loop over coordinates - X0i = xᴶ′[i] × xᴾ′[i] - ivp = @ivp(x' = activation!(x), dim = 2, x(0) ∈ X0i) - sol = RA.solve(ivp; tspan=(0.0, 1.0), alg=alg) - - # interval overapproximation of the final reach-set along - # dimension 2, which corresponds to xᴾ - xᴾ_end = sol.F.ext[:xv][end][2] - xᴾ′[i] = xᴾ_end - end - xᴾ₀ = copy(xᴾ′) - end - return CartesianProductArray([Interval(x) for x in xᴾ₀]) -end - -# ===================================== -# Methods to handle arbitrary networks -# ===================================== - -struct BlackBoxSolver <: Solver end - -struct BlackBoxController{FT} <: AbstractNeuralNetwork - f::FT -end - -function forward(solver::BlackBoxSolver, bbc::BlackBoxController, X0) - return bbc.f(X0) -end - -function forward(bbc::BlackBoxController, X0) - return bbc.f(X0) -end - -# ============================= -# Handling of singleton inputs -# ============================= - -function forward(nnet::FeedforwardNetwork, X0::AbstractSingleton) - x0 = element(X0) - x1 = forward(nnet, x0) - return Singleton(x1) -end - -for SOLVER in subtypes(Solver, true) - @eval function forward(solver::$SOLVER, nnet::FeedforwardNetwork, - X0::AbstractSingleton) - return forward(nnet, X0) - end -end diff --git a/src/setops.jl b/src/setops.jl index a3163c5d..6ee3fbd5 100644 --- a/src/setops.jl +++ b/src/setops.jl @@ -11,23 +11,6 @@ function _project_oa(R::AbstractTaylorModelReachSet, vars, t; remove_zero_genera return LazySets.project(set(Z), vars; remove_zero_generators=remove_zero_generators) end -# ======================================== -# Decomposition operations -# ======================================== - -# decompose a set into the Cartesian product of intervals - -function _decompose_1D(X0::LazySet{N}) where {N} - n = dim(X0) - out = Vector{Interval{Float64,IA.Interval{Float64}}}(undef, n) - - @inbounds for i in 1:n - eᵢ = SingleEntryVector(i, n, one(N)) - out[i] = LazySets.Interval(-ρ(-eᵢ, X0), ρ(eᵢ, X0)) - end - return CartesianProductArray(out) -end - # ============================================================== # Construct the initial states for the continuous post-operator # ============================================================== @@ -38,7 +21,7 @@ abstract type AbstractReconstructionMethod end struct CartesianProductReconstructor <: AbstractReconstructionMethod end -function _reconstruct(method::CartesianProductReconstructor, P₀::LazySet, U₀::LazySet, R, ti) +function _reconstruct(::CartesianProductReconstructor, P₀::LazySet, U₀::LazySet, R, ti) Q₀ = P₀ × U₀ return Q₀ end @@ -48,7 +31,7 @@ end end # if no Taylor model is available => use the given set P₀ -function _reconstruct(method::TaylorModelReconstructor, P₀::LazySet, U₀, R::Nothing, ti) +function _reconstruct(::TaylorModelReconstructor, P₀::LazySet, U₀, R::Nothing, ti) return _reconstruct(CartesianProductReconstructor(), P₀, U₀, R, ti) end diff --git a/src/simulate.jl b/src/simulate.jl index ade3b2fe..37867309 100644 --- a/src/simulate.jl +++ b/src/simulate.jl @@ -56,7 +56,7 @@ function simulate(cp::AbstractControlProblem, args...; kwargs...) for j in 1:trajectories x₀ = x0_vec[j] x₀ = apply(preprocessing, x₀) - network_output = forward(network, x₀) + network_output = forward(x₀, network) control_signals[j] = apply(postprocessing, network_output) end all_controls[i] = control_signals diff --git a/src/solve.jl b/src/solve.jl index 861e66cd..5eab5a9e 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -25,7 +25,7 @@ The control signals are stored in the `ext` field with each flowpipe. - Use the `T` keyword argument to specify the time horizon; the initial time is then assumed to be zero. -- Use the `alg_nn` keyword argument to specify the solver for the controller. +- Use the `alg_nn` keyword argument to specify the algorithm for the controller. - While this function is written with a neural-network controlled systems in mind, the type of the controller is arbitrary, as long as a function @@ -53,7 +53,7 @@ function solve(prob::AbstractControlProblem, args...; kwargs...) cpost = _default_cpost(ivp, tspan; kwargs...) end - solver = _get_alg_nn(args...; kwargs...) + nnpost = _get_alg_nn(args...; kwargs...) splitter = get(kwargs, :splitter, NoSplitter()) @@ -63,21 +63,21 @@ function solve(prob::AbstractControlProblem, args...; kwargs...) remove_zero_generators = get(kwargs, :remove_zero_generators, true) - sol = _solve(prob, cpost, solver, tvec, τ, splitter, input_splitter, + sol = _solve(prob, cpost, nnpost, tvec, τ, splitter, input_splitter, rec_method, remove_zero_generators) - d = Dict{Symbol,Any}(:solver => solver) + d = Dict{Symbol,Any}(:solver => nnpost) return ReachSolution(sol, cpost, d) end function _get_alg_nn(args...; kwargs...) if haskey(kwargs, :alg_nn) - solver = kwargs[:alg_nn] + nnpost = kwargs[:alg_nn] else - throw(ArgumentError("the solver for the controller `alg_nn` should be " * - "specified, but was not found")) + throw(ArgumentError("the algorithm for the controller `alg_nn` " * + "should be specified, but was not found")) end - return solver + return nnpost end # element of the waiting list: a flowpipe with corresponding iteration @@ -88,7 +88,7 @@ end function _solve(cp::ControlledPlant, cpost::AbstractContinuousPost, - solver::Solver, + nnpost::ForwardAlgorithm, tvec::AbstractVector, sampling_time::N, splitter::AbstractSplitter, @@ -135,13 +135,13 @@ function _solve(cp::ControlledPlant, # first perform an isolated analysis because of problems in TaylorSeries # (global variables need to be written once) @inbounds results[1] = _solve_one(R, first(X₀s), W₀, S, st_vars, t0, t1, - cpost, rec_method, solver, network, preprocessing, + cpost, rec_method, nnpost, network, preprocessing, postprocessing, input_splitter) # parallelize analysis of the remaining parts Threads.@threads for i in 2:length(results) @inbounds results[i] = _solve_one(R, X₀s[i], W₀, S, st_vars, t0, t1, - cpost, rec_method, solver, network, preprocessing, + cpost, rec_method, nnpost, network, preprocessing, postprocessing, input_splitter) end @@ -170,7 +170,7 @@ function _solve(cp::ControlledPlant, # parallelize analysis Threads.@threads for i in 1:length(results) @inbounds results[i] = _solve_one(R, X₀s[i], W₀, S, st_vars, t0, t1, - cpost, rec_method, solver, network, preprocessing, + cpost, rec_method, nnpost, network, preprocessing, postprocessing, input_splitter) end # collect results from all threads @@ -187,9 +187,9 @@ function _solve(cp::ControlledPlant, return MixedFlowpipe(flowpipes) end -function nnet_forward(solver, network, X, preprocessing, postprocessing) +function nnet_forward(nnpost, network, X, preprocessing, postprocessing) X = apply(preprocessing, X) - U = forward(solver, network, X) + U = forward(X, network, nnpost) U = apply(postprocessing, U) if dim(U) == 1 # simplify the control input for intervals U = overapproximate(U, Interval) @@ -197,13 +197,13 @@ function nnet_forward(solver, network, X, preprocessing, postprocessing) return U end -function _solve_one(R, X₀, W₀, S, st_vars, t0, t1, cpost, rec_method, solver, +function _solve_one(R, X₀, W₀, S, st_vars, t0, t1, cpost, rec_method, nnpost, network, preprocessing, postprocessing, splitter) # add disturbances (if any) P₀ = isnothing(W₀) ? X₀ : X₀ × W₀ # get new control inputs from the controller - U = nnet_forward(solver, network, X₀, preprocessing, postprocessing) + U = nnet_forward(nnpost, network, X₀, preprocessing, postprocessing) dt = t0 .. t1 diff --git a/src/utils.jl b/src/utils.jl index 3bee863b..3b84c770 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -132,11 +132,3 @@ function random_control_problems(n::Integer, m::Integer; ivp=nothing, period=0.1 return problems end - -# relative size between the set-based output and the (CH of) sampled output -function relative_size(X0, nsamples, controller, solver=DeepZ()) - @assert output_dim(controller) == 1 "the dimension of the output of the network needs to be 1, but is $output_dim(controller)" - o = overapproximate(forward(solver, controller, X0), Interval) - u = forward(SampledApprox(nsamples), controller, X0) - return diameter(o) / diameter(u) -end diff --git a/test/Project.toml b/test/Project.toml index 46201237..0c363327 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,10 +1,2 @@ [deps] -MAT = "23992714-dd62-5051-b70f-ba57cb901cac" -ONNX = "d0dd6a25-fac6-55c0-abf7-829e0c774d20" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" - -[compat] -MAT = "0.10" -ONNX = "0.2" -YAML = "0.4" diff --git a/test/black_box_toy_model.jl b/test/black_box_toy_model.jl deleted file mode 100644 index c3e43652..00000000 --- a/test/black_box_toy_model.jl +++ /dev/null @@ -1,53 +0,0 @@ -using ClosedLoopReachability - -@taylorize function sys!(dx, x, p, t) - dx[1] = x[2] - dx[2] = zero(x[2]) - return dx -end; - -function oscillator(X::LazySet) - if high(X)[1] > 1 - if low(X)[1] < 1 - res = Interval(-0.1, 0.1) - else - res = Interval(-0.1, -0.1) - end - else - res = Interval(0.1, 0.1) - end - return res -end; - -function oscillator(x::Vector) - if x[1] > 1 - res = [-0.1] - else - res = [0.1] - end - return res -end; -controller = BlackBoxController(oscillator); - -X₀ = Interval(0.95, 0.99); -U = ZeroSet(1); -vars_idx = Dict(:states => [1], :controls => 2); -ivp = @ivp(x' = sys!(x), dim:2, x(0) ∈ X₀ × U); -period = 0.1; -prob = ControlledPlant(ivp, controller, vars_idx, period); -alg = TMJets21b(; abstol=1e-10, orderT=8, orderQ=2, adaptive=true); -alg_nn = BlackBoxSolver(); -N = 10; -T = N * period; - -sol_raw = solve(prob; T=T, alg_nn=alg_nn, alg=alg); - -# import DifferentialEquations -# sim = simulate(prob, T=T; trajectories=10, include_vertices=true); - -# using Plots -# fig = plot() -# vars = (0, 1) -# plot!(fig, sol_raw, vars=vars, c=:yellow, lw=2.0, ls=:dash) -# plot_simulation!(fig, sim; vars=vars, color=:red, lab=""); -# savefig("oscillator_$(N)_steps") diff --git a/test/runtests.jl b/test/runtests.jl index 61eaa6e8..f5ce9fd8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1 @@ using Test, ClosedLoopReachability - -@testset "Toy model (black-box network)" begin - include("black_box_toy_model.jl") -end