Constellation constrained channel capacity with additive white Gaussian noise

Computation of the channel capacity of a given constellation with AWGN

Consider the following channel model,
Y = sqrt(snr)X + N,
where X is the (discrete) channel input such as E[X^2]<1, N is the noise with density function N(0,1) and snr the signal-to-noise ratio (linear scale). The function computes the capacity associated to any constellation in one or two dimension(s) with a given probability density function for the channel inputs. It enables to evaluate the capacity of any possible input distributions, especially classical modulation schemes (QAM, PSK, etc) with uniform inputs.
I rely on the quad2d function but it is also possible to use Monte Carlo methods to evaluate the integral that appears in the mutual information computation.