Merge pull request #14 from blegat/stochprog

Creation of AbstractStochasticProgram interface

REQUIRE

@@ -4,3 +4,4 @@ MathProgBase 0.5
 JuMP 0.17
 StructJuMP 0.0.1
 DocStringExtensions 0.2
+LightGraphs

src/StructDualDynProg.jl

@@ -13,27 +13,31 @@ import Base.show, Base.isless
 
 # Utils
 include("comp.jl")
-include("stats.jl")
 
 # Abstract components
+# Stochastic Program
+include("stochprog.jl")
+# Stats
+include("stats.jl")
 # Stopping Criterion
 include("stopcrit.jl")
 # Cut Generator
 include("cutgen.jl")
-# NLDS Model
-include("cutstore.jl")
-include("solver.jl")
-include("nlds.jl")
-# SDDP Graph
-include("node.jl")
-include("graph.jl")
 # Path Sampler
 include("sampler.jl")
+# Solution
+include("solution.jl")
 
-# SDDP algorithm
+# SDDP algorithm on top of these abstract components
 include("path.jl")
 include("sddp.jl")
 
+# Generic implementation of Stochastic Program
+include("cutstore.jl")
+include("solver.jl")
+include("nlds.jl")
+include("graph.jl")
+
 # Wait and See value
 include("waitandsee.jl")
 

src/cutgen.jl

@@ -5,20 +5,20 @@ abstract type AbstractOptimalityCutGenerator <: AbstractCutGenerator end
 
 struct FeasibilityCutGenerator <: AbstractCutGenerator
 end
-function gencut(::FeasibilityCutGenerator, g, parent, path, stats, ztol)
-    for i in 1:nchildren(g, parent)
-        if !isnull(path.childsols[i])
-            childsol = get(path.childsols[i])
-            if childsol.status == :Infeasible
+function gencut(::FeasibilityCutGenerator, sp, parent, path, stats, ztol)
+    for child in out_neighbors(sp, parent)
+        if haskey(path.childsols, child)
+            childsol = path.childsols[child]
+            if getstatus(childsol) == :Infeasible
                 stats.nfcuts += 1
-                stats.fcutstime += @_time pushfeasibilitycut!(parent.children[i], getfeasibilitycut(childsol)..., parent)
+                stats.fcutstime += @_time add_feasibility_cut!(sp, child, feasibility_cut(childsol)..., parent)
             end
         end
     end
 end
-function applycut(::FeasibilityCutGenerator, g, node)
-    for child in children(g, node)
-        applyfeasibilitycut!(child)
+function applycut(::FeasibilityCutGenerator, sp, node)
+    for child in out_neighbors(sp, node)
+        apply_feasibility_cuts!(sp, child)
     end
 end
 
@@ -27,38 +27,41 @@ struct NoOptimalityCutGenerator <: AbstractOptimalityCutGenerator
 end
 nθ(::NoOptimalityCutGenerator, proba) = 0
 needallchildsol(::NoOptimalityCutGenerator) = false
-function gencut(::NoOptimalityCutGenerator, g, parent, path, stats, ztol)
+function gencut(::NoOptimalityCutGenerator, sp, parent, path, stats, ztol)
 end
-function applycut(::NoOptimalityCutGenerator, g, node)
-    applyoptimalitycut!(node)
+function applycut(::NoOptimalityCutGenerator, sp, node)
+    apply_optimality_cuts!(sp, node)
 end
 
 # Multi Cut Generator
 struct MultiCutGenerator <: AbstractOptimalityCutGenerator
 end
 nθ(::MultiCutGenerator, proba) = length(proba)
 needallchildsol(::MultiCutGenerator) = false
-function gencut(::MultiCutGenerator, g, 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, β = getoptimalitycut(childsol)
-            @assert childsol.status == :Optimal
-            if isnull(parent.childT)
-                aT = a
+function gencut(::MultiCutGenerator, sp, parent, path, stats, ztol)
+    for (child, sol) in path.childsols
+        status = getstatus(sol)
+        if status != :Unbounded
+            a, β = optimality_cut(sol)
+            @assert status == :Optimal
+            edge = Edge(parent, child)
+            if haskey(sp.childT, edge)
+                aT = sp.childT[edge]' * a
             else
-                aT = get(parent.childT)[i]' * a
+                aT = a
             end
-            if path.sol.status == :Unbounded || _lt(path.sol.θ[i], β - dot(aT, path.sol.x), ztol)
-                stats.ocutstime += @_time pushoptimalitycutforparent!(parent.children[i], a, β, parent)
+            x = getstatevalue(path.sol)
+            θ = getθvalue(sp, parent, child, path.sol)
+            if getstatus(path.sol) == :Unbounded || _lt(θ, β - dot(aT, x), ztol)
+                stats.ocutstime += @_time add_optimality_cut_for_parent!(sp, child, a, β, parent)
                 stats.nocuts += 1
             end
         end
     end
 end
-function applycut(::MultiCutGenerator, g, node)
-    for child in children(g, node)
-        applyoptimalitycutforparent!(child)
+function applycut(::MultiCutGenerator, sp, node)
+    for child in out_neighbors(sp, node)
+        apply_optimality_cuts_for_parent!(sp, child)
     end
 end
 
@@ -67,29 +70,32 @@ struct AvgCutGenerator <: AbstractOptimalityCutGenerator
 end
 nθ(::AvgCutGenerator, proba) = 1
 needallchildsol(::AvgCutGenerator) = true
-function gencut(::AvgCutGenerator, g, parent, path, stats, ztol)
+function gencut(::AvgCutGenerator, sp, 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, β = getoptimalitycut(childsol)
-        @assert childsol.status == :Optimal
-        if isnull(parent.childT)
-            aT = a
+    avga = zeros(statedim(sp, parent))
+    avgβ = 0.
+    for (child, sol) in path.childsols
+        status = getstatus(sol)
+        @assert status != :Unbounded
+        a, β = optimality_cut(sol)
+        @assert status == :Optimal
+        edge = Edge(parent, child)
+        if haskey(sp.childT, edge)
+            aT = sp.childT[edge]' * a
         else
-            aT = get(parent.childT)[i]' * a
+            aT = a
         end
-        avga += parent.proba[i] * aT
-        avgβ += parent.proba[i] * β
+        proba = probability(sp, edge)
+        avga += proba * aT
+        avgβ += proba * β
     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)
+    x = getstatevalue(path.sol)
+    θ = getθvalue(sp, parent, path.sol)
+    if getstatus(path.sol) == :Unbounded || _lt(θ, avgβ - dot(avga, x), ztol)
+        stats.ocutstime += @_time add_optimality_cut!(sp, parent, avga, avgβ, parent)
         stats.nocuts += 1
     end
 end
-function applycut(::AvgCutGenerator, g, node)
-    applyoptimalitycut!(node)
+function applycut(::AvgCutGenerator, sp, node)
+    apply_optimality_cuts!(sp, node)
 end

src/cutstore.jl

@@ -35,7 +35,7 @@ mutable struct CutStore{S}
 
     function CutStore{S}(nvars) where {S}
         # spzeros(S, 0) -> S[] : See julia#22225
-        new{S}(spzeros(S, 0, nvars), S[], NLDS{S}[], spzeros(S, 0, nvars), S[], NLDS{S}[], Vector{Tuple{NLDS{S},Tuple{Symbol,Int64}}}(0), Vector{Bool}(0), :IfNeededElseDelete)
+        new{S}(spzeros(S, 0, nvars), S[], NLDS{S}[], spzeros(S, 0, nvars), S[], NLDS{S}[], Vector{Tuple{NLDS{S},Tuple{Symbol,Int}}}(0), Vector{Bool}(0), :IfNeededElseDelete)
     end
 end
 

src/graph.jl

@@ -1,23 +1,160 @@
-export AbstractSDDPGraph, haschildren, nchidren, children, getchild, getproba, getprobas, cutgen, numberofpaths
-abstract type AbstractSDDPGraph{S} end
+using LightGraphs
 
-mutable struct SDDPGraph{S} <: AbstractSDDPGraph{S}
-	root::SDDPNode{S}
+import MathProgBase
+
+mutable struct NodeData{S}
+    nlds::NLDS{S}
+    npath::Dict{Int, Int}
+
+    # Feasibility cuts
+    fcuts::CutStore{S}
+    # Optimality cuts
+    ocuts::CutStore{S}
+
+    function NodeData{S}(nlds::NLDS{S}, nvars_a) where S
+        new{S}(nlds, Dict{Int, Int}(), CutStore{S}(nvars_a), CutStore{S}(nvars_a))
+    end
+end
+
+NodeData(nlds::NLDS{S}, parent) where {S} = NodeData{S}(nlds, parent)
+
+function Base.show(io::IO, data::NodeData)
+    println(io, "Node of $(data.nlds.nx) variables")
 end
 
-getmaster(g::SDDPGraph) = g.root, g.root
+ET = LightGraphs.SimpleGraphs.SimpleEdge{Int}
+
+mutable struct StochasticProgram{S} <: AbstractStochasticProgram
+    graph::LightGraphs.SimpleGraphs.SimpleDiGraph{Int}
+    eid::Dict{ET, Int}
+    proba::Dict{ET, S}
+    childT::Dict{ET, AbstractMatrix{S}}
+    data::Vector{NodeData{S}}
+    function StochasticProgram{S}() where S
+        new{S}(DiGraph(), Dict{ET, Int}(), Dict{ET, S}(), Dict{ET, AbstractMatrix{S}}(), NodeData{S}[])
+    end
+end
+nodedata(sp::StochasticProgram, node::Int) = sp.data[node]
 
-# Get children scenarios
-haschildren(g::SDDPGraph, node::SDDPNode) = !isempty(node.children)
-nchildren(g::SDDPGraph, node::SDDPNode) = length(node.children)
-children(g::SDDPGraph, node::SDDPNode) = node.children
-getchild(g::SDDPGraph, node::SDDPNode, i) = node.children[i]
-# Get proba of children scenario
-getproba(g::SDDPGraph, node::SDDPNode, i) = node.proba[i]
-getprobas(g::SDDPGraph, node::SDDPNode) = node.proba
+getobjectivebound(sp::StochasticProgram, node) = getobjectivebound(nodedata(sp, node).nlds)
+setθbound!(sp::StochasticProgram, node, child, θlb) = setθbound!(nodedata(sp, node).nlds, edgeid(sp, Edge(node, child)), θlb)
+statedim(sp::StochasticProgram, node) = nodedata(sp, node).nlds.nx
 
-# Get number of paths
-numberofpaths(g::SDDPGraph, num_stages) = numberofpaths(g.root, 1, num_stages)
-numberofpaths(g::SDDPGraph, node::SDDPNode, t, num_stages) = numberofpaths(node, t, num_stages)
+# LightGraphs interface
+LightGraphs.out_neighbors(sp::StochasticProgram, node::Int) = out_neighbors(sp.graph, node)
 
-cutgen(g::SDDPGraph, node::SDDPNode) = node.nlds.cutgen
+getmaster(sp::StochasticProgram) = 1
+
+# If the graph is not a tree, this will loop if I don't use a num_stages limit
+function numberofpaths(sp::StochasticProgram, node, len)
+    @assert len >= 0
+    if iszero(len) || isleaf(sp, node)
+        1
+    else
+        npath = nodedata(sp, node).npath
+        if !(len in keys(npath))
+            npath[len] = sum(map(c -> numberofpaths(sp, c, len-1), out_neighbors(sp, node)))
+        end
+        npath[len]
+    end
+end
+
+cutgenerator(sp::StochasticProgram, node) = nodedata(sp, node).nlds.cutgen
+function setcutgenerator!(sp::StochasticProgram, node, cutgen::AbstractOptimalityCutGenerator)
+    nodedata(sp, node).nlds.cutgen = cutgen
+end
+
+function add_scenario_state!(sp::StochasticProgram, data::NodeData)
+    @assert add_vertex!(sp.graph)
+    push!(sp.data, data)
+    @assert nv(sp.graph) == length(sp.data)
+    length(sp.data)
+end
+
+function add_scenario_transition!(sp::StochasticProgram, parent, child, proba, childT=nothing)
+    edge = Edge(parent, child)
+    if !add_edge!(sp.graph, edge)
+        error("Edge already in the graph, multiple edges not supported yet")
+    end
+    @assert !haskey(sp.eid, edge)
+    @assert !haskey(sp.proba, edge)
+    @assert !haskey(sp.childT, edge)
+    data = nodedata(sp, parent)
+    sp.eid[edge] = outdegree(sp, parent)
+    sp.proba[edge] = proba
+    if childT !== nothing
+        sp.childT[edge] = childT
+    end
+    empty!(data.npath)
+    childdata = nodedata(sp, child)
+    add_scenario_transition!(data.nlds, childdata.fcuts, childdata.ocuts, proba, childT)
+    @assert length(data.nlds.childFC) == length(data.nlds.proba) == outdegree(sp, parent)
+end
+
+probability(sp::StochasticProgram, edge) = sp.proba[edge]
+
+function setprobability!(sp::StochasticProgram, edge, proba)
+    sp.proba[edge] = proba
+    data = nodedata(sp, src(edge))
+    setprobability!(data.nlds, edgeid(sp, edge), proba)
+end
+
+function edgeid(sp::StochasticProgram, edge)
+    sp.eid[edge]
+end
+
+function solve!(sp::StochasticProgram, node)
+    getsolution(nodedata(sp, node).nlds)
+end
+
+function setchildx!(sp::StochasticProgram, node, child, sol::Solution)
+    data = nodedata(sp, node)
+    edge = Edge(node, child)
+    if haskey(sp.childT, edge)
+        T = data.childT[edge]
+        x = T * sol.x
+        if sol.xuray !== nothing
+            xuray = T * sol.xuray
+        end
+    else
+        x = sol.x
+        xuray = sol.xuray
+    end
+    setparentx(nodedata(sp, child).nlds, x, xuray, sol.objvalxuray)
+end
+
+function getθvalue(sp::StochasticProgram, node, child, sol::Solution)
+    @assert length(sol.θ) == outdegree(sp, node)
+    getθvalue(sol, edgeid(sp, Edge(node, child)))
+end
+
+function getθvalue(sp::StochasticProgram, node, sol::Solution)
+    @assert length(sol.θ) == 1
+    getθvalue(sol, 1)
+end
+
+function add_feasibility_cut!(sp::StochasticProgram, node, coef, rhs, author)
+    # coef is a ray
+    # so alpha * coef is also valid for any alpha >= 0.
+    # Hence coef might have very large coefficients and alter
+    # the numerial accuracy of the master's solver.
+    # We scale it to avoid this issue
+    scaling = max(abs(rhs), maximum(abs, coef))
+    addcut(nodedata(sp, node).fcuts, coef/scaling, sign(rhs), nodedata(sp, author).nlds)
+end
+function add_optimality_cut!(sp::StochasticProgram, node, coef, rhs, author)
+    addcut(nodedata(sp, node).nlds.localOC, coef, rhs, nodedata(sp, author).nlds)
+end
+function add_optimality_cut_for_parent!(sp::StochasticProgram, node, coef, rhs, author)
+    addcut(nodedata(sp, node).ocuts, coef, rhs, nodedata(sp, author).nlds)
+end
+
+function apply_feasibility_cuts!(sp::StochasticProgram, node)
+    apply!(nodedata(sp, node).fcuts)
+end
+function apply_optimality_cuts!(sp::StochasticProgram, node)
+    apply!(nodedata(sp, node).nlds.localOC)
+end
+function apply_optimality_cuts_for_parent!(sp::StochasticProgram, node)
+    apply!(nodedata(sp, node).ocuts)
+end