tickets/DM-28474: Update to latest master

docs/0-quickstart.ipynb

@@ -47,9 +47,20 @@
     "data = np.load(\"../data/hsc_cosmos_35.npz\")\n",
     "images = data[\"images\"]\n",
     "filters = data[\"filters\"]\n",
-    "catalog = data[\"catalog\"]\n",
+    "centers = [(src['y'], src['x']) for src in data[\"catalog\"]] # Note: y/x convention!\n",
     "weights = 1/data[\"variance\"]\n",
-    "psfs = scarlet.ImagePSF(data[\"psfs\"])"
+    "psf = scarlet.ImagePSF(data[\"psfs\"])"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    ".. warning::\n",
+    "\n",
+    "    Coordinates in *scarlet* are given in the C/numpy notation (y,x) as opposed to the more conventional mathematical (x,y) ordering."
    ]
   },
   {
@@ -97,7 +108,7 @@
    "source": [
     ".. note::\n",
     "\n",
-    "    It is crictical to realize that there are two frames: one for model that reconstructs the scene on the sky, and one for any observation of the scene. They may be identical, but mostly they are not, in which case the observation is an information-reduced rendering of the model. The loss of information could be in terms of signal-to-noise ratio, spatial resolution, PSF blurring, spectral coverage, or all of the above."
+    "    It is critical to realize that there are two frames: one for model that represents the scene on the sky, and one for any observation of that scene. The frames may be identical, but usually they are not, in which case the observation is compared to a degraded rendering of the model. The degradation could be in terms of spatial resolution, PSF blurring, spectral coverage, or all of the above. If the model frame and the observation frame are the same, no attempt can be made to recover information lost by the degradation."
    ]
   },
   {
@@ -131,12 +142,12 @@
    "source": [
     "model_frame = scarlet.Frame(\n",
     "    images.shape,\n",
-    "    psfs=model_psf,\n",
+    "    psf=model_psf,\n",
     "    channels=filters)\n",
     "\n",
     "observation = scarlet.Observation(\n",
     "    images, \n",
-    "    psfs=psfs, \n",
+    "    psf=psf,\n",
     "    weights=weights, \n",
     "    channels=filters).match(model_frame)"
    ]
@@ -168,8 +179,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sky_coords = [[src['y'], src['x']] for src in catalog] # y/x!\n",
-    "scarlet.display.show_observation(observation, norm=norm, sky_coords=sky_coords, show_psf=True)\n",
+    "scarlet.display.show_observation(observation, norm=norm, sky_coords=centers, show_psf=True)\n",
     "plt.show()"
    ]
   },
@@ -180,17 +190,6 @@
     "Since we use a trivial `wcs` in this `Observation`, all coordinates are already in image pixels, otherwise RA/Dec pairs are expected as sky coordinates. Also:"
    ]
   },
-  {
-   "cell_type": "raw",
-   "metadata": {
-    "raw_mimetype": "text/restructuredtext"
-   },
-   "source": [
-    ".. warning::\n",
-    "\n",
-    "    Coordinates in *scarlet* are given in the C/numpy notation (y,x) as opposed to the more conventional mathematical (x,y) ordering."
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -201,10 +200,40 @@
     "\n",
     "* `RandomSource` fits per-band amplitude and non-parametric morphology starting from uniform random draws for both.\n",
     "* `PointSource` fits centers and per-band amplitude using the observed PSF model.\n",
-    "* `ExtendedSource` fits per-band amplitude and a non-parametric morphology (which is constrained to be monotonically decreasing from the center).\n",
+    "* `ExtendedSource` fits per-band amplitude and a non-parametric morphology (which is constrained to be monotonically decreasing from the center). It has a `compact` option, which initializes the source with the shape of a point source.\n",
     "* `ExtendedSource(K=2)` splits an `ExtendedSource` into `K` components that are initially radially separated and stacked on top of each other like a pyramid.\n",
     "\n",
-    "In our example, we assume *prior* knowledge that object 0 is a star, and object 1 should be modeled with two components. Everything else is assumed a single-component galaxy. **We generally recommend `ExtendedSource` as default** if additional information about the source is not available."
+    "We can use a robust initialization scheme, that starts by attempting to model every source with 2 components. If these components don't have enough signal-to-noise, it will fall back to 1 component. If that fails, it will work with a `compact` component. If cannot do even that, it will report the source number in `skipped`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sources, skipped = scarlet.initialization.init_all_sources(model_frame, \n",
+    "                                                           centers, \n",
+    "                                                           observation, \n",
+    "                                                           max_components=2, \n",
+    "                                                           min_snr=50,\n",
+    "                                                           thresh=1,\n",
+    "                                                           fallback=True,\n",
+    "                                                           silent=True,\n",
+    "                                                           set_spectra=True\n",
+    "                                                          )\n",
+    "\n",
+    "for k, src in enumerate(sources):\n",
+    "    print (f\"{k}: {src.__class__.__name__}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can see that this method chooses multiple (here 2) extended components for the first three sources, a single extended component for sources 3 and 5, and compact ones for the rest.\n",
+    "\n",
+    "If we know something about the scene, we might want to customize the modeling. Let's assume that we know that object 0 is a star. We could just replace that source with a `PointSource`, but for clarity we rebuild the entire source list, making one change for source 0 and accepting everything else from the default initialization above:"
    ]
   },
   {
@@ -214,22 +243,45 @@
    "outputs": [],
    "source": [
     "sources = []\n",
-    "for k,src in enumerate(catalog):\n",
+    "for k,center in enumerate(centers):\n",
     "    if k == 0:\n",
-    "        new_source = scarlet.PointSource(model_frame, (src['y'], src['x']), observation)\n",
-    "    elif k == 1:\n",
-    "        new_source = scarlet.ExtendedSource(model_frame, (src['y'], src['x']), observation, K=2)\n",
+    "        new_source = scarlet.PointSource(model_frame, center, observation)\n",
+    "    elif k in [1, 2]:\n",
+    "        new_source = scarlet.ExtendedSource(model_frame, center, observation, K=2)\n",
+    "    elif k in [3, 5]:\n",
+    "         new_source = scarlet.ExtendedSource(model_frame, center, observation, K=1)\n",
     "    else:\n",
-    "        new_source = scarlet.ExtendedSource(model_frame, (src['y'], src['x']), observation, compact=True)\n",
-    "    sources.append(new_source)"
+    "         new_source = scarlet.ExtendedSource(model_frame, center, observation, compact=True)\n",
+    "    sources.append(new_source)\n",
+    "    \n",
+    "for k, src in enumerate(sources):\n",
+    "    print (f\"{k}: {src.__class__.__name__}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "These source are initialized independently, and their spectra, i.e. the amplitudes in every channel, assume that they are isolated.\n",
+    "We can make another sweep and determine their spectra so that their superposition best matches the observation. Above, this was done for use by the option `set_specta=True`. But we can call the linear solver directly:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "scarlet.initialization.set_spectra_to_match(sources, observation)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "## Create and Fit Model\n",
-    "The `Blend` class represents the sources as a tree and has the machinery to fit all of the sources to the given images. In this example the code is set to run for a maximum of 200 iterations, but will end early if the likelihood and all of the constraints converge."
+    "\n",
+    "The `Blend` class holds the list of sources and has the machinery to fit them to the given images. In this example the code is set to run for a maximum of 100 iterations, but will end early if the likelihood and all of the constraints converge."
    ]
   },
   {
@@ -304,6 +356,7 @@
    "source": [
     "### View Source Models\n",
     "\n",
+    "\n",
     "We will now inspect the model for each source, in its original frame and in its observed frame by leveraging the `show_sources` method:"
    ]
   },
@@ -425,10 +478,10 @@
     "# add two sources at their approximate locations\n",
     "model_frame_ = sources_[0].frame\n",
     "yx = (14., 44.)\n",
-    "new_source = scarlet.ExtendedSource(model_frame_, yx, observation, shifting=True)\n",
+    "new_source = scarlet.ExtendedSource(model_frame_, yx, observation, compact=True, shifting=True)\n",
     "sources_.append(new_source)\n",
     "yx = (42., 9.)\n",
-    "new_source = scarlet.ExtendedSource(model_frame_, yx, observation, shifting=True)\n",
+    "new_source = scarlet.ExtendedSource(model_frame_, yx, observation, compact=True, shifting=True)\n",
     "sources_.append(new_source)\n",
     "\n",
     "# generate a new Blend instance\n",

docs/1-concepts.ipynb

@@ -189,7 +189,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Source 0 is a `PointSource`, parameterized by a spectrum and a center, as we have seen above. Source 1 has multiple components (each with a spectrum and an image morphology), the rest have a single spectrum and image morphology.\n",
+    "Source 0 is a `PointSource`, parameterized by a spectrum and a center, as we have seen above. Source 1 has multiple components (each with a spectrum and an image morphology), others have a single spectrum and image morphology.\n",
     "\n",
     "`Parameter` is a souped up `numpy` array. It has a value and a name, as well as additonal attributes that store a `prior` and `constraint` that were enforced during optimization; the typical `step` size during optimization; an estimate of the standard deviate `std`; and whether the parameter was held `fixed`:"
    ]

docs/tutorials/display.ipynb

@@ -36,11 +36,11 @@
     "catalog = data[\"catalog\"]\n",
     "filters = data[\"filters\"]\n",
     "weights = 1/data[\"variance\"]\n",
-    "psfs = scarlet.ImagePSF(data[\"psfs\"])\n",
+    "psf = scarlet.ImagePSF(data[\"psfs\"])\n",
     "\n",
     "observation = scarlet.Observation(\n",
     "    images, \n",
-    "    psfs=psfs, \n",
+    "    psf=psf, \n",
     "    weights=weights, \n",
     "    channels=filters)\n",
     "\n",