merge

.gitlab-ci.yml

@@ -65,11 +65,16 @@ pages:
 before_script:
   - python3 setup.py install --user -f
 
-run_ipynb:
+run_ipynb0:
   stage: demo_runs
   script:
     - jupyter nbconvert --execute --ExecutePreprocessor.timeout=None demos/getting_started_0.ipynb
 
+run_ipynb1:
+  stage: demo_runs
+  script:
+    - jupyter nbconvert --execute --ExecutePreprocessor.timeout=None demos/getting_started_4_CorrelatedFields.ipynb
+
 run_getting_started_1:
   stage: demo_runs
   script:
@@ -99,7 +104,7 @@ run_getting_started_3:
 run_getting_started_mf:
   stage: demo_runs
   script:
-    - python3 demos/getting_started_mf.py
+    - python3 demos/getting_started_5_mf.py
   artifacts:
     paths:
       - '*.png'

demos/getting_started_3.py

@@ -55,6 +55,10 @@ if __name__ == '__main__':
 
     position_space = ift.RGSpace([128, 128])
 
+    #  For a detailed showcase of the effects the parameters
+    #  of the CorrelatedField model have on the generated fields,
+    #  see 'getting_started_4_CorrelatedFields.ipynb'.
+
     cfmaker = ift.CorrelatedFieldMaker.make(
         offset_mean =      0.0,  # 0.
         offset_std_mean = 1e-3,  # 1e-3
@@ -62,22 +66,21 @@ if __name__ == '__main__':
         prefix = '')
 
     fluctuations_dict = {
-        # Amplitude of the fluctuations
+        # Amplitude of field fluctuations
         'fluctuations_mean':   2.0,  # 1.0
         'fluctuations_stddev': 1.0,  # 1e-2
 
-        # Smooth variation speed
+        # Exponent of power law power spectrum component
+        'loglogavgslope_mean': -2.0,  # -3.0
+        'loglogavgslope_stddev': 0.5,  #  0.5
+
+        # Amplitude of integrated Wiener process power spectrum component
         'flexibility_mean':   2.5,  # 1.0
         'flexibility_stddev': 1.0,  # 0.5
 
-        # How strong the ragged component of the spectrum is
-        # (Ratio of Wiener process and integrated Wiener process ?)
+        # How ragged the integrated Wiener process component is
         'asperity_mean':   0.5,  # 0.1
-        'asperity_stddev': 0.5,  # 0.5
-
-        # Slope of linear spectrum component
-        'loglogavgslope_mean': -2.0,  # -3.0
-        'loglogavgslope_stddev': 0.5  #  0.5
+        'asperity_stddev': 0.5  # 0.5
     }
     cfmaker.add_fluctuations(position_space, **fluctuations_dict)
 
@@ -135,7 +138,7 @@ if __name__ == '__main__':
         # Plot current reconstruction
         plot = ift.Plot()
         plot.add(signal(KL.position), title="reconstruction")
-        plot.add([A.force(KL.position), A.force(mock_position)], title="power")
+        plot.add([A.force(mock_position), A.force(KL.position)], title="power")
         plot.output(ny=1, ysize=6, xsize=16,
                     name=filename.format("loop_{:02d}".format(i)))
 
@@ -153,8 +156,8 @@ if __name__ == '__main__':
 
     powers = [A.force(s + KL.position) for s in KL.samples]
     plot.add(
-        powers + [A.force(KL.position),
-                  A.force(mock_position)],
+        powers + [A.force(mock_position),
+                  A.force(KL.position)],
         title="Sampled Posterior Power Spectrum",
         linewidth=[1.]*len(powers) + [3., 3.])
     plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename_res)

demos/getting_started_4_CorrelatedFields.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Notebook showcasing the NIFTy 6 Correlated Field model\n",
+    "\n",
+    "**Skip to `Parameter Showcases` for the meat/veggies ;)**\n",
+    "\n",
+    "The field model roughly works like this:\n",
+    "\n",
+    "`f = HT( A * zero_mode * xi ) + offset`\n",
+    "\n",
+    "`A` is a spectral power field which is constructed from power spectra that hold on subdomains of the target domain.\n",
+    "It is scaled by a zero mode operator and then pointwise multiplied by a gaussian excitation field, yielding\n",
+    "a representation of the field in harmonic space.\n",
+    "It is then transformed into the target real space and a offset added.\n",
+    "\n",
+    "The power spectra `A` is constructed of are in turn constructed as the sum of a power law component\n",
+    "and an integrated Wiener process whose amplitude and roughness can be set.\n",
+    "\n",
+    "## Setup code"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import nifty7 as ift\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "ift.random.push_sseq_from_seed(43)\n",
+    "\n",
+    "n_pix = 256\n",
+    "x_space = ift.RGSpace(n_pix)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plotting routine\n",
+    "def plot(fields, spectra, title=None):\n",
+    "    # Plotting preparation is normally handled by nifty7.Plot\n",
+    "    # It is done manually here to be able to tweak details\n",
+    "    # Fields are assumed to have identical domains\n",
+    "    fig = plt.figure(tight_layout=True, figsize=(12, 3.5))\n",
+    "    if title is not None:\n",
+    "        fig.suptitle(title, fontsize=14)\n",
+    "    \n",
+    "    # Field\n",
+    "    ax1 = fig.add_subplot(1, 2, 1)\n",
+    "    ax1.axhline(y=0., color='k', linestyle='--', alpha=0.25)\n",
+    "    for field in fields:\n",
+    "        dom = field.domain[0]\n",
+    "        xcoord = np.arange(dom.shape[0]) * dom.distances[0]\n",
+    "        ax1.plot(xcoord, field.val)\n",
+    "    ax1.set_xlim(xcoord[0], xcoord[-1])\n",
+    "    ax1.set_ylim(-5., 5.)\n",
+    "    ax1.set_xlabel('x')\n",
+    "    ax1.set_ylabel('f(x)')\n",
+    "    ax1.set_title('Field realizations')\n",
+    "    \n",
+    "    # Spectrum\n",
+    "    ax2 = fig.add_subplot(1, 2, 2)\n",
+    "    for spectrum in spectra:\n",
+    "        xcoord = spectrum.domain[0].k_lengths\n",
+    "        ycoord = spectrum.val_rw()\n",
+    "        ycoord[0] = ycoord[1]\n",
+    "        ax2.plot(xcoord, ycoord)\n",
+    "    ax2.set_ylim(1e-6, 10.)\n",
+    "    ax2.set_xscale('log')\n",
+    "    ax2.set_yscale('log')\n",
+    "    ax2.set_xlabel('k')\n",
+    "    ax2.set_ylabel('p(k)')\n",
+    "    ax2.set_title('Power Spectrum')\n",
+    "    \n",
+    "    fig.align_labels()\n",
+    "    plt.show()\n",
+    "\n",
+    "# Helper: draw main sample\n",
+    "main_sample = None\n",
+    "def init_model(m_pars, fl_pars):\n",
+    "    global main_sample\n",
+    "    cf = ift.CorrelatedFieldMaker.make(**m_pars)\n",
+    "    cf.add_fluctuations(**fl_pars)\n",
+    "    field = cf.finalize(prior_info=0)\n",
+    "    main_sample = ift.from_random(field.domain)\n",
+    "    print(\"model domain keys:\", field.domain.keys())\n",
+    "    \n",
+    "# Helper: field and spectrum from parameter dictionaries + plotting\n",
+    "def eval_model(m_pars, fl_pars, title=None, samples=None):\n",
+    "    cf = ift.CorrelatedFieldMaker.make(**m_pars)\n",
+    "    cf.add_fluctuations(**fl_pars)\n",
+    "    field = cf.finalize(prior_info=0)\n",
+    "    spectrum = cf.amplitude\n",
+    "    if samples is None:\n",
+    "        samples = [main_sample]\n",
+    "    field_realizations = [field(s) for s in samples]\n",
+    "    spectrum_realizations = [spectrum.force(s) for s in samples]\n",
+    "    plot(field_realizations, spectrum_realizations, title)\n",
+    "\n",
+    "def gen_samples(key_to_vary):\n",
+    "    if key_to_vary is None:\n",
+    "        return [main_sample]\n",
+    "    dct = main_sample.to_dict()\n",
+    "    subdom_to_vary = dct.pop(key_to_vary).domain\n",
+    "    samples = []\n",
+    "    for i in range(8):\n",
+    "        d = dct.copy()\n",
+    "        d[key_to_vary] = ift.from_random(subdom_to_vary)\n",
+    "        samples.append(ift.MultiField.from_dict(d))\n",
+    "    return samples\n",
+    "        \n",
+    "def vary_parameter(parameter_key, values, samples_vary_in=None):\n",
+    "    s = gen_samples(samples_vary_in)\n",
+    "    for v in values:\n",
+    "        if parameter_key in cf_make_pars.keys():\n",
+    "            m_pars = {**cf_make_pars, parameter_key: v}\n",
+    "            eval_model(m_pars, cf_x_fluct_pars, f\"{parameter_key} = {v}\", s)\n",
+    "        else:\n",
+    "            fl_pars = {**cf_x_fluct_pars, parameter_key: v}\n",
+    "            eval_model(cf_make_pars, fl_pars, f\"{parameter_key} = {v}\", s)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Before the Action: The Moment-Matched Log-Normal Distribution\n",
+    "Many properties of the correlated field are modelled as being lognormally distributed.\n",
+    "\n",
+    "The distribution models are parametrized via their means `_mean` and standard-deviations `_stddev`.\n",
+    "\n",
+    "To get a feeling of how the ratio of the `mean` and `stddev` parameters influences the distribution shape,\n",
+    "here are a few example histograms: (observe the x-axis!)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = plt.figure(figsize=(13, 3.5))\n",
+    "mean = 1.0\n",
+    "sigmas = [1.0, 0.5, 0.1]\n",
+    "\n",
+    "for i in range(3):\n",
+    "    op = ift.library.correlated_fields._LognormalMomentMatching(mean=mean,\n",
+    "                                                                sig=sigmas[i],\n",
+    "                                                                key='foo',\n",
+    "                                                                N_copies=0)\n",
+    "    op_samples = np.array(\n",
+    "        [op(s).val for s in [ift.from_random(op.domain) for i in range(10000)]])\n",
+    "\n",
+    "    ax = fig.add_subplot(1, 3, i + 1)\n",
+    "    ax.hist(op_samples, bins=50)\n",
+    "    ax.set_title(f\"mean = {mean}, sigma = {sigmas[i]}\")\n",
+    "    ax.set_xlabel('x')\n",
+    "    del op_samples\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The Neutral Field\n",
+    "\n",
+    "To demonstrate the effect of all parameters, first a 'neutral' set of parameters\n",
+    "is defined which then are varied one by one, showing the effect of the variation\n",
+    "on the generated field realizations and the underlying power spectrum from which\n",
+    "they were drawn.\n",
+    "\n",
+    "As a neutral field, a model with a white power spectrum and vanishing spectral power was chosen."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Neutral model parameters yielding a quasi-constant field\n",
+    "cf_make_pars = {\n",
+    "    'offset_mean': 0.,\n",
+    "    'offset_std_mean': 1e-3,\n",
+    "    'offset_std_std': 1e-16,\n",
+    "    'prefix': ''\n",
+    "}\n",
+    "\n",
+    "cf_x_fluct_pars = {\n",
+    "    'target_subdomain': x_space,\n",
+    "    'fluctuations_mean': 1e-3,\n",
+    "    'fluctuations_stddev': 1e-16,\n",
+    "    'flexibility_mean': 1e-3,\n",
+    "    'flexibility_stddev': 1e-16,\n",
+    "    'asperity_mean': 1e-3,\n",
+    "    'asperity_stddev': 1e-16,\n",
+    "    'loglogavgslope_mean': 0.,\n",
+    "    'loglogavgslope_stddev': 1e-16\n",
+    "}\n",
+    "\n",
+    "init_model(cf_make_pars, cf_x_fluct_pars)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Show neutral field\n",
+    "eval_model(cf_make_pars, cf_x_fluct_pars, \"Neutral Field\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Parameter Showcases\n",
+    "\n",
+    "## The `fluctuations_` parameters of `add_fluctuations()`\n",
+    "\n",
+    "determine the **amplitude of variations along the field dimension**\n",
+    "for which `add_fluctuations` is called.\n",
+    "\n",
+    "`fluctuations_mean` set the _average_ amplitude of the fields fluctuations along the given dimension,\\\n",
+    "`fluctuations_stddev` sets the width and shape of the amplitude distribution.\n",
+    "\n",
+    "The amplitude is modelled as being log-normally distributed,\n",
+    "see `The Moment-Matched Log-Normal Distribution` above for details.\n",
+    "\n",
+    "#### `fluctuations_mean`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vary_parameter('fluctuations_mean', [0.05, 0.5, 2.], samples_vary_in='xi')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### `fluctuations_stddev`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cf_x_fluct_pars['fluctuations_mean'] = 1.0\n",
+    "vary_parameter('fluctuations_stddev', [0.01, 0.1, 1.0], samples_vary_in='fluctuations')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The `loglogavgslope_` parameters of `add_fluctuations()`\n",
+    "\n",
+    "determine **the slope of the loglog-linear (power law) component of the power spectrum**.\n",
+    "\n",
+    "The slope is modelled to be normally distributed.\n",
+    "\n",
+    "#### `loglogavgslope_mean`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vary_parameter('loglogavgslope_mean', [-3., -1., 1.], samples_vary_in='xi')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### `loglogavgslope_stddev`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cf_x_fluct_pars['loglogavgslope_mean'] = -1.0\n",
+    "vary_parameter('loglogavgslope_stddev', [0.01, 0.1, 1.0], samples_vary_in='loglogavgslope')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The `flexibility_` parameters of `add_fluctuations()`\n",
+    "\n",
+    "determine **the amplitude of the integrated Wiener process component of the power spectrum**\n",
+    "(how strong the power spectrum varies besides the power-law).\n",
+    "\n",
+    "`flexibility_mean` sets the _average_ amplitude of the i.g.p. component,\n",
+    "`flexibility_stddev` sets how much the amplitude can vary.\\\n",
+    "These two parameters feed into a moment-matched log-normal distribution model,\n",
+    "see above for a demo of its behavior.\n",
+    "\n",
+    "#### `flexibility_mean`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vary_parameter('flexibility_mean', [0.1, 1.0, 3.0], samples_vary_in='spectrum')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### `flexibility_stddev`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cf_x_fluct_pars['flexibility_mean'] = 2.0\n",
+    "vary_parameter('flexibility_stddev', [0.01, 0.1, 1.0], samples_vary_in='flexibility')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The `asperity_` parameters of `add_fluctuations()`\n",
+    "\n",
+    "`asperity_` determines **how rough the integrated Wiener process component of the power spectrum is**.\n",
+    "\n",
+    "`asperity_mean` sets the average roughness, `asperity_stddev` sets how much the roughness can vary.\\\n",
+    "These two parameters feed into a moment-matched log-normal distribution model,\n",
+    "see above for a demo of its behavior.\n",
+    "\n",
+    "#### `asperity_mean`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vary_parameter('asperity_mean', [0.001, 1.0, 5.0], samples_vary_in='spectrum')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### `asperity_stddev`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cf_x_fluct_pars['asperity_mean'] = 1.0\n",
+    "vary_parameter('asperity_stddev', [0.01, 0.1, 1.0], samples_vary_in='asperity')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The `offset_mean` parameter of `CorrelatedFieldMaker.make()`\n",
+    "\n",
+    "The `offset_mean` parameter defines a global additive offset on the field realizations.\n",
+    "\n",
+    "If the field is used for a lognormal model `f = field.exp()`, this acts as a global signal magnitude offset."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Reset model to neutral\n",
+    "cf_x_fluct_pars['fluctuations_mean'] = 1e-3\n",
+    "cf_x_fluct_pars['flexibility_mean'] = 1e-3\n",
+    "cf_x_fluct_pars['asperity_mean'] = 1e-3\n",
+    "cf_x_fluct_pars['loglogavgslope_mean'] = 1e-3"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vary_parameter('offset_mean', [3., 0., -2.])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The `offset_std_` parameters of `CorrelatedFieldMaker.make()`\n",
+    "\n",
+    "Variation of the global offset of the field are modelled as being log-normally distributed.\n",
+    "See `The Moment-Matched Log-Normal Distribution` above for details.\n",
+    "\n",
+    "The `offset_std_mean` parameter sets how much NIFTy will vary the offset *on average*.\\\n",
+    "The `offset_std_std` parameters defines the with and shape of the offset variation distribution.\n",
+    "\n",
+    "#### `offset_std_mean`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vary_parameter('offset_std_mean', [1e-16, 0.5, 2.], samples_vary_in='xi')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### `offset_std_std`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cf_make_pars['offset_std_mean'] = 1.0\n",
+    "vary_parameter('offset_std_std', [0.01, 0.1, 1.0], samples_vary_in='zeromode')"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

src/library/correlated_fields.py

@@ -356,24 +356,23 @@ class CorrelatedFieldMaker:
     """Constrution helper for hirarchical correlated field models.
 
     The correlated field models are parametrized by creating
-    power spectrum operators ("amplitudes")
-    acting on their target subdomains
-    via calls to :func:`add_fluctuations`.
+    power spectrum operators ("amplitudes") via calls to
+    :func:`add_fluctuations` that act on the targeted field subdomains.
     During creation of the :class:`CorrelatedFieldMaker` via
-    :func:`make`, a global offset from zero can be added to the
-    field to be created and an operator applying gaussian fluctuations
-    around this offset needs to be parametrized.
+    :func:`make`, a global offset from zero of the field model
+    can be defined and an operator applying fluctuations
+    around this offset is parametrized.
 
     The resulting correlated field model operator has a
     :class:`~nifty7.multi_domain.MultiDomain` as its domain and
     expects its input values to be univariately gaussian.
 
     The target of the constructed operator will be a
-    :class:`~nifty7.domain_tuple.DomainTuple`
-    containing the `target_subdomains` of the added fluctuations in the
-    order of the `add_fluctuations` calls.
+    :class:`~nifty7.domain_tuple.DomainTuple` containing the
+    `target_subdomains` of the added fluctuations in the order of
+    the `add_fluctuations` calls.
 
-    Creation of the model operator is finished by calling the method
+    Creation of the model operator is completed by calling the method
     :func:`finalize`, which returns the configured operator.
 
     An operator representing an array of correlated field models
@@ -482,9 +481,9 @@ class CorrelatedFieldMaker:
         fluctuations_{mean,stddev} : float
             Total spectral energy -> Amplitude of the fluctuations