diff --git a/Deblending/scarlet_tutorial.ipynb b/Deblending/scarlet_tutorial.ipynb
index 298ecaa6..e4ae262d 100644
--- a/Deblending/scarlet_tutorial.ipynb
+++ b/Deblending/scarlet_tutorial.ipynb
@@ -6,8 +6,8 @@
    "source": [
     "# Deblending with *Scarlet*\n",
     "<br>Owner(s): **Fred Moolekamp** ([@fred3m](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@fred3m))\n",
-    "<br>Last Verified to Run: **2019-08-16**\n",
-    "<br>Verified Stack Release: **w_2019-31**\n",
+    "<br>Last Verified to Run: **2020-07-10**\n",
+    "<br>Verified Stack Release: **v20.0.0**\n",
     "\n",
     "The purpose of this tutorial is to familiarize you with the basics of using *scarlet* to model blended scenes, and how tweaking various objects and parameters affects the resulting model. A tutorial that is more specific to using scarlet in the context of the LSST DM Science Pipelines is also available.\n",
     "\n",
@@ -74,124 +74,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Display the sources\n",
-    "# Display the sources\n",
-    "def display_sources(sources, observation, norm=None, subset=None, combine=False, show_sed=True):\n",
-    "    \"\"\"Display the data and model for all sources in a blend\n",
-    "\n",
-    "    This convenience function is used to display all (or a subset) of\n",
-    "    the sources and (optionally) their SED's.\n",
-    "    \"\"\"\n",
-    "    if subset is None:\n",
-    "        # Show all sources in the blend\n",
-    "        subset = range(len(sources))\n",
-    "    for m in subset:\n",
-    "        # Load the model for the source\n",
-    "        src = sources[m]\n",
-    "        if hasattr(src, \"components\"):\n",
-    "            components = len(src.components)\n",
-    "        else:\n",
-    "            components = 1\n",
-    "        # Convolve the model with the psfs in the observation\n",
-    "        model = observation.render(src.get_model())\n",
-    "        # Extract the bounding box that contains the non-zero\n",
-    "        # pixels in the model\n",
-    "        bbox = scarlet.bbox.trim(np.sum(model, axis=0), min_value=1e-2)\n",
-    "        bb = (slice(None), *bbox.slices)\n",
-    "        # Adjust the stretch based on the maximum flux in the model for the current source\n",
-    "        if model.max() > 10 * bg_rms.max():\n",
-    "            norm = AsinhMapping(minimum=model.min(), stretch=model.max()*.05, Q=10)\n",
-    "        else:\n",
-    "            norm = LinearMapping(minimum=model.min(), maximum=model.max())\n",
-    "\n",
-    "        # Select the image patch the overlaps with the source and convert it to an RGB image\n",
-    "        img_rgb = scarlet.display.img_to_rgb(images[bb], norm=norm)\n",
-    "\n",
-    "        # Build a model for each component in the model\n",
-    "        if hasattr(src, \"components\"):\n",
-    "            rgb = []\n",
-    "            for component in src.components:\n",
-    "                # Convert the model to an RGB image\n",
-    "                _model = observation.get_model(component.get_model())\n",
-    "                _rgb = scarlet.display.img_to_rgb(_model[bb], norm=norm)\n",
-    "                rgb.append(_rgb)\n",
-    "        else:\n",
-    "            # There is only a single component\n",
-    "            rgb = [scarlet.display.img_to_rgb(model[bb], norm=norm)]\n",
-    "\n",
-    "        # Display the image and model\n",
-    "        figsize = [6,3]\n",
-    "        columns = 2\n",
-    "        # Calculate the number of columns needed and shape of the figure\n",
-    "        if show_sed:\n",
-    "            figsize[0] += 3\n",
-    "            columns += 1\n",
-    "        if not combine:\n",
-    "            figsize[0] += 3*(components-1)\n",
-    "            columns += components-1\n",
-    "        # Build the figure\n",
-    "        fig = plt.figure(figsize=figsize)\n",
-    "        ax = [fig.add_subplot(1,columns,n+1) for n in range(columns)]\n",
-    "        ax[0].imshow(img_rgb)\n",
-    "        ax[0].set_title(\"Data: Source {0}\".format(m))\n",
-    "        for n, _rgb in enumerate(rgb):\n",
-    "            ax[n+1].imshow(_rgb)\n",
-    "            if combine:\n",
-    "                ax[n+1].set_title(\"Initial Model\")\n",
-    "            else:\n",
-    "                ax[n+1].set_title(\"Component {0}\".format(n))\n",
-    "        if show_sed:\n",
-    "            if components > 1:\n",
-    "                for comp in src:\n",
-    "                    ax[-1].plot(comp.sed)\n",
-    "            else:\n",
-    "                ax[-1].plot(src.sed)\n",
-    "            ax[-1].set_title(\"SED\")\n",
-    "            ax[-1].set_xlabel(\"Band\")\n",
-    "            ax[-1].set_ylabel(\"Intensity\")\n",
-    "        # Mark the current source in the image\n",
-    "        if components > 1:\n",
-    "            y,x = src.components[0].pixel_center\n",
-    "        else:\n",
-    "            y,x = src.pixel_center\n",
-    "        ax[0].plot(x-bb[2].start, y-bb[1].start, 'rx', mew=2)\n",
-    "        plt.tight_layout()\n",
-    "        plt.show()\n",
-    "\n",
-    "def display_model_residual(images, blend, peaks, norm):\n",
-    "    \"\"\"Display the data, model, and residual for a given result\n",
-    "    \"\"\"\n",
-    "    model = blend.get_model()\n",
-    "    residual = images-model\n",
-    "    print(\"Data range: {0:.3f} to {1:.3f}\\nresidual range: {2:.3f} to {3:.3f}\\nrms: {4:.3f}\".format(\n",
-    "        np.min(images),\n",
-    "        np.max(images),\n",
-    "        np.min(residual),\n",
-    "        np.max(residual),\n",
-    "        np.sqrt(np.std(residual)**2+np.mean(residual)**2)\n",
-    "    ))\n",
-    "    # Create RGB images\n",
-    "    img_rgb = scarlet.display.img_to_rgb(images, norm=norm)\n",
-    "    model_rgb = scarlet.display.img_to_rgb(model, norm=norm)\n",
-    "    residual_rgb = scarlet.display.img_to_rgb(residual)\n",
-    "\n",
-    "    # Show the data, model, and residual\n",
-    "    fig = plt.figure(figsize=(15,5))\n",
-    "    ax = [fig.add_subplot(1,3,n+1) for n in range(3)]\n",
-    "    ax[0].imshow(img_rgb)\n",
-    "    ax[0].set_title(\"Data\")\n",
-    "    ax[1].imshow(model_rgb)\n",
-    "    ax[1].set_title(\"Model\")\n",
-    "    ax[2].imshow(residual_rgb)\n",
-    "    ax[2].set_title(\"Residual\")\n",
-    "    for k,component in enumerate(blend.components):\n",
-    "        y,x = component.pixel_center\n",
-    "        #px, py = peaks[k]\n",
-    "        ax[0].plot(x, y, \"gx\")\n",
-    "        #ax[0].plot(px, py, \"rx\")\n",
-    "        ax[1].text(x, y, k, color=\"r\")\n",
-    "    plt.show()\n",
-    "\n",
     "def show_psfs(psfs, filters, norm=None):\n",
     "    rows = int(np.ceil(len(psfs)/3))\n",
     "    columns = min(len(psfs), 3)\n",
@@ -256,6 +138,7 @@
     "peaks = data[\"peaks\"]\n",
     "psfs = data[\"psfs\"]\n",
     "filters = [\"G\", \"R\", \"I\", \"Z\", \"Y\"]\n",
+    "\n",
     "# Only a rough estimate of the background is needed\n",
     "# to initialize and resize the sources\n",
     "bg_rms = np.std(images, axis=(1,2))\n",
@@ -263,6 +146,7 @@
     "\n",
     "# Use Asinh scaling for the images\n",
     "norm = AsinhMapping(minimum=images.min(), stretch=images.max()/20, Q=10)\n",
+    "\n",
     "# Convert the image to an RGB image\n",
     "img_rgb = scarlet.display.img_to_rgb(images, norm=norm)\n",
     "plt.imshow(img_rgb)\n",
@@ -272,6 +156,29 @@
     "plt.show()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Psf display "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can also takea  look at the `psfs`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "show_psfs(psfs, filters)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -285,6 +192,15 @@
     "See https://fred3m.github.io/scarlet/user_docs.html#Frame-and-Observation for more on `Frame`s and `Observation`s."
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from functools import partial"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -292,17 +208,20 @@
    "outputs": [],
    "source": [
     "# Create a PSF image of a narrow gaussian to use as our image PSF\n",
-    "model_psf = scarlet.generate_psf_image(scarlet.gaussian, shape=(41, 41), amplitude=1, sigma=.9)\n",
+    "amplitude=1.\n",
+    "model_psf = scarlet.PSF(partial(scarlet.psf.gaussian, sigma=0.9), shape=(None,41, 41)).image * amplitude\n",
     "model_psf /= model_psf.sum()\n",
     "# Make sure that the observation PSF is normalized (otherwise the scaling in PSF matching might be off)\n",
     "psfs = psfs / psfs.sum(axis=(1, 2))[:, None, None]\n",
+    "\n",
     "# Create the initial frame (metadata for the model).\n",
     "# Note that we initialized a PSF with shape (Ny, Nx) but a frame\n",
     "# expects a PSf with shape (bands, Ny, Nx), so we have to\n",
     "# broadcast the model_psf into an extra dimension\n",
-    "frame = scarlet.Frame(images.shape, psfs=model_psf[None])\n",
+    "frame = scarlet.Frame(images.shape, psfs=model_psf, channels=filters)\n",
+    "\n",
     "# Create our observation\n",
-    "observation = scarlet.Observation(images, psfs=psfs).match(frame)"
+    "observation = scarlet.Observation(images, psfs=psfs, channels=filters, weights=weights).match(frame)"
    ]
   },
   {
@@ -315,22 +234,32 @@
     "\n",
     "The different classes that inherit from `Source` mainly differ in how they are initialized, and otherwise behave similarly during the optimization routine. This section illustrates the differences between different source initialization classes.\n",
     "\n",
-    "The simplest source is a single component intialized with only a single pixel (at the center of the object) turned on.\n",
-    "\n",
     "### <span style=\"color:red\"> *WARNING* </span>\n",
     "Scarlet accepts source positions using the numpy/C++ convention of (y,x), which is different than the astropy and LSST stack convention of (x,y)."
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Below we demonstrate the usage of `ExtendedSource`, which initializes each object as a single component with maximum flux at the peak that falls off monotonically and has 180 degree symmetry."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "sources = [scarlet.PointSource(frame, (peak[1], peak[0]), observation) for peak in peaks]\n",
+    "sources = [scarlet.ExtendedSource(frame, (peak[1], peak[0]), observation) for peak in peaks]\n",
     "\n",
     "# Display the initial guess for each source\n",
-    "display_sources(sources, observation, norm=norm)"
+    "scarlet.display.show_sources(sources,\n",
+    "                             norm=norm,\n",
+    "                             observation=observation,\n",
+    "                             show_rendered=True,\n",
+    "                             show_observed=True)\n",
+    "plt.show()"
    ]
   },
   {
@@ -339,7 +268,7 @@
    "source": [
     "## Exercise:\n",
     "\n",
-    "* Experiment with the above code by using `ExtendedSource`, which initializes each object as a single component with maximum flux at the peak that falls off monotonically and has 180 degree symmetry; and using `MultiComponentSource`, which models a source as two components (a bulge and a disk) that are each symmetric and montonically decreasing from the peak.\n",
+    "* Experiment with the above code by using ; and using `MultiComponentSource`, which models a source as two components (a bulge and a disk) that are each symmetric and montonically decreasing from the peak.\n",
     "\n",
     "# Deblending a scene\n",
     "\n",
@@ -369,29 +298,20 @@
    "outputs": [],
    "source": [
     "# Fit the data until the relative error is <= 1e-3,\n",
-    "# for a maximum of 100 iterations\n",
-    "blend.fit(100, e_rel=1e-3)\n",
-    "print(\"Deblending completed in {0} iterations\".format(blend.it))\n",
-    "display_model_residual(images, blend, peaks, norm)\n",
-    "display_sources(sources, observation, norm)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Exercises\n",
-    "\n",
-    "* Experiment by running the above code using different source models (for example `ExtendedSource`) to see how initializtion affects the belnding results."
+    "# for a maximum of 200 iterations\n",
+    "blend = scarlet.Blend(sources, observation)\n",
+    "%time blend.fit(200, e_rel = 1e-3)\n",
+    "print(\"scarlet ran for {0} iterations to logL = {1}\".format(len(blend.loss), -blend.loss[-1]))\n",
+    "plt.plot(-np.array(blend.loss))\n",
+    "plt.xlabel('Iteration')\n",
+    "plt.ylabel('log-Likelihood')"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Update Functions\n",
-    "\n",
-    "The above models used the default constraints: perfect symmetry and a weighted monotonicity that decreases from the peak. So each source is created with the update function, equivalent to the following:"
+    "There are two options for displaying the scene, using `scarlet.display.show_scene` function. This shows the model along with the observation information and the residuals defined as: `observation.images - model`. "
    ]
   },
   {
@@ -400,54 +320,20 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from scarlet import measurement, update\n",
-    "\n",
-    "class MySource(scarlet.ExtendedSource):\n",
-    "    def update(self):\n",
-    "        \"\"\"Default update parameters for a PointSource\n",
-    "        This method can be overwritten if a different set of constraints\n",
-    "        or update functions is desired.\n",
-    "        \"\"\"\n",
-    "        # Keep track of the iteration so we can skip certain updates\n",
-    "        # based on the iteration number\n",
-    "        if self._parent is None:\n",
-    "            it = 0\n",
-    "        else:\n",
-    "            it = self._parent.it\n",
-    "\n",
-    "        # Update the central pixel location (pixel_center)\n",
-    "        self.pixel_center = measurement.max_pixel(self.morph, self.pixel_center)\n",
-    "\n",
-    "        if self.symmetric:\n",
-    "            # Update the centroid position every 5th iteration\n",
-    "            if it % 5 == 0:\n",
-    "                self.pixel_center, self.shift = measurement.psf_weighted_centroid(self.morph,\n",
-    "                                                                                  self._centroid_weight,\n",
-    "                                                                                  self.pixel_center)\n",
-    "\n",
-    "            # make the morphology perfectly symmetric\n",
-    "            update.symmetric(self, algorithm=\"kspace\")\n",
-    "\n",
-    "        if self.monotonic:\n",
-    "            # make the morphology monotonically decreasing\n",
-    "            update.monotonic(self, self.pixel_center)\n",
-    "\n",
-    "        update.positive(self)  # Make the SED and morph non-negative\n",
-    "        update.normalized(self)  # Use MORPH_MAX normalization\n",
-    "        return self"
+    "scarlet.display.show_scene(sources, norm=norm,linear=True, \n",
+    "                           observation=observation, show_observed=True, \n",
+    "                           label_sources=True, \n",
+    "                           show_rendered=True, \n",
+    "                           show_residual=True\n",
+    "                          )\n",
+    "plt.show()"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Exercises\n",
-    "\n",
-    "In the code above modify the `MySource.update` method to use different constraints and run the code block following this cell for the following:\n",
-    "\n",
-    "* Add the keyword argument `use_nearest=True` to `update.monotonic` and see how that affects the results.\n",
-    "\n",
-    "*  Comment out the lines that set `self.pixel_center` so that positions are not updated and see how that changes the results."
+    "You can also do it by hand. "
    ]
   },
   {
@@ -456,22 +342,42 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "offset = np.random.rand(len(peaks), 2)/2\n",
-    "sources = [MySource(frame, (peak[1]+offset[k][0], peak[0]+offset[k][1]), observation, bg_rms) for k, peak in enumerate(peaks)]\n",
-    "blend = scarlet.Blend(sources, observation)\n",
-    "blend.fit(100, e_rel=1e-3)\n",
-    "print(\"Deblending completed in {0} iterations\".format(blend.it))\n",
-    "display_model_residual(images, blend, peaks, norm)\n",
-    "display_sources(sources, observation, norm)"
+    "# Load the model and calculate the residual\n",
+    "model = blend.get_model()\n",
+    "model_ = observation.render(model)  # adapt model to observations. \n",
+    "residual = images-model_\n",
+    "\n",
+    "# Create RGB images\n",
+    "model_rgb = scarlet.display.img_to_rgb(model_, norm=norm)\n",
+    "residual_rgb = scarlet.display.img_to_rgb(residual)\n",
+    "\n",
+    "# Show the data, model, and residual\n",
+    "fig = plt.figure(figsize=(15,5))\n",
+    "ax = [fig.add_subplot(1,3,n+1) for n in range(3)]\n",
+    "ax[0].imshow(img_rgb)\n",
+    "ax[0].set_title(\"Data\")\n",
+    "ax[1].imshow(model_rgb)\n",
+    "ax[1].set_title(\"Model\")\n",
+    "ax[2].imshow(residual_rgb)\n",
+    "ax[2].set_title(\"Residual\")\n",
+    "\n",
+    "for k,component in enumerate(blend):\n",
+    "    y,x = component.center\n",
+    "    ax[0].text(x, y, k, color=\"w\")\n",
+    "    ax[1].text(x, y, k, color=\"w\")\n",
+    "    ax[2].text(x, y, k, color=\"w\")\n",
+    "plt.show()"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Extra Credit\n",
+    "## Exercises\n",
+    "\n",
+    "* Experiment by running the above code using different source models (for example `MultiComponentSource` or `PointSource`) to see how initializtion affects the belnding results.\n",
     "\n",
-    "* Modify `MySource` to create a source that has multiple components. You might want to take a look at [MultiComponentSource](https://github.com/fred3m/scarlet/blob/master/scarlet/source.py#L478) to get an understanding about how to do this, but try your own initialization function."
+    "* Change the value of `e_rel` in the above fit and try to understand how it affects the results. "
    ]
   }
  ],
@@ -491,9 +397,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.2"
+   "version": "3.7.6"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }