Re: Speeding up Replacement Rules
- To: mathgroup at smc.vnet.net
- Subject: [mg24572] Re: Speeding up Replacement Rules
- From: Jens-Peer Kuska <kuska at informatik.uni-leipzig.de>
- Date: Tue, 25 Jul 2000 00:56:05 -0400 (EDT)
- Organization: Universitaet Leipzig
- References: <8lgqqm$1vj@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi,
> Of course, if someone has Mathematica code which takes a vector
> valued function and generates a Compiled gradient or hessian
> function, I will accept it gratefully.
Here is the code that take a system of differential equations
and make a compiled function and a compiled jacobian. It is not
exactly the same but it may give you a hint.
Clear[dirkMakeFun]
(* We must wrap a Hold[] around the Part[new, _] to
prevent that the part of a symbol is serached *)
dirkMakeVars[y_List, new_Symbol] :=
Module[{x, i},
Thread[y -> Table[Hold[new[[x]]] /. x -> i, {i, Length[y]}]]
]
DIRKNDSolve::noexpl = "No explicit solution found for equations `1`."
dirkMakeFun[eqn : {_Equal ..}, y_List, t_Symbol, compile_:True] :=
Module[{sol, deqn, f, jac, rules, newsym, fj},
deqn = Select[eqn, ! FreeQ[#, t] &];
sol = Solve[deqn, D[#, t] & /@ y];
If[sol == {} || Head[sol] === Solve,
Message[DIRKNDSolve::noexpl, eqn];
Return[$Failed];
];
f = D[y, t] /. sol; (* explicit solution *)
jac = Outer[D[#1, #2] &, #, y] & /@ f;
newsym = Unique[];
rules = dirkMakeVars[y, newsym];
fj = Transpose[{f, jac}] /. rules;
If[compile,
(* at first we add the function arguments for compile
with an evaluated newsym, than Hold[] is removed
and the Hold[part[_]] are removed with the last Replace[] *)
fj = Map[Hold,
Map[{{{t, _Real}, {newsym, _Real, 1}}, #} & , fj, {2}] , {2}]
/.
Hold[a_Part] :> a;
(* Now we can start with compile. To prevent a early
evaluation of the Part[] function in the Apply[] call an
Unevaluated[]
is needed *)
fj /. Hold[l_] :> Apply[Compile, Unevaluated[l]]
(* Else *),
ReleaseHold[Map[Function[Evaluate[{t, newsym}], #] &, fj, {2}]]
]
dirkMakeFun[] generate for the variables in y a vector newsym so that
for
y={x[t],q[t],z[t]} the new vector ist newsym[[1]], newsym[[2]] and
newsym[[3]].
The tricky part is to prevent the evaluation of newsym[[i]] before
it is in the Compile[] command.
dirkMakeVars will return something like
{x[t]->Hold[newsym[[1]]],q[t]->Hold[newsym[[2]]],z[t]->Hold[newsym[[3]]]}
than I make the argument for Compile[] but with a hold arround the full
arguments and get
Hold[{{{newsym, _Real,1}},theExpressionWithManyHoldParts}}
since a Hold[] is around the expression I can remove the Hold[a_Part]
inside the expression. Finaly I have to Apply[] the Compile command
*and* to remove the Hold[]
fj /. Hold[l_] :> Apply[Compile, Unevaluated[l]]
and I have my pair of function/jacobian compiled functions for
the right hand side of my systems of differential equations.
Hope that helps
Jens