Skip to content

Commit

Permalink
Lot of work
Browse files Browse the repository at this point in the history 
- Working BS3 implementation
- Start of event support
- Additional DP5 random testing
- Fixes major interpolant caching/re-use bug
- Adds Bisection method (used by eventing)
- Eventing may or may not work at all right now

Signed-off-by: Lamont Granquist <lamont@scriptkiddie.org>
  • Loading branch information
@lamont-granquist
lamont-granquist committed May 9, 2023
1 parent d277897 commit 51e44f9
Show file tree
Hide file tree
Showing 21 changed files with 1,001 additions and 500 deletions.
1 change: 1 addition & 0 deletions MechJeb2.sln.DotSettings
View file
Expand Up @@ -3,6 +3,7 @@
<s:String x:Key="/Default/CodeStyle/Naming/CSharpNaming/Abbreviations/=AOA/@EntryIndexedValue">AOA</s:String>
<s:String x:Key="/Default/CodeStyle/Naming/CSharpNaming/Abbreviations/=ASS/@EntryIndexedValue">ASS</s:String>
<s:String x:Key="/Default/CodeStyle/Naming/CSharpNaming/Abbreviations/=BCI/@EntryIndexedValue">BCI</s:String>
<s:String x:Key="/Default/CodeStyle/Naming/CSharpNaming/Abbreviations/=BS/@EntryIndexedValue">BS</s:String>
<s:String x:Key="/Default/CodeStyle/Naming/CSharpNaming/Abbreviations/=DP/@EntryIndexedValue">DP</s:String>
<s:String x:Key="/Default/CodeStyle/Naming/CSharpNaming/Abbreviations/=ECI/@EntryIndexedValue">ECI</s:String>
<s:String x:Key="/Default/CodeStyle/Naming/CSharpNaming/Abbreviations/=ENU/@EntryIndexedValue">ENU</s:String>
Expand Down
510 changes: 256 additions & 254 deletions MechJeb2/MechJeb2.csproj
View file

Large diffs are not rendered by default.

83 changes: 83 additions & 0 deletions MechJeb2/MechJebLib/Core/Bisection.cs
View file
@@ -0,0 +1,83 @@
/*
* Copyright Lamont Granquist, Sebastien Gaggini and the MechJeb contributors
* SPDX-License-Identifier: MIT-0 OR LGPL-2.1+ OR CC0-1.0
*/

#nullable enable

using System;
using MechJebLib.Utils;

// ReSharper disable CompareOfFloatsByEqualityOperator
namespace MechJebLib.Core
{
using RootFunc = Func<double, object?, double>;

public class Bisection
{
/// <summary>
/// Bisection search.
/// </summary>
/// <param name="f">1 dimensional function</param>
/// <param name="a">minimum bounds</param>
/// <param name="b">maximum bounds</param>
/// <param name="o">state object to pass to function</param>
/// <param name="rtol">tolerance</param>
/// <param name="preferLeft"></param>
/// <returns>value for which the function evaluates to zero</returns>
/// <exception cref="ArgumentException">guess does not brack the root</exception>
public static (double c, double fc) Solve(RootFunc f, double a, double b, object? o, double tolX = 0,
bool preferLeft = true)
{
Check.Finite(a);
Check.Finite(b);
Check.NonNegative(tolX);

if (a > b)
{
(a, b) = (b, a);
preferLeft = !preferLeft;
}

double fa = f(a, o);
double fb = f(b, o);

Check.Finite(fa);
Check.Finite(fb);

if (fa == 0)
return (a, fa);

if (fb == 0)
return (b, fb);

if (fa * fb > 0)
throw new ArgumentException("Bisection rootfinding method: guess does not bracket the root");

while (true)
{
double c = (a + b) / 2;

if (c == a || c == b || b - c < tolX)
break;

double fc = f(c, o);

if (fc * fa < 0)
{
b = c;
fb = fc;
}
else
{
a = c;
fa = fc;
}
}

if (preferLeft)
return (a, fa);
return (b, fb);
}
}
}
3 changes: 1 addition & 2 deletions MechJeb2/MechJebLib/Core/BrentMin.cs
View file
Expand Up @@ -5,14 +5,13 @@
*/

using System;
using System.Diagnostics.CodeAnalysis;
using static MechJebLib.Utils.Statics;

#nullable enable

// ReSharper disable CompareOfFloatsByEqualityOperator
namespace MechJebLib.Core
{
[SuppressMessage("ReSharper", "CompareOfFloatsByEqualityOperator")]
public static class BrentMin
{
// Brent's 1-dimensional derivative-free local minimization method
Expand Down
15 changes: 7 additions & 8 deletions MechJeb2/MechJebLib/Core/BrentRoot.cs
View file
@@ -1,23 +1,22 @@
/*
* Copyright Lamont Granquist (lamont@scriptkiddie.org)
* Dual licensed under the MIT (MIT-LICENSE) license
* and GPLv2 (GPLv2-LICENSE) license or any later version.
* Copyright Lamont Granquist, Sebastien Gaggini and the MechJeb contributors
* SPDX-License-Identifier: MIT-0 OR LGPL-2.1+ OR CC0-1.0
*/

#nullable enable

using System;
using MechJebLib.Utils;
using static MechJebLib.Utils.Statics;

// ReSharper disable CompareOfFloatsByEqualityOperator
namespace MechJebLib.Core
{
/// <summary>
/// Brent's rootfinding method.
/// </summary>
public static class BrentRoot
{
private const double EPS = 2.24e-15;

/// <summary>
/// Brent's rootfinding method, bounded search. Uses secant, inverse quadratic interpolation and bisection.
/// </summary>
Expand Down Expand Up @@ -56,7 +55,7 @@ public static class BrentRoot
(fa, fb) = (fb, fa);
}

return InternalSolve(f, a, b, fa, fb, o, maxiter, rtol, sign);
return _Solve(f, a, b, fa, fb, o, maxiter, rtol, sign);
}

/// <summary>
Expand Down Expand Up @@ -118,11 +117,11 @@ public static double Solve(Func<double, object?, double> f, double x, object? o,
return b;
}

return InternalSolve(f, a, b, fa, fb, o, maxiter, rtol, sign);
return _Solve(f, a, b, fa, fb, o, maxiter, rtol, sign);
}

// This is the actual algorithm
private static double InternalSolve(Func<double, object?, double> f, double a, double b, double fa, double fb, object? o, int maxiter,
private static double _Solve(Func<double, object?, double> f, double a, double b, double fa, double fb, object? o, int maxiter,
double rtol, int sign)
{
bool maybeOneMore = sign != 0;
Expand Down
153 changes: 111 additions & 42 deletions MechJeb2/MechJebLib/Core/ODE/AbstractIVP.cs
View file
Expand Up @@ -14,9 +14,13 @@
// ReSharper disable CompareOfFloatsByEqualityOperator
namespace MechJebLib.Core.ODE
{
using IVPFunc = Action<Vn, double, Vn>;
using IVPEvent = Func<double, Vn, Vn, (double x, bool dir, bool stop)>;
using IVPFunc = Action<IList<double>, double, IList<double>>;

// TODO:
// - Needs better MinStep based on next floating point number
// - Configurable to throw or just continue at MinStep
// - Needs better initial step guessing
// - Needs working event API
public abstract class AbstractIVP
{
/// <summary>
Expand Down Expand Up @@ -75,44 +79,51 @@ public abstract class AbstractIVP
/// <param name="events"></param>
/// <exception cref="ArgumentException"></exception>
public void Solve(IVPFunc f, IReadOnlyList<double> y0, IList<double> yf, double t0, double tf, Hn? interpolant = null,
IReadOnlyCollection<IVPEvent>? events = null)
IReadOnlyList<Event>? events = null)
{
try
{
N = y0.Count;
Ynew = Vn.Rent(N);
Dynew = Vn.Rent(N);
Dy = Vn.Rent(N);
Y = Vn.Rent(N);
Y = Y.Expand(N);
Dy = Dy.Expand(N);
Ynew = Ynew.Expand(N);
Dynew = Dynew.Expand(N);

Init();
_Solve(f, y0, yf, t0, tf, interpolant, events);
}

finally
{
Y.Dispose();
Dy.Dispose();
Ynew.Dispose();
Dynew.Dispose();
Y = Ynew = Dy = Dynew = null!;

Cleanup();
}
}

protected Vn Y = null!;
protected Vn Ynew = null!;
protected Vn Dy = null!;
protected Vn Dynew = null!;
protected double Habs;
protected int Direction;
protected double T, Tnew;
protected double MaxStep;
protected double MinStep;
private double _habsNext;
protected double[] Y = new double[1];
protected double[] Ynew = new double[1];
protected double[] Dy = new double[1];
protected double[] Dynew = new double[1];
protected double Habs;
protected int Direction;
protected double T, Tnew;
protected double MaxStep;
protected double MinStep;
private double _habsNext;
private readonly List<Event> _activeEvents = new List<Event>();

private Func<double, IList<double>, AbstractIVP, double> _eventFunc = null!;

private Func<double, object?, double> _eventFunctionDelegate => EventFuncWrapper;

private double EventFuncWrapper(double x, object? o)
{
using var yinterp = Vn.Rent(N);
Interpolate(x, yinterp);
return _eventFunc(x, yinterp, this);
}

private void _Solve(IVPFunc f, IReadOnlyList<double> y0, IList<double> yf, double t0, double tf, Hn? interpolant,
IReadOnlyCollection<IVPEvent>? events)
IReadOnlyList<Event>? events)
{
Direction = t0 != tf ? Math.Sign(tf - t0) : 1;
Habs = SelectInitialStep(t0, tf);
Expand All @@ -128,7 +139,15 @@ public abstract class AbstractIVP

interpolant?.Add(T, Y, Dy);

while (T != tf)
if (events != null)
{
for (int i = 0; i < events.Count; i++)
events[i].LastValue = events[i].F(T, Y, this);
}

bool terminate = false;

while ((Direction > 0 && T < tf) || (Direction < 0 && T > tf))
{
CancellationToken.ThrowIfCancellationRequested();

Expand All @@ -143,37 +162,57 @@ public abstract class AbstractIVP

Tnew = T + Habs * Direction;

// handle events, this assumes only one trigger per event per step
if (events != null)
{
}

// extract a low fidelity interpolant
if (interpolant != null)
{
while (interpCount < Interpnum)
for (int i = 0; i < events.Count; i++)
{
double tinterp = t0 + (tf - t0) * interpCount / Interpnum;

if (!tinterp.IsWithin(T, Tnew))
break;

using var yinterp = Vn.Rent(N);
using var finterp = Vn.Rent(N);
events[i].NewValue = events[i].F(T, Y, this);
_activeEvents.Clear();
if (IsActiveEvent(events[i]))
_activeEvents.Add(events[i]);
}

if (_activeEvents.Count > 0)
{
InitInterpolant();
Interpolate(tinterp, yinterp);
f(yinterp, tinterp, finterp);
interpolant?.Add(tinterp, yinterp, finterp);
interpCount++;

for (int i = 0; i < _activeEvents.Count; i++)
{
_eventFunc = _activeEvents[i].F;
(double tevent, _) = Bisection.Solve(_eventFunctionDelegate, T, Tnew, null, EPS);
_activeEvents[i].Time = tevent;
}

_activeEvents.Sort();

for (int i = 0; i < _activeEvents.Count; i++)
{
if (_activeEvents[i].Terminal)
{
terminate = true;
break;
}
}
}

for (int i = 0; i < events.Count; i++)
events[i].LastValue = events[i].NewValue;
}

// extract a low fidelity interpolant
if (interpolant != null)
interpCount = FillInterpolant(f, t0, tf, interpolant, interpCount);

// take a step
Y.CopyFrom(Ynew);
Dy.CopyFrom(Dynew);
T = Tnew;
Habs = _habsNext;

if (terminate)
break;

// handle max iterations
if (Maxiter > 0 && niter++ > Maxiter)
{
Expand All @@ -189,6 +228,36 @@ public abstract class AbstractIVP
Y.CopyTo(yf);
}

private bool IsActiveEvent(Event e)
{
bool up = e.LastValue <= 0 && e.NewValue >= 0;
bool down = e.LastValue >= 0 && e.NewValue <= 0;
bool either = up || down;
return (up && e.Direction > 0) || (down && e.Direction < 0) || (either && e.Direction == 0);
}

private int FillInterpolant(IVPFunc f, double t0, double tf, Hn interpolant, int interpCount)
{
while (interpCount < Interpnum)
{
double tinterp = t0 + (tf - t0) * interpCount / Interpnum;

if (!tinterp.IsWithin(T, Tnew))
break;

using var yinterp = Vn.Rent(N);
using var finterp = Vn.Rent(N);

InitInterpolant();
Interpolate(tinterp, yinterp);
f(yinterp, tinterp, finterp);
interpolant?.Add(tinterp, yinterp, finterp);
interpCount++;
}

return interpCount;
}

protected abstract (double, double) Step(IVPFunc f);
protected abstract double SelectInitialStep(double t0, double tf);
protected abstract void InitInterpolant();
Expand Down
3 changes: 1 addition & 2 deletions MechJeb2/MechJebLib/Core/ODE/AbstractRungeKutta.cs
View file
Expand Up @@ -13,8 +13,7 @@

namespace MechJebLib.Core.ODE
{
using IVPFunc = Action<Vn, double, Vn>;
using IVPEvent = Func<double, Vn, Vn, (double x, bool dir, bool stop)>;
using IVPFunc = Action<IList<double>, double, IList<double>>;

public abstract class AbstractRungeKutta : AbstractIVP
{
Expand Down

0 comments on commit 51e44f9

Please sign in to comment.