  We present a full forward-modeled $w$CDM analysis of the KiDS-1000 weak
lensing maps using graph-convolutional neural networks (GCNN). Utilizing the
$\texttt{CosmoGrid}$, a novel massive simulation suite spanning six different
cosmological parameters, we generate almost one million tomographic mock
surveys on the sphere. Due to the large data set size and survey area, we
perform a spherical analysis while limiting our map resolution to
$\texttt{HEALPix}$ $n_\mathrm{side}=512$. We marginalize over systematics such
as photometric redshift errors, multiplicative calibration and additive shear
bias. Furthermore, we use a map-level implementation of the non-linear
intrinsic alignment model along with a novel treatment of baryonic feedback to
incorporate additional astrophysical nuisance parameters. We also perform a
spherical power spectrum analysis for comparison. The constraints of the
cosmological parameters are generated using a likelihood free inference method
called Gaussian Process Approximate Bayesian Computation (GPABC). Finally, we
check that our pipeline is robust against choices of the simulation parameters.
We find constraints on the degeneracy parameter of $S_8 \equiv
\sigma_8\sqrt{\Omega_M/0.3} = 0.78^{+0.06}_{-0.06}$ for our power spectrum
analysis and $S_8 = 0.79^{+0.05}_{-0.05}$ for our GCNN analysis, improving the
former by 16%. This is consistent with earlier analyses of the 2-point
function, albeit slightly higher. Baryonic corrections generally broaden the
constraints on the degeneracy parameter by about 10%. These results offer great
prospects for full machine learning based analyses of on-going and future weak
lensing surveys.
