@CarloLucibello, Since it says [WIP], let's review and merge this pr after you've done?

it is a WIP because it's not clear to me how the output should look when the input as a complicated structure, e.g. a vector of matrices

@denizyuret

I've got a multi-argument function and was trying to add the argnum argument to hessian (and therefore jacobian), but am failing. Can you, @CarloLucibello, provide some guidance on how to do that? I've got a function from Rn -> R and just want the hessian as an nxn matrix.

If it makes it easier in the short term, I only need argnum = 1.

If I understand correctly, you have a mutli-argument function, e.g. f(x,y), and want to compute the hessian just with respect to x. In this case you can try to use the methods in this PR as follow

x0 = ...
y0 = ....
hessian(x->f(x, y0))(x0)

that creates a closure over y0, right? My second argument is changing during my iterations, so I need to input it on each call.

you should be fine, the closure is created anew each time hessian(x->f(x, y0))(x0) is executed

@CarloLucibello what is the latest state of this PR? Is it ready to merge?