# Maximising a convolution in C++ via a GSL routine

• C/++/#
Gold Member
Consider an integral of the form $$\int_{-1}^1 dx f(x)g(x).$$

I'd like to use https://www.gnu.org/software/gsl/doc/html/min.html to find the maximum of the convolution ##f(x)g(x)## in the domain ##x \in [-1,1]##. The method initiates via a double function with parameters x and a void params.

1) Is it a necessary part of the method to declare a function with such arguments? The function in the link returns cos x + 1. What is the syntax if I want to return the output of another function defined previously in my program?

DrClaude
Mentor
What you are after is not clear. The result of
$$\int_{-1}^1 dx f(x)g(x)$$
is a number, not a function of ##x##, so what is there to minimize. Also, you call ##f(x)g(x)## a convolution, but I don't see what it has to do with ##f * g##, see
http://mathworld.wolfram.com/Convolution.html
1) Is it a necessary part of the method to declare a function with such arguments? The function in the link returns cos x + 1. What is the syntax if I want to return the output of another function defined previously in my program?
You need to have a function that takes in a value of ##x## as an argument and returns a value (double) of based on that ##x##. gsl_min will then try to find a value of ##x## that is a local minimum for that function. How you actually calculate the value to return is of no importance to GSL. It may do so but calling another function you have defined,
C:
double
minimization_function (double x, void *params)
{
return my_func (x);
}

Gold Member
What you are after is not clear. The result of
$$\int_{-1}^1 dx f(x)g(x)$$
is a number, not a function of ##x##, so what is there to minimize. Also, you call ##f(x)g(x)## a convolution, but I don't see what it has to do with ##f * g##, see
http://mathworld.wolfram.com/Convolution.html
I want to maximise the convolution not the integral. f(x) and g(x) represent two kernels and the product is a convolution, language typically used in physics. I didn’t mean the mathematical sense of the word, so sorry for the confusion there.

You need to have a function that takes in a value of ##x## as an argument and returns a value (double) of based on that ##x##. gsl_min will then try to find a value of ##x## that is a local minimum for that function. How you actually calculate the value to return is of no importance to GSL. It may do so but calling another function you have defined,
C:
double
minimization_function (double x, void *params)
{
return my_func (x);
}
thanks, what is the purpose of having a pointer variable params of type void?

DrClaude
Mentor
thanks, what is the purpose of having a pointer variable params of type void?
It is used to pass any necessary parameter to the function. Declaring it void allows one to pass any type, including a struct if more than one parameter is needed. One simply needs to cast it to the correct variable type inside the function. This is used throughout GSL when user-defined functions are needed.

CAF123
Gold Member
Thanks! You mentioned struct in the above message - actually the function I wish to use for maximisation has parameters passed in via a struct, as follows
C:
struct func_params_t { classA *func1;  classB *func2; double c; };

double minimum_function (double x, void* params)
{

double object1;
double func2;
double res;
double f[2];

func_params_t &func_params = *(func_params_t *)(params);

object1 = func_params.func1->FunctionA(x);
func_params.func2->FunctionB(xfunc(x,func_params.c));
func2 = f[1];

res = object1*func2;

return res;
}

The program has #include statements for #include classA.h and #include classB.h so I think *func1 and *func2 are pointer variables to access any object (variable or function) inside classes classA and classB. Is that correct?

Could you explain the line
Code:
 func_params_t &func_params = *(func_params_t *)(params);
?

And in the following lines, object1 is declared to be FunctionA(x), a function within classA.h, here accessed via pointer variable func1. object1*func2 is my f(x)g(x) mentioned in OP.

Last edited:
DrClaude
Mentor
The program has #include statements for #include classA.h and #include classB.h so I think *func1 and *func2 are pointer variables to access any object (variable or function) inside classes classA and classB. Is that correct?
func1 and func2 are pointers to variables of type classA and classB, respectively. From the rest of the code, it does appear that these are structs. These structs have a member, FunctionA (or FunctionB) which is a pointer to a function.

Could you explain the line
Code:
 func_params_t &func_params = *(func_params_t *)(params);
?
[/code]
Is that line correct? I don't see how this could be valid code.

Mark44
Mentor
I want to maximise the convolution not the integral. f(x) and g(x) represent two kernels and the product is a convolution, language typically used in physics.
A convolution is an integral, and is a term used also in mathematics analysis. The convolution of two functions f and g is defined like this:
$$(f * g)(t) = \int_{-\infty}^\infty~f(\tau)g(t - \tau)d\tau$$
##\tau## is a dummy variable -- the convolution is a function of the parameter t.

Other limits of integration are possible.

Mark44
Mentor
BTW, this thread appears to be a continuation of what you were trying to do in your previous thread: https://www.physicsforums.com/threads/use-the-output-of-a-function-in-another-function.970298/ as well as in another thread that appeared earlier: https://www.physicsforums.com/threa...n-a-numerical-integrator.969561/#post-6158795

I'll say here what I said at the end of thread #2:
Mark44 said:
Apparently GSL is Gnu Scientific Library, which I haven't used. From a quick check of the documentation, the GSL adaptive integration routines are numeric integration techniques, meaning that if you supply a function and interval endpoints, they will return a numerical approximation to the integral. IOW, a number.
To calculate the convolution of two functions, you need an integration technique that does symbolic integration, not numeric integration.

In the last post of the first of the three threads, I gave this advice:
Mark44 said:
Two suggestions:
1. Talk to whoever assigned this work to you, and explain that you are unclear on what you are expected to do. Does the person who assigned this work to you have any idea of what all is involved? Would this person be able to do the work himself/herself? If so, ask for guidance, explaining your background. If this person wouldn't be able to do this work, it's unreasonable IMO for you to be saddled with it.
2. It's possible that Mathematica provides "hooks" into its symbolic integration/differentiation routines provided via DLLs (dynamic link libraries). If so, and I really don't know if they do, someone knowledgeable in C++ could write a program that references the relevant DLLs to use the functions in these DLLs.

Last edited:
Mark44
Mentor
Further thoughts.
I'd like to use https://www.gnu.org/software/gsl/doc/html/min.html to find the maximum of the convolution ##f(x)g(x)## in the domain ##x \in [-1,1]##.
I think this is a no-go. The example code at the link above expects a function that can be evaluated at selected points on some interval, with the example function being ##f(x) = \cos(x) + 1##.

In contrast, the function you want to maximize is the convolution of two functions, on what appears to be the interval [-1, 1]. So you're trying to find the maximum value of C(t), with C being defined as follows:
$$C(t) = (f * g)(t) = \int_{-1}^1~f(\tau)g(t - \tau)d\tau$$

So C is the function you want to maximize, but to get C(t) for some specific value of t, you have to evaluate the integral symbolically, which my reading of the GSL docs links that you have provided doesn't do.

I haven't done much checking, but if there is a GSL routine that calculates a convolution symbolically, then that's what you need before you can attempt to find the maximum value. If there isn't, you need to look elsewhere, possibly in a Mathematica DLL, as I suggested in another of your threads.

Short of that, you could divide up the interval [-1, 1] into a bunch of subintervals, find C(t), the convolution of f and g, at, say, the center of each subinterval, store it into an array, and then cycle through the array picking out the largest value. That doesn't seem to have been on your mind, but it seems to me to be the only way to do things, if you can't find a routine that calculates a convolution symbolically.

Last edited:
DrClaude
Mentor
In contrast, the function you want to maximize is the convolution of two functions
I'm not sure about that, see
I want to maximise the convolution not the integral. f(x) and g(x) represent two kernels and the product is a convolution, language typically used in physics. I didn’t mean the mathematical sense of the word, so sorry for the confusion there.

Together with the other threads, this is getting frustrating. The actual problem that needs to be solved has never been explicitly described in details, so we don't know what @CAF123 is actually after. Feels like this is an XY-problem. @Mark44 and I also seem to be pushing the solution in different directions, as we seem to be interpreting it in different ways.

Gold Member
Is that line correct? I don't see how this could be valid code.
I rechecked and it appears to be correct.
A convolution is an integral
I see, in that case, some texts are sloppy because I see things like 'convolution integral'. In any case, what I really want is to maximise the product function f(x)g(x).
BTW, this thread appears to be a continuation of what you were trying to do in your previous thread: https://www.physicsforums.com/threads/use-the-output-of-a-function-in-another-function.970298/ as well as in another thread that appeared earlier: https://www.physicsforums.com/threa...n-a-numerical-integrator.969561/#post-6158795
Yes, I had a private discussion with PeterDonis and he told me to create a new thread.

Together with the other threads, this is getting frustrating. The actual problem that needs to be solved has never been explicitly described in details, so we don't know what @CAF123 is actually after. Feels like this is an XY-problem. @Mark44 and I also seem to be pushing the solution in different directions, as we seem to be interpreting it in different ways.
Hopefully the issue about terminology is clarified now, and it's agreed that gsl routine linked in OP can definitely do what I want. Sorry for lack of details. The product f(x)g(x) is based on previous functions declared elsewhere (as mentioned, g(x) is an interpolation function) so I just wanted to understand the logic behind the gsl method presented, syntax and some of the code I've been given. I'm happy now with the set up of the minimisation function up to the line of code above.

Here is the snippet of code for the gsl minimisation routine. Few queries:

1) What if, in some situation, one doesn't know the expected maximum of the function and the value of x at which that maximum occurs?

2) Could you clarify what is meaning of
Code:
 F.function = &minimum_func;
? Without looking at the gsl header files, I guess it's just a way of accessing the relevant function.

I declare the function via a double and then put the gsl method in main function, as in example code in link. I get an error that, by referencing the parameters via
Code:
 F.params = (void*) (&func_params)
, that 'Use of unidentified identifier 'func_params'.

C:
int main (void) {
int status;
int iter = 0, max_iter = 100;
const gsl_min_fminimizer_type *T;
gsl_min_fminimizer *s;
double m = 15.0, m_expected = 0.0001;
double a = -1.0, b = 1.0;

gsl_function F;
F.function = &minimum_func;
F.params = (void*)(&func_params);

T = gsl_min_fminimizer_brent;
s = gsl_min_fminimizer_alloc (T);
gsl_min_fminimizer_set (s, &F, m, a, b);

printf ("using %s method\n",
gsl_min_fminimizer_name (s));

printf ("%5s [%9s, %9s] %9s %10s %9s\n", "iter", "lower", "upper", "min", "err", "err(est)");
printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n", iter, a, b,m, m - m_expected, b - a);

do{ iter++;
status = gsl_min_fminimizer_iterate (s);
m = gsl_min_fminimizer_x_minimum (s);
a = gsl_min_fminimizer_x_lower (s);
b = gsl_min_fminimizer_x_upper (s);
status= gsl_min_test_interval (a, b, 0.001, 0.0);

if (status == GSL_SUCCESS)
printf ("Converged:\n");
printf ("%5d [%.7f, %.7f] " "%.7f %+.7f %.7f\n",  iter, a, b, m, m - m_expected, b - a);}

while (status == GSL_CONTINUE && iter < max_iter);
gsl_min_fminimizer_free (s);

return status;
}

Last edited:
Mark44
Mentor
Together with the other threads, this is getting frustrating. The actual problem that needs to be solved has never been explicitly described in details, so we don't know what @CAF123 is actually after.
Amen.
Hopefully the issue about terminology is clarified now, and it's agreed that gsl routine linked in OP can definitely do what I want.
Since, as @DrClaude wrote, we don't understand what the actual problem is, I don't think we can dismiss the terminology issue, nor is it clear to me that the gsl routine will do what you want.
some texts are sloppy because I see things like 'convolution integral'
"Convolution" and "convolution integral" are synonymous.
In any case, what I really want is to maximise the product function f(x)g(x).
Then why did you write an integral in post #1 of this thread? If all you want to do is to maximize or minimize the product of two functions, that's a completely different situation than xxximizing an integral whose integrand is a product of two functions.
Questions:
1. Are you trying to maximize something or minimize it? The thread title mentions maximizing. Your code sample suggests that you want to minimize something.
2. What exactly are you trying to maximize (or minimize, since this isn't clear)? At various points in this and other threads, you seem to be asking about maximizing/minimizing 1) a function (e.g. f(x)g(x)), 2) a definite integral that is a number, or 3) a definite integral with a variable as a limit of integration (e.g., of the form ##g(x) = \int_a^x f(t)dt##). Choice 2 is not possible, because nothing is varying, so it can't be minimized/maximized.
3. You've mentioned the term "convolution" in several places. Unless there's a wildly different meaning for this in physics as compared to the mathematical definition (see post #7), the convolution of two functions is NOT a product - it's a definite integral that has different values, depending on the parameter, and whose integrand involves a product of two functions -- again see post #7 as well as https://en.wikipedia.org/wiki/Convolution.

DrClaude
PeterDonis
Mentor
2020 Award
I had a private discussion with PeterDonis and he told me to create a new thread.

Yes, and here is what you said you wanted in that private discussion:

I want the value of x=x_max that maximises f(x)g(x) over the interval [a,b] and to give the value of f(x_max)g(x_max).

And then, after I explained the simple way of finding the answer to your question, you said:

I want it to do this in C++ which I still regard myself as an amateur in so I posted it in the programming forum. The expression f(x)g(x) is sufficiently long that by hand it is not efficient. I think gsl numeric differentiation can do this.

But in the OP of this thread you are saying something different. So I am closing this thread because you are still not being either clear or consistent about what you want. And in any case, from the responses you have gotten, it does not seem like PF is the right place to find the answer to whatever it is you want; it sounds like you either need to read the documentation for the software you are using, or talk to whoever gave you this assignment, or both.

Mark44
Mentor
it sounds like you either need to read the documentation for the software you are using, or talk to whoever gave you this assignment, or both.
To which I would add, take a class or online training in C++, since even just using the GSL API requires considerable sophistication in the use of function pointers and other arcana.

jedishrfu
Mentor
Reopening for the moment in order for some clarifications to be posted.

Gold Member
Thanks. As in the OP, my main goal is the following:

I'd like to use https://www.gnu.org/software/gsl/doc/html/min.html to find the maximum of ##f(x)g(x)## in the domain ##x \in [-1,1]##.

I've set up the function but having a little trouble in referencing the parameters correctly. See my post #11. In the example within the link, the function is cos x + 1 without any additional parameters so the choice F.params=0; in that case.

DrClaude
Mentor
1) What if, in some situation, one doesn't know the expected maximum of the function and the value of x at which that maximum occurs?
The vast majority of times, this is exactly the case. One doesn't know the max of a function, and the numerical method is there to help find it.

2) Could you clarify what is meaning of
Code:
 F.function = &minimum_func;
? Without looking at the gsl header files, I guess it's just a way of accessing the relevant function.
The GSL function needs to call your function, so this sets a pointer to it so GSL can internally call it.

I declare the function via a double and then put the gsl method in main function, as in example code in link. I get an error that, by referencing the parameters via
Code:
 F.params = (void*) (&func_params)
, that 'Use of unidentified identifier 'func_params'.
You haven't declared the variable func_params, it doesn't exist. If you need to pass a parameter to minimum_func, you must declare a variable of the right type (int, double, some struct, ...) and then set F.params to point to it. If your function doesn't need parameters, you can simply set
C:
 F.params = NULL;

Gold Member
The vast majority of times, this is exactly the case. One doesn't know the max of a function, and the numerical method is there to help find it.
Thanks, in that case should one enter a somewhat arbitrary number for m (within domain) and m_expected within the gsl method?

You haven't declared the variable func_params, it doesn't exist. If you need to pass a parameter to minimum_func, you must declare a variable of the right type (int, double, some struct, ...) and then set F.params to point to it.
I thought it had been declared in the definition of double minimum_func. With reference to a struct, I tried F.params = (void*)(&func_params_t) but it said func_params_t does not refer to a value.

DrClaude
Mentor
Thanks, in that case should one enter a somewhat arbitrary number for m (within domain) and m_expected within the gsl method?
The minimization routine will try to find a local minimum near the initial guess m. If the function is monotonic in the specified range, then the value of m doesn't matter much. Plotting the function can be helpful in figuring out what initial value to use.

There is no m_expected in the method. This is only used in the example in the GSL documentation to check the iterations and final accuracy, since it is tested on a function where the minimum is known analytically.

I thought it had been declared in the definition of double minimum_func. With reference to a struct, I tried F.params = (void*)(&func_params_t) but it said func_params_t does not refer to a value.
The operator & means "address of." It can be applied to a variable, not to a datatype. You need to declare a variable of type func_params_t and use it in setting F.params. You should also set the proper parameters in that variable before calling gsl_min_fminimizer_set.

CAF123
Gold Member
Plotting the function can be helpful in figuring out what initial value to use.

Sorry for late reply. I would like to do the plot you suggested above. For this, probably a for loop is sufficient, something like
C:
for(int i = 0; i <=100; i++) {
z = (1-(-1))*double(i)/100 + -1;
std::cout << "z:" << z << "prod:" << minimum_function(z,void*p) << std::endl;
}

This would print out a table of z values and the corresponding values of minimum_function. There is an error about 'Expected '(' for function-style cast'. I guess this is because I am trying to print out the pointer rather than what the pointer is pointing at. I know the members of the struct are accessed via the "." operator but in terms of the code in the for loop above, what should I write in place of void*p to correctly print out minimum_function for a given z and other values of the struct?

Last edited by a moderator:
DrClaude
Mentor
If p is a pointer, simply use (z,p) as the arguments. If p is not a pointer, use (z,&p).

pbuk
Gold Member
There is an error about 'Expected '(' for function-style cast'. I guess this is because I am trying to print out the pointer rather than what the pointer is pointing at.
No, it's because
C:
void*p
is not valid syntax here. You are confusing the way a function is declared as in your post #5
C:
double minimum_function (double x, void *params)
with the way a function is called, which for the declaration above can be any of:
C:
minimum_function(a) // with only one argument of type double
minimum_function(a, b) // with two arguments where b is a pointer
minimum_function(a, &c) // with the second argument a pointer to the variable c
// or even
minimum_function((double) d , &c) // where d is of, say, integer type and is being cast to a double value

DrClaude