Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

core: compute min between 2 envelope parts and use it in ETCS #10991

Open
wants to merge 2 commits into
base: alch/core/add_svl_etcs_logic
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions core/envelope-sim/build.gradle
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,9 @@ dependencies {
// Use JUnit Jupiter Engine for testing.
testRuntimeOnly libs.junit.jupiter.engine

// Use AssertJ for testing
testImplementation libs.assertj

// for linter annotations
testFixturesCompileOnly libs.jcip.annotations
testFixturesCompileOnly libs.spotbugs.annotations
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -9,10 +9,7 @@
import fr.sncf.osrd.envelope_sim.EnvelopeProfile;
import fr.sncf.osrd.envelope_utils.ExcludeFromGeneratedCodeCoverage;
import fr.sncf.osrd.utils.SelfTypeHolder;
import java.util.Arrays;
import java.util.Collections;
import java.util.HashMap;
import java.util.Map;
import java.util.*;
import java.util.stream.Collectors;

/**
Expand Down Expand Up @@ -616,6 +613,65 @@ public EnvelopePart copyAndShift(double positionDelta, double minPosition, doubl
new HashMap<>(attrs), newPositions.toArray(), newSpeeds.toArray(), newTimeDeltas.toArray());
}

/**
* Compute the minimum EnvelopePart between 2 envelope parts with an intersecting range [a, b].
* The resulting envelope part starts at the envelope with the minimum speed at a and ends at b.
*/
public static EnvelopePart min(
EnvelopePart envelopePartA, EnvelopePart envelopePartB, Iterable<SelfTypeHolder> attrs) {
Comment on lines +616 to +621
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would have thought it would have been easier to use if it returned the whole range covered by either. And possibly dealing with Envelopes instead of Parts. But you're the caller, you probably know what's most convenient.

Copy link
Contributor Author

@Erashin Erashin Mar 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here were my thoughts on this, tell me what you think about it:

  • if by "the whole range covered by either" (probably not what you mean) you mean the whole range, even after the end of one, well the envelope part would have discontinuities, so this isn't doable.
  • if by "the whole range covered by either" you mean just the intersecting range of the two curves, then yeah that is a possibility. That would mean moving the part of the code that creates the first non intersecting part of the curve into EtcsBrakingCurves, which is clearly doable.
  • for envelopes instead of envelope parts: the pb is how we're going to manage the discontinuities between the successive envelope parts, because they're not necessarily continuous (in our case they are). Except if you have a sol for that, it would work if we assert it is continuous for a min method between 2 Envelopes. The changes would be: putting the current method into a utils, taking positions + speeds as arguments. Call it in both EnvelopePart.min and Envelope.min (2 birds with one stone). EDIT: even then, how would we be managing the different EnvelopeParts' attributes?

Copy link
Contributor

@eckter eckter Mar 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can and do have discontinuities in envelopes, it's just in the final simulation output that they're forbidden.

I meant the whole range, or possibly the whole range given by one of the inputs.

I'm just a little surprised to see a method where we don't really know what range will be covered by the output. Generally speaking we have a reference envelope and we modify parts of it, so having the solution defined on the intersection of two parts feels a little off.

It would probably be easier to understand if you describe where it would be called: what are the two parts that will be used as input? And why are we only interested in the intersection? (and why can't they be envelopes?)

Copy link
Contributor Author

@Erashin Erashin Mar 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What I meant by Envelopes with discontinuities not being possible to use the method on is that we can't use interpolateSpeed and intersectStep on them. Imagine envelopeA going from a constant speed of 100 to a constant speed of 50, and envelopeB being at 75 at the time: how do we manage the intersection between the 2 curves? Would you just cut it off, create a new EnvelopePart starting at this discontinuity? How do we manage this without the method having quite a few specific cases and ending up being quite complex?

The 2 parts used as input are just 2 random parts in my head, except they have to have a common range of position because otherwise how would we determine which of them is the min? Only alternative i see would be to return an Envelope combining them both in that case, but that would mean the method shouldn't be called min, but probably something like mergeMin or something. And it would also mean the output can be discontinuous, which leads to the following point.

Why are we only interested in the intersection: the objective of this method, and also the reason why I'm working with EnvelopeParts rather than Envelopes, is to keep the continuity property of the EnvelopePart: the output will be continuous (which is handy when creating braking curves). The only way to have the output be continuous is:

  • taking the min of the 2 curves at least in the position-range intersection
  • following the curve before the intersection if it is the min at the start of the intersection
  • following the curve after the intersection if it is the min at the end of the intersection (which i could have coded actually)

However, as you've said, a good alternative would be to have the position-range given as an input: only pb would be that, if the output is still an EnvelopePart, you'd have cases when the output's range would be smaller than the input range for continuity reasons. Which i think is ok-ish kinda.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So I feel like we're kind of dodging some important questions by only returning the intersection. The questions being what we're going to do at those borders. I wouldn't want to have the caller of this method reassemble some of the input parts to build a consistent final result, this method wouldn't be usable on its own.

What I'm really uncomfortable about is not really knowing the "domain" of the output. We expect a set of points but without knowing where they'll be, it feels like we'd need a lot of weird post-processing (and pre-processing to even have a valid overlap).

I'd really like if we could:

  1. have a clear output domain that can't be ill-defined (such as "all of the range covered by the leftmost parameter")
  2. Handle discontinuities in this method, or possibly in an auxiliary one. I don't think handling envelope parts really helps here.

Actually we've had very similar use cases with mareco. Generally we tolerate discontinuities during intermediate steps, and then add accelerations or decelerations to fix them. But we need to plan carefully to avoid having weird errors or ill-defined results, the mareco module took a lot of rewrites (and isn't even stable today).


I know I'm not really helping much here, I don't have any clear path to offer. If possible I'd like some fresh eyes (@bougue-pe 👀). Just to be clear: I'm not strongly against the current solution, just a little afraid of what's coming next

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

😅 This is a bit tough, I hope I understood the stakes.

The first use-case:

  • computing the min of 2 braking curves that are continuous, monotonous (maybe not strictly if we consider the maintain-speed part of SvL).
  • the result is expected to start from the [min(beginPositions); maxSpeed] point (both curves start -or can start- at maxSpeed).
  • the end of the result is debatable:
    • can finish at min(endPositions), then external work to add final parts: maintain-speed to reach second curve, then use the second curve to its endPosition.
    • can finish at max(endPositions): internally add a maintain-speed part to link to second curve when min(endPositions) is reached.

I see 4 valid options here (detailed by my order of preference).

Important

TLDR

  • Working on Envelopes (and assume discontinuities) has my preference (on the expected RoI also).
  • We can easily manage continuity and size of position-ranges for the described case by adding missing parts, so any solution looks valid.

1 - work on Envelopes (force one intersection, return union of position ranges)

That would mean assuming discontinuities sometimes, indeed.
But forcing monotonous (and in the same direction) inputs guarantees monotonous output on the intersecting range.
And forcing continuous inputs guarantees continuous output on the intersecting range.
Monotony and/or continuity are quite easy to obtain on the first use-case we have (adding constant parts at the beginning and the end of curves so that their position range is equal), so it wouldn't be a major issue.

That would also mean applying attributes given in inputs, which is better IMO (allows for finer result).

That would require a clear choice of what attributes are applied when strictly equal (the previous min curve, then the next min curve, then the leftmost param if no previous?).

That means that going back to working on EnvelopePart is a bit tricky if required (converting Envelope slice to a Part while checking that Part requirements are respected).

That may require some kind of "final cleanup" (optional?) to regroup parts that are "identical" to obtain a clean single braking curve when possible (in our use case, it is not mandatory to avoid mutliple consecutive braking parts in a single ETCS braking, is it?).

That may require a bit more work to implement and to fit it in current code.

From the point-of-view of an external user of this method, I think we would have more use of this: it can more safely be used to add braking and acceleration curves (the extra cost on perf may be an issue?). Superposing them blindly would just work, and the code would be much clearer IMO.

2 - work on EnvelopeParts and return a union of position ranges (force one intersection, identical monotony)

TBD if the result should be an Envelope (my preference) or an EnvelopePart, see below.

The logic of the method would actually be a bit "reversed": return a curve that's the min position for all the speed range, and keep the min speed when 2 positions lead to different speed.

That would also mean assuming discontinuities sometimes (but it would not happen for the first use-case).
That would mean adding a maintain-speed piece:

  • either a point to the EnvelopePart result if allowed (post-processing to change it to a separate EnvelopePart?)
  • or just a dedicated EnvelopePart to the Envelope result (better IMO)

The amount of work is OK I think to have it fit in the current code.

Less external "unclear" code to assemble parts, the logic is contained and a bit more easy to explain IMO.

3 - Keep the current code

Staying with the current proposition is OK IMO, but we should probably move it back to the ETCSBrakingCurve.kt to avoid any misuse, as I agree it's a bit of a trap (for 2 acceleration curves, the top part would be broken).

Or keep it here, add asserts on braking and TODO for the rest to ensure future mutualization.

4 - work on EnvelopeParts and return intersection of position ranges (force one intersection, continuity, identical monotony)

This guarantees a continuous and monotonous output.
We may also force identical attributes of inputs to use the same one in output?
In the first use-case, we could add an extra first-point (pos: 0, speed: beginSpeed) to make sure that intersection is enough, I think.

This would work for our first case, cover the concern on the output position-range and the amount of work is OK I guess.

It would keep the "unclear" external code that's assembling the parts.

var beginPosA = envelopePartA.getBeginPos();
var beginPosB = envelopePartB.getBeginPos();
var endPosA = envelopePartA.getEndPos();
var endPosB = envelopePartB.getEndPos();
assert (beginPosA < endPosB && endPosA > beginPosB);
var startIntersectingRange = Math.max(beginPosA, beginPosB);
var start = envelopePartA.interpolateSpeed(startIntersectingRange)
<= envelopePartB.interpolateSpeed(startIntersectingRange)
? beginPosA
: beginPosB;
Comment on lines +628 to +631
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(if keeping the current behaviour) If speeds are equal at the beginning of pos-intersection, we probably want the min(beginPosA, beginPosB) (or max).
Otherwise, the result won't be identical if we exchange envelopePartA and envelopePartB.

var end = Math.min(endPosA, endPosB);

TreeSet<Double> keyPositions =
Arrays.stream(envelopePartA.positions).boxed().collect(Collectors.toCollection(TreeSet::new));
Comment on lines +634 to +635
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I feel like we should really convert that file to kotlin if we write that kind of stuff 😄

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was actually in the midst of doing it and there were a few pbs while i was converting it. I'll give it another shot though, we both agree it's better.

keyPositions.addAll(Arrays.stream(envelopePartB.positions).boxed().collect(Collectors.toSet()));
keyPositions = new TreeSet<>(keyPositions.subSet(start, true, end, true));
Comment on lines +634 to +637
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
TreeSet<Double> keyPositions =
Arrays.stream(envelopePartA.positions).boxed().collect(Collectors.toCollection(TreeSet::new));
keyPositions.addAll(Arrays.stream(envelopePartB.positions).boxed().collect(Collectors.toSet()));
keyPositions = new TreeSet<>(keyPositions.subSet(start, true, end, true));
NavigableSet<Double> keyPositions =
Arrays.stream(envelopePartA.positions).boxed().collect(Collectors.toCollection(TreeSet::new));
keyPositions.addAll(Arrays.stream(envelopePartB.positions).boxed().collect(Collectors.toSet()));
keyPositions = keyPositions.subSet(start, true, end, true);

Avoids instantiating a new TreeSet from the subset

var keyPosList = new ArrayList<>(keyPositions);

var newPositions = new DoubleArrayList();
var newSpeeds = new DoubleArrayList();
for (int i = 0; i < keyPosList.size(); i++) {
var pos = keyPosList.get(i);
boolean inEnvelopePartA = pos >= beginPosA;
boolean inEnvelopePartB = pos >= beginPosB;

double speedA = inEnvelopePartA ? envelopePartA.interpolateSpeed(pos) : Double.POSITIVE_INFINITY;
double speedB = inEnvelopePartB ? envelopePartB.interpolateSpeed(pos) : Double.POSITIVE_INFINITY;
double minSpeedAtPos = Math.min(speedA, speedB);

if (i > 0) {
double prevPos = keyPosList.get(i - 1);
boolean prevInEnvelopePartA = prevPos >= beginPosA;
boolean prevInEnvelopePartB = prevPos >= beginPosB;

if (inEnvelopePartA && inEnvelopePartB && prevInEnvelopePartA && prevInEnvelopePartB) {
double prevSpeedA = envelopePartA.interpolateSpeed(prevPos);
double prevSpeedB = envelopePartB.interpolateSpeed(prevPos);

if ((prevSpeedA - prevSpeedB) * (speedA - speedB) < 0) {
var intersection = EnvelopePhysics.intersectSteps(
prevPos, prevSpeedA, pos, speedA, prevPos, prevSpeedB, pos, speedB);
// Add intersection point
newPositions.add(intersection.position);
newSpeeds.add(intersection.speed);
}
}
}
newPositions.add(pos);
newSpeeds.add(minSpeedAtPos);
}
return generateTimes(attrs, newPositions.toArray(), newSpeeds.toArray());
}

// endregion

// region EQUALS
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -126,45 +126,11 @@ private fun computeMinSvlEoaIndCurve(
releaseSpeedPositionSvl,
NATIONAL_RELEASE_SPEED
)
val startIntersectingRange =
max(slicedIndicationCurveEoa.beginPos, slicedIndicationCurveSvl.beginPos)
var refCurve = slicedIndicationCurveEoa
var otherCurve = slicedIndicationCurveSvl
if (
slicedIndicationCurveEoa.interpolateSpeed(startIntersectingRange) >
slicedIndicationCurveSvl.interpolateSpeed(startIntersectingRange)
) {
refCurve = slicedIndicationCurveSvl
otherCurve = slicedIndicationCurveEoa
}
val pointCount = refCurve.pointCount()
var indicationPositions = DoubleArray(pointCount)
var indicationSpeeds = DoubleArray(pointCount)
for (i in 0 until pointCount) {
val newPos = refCurve.getPointPos(i)
val newSpeed = refCurve.getPointSpeed(i)
if (newPos < otherCurve.beginPos) {
indicationPositions[i] = newPos
indicationSpeeds[i] = newSpeed
} else if (newPos <= otherCurve.endPos) {
val otherSpeed = slicedIndicationCurveSvl.interpolateSpeed(newPos)
indicationPositions[i] = newPos
// TODO: unneeded for now: interpolate to not approximate position at intersection.
indicationSpeeds[i] = min(otherSpeed, newSpeed)
} else {
indicationPositions[i] = otherCurve.endPos
indicationSpeeds[i] = otherCurve.endSpeed
// Clean up the last unneeded points in the arrays before exiting the loop.
indicationPositions = indicationPositions.dropLast(pointCount - 1 - i).toDoubleArray()
indicationSpeeds = indicationSpeeds.dropLast(pointCount - 1 - i).toDoubleArray()
break
}
}
val firstIndicationPart =
EnvelopePart.generateTimes(
listOf(EnvelopeProfile.BRAKING),
indicationPositions,
indicationSpeeds
EnvelopePart.min(
slicedIndicationCurveEoa,
slicedIndicationCurveSvl,
listOf(EnvelopeProfile.BRAKING)
)
if (releaseSpeedPositionSvl < releaseSpeedPositionEoa) {
val maintainReleaseSpeedPart =
Expand Down
Original file line number Diff line number Diff line change
@@ -1,14 +1,19 @@
package fr.sncf.osrd.envelope;

import static fr.sncf.osrd.envelope.EnvelopePhysics.getPartMechanicalEnergyConsumed;
import static org.assertj.core.api.Assertions.assertThat;
import static org.junit.jupiter.api.Assertions.*;

import fr.sncf.osrd.envelope.part.EnvelopePart;
import fr.sncf.osrd.envelope_sim.*;
import fr.sncf.osrd.envelope_sim.allowances.utils.AllowanceValue;
import java.util.Collections;
import java.util.List;
import org.junit.jupiter.api.Assertions;
import java.util.stream.Stream;
import org.junit.jupiter.api.Test;
import org.junit.jupiter.params.ParameterizedTest;
import org.junit.jupiter.params.provider.Arguments;
import org.junit.jupiter.params.provider.MethodSource;

class EnvelopePartTest {
@Test
Expand Down Expand Up @@ -93,7 +98,7 @@ void testGetMechanicalEnergyConsumed() {
/ 1_000_000;
break;
case 1:
Assertions.assertEquals(envelopePart.getMinSpeed(), envelopePart.getMaxSpeed());
assertEquals(envelopePart.getMinSpeed(), envelopePart.getMaxSpeed());
expectedEnvelopePartEnergy = testRollingStock.getRollingResistance(envelopePart.getBeginSpeed())
* envelopePart.getTotalDistance();
break;
Expand All @@ -109,4 +114,43 @@ void testGetMechanicalEnergyConsumed() {
assertEquals(expectedEnvelopePartEnergy, envelopePartEnergy, 0.1 * expectedEnvelopePartEnergy + 1000);
}
}

@ParameterizedTest
@MethodSource("minEnvelopePartsArgs")
void testMinEnvelopeParts(
EnvelopePart envelopePartA,
EnvelopePart envelopePartB,
Stream<Integer> intersectionIndexes,
double[] expectedPositions,
double[] expectedSpeeds) {
var minEnvelope = EnvelopePart.min(
envelopePartA, envelopePartB, Collections.singleton(envelopePartA.getAttr(EnvelopeProfile.class)));
var resultingPositions = minEnvelope.clonePositions();
var resultingSpeeds = minEnvelope.cloneSpeeds();

intersectionIndexes.forEach(intersectionIndex -> {
var intersectionPosition = resultingPositions[intersectionIndex];
var intersectionSpeed = resultingSpeeds[intersectionIndex];
assertEquals(envelopePartA.interpolateSpeed(intersectionPosition), intersectionSpeed);
assertEquals(envelopePartB.interpolateSpeed(intersectionPosition), intersectionSpeed);
});
assertThat(expectedPositions).isEqualTo(resultingPositions);
assertThat(expectedSpeeds).isEqualTo(resultingSpeeds);
}

static Stream<Arguments> minEnvelopePartsArgs() {
return Stream.of(
Arguments.of(
EnvelopeTestUtils.generateTimes(new double[] {0, 10, 20}, new double[] {100, 50, 0}),
EnvelopeTestUtils.generateTimes(new double[] {5, 14, 16}, new double[] {150, 50, 0}),
Stream.of(4),
new double[] {0, 5, 10, 14, 15, 16},
new double[] {100, 79.05694150420949, 50, 38.72983346207417, 35.35533905932738, 0}),
Comment on lines +147 to +148
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks bugged, no?
I would expect the following results (easy to check for 5 and 15):

Suggested change
new double[] {0, 5, 10, 14, 15, 16},
new double[] {100, 79.05694150420949, 50, 38.72983346207417, 35.35533905932738, 0}),
new double[] {0, 5, 10, 14, 15, 16},
new double[] {100, 75, 50, 30, 25, 0}),

Arguments.of(
EnvelopeTestUtils.generateTimes(new double[] {5, 14, 16}, new double[] {100, 50, 0}),
EnvelopeTestUtils.generateTimes(new double[] {0, 10, 20}, new double[] {150, 50, 0}),
Stream.of(1, 4),
new double[] {5, 7.142857142857143, 10, 14, 15, 16},
new double[] {100, 90.63269671749657, 50, 38.72983346207417, 35.35533905932738, 0}));
Comment on lines +153 to +154
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unexpected output here too, if using min(beginPositions) in case of equality (both speeds are 100 at position 5):

Suggested change
new double[] {5, 7.142857142857143, 10, 14, 15, 16},
new double[] {100, 90.63269671749657, 50, 38.72983346207417, 35.35533905932738, 0}));
new double[] {0, 5, 10, 14, 15, 16},
new double[] {150, 100, 50, 30, 25, 0}));

}
}
Loading