I will think about it again.  However, Since we can have different variational bounds (classes) for one likelihood (for example, logit), C++ users will use these classes when they want to use variational inference for GP.   Another way to do implement is to make the specific likelihood as a data member of the variational bound class and public some methods in the likelihood.  Do you have any suggestion? 