Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

It seems there are many similar questions on Mathematica.SE, but none addresses (to the best of google) exactly my problem. I have the following data:

L[1] = 
  {{0.0200, {-3.98, -0.902}}, {0.120, {-3.97, -0.890}}, {0.220, {-3.94, -0.860}},
   {0.320, {-3.89, -0.811}}, {0.420, {-3.82, -0.744}},{0.520, {-3.71, -0.659}}, 
   {0.620, {-3.55, -0.554}}, {0.720, {-3.32, -0.429}}, {0.820, {-2.97, -0.285}}, 
   {0.920, {-0.119}}, {1.02, {0.0666}}, {1.12, {0.272}}, {1.22, {0.497}}, {1.32, {0.740}}, 
   {1.42,{1.00, 3.30}}, {1.52, {1.28, 3.68}}, {1.62, {1.58, 4.04}}, {1.72, {1.89, 4.41}}, 
   {1.82, {2.23, 4.79}}, {1.92, {2.57,5.18}}, {2.02, {2.90, 5.58}}, {2.12, {3.38,6.00}},
   {2.22, {3.90, 6.43}}, {2.32, {4.44, 6.88}}, {2.42, {5.00,7.35}}, {2.52, {5.57, 7.84}}, 
   {2.62, {6.15,8.35}}, {2.72, {6.75, 8.87}}, {2.82, {7.36, 9.42}}, {2.92, {7.99,9.98}}};
L[2] = 
  {{0.0200, {0.377, 0.579}}, {0.120, {0.373,0.580}}, {0.220, {0.360, 0.582}}, 
   {0.320, {0.340, 0.585}}, {0.420, {0.334, 0.590}}, {0.520, {0.331,0.596}}, 
   {0.620, {0.336, 0.604}}, {0.720, {0.369,0.613}}, {0.820, {0.501,0.623}}, 
   {0.920,{0.632}}, {1.02, {0.638}}, {1.12, {0.638}},{1.22, {0.635}}, {1.32, {0.629}}, 
   {1.42, {0.622,0.974}}, {1.52, {0.614, 0.944}}, {1.62, {0.601, 0.933}}, 
   {1.72, {0.581, 0.929}}, {1.82, {0.545, 0.933}}, {1.92, {0.472, 0.937}}, 
   {2.02, {0.400, 0.936}}, {2.12, {0.864, 0.934}}, {2.22, {0.889, 0.941}}, 
   {2.32, {0.903, 0.945}}, {2.42, {0.916, 0.949}}, {2.52, {0.929, 0.958}}, 
   {2.62, {0.934, 0.958}}, {2.72, {0.935, 0.949}}, {2.82, {0.943, 0.960}}, 
  {2.92, {0.941, 0.955}}};
L[3] = 
  {{0.0200, {0, 0.0669}}, {0.120, {0, 0.0649}}, {0.220, {0,0.0601}}, 
   {0.320, {0.00730, 0.0528}}, {0.420, {0.0350, 0.0439}}, 
   {0.520, {0.0816, 0.0339}},{0.620, {0.155,0.0237}}, 
   {0.720, {0.297, 0.0141}}, {0.820, {0.727, 0.00620}}, {0.920, {0.00110}}, 
   {1.02, {0.000300}}, {1.12,{0.00460}}, {1.22, {0.0130}}, {1.32, {0.0247}}, 
   {1.42, {0.0395,2.93}}, {1.52, {0.0578, 2.57}}, {1.62, {0.0776, 2.38}}, 
   {1.72, {0.0977, 2.26}}, {1.82, {0.115, 2.18}}, {1.92, {0.123, 2.12}}, 
   {2.02, {0.135, 2.06}}, {2.12, {0.543, 2.01}}, {2.22, {0.697, 1.97}}, 
   {2.32, {0.806,1.94}}, {2.42, {0.891,1.91}}, {2.52, {0.961, 1.90}}, {2.62, 
   {1.01,1.86}}, {2.72, {1.05, 1.82}}, {2.82, {1.08, 1.81}}, {2.92, {1.10,1.77}}};

It is having the following structure

Row[MatrixForm /@ {L[1], L[2], L[3]}]

table plot

That is we have $3$ ragged arrays. Let us plot them:

Do[
 l[i] = Flatten[Table[Table[{a[[1]], b}, {b, a[[2]]}], {a, L[i]}], 1];
 gr[i] = 
   ListPlot[
     l[i], PlotRange -> {{0, 3}, All}, 
     PlotMarkers -> {Graphics@{Black, Disk[]}, 0.02}, 
     FrameLabel -> {"x", "\[CapitalDelta][" <>ToString[i] <> "]"},    
     Axes -> False, Frame -> True, 
     LabelStyle -> Directive[FontFamily -> "Helvetica", FontSize -> 14],
     PlotStyle -> Black], 
  {i, 3}]
GraphicsRow[{gr[1], gr[2], gr[3]}, ImageSize -> Full]

three list plots

These three panels describe 3 different properties of a physical system as a function of one parameter. On each panel there are $3$ curves.

My question is how to connect points belonging to each of $3$ curves. It is certainly possible to edit plots in e.g. adobe illustrator however I am more curios to properly extract the data corresponding to each curve so the parametric plots, e.g., $\Delta[2] (\Delta[1])$ can be created for each curve.

Just to avoid suggestions like Joined -> True or ListLinePlot[], this is not the expected result:

wrong result

Reformulation

Imagine you observe 3 animals in the forrest (rabbit, wolf, fox) as a function of time $(x)$. They have 3 properties (L[1] - smell, L[2] - sound, L[3] - look). Not all animals are close enough to the detector, sometimes detectors register 1, sometimes 2 animals at a time.

The question is:

  • (i) to separate signals belonging to the rabbit, wolf and fox and
  • (ii) to plot the "sound" as a function of "smell" for each of them.

Clarification

As you can see, second columns of data L[1], L[2], and L[3] are lists. In fact they are ordered lists, and the ordering is the same for L[1], L[2], and L[3]. For instance, have a look at the first raw:

{-3.98, -0.902} {0.377, 0.579} {0, 0.0669}

It is known that {-3.98, 0.377, 0} belong to one group and { -0.902, 0.579, 0.0669} to another.

Expected result

There are several very nice answers below. I really appreciate! In order to clarify what exactly I am expecting I re-colored points on the plots manually.

enter image description here

I think the right approach would be to use information from all sets.

share|improve this question
1  
@Coolwater, please read the question carefully – yarchik Feb 9 at 12:12
1  
@Edmund, please read the question carefully – yarchik Feb 9 at 12:12
2  
So I can see three distinct curves in L[1], and I could imagine adapting this answer to find them, but for L[2] and L[3] I do not see three distinct curves. – JasonB Feb 9 at 12:44
1  
Let me see if I am understanding your problem: you're measuring 3 different properties of a physical system, but all your three measurements get mixed up in one bag and you want to separate them afterwards. I m really curious about how comes you can't separate them beforehand. – Dr. belisarius Feb 9 at 13:06
5  
@Dr.belisarius, no it is not exactly how you reformulate it. Sorry for figurative language below. Imagine you observe 3 animals in the forrest (rabbit, wolf, fox) as a function of time (x). They have 3 properties (L[1] - smell, L[2] - sound, L[3] -look). Not all animals are close enough to the detector, sometimes detectors register 1, sometimes 2 animals at a time. The question is (i) to separate signals belonging to rabbit, wolf and fox and (ii) plot sound as a function of smell for each of them. – yarchik Feb 9 at 13:18
up vote 13 down vote accepted

This approach splits each data set in to a set of curves and then attempts to join curves whose end and start points are "close enough". The measure I have used is okay for the example data and it includes an element of rescaling with the curve data but YMMV with "real" data. It is simpler than some of the linked approaches.

Firstly we use a helper function to split the raw data. It enables us to group data points according to their position in the list for each ordinate point.

sowByPosn[{x_, ys_}] := MapIndexed[Sow[{x, #1}, #2[[1]]] &, ys]

We now split the data according to the number of points for each ordinate, then generate and group data points in each sub-list according to their position:

Do[l[i] = 
  Join @@ (Reap[Thread[sowByPosn[#]]][[2]] & /@ 
     SplitBy[L[i], Length[#[[2]]] &]);
 Print@ListLinePlot[l[i], PlotTheme -> "Detailed", 
   PlotLegends -> Automatic], {i, 3}]

Separate curves

There are five curves in each data set and three obviously correspond to a single curve.

To join the curves we compare the distances between the ends of each pair of lines to a limit value calculated from each curve.

The limit is calculated from the mean difference plus 5 times the standard deviation (somewhat arbitrarily):

stdLimit[x_] := 
 Plus @@ ({1, 
     5} Through[{Mean, StandardDeviation}[Norm /@ Differences[x]]])

Generate a table of differences (with mean limits) for each pair of curves and select pairs that have ends closer than the limit.

Table[p[i] = 
  Position[Table[{Norm[Last[c1] - First[c2]], 
     Mean[stdLimit /@ {c1, c2}]}, {c1, l[i]}, {c2, l[i]}], {x_, l_} /;
     x <= l], {i, 3}]

{{{2, 3}, {3, 4}}, {{2, 3}, {3, 4}}, {{2, 3}, {3, 4}}}

Use Mathematica's Graph functionality to group the components (probably we could do this without...) We only need to do this once with the example data -- all have the same connectivity.

components = 
 With[{g = 
    Graph[Join[# -> # & /@ 
       Complement[Range[Length[l[1]]], Union[Flatten[p[1]]]], 
      Rule @@@ p[1]], VertexLabels -> "Name"]}, Print[g];
  (Flatten /@ (EdgeList[Subgraph[g, #]] & /@ 
        WeaklyConnectedComponents[g] /. 
       DirectedEdge[u_, v_] :> {u, v})) /. {a___, b_, b_, c___} :> {a,
      b, c}]

graph

{{2, 3, 4}, {5}, {1}}

Finally join the lines:

Table[z[i] = Join @@ Part[l[i], #] & /@ components; 
  Print@ListLinePlot[z[i], PlotTheme -> "Detailed", 
    PlotLegends -> ToString /@ components], {i, 3}];

Final curves

share|improve this answer
    
I think your result is the most accurate one. Thanks a lot! I need to read your algorithms in order to understand if it is applicable to a larger class of problems. – yarchik Feb 9 at 18:47
    
Thanks, I can see there may be problems if eg two curves end close to a beginning. In that case I think perhaps looking at the gradient may be the way to go. Another minor point is that it reorders the curves. – MikeLimaOscar Feb 9 at 18:58

We can use the following approach: start with an end point that belongs to a path and increment the path with nearest neighbor points that are "good candidates." A point is a "good candidate" if it is not too far away from the last point and it does not produce a sharp turn in the path. To determine "too far" we a look at the distances between the last, say, 5 point pairs.

In order this approach to work the data has to be rescaled and the paths initialized.

This answer is more-or-less complete. It can be refined further with automated rescaling and path initilization.

Below is code that implements the outlined approach.

Clear[GoodPointQ]
GoodPointQ[pointPath_, candidate_, distOutliersSDFactor_: 10, 
   distNLastPoints_: 6] :=
  Block[{dists, threshold, cd, v1, v2},
   (*Path is too short.*)
   If[Length[pointPath] < 3, Return[True]];

   (*Is the distance an outlier?*)

   dists = Map[Norm[#[[1]] - #[[2]]] &, 
     Partition[Take[Reverse[pointPath], UpTo[distNLastPoints]], 2, 1]];
   (*threshold=Median[dists]+(5*Median[Abs[dists-Median[dists]]]);*)

     threshold = 
    Mean[dists] + distOutliersSDFactor*StandardDeviation[dists];
   (*Print["norm:",Norm[
   pointPath\[LeftDoubleBracket]-1\[RightDoubleBracket]-candidate]];
   Print["threshold:",threshold];*)

   If[Length[pointPath] > 4 && 
     Norm[pointPath[[-1]] - candidate] > threshold,
    Return[False]
    ];

   (*Does it make a sharp turn?*)

   v1 = pointPath[[-2]] - pointPath[[-1]];
   v2 = candidate - pointPath[[-1]];
   cd = Dot[v1/Norm[v1], v2/Norm[v2]];
   (*Print["cd:",cd];*)
   If[cd > -10^-4, Return[False]];
   True
   ];

Clear[FindNextGoodPoint]
FindNextGoodPoint[pointPath_, nnFunc_, nCandidates_: 3, 
   sameDist_: 10^-9, distOutliersSDFactor_: 10, distNLastPoints_: 6] :=

    Block[{cs, good},
   cs = nnFunc[pointPath[[-1]], nCandidates + 1];
   cs = Complement[cs, pointPath, 
     SameTest -> (Norm[#1 - #2] <= sameDist &)];
   good = 
    GoodPointQ[pointPath, #, distOutliersSDFactor, 
       distNLastPoints] & /@ cs;
   If[Or @@ good,
    {True, Append[pointPath, Pick[cs, good][[1]]]},
    {False, pointPath}]
   ];

Clear[PathByNextGoodPoint]
PathByNextGoodPoint[startPoints : {{_?NumberQ, _?NumberQ} ..}, 
   allPoints : {{_?NumberQ, _?NumberQ} ..}, nCandidates_: 3, 
   distOutliersSDFactor_: 10, distNLastPoints_: 6] :=

  Block[{pointPath, found = True, nnFunc, k = 0},
   nnFunc = Nearest[allPoints];
   pointPath = startPoints;
   While[found && k < 200,
    k++;
    {found, pointPath} = 
     FindNextGoodPoint[pointPath, nnFunc, nCandidates];
    ];
   pointPath
   ];

Steps of applying this code follow.

Select a dataset.

data = dataOrig = Flatten[Thread /@ L[1], 1];

Normalize the data. (Standardize can be used instead of Rescale.)

Do[
 data[[All, i]] = 
  Rescale[data[[All, i]], MinMax[data[[All, i]]], {0, 1}]
 , {i, {1, 2}}]

This is useful to select paths initial points.

(*Graphics[{Point[data],MapIndexed[Text[#2,#1,{-1,1}]&,data]}]*)

Initial paths of points.

startPoints = {{data[[1]], data[[3]]}, {data[[2]], 
    data[[4]]}, {data[[55]], data[[53]]}};

Find the paths.

paths = PathByNextGoodPoint[#, data, 2] & /@ startPoints;

Single points can be used too.

startPoints = {{data[[1]]}, {data[[2]]}, {data[[55]]}};
paths = PathByNextGoodPoint[#, data, 2] & /@ startPoints;

Rescale to the original data.

pathsOrig = paths;
Do[
  pathsOrig[[j, All, i]] = 
   Rescale[pathsOrig[[j, All, i]], {0, 1}, 
    MinMax[dataOrig[[All, i]]]], {j, Length[paths]}, {i, {1, 2}}];

Plot results.

Show[{ListPlot[dataOrig, 
   PlotStyle -> {PointSize[0.02], GrayLevel[0.8]}], 
  ListPlot[{Sequence @@ pathsOrig}]}, Frame -> True]

enter image description here

Results with the other two sets:

enter image description here

enter image description here

The paths on last image can be somewhat improved by tweaking the parameters of the function PathByNextGoodPoint.

enter image description here

share|improve this answer
    
I am impressed that you could get accurate results working with each set independently. Do you see a possibility to apply your method to triplets of data points? I expect it can improve the accuracy. – yarchik Feb 9 at 18:52
    
@yarchik I need to think more about the problem and my solution. I think the other solutions are much better. As you indicate in your comment, my code is a little too ad-hoc and low level. – Anton Antonov Feb 9 at 21:24

As David G. Stork has pointed out, Nearest Neighbours offer a good method of attack for this problem. Here I've not implemented a full NN chain approach but something a little more basic which gets most of the way there.

I'm using a 'dumb' NN but with a distance function that only allows points to be connected to another point one x-distance away, this works because your data are on a regular grid.

First I just pad the data into a full {x,y} format:

    p1 = Flatten[Transpose /@ ({ConstantArray[#[[1]], Length[#[[2]]]], #[[2]]} & /@ 
  L[1]), 1];
    p2 = Flatten[Transpose /@ ({ConstantArray[#[[1]], Length[#[[2]]]], #[[2]]} & /@ 
  L[2]), 1];
    p3 = Flatten[Transpose /@ ({ConstantArray[#[[1]], Length[#[[2]]]], #[[2]]} & /@ 
  L[3]), 1];

Then using the following distance function:

   distM[{u_, v_}, {x_, y_}] := With[{pen = 10}, 
     If[Round[10 (u - x)] == -1, Abs[v - y], pen]
   ];
   distP[{u_, v_}, {x_, y_}] := With[{pen = 10}, 
     If[Round[10 (u - x)] == 1, Abs[v - y], pen]
   ];

   nf1 = Nearest[p1, DistanceFunction -> distM];

   Graphics[{Point@p1, Line /@ Flatten[Transpose /@ ({ConstantArray[#[[1]], Length[#[[2]]]], #[[2]]} & /@Transpose[{p1, (nf1[#, {1, 1}] & /@ p1)}]), 1]}, opts]

enter image description here

   nf2 = Nearest[p2, DistanceFunction -> distM];

   Graphics[{Point@p2, Line /@ Flatten[Transpose /@ ({ConstantArray[#[[1]], Length[#[[2]]]], #[[2]]} & /@Transpose[{p2, (nf2[#, {1, 0.15}] & /@ p2)}]), 1]}, opts]

enter image description here

   nf3 = Nearest[p3, DistanceFunction -> distP];

   Graphics[{Point@p3, Line /@ Flatten[Transpose /@ ({ConstantArray[#[[1]], Length[#[[2]]]], #[[2]]} & /@Transpose[{p3, (nf3[#, {1, 1}] & /@ p3)}]), 1]}, opts]

enter image description here

   opts = Sequence[Frame -> True, AspectRatio -> 1/GoldenRatio, 

GridLines -> Automatic, PlotRangePadding -> None];

Imperfect but might serve as a starting point and hopefully shows a few of the possibilities.

share|improve this answer
    
Very nice! (+1) If you posted it 1h earlier, I would have not spent time to program the code for my answer. – Anton Antonov Feb 9 at 18:10
    
I think your code is the shortest one. It is also pretty accurate. Here I have the same question as for the approach of @Anton Antonov, can it be applied to the triples of data-points in order to improve the accuracy? – yarchik Feb 9 at 18:55

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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