draw_sample() function for linear operators seems to lack a precise definition.
Whereas in the case of a curvature C the
draw_sample function draws a sample from the corresponding Gaussian approximation of the underlying probability distribution with said curvature (i.e. exp(-0.5 x^+ Cx) ), in the case of a
DiagonalOperator D it seems to use the distribution wrt the inverse of it (i.e. exp(-0.5 x^+ D^-1 x) )
I propose changing the latter to match the former for consistency.
The curvature definition is more likely to be the one that someone actually wants (at least for the curvature case. For the diagonal operator case one could start arguing, but ... consistency!). Also the
draw_sample function for the
DiagonalOperator class is not that hard to change.