BIC for generalized additive models Is there any way to use BIC in model selection for gam? And if so then how to extract the penalty in the bic formula?

As long as the different observations $\{y_i, \mathbf{x}_i\}_{1 \leq i \leq T}$ that you use to fit your GAM are i.i.d. (and assuming that the number of observations $T$ is reasonably large) it is perfectly possible to compute its BIC.
Let's assume that your model has the following form:
$$
g(\mathbf{E}(y)) = \alpha + f_{\theta_1}^{(1)}(x^{(1)}) + f_{\theta_2}^{(2)}(x^{(2)})+ f_{\theta_3}^{(3)}(x^{(3)}) + \dots
$$
where $\mathbf{x} = [x^{(1)},\dots,x^{(N)}]$, and where $\theta_1,\dots,\theta_N$ are the free parameters of the functions you are fitting on your data. Then the total number of free parameters in your GAM will be
$$
1+ \sum_{j=1}^N |\theta_j|
$$
i.e. the number of free parameters in each function, plus one (for $\alpha$).
