# How to extract different curves from a mixed dataset?

Suppose I have a dataset:

testdata =
SortBy[#, #[[1]] &] &@
Flatten[{{#, #/(1 + Exp[#/20] // Log) + RandomReal[0.3]} & /@
Range[100], {#, 10 + 20/#^0.3 + RandomReal[0.3]} & /@
Range[100]}, 1];


ListPlot[testdata] looks like:

We can see it contains "two curves". If I have the testdata, then how can I separate these two curves? How can I extract them from the mixed dataset?

The dataset above is a simple example. In reality, the mixed dataset comes from the solution of a complicated dispersion equation which is hard to factorization.

• I am wondering, if the actual data points are equally dense and have similar variance per curve. Is it a problem of isolating a linear and a quadratic component or are there perhaps points that can profitably use other fitting functions? Here is the result of my first attempt but it depends on the information I am seeking through this comment.
– Syed
Oct 15, 2021 at 13:38
• If the points of your 2 curves are adjacent to each other like in the example, you may write a simple Table command to extract the curves. E.g. first curve: Table[testdata[[i]], {i, 1, Length[testdata], 2}] and second curve: Table[testdata[[i]], {i, 2, Length[testdata], 2}] Oct 15, 2021 at 14:06
• @DanielHuber No..that's cheating...the real testdata I have to deal with is unordered and contains more than two "curves". Oct 15, 2021 at 16:57
• @Syed That's really nice... The actual data points are equally dense and data error is very small. however, I don't have prior knowledge about what fitting functions should be used (Though we can guess it when I plot the dataset out). "linear and a quadratic" is just an oversimplified example. Oct 15, 2021 at 17:07
• To extract near-collinear points, it is important to have a linear component in the data. If there are oscillations or interference from other polynomials or saturation artifacts, things will only go downhill from there. Is there any chance you can load a part/picture of the data so that readers can gauge their options?
– Syed
Oct 15, 2021 at 17:19

This is a very crude, inefficient and not very Mathematica-like answer, meant only to provide an algorithmic example of approaching this problem.

This algorithm can be thought of as some kind of an "Euler-like integrator".

1. Pick a random point $$r_0$$ from the convex hull
2. Find the nearest point $$r_1$$
3. Approximate the derivative $$\Delta r \approx r_1 - r_0$$
4. Estimate the next point by $$r_3^* = r_0 + 2\Delta r$$
5. Find $$\tilde{r}$$ as the nearest point to the $$r_3^*$$:
• If this point is close enough ($$|\tilde{r}-r_3^*|<\varepsilon$$), let $$\tilde{r}$$ become $$r_3$$. Go to 3.
• Otherwise, complete the curve. Go to 1.

The interpolation could be also of higher order to provide a better estimate $$r_3^*$$.

Note 1: The algorithm doesn't complete the curve immediately, but first tries to continue it at it other end.

Note 2: There is some randomness to it, and some runs may give worse results.

Note 3: Up to 2 points can miss in the curves.

findNewStartingPoint[data_ /; Length@data <= 2] = {0, 0};
findNewStartingPoint[data_] := Module[{hull, first, second, newdata},
hull = MeshCoordinates[ConvexHullMesh[data]];
first = First@RandomSample[hull, 1];
newdata = Delete[data, Position[data, first]];
second = First@Nearest[newdata -> {"Index", "Element"}, first];
newdata = Delete[newdata, second[[1]]];
{{first, second[[2]]}, newdata}
];

findCurves[pts0_, eps_] :=
Module[{curves, currCurve, pts, pnt, direction, dr},
direction = 1;
curves = {};
{currCurve, pts} = findNewStartingPoint[pts0];

While[Length@pts > 0,
dr = currCurve[[-1]] - currCurve[[-2]];
pnt = First@Nearest[pts -> {"Index", "Element", "Distance"},
currCurve[[-2]] + 2 dr
];

If[(pnt[[3]] > eps && direction == -1) || Length@pts == 1,
(* Curve completed, start a new one *)
AppendTo[curves, currCurve];
{currCurve, pts} = findNewStartingPoint[pts];
direction = 1
,
If[pnt[[3]] > eps && direction == 1,
(* Continue at the other end *)
direction = -1;
currCurve = Reverse@currCurve;
,
(* Add point to current curve *)
AppendTo[currCurve, pnt[[2]]];
pts = Delete[pts, pnt[[1]]];
]]];
curves
];

curves = findCurves[N@testdata, 3];
ListPlot[curves]


testdata2 =
N@RandomSample[#, Length@#] &@
Catenate[
Table[{i, # + RandomReal[.5]}, {i, 0, 100}] & /@ {i - 20, i^2/200,
200/(i + 5), -i^2/180 + 50, 70 - i, 10 Sin[i/5] + 45}];
ListPlot[findCurves[N@testdata2, 3], Axes -> False,
FrameTicks -> False, Frame -> True, AspectRatio -> 1] & /@ Range[3] // Row


• This is really nice already....thank you very much! Oct 16, 2021 at 6:15
• @Domen, in step 5 of the algorithm, "If this point is close enough ($|\tilde{r}-r_3^*|<\varepsilon$), let $\tilde{r}$ become $r_3$", should it be "let $\tilde{r}$ become $r_1$ and $r_1$ become $r_0$"?
– lxy
Mar 1, 2023 at 14:40
• @Domen, btw, in your algorithm there is no $r_3$ only $r_3^\ast$
– lxy
Mar 2, 2023 at 1:36