Upgrade PRIMA.jl for the future PRIMA v0.8.0

[amontoison]

@emmt @zaikunzhang
Related issue: #28

I cross-compiled the version 0.8.0 of PRIMA and pushed the
generated PRIMA_jll.jl on my fork:
https://github.com/amontoison/PRIMA_jll.jl

I used gen/wrapper.jl and Clang.jl to regenerate src/wrappers.jl.
We have the new Julia wrappers for the updated C interface.
We can now work on src/PRIMA.jl.

Continuous integration will use my fork of PRIMA_jll.jl.
I only cross-compiled PRIMA on Linux / Windows / Mac Intel but I will add tarballs for other platforms soon (FreeBSD / Apple ARM).

[emmt]

@amontoison @zaikunzhang

I do not like the idea of interfacing in Julia the new high level API of the prima C library. It would be cleaner and simpler to update the prototypes of the different algorithms (formerly accessible via prima_bobyqa, prima_newuoa, etc.) and use Julia types to make the high-level API. First, it will involve much less changes in the Julia code (thus with less new bugs). Second, it will be easier to maintain the same Julia API for the end-user which is of uttermost importance for the other packages that depend on PRIMA.jl.

I realise that the C interface to the different base algorithms (prima_bobyqa, prima_newuoa, etc. formely in c/include/prima/prima.h) has disappeared in the main branch of libprima. This is a real concern for me as it makes writing wrappers in other languages potentially more difficult because in addition to wrapping function calls and passing simple arguments like values, pointers to functions or to arrays, the high level API imposes to deal with C structures which can be quite complex. It is not that the high level interface is bad (it really simplifies things for a C user) but the other functions (prima_bobyqa, prima_newuoa, etc.) have their utility. I have made this PR to restore the visibility of the base algorithms in the C interface.

[amontoison]

One alternative is to use a patch based on your PR @emmt when we compile the new release 0.8.0 with Yggdrasil.
The previous C functions will be visible again but only for us.
It doesn't impact the other prima users.

[emmt]

This is a very good idea and would be very nice indeed.

In the mean time, I have started the conversion of PRIMA.jl, just the wrapper code (this can easily be done by hand). You can see it as commit fc2fed9 in the eric-devel branch.

[amontoison]

Can you push on this branch? I can easily regenerate amontoison/PRIMA_jll.jl, which allow us to test the Julia interface with this PR.

[zaikunzhang]

Hi @amontoison and @emmt

One alternative is to use a patch based on your PR @emmt when we compile the new release 0.8.0 with Yggdrasil. The previous C functions will be visible again but only for us. It doesn't impact the other prima users.

The signatures of individual C functions such as prima_newuoa will not be stable. This is because not all arguments of the corresponding Fortran subroutines have been exposed (compare the signatures of newuoa_c.f90 and newuoa.f90). The signature of prima_newuoa will change if we want to expose more arguments in the future.

In contrast, the signature of the C function prima is expected to remain unchanged, namely

int prima_minimize(const prima_algorithm_t algorithm, const prima_problem_t problem, const prima_options_t options, prima_result_t *const result);

When we want to expose more arguments of the Fortran subroutines, we will just enrich the fields of options.

This is another reason why prima_newuoa etc will not be exposed in the C API (they are unstable), and why using prima_newuoa etc should be avoided.

Thanks.

[emmt]

Can you push on this branch? I can easily regenerate amontoison/PRIMA_jll.jl, which allow us to test the Julia interface with this PR.

This is the eric-devel branch. Can't you see it?

[emmt]

If more options are added to the C structure, it means that the code wrappers will have to be updated just as if other arguments are added to the base algorithms with, again, the additional difficulty to deal with a C structure. In practice, this is feasible in Julia, but not all types of C structures can be reliably directly interfaced (believe me, I would have preferred it to be possible, as it would have greatly simplified my life in a number of applications).

To solve this technical problem there is another possibility that can satisfy everyone if we can simply manage the C structure by using mutators/accessors. Typically, the C API could provide:

• a function, say prima_problem_create, to create the structure context with sensitive default parameters (for example assuming unconstrained problem);
• mutators to configure settings: one to define linear constraints, one to define non-linear constraints, etc.
• a function to solve the problem, say prima_problem_solve;
• accessors to retrieve specific output parameters: the termination status, the number of function calls, the level of constraints violation, etc. like prima_problem_get_status,
• a function, say prima_problem_delete, to release the context and associated resources.

This is commonly used by many C libraries to facilitate bindings with other langages and yet allow for the possibilities offered by the library to be enriched as the code evolves. If these functions deal with simple arguments like integer or floating-point values, pointers to arrays of such values or to opaque structures, then it is easy to write bindings and to keep a similar high level interface.

I am confident that such an API can be quickly written on top of the current C API with no or only slight changes. The structure describing the PRIMA problem can be kept public for codes that prefer a direct access. In other words any of the two types of use (accessing directly the structure or considerenig it as opaque and using accessors and mutators) can co-exist.

I thank that we can quickly converge to a solution that satisfies everyone's needs and desires and that is maintainable to the long term. I am ready to draft a first version of such a C API.

[zaikunzhang]

If more options are added to the C structure, it means that the code wrappers will have to be updated just as if other arguments are added to the base algorithms with, again, the additional difficulty to deal with a C structure.

I am not sure about this, due to my ignorance about C and Julia.

Take the structure prima_options_t as an example. It contains the options/parameters that will be used by the solvers, including rhobeg, rhoend, maxfun etc. The wrapper should use prima_init_options to initialize it, which aligns with your suggestions below. The code will be something corresponding to

    prima_options_t options;
    prima_init_options(&options);
    options.iprint = PRIMA_MSG_EXIT;
    options.rhoend = 1e-6;
    options.maxfun = 500*n;
    options.callback = &callback;

Now assume that a new field magic_parameter is added to prima_options_t, and assume that prima_init_options has been updated accordingly, initializing magic_parameter to its default value. Then we have the following.

• Even if the wrapper is not updated, the code is not broken. magic_paramter will be initialized by prima_init_options to the default value. In contrast, if we did not use a structure to provide the options but listed them one by one as arguments for prima_xxx, the wrapper code would be broken when the C binding adds a new option.
• If we want to update the wrapper to make use of options.magic_parameter, we only need to add one line (or maybe one block) of code that sets options.magic_paramater to a particular value we want. I imagine that the amount of work needed is not more than the situation where we listed all options one by one as arguments for prima_xxx.

To solve this technical problem there is another possibility that can satisfy everyone if we can simply manage the C structure by using mutators/accessors.

I agree with this methodology. See below.

• a function, say prima_problem_create, to create the structure context with sensitive default parameters (for example assuming unconstrained problem);

If my understanding is correct, you mean the prima_init_problem as follows:

https://github.com/libprima/prima/blob/ae9c559c44cda3821471485a5f24e9b8b1796f07/c/include/prima/prima.h#L184-L186

See also prima_init_options:

https://github.com/libprima/prima/blob/ae9c559c44cda3821471485a5f24e9b8b1796f07/c/include/prima/prima.h#L242-L244

• mutators to configure settings: one to define linear constraints, one to define non-linear constraints, etc.

I am not sure about this one. Is it still needed given prima_init_problem?

• a function to solve the problem, say prima_problem_solve;

I suppose it is prima_minimize as follows:

https://github.com/libprima/prima/blob/ae9c559c44cda3821471485a5f24e9b8b1796f07/c/include/prima/prima.h#L184-L186

• accessors to retrieve specific output parameters: the termination status, the number of function calls, the level of constraints violation, etc. like prima_problem_get_status,

I am not sure why we would need such a function to get the termination status etc. What would the function do other than reading directly result.status etc (see below)? Sorry for my ignorance about C and Julia.

https://github.com/libprima/prima/blob/ae9c559c44cda3821471485a5f24e9b8b1796f07/c/include/prima/prima.h#L247-L274

• a function, say prima_problem_delete, to release the context and associated resources.

There is prima_free_result:

https://github.com/libprima/prima/blob/ae9c559c44cda3821471485a5f24e9b8b1796f07/c/include/prima/prima.h#L277-L279

But we do not have prima_free_problem. Should we have it? @nbelakovski

See the following examples about the usage of these functions:

https://github.com/libprima/prima/blob/main/c/examples/newuoa/newuoa_example.c
https://github.com/libprima/prima/blob/main/c/examples/cobyla/cobyla_example.c

This is commonly used by many C libraries to facilitate bindings with other langages and yet allow for the possibilities offered by the library to be enriched as the code evolves. If these functions deal with simple arguments like integer or floating-point values, pointers to arrays of such values or to opaque structures, then it is easy to write bindings and to keep a similar high level interface.

I am confident that such an API can be quickly written on top of the current C API with no or only slight changes. The structure describing the PRIMA problem can be kept public for codes that prefer a direct access. In other words any of the two types of use (accessing directly the structure or considerenig it as opaque and using accessors and mutators) can co-exist.

I agree with the idea. However, due to my ignorance, I am not sure what "The structure describing the PRIMA problem can be kept public" and "two types of use (accessing directly the structure or considerenig it as opaque and using accessors and mutators) can co-exist" mean. Is it considered abnormal to have both types of use co-exist? Anyway, I suppose we need the help of @nbelakovski here. He developed the current version of the C binding, based on the previous one by @jschueller (many thanks to both!).

I thank that we can quickly converge to a solution that satisfies everyone's needs and desires and that is maintainable to the long term.

I hope so!

Thank you very much!

[amontoison]

If the C structures of PRIMA only contains "bit" type like int, float or opaque structure (Ptr{Cvoid}) we can easily interface them in Julia and we don't need additional setter / getter in the C interface.

The main issue is when we deal with String (vector of char) or any other vector of int / float.
Because these vectors have a specific length in C / Fortran, the equivalent type in Julia is NTuple and not Vector. An NTuple is not mutable in Julia and in that case we need a setter in the C API to modify it.

The Julia structure is a "mirror" of the C structure and must be aligned (same number of bits).

An example to illustrate the issue is when you want to solve a linear optimization problem:
min c'x s.c Ax = b and x >= 0.
If you have a C routine that take as input c, A and b, it's easy to interface because we will pass the pointer of these arrays to the C routine.
If instead the C routine takes as input a structure with (c, A, b) we need to convert the arrays (mutable) to some NTuple (immutable) and we are unable to update them in Julia direcly.

[zaikunzhang]

If the C structures of PRIMA only contains "bit" type like int, float or opaque structure (Ptr{Cvoid}) we can easily interface them in Julia and we don't need additional setter / getter in the C interface.

The main issue is when we deal with String (vector of char) or any other vector of int / float. Because these vectors have a specific length in C / Fortran, the equivalent type in Julia is NTuple and not Vector. An NTuple is not mutable in Julia and in that case we need a setter in the C API to modify it.

prima_problem_t has only int, float, and pointers:

https://github.com/libprima/prima/blob/56b6489b91e7cbd9a400da87a915d41d9bd7545e/c/include/prima/prima.h#L121-L186

prima_options_t has only int, float, and pointers:

https://github.com/libprima/prima/blob/56b6489b91e7cbd9a400da87a915d41d9bd7545e/c/include/prima/prima.h#L194-L244

prima_result_t has int, float, pointers, and a string (message):

https://github.com/libprima/prima/blob/56b6489b91e7cbd9a400da87a915d41d9bd7545e/c/include/prima/prima.h#L252-L279

[zaikunzhang]

I think we must invite @nbelakovski to the conversation because he is the one who is more familiar with the current C API. My knowledge of C is very limited.

[emmt]

That a lot of reactions! I'll try to answer to the most important points to clarify what I am suggesting.

First of all, I am not saying that the current C API is wrong. I am just suggested it has to be augmented by a few functions and perhaps modified a bit to fit our needs.

I am not sure about this one. Is it still needed given prima_init_problem?

It is almost true except that it should takes pointer of structures as arguments, not structure instances.

I am not sure why we would need such a function to get the termination status etc. What would the function do other than reading directly result.status etc (see below)? Sorry for my ignorance about C and Julia.

This is exactly what I want to avoid (and the problem is not specific to Julia) because it is not always possible/safe to deal with C structures in an another language. This requires to know exactly how the members of the structure are layered, at which memory offsets, etc. and this may depend on the compiler. Even though there maybe there is some standadized ABI, the foreign langage may not offer the flexibility to read/write at a given memory address.

If the C structures of PRIMA only contains "bit" type like int, float or opaque structure (Ptr{Cvoid}) we can easily interface them in Julia and we don't need additional setter / getter in the C interface.

Julia have structures that can be stored in memory more or less like C structures with same alignment and offsets (provided these structures are immutable) but I had issues with some such C structures for which C and Julia use different offsets (it may be due to Clang.jl incorrectly computing these offsets). This is why I suggest using mutators and accessors. This is a common solution to facilitate interfacing C libraries with other languages and this will thus be beneficial to other langage bindings than Julia. In practice, this poses less portability issues because you do not have to recompute offsets, just interface the mutator(s)/accessor(s) for the new structure member(s).

The main issue is when we deal with String (vector of char) or any other vector of int / float. Because these vectors have a specific length in C / Fortran, the equivalent type in Julia is NTuple and not Vector. An NTuple is not mutable in Julia and in that case we need a setter in the C API to modify it.

Strings are indeed another thing that may pose problem. But there is a clean solution in Julia which is fully supproted by @ccall and which is used in the current interface to the prima_get_rc_string function. This can be directly recycled, no needs to change.

The Julia structure is a "mirror" of the C structure and must be aligned (same number of bits).

As I said above, this is doable in Julia but I would rather avoid it.

prima_problem_t has only int, float, and pointers:

prima_options_t has only int, float, and pointers:

Yes this is true but it is not the types of the members that is the issue it is how they are stored in the structure.

What I propose (once again this is a very common solution) is something like:

typedef struct { ... /* members stored in the context */  } prima_context_t;

/* get address of a new instantiated context (considered as an opaque structure in what follows) */
prima_context_t* ctx = prima_context_new(... /* arguments to initialize the context */);

/* set options, and problem constraints with "mutators" */
prima_context_set_bounds(ctx, lo, up);
prima_context_set_opt1(ctx, val1);
prima_context_set_opt2(ctx, val2);
...

/* solve the problem */
prima_context_solve(ctx);

/* retrieve parameters with accessors */
prima_rc_t status = prima_context_get_status(ctx);
double fx = prima_context_get_fval(ctx);
...

/* release context structure and associated ressources */
prima_context_free(ctx);

This kind of API can sit on top of the existing C API (which is simpler for a C user) but implementing it may be the opportunity to improve some aspects:

• pass structure pointers not instances (to avoid an unnecessary copy and to share the same functions in the 2 worlds);
• group prima_problem_t and prima_options_t in a unique structure (which I called prima_context_t above but this is a just to illustrate the specific name is irrelevant).

As I proposed, I can quite quickly write this C API (with better function names) in another branch to illustrate what can be done and so that we can discuss on something concrete (as we did about the PRIMA.jl interface). We could decide later to merge.

[emmt]

I think we must invite @nbelakovski to the conversation because he is the one who is more familiar with the current C API. My knowledge of C is very limited.

Of course.

[emmt]

I have done this PR (libprima/prima#189) to serve as an example.

[nbelakovski]

Thanks for inviting me to the discussion. I have a couple questions. Keep in mind I'm not familiar with Julia.

If more options are added to the C structure, it means that the code wrappers will have to be updated just as if other arguments are added to the base algorithms with, again, the additional difficulty to deal with a C structure.

It's not clear to me what the difficultly is in dealing with C structures in Julia, could you provide a simple example to illustrate this? Keep in mind that PRIMA is meant only for the 5 algorithms already contained within it, new algorithms will not be added, and while certain implementation details may change the algorithms themselves largely (or entirely) will not, so the interface we have now is largely stable, save for adding a few more int/float options here and there.

2. I took a look at the proposed alternative C API in Alternative C API prima#189. Does this have to be merged into the main PRIMA repo or can it be implemented here? It seems that prima_context_solve just calls prima_minimize so I don't think any internal implementation details are needed (save for prima_init_result, but the code can be rearranged in such a way as to make it unnecessary to call that function explicitly, since this is ordinarily handled in prima_minimize)

[emmt]

As I said, handling C structures is doable in Julia but has a number of issues:

1. Fields may be stored at a different offset in Julia than it is in C (I had this once).
2. References to objects (e.g. pointers to arrays elements like the bounds of a constrained problem) must be carefully managed to prevent the related object to be garbage collected while its reference is in use by the called C function. This is possible with GC.@preserve but harder to make it correct (especially if there are several such objects and if some are optional). When the objects themselves are directly passed to ccall or @ccall to call a C function where a pointer is expected, not only ccall takes care of taking the relevant pointer but also of protecting the related object to not be garbage collected until the C function returns.

All these disappear when calling a function that only takes bare numerical values, address of such values, or pointers to arrays of such values. This is used by many C libraries to facilitate bindings in other languages. What I propose in libprima/prima#189 can sit next to the existing API, the two APIs can co-exist. In this PR, the only change to the existing API is to pass structures by address not by value (to save stack space), but this is not mandatory to implement the alternative API.