Merge pull request #22 from chkwon/rdeits-lcp

appveyor for julia 0.6

appveyor.yml

@@ -2,6 +2,8 @@ environment:
   matrix:
   - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.5/julia-0.5-latest-win32.exe"
   - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.5/julia-0.5-latest-win64.exe"
+  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.6/julia-0.6-latest-win32.exe"
+  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.6/julia-0.6-latest-win64.exe"
   - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
   - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"
 

src/PATHSolver.jl

@@ -10,14 +10,29 @@ else
   error("PATHSolver not properly installed. Please run Pkg.build(\"PATHSolver\")")
 end
 
-export solveMCP, options
+export solveMCP, solveLCP, options
 
 # Global function pointers for the user-supplied function and jacobian evaluators.
 const user_f = Ref(FunctionWrapper{Vector{Cdouble}, Tuple{Vector{Cdouble}}}(identity))
 # The annotated SparseMatrixCSC return type will automatically convert the 
 # jacobian into the correct sparse form for PATH
 const user_j = Ref(FunctionWrapper{SparseMatrixCSC{Cdouble, Cint}, Tuple{Vector{Cdouble}}}(identity))
 
+const cached_J = [convert(SparseMatrixCSC{Cdouble, Cint}, zeros(0, 0))]
+const cached_J_filled = Ref(false)
+
+const status = 
+   [  :Solved,                          # 1 - solved
+      :StationaryPointFound,            # 2 - stationary point found
+      :MajorIterationLimit,             # 3 - major iteration limit
+      :CumulativeMinorIterationLimit,   # 4 - cumulative minor iteration limit
+      :TimeLimit,                       # 5 - time limit
+      :UserInterrupt,                   # 6 - user interrupt
+      :BoundError,                      # 7 - bound error (lb is not less than ub)
+      :DomainError,                     # 8 - domain error (could not find a starting point)
+      :InternalError                    # 9 - internal error
+  ]
+
 count_nonzeros(M::AbstractSparseMatrix) = nnz(M)
 count_nonzeros(M::AbstractMatrix) = count(x -> x != 0, M) # fallback for dense matrices
 
@@ -47,24 +62,40 @@ function solveMCP(f_eval::Function, j_eval::Function, lb::Vector, ub::Vector, va
            Ptr{Void}, Ptr{Void}),
            n, nnz, z, f, lb, ub, var_name, con_name, f_user_cb, j_user_cb)
 
-  status =
-   [  :Solved,                          # 1 - solved
-      :StationaryPointFound,            # 2 - stationary point found
-      :MajorIterationLimit,             # 3 - major iteration limit
-      :CumulativeMinorIterationLimit,   # 4 - cumulative minor iteration limit
-      :TimeLimit,                       # 5 - time limit
-      :UserInterrupt,                   # 6 - user interrupt
-      :BoundError,                      # 7 - bound error (lb is not less than ub)
-      :DomainError,                     # 8 - domain error (could not find a starting point)
-      :InternalError                    # 9 - internal error
-  ]
-
   remove_option_file()
   return status[t], z, f
+end
 
+function solveLCP(f_eval::Function, lb::AbstractVector, ub::AbstractVector, 
+                  var_name=C_NULL, con_name=C_NULL)
+  J = ForwardDiff.jacobian(f_eval, lb)
+  solveLCP(f_eval, J, lb, ub, var_name, con_name)
 end
 
+function solveLCP(f_eval::Function, M::AbstractMatrix, 
+                  lb::AbstractVector, ub::AbstractVector,
+                  var_name=C_NULL, con_name=C_NULL)
+  user_f[] = f_eval
+  cached_J[] = M
+  cached_J_filled[] = false
+  f_user_cb = cfunction(f_user_wrap, Cint, (Cint, Ptr{Cdouble}, Ptr{Cdouble}))
+  j_user_cb = cfunction(cached_j_user_wrap, Cint, (Cint, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}))
 
+  n = length(lb)
+  z = copy(lb)
+  f = zeros(n)
+
+  nnz = count_nonzeros(M)
+  t = ccall( (:path_main, "libpath47julia"), Cint,
+          (Cint, Cint,
+           Ptr{Cdouble}, Ptr{Cdouble},
+           Ptr{Cdouble}, Ptr{Cdouble},
+           Ptr{Ptr{Cchar}}, Ptr{Ptr{Cchar}},
+           Ptr{Void}, Ptr{Void}),
+           n, nnz, z, f, lb, ub, var_name, con_name, f_user_cb, j_user_cb)
+  remove_option_file()
+  return status[t], z, f
+end
 
 function remove_option_file()
   if isfile("path.opt")
@@ -104,7 +135,23 @@ function j_user_wrap(n::Cint, expected_nnz::Cint, z_ptr::Ptr{Cdouble},
   if nnz(J) > expected_nnz
     error("Evaluated jacobian has more nonzero entries than were initially provided in solveMCP()")
   end
+  load_sparse_matrix(J, n, expected_nnz, col_start_ptr, col_len_ptr, row_ptr, data_ptr)
+  return Cint(0)
+end
+
+function cached_j_user_wrap(n::Cint, expected_nnz::Cint, z_ptr::Ptr{Cdouble},
+                     col_start_ptr::Ptr{Cint}, col_len_ptr::Ptr{Cint}, 
+                     row_ptr::Ptr{Cint}, data_ptr::Ptr{Cdouble})
+  if !(cached_J_filled[])
+    load_sparse_matrix(cached_J[], n, expected_nnz, col_start_ptr, col_len_ptr, row_ptr, data_ptr)
+    cached_J_filled[] = true
+  end
+  return Cint(0)
+end
 
+function load_sparse_matrix(J::SparseMatrixCSC, n::Cint, expected_nnz::Cint,
+                            col_start_ptr::Ptr{Cint}, col_len_ptr::Ptr{Cint}, 
+                            row_ptr::Ptr{Cint}, data_ptr::Ptr{Cdouble})
   # Transfer data from the computed jacobian into the sparse format that PATH
   # expects. Fortunately, PATH uses a compressed-sparse-column storage which
   # is compatible with Julia's default SparseMatrixCSC format.
@@ -130,7 +177,8 @@ function j_user_wrap(n::Cint, expected_nnz::Cint, z_ptr::Ptr{Cdouble},
     row[i] = rv[i]
     data[i] = nz[i]
   end
-  return Cint(0)
 end
 
+
+
 end # Module

test/lcp.jl

@@ -0,0 +1,99 @@
+@testset "linear complementarity problem" begin
+    M = [0  0 -1 -1 ;
+         0  0  1 -2 ;
+         1 -1  2 -2 ;
+         1  2 -2  4 ]
+
+    q = [2; 2; -2; -6]
+
+    myfunc(x) = M*x + q
+
+    n = 4
+    lb = zeros(n)
+    ub = 100*ones(n)
+
+    status, z, f = solveLCP(myfunc, M, lb, ub)
+    @show status
+    @show z
+    @show f
+    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
+    @test status == :Solved
+
+    status, z, f = solveLCP(myfunc, lb, ub)
+    @show status
+    @show z
+    @show f
+    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
+    @test status == :Solved
+
+
+    println("-------------------------------------------------------")
+
+
+    M = [0  0 -1 -1 ;
+         0  0  1 -2 ;
+         1 -1  2 -2 ;
+         1  2 -2  4 ]
+
+    q = [2; 2; -2; -6]
+
+    myfunc(x) = M*x + q
+
+    n = 4
+    lb = zeros(n)
+    ub = 100*ones(n)
+
+    options(convergence_tolerance=1e-2, output=:no, time_limit=3600)
+
+    status, z, f = solveLCP(myfunc, M, lb, ub)
+    @show status
+    @show z
+    @show f
+    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
+    @test status == :Solved
+
+    status, z, f = solveLCP(myfunc, lb, ub)
+    @show status
+    @show z
+    @show f
+    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
+    @test status == :Solved
+
+
+
+    println("-------------------------------------------------------")
+
+
+    function elemfunc(x)
+        val = similar(x)
+        val[1] = -x[3]-x[4] + q[1]
+        val[2] = x[3] -2x[4] + q[2]
+        val[3] = x[1]-x[2]+2x[3]-2x[4] + q[3]
+        val[4] = x[1]+2x[2]-2x[3]+4x[4] + q[4]
+        return val
+    end
+
+    n = 4
+    lb = zeros(n)
+    ub = 100*ones(n)
+
+    var_name = ["x1", "x2", "x3", "x4"]
+    con_name = ["F1", "F2", "F3", "F4"]
+
+    options(convergence_tolerance=1e-2, output=:yes, time_limit=3600)
+
+    status, z, f = solveLCP(elemfunc, M, lb, ub)
+    status, z, f = solveLCP(elemfunc, M, lb, ub, var_name)
+    status, z, f = solveLCP(elemfunc, M, lb, ub, var_name, con_name)
+    status, z, f = solveLCP(elemfunc, lb, ub)
+    status, z, f = solveLCP(elemfunc, lb, ub, var_name)
+    status, z, f = solveLCP(elemfunc, lb, ub, var_name, con_name)
+
+    @show status
+    @show z
+    @show f
+
+
+    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
+    @test status == :Solved
+end

test/mcp.jl

@@ -0,0 +1,131 @@
+@testset "mixed complementarity problem" begin
+    M = [0  0 -1 -1 ;
+         0  0  1 -2 ;
+         1 -1  2 -2 ;
+         1  2 -2  4 ]
+
+    q = [2; 2; -2; -6]
+
+    myfunc(x) = M*x + q
+
+    n = 4
+    lb = zeros(n)
+    ub = 100*ones(n)
+
+    status, z, f = solveMCP(myfunc, lb, ub)
+    @show status
+    @show z
+    @show f
+
+
+    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
+    @test status == :Solved
+
+
+    println("-------------------------------------------------------")
+
+
+    M = [0  0 -1 -1 ;
+         0  0  1 -2 ;
+         1 -1  2 -2 ;
+         1  2 -2  4 ]
+
+    q = [2; 2; -2; -6]
+
+    myfunc(x) = M*x + q
+
+    n = 4
+    lb = zeros(n)
+    ub = 100*ones(n)
+
+    options(convergence_tolerance=1e-2, output=:no, time_limit=3600)
+
+    status, z, f = solveMCP(myfunc, lb, ub)
+
+    @show status
+    @show z
+    @show f
+
+    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
+    @test status == :Solved
+
+
+
+    println("-------------------------------------------------------")
+
+
+    function elemfunc(x)
+        val = similar(x)
+        val[1] = -x[3]-x[4] + q[1]
+        val[2] = x[3] -2x[4] + q[2]
+        val[3] = x[1]-x[2]+2x[3]-2x[4] + q[3]
+        val[4] = x[1]+2x[2]-2x[3]+4x[4] + q[4]
+        return val
+    end
+
+    n = 4
+    lb = zeros(n)
+    ub = 100*ones(n)
+
+    var_name = ["x1", "x2", "x3", "x4"]
+    con_name = ["F1", "F2", "F3", "F4"]
+
+    options(convergence_tolerance=1e-2, output=:yes, time_limit=3600)
+
+
+    jacfunc(x) = M
+
+    status, z, f = solveMCP(elemfunc, lb, ub)
+    status, z, f = solveMCP(elemfunc, lb, ub, var_name)
+    status, z, f = solveMCP(elemfunc, lb, ub, var_name, con_name)
+    status, z, f = solveMCP(elemfunc, jacfunc, lb, ub)
+    status, z, f = solveMCP(elemfunc, jacfunc, lb, ub, var_name)
+    status, z, f = solveMCP(elemfunc, jacfunc, lb, ub, var_name, con_name)
+
+
+    @show status
+    @show z
+    @show f
+
+
+    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
+    @test status == :Solved
+
+    function test_in_local_scope()
+        # Verify that we can solve MCPs in local scope. Surprisingly, this is 
+        # is relevant because it affects the way closures are generated. To be 
+        # specific, you can do the following in global scope:
+        # 
+        # julia> y = [1]
+        # julia> cfunction(x -> x + y[1], Int, (Int,))
+        # 
+        # but running the same code inside a function will fail with:
+        #   ERROR: closures are not yet c-callable
+
+        M = [0  0 -1 -1 ;
+             0  0  1 -2 ;
+             1 -1  2 -2 ;
+             1  2 -2  4 ]
+
+        q = [2; 2; -2; -6]
+
+        myfunc(x) = M*x + q
+
+        n = 4
+        lb = zeros(n)
+        ub = 100*ones(n)
+
+        options(convergence_tolerance=1e-2, output=:no, time_limit=3600)
+
+        status, z, f = solveMCP(myfunc, lb, ub)
+
+        @show status
+        @show z
+        @show f
+
+        @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
+        @test status == :Solved
+    end
+
+    test_in_local_scope()
+end