7
$\begingroup$

I'm plotting several lists of data (constant lines of density, temperature quality, etc. on a pressure-enthalpy diagram), none of which are fitted well by Interpolation or Fit. Since each table/line represents a constant density (in this case), I want to determine what density is at any selected point based off of the two nearest density lines (blue).

enter image description here

  • Blue points are data points in rho lists.
  • Orange points are all points from my function TotalClose which returns one point from each rho list.
  • Green points are from my function PointClose.
  • Largest black point is set with Manipulate controls.

The Locator is on a log scale so pt makes it linear for calculations:

Control[{{pt0,{2295,Log@145}},Locator}]

pt={First@pt0,Exp@Last@pt0};

In my attempt I found the point from each rho list closest to the point pt selected with the Locator:

TotalClose[list_]:=First@Nearest[list[[#]],pt]&/@Range@Length@list;

PointClose[list_]:=Nearest[TotalClose[list],pt,2];

Functions are evaluated at {rho1000,rho700,rho50,rho2,rho02}

The green points PointClose returned are not correct here. It needs to return points from the rho700 and rho50 lists. This would be easier if I could fit these lists to a function but they won't.

I was planning on using the lever rule to determine what the density would be at the selected point. I am stuck as to how to get which two rho lists are either immediately to the left and right or above and below pt.

I tried changing the DistanceFunction in Nearest with no luck:

 DistanceFunction->(((#1[[1]]-#2[[1]])^2+(Log@#1[[2]]-Log@#2[[2]])^2)^1/2&)

Here are the rho lists that make the blue constant density lines:

rho1000 = {{2674.112, 10000.}, {2442.182, 9000.}, {2211.362, 
    8000.}, {1980.958, 7000.}, {1750.036, 6000.}, {1517.27, 
    5000.}, {1280.666, 4000.}, {1036.931, 3000.}, {779.806, 
    2000.}, {493.584, 1000.}, {461.831, 900.}, {429.132, 
    800.}, {395.296, 700.}, {360.063, 600.}, {323.059, 
    500.}, {283.706, 400.}, {241.029, 300.}, {193.124, 200.}, {135.03,
     100.}, {120.986, 80.}, {105.405, 60.}, {87.389, 40.}, {64.549, 
    20.}, {48.572, 10.}, {44.443, 8.}, {39.658, 6.}, {34.028, 
    4.098}, {29.72, 2.984}, {25.44, 2.181}, {21.188, 1.692}, {16.964, 
    1.521}, {16.964, 0.05}};

rho700 = {{5000., 8801.33}, {4642.4, 8000.}, {4207.4, 7000.}, {3784.3,
     6000.}, {3371.92, 5000.}, {2968.32, 4000.}, {2570.4, 
    3000.}, {2173.14, 2000.}, {1767.73, 1000.}, {1726.18, 
    900.}, {1684.37, 800.}, {1642.24, 700.}, {1599.79, 
    600.}, {1556.95, 500.}, {1513.7, 400.}, {1469.97, 300.}, {1425.7, 
    200.}, {1380.8, 100.}, {1377.65, 93.070}};

rho50 = {{5000., 331.29}, {4673.64, 300.}, {3703.16, 200.}, {2815.85, 
    100.}, {2739.86, 91.84}, {2444, 71}, {2100, 50}, {1900, 
    37}, {1700, 26}, {1490, 15.5}, {1300, 9.5}, {1050, 4.1}, {820, 
    1.6}, {590, 0.5}, {425, 0.2}, {325, 0.1}, {250, 0.05}};

rho2 = {{5000., 13.07}, {4184.58, 10.}, {3693.09, 8.}, {3234.16, 
    6.}, {2802.61, 4.}, {2734.26, 3.68}, {2500, 3.3}, {2150, 
    2.5}, {1950, 2}, {1750, 1.5}, {1500, 1}, {1230, 0.6}, {1020, 
    0.38}, {815, 0.2}, {660, 0.115}, {450, 0.05}};

rho02 = {{5000., 1.31}, {4184.66, 1.}, {3692.28, 0.8}, {3232.73, 
    0.6}, {2803.20, 0.4}, {2626.43, 0.31}, {2500, 0.32}, {2250, 
    0.32}, {2000, 0.3}, {1840, 0.28}, {1660, 0.24}, {1480, 
    0.2}, {1250, 0.15}, {1000, 0.105}, {810, 0.08}, {550, 0.05}};
$\endgroup$
  • 1
    $\begingroup$ Why can't you use Interpolation of the data sets? You should be able to. $\endgroup$ – march Nov 11 '15 at 17:29
  • $\begingroup$ For instance, this was generated by interpolating your points. $\endgroup$ – march Nov 11 '15 at 17:35
  • $\begingroup$ Wow I could never get the curves to look that nice! How exactly did you get that? Some of the other data I'm plotting looks even worse with Interpolation, that and some don't pass the vertical line test $\endgroup$ – baumannr Nov 11 '15 at 18:01
  • $\begingroup$ f = Interpolation /@ Transpose @ rho1000, for instance. This returns a list of two pure interpolating functions that you can plot parametrically: ParametricPlot[Through@f@t, Evaluate@Flatten@{t, f[[1]]["Domain"]}, AspectRatio -> 1]. I've been trying to figure out how to do the calculation of the nearest points on the curves, but haven't succeeded yet. Perhaps you can have a go with these interpolations. $\endgroup$ – march Nov 11 '15 at 18:04
  • 1
    $\begingroup$ @march : VoronoiMesh[] should give you the nearest point (which can be mapped to its curve via an association). Dropping out the winning curve's points and VoronoiMesh[]ing the surviving points will give the closest point on the second closest curve. There are only five curves, so precomputing and storing the six meshes should be pretty easy... (I'm far from a Mathematica currently, or I'd just whack this up myself.) $\endgroup$ – Eric Towers Nov 12 '15 at 0:31
6
$\begingroup$

Update

I don't know what the "lever rule" is, even after visiting the Wikipedia page on it, but here's the routine wrapped up into a Module.

Clear[interpolateDensity]
interpolateDensity[pnt : {_, _}, fs : {{_InterpolatingFunction, _InterpolatingFunction} ..}, dens_List] /; Length@dens == Length@fs := Module[
  {mins}
  , mins = NMinimize[{#.# &@(pt - Through@#@t), 1. <= t <= #[[1]]["Domain"][[1, 2]]}, t] & /@ fs
  ; Dot @@ Transpose@#/Plus @@ #[[All, 1]] &@ Sort[Transpose[{mins[[All, 1]], dens}]][[1 ;; 2]]
 ]

This calculates the density by finding the nearest points on the nearest curves and using the relative distances as weights for the densities of those curves. Usage:

Clear[pt, fs, dens]
pt = {3604, 693};
fs = Interpolation /@ Transpose@# & /@ {rho02, rho2, rho50, rho700, rho1000};
dens = {0.2, 2, 50, 700, 1000};
interpolateDensity[pt, fs, dens]
(* 18.7113 *)

Alternatively, here is a routine that returns the list of the nearest points along with the densities corresponding to those curves (which you can then use to get the density via the lever rule?):

nearestPoints[pnt : {_, _}, fs : {{_InterpolatingFunction, _InterpolatingFunction} ..}, dens_List] /; Length@dens == Length@fs := Module[
  {mins}
  , mins = NMinimize[{#.# &@(pt - Through@#@t), 1. <= t <= #[[1]]["Domain"][[1, 2]]}, t] & /@ fs
  ; Sort[MapThread[{#2[[1]], Through@#1@t /. #2[[2]], #3} &, {fs, mins, dens}]][[1 ;; 2]][[All, 2 ;; 3]]
  ]

Usage:

nearestPoints[pt, fs, dens]
(* {{{3658.84, 195.126}, 50}, {{3606.93, 7.63443}, 2}} *)

Original Post

For starters, I will not re-paste in your data, but I will define a list of densities:

dens = {0.2, 2, 50, 700, 1000};

Here's something that I think will work. We first carefully Interpolate the data by interpolating the x and y lists separately to create a parametric interpolating function, like so:

fs = Interpolation /@ Transpose@# & /@ {rho02, rho2, rho50, rho700, rho1000};

These functions will be defined on the intervals 1 to Length@rho:

domains = Flatten /@ Through[Transpose[fs][[1]]@"Domain"]
(* {{1., 16.}, {1., 16.}, {1., 17.}, {1., 19.}, {1., 32.}} *)

We can plot the InterpolatingFunctions to see that they look nice:

p1 = Show[
 ListLogPlot[{rho02, rho2, rho50, rho700, rho1000}]
 , ParametricPlot[Evaluate[{#1, Log@#2} & @@ (Through@#@t)], Evaluate@Flatten@{t, #[[1]]["Domain"]}] & /@ fs
]

enter image description here

We now go through the process of minimizing the Euclidean distance between some point

pt = {RandomReal[{0, 5000}], RandomReal[{10, 10^3}]}
(* {3604.77,693.335} *)

and the curves using the code:

mins = NMinimize[{#.# &@(pt - Through@#@t), 1. <= t <= #[[1]]["Domain"][[1, 2]]}, t] & /@ fs;

We are minimizing the square of the Euclidean distance here, but putting the square root around #.# doesn't seem to spit out any errors, so it should work just as well. Obviously, if you think some other distance function is more relevant, then here's where to modify the code. For illustration purposes, we extract the points on the curves as

pts = MapThread[Through@#1@t /. Last@#2 &, {fs, mins}]
(* {{3605.06,0.763082},{3607.7,7.63772},{3659.63,195.214},{1921.07,1376.9},{649.59,1537.91}} *)

and plot everything together:

Show[
  p1
  , ListLogPlot[{pt}, PlotStyle -> {Black, PointSize[0.02]}]
  , ListLogPlot[pts, PlotStyle -> {Red, PointSize[0.02]}]
 ]

enter image description here

Finally, we choose the closest two points, keeping their distances and their densities:

weights = Sort[Transpose[{mins[[All, 1]], dens}]][[1 ;; 2]]
(* {{251135., 50}, {470189., 2}} *)

And finally, we linearly interpolate the densities by weighting the densities by the relative distances:

density = Dot @@ Transpose@#/Plus @@ #[[All, 1]] &@weights
(* 18.7116 *)

All of this can be automated into a neat little Module, and of course the code as written isn't the nicest to look at, but I think this works! Take it and run with it.

$\endgroup$
4
$\begingroup$
fs = Interpolation /@ {Mean /@ GatherBy[rho1000, First], rho700, rho50, rho2, rho02};



Show[
 ListLogPlot[{Mean /@ GatherBy[rho1000, First], rho700, rho50, rho2, 
   rho02}
  , PlotRange -> All
  , PlotTheme -> "Scientific"
  ]
 , Sequence @@ Table[
   With[{f = Evaluate[fs[[k]]]},
    LogPlot[f[t], Evaluate[{t, Sequence @@ First[f["Domain"]]}]]
    ], {k, Length[fs]}]
 ]

Mathematica graphics

minpf[pnt_, func_] := Block[{r, x, xr, y},
  {r, xr} = NMinimize[
    {
     Norm[{pnt, {t, func[t]}}]
     , Between[t, func["Domain"]]
     }, t];
  x = Part[xr, 1, 2];
  {r, x, func[x]}
  ]

minpf[pnt_, func_] := Block[{r, x, xr, y},
  {r, xr} = NMinimize[
    {
     Norm[{pnt, {t, func[t]}}]
     , Between[t, func["Domain"]]
     }, t];
  x = Part[xr, 1, 2];
  {r, x, func[x]}
  ]

findpnt[pnt_] := Rest@First@MinimalBy[Map[minpf[pnt, #] &, fs], First]

findpnt[{100, 1}]
(* {16.964, 0.7855} *)
$\endgroup$
4
$\begingroup$

Following my OP comment, the first thing that worked. (Optimized? Heck, it barely works. I considered using the rho-sets as their labels, but that just seemed like nuking flies, so string labels.)

Module[{rhos, logRhos, minX, maxX, minY, maxY},
  rhos = {rho1000, rho700, rho50, rho2, rho02};
  logRhos = {#[[1]], Log[10, #[[2]]]} & /@ # & /@ rhos;
  {minX, maxX} = MinMax[(Join @@ logRhos)[[All, 1]]];
  {minY, maxY} = MinMax[(Join @@ logRhos)[[All, 2]]];
  {lrs1000, lrs700, lrs50, lrs2, lrs02} = (# - {minX, minY})/({maxX, maxY} - {minX, minY}) & /@ # & /@ logRhos;
  lrsPointToCurve = Association[Sequence[
    (# -> "lrs1000" &) /@ lrs1000,
    (# -> "lrs700" &) /@ lrs700,
    (# -> "lrs50" &) /@ lrs50,
    (# -> "lrs2" &) /@ lrs2,
    (# -> "lrs02" &) /@ lrs02]];
  logRhosScaled = Join[{lrs1000, lrs700, lrs50, lrs2, lrs02}];
  pointToFirstLRSPoint = Nearest[Join @@ logRhosScaled];
  knockout["lrs1000"] = Join[{lrs700, lrs50, lrs2, lrs02}];
  knockout["lrs700"] = Join[{lrs1000, lrs50, lrs2, lrs02}];
  knockout["lrs50"] = Join[{lrs1000, lrs700, lrs2, lrs02}];
  knockout["lrs2"] = Join[{lrs1000, lrs700, lrs50, lrs02}];
  knockout["lrs02"] = Join[{lrs1000, lrs700, lrs50, lrs2}];
  pointToSecondLRSPoint[pt_] := Block[{},
    Nearest[Join @@ knockout[
      lrsPointToCurve[pointToFirstLRSPoint[pt][[1]]]]][pt]];
]

ListPlot[logRhosScaled]
(* Point plot.  Happily, where they're expected to be. *)

VoronoiMesh[Join @@ logRhosScaled]
(* Pretty blue cellulation. *)

pointToFirstLRSPoint[{0.5, 0.5}]
lrsPointToCurve[pointToFirstLRSPoint[{0.5, 0.5}][[1]]]
pointToSecondLRSPoint[{0.5, 0.5}]
lrsPointToCurve[pointToSecondLRSPoint[{0.5, 0.5}][[1]]]
(* {{0.48706, 0.594656}} *)
(* lrs50 *)
(* {{0.559026, 0.359004}} *)
(* lrs2 *)

This strictly meets the request : go from an arbitrary point to the two nearest curves. The same sort of idea (use an Association[]) will let us go from logRhoScaled points to original points if that was actually intended.

$\endgroup$
  • $\begingroup$ I like it. If for no other reason than for the fact that it's a very different take on the problem. $\endgroup$ – march Nov 12 '15 at 6:26

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.