memilio icon indicating copy to clipboard operation
memilio copied to clipboard

Use floating point types consistently

Open reneSchm opened this issue 9 months ago • 1 comments

Motivation / Current Behaviour

We did not consistently use ScalarType for some time, but now we added a FP ("floating point") template parameter to several classes. Instead of fixing inconsistencies we now have three contenders for a floating point type, each of which being almost exclusively used as double - hence it is not a current problem, but could definitely become one in the future. Below is an example where each type is present, but probably only one (FP) should be used.

We probably want to use FP where available and ScalarType where it is not. The difficulty is, that there may be components that require double explicitly.

Enhancement description

Decide which type to use where. Then make usage of that type consistent.

Additional context

template <typename FP = ScalarType, class... Categories>
class Populations : public CustomIndexArray<UncertainValue<FP>, Categories...>
{
public:

...

    void set_total(ScalarType value)
    {
        double current_population = get_total();
        if (fabs(current_population) < 1e-12) {
            double ysize = double(this->array().size());
            for (size_t i = 0; i < this->get_num_compartments(); i++) {
                this->array()[(Eigen::Index)i] = value / ysize;
            }
        }
        else {
            for (size_t i = 0; i < this->get_num_compartments(); i++) {
                this->array()[(Eigen::Index)i] *= value / current_population;
            }
        }
    }

...

};

Checklist

  • [X] Attached labels, especially loc:: or model:: labels.
  • [X] Linked to project

reneSchm avatar Apr 25 '24 13:04 reneSchm

Here I will try to list problems of or around as well as implicit requirements on ScalarType or FP. Solutions or alternatives listed are to be understood as ideas, and should be discussed.

  • Requirement: Compatibility with double: We have a lot of double literals in our code (e.g. sin(3.141592653589793 * ((parameters.get<StartDay>() + t) / 182.5 + 0.5));). Hence, we require that ScalarType implements operators for double * ScalarType etc.
    • ~~Alternative: Require that ScalarType can be implicitly cast to and constructed from double~~ This will not work, since it makes binary operators ambiguous (and it may cause loss of data)
    • Alternative: Never use double. Or: Only ever use ScalarType, e.g. via casts or a user-defined literal
  • Problem: Mismatch by defaults: For example, in examples/euler_test.cpp (weird name?) we explicitly use double, except for the EulerIntegratorCore, which defaults to ScalarType. This, in theory, can cause issues.
    • Solution: Do not use defaults for FP
    • Solution: Require perfect Compatibility with double, including Eigen Matrixes and other types using FP or ScalarType internally

reneSchm avatar May 17 '24 14:05 reneSchm