Merge pull request #1 from kibaekkim/dev

Dev

README.md

@@ -13,29 +13,34 @@ Pkg.clone("https://github.com/kibaekkim/Dsp.jl")
 DSP can read and solve model from JuMP:
 
 ```julia
-using Dsp, JuMP
+using JuMP, Dsp, MPI
+
+# Comment out this line if you want to run in serial
+MPI.Init()
 
 xi = [[7,7] [11,11] [13,13]]
 
-m = Model()
+# create JuMP.Model with number of blocks
+m = Model(3)
+
 @variable(m, 0 <= x[i=1:2] <= 5, Int)
 @objective(m, Min, -1.5 * x[1] - 4 * x[2])
 for s = 1:3
-    blk = Model()
+    # create a JuMP.Model block linked to m with id s and probability 1/3
+    blk = Model(m, s, 1/3)
     @variable(blk, y[j=1:4], Bin)
     @objective(blk, Min, -16 * y[1] + 19 * y[2] + 23 * y[3] + 28 * y[4])
     @constraint(blk, 2 * y[1] + 3 * y[2] + 4 * y[3] + 5 * y[4] <= xi[1,s] - x[1])
     @constraint(blk, 6 * y[1] + y[2] + 3 * y[3] + 2 * y[4] <= xi[2,s] - x[2])
-
-    # Model m has block blk indexed by s with objective function weight of 1/3.
-    @block(m, blk, s, 1/3)
-
 end
 
 solve_types = [:Dual, :Benders, :Extensive]
 solve(m, solve_type = solve_types[1], param = "myparam.txt")
 
 getobjectivevalue(m)
+
+# Comment out this line if you want to run in serial
+MPI.Finalize()
 ```
 
 or, it can also read and solve model from SMPS files:

examples/farmer.jl

@@ -20,15 +20,14 @@ Yield    = [3.0 3.6 24.0;
 Minreq   = [200 240 0]     # minimum crop requirement
 
 # JuMP model
-# m = Model(solver = DspSolver(method = :DD, param = "myparams.txt"))
-m = Model()
+m = Model(NS)
 
 @variable(m, x[i=CROPS] >= 0, Int)
 @objective(m, Min, sum{Cost[i] * x[i], i=CROPS})
 @constraint(m, const_budget, sum{x[i], i=CROPS} <= Budget)
 
 for s in 1:NS
-    blk = Model()
+    blk = Model(m, s, probability[s])
 
     @variable(blk, y[j=PURCH] >= 0)
     @variable(blk, w[k=SELL] >= 0)
@@ -38,7 +37,4 @@ for s in 1:NS
     @constraint(blk, const_minreq[j=PURCH], Yield[s,j] * x[j] + y[j] - w[j] >= Minreq[j])
     @constraint(blk, const_minreq_beets, Yield[s,3] * x[3] - w[3] - w[4] >= Minreq[3])
     @constraint(blk, const_aux, w[3] <= 6000)
-
-    # add blk to m as a block
-    @block(m, blk, s, probability[s])
 end

examples/farmer_mpi.jl

@@ -0,0 +1,46 @@
+# Kibaek Kim - ANL MCS 2016
+# Farmer example from Birge and Louveaux book.
+#
+# To use MPI, add MPI.Init() and MPI.Finalize() around the scritp lines.
+
+using MPI, JuMP, Dsp
+
+MPI.Init()
+
+NS = 3;                        # number of scenarios
+probability = [1/3, 1/3, 1/3]; # probability
+
+CROPS  = 1:3 # set of crops (wheat, corn and sugar beets, resp.)
+PURCH  = 1:2 # set of crops to purchase (wheat and corn, resp.)
+SELL   = 1:4 # set of crops to sell (wheat, corn, sugar beets under 6K and those over 6K)
+
+Cost     = [150 230 260]   # cost of planting crops
+Budget   = 500             # budget capacity
+Purchase = [238 210];      # purchase price
+Sell     = [170 150 36 10] # selling price
+Yield    = [3.0 3.6 24.0;
+            2.5 3.0 20.0;
+            2.0 2.4 16.0]
+Minreq   = [200 240 0]     # minimum crop requirement
+
+# JuMP model
+m = Model(NS)
+
+@variable(m, x[i=CROPS] >= 0, Int)
+@objective(m, Min, sum{Cost[i] * x[i], i=CROPS})
+@constraint(m, const_budget, sum{x[i], i=CROPS} <= Budget)
+
+for s in 1:NS
+    blk = Model(m, s, probability[s])
+
+    @variable(blk, y[j=PURCH] >= 0)
+    @variable(blk, w[k=SELL] >= 0)
+
+    @objective(blk, Min, sum{Purchase[j] * y[j], j=PURCH} - sum{Sell[k] * w[k], k=SELL})
+
+    @constraint(blk, const_minreq[j=PURCH], Yield[s,j] * x[j] + y[j] - w[j] >= Minreq[j])
+    @constraint(blk, const_minreq_beets, Yield[s,3] * x[3] - w[3] - w[4] >= Minreq[3])
+    @constraint(blk, const_aux, w[3] <= 6000)
+end
+
+MPI.Finalize()

src/Dsp.jl

@@ -1,4 +1,4 @@
-__precompile__()
+# __precompile__()
 
 module Dsp
 
@@ -8,21 +8,60 @@ include("block.jl")
 using Dsp.DspCInterface
 
 import JuMP
-export readSmps, getblocksolution, getdualobjval
+export
+    readSmps, 
+    getblocksolution, 
+    optimize, 
+    getprimobjval, 
+    getdualobjval, 
+    getprimvalue,
+    getdualvalue,
+    getsolutiontime, 
+    blockids
 
 # DspModel placeholder
 model = DspModel()
 
 ###############################################################################
-# JuMP.solvehook
+# Override JuMP.Model
 ###############################################################################
 
-# This function is hooked by JuMP (see block.jl)
-function dsp_solve(m::JuMP.Model; suppress_warnings = false, options...)
-
+# This is for the master problem.
+function JuMP.Model(nblocks::Integer)
+    # construct model
+    m = JuMP.Model()
     # free DspModel
     DspCInterface.freeModel(Dsp.model)
+    # set block ids
+    DspCInterface.setBlockIds(Dsp.model, nblocks)
+    # set extension
+    m.ext[:DspBlocks] = BlockStructure(nothing, Dict{Int,JuMP.Model}(), Dict{Int,Float64}())
+    # set solvehook
+    JuMP.setsolvehook(m, Dsp.dsp_solve)
+    # return model
+    m
+end
+
+# This is for the subproblems.
+function JuMP.Model(parent::JuMP.Model, id::Integer, weight::Float64)
+    # construct model
+    m = JuMP.Model()
+    # set extension
+    m.ext[:DspBlocks] = BlockStructure(nothing, Dict{Int,JuMP.Model}(), Dict{Int,Float64}())
+    # set block
+    m.ext[:DspBlocks].parent = parent
+    setindex!(parent.ext[:DspBlocks].children, m, id)
+    setindex!(parent.ext[:DspBlocks].weight, weight, id)
+    # return model
+    m
+end
+
+###############################################################################
+# JuMP.solvehook
+###############################################################################
 
+# This function is hooked by JuMP (see block.jl)
+function dsp_solve(m::JuMP.Model; suppress_warnings = false, options...)
     # parse options
     for (optname, optval) in options
         if optname == :param
@@ -60,37 +99,11 @@ function dsp_solve(m::JuMP.Model; suppress_warnings = false, options...)
     end
 
     if !(stat == :Infeasible || stat == :Unbounded)
-        try
-            getDspSolution(suppress_warnings)
-
-            # Don't corrupt the answers if one of the above two calls fails
-            m.objVal = Dsp.model.primVal
-            # parse solution to each block
-            n_start = 1
-            n_end = m.numCols
-            m.colVal = copy(Dsp.model.colVal[n_start:n_end])
-            n_start += m.numCols
-            if haskey(m.ext, :DspBlocks) == true
-                blocks = m.ext[:DspBlocks].children
-                for s in keys(blocks)
-                    n_end += blocks[s].numCols
-                    blocks[s].colVal = Dsp.model.colVal[n_start:n_end]
-                    n_start += blocks[s].numCols
-                end
-            end
-        end
-    end
-
-    # maximization?
-    if m.objSense == :Max
-        m.objVal *= -1
-        dsp.primVal *= -1
-        dsp.dualVal *= -1
+        getDspSolution(m)
     end
 
     # Return the solve status
     stat
-
 end
 
 ###############################################################################
@@ -130,7 +143,7 @@ function optimize(;suppress_warnings = false, options...)
     end
 
     if !(stat == :Infeasible || stat == :Unbounded)
-        getDspSolution(suppress_warnings)
+        getDspSolution()
     end
 
     # Return the solve status
@@ -140,10 +153,8 @@ JuMP.solve(;suppress_warnings = false, options...) = optimize(;suppress_warnings
 
 # Read model from SMPS files
 function readSmps(filename::AbstractString)
-
     # free DspModel
     DspCInterface.freeModel(Dsp.model)
-
     # read Smps file
     DspCInterface.readSmps(Dsp.model, filename)
 end
@@ -170,8 +181,14 @@ end
 getprimobjval() = Dsp.model.primVal
 # get dual objective value
 getdualobjval() = Dsp.model.dualVal
+# get primal value
+getprimvalue() = Dsp.model.colVal
 # get dual value
-JuMP.getdual() = Dsp.model.rowVal
+getdualvalue() = Dsp.model.rowVal
+# get solution time
+getsolutiontime() = DspCInterface.getWallTime(Dsp.model)
+# get block ids
+blockids() = Dsp.model.block_ids
 
 function parseStatusCode(statcode::Integer)
     stat = :NotSolved
@@ -199,15 +216,45 @@ function parseStatusCode(statcode::Integer)
     stat
 end
 
-function getDspSolution(suppress_warnings)
+function getDspSolution(m::JuMP.Model = nothing)
     Dsp.model.primVal = DspCInterface.getPrimalBound(Dsp.model)
     Dsp.model.dualVal = DspCInterface.getDualBound(Dsp.model)
 
     if Dsp.model.solve_type == :Dual
-        suppress_warnings || warn("Primal solution is not available in dual decomposition")
+        Dsp.model.rowVal = DspCInterface.getDualSolution(Dsp.model)
     else
         Dsp.model.colVal = DspCInterface.getSolution(Dsp.model)
+        if m != nothing
+            # parse solution to each block
+            n_start = 1
+            n_end = m.numCols
+            m.colVal = Dsp.model.colVal[n_start:n_end]
+            n_start += m.numCols
+            if haskey(m.ext, :DspBlocks) == true
+                numBlockCols = DspCInterface.getNumBlockCols(Dsp.model, m)
+                blocks = m.ext[:DspBlocks].children
+                for i in 1:Dsp.model.nblocks
+                    n_end += numBlockCols[i]
+                    if haskey(blocks, i)
+                        # @show b
+                        # @show n_start
+                        # @show n_end
+                        blocks[i].colVal = Dsp.model.colVal[n_start:n_end]
+                    end
+                    n_start += numBlockCols[i]
+                end
+            end
+        end
+    end
+    if m != nothing
+        m.objVal = Dsp.model.primVal
+        # maximization?
+        if m.objSense == :Max
+            m.objVal *= -1
+            dsp.primVal *= -1
+            dsp.dualVal *= -1
+        end
     end
 end
 
-end # module
+end # module