Maximising a convolution in C++ via a GSL routine

In summary: To use it, you first need to alloc it using malloc, and then you can access its members using the & operator. In this particular line, the & operator is used to point to the func_params_t structure, and then the = operator is used to assign func_params_t's members to the params variable.
  • #1
CAF123
Gold Member
2,948
88
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?
 
Technology news on Phys.org
  • #2
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
CAF123 said:
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);
}
 
  • #3
DrClaude said:
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?
 
  • #4
CAF123 said:
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.
 
  • Like
Likes CAF123
  • #5
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:
  • #6
CAF123 said:
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.

CAF123 said:
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.
 
  • #7
CAF123 said:
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.
 
  • #8
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:
  • #9
Further thoughts.
CAF123 said:
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:
  • #10
Mark44 said:
In contrast, the function you want to maximize is the convolution of two functions
I'm not sure about that, see
CAF123 said:
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.
 
  • #11
DrClaude said:
Is that line correct? I don't see how this could be valid code.
I rechecked and it appears to be correct.
Mark44 said:
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).
Mark44 said:
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.

DrClaude said:
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:
  • #12
DrClaude said:
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.
CAF123 said:
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.
CAF123 said:
some texts are sloppy because I see things like 'convolution integral'
"Convolution" and "convolution integral" are synonymous.
CAF123 said:
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.
 
  • Like
Likes DrClaude
  • #13
CAF123 said:
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:

CAF123 said:
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:

CAF123 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.
 
  • #14
PeterDonis said:
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.
 
  • #15
Reopening for the moment in order for some clarifications to be posted.

Thanks for your patience.
 
  • #16
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.
 
  • #17
CAF123 said:
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.

CAF123 said:
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.

CAF123 said:
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;
 
  • #18
DrClaude said:
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.
 
  • #19
CAF123 said:
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.

CAF123 said:
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.
 
  • Like
Likes CAF123
  • #20
DrClaude said:
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:
  • #21
If p is a pointer, simply use (z,p) as the arguments. If p is not a pointer, use (z,&p).
 
  • #22
CAF123 said:
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
 
  • Like
Likes DrClaude

1. How can I use the GSL routine to maximize a convolution in C++?

To use the GSL (GNU Scientific Library) routine for maximizing a convolution in C++, you will first need to include the appropriate header file (gsl/gsl_multimin.h). Then, you can use the gsl_multimin_fminimizer_set() function to set up the minimizer and the gsl_multimin_fminimizer_minimize() function to perform the minimization.

2. What is a convolution and why is it important to maximize it?

A convolution is a mathematical operation that combines two functions to create a third function. In the context of signal processing and data analysis, convolutions are frequently used to smooth or filter data. Maximizing a convolution can help improve the accuracy and reliability of data analysis results.

3. Are there any alternatives to using the GSL routine for maximizing a convolution in C++?

Yes, there are other libraries and routines available for maximizing a convolution in C++, such as the Numerical Recipes library and the Boost C++ libraries. However, the GSL routine is a popular and well-documented option that is specifically designed for scientific computing.

4. Can the GSL routine be used for multidimensional convolutions?

Yes, the GSL routine supports multidimensional convolutions. You can specify the dimensionality of the convolution in the gsl_multimin_fminimizer_set() function by setting the size of the initial guess vector and the function to be minimized.

5. How can I test the performance of the GSL routine for maximizing a convolution in C++?

The GSL library includes a set of example programs that demonstrate the usage and performance of the routines, including the gsl_multimin_fminimizer. You can also write your own test code to compare the results of the GSL routine with other methods or to test its performance on different input data sets.

Similar threads

  • Programming and Computer Science
Replies
7
Views
1K
  • Programming and Computer Science
Replies
11
Views
947
  • Programming and Computer Science
Replies
1
Views
2K
  • Programming and Computer Science
Replies
19
Views
2K
  • Programming and Computer Science
Replies
1
Views
662
  • Programming and Computer Science
Replies
17
Views
1K
  • Calculus and Beyond Homework Help
Replies
10
Views
1K
  • Programming and Computer Science
Replies
3
Views
1K
  • Programming and Computer Science
Replies
1
Views
1K
  • Programming and Computer Science
Replies
1
Views
613
Back
Top