  I present a new approach to recover the primordial density fluctuations and
the cosmic web structure underlying a galaxy distribution. The method is based
on sampling Gaussian fields which are compatible with a galaxy distribution and
a structure formation model. This is achieved by splitting the inversion
problem into two Gibbs-sampling steps: the first being a Gaussianisation step
transforming a distribution of point sources at Lagrangian positions -which are
not a priori given- into a linear alias-free Gaussian field. This step is based
on Hamiltonian sampling with a Gaussian-Poisson model. The second step consists
on a likelihood comparison in which the set of matter tracers at the initial
conditions is constrained on the galaxy distribution and the assumed structure
formation model. For computational reasons second order Lagrangian Perturbation
Theory is used. However, the presented approach is flexible to adopt any
structure formation model. A semi-analytic halo-model based galaxy mock catalog
is taken to demonstrate that the recovered initial conditions are closely
unbiased with respect to the actual ones from the corresponding N-body
simulation down to scales of a ~ 5 Mpc/h. The cross-correlation between them
shows a substantial gain of information, being at k ~ 0.3 h/Mpc more than
doubled. In addition the initial conditions are extremely well Gaussian
distributed and the power-spectra follow the shape of the linear power-spectrum
being very close to the actual one from the simulation down to scales of k ~ 1
h/Mpc.
