Interpolating four-dimensional data

Question

In this case it is not easy to work with a MWE, so I shall use the actual data which contain initial conditions $(x,y,p_x,p_y)$.

Let's see how the data look on the $(x,y)$ plane

data = Import["4D_3200.dat", "Table"];
d0 = data[[All, {1, 2}]];
L0 = ListPlot[d0, PlotStyle -> {Blue, PointSize[0.001]}]

It is seen that the data are not random but they fill an oval region. The data file contains about 21000 initial conditions.

Question #1: How can I generate, let's say 100000, random $(x,y)$ initial conditions inside the oval region of the above plot?

Question #2: Assuming that we have generated 100000 $(x,y)$ initial conditions, inside the oval region, how can we exploit the original data in order to interpolate the momenta $p_x$ and $p_y$ of the additional initial conditions?

In other words, we have an initial distribution of $N$ four-dimensional initial conditions and we want to generate more $N'$.

Any suggestions?

Many thanks in advance!

Are px and py dependent variables? What is the relationship between them? — MikeY, Apr 24 '17 at 18:14
@MikeY There is no known relationship between the momenta. The values of the initial data file should be used for the interpolation. — Vaggelis_Z, Apr 24 '17 at 18:18

Jason B. · Accepted Answer · 2017-04-24T20:30:19.363

Before doing anything, you have to remove any duplicate $(x,y)$ pairs from your data. Interpolating functions don't like double-valued functions. 10 of your points are like these two:

data[[{3197, 3250}]]
(* {{10.20913`, -0.18115`, -0.15105`, 50.42686`}, 
    {10.20913`, -0.18115`, -0.15105`, 50.42685`}} *)

So just throw out the duplicates via

data = DeleteDuplicatesBy[data, #[[;; 2]] &];
xy = data[[All, ;; 2]];

How can I generate, let's say 100000, random (x,y) initial conditions inside the oval region of the above plot?

You can define the oval region using ConvexHullMesh, and then use RandomPoint to generate your points.

xy = data[[All, ;; 2]];
pts = RandomPoint[ConvexHullMesh[xy], 100000];
Graphics[{Point[pts], Red, Point[xy]}]

Mathematica graphics

Assuming that we have generated 100000 (x,y) initial conditions, inside the oval region, how can we exploit the original data in order to interpolate the momenta $p_x$ and $p_y$ of the additional initial conditions?

Now just make your interpolation functions. You need to set the InterpolationOrder to 1 because the $(x,y)$ points don't lie on a rectangular grid. Check the documentation for Interpolation to see why you have to rearrange the list structure.

px = Interpolation[{{#1, #2}, #3} & @@@ data, 
  InterpolationOrder -> 1]
py = Interpolation[{{#1, #2}, #4} & @@@ data, 
   InterpolationOrder -> 1];

and generate the new data via

newData = {#1, #2, px[##], py[##]} & @@@ pts;

Here is the new data (in blue) along with the old data (in red) for the px variable

Graphics3D[{
  Red, PointSize@Medium, Point[data[[All, ;; 3]]],
  Blue, PointSize@Small, Point[
   newData[[All, ;; 3]]]}]

Mathematica graphics

Edit If you are working in an older version of Mathematica, before there was a ConvexHullMesh or RandomPoint, you can still make this work, but the only way I could think was to use rejection sampling to get the new data points.

positionDuplicates[list_] := 
    Flatten[Rest /@ GatherBy[Range@Length[list], list[[#]] &]]; 
data = ReplacePart[data, 
    Thread[positionDuplicates[data[[All, ;; 2]]] -> Sequence[]]]; 
xy = data[[All, ;; 2]]; 
px = Interpolation[{{#1, #2}, #3} & @@@ data, InterpolationOrder -> 1]; 
py = Interpolation[{{#1, #2}, #4} & @@@ data, InterpolationOrder -> 1]; 
range = {Min@#, Max@#} & /@ Transpose[xy]; 
pointsBeta = Transpose[{ 
    RandomReal[First@range, 100000], 
    RandomReal[Last@range, 100000] 
    } 
]; 
pgon = Polygon[ 
    data[[Graphics`Mesh`ConvexHull[xy], ;; 2]] 
]; 
pts = Select[pointsBeta, Graphics`Mesh`InPolygonQ[pgon, #] &];
newData = {#1, #2, px[##], py[##]} & @@@ pts; 

Graphics3D[{ 
    Red, PointSize@Medium, Point[data[[All, ;; 3]]], 
    Blue, PointSize@Small, Point[ newData[[All, ;; 3]]]
    }, Method -> {"ShrinkWrap" -> True}]

should give the same result, but only giving 80,000 points instead of 100,000. Code borrowed and adapted from this answer and this answer.

It works like a charm! Many thanks! Just out of curiosity, is this task doable using v9? — Vaggelis_Z, Apr 24 '17 at 19:43
@Vaggelis_Z - the first thing I notice is that you need to use the much slower DeleteDuplicates[data, #1[[;; 2]] === #2[[;; 2]] &] since you don't have DeleteDuplicatesBy — Jason B., Apr 24 '17 at 19:50
checking now - to be honest, using that form of DeleteDuplicates on a list this big was taking so long I abandoned it. — Jason B., Apr 24 '17 at 19:53
RandomPoint and ConvexHullMesh are not recognized in v9. — Vaggelis_Z, Apr 24 '17 at 19:59
@Vaggelis_Z - you can get rid of the duplicates efficiently using positionDuplicates[list_] := Flatten[Rest /@ GatherBy[Range@Length[list], list[[#]] &]]; data = ReplacePart[data, Thread[positionDuplicates[data[[All, ;; 2]]] -> Sequence[]]]; - and the interpolating functions generate just fine. But you don't have access to RandomPoint or ConvexHullMesh, so you need another way to generate the points. — Jason B., Apr 24 '17 at 20:00
@Vaggelis, you're already two versions behind; next time, please mention at the outset that you're using an older version. — J. M.'s missing motivation, Apr 25 '17 at 02:28

george2079 · Answer 2 · 2017-04-25T09:49:47.220

For the random generation part, once you have an interpolation function you can actually just let the interpolation function tell you when points are out of the region:

points = Select[ RandomReal[{-2, 2}, {200, 3}], 
                Norm[#[[1 ;; 2]]] < 1 &];
int = Interpolation[points, InterpolationOrder -> 1]
random = DeleteCases[
   Quiet@Check[Append[#, int @@ #], "Fail"] & /@
    RandomReal[{Min@Flatten[points], Max@Flatten[points]}, {2000, 2}],
     "Fail"];
Show[{
  ConvexHullMesh[points[[All, 1 ;; 2]]],
  Graphics@Point[random[[All, 1 ;; 2]]]}]

note I used ConvexHullMesh here only for the illustration graphic.

Interpolating four-dimensional data

2 Answers2