RE: compile NIntegrate
- To: mathgroup at smc.vnet.net
- Subject: [mg17435] RE: [mg17402] compile NIntegrate
- From: "Ersek, Ted R" <ErsekTR at navair.navy.mil>
- Date: Sun, 9 May 1999 04:43:47 -0400
- Sender: owner-wri-mathgroup at wolfram.com
Francesco Siano wrote:
-------------------------
Hi, does anybody know if it is possible to compile efficiently a
function which, among other things, performs a numerical integration ?
Say, my function is
f[x_]=...NIntegrate[g[x-y],{y,-30,30}]...
Any help would be greatly appreciated. Thanks.
-------------------------
Note NIntegrate uses CompiledFunction as long as the option below is true
and the function you want to integrate can be converted to opcode. If the
function you want to integrate can't be converted to opcode by the
CompiledFunction interpreter standard Mathematica evaluation is used.
In[1]:=
Options[NIntegrate,Compiled]
Out[1]=
{Compiled\[Rule]True}
On the other hand maybe you wanted to do:
f[x_]:=Compile[{x}, ..... NIntegrate[g[x-y],{y,-30,30}], ..... ]
The CompiledFunction interpreter can't express NIntegrate using opcode, so
the use of Compile may not help. You might be better off using a
CompiledFunction to do everything (or almost everything) except NIntegrate.
You can then evaluate NIntegrate outside the CompiledFunction. (see below)
cf=Compile[{x}, ..... ]
f[x_]:= (cf[x];NIntegrate[g[x-y],{y,-30,30}])
Note:
Features the CompiledFunction interpreter can convert to opcode are listed
at
http://www.wolfram.com/support/Kernel/Symbols/System/Compile.html
Earlier I posted the discussion below and it may interest you.
------------------
As another member of the group indicated a CompiledFunction doesn't run very
fast if it can't be represented using all op-code instructions. If you have
a CompiledFunction named (func) you can see if any part of it didn't get
converted to op-code by evaluating func[[3]]. What you get back will mostly
be a list of lists. Most of the sub-lists will be vectors of integers.
The parts that don't make it into op-code look like
{31,Function[{x},......]}.
What you see inside Function will be a part of your program the interpreter
has a problem with.
The function 'CheckCompiled' below takes a CompiledFunction and returns a
list of all the expressions that don't get converted to op-code. If it
gives you an empty list you have an efficient CompiledFunction.
In[6]:=
CheckCompiled[func_CompiledFunction]:=
Cases[func[[3]],
bad_Function:>Extract[bad,2,HoldForm],{2}
]
______________________________
But even if every symbol you use is included on the web page above your
program still may not convert to op-code. I hear this can happen for
lots of reasons, but I only know of a few (each discussed below).
_____________________
A global variable can't be represented using op-code. The compiled function
(f2) below isn't converted to op-code because the global variable (w) is
used in the function. The best way around this that I know is to provide
the global variable to the function as an argument.
In[1]:=
w = 0.25;
In[2]:=
f2 = Compile[{x}, (x^2 + Sqrt[Exp[x]]+w)];
Evaluate f2[[3]] and you will see the op-code translation of this program.
______________________
You can't use opcode to change the value of a function's argument. The
CompiledFunction below doesn't convert to op-code because the function
changes the value of (x) (one of the function's arguments).
In[3]:=
f3=Compile[{x,y},
x=Floor[y];
x+y]
_____________________
A little more than a year ago a relatively long program using Compile
had the same problem as the short program below. The solution with
extensive discussion was provided by Dave Withoff (of Wolfram Research).
The program below isn't converted to op-code because of the line
If[x<0, a=1.0]. When x<0 the If returns a Real. When x>=0 the If
returns Null (Mathematica jargon for nothing). This is unacceptable
for the interpreter so it doesn't get converted to op-code.
In[4]:=
f4=Compile[{x},
Module[{a=2.0},
If[x<0,a=1.0];
x+a]
]
To fix the program above change (a=1.0) to (a=1.0;) as below. That way the
interpreter knows the If will always return 'Null'.
In[5]:=
(* This will convert to op-code *)
f4=Compile[{x},
Module[{a=2.0},
If[x<0,a=1.0;];
x+a]
]
______________________
Regards,
Ted Ersek