8
$\begingroup$

Question How can I find the closest point to a parametric curve produced by a numerical method.

General Context I have a parametric curve (which I will call $g$) produced from a averaging many solutions of an ODE system. As part of further analysis I would like to construct a function that given a point outside the curve, it gives the closest point on the parametric curve $g$.

i.e. I am looking for a function $f:x_0\to p$ where $x_0\in \mathbb{R}^3_+$, and $p$ is a point on the parametric curve $g:t \to (x,y,z)$ such that $||x_0-p||$ is minimized (Where the norm is the standard euclidean norm).

Why do I want such a function? I want to compare the distance of another parametric curve (which I will call $h$) to $g$ at a given $t$. Or put another way for a time $t$ what is $||h(t)-p||$. (Note that I am not trying to do what it done in this post which is trying to find $||h(t)-g(t)||$. I am looking for the point on $g$ which is closest to $h(t)$, anywhere on $g$.)

My biggest problem is using Mathematica to construct an appropriate $f$.

Problems Using Mathematica

There are quite a few Mathematica functions that I have looked at, and I believe that the the appropriate answer is just a question of calling the right functions with proper syntax. Possible Mathematica functions include

  • RegionNearest
  • NMinimize
  • ParametricRegion
  • ImplicitRegion

the idea being that we can use either ImplicitRegion or ParametricRegion to define a region which is all of $g$, then use either RegionNearest or NMinimize. As an example

RegionNearest[ParametricRegion[{Cos[theta], Sin[theta]}, {{theta, 0, 2 \[Pi]}}], {2,2}]

enter image description here

I haven't been apply to figure the right combination functions and syntax though.

Minimum Working Example Note there are two separate sections. One where $g$ is generated and another where I test different options. I have included the way in which $g$ is generated is as there may be a syntax problems, in the way $g$ is generated. Otherwise you may consider $g$ a blackbox.

(*Simulation Parameters and definitions*)
Clear[f]
Clear[i, P, B]
f[P_, B_] := 1/2 P + 10 B/(1 + B);
tmax = 20;
randNum = 50;
A = {{1/20, 1/4, 1/50}, {1/4, 1/26, 1/40}};
point = {16.666666666666735`, 0.`, 8.333333333333345`};

(*ODE System*)
ODEsys = {i'[t] == f[P[t], B[t]] - i[t],
   P'[t] == 
    P[t] (1 - A[[1, 1]] P[t] - A[[1, 2]] B[t] - A[[1, 3]] i[t]),
   B'[t] == B[t] (1 - A[[2, 2]] B[t] - A[[2, 3]] i[t])};

(*Generating parametric g curve.*)
(*--------------------------*)
eps = 0.001;
set = List[];
While[Length[set] < randNum,
 holdSet = 
  Join[set, Map[point + # &, RandomReal[{-eps, eps}, {randNum, 3}]]];
 set = Select[holdSet, #[[2]] >= 0 &];
 ]

set = Drop[set, -(randNum - Length[set])];

(* Simulation *)
sol = ParametricNDSolveValue[{ODEsys, {P[0] == init1, B[0] == init2, 
     i[0] == init0}}, {P, B, i}, {t, 0, tmax}, {init1, init2, init0}];

(* Averging solution over multiple inital conditions. *)
gCurve[t_] := 
  Evaluate[Mean[
    Through /@ (sol[#[[1]], #[[2]], #[[3]]][t] & /@ set)]] ;
(* Below is a test that gCorve works as intended. If gcurve works as \
intended we should get a single curve in 3D *)
simplexPlot = 
 ParametricPlot3D[gCurve[t], {t, 0, tmax}, PlotRange -> All, 
  ImageSize -> Large, PlotStyle -> Black] 

(*-------------------------*)
(*Attempt to solve the problem*)
(*-------------------------*)

(* Attempt *)
(* Problem: gives out a function and not a number *)
nearestPoint[{x_, y_, z_}] := 
  Evaluate[RegionNearest[
    ParametricRegion[gCurve[t], {{t, 0, tmax}}], {x, y, z}]];
nearestPoint[{2, 2, 2}]

Notes

If you need any clarification don't be afraid to ask.

Some previous posts that I've seen but have answered my question.

$\endgroup$
  • 3
    $\begingroup$ Nicely written question. randNum and simplexSol are not defined as far as I can tell though. $\endgroup$ – C. E. Feb 7 at 18:30
  • $\begingroup$ Apologies I copied the code from a match larger document. Forgot to clear memory to test it. The new edited code should work as intended. $\endgroup$ – AzJ Feb 7 at 18:37
5
$\begingroup$

It seems like the most straightforward way to handle this is to just use NMinimize directly on the norm. I do this here by defining a function to by find the t value which minimizes the distance between the function and the point:

nearestT[f_, pt_, tmin_, tmax_] := 
  Module[{t}, 
   t /. NMinimize[{Norm[f[t] - pt], tmin <= t, t <= tmax}, t][[2]]
  ];

And then subsequently a function which inserts that t value back into the original function to find the nearest point on the curve:

nearestPt[f_, pt_, tmin_, tmax_] := f[nearestT[f, pt, tmin, tmax]];

For gCurve, this is used as:

{nearestT[gCurve, {2, 2, 2}, 0, tmax], nearestPt[gCurve, {2, 2, 2}, 0, tmax]}

{12.9809, {4.25274, 6.24997, 10.07}}

If it's too slow and you suspect that gCurve will be relatively well behaved with regards to the point you're minimizing at, you may wish to use FindMinimum instead. However, FindMinimum does not make nearly as many attempts to find global minima, so it may not always provide the best answer.

$\endgroup$
  • $\begingroup$ I tested your code a little bit for robustness and it seems to work. I did not know that FindMinimum would work given that the function gcurve is not a standard function. I really like your approach. I am keeping the question open for a while to see if people have any alternative answers. If no one answers this question in 24 hours I will accept your answer. $\endgroup$ – AzJ Feb 7 at 19:23
  • 1
    $\begingroup$ @AzJ It is true that there can be some odd issues here and there when using non-standard function forms, but the numerical functions tend to be pretty uncaring about the underlying nature of the function so long as the objective function you're working with is numerical by the time it's fully evaluated. With some functions, that might involve redefining things to use x_?NumericQ and such, but once that's done they're usually quite reliable. $\endgroup$ – eyorble Feb 7 at 19:30
7
$\begingroup$
nF = RegionNearest[Cases[simplexPlot, _Line, All][[1]]];
pnt = {2, 2, 2};
npnt = nF @ pnt

{4.56162,6.16536,10.9321}

EuclideanDistance[pnt, nF@pnt]

10.1831

Show[simplexPlot, 
 Graphics3D[{Red, Sphere[pnt, .2], Green, Sphere[npnt, .2], 
   Gray, Dashed, Arrow[{pnt, npnt}]}]]

enter image description here

$\endgroup$
  • $\begingroup$ I like how elegant your method is. Could you explain a little bit how Cases[simplexPlot, _Line, All] works? $\endgroup$ – AzJ Feb 7 at 23:24
  • 1
    $\begingroup$ @AzJ, Cases[expr, _h, All] extracts objects with Head h that appear at any Level in expr. If you look at InputForm[simplexPlot] it is a collection of directives and a single primitive Line[...] and Cases[simplexPlot, _Line, All] gets that line expression. $\endgroup$ – kglr Feb 7 at 23:38
  • $\begingroup$ Okay that makes sense. Does this mean that the accuracy of RegionNearest[Cases[simplexPlot, _Line, All][[1]]] depend on how accurate the graphing method is? i..e. Changing MaxRecursion , PerformaceGoal,WorkingPrecision changes the accuracy of the method. $\endgroup$ – AzJ Feb 7 at 23:50
  • 1
    $\begingroup$ @AzJ, good point re accuracy. In addition to the options you mentioned PlotPoints also affects the Line produced by the plot function. $\endgroup$ – kglr Feb 7 at 23:55

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.