Merge bed012f into f4519cd

.travis.yml

@@ -1,23 +1,23 @@
-language: cpp
-compiler:
-    - gcc
+language: julia
+os:
+    - linux
+julia:
+    - nightly
 notifications:
     email: false
-env:
-    - JULIAVERSION="juliareleases"
-    - JULIAVERSION="julianightlies"
-before_install:
-    - sudo add-apt-repository ppa:staticfloat/julia-deps -y
-    - sudo add-apt-repository ppa:staticfloat/${JULIAVERSION} -y
-    - sudo apt-get update -qq -y
-    - sudo apt-get install libpcre3-dev libgmp-dev julia -y
-    - git config --global user.name "Travis User"
-    - git config --global user.email "travis@example.net"
-    - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
+sudo: false
+addons:
+    apt_packages:
+        - gfortran
+        - liblapack-dev
+        - libgmp-dev
+        - libglpk-dev
 script:
-    - julia -e 'versioninfo(); Pkg.init(); Pkg.clone(pwd())'
+    - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
+    - julia -e 'Pkg.add("JuMP"); Pkg.checkout("JuMP")'
+    - julia -e 'Pkg.clone(pwd())'
     - julia -e 'Pkg.test("JuMPeR", coverage=true)'
-    - if [ $JULIAVERSION = "juliareleases" ]; then julia -e 'Pkg.clone("https://github.com/vgupta1/DDUS.jl.git")'; fi
-    - if [ $JULIAVERSION = "juliareleases" ]; then julia -e 'Pkg.test("DDUS")'; fi
 after_success:
-    - julia -e 'cd(Pkg.dir("JuMPeR")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'
+    - julia -e 'Pkg.add("Coverage")'
+    - julia -e 'cd(Pkg.dir("JuMPeR")); using Coverage; Coveralls.submit(process_folder())'
+    - julia -e 'cd(Pkg.dir("JuMPeR")); using Coverage; Codecov.submit(process_folder())'

README.md

@@ -1,6 +1,6 @@
 ![JuMPeR Logo](http://iainnz.github.io/JuMPeR.jl/logo.svg)
 
-[![Build Status](https://travis-ci.org/IainNZ/JuMPeR.jl.png?branch=master)](https://travis-ci.org/IainNZ/JuMPeR.jl)
+[![Build Status](https://travis-ci.org/IainNZ/JuMPeR.jl.svg?branch=master)](https://travis-ci.org/IainNZ/JuMPeR.jl)
 [![Coverage Status](https://img.shields.io/coveralls/IainNZ/JuMPeR.jl.svg)](https://coveralls.io/r/IainNZ/JuMPeR.jl?branch=master)
 [![JuMPeR](http://pkg.julialang.org/badges/JuMPeR_release.svg)](http://pkg.julialang.org/?pkg=JuMPeR&ver=release)
 

REQUIRE

@@ -1,3 +1,3 @@
-julia 0.3
-JuMP 0.9
+julia 0.4-
+JuMP 0.10 0.11
 Compat

example/portfolio.jl

@@ -4,39 +4,43 @@
 # See http://github.com/IainNZ/JuMPeR.jl
 #############################################################################
 # portfolio.jl
-# Use the CovarOracle defined in oracle_covar.jl to solve an asset
-# allocation problem. Solves for multiple values of Γ and shows the
-# distribution of returns of the optimal portfolio for each value.
+# Solve a robust portfolio optimization problem. Can choose between:
+# - Polyhedral or ellipsoidal uncertainty sets
+# - Reformulation or cutting planes
 #############################################################################
 
-using JuMPeR
+using JuMP, JuMPeR
 using Distributions
-using Gurobi
-include("oracle_covar.jl")
+using Compat
 
-const NUM_ASSET = 10
+# Number of stocks
+const NUM_ASSET = 5
+# Number of samples of past returns
+const NUM_SAMPLES = 1000
+# Uncertainty set type (:Polyhedral or :Ellipsoidal)
+const UNCERTAINY_SET = :Polyhedral
 
 # Seed the RNG for reproducibilities sake
-srand(10)
+srand(1234)
 
 #############################################################################
 # generate_data
 # Sample returns of the assets. The parameters of the assets are fixed, we
 # simply make num_samples draws from the joint distribution.
-# Returns matrix with the samples in rows, each column is an assets
-function generate_data(num_samples::Int)
-    data = zeros(num_samples, NUM_ASSET)
+# Returns matrix with the samples in rows, each column is an asset
+function generate_data()
+    data = zeros(NUM_SAMPLES, NUM_ASSET)
 
     # Linking factors
     beta = [(i-1.0)/NUM_ASSET for i = 1:NUM_ASSET] 
 
-    for sample_ind = 1:num_samples
+    for sample_ind in 1:NUM_SAMPLES
         # Common market factor, mean 3%, sd 5%, truncate at +- 3 sd
         z = rand(Normal(0.03, 0.05))
         z = max(z, 0.03 - 3*0.05)
         z = min(z, 0.03 + 3*0.05)
 
-        for asset_ind = 1:NUM_ASSET
+        for asset_ind in 1:NUM_ASSET
             # Idiosyncratic contribution, mean 0%, sd 5%, truncated at +- 3 sd
             asset = rand(Normal(0.00, 0.05))
             asset = max(asset, 0.00 - 3*0.05)
@@ -51,42 +55,65 @@ end
 
 #############################################################################
 # solve_portfolio
-# Solve the robust portfolio problem given a matrix of past returns and
-# a degree of conservatism, Γ. The particular problem solved is
+# Solve the robust portfolio problem given a matrix of past returns.
+# The particular problem solved is
 #   max   z
-#   s.t.  ∑ xi    = 1
-#         ∑ ri*xi ≥ z  ∀ r ∈ R
+#   s.t.  ∑ ri*xi ≥ z  ∀ r ∈ R      # Worst-case return
+#         ∑ xi    = 1               # Is a portfolio
 #         xi ≥ 0
-#   where R is CovarOracle's set, an ellipse centered on the mean return
-#         and tilted using the covariance of the returns
+# where R is the uncertainty set:
+# R = { (r,z) | r = μ + Σ^0.5 z,
+#        Polyhedral:  ‖z‖₁ ≤ Γ, ‖z‖∞ ≤ 1
+#       Ellipsoidal:  ‖z‖₂ ≤ Γ
+#     }
 function solve_portfolio(past_returns, Γ)
     # Setup the robust optimization model
-    m = RobustModel(solver=GurobiSolver(OutputFlag=0))
+    m = RobustModel()
 
-    # Create the CovarOracle
-    setDefaultOracle!(m, CovarOracle(past_returns, Γ))
-
-    # Variables
-    @defVar(m, obj)  # Put objective as constraint (epigraph form)
-    @defVar(m, 0 <= x[1:NUM_ASSET] <= 1)  # Fractional allocation
-
-    # Uncertainties
-    @defUnc(m, r[1:NUM_ASSET])  # The returns
+    # Each asset is a share of the money to invest...
+    @defVar(m, 0 <= x[1:NUM_ASSET] <= 1)
+    # ... and we must allocate all the money.
+    @addConstraint(m, sum(x) == 1)
 
+    # JuMPeR doesn't support uncertain objectives directly. To have an
+    # uncertain objective, put objective as constraint (epigraph form)
+    @defVar(m, obj)
     @setObjective(m, Max, obj)
 
-    # Portfolio constraint
-    @addConstraint(m, sum(x) == 1)
+    # The uncertain parameters are the returns for each asset
+    @defUnc(m, r[1:NUM_ASSET])
+
+    # The uncertainty set requires the "square root" of the covariance
+    μ = vec(mean(past_returns, 1))  # want column vector 
+    Σ = cov(past_returns)
+    #L = rand(NUM_ASSET,NUM_ASSET)
+    L = full(@compat chol(Σ, Val{:L}))
+    # Define auxiliary uncertain parameters to model underlying factors
+    @defUnc(m, z[1:NUM_ASSET])
+    # Link r with z
+    @addConstraint(m, r .== L*z + μ)
+    if UNCERTAINY_SET == :Polyhedral
+        # ‖z‖∞ ≤ 1
+        for i in 1:NUM_ASSET
+            @addConstraint(m, z[i] <=  1)
+            @addConstraint(m, z[i] >= -1)
+        end
+        # ‖z‖₁ ≤ Γ
+        @defUnc(m, zabs[1:NUM_ASSET])
+        for i in 1:NUM_ASSET
+            @addConstraint(m, zabs[i] >=  z[i])
+            @addConstraint(m, zabs[i] >= -z[i])
+        end
+        @addConstraint(m, sum(zabs) <= Γ)
+    end
 
-    # The objective constraint - uncertain
-    @addConstraint(m, obj - dot(r, x) <= 0)
+    # The objective function: the worst-case return
+    @addConstraint(m, obj <= dot(r, x))
 
-    # Solve it, report statistics on number of cutting planes etc.
-    #println("Solving model for Γ=$Γ")
-    #solveRobust(m, report=true)
-    #println("Objective value: ", getValue(obj))
-    solveRobust(m)
+    # Solve it
+    solve(m)
 
+    # Return the allocation and the worst-case return
     return getValue(x), getValue(obj)
 end
 
@@ -95,27 +122,27 @@ end
 # simulate
 # Generate some past returns, then obtain portfolios for different values
 # of Γ to see the effect of robustness on the distribution of returns
-function simulate(num_past, num_future)
+function simulate()
     # Generate the simulated data
-    past_returns   = generate_data(num_past)
-    future_returns = generate_data(num_future)
+    past_returns   = generate_data()
+    future_returns = generate_data()
 
     # Print table header
-    println("   Γ|   Min|   10%|   20%|  Mean|   Max")
-    for Γ in [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
+    println("   Γ |    Min |    10% |    20% |   Mean |    Max")
+    for Γ in 0:NUM_ASSET
         # Solve for portfolio
         x, obj = solve_portfolio(past_returns, Γ)
         # Evaluate portfolio returns in future
-        future_z = future_returns*x[:]
+        future_z = future_returns*x
         sort!(future_z)
         min_z    = future_z[1]*100
-        ten_z    = future_z[int(num_future*0.1)]*100
-        twenty_z = future_z[int(num_future*0.2)]*100
+        ten_z    = future_z[div(NUM_SAMPLES,10)]*100
+        twenty_z = future_z[div(NUM_SAMPLES, 5)]*100
         mean_z   = mean(future_z)*100
         max_z    = future_z[end]*100
-        @printf(" %3.1f| %5.1f| %5.1f| %5.1f| %5.1f| %5.1f\n",
+        @printf(" %3.1f | %6.2f | %6.2f | %6.2f | %6.2f | %6.2f\n",
                     Γ, min_z, ten_z, twenty_z, mean_z, max_z)
     end
 end
 
-simulate(1000, 1000)
+simulate()

src/JuMPeR.jl

@@ -1,20 +1,29 @@
-#############################################################################
-# JuMPeR
-# Julia for Mathematical Programming - extension for Robust Optimization
-# See http://github.com/IainNZ/JuMPeR.jl
-#############################################################################
+#-----------------------------------------------------------------------
+# JuMPeR  --  JuMP Extension for Robust Optimization
+# http://github.com/IainNZ/JuMPeR.jl
+#-----------------------------------------------------------------------
+# Copyright (c) 2015: Iain Dunning
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at http://mozilla.org/MPL/2.0/.
+#-----------------------------------------------------------------------
 
 module JuMPeR
 
 using Compat
 
+importall Base.Operators
+
 # Import everything we need from JuMP, so we can build on it
 importall JuMP
-import JuMP.GenericAffExpr, JuMP.JuMPConstraint, JuMP.GenericRangeConstraint
-import JuMP.sense, JuMP.rhs
+import JuMP: GenericAffExpr, GenericRangeConstraint
+import JuMP: sense, rhs
 import JuMP.IndexedVector, JuMP.addelt!, JuMP.isexpr
-import JuMP.JuMPContainer, JuMP.JuMPDict, JuMP.JuMPArray
+import JuMP: JuMPContainer, JuMPDict, JuMPArray
 import JuMP.@gendict
+import JuMP: assert_validmodel, validmodel, esc_nonconstant
+import JuMP: getloopedcode, buildrefsets, getname
+import JuMP: addVectorizedConstraint
 
 # JuMPeRs exported interface
 export RobustModel, getNumUncs
@@ -92,7 +101,7 @@ getNumUncs(m::Model) = getRobust(m).numUncs
 #############################################################################
 # Uncertain
 # Similar to JuMP.Variable, has an reference back to the model and an id num
-type Uncertain
+type Uncertain <: JuMP.AbstractJuMPScalar
     m::Model
     unc::Int
 end
@@ -107,8 +116,15 @@ function Uncertain(m::Model, lower::Number, upper::Number, cat::Symbol, name::St
 end
 Uncertain(m::Model, lower::Number, upper::Number, cat::Symbol) = Uncertain(m,lower,upper,cat,"")
 getName(u::Uncertain) = unc_str(REPLMode, u.m, u.unc)
+Base.zero(::Type{Uncertain}) = UAffExpr()
+Base.zero(::Uncertain) = zero(Uncertain)
+Base.one(::Type{Uncertain}) = UAffExpr(1)
+Base.one(::Uncertain) = one(Uncertain)
 Base.isequal(u1::Uncertain, u2::Uncertain) = isequal(u1.unc, u2.unc)
-Base.promote_rule{T<:Real}(::Type{Uncertain},::Type{T}) = UAffExpr
+# Generic matrix operators in Base expect to be able to find
+# a type to use for the result, but that type needs to have a zero
+# so we can't leave it as Any
+#Base.promote_rule{T<:Real}(::Type{Uncertain},::Type{T}) = UAffExpr
 
 #############################################################################
 # Uncertain Affine Expression class
@@ -117,11 +133,7 @@ typealias UAffExpr GenericAffExpr{Float64,Uncertain}
 UAffExpr() = UAffExpr(Uncertain[],Float64[],0.)
 UAffExpr(c::Real) = UAffExpr(Uncertain[],Float64[],float(c))
 UAffExpr(u::Uncertain) = UAffExpr([u],[1.],0.)
-UAffExpr(u::Uncertain, c::Real) = UAffExpr([u],[float(c)],0.)
-# Next is a bit weird - its basically the vectorized version of UAffExpr(c) 
-UAffExpr(coeffs::Array{Float64,1}) = [UAffExpr(c) for c in coeffs]
-Base.zero(a::Type{UAffExpr}) = UAffExpr()  # For zeros(UAffExpr, dims...)
-Base.zero(a::UAffExpr) = zero(typeof(a))
+UAffExpr(c::Real,u::Uncertain,constant=0) = UAffExpr([u],[float(c)],constant)
 Base.convert(::Type{UAffExpr}, u::Uncertain) = UAffExpr(u)
 Base.convert(::Type{UAffExpr}, c::Number) = UAffExpr(c)
 
@@ -136,6 +148,8 @@ FullAffExpr() = FullAffExpr(Variable[], UAffExpr[], UAffExpr())
 Base.zero(a::Type{FullAffExpr}) = FullAffExpr()
 Base.zero(a::FullAffExpr) = zero(typeof(a))
 Base.convert(::Type{FullAffExpr}, x::Variable) = FullAffExpr([x],[UAffExpr(1)], UAffExpr())
+Base.convert(::Type{FullAffExpr}, aff::AffExpr) =
+    FullAffExpr(aff.vars,map(UAffExpr,aff.coeffs), UAffExpr(aff.constant))
 function Base.push!(faff::FullAffExpr, new_coeff::UAffExpr, new_var::Variable)
     push!(faff.vars, new_var)
     push!(faff.coeffs, new_coeff)
@@ -149,6 +163,9 @@ end
 # UncSetConstraint      Just uncertainties
 typealias UncSetConstraint GenericRangeConstraint{UAffExpr}
 addConstraint(m::Model, c::UncSetConstraint) = push!(getRobust(m).uncertaintyset, c)
+addConstraint(m::Model, c::Array{UncSetConstraint}) =
+    error("The operators <=, >=, and == can only be used to specify scalar constraints. If you are trying to add a vectorized constraint, use the element-wise dot comparison operators (.<=, .>=, or .==) instead")
+addVectorizedConstraint(m::Model, v::Array{UncSetConstraint}) = map(c->addConstraint(m,c), v)
 
 # UncConstraint         Mix of variables and uncertains
 typealias UncConstraint GenericRangeConstraint{FullAffExpr}
@@ -167,6 +184,9 @@ function addConstraint(m::Model, c::UncConstraint, w=nothing)
     push!(rd.activecuts, Any[])
     return ConstraintRef{UncConstraint}(m,length(rd.uncertainconstr))
 end
+addConstraint(m::Model, c::Array{UncConstraint}) =
+    error("The operators <=, >=, and == can only be used to specify scalar constraints. If you are trying to add a vectorized constraint, use the element-wise dot comparison operators (.<=, .>=, or .==) instead")
+addVectorizedConstraint(m::Model, v::Array{UncConstraint}) = map(c->addConstraint(m,c), v)
 
 
 #############################################################################

src/oracle.jl

@@ -93,8 +93,8 @@ function build_certain_constraint(  master::Model,
                                         coeff.coeffs[unc_ind]
         end
 
-    return sense(unc_con) == :(<=) ? new_lhs <= unc_con.ub :
-                                     new_lhs >= unc_con.lb
+    return sense(unc_con) == :(<=) ? @LinearConstraint(new_lhs <= unc_con.ub) :
+                                     @LinearConstraint(new_lhs >= unc_con.lb)
 end
 
 function build_certain_constraint(  master::Model,

src/oracle_gen_graph.jl

@@ -83,7 +83,7 @@ function setup(gen::GeneralGraphOracle, rm::Model, prefs)
             m.colNames = rd.uncNames[comp_to_unc[comp]]
             m.colLower = rd.uncLower[comp_to_unc[comp]]
             m.colUpper = rd.uncUpper[comp_to_unc[comp]]
-            m.colCat   = rd.uncCat  [comp_to_unc[comp]]
+            m.colCat   = rd.uncCat[  comp_to_unc[comp]]
             cut_vars   = [Variable(m, i) for i = 1:length(comp_to_unc[comp])]
             # Polyhedral constraints only right now
             for con_ind in 1:length(rd.uncertaintyset)
@@ -330,8 +330,8 @@ function generateCut(gen::GeneralGraphOracle, master::Model, rm::Model, inds::Ve
         end
 
         @assert comp != 0
-        cut_model       = gen.comp_cut_model[comp]
-        cut_vars        = gen.comp_cut_vars [comp]
+        cut_model       = gen.comp_cut_model[ comp]
+        cut_vars        = gen.comp_cut_vars[  comp]
         unc_to_comp_unc = gen.unc_to_comp_unc[comp]
 
         # Update the cutting plane problem's objective, and solve

src/print.jl

@@ -117,7 +117,7 @@ function fill_unc_names(mode, uncNames, u::JuMPArray{Uncertain})
     for (ind,unc) in enumerate(u.innerArray)
         idx_strs = [string( idxsets[1][mod1(ind,lengths[1])] )]
         for i = 2:N
-            push!(idx_strs, string(idxsets[i][int(ceil(mod1(ind,cprod[i]) / cprod[i-1]))]))
+            push!(idx_strs, string(idxsets[i][@compat Int(ceil(mod1(ind,cprod[i]) / cprod[i-1]))]))
         end
         #if mode == IJuliaMode
         #    uncNames[unc.unc] = string(name, "_{", join(idx_strs,",") , "}")