Same code in C++ and Python calculates diff values after a lot of loops, is it the nature of float/double?

Question

I am writing a driver model using C++ and Python to compare the performance. The simulation gives data like width, position, speed, etc. and the driver model does some calculations to decide if it needs to brake or not. Both models have the same variables and calculations, but after looping over 500 times first divergence surface and in the end, result in different brake responses. I am aware that floating point error is a thing, but shouldn't it be the same for both languages? I checked C++ and Python if they have IEEE 754 and it seems to be the case. I will attach my check at the end of the question.

Is there something I don't know about float calculations or I do sth wrong in my codes?

What did I expect: That I get the same response for both driver models.

C++ Code:

void DM::runAccumulator()
{    // storing values from previous run
    double prevVTheta = vTheta;
    double prevVActivation = vActivation;

    // shorter handle for convinience and to stay true to given Modell
    double egoSpeed    = data.egoSpeed;
    double egoAcc      = data.egoAcc;
    double egoLength = data.egoLength;
    double targetSpeed = data.targetSpeed;
    double targetAcc   = data.targetAcc;
    double targetWidth = data.targetWidth;
    double targetLength = data.targetLength;
    
    // calculate positional data
    double targetPosRear = data.targetPosX - targetLength/2.0;
    double egoPosFront = data.egoPosX + egoLength/2.0;
    double distEgotoTarget = targetPosRear - egoPosFront;

    // optical size of target vehicle
    vTheta = 2.0 * std::atan(targetWidth / (2.0 * prevDistEgoToTarget));
    prevDistEgoToTarget = distEgotoTarget;
    vThetaDot = (vTheta - prevVTheta)/timestep;
    vp = vThetaDot/vTheta;
   
    if (log = true)
    {
        this -> file << std::fixed << std::setprecision(2) << data.timestamp << ',' << std::setprecision(6) << targetPosRear << ',' << egoPosFront << ',' << prevDistEgoToTarget << ',' << vTheta << ',' << vThetaDot << ',' << vp << '\n';

    }

Python:

def runAccumulator(self):
        #storing values from previous run
        prevVTheta = self.vTheta
        prevVActivation = self.vActivation

        #shorter handle for convinience and to stay true to given Modell
        egoSpeed    = self.data.egoSpeed
        egoAcc      = self.data.egoAcc
        egoLength = self.data.egoLength
        targetSpeed = self.data.targetSpeed
        targetAcc   = self.data.targetAcc
        targetWidth = self.data.targetWidth
        targetLength = self.data.targetLength

        # calculate positional data
        targetPosRear = self.data.targetPosX - targetLength/2.0
        egoPosFront = self.data.egoPosX + egoLength/2.0
        distEgotoTarget = targetPosRear - egoPosFront

        # optical size target vehicle and looming
        self.vTheta = 2.0*arctan(targetWidth/(2.0*self.prevDistEgoToTarget))
        self.prevDistEgoToTarget = distEgotoTarget
        self.vThetaDot = (self.vTheta - prevVTheta)/self.timestep
        self.vp = self.vThetaDot/self.vTheta

        # Driver Model update
        self.vEpsilon = self.vp - self.vpp 
        activationChange = (self.accuK * self.vEpsilon - self.accuM - self.accuC * prevVActivation) * self.timestep
        self.vActivation = max(0.0, prevVActivation-activationChange)


        # log some parameters
        if self.log == True:
            logString = str(self.data.timestamp) + ','
            logString += str(targetPosRear) + ',' #'{:.7f}'.format(targetPosRear) + ','
            logString += str(egoPosFront) + ',' #'{:.7f}'.format(egoPosFront) + ','
            logString += str(self.prevDistEgoToTarget) + ',' #'{:.7f}'.format(targetPosRear - egoPosFront) + ',';
            logString += str(self.vTheta) + ',' #'{:.7f}'.format(self.vTheta) + ',' 
            logString +=str(self.vThetaDot) + ',' # '{:.7f}'.format(self.vThetaDot) + ','
            logString += str(self.vp) + '\n' #'{:.7f}'.format(self.vp) + '\n'
            self.file.write(logString) 

C++ IEEE 754 Check:

#include <cfloat>
#include <iomanip>
#include <iostream>
 
int main()
{
    int w = 16;
    std::cout << std::left; // std::cout << std::setprecision(53);
#   define COUT(x) std::cout << std::setw(w) << #x << " = " << x << '\n'
 
    COUT( FLT_RADIX        );
    COUT( DECIMAL_DIG      );
    COUT( FLT_DECIMAL_DIG  );
    COUT( DBL_DECIMAL_DIG  );
    COUT( LDBL_DECIMAL_DIG );
    COUT( FLT_MIN          );
    COUT( DBL_MIN          );
    COUT( LDBL_MIN         );
    COUT( FLT_TRUE_MIN     );
    COUT( DBL_TRUE_MIN     );
    COUT( LDBL_TRUE_MIN    );
    COUT( FLT_MAX          );
    COUT( DBL_MAX          );
    COUT( LDBL_MAX         );
    COUT( FLT_EPSILON      );
    COUT( DBL_EPSILON      );
    COUT( LDBL_EPSILON     );
    COUT( FLT_DIG          );
    COUT( DBL_DIG          );
    COUT( LDBL_DIG         );
    COUT( FLT_MANT_DIG     );
    COUT( DBL_MANT_DIG     );
    COUT( LDBL_MANT_DIG    );
    COUT( FLT_MIN_EXP      );
    COUT( DBL_MIN_EXP      );
    COUT( LDBL_MIN_EXP     );
    COUT( FLT_MIN_10_EXP   );
    COUT( DBL_MIN_10_EXP   );
    COUT( LDBL_MIN_10_EXP  );
    COUT( FLT_MAX_EXP      );
    COUT( DBL_MAX_EXP      );
    COUT( LDBL_MAX_EXP     );
    COUT( FLT_MAX_10_EXP   );
    COUT( DBL_MAX_10_EXP   );
    COUT( LDBL_MAX_10_EXP  );
    COUT( FLT_ROUNDS       );
    COUT( FLT_EVAL_METHOD  );
    COUT( FLT_HAS_SUBNORM  );
    COUT( DBL_HAS_SUBNORM  );
    COUT( LDBL_HAS_SUBNORM );
}

EDIT: Minimal Working examples (Removed the logging but left the parsing, since I believe it is good to have he actual data from the sim)

It has been some time, since I tried different ways to work around this issue. Sadly, everything I tried didn't work out, so here I am finally following up with the minimal working example. There is one .csv file with the data from the simulation and a minimal version for both, C++ and Python. To clarify, I have the driver model and would like them to break at the same time. For this, vp, vTheta and vThetaDot should be roughly the same. As it is now, vp is 0.0825 in Pyhton and 0.00054 in C++...

minimal working examples: accuX are constants liek 1.5 Python

from numpy import nan, arctan2
from csv import DictReader


class DM():
    def __init__(self):
        self.vTheta: float = nan
        self.vThetaDot: float = nan
        self.prevDistEgoToTarget: float = nan
        self.timestep: float = 0.01
        self.vpp: float = nan

    def runAccumulator(self, data):
            #storing values from previous run
            prevVTheta = self.vTheta
            
            length = 5.4
            width = 2.0

            # calculate positional data
            targetPosRear = float(data["targetPosX"]) - length/2.0 # adapt arctan2 input
            egoPosFront = float(data["egoPosX"]) + length/2.0
            distEgotoTarget = targetPosRear - egoPosFront

            # optical size target vehicle and looming
            self.vTheta = 2.0*arctan2(width,2.0*self.prevDistEgoToTarget)
            self.vEpsilon = self.vp -self.vpp

            
            self.vThetaDot = (self.vTheta - prevVTheta)/self.timestep
            self.vp = self.vThetaDot/self.vTheta 
            self.prevDistEgoToTarget = distEgotoTarget
            activationChange = (self.accuK * self.vEpsilon - self.accuM - self.accuC * prevVActivation) * self.timestep
            self.vActivation = max(0.0, prevVActivation + activationChange)
            


def main():
    with open("../log/sample.csv", 'r') as f:
        dictReader = DictReader(f)     
        listOfDict = list(dictReader)
    
    dm = DM()
    for i in range(700):
        dm.runAccumulator(listOfDict[i])

if __name__ == "__main__":
    main()

C++ minimal working example

#include <limits>
#include <cmath>
#include <fstream>
#include <sstream>
#include <iostream>
#include <vector>

class DM
{
private:
    /* data */
public:
    DM(/* args */);
    ~DM();
    double vTheta = std::numeric_limits<double>::quiet_NaN();
    double vThetaDot = std::numeric_limits<double>::quiet_NaN();
    double prevDistEgoToTarget = std::numeric_limits<double>::quiet_NaN();
    double timestep = 0.01;
    double vp = std::numeric_limits<double>::quiet_NaN();
    void runAccumulator(std::vector<std::string> data);

    std::ifstream sampleCSV;
    void readFileIntoString(const std::string& path);
};

DM::DM(/* args */)
{
}

DM::~DM()
{
}

void DM::runAccumulator(std::vector<std::string> data)
{
    double prevVTheta = vTheta;
    double length = 5.4;
    double width = 2.0;

    // calculate positional data
    double targetPosX = stod(data[3]) - length/2.0;
    double egoPosX = stod(data[2]) + length/2.0;
    double distEgoToTarget = targetPosX - egoPosX;

    vTheta = 2.0* std::atan2(width,distEgoToTarget);
    vThetaDot = (vTheta-prevVTheta)/timestep;
    vp = vThetaDot/vTheta;
    prevDistEgoToTarget = distEgoToTarget;
    vEpsilon = vp -vpp
    std::cout << vp << std::endl;
    auto activationChange = (accuK * vEpsilon - accuM - accuC*prevVActivation) * timestep;
    vActivation = std::max(0.0, prevVActivation + activationChange);

}

void DM::readFileIntoString(const std::string& path)
{
    this -> sampleCSV.open(path);
    if (!sampleCSV.is_open()) {
        std::cerr << "Could not open the file - '"
             << path << "'" << std::endl;
        exit(EXIT_FAILURE);
    }
    std::string temp;
    getline(this->sampleCSV, temp);
}

std::vector<std::string> split (std::string s, std::string delimiter) {
    size_t pos_start = 0, pos_end, delim_len = delimiter.length();
    std::string token;
    std::vector<std::string> res;

    while ((pos_end = s.find (delimiter, pos_start)) != std::string::npos) {
        token = s.substr (pos_start, pos_end - pos_start);
        pos_start = pos_end + delim_len;
        
        res.push_back(token);
    }

    res.push_back(s.substr (pos_start));
    return res;
}

int main(int argc, char const *argv[])
{
    DM dm;
    dm.readFileIntoString("../log/sample.csv");
    std::string buffer;
    std::vector<std::string> linevec;
    
    
    for (size_t i = 0; i < 700; i++)
    {
        getline(dm.sampleCSV, buffer);
        linevec = split(buffer, ",");
        dm.runAccumulator(linevec);
    }
    
    
    
    return 0;
}

CSV: https://pastebin.com/8gBm5z8q

Answer 1

As a rule of thumb, if you every find yourself in a situation, where your code does not work because of floating point arithmetic, it is very likely that at least one of the following sentences applies to you:

• You work in a very niche field of research.
• You have a bug in your code.
• You are dealing with a mathematically ill-posed problem.

In your case, you have a bug in your code. Your C++ code and your Python code do not do the same thing:

• You have an index error (you should get targetPosX from data[4] not from data[3])
• In your Python code, the second argument to the arctan2 function is 2.0*self.prevDistEgoToTarget, while in the C++ code it's distEgoToTarget.
• Your Python code does one more iteration than your C++ code.

These errors can easily be found by printing all intermediate values in just one iteration. You don't even need a debugger for this. If you fix these issues, you get a vp value of 0.08624797012674965 from Python and 0.086247970126749646 from C++ after 700 iterations. I would argue that the difference is negligible.

Answer 2

IEEE754 guarantees a few things such as:

• The format of FP values
• Basic operations (+,-,*,/ and sqrt) are exact. The result is as if you did the calculation with infinite precision, and then rounded the result.

What IEEE754 does not guarantee is that other functions such as arctan are rounded in the same way. And therefore, Python and C++ can reasonably differ by 1 ULP.

[background] The problem with "exact rounding" is that you need an algorithm to do it correctly. For the 5 basic functions, there were well-known, efficient algorithms at the time IEEE754 was drafted. Since then, efficient algorithms have been discovered for more functions, but there are many, many functions that still do not have such algorithms. In particular, many functions still use tables and interpolation, and that's almost automatically inexact.