GCFM miscalculations when not initialised properly #1438

chraibi · 2025-01-16T19:18:37Z

I found a bug in the calculation of the ellipse distances.
Since orientation is a unit vector, it directly provides the cosine and sine of the rotation angle.
However, at first the orientation might be a zero vector.

In this PR, I will check all equations related to ellipses and update their documentation.

Hopefully, by then we will have a correctly performing gcfm.

libsimulator/src/Point.cpp

Ozaq

Hey @chraibi if I understand it correctly the issue only manifests when adding agents, when their orientation vector is erroneously [0,0], right?

If the above is the case, I think it would make more sense to require users to provide a unit length orientation vector as part of the construction? And fail on construction, i.e as early as possible.

Further, the change you made would still allow users to set an initial orientation of [0,0] but would abort(!) the program at runtime. Not what we generally want as this kills the whole process. Instead we should throw an exception and always do the check, this way we can properly inform the user.

Also: could you please add a regression test that exhibits the problem? Let me know if I should lend you a hand in that.

Co-authored-by: Kai Kratz <[email protected]>

chraibi · 2025-01-17T14:50:29Z

@Ozaq

If the above is the case, I think it would make more sense to require users to provide a unit length orientation vector as part of the construction? And fail on construction, i.e as early as possible.

Yes, that's what I just realized. However, GCFM is the only model that works with Ellipses. Since, we initialize the orientation to (0,0) and this means that cos(phi)=0 and sin(phi)=0 at the same time.

You are right. The problem can easily be solved by having the user initialize the orientation.

Without initialization

With initialization

I personally prefer not to make these calls and rather validate the function input like so:

{
    if((std::abs(cphi) < 1e-6 && std::abs(sphi) < 1e-6) ||
       std::abs(cphi * cphi + sphi * sphi - 1.0) > 1e-6) {
        throw std::invalid_argument(
            "Invalid rotation inputs: vector is zero or not normalized (cphi^2 + sphi^2 != 1). "
            "Values: cphi=" +
            std::to_string(cphi) + ", sphi=" + std::to_string(sphi));
    }

Further, the change you made would still allow users to set an initial orientation of [0,0] but would abort(!) the program at runtime. Not what we generally want as this kills the whole process. Instead we should throw an exception and always do the check, this way we can properly inform the user.

I did not try yet to solve the problem.
In this Draft PR, I'm just trying, very cautiously, to go through the functions and update their documentations.

Also: could you please add a regression test that exhibits the problem?

Yes, I updated one test by adding a second agent. This already produces the problem.

When using GCFM, agent orientation represents sin/cos values for ellipse orientation, but defaults to (0,0) in constructor. This causes issues since (0,0) is an invalid sin/cos pair. If orientation remains (0,0), assume it wasn't set by user and automatically initialize it to point towards agent's target position.

Added conditional checks to set the scaling factor to 1 when model1.v0 or model2.v0 is zero.

in the case of two agents on the line: One moving toward the other, who is not moving, we have the following situation: - The driving force and repulsive force cancel out but don’t completely stop the agent, who's speed is very small but not 0. - So, Agent 1 (very) slowly pushes through Agent 2, which is incorrect. An alternative for choosing such a big value would be to introduce an explicit stopping mechanism e.g: if (netForceMagnitude < 1e-3 && model.speed < 1e-3) { update.velocity = Point(0, 0); } else { update.velocity = (agent.orientation * model.speed) + acc * dT; }

Ozaq · 2025-01-20T20:16:19Z

However, GCFM is the only model that works with Ellipses. Since, we initialize the orientation to (0,0) and this means that cos(phi)=0 and sin(phi)=0 at the same time.

True, but an orientation of (0,0) is also a constraint violation because we say "Orientation is a unit length vector", we could also have used an angle in the API when creating the agent and then convert it into a vector internally. This might be a useful API change for a potential 2.0. What do you think?

TODO: Check for invalid orientations from user during simulation

chraibi · 2025-01-21T19:01:47Z

@Ozaq I find orientation is more intuitive and easier to use than angle. For now, I check for it during the initialization phase.
Just need to figure out how to check for it again during a running simulation when the user decides to change it again.

I think this PR solves already the problem I had before.

Ozaq

Also you did forget to actually change the defaults on the python layer, the change in the C-API layer would always be override by the python layer.

Ozaq · 2025-01-22T08:19:16Z

systemtest/test_model_properties.py

+    agent2 = jps.GeneralizedCentrifugalForceModelAgentParameters(
+        journey_id=journey_id,
+        stage_id=wp,
+        position=(3, 1),
+    )
+
    agent_id = sim.add_agent(agent)
+    sim.add_agent(agent2)


What do you want to test here? You add an additional Agent and then no new checks are added, this looks incomplete.

By placing two pedestrians the test fails.
This reproduces the problem with 0 orientation.

Ozaq · 2025-01-22T08:22:03Z

python_modules/jupedsim/jupedsim/models/generalized_centrifugal_force.py

@@ -29,7 +29,7 @@ class GeneralizedCentrifugalForceModel:
    max_geometry_interaction_distance: float = 2
    max_neighbor_interpolation_distance: float = 0.1
    max_geometry_interpolation_distance: float = 0.1
-    max_neighbor_repulsion_force: float = 9
+    max_neighbor_repulsion_force: float = 30


We should not change the defaults of our models. We have a different issue open with a similar request for change, see #1391. I think that we should make it obvious that defaults change and provide them as selectable presets.

This change should be removed.

Ozaq · 2025-01-22T08:57:39Z

libsimulator/src/Simulation.cpp

+    Point targetPoint = std::get<0>(_journeys.at(agent.journeyId)->Target(agent));
+    if((targetPoint - agent.pos).Norm() < 1e-6) {
+        LOG_WARNING(
+            "Agent's position ({}, {})  is the same as the target ({}, {}). Set defaul orientation "
+            "to (1,0)",
+            agent.pos.x,
+            agent.pos.y,
+            targetPoint.x,
+            targetPoint.y);
+        agent.orientation = Point(1.0, 0.0);
+    }
+
+    if(agent.orientation.isInvalidOrientation()) {
+        Point direction = (targetPoint - agent.pos).Normalized();
+        agent.orientation = direction;
+    }
+    agent.orientation = agent.orientation.Normalized();
+    _operationalDecisionSystem.ValidateAgent(agent, _neighborhoodSearch, *_geometry);
+


I would this whole block to look like this:

if(isZeroLength(agent.orientation)) { throw SimulationError("Invalid orientation: {}", agent.orientation); } agent.orientation = agent.orientation.Normalized();

Note: I am aware that we do not have a 'isZeroLength' function (yet).

Ozaq · 2025-01-22T09:14:29Z

libsimulator/src/Simulation.cpp

@@ -8,6 +8,7 @@
 #include "OperationalModel.hpp"
 #include "Stage.hpp"
 #include "Visitor.hpp"
+#include <Logger.hpp>


Ozaq · 2025-01-22T09:14:45Z

libsimulator/src/Simulation.cpp

@@ -177,17 +178,33 @@ BaseStage::ID Simulation::AddStage(const StageDescription stageDescription)

 GenericAgent::ID Simulation::AddAgent(GenericAgent&& agent)
 {
+


Ozaq · 2025-01-22T09:25:41Z

libsimulator/src/Point.cpp

@@ -1,12 +1,20 @@
 // Copyright © 2012-2024 Forschungszentrum Jülich GmbH
 // SPDX-License-Identifier: LGPL-3.0-or-later
 #include "Point.hpp"
-
-#include "Macros.hpp"
+#include <cassert>


This is unused.

Ozaq · 2025-01-22T09:27:15Z

libsimulator/src/Ellipse.cpp

+/**
+ * @brief Computes the point on the ellipse boundary along the line from the ellipse center to point
+ * P.
+ *
+ * Given a point \( P \) in the local coordinate system of the ellipse, this function finds the
+ * corresponding point on the ellipse that lies on the same line extending from the center through
+ * \( P \).
+ *
+ * **Behavior:**
+ * - If \( P \) is very close to the ellipse center, it defaults to returning the point \((a, 0)\)
+ * on the ellipse.
+ * - Otherwise, it scales the direction from the center to \( P \) by the ellipse’s semi-major and
+ * semi-minor axes.
+ *
+ * @return Point on the ellipse boundary, transformed into the global coordinate system.
+ */


Please follow the established documentation style and use /// to indicate doc-style comments. Please omit the brief because we usually use Doxygen with auto-brief, i.e. First line / first sentence is the brief.

Ozaq · 2025-01-22T09:27:55Z

libsimulator/src/Ellipse.cpp

+/**
+ * @brief Calculates the effective distance between two ellipses.
+ *
+ * This function computes the shortest distance between two ellipses. It does so by:
+ *
+ * 1. **Coordinate Transformation:**
+ *    Transforms the center of each ellipse into the local coordinate system of the other
+ *    using their orientation vectors.
+ *
+ * 2. **Closest Point Determination:**
+ *    Identifies the closest point on each ellipse to the other ellipse in their respective
+ *    coordinate systems.
+ *
+ * 3. **Effective Distance Calculation:**
+ *    Calculates the Euclidean distance between these two closest points on the ellipses.
+ */


Please follow the established documentation style and use /// to indicate doc-style comments.

Ozaq · 2025-01-22T09:28:20Z

libsimulator/src/Ellipse.cpp

+/**
+ * @brief Calculates the semi-minor axis (EB) of an ellipse orthogonal to the direction of velocity.
+ *
+ * This method computes the ellipse's semi-axis length perpendicular to the object's movement,
+ * allowing the ellipse to contract or expand based on the provided scaling factor.
+ *
+ * @param scale A scaling factor typically in the range \f$ [0, 1] \f$.
+ *        - \f$ \text{scale} = 0 \f$ → returns \f$ B_{\text{max}} \f$.
+ *        - \f$ \text{scale} = 1 \f$ → returns \f$ B_{\text{min}} \f$.
+ * @return The computed semi-axis length orthogonal to the velocity direction.
+ *
+ * @note Values of `scale` outside the \f$ [0, 1] \f$ range may produce unexpected results.
+ * @warning No input validation is performed on the `scale` parameter.
+ */


Please follow the established documentation style and use /// to indicate doc-style comments.

Ozaq · 2025-01-22T09:28:39Z

libsimulator/src/Ellipse.cpp

+/**
+ * @brief Calculates the semi-major axis (EA) of an ellipse based on the given speed.
+ *
+ * The ellipse adapts dynamically depending on the agent's speed. This method computes
+ * the semi-axis length in the direction of movement, reflecting how the ellipse
+ * elongates as speed increases.
+ *
+ * @param speed The current speed of the object (must be non-negative).
+ * @return The computed semi-axis length in the velocity direction.
+ *
+ */


Please follow the established documentation style and use /// to indicate doc-style comments.

Update documentation of Point functions

45eb9fc

chraibi changed the title ~~Update documentation of Point functions~~ Fix Buggy GCFM Jan 16, 2025

Ozaq reviewed Jan 17, 2025

View reviewed changes

libsimulator/src/Point.cpp Outdated Show resolved Hide resolved

Ozaq requested changes Jan 17, 2025

View reviewed changes

chraibi and others added 5 commits January 17, 2025 14:51

Update libsimulator/src/Point.cpp

2becf11

Co-authored-by: Kai Kratz <[email protected]>

Simplify calculation of Effective distance to ellipse

e938bf7

Format documentation of Ellipse functions

cc253ca

Check angle value before rotating

3022bcd

Update GCFM test by adding two agents

0a39ffb

chraibi changed the title ~~Fix Buggy GCFM~~ GCFM miscalculations when not initialised properly Jan 17, 2025

chraibi added 6 commits January 17, 2025 16:16

Uncomment throw for better debugging

e6cceb0

Fix division by zero in AgentToAgentSpacing by handling zero v0

fe3299a

Added conditional checks to set the scaling factor to 1 when model1.v0 or model2.v0 is zero.

remove not use return value

c3bca76

write a warning instead of throwing an error

bdce3ce

chraibi marked this pull request as ready for review January 17, 2025 20:39

chraibi added 2 commits January 21, 2025 14:30

Check valid orientation when adding agents

dc21dfe

TODO: Check for invalid orientations from user during simulation

default orientation to 1 0

9573d01

Change test since we changed default orientation to 1 0

2793271

Ozaq added this to the v1.3.0 milestone Jan 22, 2025

Set default orientation in python layer

8f3a02b

Ozaq requested changes Jan 22, 2025

View reviewed changes

chraibi added 4 commits January 22, 2025 16:02

Throw explicit error when orientation is not properly set.

92d20f3

unchange default parameter

e755391

check for zero orientation

a23d141

Change commenting style

3d6e84e

chraibi added 4 commits January 22, 2025 17:20

remove unused include

473025d

change iszerolength check position

493e06b

Test orientation only for

7fa88ed

remove unnecessary initialisation

87a1d11

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

GCFM miscalculations when not initialised properly #1438

GCFM miscalculations when not initialised properly #1438

chraibi commented Jan 16, 2025 •

edited

Loading

Ozaq left a comment

chraibi commented Jan 17, 2025 •

edited

Loading

Ozaq commented Jan 20, 2025

chraibi commented Jan 21, 2025

Ozaq left a comment

Ozaq Jan 22, 2025

chraibi Jan 22, 2025

Ozaq Jan 22, 2025

Ozaq Jan 22, 2025

Ozaq Jan 22, 2025

Ozaq Jan 22, 2025

Ozaq Jan 22, 2025

Ozaq Jan 22, 2025

Ozaq Jan 22, 2025

Ozaq Jan 22, 2025

Ozaq Jan 22, 2025

		@@ -177,17 +178,33 @@ BaseStage::ID Simulation::AddStage(const StageDescription stageDescription)

		GenericAgent::ID Simulation::AddAgent(GenericAgent&& agent)
		{

GCFM miscalculations when not initialised properly #1438

Are you sure you want to change the base?

GCFM miscalculations when not initialised properly #1438

Conversation

chraibi commented Jan 16, 2025 • edited Loading

Ozaq left a comment

Choose a reason for hiding this comment

chraibi commented Jan 17, 2025 • edited Loading

Ozaq commented Jan 20, 2025

chraibi commented Jan 21, 2025

Ozaq left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

chraibi commented Jan 16, 2025 •

edited

Loading

chraibi commented Jan 17, 2025 •

edited

Loading