Merge pull request #3 from blegat/refact-stopcrit

Refact Stopping Criterion

src/StochasticDualDynamicProgramming.jl

@@ -13,6 +13,7 @@ import Base.show, Base.isless
 
 # Utils
 include("mycomp.jl")
+include("stats.jl")
 
 # Abstract components
 # Cut Manager

src/sddp.jl

@@ -5,64 +5,6 @@ type SDDPSolution
     attrs
 end
 
-macro mytime(x)
-    quote
-        y = @timed $(esc(x))
-        # y[1] is returned value
-        # y[2] is time in seconds
-        y[2]
-    end
-end
-
-type SDDPStats
-    nsolved::Int
-    solvertime::Float64
-    nmerged::Int
-    mergetime::Float64
-    nfcuts::Int
-    fcutstime::Float64
-    nocuts::Int
-    ocutstime::Float64
-    nsetx::Int
-    setxtime::Float64
-end
-
-SDDPStats() = SDDPStats(0,.0,0,.0,0,.0,0,.0,0,.0)
-
-import Base: +, show
-
-function +(a::SDDPStats, b::SDDPStats)
-    SDDPStats(a.nsolved + b.nsolved, a.solvertime + b.solvertime,
-    a.nmerged + b.nmerged, a.mergetime  + b.mergetime,
-    a.nfcuts  + b.nfcuts , a.fcutstime  + b.fcutstime,
-    a.nocuts  + b.nocuts , a.ocutstime  + b.ocutstime,
-    a.nsetx   + b.nsetx  , a.setxtime   + b.setxtime)
-end
-
-function showtime(t::Float64)
-    if !isfinite(t)
-        "   min    s    ms    μs"
-    else
-        s = Int(floor(t))
-        t = (t - s) * 1000
-        min = div(s, 60)
-        s = mod(s, 60)
-        ms = Int(floor(t))
-        t = (t - ms) * 1000
-        μs = Int(floor(t))
-        @sprintf "%3dmin %3ds %3dms %3dμs" min s ms μs
-    end
-end
-
-function Base.show(io::IO, stat::SDDPStats)
-    println("                        |     Total time [s]      |  Number  | Average time [s]")
-    @printf "        Solving problem | %s | %8d | %s\n" showtime(stat.solvertime) stat.nsolved showtime(stat.solvertime / stat.nsolved)
-    @printf "          Merging paths | %s | %8d | %s\n" showtime(stat.mergetime ) stat.nmerged showtime(stat.mergetime  / stat.nmerged)
-    @printf "Adding feasibility cuts | %s | %8d | %s\n" showtime(stat.fcutstime ) stat.nfcuts  showtime(stat.fcutstime  / stat.nfcuts )
-    @printf "Adding  optimality cuts | %s | %8d | %s\n" showtime(stat.ocutstime ) stat.nocuts  showtime(stat.ocutstime  / stat.nocuts )
-    @printf "Setting parent solution | %s | %8d | %s\n" showtime(stat.setxtime  ) stat.nsetx   showtime(stat.setxtime   / stat.nsetx  )
-    @printf "                        | %s |" showtime(stat.solvertime+stat.fcutstime+stat.ocutstime+stat.setxtime)
-end
 
 function meanstdpaths(z::Vector{Float64}, proba::Vector{Float64}, npaths::Vector{Int}, Ktot)
     if Ktot != -1
@@ -176,6 +118,7 @@ function iteration{S}(root::SDDPNode{S}, Ktot::Int, num_stages, verbose, pathsel
 
     stats.solvertime += @mytime rootsol = loadAndSolve(root)
     stats.nsolved += 1
+    stats.niterations += 1
     infeasibility_detected = rootsol.status == :Infeasible
     if infeasibility_detected
         pathsd = Dict{SDDPNode, Vector{SDDPPath}}()
@@ -344,7 +287,14 @@ function iteration{S}(root::SDDPNode{S}, Ktot::Int, num_stages, verbose, pathsel
         end
         z_UB, σ = meanstdpaths(endedpaths, Ktot)
     end
-    rootsol, stats, z_UB, σ
+
+    # update stats
+    stats.upperbound = z_UB
+    stats.σ_UB = σ
+    stats.npaths = Ktot
+    stats.lowerbound = rootsol.objval
+
+    rootsol, stats
 end
 
 """
@@ -361,7 +311,6 @@ function SDDP(root::SDDPNode, num_stages; K::Int=25, stopcrit::AbstractStoppingC
         error("Invalid pathsel")
     end
     rootsol = nothing
-    totaltime = 0
     totalstats = SDDPStats()
 
     z_UB = Inf
@@ -371,17 +320,18 @@ function SDDP(root::SDDPNode, num_stages; K::Int=25, stopcrit::AbstractStoppingC
     nfcuts = 0
     nocuts = 0
 
-    while (rootsol === nothing || rootsol.status != :Infeasible) && !stop(stopcrit, iter, nfcuts, nocuts, K, z_LB, z_UB, σ)
-        itertime = @mytime rootsol, stats, z_UB, σ = iteration(root, K, num_stages, verbose, pathsel, ztol)
+    while (rootsol === nothing || rootsol.status != :Infeasible) && !stop(stopcrit, totalstats)
+        itertime = @mytime rootsol, stats = iteration(root, K, num_stages, verbose, pathsel, ztol)
         z_LB = rootsol.objval
         iter += 1
         nfcuts = stats.nfcuts
         nocuts = stats.nocuts
 
-        totaltime  += itertime
+        totalstats.time += itertime
+
         totalstats += stats
         if verbose >= 2
-            println("Iteration $iter completed in $itertime s (Total time is $totaltime)")
+            println("Iteration $iter completed in $itertime s (Total time is $(stats.time))")
             println("Status: $(rootsol.status)")
             println("z_UB: $(z_UB)")
             println("z_LB: $(z_LB)")

src/stats.jl

@@ -0,0 +1,88 @@
+export AbstractSDDPStats
+
+abstract AbstractSDDPStats
+
+type SDDPStats <: AbstractSDDPStats
+    # number of calls to solver
+    nsolved::Int
+    # total time passed inside solver
+    solvertime::Float64
+    # number of cuts merged
+    nmerged::Int
+    # time required to merge
+    mergetime::Float64
+    # number of feasibility cuts
+    nfcuts::Int
+    fcutstime::Float64
+    # number of optimality cuts
+    nocuts::Int
+    ocutstime::Float64
+    nsetx::Int
+    setxtime::Float64
+
+    # number of iterations
+    niterations::Int
+    # n forwards passes of last computation of upper-bound:
+    npaths::Int
+    # current lower bound
+    lowerbound::Float64
+    # current lower bound
+    upperbound::Float64
+    # upper-bound std:
+    σ_UB::Float64
+    # total time
+    time::Float64
+end
+
+SDDPStats() = SDDPStats(0,.0,0,.0,0,.0,0,.0,0,.0, 0, 0, 0, 0, 0, 0)
+
+import Base: +, show
+
+# Add two `SDDPStats` objects and return a new `SDDPStats`
+# The second SDDPStats is supposed to be more up to date that the first one
+# WARNING: the addition is currently non commutative !
+function +(a::SDDPStats, b::SDDPStats)
+    SDDPStats(a.nsolved + b.nsolved, a.solvertime + b.solvertime,
+    a.nmerged + b.nmerged, a.mergetime  + b.mergetime,
+    a.nfcuts  + b.nfcuts , a.fcutstime  + b.fcutstime,
+    a.nocuts  + b.nocuts , a.ocutstime  + b.ocutstime,
+    a.nsetx   + b.nsetx  , a.setxtime   + b.setxtime,
+    a.niterations   + b.niterations  , b.npaths,
+    b.lowerbound, b.upperbound,
+    b.σ_UB, a.time + b.time)
+end
+
+macro mytime(x)
+    quote
+        y = @timed $(esc(x))
+        # y[1] is returned value
+        # y[2] is time in seconds
+        y[2]
+    end
+end
+
+
+function showtime(t::Float64)
+    if !isfinite(t)
+        "   min    s    ms    μs"
+    else
+        s = Int(floor(t))
+        t = (t - s) * 1000
+        min = div(s, 60)
+        s = mod(s, 60)
+        ms = Int(floor(t))
+        t = (t - ms) * 1000
+        μs = Int(floor(t))
+        @sprintf "%3dmin %3ds %3dms %3dμs" min s ms μs
+    end
+end
+
+function Base.show(io::IO, stat::SDDPStats)
+    println("                        |     Total time [s]      |  Number  | Average time [s]")
+    @printf "        Solving problem | %s | %8d | %s\n" showtime(stat.solvertime) stat.nsolved showtime(stat.solvertime / stat.nsolved)
+    @printf "          Merging paths | %s | %8d | %s\n" showtime(stat.mergetime ) stat.nmerged showtime(stat.mergetime  / stat.nmerged)
+    @printf "Adding feasibility cuts | %s | %8d | %s\n" showtime(stat.fcutstime ) stat.nfcuts  showtime(stat.fcutstime  / stat.nfcuts )
+    @printf "Adding  optimality cuts | %s | %8d | %s\n" showtime(stat.ocutstime ) stat.nocuts  showtime(stat.ocutstime  / stat.nocuts )
+    @printf "Setting parent solution | %s | %8d | %s\n" showtime(stat.setxtime  ) stat.nsetx   showtime(stat.setxtime   / stat.nsetx  )
+    @printf "                        | %s |" showtime(stat.solvertime+stat.fcutstime+stat.ocutstime+stat.setxtime)
+end

src/stopcrit.jl

@@ -1,5 +1,5 @@
 import Base.|, Base.&
-export stop, AbstractStoppingCriterion, OrStoppingCriterion, AndStoppingCriterion, IterLimit, Pereira, CutLimit
+export stop, AbstractStoppingCriterion, OrStoppingCriterion, AndStoppingCriterion, IterLimit, Pereira, CutLimit, TimeLimit
 
 abstract AbstractStoppingCriterion
 
@@ -11,7 +11,7 @@ If `iter` is 0, no iteration has already been done, otherwise, the `iter`th iter
 This iteration used `K` paths and generated `nfcuts` (resp. `nocuts`) new feasibility (resp. optimality) cuts.
 The lower bound is now `z_LB` and the upper bound has mean `z_UB` and variance `σ`.
 """
-function stop(s::AbstractStoppingCriterion, iter, nfcuts, nocuts, K, z_LB, z_UB, σ)
+function stop(s::AbstractStoppingCriterion, stats::AbstractSDDPStats)
     error("`stop' function not defined for $(typeof(s))")
 end
 
@@ -25,8 +25,8 @@ type OrStoppingCriterion <: AbstractStoppingCriterion
     rhs::AbstractStoppingCriterion
 end
 
-function stop(s::OrStoppingCriterion, iter, nfcuts, nocuts, K, z_LB, z_UB, σ)
-    stop(s.lhs, iter, nfcuts, nocuts, K, z_LB, z_UB, σ) || stop(s.rhs, iter, nfcuts, nocuts, K, z_LB, z_UB, σ)
+function stop(s::OrStoppingCriterion, stats::AbstractSDDPStats)
+    stop(s.lhs, stats) || stop(s.rhs, stats)
 end
 
 function (|)(lhs::AbstractStoppingCriterion, rhs::AbstractStoppingCriterion)
@@ -43,8 +43,8 @@ type AndStoppingCriterion <: AbstractStoppingCriterion
     rhs::AbstractStoppingCriterion
 end
 
-function stop(s::AndStoppingCriterion, iter, nfcuts, nocuts, K, z_LB, z_UB, σ)
-    stop(s.lhs, iter, nfcuts, nocuts, K, z_LB, z_UB, σ) && stop(s.rhs, iter, nfcuts, nocuts, K, z_LB, z_UB, σ)
+function stop(s::AndStoppingCriterion, stats::AbstractSDDPStats)
+    stop(s.lhs, stats) && stop(s.rhs, stats)
 end
 
 function (&)(lhs::AbstractStoppingCriterion, rhs::AbstractStoppingCriterion)
@@ -60,8 +60,8 @@ type IterLimit <: AbstractStoppingCriterion
     limit::Int
 end
 
-function stop(s::IterLimit, iter, nfcuts, nocuts, K, z_LB, z_UB, σ)
-    iter >= s.limit
+function stop(s::IterLimit, stats::AbstractSDDPStats)
+    stats.niterations >= s.limit
 end
 
 """
@@ -74,10 +74,26 @@ type CutLimit <: AbstractStoppingCriterion
     limit::Int
 end
 
-function stop(s::CutLimit, iter, nfcuts, nocuts, K, z_LB, z_UB, σ)
-    iter > 0 && nfcuts + nocuts <= s.limit
+function stop(s::CutLimit, stats::AbstractSDDPStats)
+    stats.niterations > 0 && stats.nfcuts + stats.nocuts <= s.limit
 end
 
+
+"""
+$(TYPEDEF)
+
+Stops if total time of execution is greater than the time limit specified.
+For instance, `TimeLimit(100)` stops after 100s.
+"""
+type TimeLimit <: AbstractStoppingCriterion
+    timelimit::Float64
+end
+
+function stop(s::TimeLimit, stats::AbstractSDDPStats)
+    stats.niterations > 0 && stats.time > s.timelimit
+end
+
+
 """
 $(TYPEDEF)
 
@@ -90,8 +106,13 @@ type Pereira <: AbstractStoppingCriterion
     Pereira(α=2.0, β=0.05) = new(α, β)
 end
 
-function stop(s::Pereira, iter, nfcuts, nocuts, K, z_LB, z_UB, σ)
-    if iter > 0
+function stop(s::Pereira, stats::AbstractSDDPStats)
+    z_UB = stats.upperbound
+    z_LB = stats.lowerbound
+    K = stats.npaths
+    σ = stats.σ_UB
+
+    if stats.niterations > 0
         @assert K >= 0
         σ1 = σ / √K
         σ2 = s.α * σ1

test/stopcrit.jl

@@ -5,7 +5,21 @@ K = 42
 z_LB = 1
 z_UB = 1
 σ = 1
-@test_throws ErrorException stop(InvalidStoppingCriterion(), 1, 2, 3, K, z_LB, z_UB, σ)
-@test stop(AndStoppingCriterion(IterLimit(8), CutLimit(8)), 8, 2, 6, K, z_LB, z_UB, σ)
-@test !stop(AndStoppingCriterion(IterLimit(8), CutLimit(8)), 7, 2, 6, K, z_LB, z_UB, σ)
-@test !stop(AndStoppingCriterion(IterLimit(8), CutLimit(8)), 7, 2, 5, K, z_LB, z_UB, σ)
+stats = StochasticDualDynamicProgramming.SDDPStats()
+
+stats.upperbound = z_UB
+stats.niterations = 8
+stats.nocuts = 8
+stats.lowerbound = z_LB
+stats.npaths = K
+stats.σ_UB = σ
+
+@test_throws ErrorException stop(InvalidStoppingCriterion(), stats)
+@test stop(AndStoppingCriterion(IterLimit(8), CutLimit(8)), stats)
+stats.niterations = 7
+stats.nocuts = 2
+stats.nfcuts = 6
+@test !stop(AndStoppingCriterion(IterLimit(8), CutLimit(8)), stats)
+stats.nocuts = 2
+stats.nfcuts = 5
+@test !stop(AndStoppingCriterion(IterLimit(8), CutLimit(8)), stats)