Math 4500 Final Project 2013
Navigation on a Landscape

Defining the model.

For this project, we’re going to try to find geodesics on a landscape z = f(x,y) using numerical methods in a few different ways. The first problem is finding a way to specify the path. The idea is that we’ll let the variables in the problem be a set of waypoints in the form

{x1,y1}, {x2,y2}, ... , {xn,yn}

We need a way to specify a path through these points which 

a) passes through the start (at t=0), all these points, and the end (at t=1)
b) is differentiable

We’ll use Mathematica’s Interpolation to do this because it’s actually relatively easy. Other options would include a piecewise-linear path (just connect the dots with straight line segments) and the BezierFunction call. There’s a real hazard here!  Because of the way we’re numerically integrating, the relationship between the parametrization of the curve (ideally, by arclength), and the parameter (t) has to be constructed a bit carefully. We do so by spacing the t values according to the Euclidean distance between the input points.

The first step here is to figure out a list of (t,x) coordinates and a list of (t,y) coordinates to interpolate.

PathFunction[Start_,Points_,End_] :=  Module[{AllPoints,tVals,Distances},
AllPoints =Join[{Start},Points,{End}];
Distances = (Norm /@ (AllPoints - RotateRight[AllPoints]))[[2;;]];
tVals =Join[{0},Accumulate[Distances]/ Total[Distances]];
AllPoints = Transpose[AllPoints];
Interpolation[#,InterpolationOrder->1]& /@ (Transpose[{tVals,#}]& /@ AllPoints)
];

Let’s do a quick example of a path from {1,1} to {10,10}, where we can see the nodes and move them.

PathFunction[{1,1},{{2,4},{7,5}},{10,10}]

Manipulate[
MyPath=PathFunction[{1,1},pts,{10,10}];
γ[PathFunction_,t_] := {PathFunction[[1]][t],PathFunction[[2]][t]};	
ParametricPlot[γ[MyPath,t],{t,0,1},AxesOrigin->{0,0},PlotRange->{{-1,11},{-1,11}}],
{{pts,{{3,6},{5,6},{7,2}}},Locator,LocatorAutoCreate->True}
]

Now we need a demonstration landscape:

f[x_,y_] := 1 - 0.3(x^2 + y^2);
gradF[{x_,y_}] := {-0.6x,-0.6y};

We are now ready to work on the function which gives you the length of the path from A to B.

PathLength[f_,gradF_,A_,Points_, B_] := Module[{γ,PathFunc,dp1,dp2,dγ},
PathFunc = PathFunction[A,Points,B];
dp1 = Derivative[1][PathFunc[[1]]]; dp2 = Derivative[1][PathFunc[[2]]];
γ[t_] := {PathFunc[[1]][t],PathFunc[[2]][t]};
dγ[t_] := {dp1[t],dp2[t]};
NIntegrate[Sqrt[dγ[t].dγ[t] + (dγ[t].gradF[γ[t]])^2],{t,0.0,1.0},AccuracyGoal->4,PrecisionGoal->4,Method->"LocalAdaptive"]
];

PathLength[f,gradF,{-2,-2},{{1,1},{2,0},{2,2.5}},{3,3}] // Timing

For debugging, it’s actually helpful to be able to look at the profile of a curve along a path as well. This is easiest to encapsulate, since we have to think about the order of operations for performance reasons.

PathPlot[f_,A_,Points_,B_] := Module[{γ,t,PathFunc},
PathFunc = PathFunction[A,Points,B];
γ [t_] := {PathFunc[[1]][t],PathFunc[[2]][t]};
Plot[f @@ γ[t],{t,0,1},Filling->Axis]
];

PathPlot[f,{-2,-2},{{1,1},{2,0},{2,2.5}},{3,3}]

PathIntegrandPlot[gradF_,A_,Points_,B_] := Module[{γ,t,PathFunc},
PathFunc = PathFunction[A,Points,B];
γ[t_] := {PathFunc[[1]][t],PathFunc[[2]][t]};
Plot[Sqrt[γ'[t].γ'[t] + (γ'[t].gradF[γ[t]])^2],{t,0.0,1.0},Filling->Axis]
];

PathIntegrandPlot[gradF,{-2,-2},{{1,1},{2,0},{2,2.5}},{3,3}]

The Penalty Function

The waypoints seem to bunch up in a very weird way during the minimization. In order to prevent this from happening, we’ll add an artificial “penalty term” in order to keep the waypoints apart from one another.

PenaltyFunction[f_,A_,Pts_,B_] := Module[{AllPoints,Distances,tVals},
AllPoints = Join[{A},Pts,{B}];
AllPoints = {#[[1]],#[[2]],f[#[[1]],#[[2]]]} & /@ AllPoints;
Distances = (Norm /@ (AllPoints - RotateRight[AllPoints]))[[2;;]];
Distances /= Total[Distances];
(0.1*Norm[Table[1/Length[Distances],{Length[Distances]}] - Distances])^5
];

PenaltyFunction[Sin[#1 * #2]&,{1,1},{{2,2},{3,3}},{4,4}]

PenaltyFunction[{1,1},{{2,5},{3,3}},{4,4}]

Feeding the model to NMinimize.

NMinimize expects the input function to simply take a collection of arguments and for those arguments to be variables with explicit names. This would be a pain to write out by hand, but we can use the amazing ability of Mathematica to generate Mathematica code to save a lot of time.

Basically, we want lists of variables in the form 

{x1,y2},{x2,y2}, ..., {xn,yn} 

for any n. We’ll generate these as strings using a Table:

Vars[nVectors_] := "{" <> StringJoin[Riffle[Table["{" <> Riffle[# <> ToString[i] &  /@ {"x","y"},","] <> "}",{i,1,nVectors}],","]] <> "}";

Vars[4]

We’ll also need to be able to write these as (flat) inputs:

VarPattern[nVectors_] := StringJoin[Riffle[Flatten[Table[# <> ToString[i] <> "_" &  /@ {"x","y"},{i,1,nVectors}]],","]];

VarPattern[4]

VarsFlat[nVectors_] := StringJoin[Riffle[Flatten[Table[# <> ToString[i] &  /@ {"x","y"},{i,1,nVectors}]],","]];

VarsFlat[5]

And as variables with bounds given by a search rectangle {{xmin,xmax},{ymin,ymax}}.

NMinPattern[n_,SearchRect_] := "{" <> StringJoin[Riffle[Flatten[Table[("{" <> #[[1]] <> ToString[i] <> "," <> ToString[SearchRect[[#[[2]],1]]] <> "," <> ToString[SearchRect[[#[[2]],2]]] <> "}" & /@ {{"x",1},{"y",2}}),{i,1,n}]],","]] <> "}";

NMinPattern[2,{{1,3},{2,4}}]

Now generating and running the NMinimize call is a matter of putting all of this together. First, we need to define the “dummy” function that we’ll pass to NMinimize. 

The point is that the dummy function must hide the dependence on the endpoints A and B and on the functions f and gradF, so we’ll need to make an internal definition here. 

GOTCHA: Since these are the names of local variables in the Module, Mathematica is actually renaming them with some completely weird name like “A$1322” which we can’t know in advance. So we have to refer to them by ToString[A] instead of just putting “A” in the string we’re creating, in order to force Mathematica to substitute whatever weird string it’s using internally into the right place in our function definition. We have to do the same thing with “f” and “gradF” in order to make sure we’re referencing the LOCAL copies of f and gradF passed in as inputs instead of any globals with the same names.

LandscapePath[f_,gradF_,A_,B_,SearchRectangle_,n_] := Module[{DefineDummyFunction,DefineNMinCall,NMinCall,InitialPath,t,DummyFunction,Solution,TestResults,RawSol},

DefineDummyFunction = ToString[DummyFunction]<>"[" <> VarPattern[n] <> "] := PathLength[" <> ToString[f] <> ","<> ToString[gradF]<>"," <> ToString[A] <>"," <> Vars[n] <> "," <> ToString[B] <> "] + PenaltyFunction[" <> ToString[f] <> ","<>ToString[A] <>","<> Vars[n] <> "," <> ToString[B] <>"];";

Print[DefineDummyFunction];

InitialPath = Flatten[Table[(1-t)A + t B,{t,1.0/(n+1),1-(1.0/(n+1)),1.0/(n+1)}]];

ToExpression[DefineDummyFunction];

TestResults = DummyFunction @@InitialPath // Timing;
Print["Testing Dummy Function: got ",TestResults[[2]] ," for initial (straight line) path in ",TestResults[[1]]," seconds"];

DefineNMinCall = "NMinimize[Hold[" <> ToString[DummyFunction] <> "[" <>VarsFlat[n] <> "]],"
<> NMinPattern[n,SearchRectangle] <> ",Method->\"NelderMead\",MaxIterations->5000,AccuracyGoal->4] // Timing" ;

Print[DefineNMinCall];

RawSol = ToExpression[DefineNMinCall];

Print[RawSol];
Print["Completed minimization, computing path length, penalty, timing..."];

{"PathLength" -> ToExpression["PathLength[" <> ToString[f] <> ","<> ToString[gradF]<>"," <> ToString[A] <>"," <> ToString[(ToExpression[Vars[n]] /. RawSol[[2,2]])] <> "," <> ToString[B] <> "]"] ,
"PenalizedValue" -> RawSol[[2,1]],
"Timing" -> RawSol[[1]],
 "Solution"->RawSol[[2,2]]} 

];

LandscapePath[f,gradF,{-2,1},{5,3},{{-10,10},{-10,10}},3]

Solution = {"PathLength"->14.20952917089689`,"PenalizedValue"->14.945435751848088`,"Timing"->109.34066599999999880310497246682643890381`8.059381628154753,"Solution"->{x1->-0.01715657071662758`,y1->3.34728013684471`,x2->1.445799729115478`,y2->3.8314866956311815`,x3->4.3396515067621335`,y3->1.932136442423856`}};

Plotting the solution

Now in order to make any kind of sense of the computation, we need to plot the resulting solution. To do so, we’ll make two plots: one overlaying our track on a contour plot of the function f, and one 3d graph showing the actual curve in space on a transparent plot of the surface.

We’ll need to first interpret the solution as a PathFunction, then plot the PathFunction on top of a ContourPlot for the landscape.

SolutionContourPlot[f_,A_,B_,SearchRect_,Sol_] := Module[{PathCurve,Points,γ,n,Solution},
Solution = "Solution" /. Sol;
n = Length[Solution]/2;
Points = ToExpression[Vars[n]] /. Solution;
PathCurve = PathFunction[A,Points,B];
γ[t_] := {PathCurve[[1]][t],PathCurve[[2]][t]};
ContourPlot[f[x,y],{x,SearchRect[[1,1]],SearchRect[[1,2]]},{y,SearchRect[[2,1]],SearchRect[[2,2]]},Epilog->{List @@ ParametricPlot[γ[t],{t,0,1}]}]
];

SolutionContourPlot[f,{-2,1},{5,3},{{-10,10},{-10,10}},Solution]

The 3d version of the plot involves defining the space curve in question, then plotting it on top of a 3d model of the landscape in question.

Solution3DPlot[f_,A_,B_,SearchRect_,Sol_] := Module[{PathCurve,Points,γ,n,Solution},
Solution ="Solution" /. Sol;
n = Length[Solution]/2;
Points = ToExpression[Vars[n]] /. Solution;
PathCurve = PathFunction[A,Points,B];
γ[t_] := {PathCurve[[1]][t],PathCurve[[2]][t],f[PathCurve[[1]][t],PathCurve[[2]][t]]};

Print[γ[0.2]];

Show[ParametricPlot3D[γ[t],{t,0,1},PlotStyle->Directive[Thick,Black]],Plot3D[f[x,y],{x,SearchRect[[1,1]],SearchRect[[1,2]]},{y,SearchRect[[2,1]],SearchRect[[2,2]]},PlotStyle->Directive[Opacity[0.7]],MeshStyle->None],PlotRange->Join[SearchRect,{{-50,2}}]]
];

Solution3DPlot[f,{-2,1},{5,3},{{-10,10},{-10,10}},Solution]

5. Loading a real landscape.

We now want to try this out with an actual landscape. We’ll use the package from the webpage:

AppendTo[$Path,NotebookDirectory[]]
<<TerrainPackage`

If all has gone well, we now have functions z[x,y] and GradZ[x,y]. Let’s check

z[10.1,12.1]
GradZ[10.1,12.1]

So far, so good! Let’s go and and plot the landscape in various ways to make sure we understand it.

CPlot = ContourPlot[z[x,y],{x,0,25},{y,0,25},Contours->15,PlotPoints->40,ColorFunction->(Lighter[ColorData["LightTerrain"][#+0.2],0.2] &),ImageSize->600,Epilog->{{Disk[{9.57,10.23},0.1],Disk[{20.12,22.3},0.1]},Text[Style["El Chato",FontFamily->"Helvetica",FontSize->14],{9.57,10.75}],
Text[Style["Treasure",FontFamily->"Helvetica",FontSize->14],{20.12,22.9}]}]

myZ[x_,y_] := z[x,y];
myGradZ[{x_,y_}] := GradZ[x,y];

PathLength[myZ,myGradZ,{9.57,10.23},{},{20.12,22.3}]

NMLongSol = LandscapePath[myZ,myGradZ,{9.57,10.23},{20.12,22.3},{{0,25},{0,25}},20]

ReallyLousySol = {54.10947266552972`,{x1->8.11527755494154`,y1->9.596651391609027`,x2->8.323377221850878`,y2->5.0204277023393304`,x3->16.341373807907253`,y3->12.377109123515295`,x4->17.115333513023568`,y4->10.587183479823805`,x5->14.107917943668621`,y5->11.89948357424181`,x6->13.700347223933298`,y6->10.974107183466026`,x7->13.674047690021194`,y7->13.39997423277222`,x8->15.923290840401576`,y8->16.685704600724502`,x9->11.442454431368525`,y9->20.354139903343654`,x10->12.006359615740898`,y10->23.795490560349155`,x11->15.698400638779773`,y11->21.537085054427013`,x12->20.724735005110226`,y12->20.645739641194236`}};

LongSASol = {16.854082885761663`,{x1->9.625828240652673`,y1->10.417473454646672`,x2->9.667063084750122`,y2->10.680450706756071`,x3->9.665177048582855`,y3->10.764685364808717`,x4->9.687829607732585`,y4->10.868328430854103`,x5->9.692596564524544`,y5->10.945225291627048`,x6->9.703769671400744`,y6->10.991780855203048`,x7->9.747432263246976`,y7->11.140536084282722`,x8->9.759856789130515`,y8->11.177296235811019`,x9->9.777984639104947`,y9->11.249872601292891`,x10->9.768186783219816`,y10->11.300182897341134`,x11->9.76358219772571`,y11->11.38925841298366`,x12->9.770667609379874`,y12->11.485810644643914`,x13->9.796586271207214`,y13->11.656052972679412`,x14->9.845236057396262`,y14->11.767541320852747`,x15->9.877164074695573`,y15->11.915726557828027`,x16->9.92931582121598`,y16->12.106337851821422`,x17->9.98225628705648`,y17->12.23368886542739`,x18->10.088952470189902`,y18->12.372233971371333`,x19->10.181722731356572`,y19->12.504887220137656`,x20->10.260688137773629`,y20->12.63732266140134`,x21->10.356118495697512`,y21->12.791494740073261`,x22->10.368388170122246`,y22->12.847829001159011`,x23->10.423032430429261`,y23->12.88623207219203`,x24->10.437040602767299`,y24->12.97266243390815`,x25->10.470485182250371`,y25->13.067118201829155`,x26->10.603144433357405`,y26->13.215513969421806`,x27->10.643409269186725`,y27->13.248433759274377`,x28->10.643032331731357`,y28->13.254266274865744`,x29->10.661968985369706`,y29->13.299724098984285`,x30->10.84426805465367`,y30->13.453793535237974`,x31->10.89752610256596`,y31->13.517868824156693`,x32->10.96480660325899`,y32->13.57900931491072`,x33->10.980171998349878`,y33->13.57376985857335`,x34->10.986505662299606`,y34->13.584778703528444`,x35->11.006424796926826`,y35->13.609831729096603`,x36->11.015418440923016`,y36->13.62910768546854`,x37->11.036978339064667`,y37->13.67693823949964`,x38->11.117090803265619`,y38->13.71600953385222`,x39->11.125107861940434`,y39->13.72319797272263`,x40->11.255258973712634`,y40->13.81908869284822`,x41->11.28956242560226`,y41->13.871911887543122`,x42->11.34672173257365`,y42->13.914719022759654`,x43->11.401987306439363`,y43->13.94064499965686`,x44->11.46296727182298`,y44->13.990780644156827`,x45->11.550473751654156`,y45->14.055729508975073`,x46->11.668990855736292`,y46->14.16137062798129`,x47->11.874946076714869`,y47->14.355880626652622`,x48->12.595800643140546`,y48->14.991755975929284`,x49->14.906560710383218`,y49->17.140955885813863`,x50->16.847441830202737`,y50->19.076838800686385`}};

RottenSol = {17.689817085806766`,{x1->9.840277043539475`,y1->10.216891912728917`,x2->10.014838174369352`,y2->10.194186118359548`,x3->10.098010964458478`,y3->10.13960161936455`,x4->10.190705193708204`,y4->10.125338969463655`,x5->10.32169018089812`,y5->10.136187391857751`,x6->10.469979979957`,y6->10.126072911325823`,x7->10.53631713962968`,y7->10.131574009116365`,x8->10.65809791557109`,y8->10.19909097711026`,x9->10.816575172136814`,y9->10.257582091673095`,x10->10.892937815105867`,y10->10.282720983673853`,x11->10.906593136602051`,y11->10.28470386303373`,x12->10.96418632632281`,y12->10.337534560812543`,x13->11.035496175954052`,y13->10.369106816496089`,x14->11.360167915450896`,y14->10.496239846655655`,x15->11.442990888836597`,y15->10.623115157810112`,x16->11.741533949343456`,y16->10.989296686432125`,x17->11.79252763017392`,y17->11.055666056901062`,x18->11.926794333474705`,y18->11.3203931855274`,x19->11.981773126792516`,y19->11.404710681962921`,x20->11.954686784671363`,y20->11.590597037972`,x21->11.940352124374506`,y21->11.624355935947916`,x22->11.941874747236707`,y22->11.64811130707023`,x23->11.925241053473991`,y23->11.718078790979579`,x24->11.910229344033445`,y24->11.744796751953917`,x25->11.901373318042365`,y25->11.748567613294183`,x26->11.89212525171909`,y26->11.779640734945662`,x27->11.864473631760951`,y27->11.839869072425829`,x28->11.84501217504143`,y28->11.89167572529095`,x29->11.785064957349816`,y29->12.00405389193363`,x30->11.77201744630628`,y30->12.008616577584723`,x31->11.749757505785158`,y31->12.039757188404321`,x32->11.743490782372092`,y32->12.05054298910826`,x33->11.74157753285987`,y33->12.054793714128111`,x34->11.746504474818769`,y34->12.061665847070337`,x35->11.739393099500576`,y35->12.072929684298433`,x36->11.73873844616448`,y36->12.092076589861598`,x37->11.738686475612345`,y37->12.09447894190077`,x38->11.755054944632697`,y38->12.110488427240817`,x39->11.774036915299025`,y39->12.134648866796516`,x40->11.786820141521504`,y40->12.149095240365455`,x41->11.840730598035698`,y41->12.204141623229185`,x42->11.925785145203301`,y42->12.397956346653517`,x43->12.082879467276575`,y43->12.748855040473385`,x44->12.276345352050114`,y44->12.907686890308867`,x45->12.3307656928498`,y45->12.937349501883162`,x46->12.92822709446221`,y46->13.266844002017866`,x47->12.95305606915957`,y47->13.335883606886595`,x48->13.40639673541439`,y48->14.130116220994772`,x49->14.790554299624178`,y49->16.085721472293635`,x50->16.271835658153616`,y50->17.1900633322529`}}

Sol // InputForm

Sol = {16.340466910627192`,{x1->10.620503826627484`,y1->11.313539262092313`,x2->11.286361159506763`,y2->11.992152414710015`,x3->13.007758965199603`,y3->14.021176927561392`,x4->15.281192401394764`,y4->16.609896999439623`,x5->16.430809761446493`,y5->17.831447370570675`}};

SolutionLandscapePlot[A_,B_,SearchRect_,Sol_] := Module[{PathCurve,Points,γ,n,Solution},
Solution = "Solution" /. Sol;
n = Length[Solution]/2;
Points = ToExpression[Vars[n]] /. Solution;
PathCurve = PathFunction[A,Points,B];
γ[t_] := {PathCurve[[1]][t],PathCurve[[2]][t]};
Show[{CPlot,ParametricPlot[γ[t],{t,0,1}],Graphics[Disk[#,0.1]] & /@ Points}]
];

SolutionLandscapePlot[{9.57,10.23},{20.12,22.3},{{0,25},{0,25}},InterpOrder1Sol]

LandscapeContourPlot[Sol_] := Module[{CPoints,FPlot,BCurve},
CPoints = Join[{{-115,35}},{{x1,y1},{x2,y2},{x3,y3},{x4,y4}},{{-119,32}}] /. Sol[[2,2]];
BCurve = BezierFunction[CPoints];
ContourPlot[z[x,y],{x,-120,-110},{y,30,40},ColorFunction->ColorData["LightTerrain"],Epilog->{List @@ ParametricPlot[BCurve[t],{t,0,1}]}]
];

LandscapeContourPlot[Solution]

LandscapeSolution3DPlot[Sol_] := Module[{CPoints,FPlot,BCurve,SpaceCurve},
CPoints = Join[{{-115,35}},{{x1,y1},{x2,y2},{x3,y3},{x4,y4}},{{-119,32}}] /. Sol[[2,2]];
BCurve = BezierFunction[CPoints];
SpaceCurve[t_,BCurve_] := Module[{x,y}, {x,y} = BCurve[t]; {x,y,z[x,y]}];
Show[ParametricPlot3D[SpaceCurve[t,BCurve],{t,0,1},PlotStyle->Directive[Thick,Black]],Plot3D[z[x,y],{x,-120,-110},{y,30,40},PlotStyle->Directive[Opacity[0.9]],MeshStyle->None,ColorFunction->ColorData["LightTerrain"]],PlotRange->{{-120,-110},{30,40},{-3000,4000}},BoxRatios->{1,1,0.1}]
];

LandscapeSolution3DPlot[Solution]