Deprecating SOFA_FLOAT ?

damienmarchal commented 7 years ago

From my perspective maintaining three code base (SOFA_FLOAT, SOFA_DOUBLE, SOFA_FLOAT+SOFA_DOUBLE) and insuring everything is ok on the three code path, costs a lot (amount of code to write, maintain, read (all those #ifdef all around) and compilation time.

So I wonder if it would be acceptable to drop the SOFA_FLOAT part of sofa.

PS: I actually also have the same question for the PS3 #ifdef.

Example of code "readability" improvement:

// Register in the Factory
int ImplicitSurfaceMappingClass = core::RegisterObject("Compute an iso-surface from a set of particles")
#ifndef SOFA_FLOAT
        .add< ImplicitSurfaceMapping< Vec3dTypes, Vec3dTypes > >()
        .add< ImplicitSurfaceMapping< Vec3dTypes, ExtVec3fTypes > >()
#endif
#ifndef SOFA_DOUBLE
        .add< ImplicitSurfaceMapping< Vec3fTypes, Vec3fTypes > >()
        .add< ImplicitSurfaceMapping< Vec3fTypes, ExtVec3fTypes > >()
#endif

#ifndef SOFA_FLOAT
#ifndef SOFA_DOUBLE
        .add< ImplicitSurfaceMapping< Vec3fTypes, Vec3dTypes > >()
        .add< ImplicitSurfaceMapping< Vec3dTypes, Vec3fTypes > >()
#endif
#endif
        ;

#ifndef SOFA_FLOAT
template class SOFA_SOFAIMPLICITFIELD_API ImplicitSurfaceMapping< Vec3dTypes, Vec3dTypes >;
template class SOFA_SOFAIMPLICITFIELD_API ImplicitSurfaceMapping< Vec3dTypes, ExtVec3fTypes >;
#endif
#ifndef SOFA_DOUBLE
template class SOFA_SOFAIMPLICITFIELD_API ImplicitSurfaceMapping< Vec3fTypes, Vec3fTypes >;
template class SOFA_SOFAIMPLICITFIELD_API ImplicitSurfaceMapping< Vec3fTypes, ExtVec3fTypes >;
#endif

#ifndef SOFA_FLOAT
#ifndef SOFA_DOUBLE
template class SOFA_SOFAIMPLICITFIELD_API ImplicitSurfaceMapping< Vec3dTypes, Vec3fTypes >;
template class SOFA_SOFAIMPLICITFIELD_API ImplicitSurfaceMapping< Vec3fTypes, Vec3dTypes >;
#endif
#endif

Would become

int ImplicitSurfaceMappingClass = core::RegisterObject("Compute an iso-surface from a set of particles")
 .add< ImplicitSurfaceMapping< Vec3dTypes, Vec3dTypes > >()
 .add< ImplicitSurfaceMapping< Vec3dTypes, ExtVec3fTypes > >();

template class SOFA_SOFAIMPLICITFIELD_API ImplicitSurfaceMapping< Vec3dTypes, Vec3dTypes >;
template class SOFA_SOFAIMPLICITFIELD_API ImplicitSurfaceMapping< Vec3dTypes, ExtVec3fTypes >;

Suggested labels:

maxime-tournier commented 7 years ago

Even shorter ;-)

int ImplicitSurfaceMappingClass = core::RegisterObject("Compute an iso-surface from a set of particles")
 .add< ImplicitSurfaceMapping< Vec3Types, Vec3Types > >()
 .add< ImplicitSurfaceMapping< Vec3Types, ExtVec3fTypes > >();

template class SOFA_SOFAIMPLICITFIELD_API ImplicitSurfaceMapping< Vec3Types, Vec3Types >;
template class SOFA_SOFAIMPLICITFIELD_API ImplicitSurfaceMapping< Vec3Types, ExtVec3fTypes >;

(the point being to keep SReal customizable instead of forcing it to double)

matthieu-nesme commented 7 years ago

To me, the real question is do we want to keep float+double at the same time? I.e. do we want to be able to have double dofs mapped from float dofs themselves mapped from double dofs?

A simple typedef for SReal would be so easy...

Otherwise the graal would rather be to keep only float everywhere, and having just a few double where it is really necessary for numerical precision.

maxime-tournier commented 7 years ago

@matthieu-nesme that's precisely what I'm hinting to: having a single typedef for SReal and keeping only the needed exact types for e.g. ExtVec3fTypes.

damienmarchal commented 7 years ago

So if I understand you maxime the idea is to have: SReal to be switched via a define Vec3 that is using SReal Vec3d that is using double Vec3f that is using float

And at the same time dropping the SOFA_FLOAT/SOFA_DOUBLE if def that are there to, as mathieu said, mix in the same scenes component in Vec3f & Vec3d.

maxime-tournier commented 7 years ago

@damienmarchal I would drop Vec3d/Vec3f altogether and only use Vec3. There would be no mixing of Vec3f and Vec3d since they would not even exist!

When float or double is needed for a precise use-case one can always use ExtVec3f/ExtVec3d instead (e.g. opengl). It seems to cover all the needs, but maybe I am missing something?

matthieu-nesme commented 7 years ago

ExtVec3f/ExtVec3d is not enough. In the most general case, you could want a part of the simulation in float and another one in double (where the dofs can be anything, vec1, deformation gradients...) Optionally each one with their own solvers, but not necessarily, mappings can map from float to double. I think remembering that this was important for GPU/CUDA simulations to be able to map to float.

It is easy enough to let the code as now, and compiling only in double, or to remove the extra template instantiations to keep only the one based on SReal as suggested by Max with a possibility to let the other ones that are really required.

I repeat the ultimate goal would be to work only with single float, they are more efficient (using less bandwidth, serialization will handle more computations at the same time...) and they are numerically enough in 99.999%. cases We are using double because we never tried to tackle the numerical challenges. For instance https://www.youtube.com/watch?v=etAYNa6DLKw is only using float.

matthieu-nesme commented 7 years ago

@damienmarchal

SReal to be switched via a define Vec3 that is using SReal Vec3d that is using double Vec3f that is using float

That is exactly the case right now! Vec3fTypes/Vec3dTypes/Vec3Types - Rigid3fTypes/Rigid3dTypes/Rigid3Types

maxime-tournier commented 7 years ago

@matthieu-nesme

In the most general case, you could want a part of the simulation in float and another one in double (where the dofs can be anything, vec1, deformation gradients...)

But do we actually? Did we ever use this feature in the past, apart from mapping external vectors (as I said, for Cuda, OpenGL and the likes)?

matthieu-nesme commented 7 years ago

Who is we? The feature was used ; removing it has already been discussed many times and refused.

I think the problem was: our code (e.g. CGSolver) is not numerical robust in single precision so we need to use double precision implementations, but some part of the code is (was) single precision only (e.g. Cuda forcefields).

Once again the real point is why using double precision when single float should be enough? So if you want to save compilation time and have lighter code => instantiate components only on SReal-based types (what it is already done in Flexible for example). If you really want to make a change => clean every components that need to be robust in single float and compile sofa always in float only!

vmagno commented 7 years ago

I just did a quick test in a simple scene (horizontal beam with downward force at one end), and the Newton static solver does not converge with the TetrahedronFEMForceField in single precision, using a CGSolver in double precision. The problem may be with the addDforce function used by CG, but it would need some more investigating. So basically, we can't simply switch to single precision without checking everything that could be sensitive to numerical errors (but I agree that it would be better if we could use mostly single precision).

damienmarchal commented 7 years ago

Thanks all,

Matthieu pointed that "The feature was used ; removing it has already been discussed many times and refused" but having something discussed or used (a long time ago) shouldn't be an argument not to "refresh" the discussion (in the last 10 years things have changed in term of hardware as well as we have now more insight on the intrusiveness of the selected approach). Having the featured re-discussed again is probably an indication that something is somehow problematic :)

As I said before there is fact two different issues...one is about mixing in the same scene object/solver with different floating point representation while the other is about having SReal mapped to float or double by default.

The current approach, despite it does not match individual preferences, makes a relative consensus: SReal to be switched via a define Vec3 that is using SReal Vec3d that is using double Vec3f that is using float We use vec3d or vec3f when we want to be explicit on the type or we use Vec3 if we don't. And of course life would be easier if everyone agreed to use double only code (or according to Matthieu's opinion float only) but I don't think this will happen soon ;)

To me the real problem is not there, it is in how we have implemented the mixing in the same scene object/solver with different floating point representation. Our current implement is based on instantiating templates in the factory for each types so that when the types matches the objects can work together without conversion cost. This design is implemented with the conditional #ifdef SOFA_FLOAT/SOFA_DOUBLE and is very intrusive (our code is harder to write/read/understand, easy to forget things and/or have 'hidden' or subtle bugs (confirmed by the quick test made by @vmagno)) while there is only very specific gain in certain limited use case.

So we are trading run-time speed for a development cost. But the use case is very rare that we may wonder if other approach shouldn't be preferable (eg: a conversion layer instead) ?