JuliaCon2016-EmbeddingReverseCommunicationSolvers.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Enabling and embedding reverse communication solvers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia Computing regularly assists clients in developing custom Julia functionality tailored to their specific application and operational needs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A recent client engagement involved the following requirements:\n",
    "\n",
    "* Integrate Julia in a large C/C++ batch processing workflow\n",
    "* Call Julia as a library from C\n",
    "* Provide Julian alternatives to existing optimization and root finding solvers\n",
    "* Wrap client's existing solvers with a Julia interface\n",
    "* Provide a reverse communication interface to Julia solvers, including the wrapped solvers\n",
    "  * Objective & derivative functions cannot be passed as arguments to Julia solver (no wrapping of cost function)\n",
    "  * Objective & derivative functions should be executed in C/C++ \n",
    "* Enable future infrastructure for easy migration of algorithms from front office to back office without recoding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is the difference between forward and reverse communication solvers?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forward communication solvers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Forward communication is the strategy employed by most Julia optimization and root finding packages (Roots.jl, Optim.jl, NLsolve.jl, NLopt.jl, etc.)\n",
    "\n",
    "In each of these packages, an objective function is passed as one of the arguments to a root finder or optimizer.  \n",
    "\n",
    "The solver directly invokes the passed in objective function as necessary.\n",
    "\n",
    "Use of forward communication can allow for incorporating powerful techniques, such as automatic differentiation.\n",
    "\n",
    "Example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Optim, NLsolve\n",
    "\n",
    "###################################################\n",
    "# Rosenbrock objective and derivative functions defined for use by Optim.jl\n",
    "\n",
    "function rosenbrock_optim(x::Vector)\n",
    "    (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2\n",
    "end \n",
    "\n",
    "function rosenbrock_optim_gradient!(x::Vector, storage::Vector)\n",
    "    storage[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]\n",
    "    storage[2] = 200.0 * (x[2] - x[1]^2)\n",
    "end\n",
    "\n",
    "function rosenbrock_optim_hessian!(x::Vector, storage::Matrix)\n",
    "    storage[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2\n",
    "    storage[1, 2] = -400.0 * x[1]\n",
    "    storage[2, 1] = -400.0 * x[1]\n",
    "    storage[2, 2] = 200.0\n",
    "end\n",
    "\n",
    "f = rosenbrock_optim\n",
    "g = rosenbrock_optim_gradient!\n",
    "h = rosenbrock_optim_hessian!\n",
    "\n",
    "####################################################\n",
    "# Rosenbrock objective and derivative functions defined for use by NLsolve.jl\n",
    "\n",
    "function rosenbrock_nlsolve_f!(x::Vector, fvec::Vector)\n",
    "    fvec[1] = (1 - x[1])^2\n",
    "    fvec[2] = 100.0 * (x[2] - x[1]^2)^2\n",
    "end\n",
    "\n",
    "function rosenbrock_nlsolve_g!(x::Vector, fjac::Matrix)\n",
    "    fjac[1,1] = -2.0 * (1 - x[1])\n",
    "    fjac[1,2] = 0.0\n",
    "    fjac[2,1] = -400.0 * (x[2] - x[1]^2) * x[1]\n",
    "    fjac[2,2] = 200.0 * (x[2] - x[1]^2)\n",
    "end\n",
    "\n",
    "function rosenbrock_nlsolve_fg!(x::Vector, fvec::Vector, fjac::Matrix)\n",
    "    fvec[1]   = (1.0 - x[1])^2\n",
    "    fvec[2]   = 100.0 * (x[2] - x[1]^2)^2\n",
    "    fjac[1,1] = -2.0 * (1 - x[1])\n",
    "    fjac[1,2] = 0.0\n",
    "    fjac[2,1] = -400.0 * (x[2] - x[1]^2) * x[1]\n",
    "    fjac[2,2] = 200.0 * (x[2] - x[1]^2)\n",
    "end\n",
    "\n",
    "f! = rosenbrock_nlsolve_f!\n",
    "g! = rosenbrock_nlsolve_g!\n",
    "fg! = rosenbrock_nlsolve_fg!\n",
    "\n",
    "df = NLsolve.DifferentiableMultivariateFunction(f!, g!, fg!)\n",
    "\n",
    ";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Passing cost/derivative functions as an input arguments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "results = Optim.optimize(f,[-1.2,1.0],LBFGS()) = Results of Optimization Algorithm\n",
      " * Algorithm: L-BFGS\n",
      " * Starting Point: [-1.2,1.0]\n",
      " * Minimizer: [0.9999999926455753,0.999999985283916]\n",
      " * Minimum: 0.000000\n",
      " * Iterations: 30\n",
      " * Convergence: true\n",
      "   * |x - x'| < 1.0e-32: false\n",
      "   * |f(x) - f(x')| / |f(x)| < 1.0e-08: false\n",
      "   * |g(x)| < 1.0e-08: true\n",
      "   * Reached Maximum Number of Iterations: false\n",
      " * Objective Function Calls: 128\n",
      " * Gradient Calls: 128\n",
      "results = Optim.optimize(f,g,[-1.2,1.0],LBFGS()) = Results of Optimization Algorithm\n",
      " * Algorithm: L-BFGS\n",
      " * Starting Point: [-1.2,1.0]\n",
      " * Minimizer: [1.000000000022946,1.000000000043234]\n",
      " * Minimum: 0.000000\n",
      " * Iterations: 26\n",
      " * Convergence: true\n",
      "   * |x - x'| < 1.0e-32: false\n",
      "   * |f(x) - f(x')| / |f(x)| < 1.0e-08: false\n",
      "   * |g(x)| < 1.0e-08: true\n",
      "   * Reached Maximum Number of Iterations: false\n",
      " * Objective Function Calls: 115\n",
      " * Gradient Calls: 115\n",
      "results = Optim.optimize(f,g,h,[-1.2,1.0],Newton()) = Results of Optimization Algorithm\n",
      " * Algorithm: Newton's Method\n",
      " * Starting Point: [-1.2,1.0]\n",
      " * Minimizer: [0.9999999999999988,0.9999999999999962]\n",
      " * Minimum: 0.000000\n",
      " * Iterations: 23\n",
      " * Convergence: true\n",
      "   * |x - x'| < 1.0e-32: false\n",
      "   * |f(x) - f(x')| / |f(x)| < 1.0e-08: false\n",
      "   * |g(x)| < 1.0e-08: true\n",
      "   * Reached Maximum Number of Iterations: false\n",
      " * Objective Function Calls: 97\n",
      " * Gradient Calls: 97\n",
      "results = NLsolve.nlsolve(f!,[-1.2,1.0],method=:newton) = Results of Nonlinear Solver Algorithm\n",
      " * Algorithm: Newton with line-search\n",
      " * Starting Point: [-1.2,1.0]\n",
      " * Zero: [0.9999979019165038,0.9999907684719761]\n",
      " * Inf-norm of residuals: 0.000000\n",
      " * Iterations: 20\n",
      " * Convergence: true\n",
      "   * |x - x'| < 0.0e+00: false\n",
      "   * |f(x)| < 1.0e-08: true\n",
      " * Function Calls (f): 21\n",
      " * Jacobian Calls (df/dx): 20\n",
      "results = NLsolve.nlsolve(f!,g!,[-1.2,1.0],method=:trust_region) = Results of Nonlinear Solver Algorithm\n",
      " * Algorithm: Trust-region with dogleg and autoscaling\n",
      " * Starting Point: [-1.2,1.0]\n",
      " * Zero: [0.9999734745492139,0.9999379198581341]\n",
      " * Inf-norm of residuals: 0.000000\n",
      " * Iterations: 37\n",
      " * Convergence: true\n",
      "   * |x - x'| < 0.0e+00: false\n",
      "   * |f(x)| < 1.0e-08: true\n",
      " * Function Calls (f): 38\n",
      " * Jacobian Calls (df/dx): 25\n",
      "results = NLsolve.nlsolve(df,[-1.2,1.0]) = Results of Nonlinear Solver Algorithm\n",
      " * Algorithm: Trust-region with dogleg and autoscaling\n",
      " * Starting Point: [-1.2,1.0]\n",
      " * Zero: [0.9999734745492139,0.9999379198581341]\n",
      " * Inf-norm of residuals: 0.000000\n",
      " * Iterations: 37\n",
      " * Convergence: true\n",
      "   * |x - x'| < 0.0e+00: false\n",
      "   * |f(x)| < 1.0e-08: true\n",
      " * Function Calls (f): 38\n",
      " * Jacobian Calls (df/dx): 25\n"
     ]
    }
   ],
   "source": [
    "@show results = Optim.optimize(f, [-1.2, 1.0], LBFGS())\n",
    "\n",
    "@show results = Optim.optimize(f, g, [-1.2, 1.0], LBFGS())\n",
    "\n",
    "@show results = Optim.optimize(f, g, h, [-1.2, 1.0], Newton())\n",
    "\n",
    "@show results = NLsolve.nlsolve(f!, [-1.2, 1.0], method = :newton)\n",
    "\n",
    "@show results = NLsolve.nlsolve(f!, g!, [-1.2, 1.0], method = :trust_region)\n",
    "\n",
    "@show results = NLsolve.nlsolve(df, [-1.2, 1.0])\n",
    ";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objective and derivative functions can take required parameters as input arguments, or some parameters can be passed via closures.\n",
    "\n",
    "The important point is that the forward communication solver directly [invokes the cost function](https://github.com/JuliaOpt/Optim.jl/blob/b01425755c653ff8fe11cc1e3aa41b648331f6dd/src/l_bfgs.jl#L97) within its own internal logic.  \n",
    "\n",
    "Hence the objective and derivative functions need input arguments structured in a well defined manner for being called directly within the solver."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reverse communication solvers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reverse communication is a classic technique in which objective and derivative functions are not passed directly as inputs to the solver.  \n",
    "\n",
    "Solver returns to its calliing environment whenever an objective or derivative function value is required.  \n",
    "\n",
    "After execution of the objective or derivative function in the calling environment, the solver is called, advancing execution from its prior location."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Traditional methods of transforming forward communication to reverse communication"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transforming a forward communication solver into a reverse communication solver can be performed by:\n",
    "* Tracing execution of the forward solver\n",
    "* Marking the locations of each objective or derivative function evaluation\n",
    "* Defining a stateful input parameter for storing current objective and derivative function values and location within the solver\n",
    "* Adding logic to a modified solver to resume execution at desired locations\n",
    "* Calling the new reverse communication solver repeatedly until reaching an exit condition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example of trapezoidal integration by forward communication, marking locations of interest for reverse communication\n",
    "\n",
    "(Adapted from example in [\"Writing Scientific Software: A Guide to Good Style\"](http://www.cambridge.org/us/academic/subjects/computer-science/scientific-computing-scientific-software/writing-scientific-software-guide-good-style))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Trapezoidal integration via foward communication \n",
    "```julia\n",
    "function integrate(f::Function, a, b, n) # f is a function taking one scalar input\n",
    "    h = (b-a)/n\n",
    "    # Position 1\n",
    "    s = 0.5*f(a)\n",
    "    # Position 2\n",
    "    s += 0.5*f(b)\n",
    "    for i = 1:n-1\n",
    "        # Position 3\n",
    "        s += f(a+i*h)\n",
    "    end\n",
    "    # Position 4\n",
    "    val = s*h\n",
    "end\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Reverse communication construction via a mechanical transformation process\n",
    "```julia\n",
    "type InternalState\n",
    "    position::Int\n",
    "    i::Int\n",
    "    n::Int\n",
    "    h::Float64\n",
    "    v::Float64\n",
    "    flag::Int\n",
    "end\n",
    "\n",
    "# Controller flag argument\n",
    "# flag = 0   Initialize on entry,           success on exit\n",
    "# flag = 1    fx = f(x) on entry, Evaluate fx= f(x) on exit\n",
    "# flag = -1       Error on entry,           failure on exit\n",
    "\n",
    "function integrate_rc(s::InternalState, x::Ptr{Float64}, fx::Float64, a::Float64, b::Float64, n::Int)\n",
    "    val = NaN\n",
    "    if s.flag == 0\n",
    "        s.position = 1\n",
    "        s.n = n\n",
    "        x[] = a\n",
    "        s.flag = 1\n",
    "    elseif s.flag == 1\n",
    "        if s.position == 1\n",
    "            s.v = 0.5*fx\n",
    "            x[] = b\n",
    "            s.flag = 1\n",
    "            s.position = 2\n",
    "        elseif s.position == 2\n",
    "            s.v += 0.5*fx\n",
    "            s.i = 1\n",
    "            x[] = a + s.i*s.h\n",
    "            s.flag = 1\n",
    "            s.position = 3\n",
    "        elseif s.position == 3\n",
    "            s.v += fx\n",
    "            s.i += 1\n",
    "            x[] = a + s.i*s.h\n",
    "            s.flag = 1\n",
    "            if s.i == n-1\n",
    "                s.position = 4\n",
    "            end\n",
    "        elseif s.position == 4\n",
    "            s.v += fx\n",
    "            val = s.v*s.h\n",
    "            s.flag = 0\n",
    "        else\n",
    "            s.flag = -1\n",
    "        end\n",
    "        s.flag = -1\n",
    "    end\n",
    "    return val\n",
    "end\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Performing this type of transformation can be a complicated process.\n",
    "\n",
    "Especially for solvers that use multiple layers of different functions, files and dependent packages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using Julia tasks to transform an existing forward communication solver"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using Julia coroutines and task scheduling, a forward solver executing as a task can be paused to obtain objective or derivative values in a separate task.\n",
    "\n",
    "Upon obtaining the desired function value, a task switch can resume the original solver execution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Structuring this transformation incorporates the following components:\n",
    "\n",
    "* A task for the main execution environment\n",
    "  * Objective and derivative functions are evaluated in the scope of this main task\n",
    "* A task for the solver\n",
    "* Stub function(s) as input arguments to the forward communication solver\n",
    "  * Each stub function:\n",
    "    * Defines a value, or a tuple of values, for passing to the main calling environment\n",
    "    * Performs a task switch to the main calling environment\n",
    "* A driver while loop\n",
    "  * Controls the solver and function execution\n",
    "  * Checks the state of the solver task\n",
    "  * Contains logic for execution of the actual objective and derivative functions\n",
    "  * Performs a task switch to the solver, resuming execution within a stub function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reverse communication example using only an objective function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "result = Results of Optimization Algorithm\n",
      " * Algorithm: Nelder-Mead\n",
      " * Starting Point: [-1.2,1.0]\n",
      " * Minimizer: [1.0000435704852513,1.0000890072130664]\n",
      " * Minimum: 0.000000\n",
      " * Iterations: 79\n",
      " * Convergence: true\n",
      "   * |x - x'| < NaN: false\n",
      "   * |f(x) - f(x')| / |f(x)| < 1.0e-08: true\n",
      "   * |g(x)| < NaN: false\n",
      "   * Reached Maximum Number of Iterations: false\n",
      " * Objective Function Calls: 149\n",
      " * Gradient Calls: 0\n"
     ]
    }
   ],
   "source": [
    "# Define the main task\n",
    "maintask = current_task()\n",
    "\n",
    "# Define the stub function\n",
    "function f_t(x::Vector)\n",
    "    yieldto(maintask, x)\n",
    "end\n",
    "\n",
    "# Define an initial guess\n",
    "guess = [-1.2; 1.0]\n",
    "\n",
    "# Create a task for the solver\n",
    "solver = @task Optim.optimize(f_t, guess, NelderMead())\n",
    "\n",
    "# Driving loop\n",
    "next_x = copy(guess)\n",
    "while !istaskdone(solver)\n",
    "    fx = f(next_x)                 # Objective function evaluation\n",
    "    next_x = yieldto(solver, fx)   # Yielding back to the solver task\n",
    "end\n",
    "result = next_x\n",
    "\n",
    "@show result\n",
    ";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reverse communication example using objective and derivative functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "vals = Results of Nonlinear Solver Algorithm\n",
      " * Algorithm: Trust-region with dogleg and autoscaling\n",
      " * Starting Point: [-1.2,1.0]\n",
      " * Zero: [0.9999734745492139,0.9999379198581341]\n",
      " * Inf-norm of residuals: 0.000000\n",
      " * Iterations: 37\n",
      " * Convergence: true\n",
      "   * |x - x'| < 0.0e+00: false\n",
      "   * |f(x)| < 1.0e-08: true\n",
      " * Function Calls (f): 38\n",
      " * Jacobian Calls (df/dx): 25\n"
     ]
    }
   ],
   "source": [
    "# Main task\n",
    "maintask = current_task()\n",
    "\n",
    "# Stub functions\n",
    "# Include a symbol in \"vals\" for detecting which stub function was called within while loop\n",
    "function f_t!(x::Vector, fval::Vector)\n",
    "    vals = (:f, x, fval)\n",
    "    yieldto(maintask, vals)\n",
    "end\n",
    "function g_t!(x::Vector, fjac::Matrix)\n",
    "    vals = (:g, x, fjac)\n",
    "    yieldto(maintask, vals)\n",
    "end\n",
    "function fg_t!(x::Vector, fval::Vector, fjac::Matrix)\n",
    "    vals = (:fg, x, fval, fjac)\n",
    "    yieldto(maintask, vals)\n",
    "end\n",
    "\n",
    "# Initial guess and fvec allocations\n",
    "guess = [-1.2; 1.0]\n",
    "fval = [Inf; Inf]\n",
    "\n",
    "df = NLsolve.DifferentiableMultivariateFunction(f_t!, g_t!, fg_t!)\n",
    "\n",
    "# Create a task for the solver\n",
    "solver = @task NLsolve.nlsolve(df, guess)\n",
    "\n",
    "vals = (:f, guess, fval)\n",
    "while !istaskdone(solver)\n",
    "    if vals[1] == :f\n",
    "        guess = vals[2]\n",
    "        fval  = vals[3]\n",
    "        f!(guess, fval)              # Evaluate the objective function\n",
    "        vals = yieldto(solver, vals)\n",
    "    elseif vals[1] == :g\n",
    "        guess    = vals[2]\n",
    "        jacobian = vals[3]\n",
    "        g!(guess, jacobian)          # Evaluate the derivative function\n",
    "        vals = yieldto(solver, vals)\n",
    "    elseif vals[1] == :fg\n",
    "        guess    = vals[2]\n",
    "        fval     = vals[3]\n",
    "        jacobian = vals[4]\n",
    "        fg!(guess, fval, jacobian)   # Evaluate the combined objective and derivative function\n",
    "        vals = yieldto(solver, vals)\n",
    "    end\n",
    "end\n",
    "\n",
    "@show vals\n",
    ";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Embedding the reverse communication solver in a C calling environment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we replicate the above logic using Julia's C API, allowing for Julia to be embedded in a C application.\n",
    "\n",
    "In so doing, the objective and derivative functions can be defined as C functions instead of being Julia functions, meeting one of the client's main objectives."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```C\n",
    "//rc.c\n",
    "\n",
    "#include <julia.h>\n",
    "#include <stdio.h>\n",
    "#include <math.h>\n",
    "\n",
    "#ifdef _OS_WINDOWS_\n",
    "__declspec(dllexport) __cdecl\n",
    "#endif\n",
    "\n",
    "int r_f(jl_array_t *x, jl_array_t *f)\n",
    "{\n",
    "    double* xd = (double*)jl_array_data(x);\n",
    "    double* fd = (double*)jl_array_data(f);\n",
    "\n",
    "    fd[0] = pow(1.0 - xd[0], 2.0);\n",
    "    fd[1] = 100.0 * pow(xd[1] - pow(xd[0], 2.0), 2.0);\n",
    "\n",
    "    return 0;\n",
    "}\n",
    "\n",
    "int r_g(jl_array_t *x, jl_array_t *jac)\n",
    "{\n",
    "    double* xd = (double*)jl_array_data(x);\n",
    "    double* jd = (double*)jl_array_data(jac);\n",
    "    \n",
    "    jd[0] = -2.0 * (1 - xd[0]);                         //1,1 element\n",
    "    jd[1] = -400.0 * (xd[1] - pow(xd[0], 2.0)) * xd[0]; //2,1 element\n",
    "    jd[2] = 0.0;                                        //1,2 element\n",
    "    jd[3] = 200.0 * (xd[1] - pow(xd[0], 2.0));          //2,2 element\n",
    "\n",
    "    return 0;\n",
    "}\n",
    "\n",
    "void r_fg(jl_array_t *x, jl_array_t *f, jl_array_t *jac)\n",
    "{\n",
    "    double* xd = (double*)jl_array_data(x);\n",
    "    double* fd = (double*)jl_array_data(f);\n",
    "    double* jd = (double*)jl_array_data(jac);\n",
    "    \n",
    "    fd[0] = pow(1.0 - xd[0], 2.0);\n",
    "    fd[1] = 100.0 * pow(xd[1] - pow(xd[0], 2.0), 2.0);\n",
    "\n",
    "    jd[0] = -2.0 * (1 - xd[0]);                         //1,1 element\n",
    "    jd[1] = -400.0 * (xd[1] - pow(xd[0], 2.0)) * xd[0]; //2,1 element\n",
    "    jd[2] = 0.0;                                        //1,2 element\n",
    "    jd[3] = 200.0 * (xd[1] - pow(xd[0], 2.0));          //2,2 element\n",
    "}\n",
    "\n",
    "jl_value_t* yieldto(jl_value_t *solver, jl_value_t *vals) {\n",
    "    jl_function_t *yieldto_func = jl_get_function(jl_base_module, \"yieldto\");\n",
    "    return jl_call2(yieldto_func, solver, vals);\n",
    "}\n",
    "\n",
    "int istaskdone(jl_value_t *solver)\n",
    "{\n",
    "    jl_function_t *taskdone = jl_get_function(jl_base_module, \"istaskdone\");\n",
    "    return jl_call1(taskdone, solver) == jl_true;\n",
    "}\n",
    "\n",
    "#define __noinline __attribute__((noinline))\n",
    "void __noinline real_main()\n",
    "{\n",
    "    int n = 2;\n",
    "\n",
    "    // Represent arrays that will be created and owned by C\n",
    "    double *guess    = (double*)malloc(sizeof(double)*n);\n",
    "    double *fval     = (double*)malloc(sizeof(double)*n);\n",
    "    double *jacobian = (double*)calloc(n*n, sizeof(double));\n",
    "\n",
    "    // Assign initial values into the above arrays\n",
    "    guess[0] = -1.2;\n",
    "    guess[1] =  1.0;\n",
    "    fval[0]  =  1.0/0.0;\n",
    "    fval[1]  =  1.0/0.0;\n",
    "\n",
    "    // Root objects with the GC\n",
    "    jl_value_t **args;\n",
    "    // args[0] - NLsolve.DifferentiableMultivariateFunction(f!, g!, fg!)\n",
    "    // args[1] - :f\n",
    "    // args[2] - :g\n",
    "    // args[3] - :fg\n",
    "    // args[4] - array_type\n",
    "    // args[5] - guess::Vector{Float64}\n",
    "    // args[6] - fval::Vector{Float64}\n",
    "    // args[7] - jacobian::Matrix{Float64}\n",
    "    // args[8] - tuple function\n",
    "    // args[9] - solver_task function\n",
    "    \n",
    "    JL_GC_PUSHARGS(args,10);\n",
    "\n",
    "    jl_eval_string(\"using NLsolve\");\n",
    "    jl_eval_string(\"maintask = current_task()\");\n",
    "    jl_eval_string(\"f_t!(x::Vector, fvec::Vector) = yieldto(maintask, (:f, x, fvec))\");\n",
    "    jl_eval_string(\"g_t!(x::Vector, fjac::Matrix) = yieldto(maintask, (:g, x, fjac))\");\n",
    "    jl_eval_string(\"fg_t!(x::Vector, fvec::Vector, fjac::Matrix) = yieldto(maintask, (:fg, x, fvec, fjac))\");\n",
    "    \n",
    "    args[0] = jl_eval_string(\"NLsolve.DifferentiableMultivariateFunction(f_t!, g_t!, fg_t!)\");\n",
    "    args[1] = (jl_value_t*) jl_symbol(\"f\");\n",
    "    args[2] = (jl_value_t*) jl_symbol(\"g\");\n",
    "    args[3] = (jl_value_t*) jl_symbol(\"fg\");\n",
    "    args[4] = (jl_value_t*) jl_apply_array_type(jl_float64_type, 1 );\n",
    "    // Wrap the C arrays as Julia arrays. The final \"0\" argument means\n",
    "    // that Julia does not take ownership of the array data for GC purposes.\n",
    "    args[5] = (jl_value_t*) jl_ptr_to_array_1d(args[4], guess, n, 0);      \n",
    "    args[6] = (jl_value_t*) jl_ptr_to_array_1d(args[4], fval, n, 0);       \n",
    "    args[7] = (jl_value_t*) jl_ptr_to_array_1d(args[4], jacobian, n*n, 0);\n",
    "    args[8] = (jl_value_t*) jl_get_function(jl_base_module, \"tuple\");\n",
    "    args[9] = jl_eval_string(\"solver_task(df, guess) = @task nlsolve(df, guess)\");\n",
    "\n",
    "    { \n",
    "        jl_value_t* vals   = jl_call3((jl_function_t*) args[8], args[1], args[5], args[6]);  // \"vals\" tuple\n",
    "        jl_value_t* solver = jl_call2((jl_function_t*) args[9], args[0], args[5]);           // solver task\n",
    "        \n",
    "        JL_GC_PUSH2(&vals, &solver);\n",
    "\n",
    "        jl_yield();\n",
    "        while (!istaskdone(solver)) {\n",
    "            if (jl_get_nth_field(vals, 0) == args[1]) {        // :f case\n",
    "                args[5] = jl_get_nth_field(vals, 1);\n",
    "                args[6] = jl_get_nth_field(vals, 2);\n",
    "                r_f((jl_array_t*) args[5], (jl_array_t*) args[6]);\n",
    "                vals = yieldto(solver, vals);\n",
    "            } else if (jl_get_nth_field(vals, 0) == args[2]) { // :g case\n",
    "                args[5] = jl_get_nth_field(vals, 1);\n",
    "                args[7] = jl_get_nth_field(vals, 2);\n",
    "                r_g((jl_array_t*) args[5], (jl_array_t*) args[7]); \n",
    "                vals = yieldto(solver, vals);\n",
    "            } else if (jl_get_nth_field(vals, 0) == args[3]) { // :fg case\n",
    "                args[5] = jl_get_nth_field(vals, 1);\n",
    "                args[6] = jl_get_nth_field(vals, 2);\n",
    "                args[7] = jl_get_nth_field(vals, 3);\n",
    "                r_fg((jl_array_t*) args[5], (jl_array_t*) args[6], (jl_array_t*) args[7]); \n",
    "                vals = yieldto(solver, vals);\n",
    "            }\n",
    "        }\n",
    "        jl_show(jl_stderr_obj(), vals);\n",
    "        jl_eval_string(\"println(\\\"\\n\\\")\");\n",
    "        JL_GC_POP();\n",
    "    }\n",
    "    JL_GC_POP();\n",
    "    \n",
    "    free(guess);\n",
    "    free(fval);\n",
    "    free(jacobian);\n",
    "}\n",
    "\n",
    "void __noinline real_main_wrapper()\n",
    "{\n",
    "    real_main();\n",
    "}\n",
    "\n",
    "int main()\n",
    "{\n",
    "    jl_init(NULL);\n",
    "\n",
    "    real_main(); // this should be a separate function from main()\n",
    "\n",
    "    int ret = 0;\n",
    "    jl_atexit_hook(ret);\n",
    "    return ret;\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example gcc compilation\n",
    "\n",
    "```bash\n",
    "$ export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/julia-0.4.5/lib/julia\n",
    "$ export JULIA_HOME=/opt/julia-0.4.5/bin\n",
    "$ export JULIA_DIR=/opt/julia-0.4.5 \n",
    "$ gcc -o rc -fPIC -I$JULIA_DIR/include/julia -L$JULIA_DIR/lib/julia rc.c -ljulia $JULIA_DIR/lib/julia/libstdc++.so.6 -lm\n",
    "```\n",
    "\n",
    "Also see the [Embedding Julia documentation](http://docs.julialang.org/en/release-0.4/manual/embedding/) as well as the [embedding.c](https://github.com/JuliaLang/julia/blob/master/examples/embedding.c) and [Makefile](https://github.com/JuliaLang/julia/blob/master/examples/Makefile) examples in the Julia source code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.5",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}