Merge 225df98 into 8d64882

src/StructDualDynProg.jl

@@ -20,6 +20,8 @@ include("stats.jl")
 # Abstract components
 # Stopping Criterion
 include("stopcrit.jl")
+# Cut Generator
+include("cutgen.jl")
 # NLDS Model
 include("cutstore.jl")
 include("solver.jl")

src/comp.jl

@@ -4,16 +4,15 @@ _isapproxzero{T<:Real}(x::T, ztol) = x == zero(T)
 _isapproxzero{T<:AbstractFloat}(x::T, ztol) = -ztol < x < ztol
 _isapproxzero{T<:Real}(x::Vector{T}, ztol) = _isapproxzero(sum(abs, x), ztol)
 
-_isapprox{T<:Real}(x::T, y::T, ztol) = x == y
+_isapprox(x::Real, y::Real, ztol) = x == y
 # I check with zero because isapprox(0, 1e-100) is false...
 # but isapprox(1e-100, 2e-100) should be false
-_isapprox{T<:AbstractFloat}(x::T, y::T, ztol) = (x == zero(T) ? _isapproxzero(y, ztol) : (y == zero(T) ? _isapproxzero(x, ztol) : isapprox(x, y)))
-_isapprox{T<:Real, S<:Real}(x::T, y::S) = _isapprox(promote(x, y)...)
-_isapprox{T<:AbstractFloat}(x::Vector{T}, y::Vector{T}, ztol) = (x == zero(x) ? _isapproxzero(y, ztol) : (y == zero(y) ? _isapproxzero(x, ztol) : isapprox(x, y)))
+_isapprox(x::AbstractFloat, y::AbstractFloat, ztol) = (_isapproxzero(x, ztol) ? _isapproxzero(y, ztol) : (_isapproxzero(y, ztol) ? false : isapprox(x, y)))
+_isapprox{S<:AbstractFloat, T<:AbstractFloat}(x::Vector{S}, y::Vector{T}, ztol) = (x == zero(x) ? _isapproxzero(y, ztol) : (y == zero(y) ? _isapproxzero(x, ztol) : isapprox(x, y)))
 
 _lt{T<:Real}(x::T, y::T, ztol) = x < y
-_lt{T<:AbstractFloat}(x::T, y::T, ztol) = x < y && !_isapprox(x, y, ztol)
 _lt{S<:Real,T<:Real}(x::S, y::T, ztol) = _lt(promote(x, y)..., ztol)
+_lt(x::AbstractFloat, y::AbstractFloat, ztol) = x < y && !_isapprox(x, y, ztol)
 # _gt{S<:Real, T<:Real}(x::S, y::T) = _lt(y, x)
 # _leq{T<:Real}(x::T, y::T) = x <= y
 # _leq{T<:AbstractFloat}(x::T, y::T) = x <= y || _isapprox(x, y)

src/cutgen.jl

@@ -0,0 +1,79 @@
+export AbstractCutGenerator, NoOptimalityCutGenerator, MultiCutGenerator, AvgCutGenerator
+export nθ, solveallchildren, needallchildsol
+@compat abstract type AbstractCutGenerator end
+
+# No Optimality Cut Generator
+immutable NoOptimalityCutGenerator <: AbstractCutGenerator
+end
+nθ(::NoOptimalityCutGenerator, proba) = 0
+needallchildsol(::NoOptimalityCutGenerator) = false
+function gencut(::NoOptimalityCutGenerator, parent, path, stats, ztol)
+end
+function applycut(::NoOptimalityCutGenerator, node, g)
+    applyoptimalitycut!(node)
+end
+
+# Multi Cut Generator
+immutable MultiCutGenerator <: AbstractCutGenerator
+end
+nθ(::MultiCutGenerator, proba) = length(proba)
+needallchildsol(::MultiCutGenerator) = false
+function gencut(::MultiCutGenerator, parent, path, stats, ztol)
+    for i in 1:length(parent.children)
+        if !isnull(path.childsols[i]) && get(path.childsols[i]).status != :Unbounded
+            childsol = get(path.childsols[i])
+            a = childsol.πT
+            β = childsol.πh + childsol.σd
+            @assert childsol.status == :Optimal
+            if isnull(parent.childT)
+                aT = a
+            else
+                aT = get(parent.childT)[i]' * a
+            end
+            β += childsol.ρe
+            if path.sol.status == :Unbounded || _lt(path.sol.θ[i], β - dot(aT, path.sol.x), ztol)
+                stats.ocutstime += @_time pushoptimalitycutforparent!(parent.children[i], a, β, parent)
+                stats.nocuts += 1
+            end
+        end
+    end
+end
+function applycut(::MultiCutGenerator, node, g)
+    for child in children(g, node)
+        applyoptimalitycutforparent!(child)
+    end
+end
+
+# Average Cut Generator
+immutable AvgCutGenerator <: AbstractCutGenerator
+end
+nθ(::AvgCutGenerator, proba) = 1
+needallchildsol(::AvgCutGenerator) = true
+function gencut(::AvgCutGenerator, parent, path, stats, ztol)
+    (!path.childs_feasible || !path.childs_bounded) && return
+    avga = zeros(parent.nlds.nx)
+    avgβ = 0
+    for i in 1:length(parent.children)
+        @assert !isnull(path.childsols[i])
+        @assert get(path.childsols[i]).status != :Unbounded
+        childsol = get(path.childsols[i])
+        a = childsol.πT
+        β = childsol.πh + childsol.σd
+        @assert childsol.status == :Optimal
+        if isnull(parent.childT)
+            aT = a
+        else
+            aT = get(parent.childT)[i]' * a
+        end
+        β += childsol.ρe
+        avga += parent.proba[i] * aT
+        avgβ += parent.proba[i] * β
+    end
+    if path.sol.status == :Unbounded || _lt(path.sol.θ[1], avgβ - dot(avga, path.sol.x), ztol)
+        stats.ocutstime += @_time pushoptimalitycut!(parent, avga, avgβ, parent)
+        stats.nocuts += 1
+    end
+end
+function applycut(::AvgCutGenerator, node, g)
+    applyoptimalitycut!(node)
+end

src/interface.jl

@@ -1,10 +1,10 @@
-export model2lattice, SDDP, SDDPclear
+export model2lattice, getSDDPNode, SDDPclear
 
 type SDDPModelData
     nodes::Vector{Nullable{SDDPNode}}
 end
 
-function getSDDPNode(allnodes, m::Model, t, num_stages, solver, parent, pruningalgo::AbstractCutPruningAlgo, cutmode::Symbol, detectlb::Bool, newcut::Symbol)
+function getSDDPNode(m::Model, t, num_stages, solver, parent, pruningalgo::AbstractCutPruningAlgo, cutgen::AbstractCutGenerator, detectlb::Bool=true, newcut::Symbol=:InvalidateSolver)
     if !(:SDDP in keys(m.ext))
         nodes = Vector{Nullable{SDDPNode}}(num_stages)
         fill!(nodes, nothing)
@@ -16,17 +16,16 @@ function getSDDPNode(allnodes, m::Model, t, num_stages, solver, parent, pruninga
         c, T, W, h, C, K, _ = StructJuMP.conicconstraintdata(m)
         newnode = SDDPNode(NLDS{Float64}(W,h,T,K,C,c,solver,pruningalgo, newcut), parent)
         nodes[t] = newnode
-        push!(allnodes[t], newnode)
         struc = getStructure(m)
         if t < num_stages
             num_scen = length(struc.children)
             children = Vector{SDDPNode{Float64}}(num_scen)
             probability = Vector{Float64}(num_scen)
             for (i, id) in enumerate(keys(struc.children))
-                children[i] = getSDDPNode(allnodes, struc.children[id], t+1, num_stages, solver, newnode, pruningalgo, cutmode, detectlb, newcut)
+                children[i] = getSDDPNode(struc.children[id], t+1, num_stages, solver, newnode, pruningalgo, cutgen, detectlb, newcut)
                 probability[i] = struc.probability[id]
             end
-            setchildren!(newnode, children, probability, cutmode)
+            setchildren!(newnode, children, probability, cutgen)
             if detectlb
                 setθlb!(newnode, map(getobjlb, children))
             end
@@ -41,15 +40,10 @@ $(SIGNATURES)
 Transforms a [StructJuMP](https://github.com/StructJuMP/StructJuMP.jl) model into a lattice that can be used by the SDDP algorithm.
 The master problem is assumed to have model `m` and the scenarios are considered up to `num_stages` stages.
 The `pruningalgo` is as defined in [CutPruners](https://github.com/JuliaPolyhedra/CutPruners.jl).
-If `cutmode` is `:MultiCut`, one variable `θ_i` is created for each scenario. Otherwise, if `cutmode` is `:AveragedCut`, only one variable `θ` is created and it represents the expected value of the objective value of the scenarios. If `cutmode` is `:NoOptimalityCut` then no `θ` is created, only use this option if the objective of all models is zero except fo the master model.
+If `cutgen` is `MultiCutGenerator`, one variable `θ_i` is created for each scenario. Otherwise, if `cutgen` is `AveragedCutGenerator`, only one variable `θ` is created and it represents the expected value of the objective value of the scenarios. If `cutgen` is `NoOptimalityCut` then no `θ` is created, only use this option if the objective of all models is zero except fo the master model.
 """
-function model2lattice(m::Model, num_stages, solver, pruningalgo::AbstractCutPruningAlgo, cutmode::Symbol=:MultiCut, detectlb::Bool=true, newcut::Symbol=:InvalidateSolver)
-    nodes = Vector{Vector{SDDPNode}}(num_stages)
-    for i in 1:num_stages
-        nodes[i] = SDDPNode[]
-    end
-
-    root = getSDDPNode(nodes, m, 1, num_stages, solver, nothing, pruningalgo, cutmode, detectlb, newcut)
+function model2lattice(m::Model, num_stages, solver, pruningalgo::AbstractCutPruningAlgo, cutgen::AbstractCutGenerator=MultiCutGenerator, detectlb::Bool=true, newcut::Symbol=:InvalidateSolver)
+    root = getSDDPNode(m, 1, num_stages, solver, nothing, pruningalgo, cutgen, detectlb, newcut)
 	GraphSDDPTree(root)
 end
 

src/nlds.jl

@@ -24,14 +24,14 @@ end
 
 # Primal
 # min c x
-#     + θ                         -> :AveragedCut
-#     + sum θ_i                   -> :MultiCut
+#     + θ                         -> AvgCutGenerator
+#     + sum θ_i                   -> MultiCutGenerator
 #     + λ c_a x_a                 -> parent unbounded
 #     h -   Tx_a - Wx ∈ K         -> parent bounded
 #     h - λ Tx_a - Wx ∈ K         -> parent unbounded
 #     Dx >= d                     (feasibility cuts)
-#     Ex + θ   >= e               -> :AveragedCut
-#     Ex + θ_i >= e               -> :MultiCut
+#     Ex + θ   >= e               -> AvgCutGenerator
+#     Ex + θ_i >= e               -> MultiCutGenerator
 #     x in C
 #     λ >= 0
 
@@ -67,7 +67,7 @@ type NLDS{S}
     θlb::Vector{Float64}
     θC
     childT::Nullable{Vector{AbstractMatrix{S}}}
-    cutmode::Symbol
+    cutgen::AbstractCutGenerator
 
     nx::Int
     nθ::Int
@@ -105,7 +105,7 @@ type NLDS{S}
         else
             model = MathProgBase.LinearQuadraticModel(solver)
         end
-        nlds = new{S}(W, h, T, K, C, c, S[], nothing, nothing, CutStore{S}[], CutStore{S}[], localFC, localOC, nothing, Float64[], [], nothing, :NoOptimalityCut, nx, nθ, nπ, 1:nπ, 0, Int[], Int[], Vector{Int}[], model, false, false, nothing, newcut, pruningalgo, FCpruner, OCpruners)
+        nlds = new{S}(W, h, T, K, C, c, S[], nothing, nothing, CutStore{S}[], CutStore{S}[], localFC, localOC, nothing, Float64[], [], nothing, NoOptimalityCutGenerator(), nx, nθ, nπ, 1:nπ, 0, Int[], Int[], Vector{Int}[], model, false, false, nothing, newcut, pruningalgo, FCpruner, OCpruners)
         addfollower(localFC, (nlds, (:Feasibility, 0)))
         addfollower(localOC, (nlds, (:Optimality, 1)))
         nlds
@@ -116,17 +116,11 @@ function (::Type{NLDS{S}}){S}(W::AbstractMatrix, h::AbstractVector, T::AbstractM
     NLDS{S}(AbstractMatrix{S}(W), AbstractVector{S}(h), AbstractMatrix{S}(T), K, C, AbstractVector{S}(c), solver, pruningalgo, newcut)
 end
 
-function setchildren!{S}(nlds::NLDS{S}, childFC, childOC, proba, cutmode, childT)
-    nlds.cutmode = cutmode
+function setchildren!{S}(nlds::NLDS{S}, childFC, childOC, proba, cutgen::AbstractCutGenerator, childT)
+    nlds.cutgen = cutgen
     @assert length(childFC) == length(childOC) == length(proba)
     nlds.proba = proba
-    if cutmode == :MultiCut
-        nlds.nθ = length(proba)
-    elseif cutmode == :AveragedCut
-        nlds.nθ = 1
-    else
-        nlds.nθ = 0
-    end
+    nlds.nθ = nθ(nlds.cutgen, proba)
     nlds.θlb = zeros(nlds.nθ)
     nlds.θC = [(:Free, collect(nlds.nx+(1:nlds.nθ)))]
     nlds.nρ = zeros(Int, nlds.nθ)
@@ -193,9 +187,7 @@ function appendchildren!{S}(nlds::NLDS{S}, childFC, childOC, proba, childT)
     @assert length(proba) == length(nlds.childOC) + length(childOC) == length(nlds.childFC) + length(childFC)
     oldnθ = nlds.nθ
     nlds.proba = proba
-    if nlds.cutmode == :MultiCut
-        nlds.nθ = length(proba)
-    end
+    nlds.nθ = nθ(nlds.cutgen, proba)
     Δθ = nlds.nθ - oldnθ
     append!(nlds.nρ, zeros(Int, Δθ))
     append!(nlds.ρs, [Int[] for i in 1:Δθ])

src/node.jl

@@ -43,7 +43,7 @@ function Base.show(io::IO, node::SDDPNode)
     println(io, "ode of $(node.nvars) variables and outdegree of $(length(node.children)) with proba: $(node.proba)")
 end
 
-function setchildren!(node::SDDPNode, children, proba, cutmode, childT=nothing)
+function setchildren!(node::SDDPNode, children, proba, cutgen::AbstractCutGenerator, childT=nothing)
     @assert length(children) == length(proba)
     node.children = children
     node.proba = proba
@@ -52,7 +52,7 @@ function setchildren!(node::SDDPNode, children, proba, cutmode, childT=nothing)
     childOC = map(child -> child.ocuts, children)
     node.childT = childT
     empty!(node.npath)
-    setchildren!(node.nlds, childFC, childOC, proba, cutmode, childT)
+    setchildren!(node.nlds, childFC, childOC, proba, cutgen, childT)
 end
 
 getobjlb(node::SDDPNode) = getobjlb(node.nlds)

src/sddp.jl

@@ -1,11 +1,12 @@
+export SDDP
+
 type SDDPSolution
     status
     objval
     sol
     attrs
 end
 
-
 function meanstdpaths(z::Vector{Float64}, proba::Vector{Float64}, npaths::Vector{Int}, Ktot)
     if Ktot != -1
         @assert sum(npaths) == Ktot
@@ -141,7 +142,7 @@ function iteration{S}(g::AbstractSDDPTree{S}, Ktot::Int, num_stages, verbose, pa
                     npaths = samplepaths(pathsampler, g, state, path.K, t, num_stages)
 					childocuts = Array{Any}(nchildren(g, state))
                     for i in 1:nchildren(g, state)
-						if t == 2 || sum(npaths[i]) != 0 || cutmode(g, state) == :AveragedCut
+                        if t == 2 || sum(npaths[i]) != 0 || needallchildsol(cutgen(g, state))
 							addjob!(jobsd, getchild(g, state, i), SDDPJob(path.proba * getproba(g, state, i), npaths[i], state, path, i))
                         end
                     end
@@ -178,68 +179,31 @@ function iteration{S}(g::AbstractSDDPTree{S}, Ktot::Int, num_stages, verbose, pa
         # e = π h + ρ e + σ d
         for (parent, paths) in pathsd
             for path in paths
-                if parent.nlds.cutmode == :AveragedCut && path.childs_feasible && path.childs_bounded
-                    avga = zeros(parent.nlds.nx)
-                    avgβ = 0
-                end
-                for i in 1:length(parent.children)
-                    if isnull(path.childsols[i])
-                        @assert !path.childs_feasible || parent.nlds.cutmode != :AveragedCut
-                    elseif get(path.childsols[i]).status != :Unbounded
-                        childsol = get(path.childsols[i])
-                        a = childsol.πT
-                        β = childsol.πh + childsol.σd
-                        if !path.childs_feasible
+                if path.childs_feasible
+                    gencut(cutgen(g, parent), parent, path, stats, ztol)
+                else
+                    for i in 1:length(parent.children)
+                        if !isnull(path.childsols[i])
+                            childsol = get(path.childsols[i])
                             if childsol.status == :Infeasible
+                                a = childsol.πT
+                                β = childsol.πh + childsol.σd
                                 infeasibility_detected = true
                                 stats.nfcuts += 1
                                 stats.fcutstime += @_time pushfeasibilitycut!(parent.children[i], a, β, parent)
                             end
-                        elseif !(parent.nlds.cutmode == :AveragedCut && (!path.childs_bounded || !path.childs_feasible))
-                            @assert childsol.status == :Optimal
-                            if isnull(parent.childT)
-                                aT = a
-                            else
-                                aT = get(parent.childT)[i]' * a
-                            end
-                            β += childsol.ρe
-                            # This assumption is dropped now
-                           #if _lt(β - dot(aT, path.sol.x), .0, ztol)
-                           #    error("The objectives are supposed to be nonnegative")
-                           #end
-                            if parent.nlds.cutmode == :MultiCut
-                                if path.sol.status == :Unbounded || _lt(path.sol.θ[i], β - dot(aT, path.sol.x), ztol)
-                                    stats.ocutstime += @_time pushoptimalitycutforparent!(parent.children[i], a, β, parent)
-                                    stats.nocuts += 1
-                                end
-                            elseif parent.nlds.cutmode == :AveragedCut
-                                avga += parent.proba[i] * aT
-                                avgβ += parent.proba[i] * β
-                            end
                         end
                     end
                 end
-                if parent.nlds.cutmode == :AveragedCut && path.childs_feasible && path.childs_bounded
-                    if path.sol.status == :Unbounded || _lt(path.sol.θ[1], avgβ - dot(avga, path.sol.x), ztol)
-                        stats.ocutstime += @_time pushoptimalitycut!(parent, avga, avgβ, parent)
-                        stats.nocuts += 1
-                    end
-                end
             end
         end
 
         # Apply cut addition
         for (state, paths) in pathsd
-			for child in children(g, state)
+            for child in children(g, state)
                 applyfeasibilitycut!(child)
             end
-			if cutmode(g, state) == :MultiCut
-				for child in children(g, state)
-                    applyoptimalitycutforparent!(child)
-                end
-			elseif cutmode(g, state) == :AveragedCut
-                applyoptimalitycut!(state)
-            end
+            applycut(cutgen(g, state), state, g)
         end
 
         # Jobs -> Paths

src/sddptree.jl

@@ -1,4 +1,4 @@
-export AbstractSDDPTree, haschildren, nchidren, children, getchild, getproba, getprobas, cutmode, numberofpaths
+export AbstractSDDPTree, haschildren, nchidren, children, getchild, getproba, getprobas, cutgen, numberofpaths
 @compat abstract type AbstractSDDPTree{S} end
 
 type GraphSDDPTree{S} <: AbstractSDDPTree{S}
@@ -20,4 +20,4 @@ getprobas(g::GraphSDDPTree, node::SDDPNode) = node.proba
 numberofpaths(g::GraphSDDPTree, num_stages) = numberofpaths(g.root, 1, num_stages)
 numberofpaths(g::GraphSDDPTree, node::SDDPNode, t, num_stages) = numberofpaths(node, t, num_stages)
 
-cutmode(g::GraphSDDPTree, node::SDDPNode) = node.nlds.cutmode
+cutgen(g::GraphSDDPTree, node::SDDPNode) = node.nlds.cutgen

test/hydro_thermal_scheduling.jl

@@ -34,7 +34,7 @@
 		end
 	end
 
-    for cutmode in [:MultiCut, :AveragedCut]
+    for cutmode in [MultiCutGenerator(), AvgCutGenerator()]
         for detectlb in [false, true]
             !detectlb && iscpx(solver) && continue
             lattice = model2lattice(models[1], num_stages, solver, AvgCutPruningAlgo(-1), cutmode, detectlb)