8
$\begingroup$

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]}]

enter image description here

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!

$\endgroup$
  • $\begingroup$ Are px and py dependent variables? What is the relationship between them? $\endgroup$ – MikeY Apr 24 '17 at 18:14
  • $\begingroup$ @MikeY There is no known relationship between the momenta. The values of the initial data file should be used for the interpolation. $\endgroup$ – Vaggelis_Z Apr 24 '17 at 18:18
  • $\begingroup$ OK, but is it px = px(x,y) and py = py(x,y)? $\endgroup$ – MikeY Apr 24 '17 at 18:25
  • $\begingroup$ @MikeY Perhaps, but I am not sure! $\endgroup$ – Vaggelis_Z Apr 24 '17 at 18:29
8
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ It works like a charm! Many thanks! Just out of curiosity, is this task doable using v9? $\endgroup$ – Vaggelis_Z Apr 24 '17 at 19:43
  • $\begingroup$ @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 $\endgroup$ – Jason B. Apr 24 '17 at 19:50
  • $\begingroup$ Is this the only change in v9? $\endgroup$ – Vaggelis_Z Apr 24 '17 at 19:52
  • $\begingroup$ checking now - to be honest, using that form of DeleteDuplicates on a list this big was taking so long I abandoned it. $\endgroup$ – Jason B. Apr 24 '17 at 19:53
  • $\begingroup$ RandomPoint and ConvexHullMesh are not recognized in v9. $\endgroup$ – Vaggelis_Z Apr 24 '17 at 19:59
2
$\begingroup$

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]]]}]

enter image description here

note I used ConvexHullMesh here only for the illustration graphic.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.