# Numerical integration of a numeric data available as a nested list

I have some numerical data in the form of a list with the following structure: {...{x,y,z},...} defining a surface z=z(x,y) in a 3D space (x,y,z). The data came from a simulation, and I am post-processing it within Mathematica. The precise numbers entering this list are not very important. One can think of this as of a small simplified example:

 lst = RandomReal[{-2, 2}, {200, 3}] /. {x_, y_, z_} -> {x, y,
1.5 Exp[-(x^2 + y^2) - x] + z/10};


However, my question is more broad than that. I am aware that there is a Mathematica package enabling one to apply derivatives to such numerical data in order to, say, numerically calculate gradients of z.

My question is, if there are some standard Mathematica approaches to numerically calculate integrals either of z=z(x,y), or of some function of z: f[z(x,y)] over some domain in the (x,y) plane.

Let me express my interest more clearly, I would be grateful for any proposal from you, but my primary interest is if a standard function or a standard package for this purpose exists.

-
I thank very much for the answers. I will work them through during next few days. It is difficult to decide which answer is closer to my aim. –  Alexei Boulbitch Feb 28 '13 at 13:23
There is, however, one quite clear outcome: there is no such built-in and ready to use Mathematica function making the job. Ok, good to know. –  Alexei Boulbitch Feb 28 '13 at 13:24
Probably it is already a technical question. In the data that I have in mind the area in the (x,y) plane is very large. Probably in such a case a good idea would be to split it into parts and make one of the procedures proposed in your answers for each part? In general, are there any restrictions on the size of the (x,y) domain? –  Alexei Boulbitch Feb 28 '13 at 13:37

You can install the Obtuse package for unstructured grid interpolation from here. Then interpolating your given data as a function func is simple.

Interpolation Part:

Off[General::compat];
Needs["Obtuse"];
On[General::compat];
lst=RandomReal[{-2,2},{200,3}];
cpts=lst[[All,1;;2]];
cptab=lst/.{x_,y_,z_}->{{x,y},1.5 Exp[-(x^2+y^2)-x]+z/10};
func=Interpolation[cptab,Method->"Delaunay"];
(*Fartest x,y coordinates in data *)
{xBD, yBD} = MapThread[Prepend, {{Min[#], Max[#]} & /@ Transpose[cpts], {x, y}}]


{{x, -1.97466, 1.99827}, {y, -1.97415, 1.98924}}

So the plot..

I have chosen the Delaunay interpolation method based on convex hull. You can see the convex-hull in transparent white.

Integration Part:

Lets try to integrate this function. Take help of the Boole function trick from @Mark McClure.

GraphicsMeshMeshInit[];
lst2D = Most /@ lst;
poly = lst2D[[ConvexHull[lst2D]]];
InsideTest[x_?NumericQ, y_?NumericQ] :=
Boole[InPolygonQ[poly, {x, y}]];
test[x_?NumericQ, y_?NumericQ] := InsideTest[x, y]*func[{x, y}];
NIntegrate[test[x, y], Evaluate@xBD, Evaluate@yBD,
Method -> {"AdaptiveQuasiMonteCarlo", "SymbolicProcessing" -> 0},
AccuracyGoal -> 3]


5.89323

Code for visualization:

Mainly copy-paste from Obtuse documentation!

Row@{Show[
ContourPlot[func[{x, y}], Evaluate@xBD, Evaluate@yBD,
PlotPoints -> 100, PerformanceGoal -> "Speed",
ColorFunction -> "Rainbow", PlotRange -> Full,
ImageSize -> {400, 400}],
ListPlot[cpts, PlotStyle -> {{Black, PointSize[Medium]}}],
ListPlot[
Append[cpts[[ConvexHull[cpts]]], cpts[[ConvexHull[cpts]]][[1]]],
Joined -> True,
PlotStyle -> Directive[White, Thickness[.011], Opacity[.5]]],
ImageSize -> 500],
Show[Plot3D[Evaluate@func[{x, y}], Evaluate@xBD, Evaluate@yBD,
ColorFunction -> "Rainbow", PlotRange -> Full],
ListPointPlot3D[Flatten /@ cptab,
PlotStyle -> {{PointSize[Large], Black}}, Boxed -> False,
Axes -> None], ImageSize -> 500]}

-
+1 I had no idea that package existed - very useful! It seems similar to this one –  gpap Feb 27 '13 at 16:23
Thank you. It is very interesting. I will try the method. There is however, the question concerning the domain size. Please have a look right below my question above. –  Alexei Boulbitch Feb 28 '13 at 13:42

I look forward to reading other answers to your question because I would like to know more myself.

When I faced the same problem, I created an interpolation function over an $(x,y)$-grid. This will effectively mean that your data will have been converted to a continuous (and, depending on the options you use in interpolation, even smooth) function $z(x,y)$ in the $(x,y)$ domain which you can integrate, differentiate, use as input to other functions etc.

The catch is that Interpolation in higher dimensions works only in regular grids. The problem of having non-regular grids has been addressed several times here: i.e. see interpolation of 3D data, Improved interpolation of mostly-structured 3d data or Interpolation of mostly-structured 3D data but it is mostly treated on a case by case basis. I would assume a higher order Interpolation function doesn't work on non-regular grids BECAUSE the non-regularity of a grid affects the approach one should be taking and I can't see how to automate this.

In any case, if your datasets are of the general form you have posted (i.e. come from a Poisson-sprinkling), there is a few functions in Haneberg's book (if you have access to a scientific library) that do this thing quite well. I can't post the code because I assume it's proprietary but it is not something I couldn't have written myself so, given your reputation, you'd probably find it easy to implement.

The idea is to create a regular grid by dividing up your range in $x$ and $y$, then create a nearest function that assigns $z$ values to the points in the grid nearest to your given points and then assign $z$ values to the rest of the points in the regular grid using some distance function of your choice. Once your regular grid has been populated with $z$ values, you can use Interpolation as usual.

InverseDistanceGrid2D from that book is the one I tend to prefer and it works quite well on the dataset you posted - I expect it to work well in any dataset with the same kind of sprinkling of points.

lst = RandomReal[{-2, 2}, {200, 3}] /. {x_, y_, z_} -> {x, y,
1.5 Exp[-(x^2 + y^2) - x] + z/10};
grid = Flatten[Table[{i, j}, {i, -2, 2, 0.2}, {j, -2, 2, 0.2}], 1];
(* grid chosen so it is roughly twice as dense as your original set*)
f = InverseDistanceGrid2D[lst, 3, 2, grid];
(*calling of magical function whose code I can't post*)
Show[Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}, PlotPoints -> 80, Lighting -> "Neutral"],
ListPointPlot3D[lst, PlotStyle -> Directive[{Red, PointSize[0.015]}]]]


I assume your data are more unpredictable than a Gaussian but if -like here- you know what function to expect, then FindFit is what you need.

-
Thank you very much. The book you mentioned is not available, but I have got the direction. Without the book it may be even better. –  Alexei Boulbitch Feb 28 '13 at 13:39
If you want to see the routines in Haneberg's book, you can download them here. It seems to be called ReciprocalDistanceGrid[] in the book, tho. Looking at the implementation, it seems to me that one can now implement this in a less cumbersome manner using Nearest[]. –  Guess who it is. Apr 24 '13 at 0:09
Thanks! The implementation I had was much quicker (and using Nearest) but I think that the obtuse package is a superset of that. It includes other algorithms and seems to be much quicker. Having looked around, I'd say it's the most complete collection of gridding algorithms for mathematica for irregularly spaced data. Inverse distance is definitely not as good as thin plate splines but much quicker for larger datasets. –  gpap Jun 4 '13 at 11:07

You could interpolate and then integrate the interpolating function. I don't know that this is perfect, but I think any approach you take is likely to error prone, as you presumably already have error in the data.

At any rate, here's how to use interpolation.

SeedRandom[1];
lst = RandomReal[{-2, 2}, {200, 3}] /. {x_, y_, z_} -> {x, y,
1.5 Exp[-(x^2 + y^2) - x] + z/10};
interpF = Interpolation[lst]


We've generated a function defined on most of $[-2,2]\times[-2,2]$. You should understand of course, that we shouldn't trust the function evaluated outside of this domain.

{interpF[0, 0], interpF[2, 2]}


To integrate this over an $xy$ domain, let's set up a funcion test that returns 0 or 1, depending on whether the point $(x,y)$ is inside the convex hull of the projection of your points onto the $xy$-plane.

GraphicsMeshMeshInit[];
lst2D = Most /@ lst;
poly = lst2D[[ConvexHull[lst2D]]];
test[x_?NumericQ, y_?NumericQ] := Boole[InPolygonQ[poly, {x, y}]];
Plot3D[test[x, y]*interpF[x, y], {x, -2, 2}, {y, -2, 2}]


Note that we generate a message, but test is zero for those points where extrapolation was used and the message generated. Similarly, we can NIntegrate:

NIntegrate[test[x, y]*interpF[x, y], {x, -2, 2}, {y, -2, 2},
Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}]

(* Out: 5.83902 *)
`
-
Thank you. I will try. –  Alexei Boulbitch Feb 28 '13 at 13:40

I know this is a Mathematica site, and my answer might be considered blasphemy, but it might solve your problem (which is the purpose of all this, I suppose). MATLAB has a function that does exactly that: it's called griddata.

Since this is post-processing, I thought you might as well pass your data through MATLAB before handing it to Mathematica. If you don't have a MATLAB license, there's an open source alternative, which I never tried but people say it's working.

-
Thank you, but unfortunately I am illiterate in MATLAB. I prefer, therefore, to stay within Mathematica, even if there are better solutions. –  Alexei Boulbitch Feb 28 '13 at 13:31