Finding ioctls with Clang and Cxx.jl

By: Julia Computing, Inc.

Re-posted from: http://juliacomputing.com/blog/2017/02/21/finding-ioctls-with-cxx.html

Among the more popular tools in the linux debugging toolbox is strace, which allows users to easily trace and print out all system calls made by a program (as well as the arguments to these system calls). I recently found myself writing a similar tool to trace all the requests made to a proprietary network driver. I had access to the sources of the userspace-facing API for this driver,
but strace proper did not know about it. I thus took this opportunity to write a general tool to extract ioctls from header files. The result is a compelling, but nevertheless compact enough for a blog post, application of Cxx.jl, the Julia-C++ FFI and Clang, LLVM’s C/C++ compiler project. In this blog post, I will walk you through my approach to this problem, highlighting both how to use Cxx.jl, and how to use the Clang C++ API. I will be focusing solely on extracing this data from header files. How to use it to write an strace like tool is a topic for another time.

Aside: About Cxx.jl

If you already know about Cxx.jl, feel free to move on to the next section. Cxx.jl is a julia package (available though the julia package manager), that allows julia to seamlessly interact with C++ code. It does this by taking advantage of Julia’s capabilities for staged code generation to put an actual C++ compiler into julia’s compilation pipeline. This looks roughtly as follows:

1. Julia parses code in julia syntax
2. The Cxx.jl package provides macros that translate from Julia to
   C++ code (either by punning on julia syntax, or by providing a string macro
   that allows the user to write C++ code directly). The package remembers
   what C++ code the user wants to run and leaves behind a generated function
   (basically a way for to ask the compiler to call back into the package
   when it wants to generate code for this particular function).
3. Later when Julia wants to run a function that includes, C++ code, it sees
   the generated function, calls back into Cxx.jl, which then performs
   proper type translation (of any julia values used in the C++ code and
   back to julia from C++ code), and creates a C++ AST which it then passes
   to Clang to compile. Clang compiles this AST and hands back LLVM IR, which
   Cxx.jl can then declare to julia is the implementation of the generated
   function.

Note that we could have generated native code in step 3, instead of LLVM IR,
but using LLVM IR, allows cross-language LTO-style optimization.

The easiest way to interact with Cxx.jl package, is through the C++ REPL that comes with the package. After the Cxx package is loaded, this mode is automatically added to the julia REPL and can be accessed by pressing the ‘<’ key:

julia> using Cxx

julia> # Press '<' here

C++ > #include <iostream>
true

C++ > std::cout << "Hello World" << std::endl;
Hello World

The problem

Before getting into the code, let’s first carefully understand the problem at hand.
We’re interested in ioctl requests made to a certain driver. ioctl is essentially the catch-all system call for all requests to drivers that don’t fit anywhere else.
Such requests generally look like so:

int result = ioctl(fd, REQUEST_CODE, argument);

Where fd is a file descriptor associated with a resource managed by the driver, and argument is generally either an integer or (more commonly) a pointer to a more complicated structure of arguments for this request. REQUEST_CODE is a driver-specific code that specified what kind of request to make. In practice, there are exceptions to these rules, for a variety of reasons, but that’s the general structure. So let’s look at how the possible ioctl requests (I’ll just call them ioctls for short, even though there’s only one ioctl system call) are declared in the kernel headers. To be concrete,
I’ll pick out the USB ioctls, but the discussion applies generally. Let’s look at an excerpt from the linked file:

#define USBDEVFS_SETCONFIGURATION  _IOR('U', 5, unsigned int)
#define USBDEVFS_GETDRIVER         _IOW('U', 8, struct usbdevfs_getdriver)
#define USBDEVFS_SUBMITURB32       _IOR('U', 10, struct usbdevfs_urb32)
#define USBDEVFS_DISCARDURB        _IO('U', 11)

Each of these lines defines an ioctl request (what I called request code above). Regular ioctls (ioctls defined like the ones above) have their request code split up as follows:

0xAAAABBCC
  \__/ | |
  Size | Code
    Category

which are the three values encoded by the #define above. The category is a (in theory unique by driver) 8-bit value that identifies to the kernel which driver to route the request to. The code is then used by the driver to identify the requested function. The size portion of the ioctl is ignored by the kernel, but may be used by the driver for backwards compatibility purposes. In the above define, the category is always 'U' (an ASCII-printable value is often chosen, but this is not a requirement), the numerical code follows, and lastly, we have the argument struct, which is used to compute the size.

For our ioctl dumper, we want four pieces of information:
1. The name of the ioctl
2. The category
3. The code
4. The argument struct (for size, as well as to extract field names such that we can print the argument structures in our ioctl dumper).

With a clear understand of what our goal is, let’s get to work!

Playing with the Preprocessor

It is probably possible to accomplish a lot of this using regexes or similar text processing, but there is a few distinct advantages to using a proper C compiler, such as clang for the task:
1. It has a correct preprocessor, so we can see though any defines, as well as making sure to ignore anything not reachable due to ifdef or similar
2. It is easier to use it to automatically extract the fieldnames/offset etc, while seeing through typedefs and anything else that might make it hard for a text processor to understand what’s going on.

So, let’s get us a Clang instance. Setting one up from the C++ API requires a bit of
boilerplate, but luckily for us, the Cxx.jl package, comes with the ability to create
separate Clang instances from the one it using to process C++ code:

julia> CCompiler = Cxx.new_clang_instance(
    false #= don't julia definitions =#,
    true #= C mode (as opposed to C++) =#)
Cxx.CxxInstance{2}()

Now, let’s use that instance, to load up the header file we discussed above:

julia> Cxx.cxxinclude(CCompiler, "linux/usbdevice_fs.h")
true

To achive our goal, we’ll need to manually work with Clang’s Parser and
Preprocessor objects, so let’s extract those for easy reference:

julia> PP = icxx"&$(Cxx.active_instances[2].CI)->getPreprocessor();"
(class clang::Preprocessor *) @0x000055ff269d2380

julia> P  = Cxx.active_instances[2].Parser
(class clang::Parser *) @0x000055ff26338870

Before, going on, let’s see what the preprocessor knows about one of the
macros we were looking at above:

C++ > $PP->dumpMacroInfo($PP->getIdentifierInfo("USBDEVFS_SETCONFIGURATION"))
MacroState 0x55ff26a9d3c8 USBDEVFS_SETCONFIGURATION
 DefMacroDirective 0x55ff1acbb980
  MacroInfo 0x55ff1acbb880
    #define <macro> _IOR('U', 5, unsigned int)

Ok, great. We have confirmed that the compiler parsed the header file and that it knows about our macro of interest. Let’s see where we can go from there. Consulting the clang documentation we find out about clang::Preprocessor::getMacroInfo and clang::MacroInfo::tokens, which would give us what we want. Let’s encode this into some julia functions for easy reference:

getIdentifierInfo(PP, name) = icxx"$PP->getIdentifierInfo($name);"
getMacroInfo(PP, II::pcpp"clang::IdentifierInfo") = icxx"$PP->getMacroInfo($II);"
getMacroInfo(PP, name::String) = getMacroInfo(PP, getIdentifierInfo(PP, name))
tokens(MI::pcpp"clang::MacroInfo") = icxx"$MI->tokens();"

We can now do:

julia> tokens(getMacroInfo(PP, "USBDEVFS_SETCONFIGURATION"))
(class llvm::ArrayRef<class clang::Token>) {
 .Data = (const class clang::Token *&) (const class clang::Token *) @0x000055ff269cec60
 .Length = (unsigned long &) 9
}

So we have our array of tokens. Of course, this is not very useful to us in this form, so let’s do two things. First, we’ll teach julia how to properly display Tokens:

# Convert Tokens that are identifiers to strings, we'll use these later
tok_is_identifier(Tok) = icxx"$Tok.is(clang::tok::identifier);"
Base.String(II::pcpp"clang::IdentifierInfo") = unsafe_string(icxx"$II->getName().str();")
function Base.String(Tok::cxxt"clang::Token")
    @assert tok_is_identifier(Tok)
    II = icxx"$Tok.getIdentifierInfo();"
    @assert II != C_NULL
    String(II)
end
getSpelling(PP, Tok) = unsafe_string(icxx"$PP->getSpelling($Tok);")
function Base.show(io::IO, Tok::Union{cxxt"clang::Token",cxxt"clang::Token&"})
    print(io, unsafe_string(icxx"clang::tok::getTokenName($Tok.getKind());"))
    print(io, " '", getSpelling(PP, Tok), "'")
end

Which’ll looks something like this (not I used the pointer from above)[^1]:

[^1] The astute reader may complain that I’m using the global PP instance to print this value. That is a valid complaint, and in the actual code, I made it an IOContext property, but I did not want to complicate this blog post with that discussion.

C++> *(clang::Token*) 0x000055ff269cec60
identifier '_IOR'

Great, we’re on the right track. Let’s also teach julia how to automatically iterate over llvm::ArrayRefs:

# Iteration for ArrayRef
import Base: start, next, length, done
const ArrayRef = cxxt"llvm::ArrayRef<$T>" where T
start(AR::ArrayRef) = 0
function next(AR::cxxt"llvm::ArrayRef<$T>", i) where T
    (icxx"""
        // Force a copy, otherwise we'll retain reference semantics in julia
        // which is not what people expect.
        $T element = ($AR)[$i];
        return element;
    """, i+1)
end
length(AR::ArrayRef) = icxx"$AR.size();"
done(AR::ArrayRef, i) = i >= length(AR)

Even though this may looks a bit complicated, all this is saying is that arrayrefs are indexed from one to AR.size(); and we can use the C++ bracket operator to access elements . With this defined, we get:

julia> collect(tokens(getMacroInfo(PP, "USBDEVFS_SETCONFIGURATION")))
9-element Array{Any,1}:
 identifier '_IOR'
 l_paren '('
 char_constant ''U''
 comma ','
 numeric_constant '5'
 comma ','
 unsigned 'unsigned'
 int 'int'
 r_paren ')'

We’re off to a great start.

Getting all the ioctls

As the previous section may have indicated, defining iteration over an object, is an enourmously powerful way to work with said object. Because everything in julia is generation, enabling iteration over an object, immidiately allows us to use any of the standard iteration tools (e.g. filters, maps, etc) to work with our objects.

With this, in mind, let’s see what we want to do. We know that ioctls are introduced by a macro that expands to _IO(...), _IOR(...), _IOW(...) or _IOWR(...). So let’s
define iteration over Clang’s identifier table and write down exactly that:

start(tab::rcpp"clang::IdentifierTable") = icxx"$tab.begin();"
next(tab::rcpp"clang::IdentifierTable", it) = (icxx"$it->second;", icxx"++$it;")
done(tab::rcpp"clang::IdentifierTable", it) = icxx"$it == $tab.end();"
length(tab::rcpp"clang::IdentifierTable") = icxx"$tab.size();"
# Get all identifier that are macros
macros = Iterators.filter(II->icxx"$II->hasMacroDefinition();", icxx"$PP->getIdentifierTable();")
# Expand into tuples of (II, tokens)
IItokens = map(II->(II, collect(tokens(getMacroInfo(PP, II)))), macros)
# Now filter down to the ones we're interested in
ioctl_defs = filter(IItokens) do x
      II, tokens = x
      isempty(tokens) && return false
      tok_is_identifier(tokens[1]) && String(tokens[1]) in ["_IO","_IOR","_IOW","_IOWR"]
  end;

And if all worked well, we end up with:

julia> map(x->(String(x[1]),x[2]), ioctl_defs)
34-element Array{Tuple{String,Array{Any,1}},1}:
 ("USBDEVFS_FREE_STREAMS", Any[identifier '_IOR', l_paren '(', char_constant ''U'', comma ',', numeric_constant '29', comma ',', struct 'struct', identifier 'usbdevfs_streams', r_paren ')'])
 ("USBDEVFS_BULK32", Any[identifier '_IOWR', l_paren '(', char_constant ''U'', comma ',', numeric_constant '2', comma ',', struct 'struct', identifier 'usbdevfs_bulktransfer32', r_paren ')'])
 ("USBDEVFS_DISCARDURB", Any[identifier '_IO', l_paren '(', char_constant ''U'', comma ',', numeric_constant '11', r_paren ')'])

Extracting the fields from the structures

Now, it’s fairly simple to do any any post-processing we want here, and what to do exactly will depend on our intended application, but I do want to highlight how to extract the fields. At first I attempted to simply use the second to last token as the type name, but that doesn’t work very well, because some types are multiple tokens (e.g. unsigned int) and some others are only defined via (sometimes complicated preprocessor rules). Instead, the right way to do this, is to simply feed those tokens back through the parser. We’ll use a couple of definitions

"Given an array of tokens, queue them up for parsing"
function EnterTokenStream(PP::pcpp"clang::Preprocessor", tokens::Vector{cxxt"clang::Token"})
  # Vector memory layout is incompatible, convert to clang::Token**
  toks = typeof(tokens[1].data)[x.data for x in tokens]
  icxx"$PP->EnterTokenStream(llvm::ArrayRef<clang::Token>{
    (clang::Token*)$(pointer(toks)),
    (size_t)$(length(toks))
  },false);"
end
"Advance the parse if it's currently at EOF. This happens in incremental parsing mode and should be called before parsing."
function AdvanceIfEof(P)
  icxx"""
  if ($P->getPreprocessor().isIncrementalProcessingEnabled() &&
    $P->getCurToken().is(clang::tok::eof))
      $P->ConsumeToken();
  """
end
"Parse a type name"
function ParseTypeName(P::pcpp"clang::Parser")
  AdvanceIfEof(P)
  res = icxx"$P->ParseTypeName(nullptr, clang::Declarator::TypeNameContext);"
  !icxx"$res.isUsable();" && error("Parsing failed")
  Cxx.QualType(icxx"clang::Sema::GetTypeFromParser($res.get());")
end
"Parse a constant expression"
function ParseConstantExpression(P::pcpp"clang::Parser")
  AdvanceIfEof(P)
  res = icxx"$P->ParseConstantExpression();"
  !icxx"$res.isUsable();" && error("Parsing failed")
  e = icxx"$res.get();"
  e
end
"Convert a parsed constant literal to a julia Char (asserts on failure)"
function CharFromConstExpr(e)
  Char(icxx"""
    clang::cast<clang::CharacterLiteral>($e)->getValue();
  """)
end
tok_is_comma(Tok) = icxx"$Tok.is(clang::tok::comma);"
tok_is_numeric(Tok) = icxx"$Tok.is(clang::tok::numeric_constant);"

With these definitions:

julia> ioctl_tokens = first(ioctl_defs)[2]
9-element Array{Any,1}:
 identifier '_IOR'
 l_paren '('
 char_constant ''U''
 comma ','
 numeric_constant '29'
 comma ','
 struct 'struct'
 identifier 'usbdevfs_streams'
 r_paren ')'

julia> typename_tokens = Vector{cxxt"clang::Token"}(ioctl_tokens[findlast(tok_is_comma, ioctl_tokens)+1:end-1])
2-element Array{Cxx.CppValue{Cxx.CxxQualType{Cxx.CppBaseType{Symbol("clang::Token")},(false, false, false)},N} where N,1}:
 struct 'struct'
 identifier 'usbdevfs_streams'

julia> EnterTokenStream(PP, typename_tokens); QT = Cxx.desugar(ParseTypeName(P))
Cxx.QualType(Ptr{Void} @0x000055ff24a13960)

C++ > ((clang::RecordType*)&*$QT)->getDecl()->dump()
RecordDecl 0x55951860c240 </usr/lib/gcc/x86_64-linux-gnu/6.2.0/../../../../include/linux/usbdevice_fs.h:153:1, line:157:1> line:153:8 struct usbdevfs_streams definition
|-FieldDecl 0x55951860c300 <line:154:2, col:15> col:15 num_streams 'unsigned int'
|-FieldDecl 0x55951860c358 <line:155:2, col:15> col:15 num_eps 'unsigned int'
`-FieldDecl 0x55951860c418 <line:156:2, col:21> col:16 eps 'unsigned char [0]'

We could process these for exmaple as such. Here I’ll be using a different approach, where instead of using julia to do the iteration, I’ll just write most of the function in C++ and only call back into julia once at the end:

# C structure to julia array of fields
function inspectStruct(CC, S)
    CC = Cxx.instance(CC)
    ASTCtx = icxx"&$(CC.CI)->getASTContext();"
    fields = Any[]
    icxx"""
    auto &ARL = $ASTCtx->getASTRecordLayout($S);
    for (auto field : ($S)->fields()) {
      unsigned i = field->getFieldIndex();
      // Skip these for now
      if (field->isImplicit())
        continue;
      if (field->getType()->isUnionType())
        continue;
      if (field->getType()->isArrayType())
        continue;
      if (field->getType()->isRecordType() ||
          field->getType()->isEnumeralType())
        continue;
      if (field->getType()->isPointerType() &&
          field->getType()->getPointeeOrArrayElementType()->isRecordType())
        continue;
      $:(begin
        QT = Cxx.QualType(icxx"return field->getType();")
        push!(fields, (
          String(icxx"return field;"),
          Cxx.juliatype(QT),
          icxx"return $ASTCtx->toCharUnitsFromBits(ARL.getFieldOffset(i)).getQuantity();"
        ))
      nothing
    end);
    }
    """
    fields
end
julia> inspectStruct(CCompiler, icxx"((clang::RecordType*)&*$QT)->getDecl();")
2-element Array{Any,1}:
 ("num_streams", UInt32, 0)
 ("num_eps", UInt32, 4)

Conclusion

With the above code, we can easily extract and work with the definitions of ioctls in the linux headers. I hope this blog post has given you an idea of both how to use the Clang C++ API to do some C introspection, as well as some idea, of how to use some of the generic programming features in julia. The above is a pretty decent summary some of the first things I do when working with new data sources in julia:
1. Defining printing method for the relevant types
2. Define iteration on any container data structures
3. Use julia’s iteration tools to write whatever query I’m interested in
Following this strategy usually gets one pretty far. In this case, it was essentially sufficient to solve our problem and provide a useful list of ioctls and the fields of their arguments to use in our ioctl dumping tool.

JuliaPro Featured in Danske Bank’s Business Analytics Challenge 2017

By: Julia Computing, Inc.

Re-posted from: http://juliacomputing.com/press/2017/02/21/danske.html

Copenhagen, Denmark – Danske Bank, Denmark’s largest bank, announced that JuliaPro will be available on Microsoft Azure’s Data Science Virtual Machine (DSVM) for participants in the Business Analytics Challenge 2017.

The Business Analytics Challenge 2017 is sponsored by Danske Bank, Microsoft and KMD. The competition is open to all undergraduate and master’s degree students in Denmark and the first prize is 75 thousand kroner. Registration is open until March 31.

This announcement comes two months after the release of JuliaPro and one month after JuliaPro launched on Microsoft Azure’s Data Science Virtual Machine (DSVM).

Viral Shah, Julia Computing CEO says, “We are thrilled that Julia adoption is accelerating so rapidly during the first quarter of 2017. In the last three months, we introduced the new JuliaPro and launched it on the world’s two largest cloud environments: Amazon’s AWS and Microsoft Azure’s Data Science Virtual Machine (DSVM). Julia Computing wishes the best of luck to all contestants in the Danske Bank Business Analytics Challenge 2017.”

About Julia Computing and Julia

Julia Computing (JuliaComputing.com) was
founded in 2015 by the co-creators of the Julia language to provide
support to businesses and researchers who use Julia.

Julia is the fastest modern high performance open source computing language for data and analytics. It combines the functionality and ease of use of Python, R, Matlab, SAS and Stata with the speed of Java and C++. Julia delivers dramatic improvements in simplicity, speed, capacity and productivity.

  1. Julia is lightning fast. Julia provides speed improvements up to
    1,000x for insurance model estimation, 225x for parallel
    supercomputing image analysis and 11x for macroeconomic modeling.

  2. Julia is easy to learn. Julia’s flexible syntax is familiar and
    comfortable for users of Python, R and Matlab.

  3. Julia integrates well with existing code and platforms. Users of
    Python, R, Matlab and other languages can easily integrate their
    existing code into Julia.

  4. Elegant code. Julia was built from the ground up for
    mathematical, scientific and statistical computing, and has advanced
    libraries that make coding simple and fast, and dramatically reduce
    the number of lines of code required – in some cases, by 90%
    or more.

  5. Julia solves the two language problem. Because Julia combines
    the ease of use and familiar syntax of Python, R and Matlab with the
    speed of C, C++ or Java, programmers no longer need to estimate
    models in one language and reproduce them in a faster
    production language. This saves time and reduces error and cost.

Employers looking to hire Julia programmers in 2017 include: Google, Apple, Amazon, Facebook, IBM, BlackRock, Capital One, PricewaterhouseCoopers, Ford, Oracle, Comcast, Massachusetts General Hospital, NaviHealth, Harvard University, Columbia University, Farmers Insurance, Pilot Flying J, Los Alamos National Laboratory, Oak Ridge National Laboratory and the National Renewable Energy Laboratory.

Julia users and partners include: Amazon, IBM, Intel, Microsoft, DARPA, Lawrence Berkeley National Laboratory, National Energy Research Scientific Computing Center (NERSC), Federal Aviation Administration (FAA), MIT Lincoln Labs, Moore Foundation, Nobel Laureate Thomas J. Sargent, Federal Reserve Bank of New York (FRBNY), Capital One, Brazilian National Development Bank (BNDES), BlackRock, Conning, Berkery Noyes, BestX, Path BioAnalytics, Invenia, AOT Energy, AlgoCircle, Trinity Health, Gambit, Augmedics, Tangent Works, Voxel8, UC Berkeley Autonomous Race Car (BARC) and many of the world’s largest investment banks, asset managers, fund managers, foreign exchange analysts, insurers, hedge funds and regulators.

Universities and institutes using Julia include: MIT, Caltech, Stanford, UC Berkeley, Harvard, Columbia, NYU, Oxford, NUS, UCL, Nantes, Alan Turing Institute, University of Chicago, Cornell, Max Planck Institute, Australian National University, University of Warwick, University of Colorado, Queen Mary University of London, London Institute of Cancer Research, UC Irvine, University of Kaiserslautern.

Julia is being used to: analyze images of the universe and research dark matter, drive parallel computing on supercomputers, diagnose medical conditions, provide surgeons with real-time imagery using augmented reality, analyze cancer genomes, manage 3D printers, pilot self-driving racecars, build drones, improve air safety, manage the electric grid, provide analytics for foreign exchange trading, energy trading, insurance, regulatory compliance, macroeconomic modeling, sports analytics, manufacturing and much, much more.

DynProg Class – Week 2

By: pkofod

Re-posted from: http://www.pkofod.com/2017/02/18/dynprog-class-week-2/

This post, and other posts with similar tags and headers, are mainly directed at the students who are following the Dynamic Programming course at Dept. of Economics, University of Copenhagen in the Spring 2017. The course is taught using Matlab, but I will provide a few pointers as to how you can use Julia to solve the same problems. So if you are an outsider reading this I welcome you, but you won’t find all the explanations and theory here. If you want that, you’ll have to come visit us at UCPH and enroll in the class!

This week we continue with a continuous choice model. This means we have to use interpolation and numerical optimization.

A (slightly less) simple model

Consider a simple continuous choice consumption-savings model:

\(V_t(M_t) = \max_{C_t}\sqrt{C_t}+\mathbb{E}[V_{t+1}(M_{t+1})|C_t, M_t]\)
subject to
\(M_{t+1} = M_t – C_t+R_t\\
C_t\leq M_t\\
C_t\in\mathbb{R}_+\)

where \(R_t\) is 1 with probability \(\pi\) and 0 with probability \(1-\pi\), \(\beta=0.9\), and \(\bar{M}=5\)

Last week the maximization step was merely comparing values associated with the different discrete choices. This time we have to do continuous optimization in each time period. Since the problem is now continuous, we cannot solve for all \(M_t\). Instead, we need to solve for \(V_t\) at specific values of \(M_t\), and interpolate in between. Another difference to last time is the fact that the transitions are stochastic, so we need to form (very simple) expectations as part of the Bellman equation evaluations.

Interpolation

It is of course always possible to make your own simple interpolation scheme, but we will use the functionality provided by the Interpolations.jl package. To perform interpolation, we need a grid to be used for interpolation \(\bar{x}\), and calculate the associated function values.

f(x) = (x-3)^2
x̄ = linspace(1,5,5)
fx̄ = f.()

Like last time, we remember that the dot after the function name and before the parentheses represent a “vectorized” call, or a broadcast – that is we call f on each element of the input. We now use the Interpolations.jl package to create an interpolant \(\hat{f}\).

using Interpolations
f̂ = interpolate((collect(),), fx̄, Gridded(Linear()))

We can now index into \(\hat{f}\) as if it was an array

[-3] #returns 16.0

We can also plot it

using Plots
plot()

which will output

Solving the model

Like last time, we prepare some functions, variables, and empty containers (Arrays)

# Model again
u(c) = sqrt(c)
T = 10; β = 0.9
π = 0.5 ;M₁ = 5
# Number of interpolation nodes
Nᵐ = 50 # number of grid points in M grid
Nᶜ  = 50 # number of grid points in C grid
 
M = Array{Vector{Float64}}(T)
V = Array{Any}(T)
C = Array{Any}(T)

The V and C arrays are allocated using the type “Any”. We will later look at how this can hurt performance, but for now we will simply do the convenient thing. Then we solve the last period

M[T] = linspace(0,M₁+T,Nᵐ)
C[T] = M[T]
V[T] = interpolate((M[T],), u.(C[T]), Gridded(Linear()))

This new thing here is that we are not just saving V[T] as an Array. The last element is the interpolant, such that we can simply index into V[T] as if we had the exact solution at all values of M (although we have to remember that it is an approximation). For all periods prior to T, we have to find the maximum as in the Bellman equation from above. To solve this reduced “two-period” problem (sum of utility today and discounted value tomorrow), we need to form expectations over the two possible state transitions given an M and a C today, and then we need to find the value of C that maximizes current value. We define the following function to handle this

# Create function that returns the value given a choice, state, and period
v(c, m, t, V) = u(c)*(π*V[t+1][m-c+1]+(1)*V[t+1][m-c])

Notice how convenient it is to simply index into V[t] using the values we want to evaluate tomorrow’s value function at. We perform the maximization using grid search on a predefined grid from 0 to the particular M we’re solving form. If we abstract away from the interpolation step, this is exactly what we did last time.

for t = T-1:-1:1
    M[t] = linspace(0,M₁+t,Nᵐ)
    C[t] = zeros(M[t])
    Vt = fill(-Inf, length(M[t]))
    for (iₘ, m) = enumerate(M[t])
        for c in linspace(0, m, Nᶜ)
            _v = v(c, m, t, V)
            if _v >= Vt[iₘ]
                Vt[iₘ] = _v
                C[t][iₘ] = c
            end
        end
    end
    V[t] = interpolate((M[t],), Vt, Gridded(Linear()))
end

Then we can plot the value functions to verify that they look sensible

Nicely increasing in time and in M.

Using Optim.jl for optimization

The last loop could just as well have been done using a proper optimization routine. This will in general be much more robust, as we don’t confine ourselves to a certain amount of C-values. We use the one of the procedures in Optim.jl. In Optim.jl, constrained, univariate optimization is available as either Brent’s method or Golden section search. We will use Brent’s method. This is the standard method, so an optimization call simply has the following syntax

using Optim
f(x) = x^2
optimize(f, -1.0, 2.0)

Unsurprisingly, this will return the global minimizer 0.0. However, if we constrain ourselves to a strictly positive interval

optimize(f, 1.0, 2.0)

we get a minimizer of 1.0. This is not the unconstrained minimizer of the square function, but it is minimizer given the constraints. Then, it should be straight forward to see how the grid search loop can be converted to a loop using optimization instead.

for t = T-1:-1:1
    update_M!(M, M₁, t, Nᵐ)
    C[t] = zeros(M[t])
    Vt = fill(-Inf, length(M[t]))
    for (iₘ, m) = enumerate(M[t])
        if m == 0.0
            C[t][iₘ] = m
            Vt[iₘ] = v(m, m, t, V)
        else
            res = optimize(c->-v(c, m, t, V), 0.0, m)
            Vt[iₘ] = -Optim.minimum(res)
            C[t][iₘ] = Optim.minimizer(res)
        end
    end
    V[t] = interpolate((M[t],), Vt, Gridded(Linear()))
end

If our agent has no resources at the beginning of the period, the choice set has only one element, so we skip the optimization step. We also have to negate the minimum to get the maximum we’re looking for. The main advantage of using a proper optimization routine is that we’re not restricting C to be in any predefined grid. This increases precision. If we look at the number of calls to “v” (using Optim.f_calls(res)), we see that it generally takes around 10-30 v calls. With only 10-30 grid points from 0 up to M, we would generally get an approximate solution of much worse quality.

Julia bits and pieces

This time we used a few different package from the Julia ecosystem: Plots.jl, Interpolations.jl, and Optim.jl. These are based on my personal choices (full disclosure: I’ve contributed to the former and the latter), but there are lots of packages to explore. Visit the JuliaLang discourse forum or gitter channel to discuss Julia and the various packages with other users.