Cleanup and performance improvements

.travis.yml

@@ -6,6 +6,7 @@ os:
     - osx
 julia:
     - 0.5
+    - 0.6
     - nightly
 notifications:
     email: false

REQUIRE

@@ -1,3 +1,5 @@
 julia 0.5
 ForwardDiff 0.2
 BinDeps
+FunctionWrappers 0.1
+Compat 0.40

src/PATHSolver.jl

@@ -1,54 +1,34 @@
 module PATHSolver
 
 using ForwardDiff
+using FunctionWrappers: FunctionWrapper
+using Compat: pairs
 
 if isfile(joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl"))
   include(joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl"))
 else
   error("PATHSolver not properly installed. Please run Pkg.build(\"PATHSolver\")")
 end
 
-
-
-
 export solveMCP, 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))
 
+count_nonzeros(M::AbstractSparseMatrix) = nnz(M)
+count_nonzeros(M::AbstractMatrix) = count(x -> x != 0, M) # fallback for dense matrices
 
-
-function solveMCP(f_eval::Function, lb::Vector, ub::Vector)
-  var_name = C_NULL
-  con_name = C_NULL
-  return solveMCP(f_eval, lb, ub, var_name, con_name)
-end
-
-function solveMCP(f_eval::Function, lb::Vector, ub::Vector, var_name)
-  con_name = C_NULL
-  return solveMCP(f_eval, lb, ub, var_name, con_name)
-end
-
-function solveMCP(f_eval::Function, lb::Vector, ub::Vector, var_name, con_name)
+function solveMCP(f_eval::Function, lb::Vector, ub::Vector, var_name=C_NULL, con_name=C_NULL)
   j_eval = x -> ForwardDiff.jacobian(f_eval, x)
   return solveMCP(f_eval, j_eval, lb, ub, var_name, con_name)
 end
 
-function solveMCP(f_eval::Function, j_eval::Function, lb::Vector, ub::Vector)
-  var_name = C_NULL
-  con_name = C_NULL
-  return solveMCP(f_eval, j_eval, lb, ub, var_name, con_name)
-end
-
-function solveMCP(f_eval::Function, j_eval::Function, lb::Vector, ub::Vector, var_name)
-  con_name = C_NULL
-  return solveMCP(f_eval, j_eval, lb, ub, var_name, con_name)
-end
-
-
-
-function solveMCP(f_eval::Function, j_eval::Function, lb::Vector, ub::Vector, var_name, con_name)
-  global user_f = f_eval
-  global user_j = j_eval
-
+function solveMCP(f_eval::Function, j_eval::Function, lb::Vector, ub::Vector, var_name=C_NULL, con_name=C_NULL)
+  user_f[] = f_eval
+  user_j[] = j_eval
   f_user_cb = cfunction(f_user_wrap, Cint, (Cint, Ptr{Cdouble}, Ptr{Cdouble}))
   j_user_cb = cfunction(j_user_wrap, Cint, (Cint, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}))
 
@@ -57,8 +37,7 @@ function solveMCP(f_eval::Function, j_eval::Function, lb::Vector, ub::Vector, va
   f = zeros(n)
 
   J0 = j_eval(z)
-  s_col, s_len, s_row, s_data = sparse_matrix(J0)
-  nnz = length(s_data)
+  nnz = count_nonzeros(J0)
 
   t = ccall( (:path_main, "libpath47julia"), Cint,
           (Cint, Cint,
@@ -96,7 +75,7 @@ end
 function options(;kwargs...)
   opt_file = open("path.opt", "w")
   println(opt_file, "* Generated by PATHSolver.jl. Do not edit.")
-  for (key, value) in kwargs
+  for (key, value) in pairs(kwargs)
     println(opt_file, key, " ", value)
   end
   close(opt_file)
@@ -109,112 +88,49 @@ end
 # static int (*f_eval)(int n, double *z, double *f);
 # static int (*j_eval)(int n, int nnz, double *z, int *col_start, int *col_len,
       # int *row, double *data);
-
-function f_user_wrap(n::Cint, z::Ptr{Cdouble}, f::Ptr{Cdouble})
-  F = user_f(unsafe_wrap(Array{Float64}, z, Int(n), false))
-  unsafe_store_vector!(f, F)
-  return Cint(0)
-end
-
-function j_user_wrap(n::Cint, nnz::Cint, z::Ptr{Cdouble},
-  col_start::Ptr{Cint}, col_len::Ptr{Cint}, row::Ptr{Cint}, data::Ptr{Cdouble})
-
-  J = user_j(unsafe_wrap(Array{Float64}, z, Int(n), false))
-
-  s_col, s_len, s_row, s_data = sparse_matrix(J)
-
-  unsafe_store_vector!(col_start, s_col)
-  unsafe_store_vector!(col_len, s_len)
-  unsafe_store_vector!(row, s_row)
-  unsafe_store_vector!(data, s_data)
-
+function f_user_wrap(n::Cint, z_ptr::Ptr{Cdouble}, f_ptr::Ptr{Cdouble})
+  z = unsafe_wrap(Array{Cdouble}, z_ptr, Int(n), false)
+  f = unsafe_wrap(Array{Cdouble}, f_ptr, Int(n), false)
+  f .= user_f[](z)
   return Cint(0)
 end
-###############################################################################
-
-
 
+function 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})
 
-###############################################################################
-# Converting the Jacobian matrix to the sparse matrix format of the PATH Solver
-###############################################################################
-function sparse_matrix(A::AbstractSparseArray)
-  m, n = size(A)
-  @assert m==n
-
-  col_start = Array{Int}(n)
-  col_len = Array{Int}(n)
-  row = Array{Int}(0)
-  data = Array{Float64}(0)
-  for j in 1:n
-    if j==1
-      col_start[j] = 1
-    else
-      col_start[j] = col_start[j-1] + col_len[j-1]
-    end
-
-    col_len[j] = 0
-    for i in 1:n
-      if A[i,j] != 0.0
-        col_len[j] += 1
-        push!(row, i)
-        push!(data, A[i,j])
-      end
-    end
+  z = unsafe_wrap(Array{Cdouble}, z_ptr, Int(n), false)
+  J::SparseMatrixCSC{Cdouble, Cint} = user_j[](z)
+  if nnz(J) > expected_nnz
+    error("Evaluated jacobian has more nonzero entries than were initially provided in solveMCP()")
   end
 
-  return col_start, col_len, row, data
-end
-
-function sparse_matrix(A::Matrix)
-  return sparse_matrix(sparse(A))
-  # m, n = size(A)
-  # @assert m==n
-  #
-  # col_start = Array{Int}(n)
-  # col_len = Array{Int}(n)
-  # row = Array{Int}(0)
-  # data = Array{Float64}(0)
-  # for j in 1:n
-  #   if j==1
-  #     col_start[j] = 1
-  #   else
-  #     col_start[j] = col_start[j-1] + col_len[j-1]
-  #   end
-  #
-  #   col_len[j] = 0
-  #   for i in 1:n
-  #     if A[i,j] != 0.0
-  #       col_len[j] += 1
-  #       push!(row, i)
-  #       push!(data, A[i,j])
-  #     end
-  #   end
-  # end
-  #
-  # return col_start, col_len, row, data
-end
-###############################################################################
-
-
-
-
-
-
-function unsafe_store_vector!(x_ptr::Ptr{Cint}, x_val::Vector)
-  for i in 1:length(x_val)
-    unsafe_store!(x_ptr, x_val[i], i)
+  # 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.
+
+  # col_start in PATH corresponds to J.colptr[1:end-1]
+  col_start = unsafe_wrap(Array{Cint}, col_start_ptr, Int(n), false)
+  # col_len in PATH corresponds to diff(J.colptr)
+  col_len = unsafe_wrap(Array{Cint}, col_len_ptr, Int(n), false)
+  # row in PATH corresponds to rowvals(J)
+  row = unsafe_wrap(Array{Cint}, row_ptr, Int(expected_nnz), false)
+  # data in PATH corresponds to nonzeros(J)
+  data = unsafe_wrap(Array{Cdouble}, data_ptr, Int(expected_nnz), false)
+
+  @inbounds for i in 1:n
+    col_start[i] = J.colptr[i]
+    col_len[i] = J.colptr[i + 1] - J.colptr[i]
   end
-  return
-end
 
-function unsafe_store_vector!(x_ptr::Ptr{Cdouble}, x_val::Vector)
-  for i in 1:length(x_val)
-    unsafe_store!(x_ptr, x_val[i], i)
+  rv = rowvals(J)
+  nz = nonzeros(J)
+  num_nonzeros = nnz(J)
+  @inbounds for i in 1:num_nonzeros
+    row[i] = rv[i]
+    data[i] = nz[i]
   end
-  return
+  return Cint(0)
 end
 
-
-
 end # Module

test/runtests.jl

@@ -2,6 +2,8 @@ using PATHSolver
 using ForwardDiff
 using Base.Test
 
+include("sparse_matrix.jl")
+
 M = [0  0 -1 -1 ;
      0  0  1 -2 ;
      1 -1  2 -2 ;
@@ -93,3 +95,41 @@ status, z, f = solveMCP(elemfunc, jacfunc, lb, ub, var_name, con_name)
 
 @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()

test/sparse_matrix.jl

@@ -0,0 +1,60 @@
+using PATHSolver
+using Base.Test
+
+
+# Verify that pulling data directly out of the SparseMatrixCSC form gives
+# exactly what the previous hand-coded algorithm gave.
+@testset "sparse matrix" begin
+  function sparse_matrix_reference(A::AbstractSparseMatrix)
+    m, n = size(A)
+    @assert m==n
+
+    col_start = Array{Int}(n)
+    col_len = Array{Int}(n)
+    row = Array{Int}(0)
+    data = Array{Float64}(0)
+    for j in 1:n
+      if j==1
+        col_start[j] = 1
+      else
+        col_start[j] = col_start[j-1] + col_len[j-1]
+      end
+
+      col_len[j] = 0
+      for i in 1:n
+        if A[i,j] != 0.0
+          col_len[j] += 1
+          push!(row, i)
+          push!(data, A[i,j])
+        end
+      end
+    end
+    return col_start, col_len, row, data
+  end
+
+  sparse_matrix_reference(A::Matrix) = sparse_matrix_reference(sparse(A))
+
+  function sparse_matrix_csc(A::SparseMatrixCSC)
+    m, n = size(A)
+    @assert m==n
+    col_start = A.colptr[1:end-1]
+    col_len = diff(A.colptr)
+    row = rowvals(A)
+    data = nonzeros(A)
+    return col_start, col_len, row, data
+  end
+
+  sparse_matrix_csc(A::Matrix) = sparse_matrix_csc(convert(SparseMatrixCSC, A))
+
+  srand(42)
+  for i in 1:100
+    n = rand(5:20)
+    M = sprandn(n, n, 0.1)
+    @test sparse_matrix_csc(M) == sparse_matrix_reference(M)
+
+    Mf = full(M)
+    @test sparse_matrix_csc(Mf) == sparse_matrix_reference(Mf)
+    @test sparse_matrix_csc(Mf) == sparse_matrix_csc(M)
+  end
+end
+