Skip to content

Commit

Permalink
Add macro with options
Browse files Browse the repository at this point in the history
  • Loading branch information
ChrisRackauckas committed Oct 19, 2016
1 parent 02aa0e2 commit 6484e98
Show file tree
Hide file tree
Showing 2 changed files with 105 additions and 65 deletions.
27 changes: 26 additions & 1 deletion README.md
Expand Up @@ -42,6 +42,26 @@ computed, a warning is thrown and only the function itself is usable. The boolea
`f.jac_exists` and `f.invjac_exists` can be used to see whether the Jacobian
and the function for its inverse exist.

#### Extra Options

In most cases the `@ode_def` macro should be sufficient. This is because by default
the macro will simply calculate each function symbolically, and if it can't it
will simply throw a warning and move on. However, in extreme cases the symbolic
calculations may take a long time, in which case it is necessary to turn them
off. To do this, use the `ode_def_opts` function. The `@ode_def` macro simply defines the specifiable options:

```julia
opts = Dict{Symbol,Bool}(
:build_Jac => true,
:build_InvJac => true,
:build_Hes => true,
:build_InvHes => true,
:build_dpfuncs => true)
```

and calls the function `ode_def_opts(name::Symbol,opts,ex::Expr,params)`. Note that
params is an iterator holding expressions for the parameters.

### Finite Element PDEs

Similar macros for finite element method definitions also exist. For the finite
Expand Down Expand Up @@ -76,7 +96,7 @@ The ParameterizedFunction Interface is as follows:
Hessians, Inverse Jacobians, Inverse Hessians, explicit parameter functions,
and parameter derivatives.
- The standard call `(p::TypeName)(t,u,du)` must be overloaded for the function
calculation.
calculation. All other functions are optional.

Solvers can interface with ParameterizedFunctions as follows:

Expand All @@ -96,6 +116,11 @@ It is requested that solvers should only use the explicit functions when they ex

## Manually Defining `ParameterizedFunction`s

It's recommended that for simple uses you use the macros. However, in many cases
the macros will not suffice, but you may still wish to provide Jacobians to the
solvers. This shows how to manually build a ParameterizedFunction to give to
a solver.

### Template

An example of explicitly defining a parameterized function is as follows. This serves
Expand Down
143 changes: 79 additions & 64 deletions src/ParameterizedFunctions.jl
Expand Up @@ -9,7 +9,7 @@ getindex{s}(p::ParameterizedFunction,::Val{s}) = getfield(p,s) ## Val for type-s

### Macros

macro ode_def(name,ex,params...)
function ode_def_opts(name::Symbol,opts::Dict{Symbol,Bool},ex::Expr,params...)
origex = ex # Save the original expression

## Build independent variable dictionary
Expand Down Expand Up @@ -42,50 +42,55 @@ macro ode_def(name,ex,params...)
hes_exists = false
invHex = :(error("Inverse Hessian Does Not Exist"))
invhes_exists = false
try #Jacobians and Hessian
# Build the Jacobian Matrix of SymEngine Expressions
numsyms = length(symtup)
symjac = Matrix{SymEngine.Basic}(numsyms,numsyms)
for i in eachindex(funcs)
funcex = funcs[i]
symfunc = @eval $funcex
for j in eachindex(symtup)
symjac[i,j] = diff(SymEngine.Basic(symfunc),symtup[j])
if opts[:build_Jac]
try #Jacobians and Hessian
# Build the Jacobian Matrix of SymEngine Expressions
numsyms = length(symtup)
symjac = Matrix{SymEngine.Basic}(numsyms,numsyms)
for i in eachindex(funcs)
funcex = funcs[i]
symfunc = @eval $funcex
for j in eachindex(symtup)
symjac[i,j] = diff(SymEngine.Basic(symfunc),symtup[j])
end
end
end

# Build the Julia function
Jex = build_jac_func(symjac,indvar_dict,param_dict,inline_dict)
jac_exists = true

try # Jacobian Inverse
invjac = inv(symjac)
invJex = build_jac_func(invjac,indvar_dict,param_dict,inline_dict)
invjac_exists = true
catch err
warn("Jacobian could not invert")
end

try # Hessian
symhes = Matrix{SymEngine.Basic}(numsyms,numsyms)
for i in eachindex(funcs), j in eachindex(symtup)
symhes[i,j] = diff(symjac[i,j],symtup[j])
end
# Build the Julia function
Hex = build_jac_func(symhes,indvar_dict,param_dict,inline_dict)
hes_exists = true

try # Hessian Inverse
invhes = inv(symhes)
invHex = build_jac_func(invhes,indvar_dict,param_dict,inline_dict)
invhes_exists = true
catch err
warn("Hessian could not invert")
Jex = build_jac_func(symjac,indvar_dict,param_dict,inline_dict)
jac_exists = true

if opts[:build_InvJac]
try # Jacobian Inverse
invjac = inv(symjac)
invJex = build_jac_func(invjac,indvar_dict,param_dict,inline_dict)
invjac_exists = true
catch err
warn("Jacobian could not invert")
end
end
if opts[:build_Hes]
try # Hessian
symhes = Matrix{SymEngine.Basic}(numsyms,numsyms)
for i in eachindex(funcs), j in eachindex(symtup)
symhes[i,j] = diff(symjac[i,j],symtup[j])
end
# Build the Julia function
Hex = build_jac_func(symhes,indvar_dict,param_dict,inline_dict)
hes_exists = true
if opts[:build_InvHes]
try # Hessian Inverse
invhes = inv(symhes)
invHex = build_jac_func(invhes,indvar_dict,param_dict,inline_dict)
invhes_exists = true
catch err
warn("Hessian could not invert")
end
end
end
end
catch err
warn("Failed to build the Jacoboian. This means the Hessian is not built as well.")
end

catch err
warn("Failed to build the Jacoboian. This means the Hessian is not built as well.")
end

# Parameter function calculations
Expand All @@ -101,32 +106,32 @@ macro ode_def(name,ex,params...)
end
pfuncs = build_p_funcs(paramfuncs,paramtup,indvar_dict,param_dict,inline_dict)

local d_pfuncs
local pderiv_exists
try # Parameter Gradients
d_paramfuncs = Vector{Vector{Expr}}(numparams)
for i in eachindex(paramtup)
tmp_dpfunc = Vector{Expr}(length(funcs))
for j in eachindex(funcs)
funcex = funcs[j]
symfunc = @eval $funcex
symfunc_str = parse(string(diff(SymEngine.Basic(symfunc),paramtup[i])))
if typeof(symfunc_str) <: Number
tmp_dpfunc[j] = :(1*$symfunc_str)
elseif typeof(symfunc_str) <: Symbol
tmp_dpfunc[j] = :(1*$symfunc_str)
else
tmp_dpfunc[j] = symfunc_str
d_pfuncs = Vector{Expr}(0)
pderiv_exists = false
if opts[:build_dpfuncs]
try # Parameter Gradients
d_paramfuncs = Vector{Vector{Expr}}(numparams)
for i in eachindex(paramtup)
tmp_dpfunc = Vector{Expr}(length(funcs))
for j in eachindex(funcs)
funcex = funcs[j]
symfunc = @eval $funcex
symfunc_str = parse(string(diff(SymEngine.Basic(symfunc),paramtup[i])))
if typeof(symfunc_str) <: Number
tmp_dpfunc[j] = :(1*$symfunc_str)
elseif typeof(symfunc_str) <: Symbol
tmp_dpfunc[j] = :(1*$symfunc_str)
else
tmp_dpfunc[j] = symfunc_str
end
end
d_paramfuncs[i] = tmp_dpfunc
end
d_paramfuncs[i] = tmp_dpfunc
d_pfuncs = build_p_funcs(d_paramfuncs,paramtup,indvar_dict,param_dict,inline_dict)
pderiv_exists = true
catch err
warn("Failed to build the parameter derivatives.")
end
d_pfuncs = build_p_funcs(d_paramfuncs,paramtup,indvar_dict,param_dict,inline_dict)
pderiv_exists = true
catch err
warn("Failed to build the parameter derivatives.")
d_pfuncs = Vector{Expr}(0)
pderiv_exists = false
end

# Build the type
Expand Down Expand Up @@ -461,7 +466,17 @@ end

const FEM_SYMBOL_DICT = Dict{Symbol,Expr}(:x=>:(x[:,1]),:y=>:(x[:,2]),:z=>:(x[:,3]))

export ParameterizedFunction, @ode_def, @fem_def
macro ode_def(name,ex,params...)
opts = Dict{Symbol,Bool}(
:build_Jac => true,
:build_InvJac => true,
:build_Hes => true,
:build_InvHes => true,
:build_dpfuncs => true)
ode_def_opts(name,opts,ex,params...)
end

export ParameterizedFunction, @ode_def, @fem_def, ode_def_opts
end # module


Expand Down

0 comments on commit 6484e98

Please sign in to comment.