Reimplementing constraint solver

[abdullah]

Trying to rebase, having some difficulties resolving conflicts. See gitter chat.

Sure I help!!!

Oh, I see you have helped yourself

I think that one is the only conflicting change in Sketch.cpp.

The code for setting the bool variable in FCSSketch.cpp is missing, but I would not worry about that now.

Are you planning to restart development?

[DeepSOIC]

Are you planning to restart development?

Well, I hope to.

Yesterday I tried incorporating the malformedConstraint all the way through commit history, but after failing two times, I give up, so I'll just commit the fixes. This is a bit of a bummer, FreeCAD won't compile for quite a range of commits now.

[abdullah]

Are you planning to restart development?

Well, I hope to.

Yesterday I tried incorporating the malformedConstraint all the way through commit history, but after failing two times, I give up, so I'll just commit the fixes. This is a bit of a bummer, FreeCAD won't compile for quite a range of commits now.

Oh! I misunderstood the situation. I saw your commit about the malformedConstraint:
https://github.com/DeepSOIC/FreeCAD-ell ... 2dddbab378

and somehow I misunderstood you had successfully rebased (I am too tired. It is being a looong year and it is only mid-year ).

I try to rebase now (it takes a while to compile when things are not in ccache) ...

[abdullah]

Are you planning to restart development?

Well, I hope to.

Yesterday I tried incorporating the malformedConstraint all the way through commit history, but after failing two times, I give up, so I'll just commit the fixes. This is a bit of a bummer, FreeCAD won't compile for quite a range of commits now.

I have forced-push the branch. It should be correctly rebased now...

[DeepSOIC]

It should be correctly rebased now...

Thanks!

[DeepSOIC]

I am currently fighting with templates and hyperdual numbers. I really want that for splines.
The idea is:
HyperDual is defined as a dual number with template types for real and dual parts. In practice, real part will be Base::DualNumber, and the dual part is either another DualNumber, or a hyperdual number itself.

The main use is to implement a ParaCurve::value method with hyperdual numbers, in order to support automatic computation of derivatives by U. Base::DualNumbers are for tracking solver-parameter derivatives, while the hyperdual's dual part is for tracking derivatives by U.

I'll simplify that, let's forget about derivatives by solver parameters for a moment.
If we compute value(U) with doubles, we get just the value, no derivatives.
If we compute value(HyperDual<double, double>(U, 1.0)), we get first derivative of the value function, in addition to the value of the function.
If we compute value(HyperDual<double, HyperDual<double, double> >(U, HyperDual<..>(1.0, 1.0))), we also get second derivative by U, as value(...).du.du.
And so on, we can implicitly compute higher derivatives, although i think it will blow up into so much code that it gets inefficient.

Now, if I replace all double's here with Base::DualNumber, I magically combine derivatives by U with derivatives by solver parameters.

Now comes the trick question. How do I implement math and operators, that support combinations between HyperDual<>, DualNumber and double? Combinations between different flavors of HyperDuals should not be allowed, but all other combinations are. HyperDuals will always be built from Base::DualNumber, there will never be a double in hyperduals in practice.

I don't have any specific questions "how do I ...", i'm just trying to get a some comprehension of how to implement it. I haven't even ruled out manual definition of the hyperdual classes (it is certainly more writing and a lot of copy-pasting, but feels like it is going to be more readable)



Okay, while I was writing that down, i realized i don't need the type of real part of hyperduals be templated. Only the dual parts.

One specific question that I have is how to prevent C++ from trying to recursively call the function being implemented in situations like this:

Code: Select all

inline DualNumber sqrt(DualNumber x){
    double re = ::sqrt(x.re);
    return DualNumber(re, x.du * 0.5 / re);
}

This time, ::sqrt won't cut it, as I would need to call FCS::sqrt, possibly a (different) template version of one.

[DeepSOIC]

One specific question that I have is how to prevent C++ from trying to recursively call the function being implemented in situations like this:

Code: Select all

inline DualNumber sqrt(DualNumber x){
    double re = ::sqrt(x.re);
    return DualNumber(re, x.du * 0.5 / re);
}

This time, ::sqrt won't cut it, as I would need to call FCS::sqrt, possibly a (different) template version of one.

https://stackoverflow.com/a/23479802 suggests using is the solution, something like this:

Code: Select all

inline DualNumber sqrt(DualNumber x){
    using ::sqrt;
    using FCS::sqrt;
    double re = sqrt(x.re);
    return DualNumber(re, x.du * 0.5 / re);
}

I'll see when I actually write the stuff.

[abdullah]

One specific question that I have is how to prevent C++ from trying to recursively call the function being implemented in situations like this:

Code: Select all

inline DualNumber sqrt(DualNumber x){
    double re = ::sqrt(x.re);
    return DualNumber(re, x.du * 0.5 / re);
}

This time, ::sqrt won't cut it, as I would need to call FCS::sqrt, possibly a (different) template version of one.

https://stackoverflow.com/a/23479802 suggests using is the solution, something like this:

Code: Select all

inline DualNumber sqrt(DualNumber x){
    using ::sqrt;
    using FCS::sqrt;
    double re = sqrt(x.re);
    return DualNumber(re, x.du * 0.5 / re);
}

I'll see when I actually write the stuff.

If the solution above does not work properly to your specific case, and I understood it correctly that you are going to have different "T" types for the recursive exit condition (e.g. T=HyperDual for the first N-1 recursions and T=DualNumber for the last one) you can write two sqrt template functions and use the type T to change the implementation when T=DualNumber. This is there is exactly the same template function definition for both functions, differentiating only in which template is enabled based on the actual type of T.

Something like this (look just at the ! in the "is_base_of" part) is a generalisation of the problem:

Code: Select all

    template <  typename ChildType,
                typename = typename std::enable_if<
                    std::is_base_of<ParaObject, typename std::decay<ChildType>::type>::value
             >::type
    >
    void tieAttr_Child(UnsafePyHandle<ChildType>& ref, std::string name, bool make = false, bool required = true, bool writeOnce = false);

    template <  typename ChildType,
                typename = typename std::enable_if<
                    !std::is_base_of<ParaObject, typename std::decay<ChildType>::type>::value
             >::type
    >
    void tieAttr_Child(UnsafePyHandle<ChildType>& ref, std::string name, bool make = false, bool required = true, bool writeOnce = false);

If you do not really need T type inheritance checking (is_base_of) and (is_same) would do, then you may just use a single templated function removing the type you do not want to have recursively called and then provide a regular member function for that specific type:

Code: Select all

    template <  typename ChildType,
                typename = typename std::enable_if<
                    ! std::is_same<ParaObject, typename std::decay<ChildType>::type>::value
             >::type
    >
    void tieAttr_Child(UnsafePyHandle<ChildType>& ref, std::string name, bool make = false, bool required = true, bool writeOnce = false);

    void tieAttr_Child(UnsafePyHandle<ParaObject>& ref, std::string name, bool make = false, bool required = true, bool writeOnce = false);

So here the template generates all possible required types where T is not ParaObject (in your case this may be something like if it is not HyperDual), and the second function provides the actual implementation when it is indeed HyperDual.

I am a fan of these conditional templates, because they remove the greed of templates to make T everything they can throw in. But maybe it does not match your use case.

Maybe this is not what you were asking, but I hope it is

I am curious what you will come to

[DeepSOIC]

After some thinking, I decided to implement a hyperdual number with a hard-coded fixed number of derivatives baked right into it.

The main reasons for the decision are:
* easy to read and understand, although it's quite heavy math-wise.
* performance: number of operations should be reduced to the absolute minimum required.
* easy to write and debug. Templates are difficult for me. Hyperdual numbers are already a tricky thing, and i don't want the extra headache of template programming right now. Template conversion might be done later, when the thing is proven to work.

So far, looks like this. Not functional yet.

Code: Select all

class FCSExport HyperDualNumber
{
public:
    DualNumber data[HYPERDEPTH] = {0};
    DualNumber& re = data[0];
    DualNumber& du1 = data[1];
public:
    HyperDualNumber(){}
    HyperDualNumber(DualNumber re) {this->re = re;}
    HyperDualNumber(DualNumber re, DualNumber* duals);
    
    HyperDualNumber operator-() const;

    PyObject* getPyObject() const;
    std::string repr() const;
    
    static HyperDualNumber compoundDeriv(HyperDualNumber g_derivs, HyperDualNumber f_derivs);
};


HyperDualNumber operator+(HyperDualNumber a, HyperDualNumber b) {
    HyperDualNumber ret;
    for(int i = 0; i < HYPERDEPTH; ++i)
        ret.data[i] = a.data[i] + b.data[i];
    return ret;
}
inline HyperDualNumber operator+(HyperDualNumber a, DualNumber b) {
    return HyperDualNumber(a.re + b, &a.du1);
}
inline HyperDualNumber operator+(DualNumber b, HyperDualNumber a) {
    return HyperDualNumber(a.re + b, &a.du1);
}

HyperDualNumber operator-(HyperDualNumber a, HyperDualNumber b) {
    HyperDualNumber ret;
    for(int i = 0; i < HYPERDEPTH; ++i)
        ret.data[i] = a.data[i] - b.data[i];
    return ret;
}
inline HyperDualNumber operator-(HyperDualNumber a, DualNumber b) {
    return HyperDualNumber(a.re - b, &a.du1);
}
inline HyperDualNumber operator-(DualNumber b, HyperDualNumber a) {
    return b + (-a);
}

//operator*


//operator/

[DeepSOIC]

As it turns out, nested template-based dual numbers are already done:
https://gitlab.com/tesch1/cppduals/-/tree/master/
I am not sure if I want to use the library, or rather write my own adaptation. If I want to use it, the license question comes up: is Mozilla Public License 2 compatible with FreeCAD?
There are reasons to use the library, such as bugfixes, and vectorization support (++performance).
There is a big reason to not use the library: it does not offer different types for real and dual part (which i think can improve performance, not really sure...).

After quite a bit of digging the math of chain rule of higher order derivatives, I come to a conclusion that nested duals are probably not that bad in performance. Additionally, writing math functions such as atan2 is easier, as I don't have to compute n derivatives of the function by hand (and when it comes to 2-argument functions like atan2, I don't even know how to write one at all). So I am returning the primary focus to nested duals.