Update: I've found a significantly cleaner way to do this which avoids
FastAnonymous entirely and runs faster for me. Now in 0.4 arbitrary objects
can be overloaded to use the f() syntax using Base.call. I create an
accessory type for each function that I want to pass to a non-linear solver
or optimization routine. This accessory type encapsulates the parameters of
the function, Then I define a new method for Base.call which acts like a
single-variable function, and simply redirects the parameters stored in my
type into the function. It looks like this
function
hoursFOC(UF::UtilityFunction,hours,state::State{IdioState1,AggState1},
aprime)
w = state.as.w
consump = budget_consump_givenhours(state, aprime, hours)
w*uc(UF,consump,hours) - ul(UF,consump,hours) # Is w*uc = ul
end
immutable pars_hoursFOC{TUF <: UtilityFunction, TIS<:IdioState,
TAS<:AggState}
UF::TUF
state::State{TIS,TAS}
aprime::Float64
end
Base.call(f::pars_hoursFOC, h) = hoursFOC(f.UF, h, f.state, f.aprime)
When I need the nonlinear equation solver, I do
f = pars_hoursFOC(UF,state,aprime)
j = pars_JACOBhoursFOC(UF, state, aprime) # jacobian, same idea
hours = myNewton(f, j, hoursguess)
f and j are treated like functions, so this works as expected.
This is much cleaner than my previous method where I was using @anon to
store an anonymous function, then passing that function around everywhere.
Also it's faster for some reason. My new code runs in 14s, the version
using @anon was about 20s, and a version with nested functions(or regular
anonymous functions) would be >30s.
I see mentioned in #8712 <https://github.com/JuliaLang/julia/pull/8712> that
there may be a Callable type. That would be useful, since currently these
callable types don't work with the Optim functions or anything else that
asks for a ::Function argument.
On Thursday, July 16, 2015 at 9:20:51 PM UTC-4, Andrew wrote:
>
> fzero(f, j, guess) works for me when f and j are functions. f(Af, guess)
> works for me now when Af is an @anon function.
>
> On Tuesday, July 7, 2015 at 7:34:39 PM UTC-4, j verzani wrote:
>>
>> Okay, this just got fixed as much as I could with v"0.1.15" (there is no
>> fzero(f,j,guess) signature).
>>
>> On Tuesday, July 7, 2015 at 4:38:41 PM UTC-4, Andrew wrote:
>>>
>>> Just checked. So, Roots.fzero(f, guess) does work. However,
>>> Roots.fzero(f, j, guess) doesn't work, and neither does Roots.newton(f, j,
>>> guess).
>>>
>>> I looked at the Roots.jl source and I see ::Function annotations on the
>>> methods with the jacobian, but not the regular one.
>>>
>>> On Tuesday, July 7, 2015 at 4:22:17 PM UTC-4, j verzani wrote:
>>>>
>>>> It isn't your first choice, but `Roots.fzero` can have `@anon`
>>>> functions passed to it, unless I forgot to tag a new version after making
>>>> that change on master not so long ago.
>>>>
>>>> On Tuesday, July 7, 2015 at 2:29:51 PM UTC-4, Andrew wrote:
>>>>>
>>>>> I'm writing this in case other people are trying to do the same thing
>>>>> I've done, and also to see if anyone has any suggestions.
>>>>>
>>>>> Recently I have been writing some code that requires solving lots(tens
>>>>> of thousands) of simple non-linear equations. The application is
>>>>> economics,
>>>>> I am solving an intratemporal first order condition for optimal labor
>>>>> supply given the state and a savings decision. This requires solving the
>>>>> same equation many times, but with different parameters.
>>>>>
>>>>> As far as I know, the standard ways to do this are to either define a
>>>>> nested function which by the lexical scoping rules inherits the
>>>>> parameters
>>>>> of the outer function, or use an anonymous function. Both these methods
>>>>> are
>>>>> slow right now because Julia can't inline those functions. However, the
>>>>> FastAnonymous package lets you define an anonymous "function", which
>>>>> behaves exactly like a function but isn't type ::Function, which is fast.
>>>>> Crucially for me, in Julia 0.4 you can modify the parameters of the
>>>>> function you get out of FastAnonymous. I rewrote some code I had which
>>>>> depended on solving a lot of non-linear equations, and it's now 3 times
>>>>> as
>>>>> fast, running in 2s instead of 6s.
>>>>>
>>>>> Here I'll describe a simplified version of my setup and point out a
>>>>> few issues.
>>>>>
>>>>> 1. I store the anonymous function in a type that I will pass along to
>>>>> the function which needs to solve the nonlinear equation. I use a
>>>>> parametric type here since the type of an anonymous function seems to
>>>>> vary
>>>>> with every instance. For example,
>>>>>
>>>>> typeof(UF.fhoursFOC)
>>>>> FastAnonymous.##Closure#11431{Ptr{Void}
>>>>> @0x00007f2c2eb26e30,0x10e636ff02d85766,(:h,)}
>>>>>
>>>>>
>>>>> To construct the type,
>>>>>
>>>>> immutable CRRA_labor{T1, T2} <: LaborChoice # <: means "subtype of"
>>>>> sigmac::Float64
>>>>> sigmal::Float64
>>>>> psi::Float64
>>>>> hoursmax::Float64
>>>>> state::State # Encodes info on how to solve itself
>>>>> fhoursFOC::T1
>>>>> fJACOBhoursFOC::T2
>>>>> end
>>>>>
>>>>> To set up the anonymous functions fhoursFOC and fJACOBhoursFOC (the
>>>>> jacobian), I define a constructor
>>>>>
>>>>> function CRRA_labor(sigmac,sigmal,psi,hoursmax,state)
>>>>> fhoursFOC = @anon h -> hoursFOC(CRRA_labor(sigmac,sigmal,psi,
>>>>> hoursmax,state,0., 0.) , h, state)
>>>>> fJACOBhoursFOC = @anon jh -> JACOBhoursFOC(CRRA_labor(sigmac,
>>>>> sigmal,psi,hoursmax,state,0., 0.) , jh, state)
>>>>> CRRA_labor(sigmac,sigmal,psi,hoursmax,state,fhoursFOC,
>>>>> fJACOBhoursFOC)
>>>>> end
>>>>>
>>>>> This looks a bit complicated because the nonlinear equation I need to
>>>>> solve, hoursFOC, relies on the type CRRA_labor, as well as some aggregate
>>>>> and idiosyncratic state info, to set up the problem. To encode this
>>>>> information, I define a dummy instance of CRRA_labor, where I supply 0's
>>>>> in
>>>>> place of the anonymous functions. I tried to make a self-referential type
>>>>> here as described in the documentation, but I couldn't get it to work, so
>>>>> I
>>>>> went with the dummy instance instead.
>>>>>
>>>>> @anon sets up the anonymous function. This means that code like
>>>>> fhoursFOC(0.5) will return a value.
>>>>>
>>>>> 2. Now that I have my anonymous function taking only 1 variable, I can
>>>>> use the nonlinear equation solver. Unfortunately, the existing nonlinear
>>>>> equation solvers like Roots.fzero and NLsolve ask the argument to be of
>>>>> type ::Function. Since anonymous functions work like functions but are
>>>>> actually some different type, they wouldn't accept my argument. Instead,
>>>>> I
>>>>> wrote my own Newton method, which is like 5 lines of code, where I don't
>>>>> restrict the argument type.
>>>>>
>>>>> I think it would be very straightforward to make this a multivariate
>>>>> Newton method.
>>>>>
>>>>> function myNewton(f, j, x)
>>>>> for n = 1:100
>>>>> fx , jx = f(x), j(x)
>>>>> abs(fx) < 1e-6 && return x
>>>>> d = fx/jx
>>>>> x = x - d
>>>>> end
>>>>> println("Too many iterations")
>>>>> return NaN
>>>>> end
>>>>>
>>>>> 3. The useful thing here in 0.4 is that you can edit the parameters of
>>>>> the anonymous function. The parameters are encoded in a custom type
>>>>> state::State, and I update the state. Then I call my nonlinear equation
>>>>> solver
>>>>>
>>>>> UF.fhoursFOC.state, UF.fJACOBhoursFOC.state = state, state
>>>>> f = UF.fhoursFOC
>>>>> j = UF.fJACOBhoursFOC
>>>>> hours = myNewton(f, j, hoursguess)
>>>>>
>>>>> This runs much faster than my old version which used NLsolve, which
>>>>> itself ran faster than a version using Roots.fzero.
>>>>>
>>>>> Issues:
>>>>>
>>>>> 1. Since the type of the anonymous function isn't ::Function, I had to
>>>>> write my own solver. I'm pretty sure a 1-line edit to Roots.fzero where I
>>>>> just remove the ::Function type annotation would let it work there, but
>>>>> I'm
>>>>> not aware of another workaround.
>>>>>
>>>>> 2. I would rather use NLsolve, which uses in-place updating of its
>>>>> arguments ( f!(input::Array, output::Array) ), but I've tried
>>>>> constructing
>>>>> an anonymous function that does that, and @anon didn't work. Perhaps
>>>>> there
>>>>> is a workaround.
>>>>>
>>>>> 3. Since I'm using an anonymous function, I have to explicitly pass it
>>>>> around. Encoding it into the type CRRA_labor wasn't really hard though.
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>