Merge pull request #120 from sisl/bigM

automate calculating M

src/optimization/nsVerify.jl

@@ -25,23 +25,46 @@ Sound and complete.
 """
 @with_kw struct NSVerify <: Solver
     optimizer = GLPK.Optimizer
-    m::Float64 = 1000.0  # The big M in the linearization
+    m = nothing
 end
 
 function solve(solver::NSVerify, problem::Problem)
+    # Set M automatically if not already set.
+    if isnothing(solver.m)
+        M = set_automatic_M(problem)
+    else
+        @warn "M should be chosen carefully. An M which is too small will cause
+        the problem to be solved incorrectly. Not setting the `m` keyword, or setting
+        it equal to `nothing` will cause a safe value to be calculated automatically.
+        E.g. NSVerify()." maxlog = 1
+        M = solver.m
+    end
+
+    # Set up and solve the problem
     model = Model(solver)
     neurons = init_neurons(model, problem.network)
     deltas = init_deltas(model, problem.network)
     add_set_constraint!(model, problem.input, first(neurons))
     add_complementary_set_constraint!(model, problem.output, last(neurons))
-    encode_network!(model, problem.network, neurons, deltas, MixedIntegerLP(solver.m))
+    encode_network!(model, problem.network, neurons, deltas, MixedIntegerLP(M))
     feasibility_problem!(model)
     optimize!(model)
+
+
     if termination_status(model) == OPTIMAL
         return CounterExampleResult(:violated, value.(first(neurons)))
-    end
-    if termination_status(model) == INFEASIBLE
+
+    elseif termination_status(model) == INFEASIBLE
         return CounterExampleResult(:holds)
+    else
+        return CounterExampleResult(:unknown)
     end
-    return CounterExampleResult(:unknown)
+end
+
+function set_automatic_M(problem)
+    # Compute the largest pre-activation bound absolute value.
+    # M must be larger than any value a variable in the problem
+    # can take.
+    bounds = get_bounds(problem, false)
+    M = maximum(abs, Iterators.flatten(abs.(hr.center) + hr.radius for hr in bounds))
 end

src/optimization/utils/constraints.jl

@@ -184,8 +184,8 @@ function encode_layer!(MIP::MixedIntegerLP,
         @constraints(model, begin
                                 z_next[j] >= ẑ[j]
                                 z_next[j] >= 0.0
-                                z_next[j] <= ẑ[j] + m * δ[j]
-                                z_next[j] <= m - m * δ[j]
+                                z_next[j] <= ẑ[j] + m * (1 - δ[j])
+                                z_next[j] <= m * δ[j]
                             end)
     end
 end

src/utils/util.jl

@@ -331,27 +331,30 @@ end
 
 """
     get_bounds(problem::Problem)
-    get_bounds(nnet::Network, input::Hyperrectangle)
+    get_bounds(nnet::Network, input::Hyperrectangle, [true])
 
-This function calls maxSens to compute node-wise bounds given a input set.
+Computes node-wise bounds given a input set. The optional last 
+argument determines whether the bounds are pre- or post-activation.
 
 Return:
-- `Vector{Hyperrectangle}`: bounds for all nodes **after** activation. `bounds[1]` is the input set.
+- `Vector{Hyperrectangle}`: bounds for all nodes. `bounds[1]` is the input set.
 """
 function get_bounds(nnet::Network, input::Hyperrectangle, act::Bool = true) # NOTE there is another function by the same name in convDual. Should reconsider dispatch
-    if act
-        solver = MaxSens(0.0, true)
-        bounds = Vector{Hyperrectangle}(undef, length(nnet.layers) + 1)
-        bounds[1] = input
-        for (i, layer) in enumerate(nnet.layers)
-            bounds[i+1] = forward_layer(solver, layer, bounds[i])
+    bounds = Vector{Hyperrectangle}(undef, length(nnet.layers) + 1)
+    bounds[1] = input
+    b = input
+    for (i, layer) in enumerate(nnet.layers)
+        if act
+            b = approximate_affine_map(layer, bounds[i])
+            bounds[i+1] = approximate_act_map(layer, b)
+        else
+            bounds[i+1] = approximate_affine_map(layer, b)
+            b = approximate_act_map(layer, bounds[i+1])
         end
-        return bounds
-    else
-       error("before activation bounds not supported yet.")
     end
+    return bounds
 end
-get_bounds(problem::Problem) = get_bounds(problem.network, problem.input)
+get_bounds(problem::Problem, args...) = get_bounds(problem.network, problem.input, args...)
 
 """
     affine_map(layer, x)
@@ -375,6 +378,23 @@ function approximate_affine_map(layer::Layer, input::Hyperrectangle)
     return Hyperrectangle(c, r)
 end
 
+"""
+   approximate_act_map(layer, input::Hyperrectangle)
+
+Returns a Hyperrectangle overapproximation of the activation map of the input.
+`act`must be monotonic.
+"""
+function approximate_act_map(act::ActivationFunction, input::Hyperrectangle)
+    β    = act.(input.center)
+    βmax = act.(high(input))
+    βmin = act.(low(input))
+    c    = (βmax + βmin)/2
+    r    = (βmax - βmin)/2
+    return Hyperrectangle(c, r)
+end
+
+approximate_act_map(layer::Layer, input::Hyperrectangle) = approximate_act_map(layer.activation, input)
+
 
 """
     split_interval(dom, i)