Questions about the C interface

[zaikunzhang]

Hi @nbelakovski , I have got two questions regarding the C interface.

1. Why do we copy problem->x0 into result->x in

prima/c/prima.c

Lines 120 to 121 in a617267

	// We copy problem->x0 into result->x so that problem->x0 does not get overwritten by the solver.
	memcpy(result->x, problem->x0, problem->n * sizeof(double));

2. Why do we define use_constr in prima_minimize and pass it to prima_check_problem in

prima/c/prima.c

Line 156 in a617267

int use_constr = (algorithm == PRIMA_COBYLA);

Since prima_check_problem is the only place that uses use_constr, why don't we define it in prima_check_problem?

3. Why do we pass problem and options by pointers to prima_minimize?

prima/c/include/prima/prima.h

Line 241 in 19df480

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

In my opinion, they should be passed by value, as we do not intend to change them within the function --- if that is not our intention, then the possibility should be avoided as much as possible, to avoid accidental mistakes. Is there any disadvantage to pass by values here?

Thank you for taking a look.

[nbelakovski]

I'm just seeing this.

1. From your comments it sounds like you didn't read the comment that's right above the line you're questioning?
2. This was left over from how Julien wrote the code.
3. Again this wasn't my code initially. Looking through the code I see now reason why we need to pass them by pointer, in fact passing them by value is probably a good idea. We can get rid of the null checks.

[zaikunzhang]

1. From your comments it sounds like you didn't read the comment that's right above the line you're questioning?

The comment is exactly my question: why is it an issue if x0 is overwritten? Does this mean that we should check result.x instead of problem.x0 if we need x0 sometimes?

I hope this is not the case in our code because the logic is wrong. Each variable/value should have one and only one meaning. result.x should never mean x0 and tell the user “I am sometimes the value of the returned x, but sometime x0; rely on me for x0 if you know what ‘sometimes’ means.”

Logically, the result is not defined yet after the “initialisation” function is called. Therefore, each component should be set to a value indicating “I am not defined, so don’t use me”. To this end, we should set

• any pointer to NULL,
• any real scalar to NAN,
• and any integer to a strange value, for which I suggest the most negative value that can be contained by the integer, which is what I would do in the Fortran code if needed.

Note that the logic is different from the initialisation of the options, even though they have some superficial similarities:

• The options are (optionally) provided by the user, and its initial status should be “default”;
• The result is to be defined by the solver, and its initial status should be “undefined”.

What do you think?

[nbelakovski]

I'm not sure what you propose? Should we remove x0 from problem_t and have the user input the actual x0 into result.x? That would both eliminate any ambiguity, and also deal with the problem over not overwriting the user provided x0. My main goal in not overwriting problem.x0 was that it wouldn't make a lot of sense if the user set problem.x0 and then after they called the solver x0 had changed and contained the result, that just seemed confusing and ugly.

result.x should never mean x0 and tell the user “I am sometimes the value of the returned x, but sometime x0; rely on me for x0 if you know what ‘sometimes’ means.”

This I don't understand. All the algorithms take x as an argument and expect that it is initially equal to x0. They don't provide an input x0 and a separate output x parameter. It's the same with result.x in the current code.

Logically, the result is not defined yet after the “initialisation” function is called.

This I also don't understand. I think I don't know what you mean by "defined" in this context. The result structure is defined in the header file. It is then initialized by prima_init_result.

I think it should be pointed out that prima_init_result is an internal function not exposed to the user. It is called by prima_minimize and it is not defined in prima.h (meaning it is not possible to call it as a user without a custom build of libprima or by manually declaring it).

[zaikunzhang]

@nbelakovski Reading your response, I believe that we are talking about completely different things ...

result means the result (to be) produced by the solvers. Is this true to you?

[zaikunzhang]

What I will do in prima_init_result is the following.

result->x = NULL
result->f = NAN
result->nlconstr = NULL
result->cstrv = NAN
result->nf = -999999999999999999999999999  // a huge negative value
// ... ... 
// others are similar

For me, this is the only reasonable "initialization" of result. Anything else is unreasonable.

Before the solver finishes its job, the result is not yet defined (mathematically or logically; not in the coding sense). Or in other words, the result does not exist. Its component should not be set to any particular values like x0 or 0. Why should some vector that is not yet defined (mathematically) take the value of x0? Why should some number that is not yet defined (mathematically) take the value 0? Then why not 1? What is particular about 0?

Should we remove x0 from problem_t and have the user input the actual x0 into result.x?

No. This sounds totally illogical to me.

[nbelakovski]

I think this is all very philosophical and if you'd like to change the initialization of result to what you've outlined, feel free to do so. I don't think it will make any difference.

[zaikunzhang]

I think this is all very philosophical

Sure, but this does not make it unimportant. Viewpoints are as important as implementations. Incorrect viewpoints will likely lead to wrong or improper implementations.

and if you'd like to change the initialization of result to what you've outlined, feel free to do so.

I will do it.

I don't think it will make any difference.

I agree, assuming that everything works and there is no bug hidden somewhere. However, it will reduce the potential of having bugs/confusion due to the initialization and make bugs much easier to spot if they occur. This has already happened in one of my tests --- no solver was invoked at all due to a mismatch of the problem of the solver, but result.x and result.f received values that looked totally decent (of course, this is related to #150, which I will address later).

Thank you.

[nbelakovski]

no solver was invoked at all due to a mismatch of the problem of the solver, but result.x and result.f received values that looked totally decent

This sounds like prima_init_result was called due to prima_check_problem returning a 0 error code when it should have returned a non-zero exit code, and this would be the real bug.

This problem currently exists in Fortran as well for x/x0. Since x is set to x0 before calling a solver, it will seem to have a reasonable value even if the solver exit without a clean exit code.

Sure, but this does not make it unimportant.

I think there may be more consequential things to explore, such as setting up a profiler to take a look at the issue raised a couple weeks ago about the performance of PRIMA being substantially worse than another implementation of Powell's algorithms, even after expecting poorer performance due to the focus on code quality. Or working towards a unified interface in Fortran, or working on documentation. The focus on prima_init_result seems like it has less impact than some of those other items.

[zaikunzhang]

This problem currently exists in Fortran as well for x/x0. Since x is set to x0 before calling a solver, it will seem to have a reasonable value even if the solver exit without a clean exit code.

This is not true. There is no solver selection in Fortran. There is no chance for the user to call a solver that does not match the problem.

However, a similar problem does exist under some very marginal case for some other reasons, which I do not want to go in details here.

I think there may be more consequential things to explore, such as setting up a profiler to take a look at the issue raised a couple weeks ago about the performance of PRIMA being substantially worse than another implementation of Powell's algorithms, even after expecting poorer performance due to the focus on code quality. Or working towards a unified interface in Fortran, or working on documentation.

I agree that there are many things to do. The most urgent one is the Python binding and the integration to SciPy. (BTW, for your information, we are getting COBYQA into SciPy).

Or working towards a unified interface in Fortran,

I do not think this is necessary. I have no plan for that. I am super reluctant to make any change to the Fortran code anymore.

working on documentation.

Sure!

the issue raised a couple weeks ago about the performance of PRIMA being substantially worse than another implementation of Powell's algorithms

I suppose you meant #130.

This is not really an issue. In short, we are not concerned about the CPU time in this kind of tests, where the function evaluations take no time. PRIMA is designed for black-box optimization problems where each function evaluation takes hours or days (this is the case for the real applications in Intel, Airbus, and Alibaba, for example). Thus we are only concerned about the number of function evaluations. I hope you can have a look at my response to the issue the explanations here: https://github.com/orgs/libprima/discussions/145

The focus on prima_init_result seems like it has less impact than some of those other items.

It takes a short time as long as I can confirm that what I proposed is not wrong. Again, I do not agree these kinds of changes are less important than those you listed. It avoids potential bugs and the importance is high.

[zaikunzhang]

and if you'd like to change the initialization of result to what you've outlined, feel free to do so.

I will do it.

Done: 488abda , 8d30ef4

[nbelakovski]

I made a comment in the relevant commits. You've introduced a possible memory leak and ignored x0 for all problems. So that's 2 bugs introduced and 0 fixed so far, if you don't mind I'll keep score until we close this issue :)

It might be easier to do this in the form of a PR, instead of an an issue in which you make commits directly to main.