Repository: ICB-DCM/pyPESTO

Issue:
Save optimizer trace
Currently, no optimizer trace is/can be saved. This could be done via the callback-input in the minimize-framework...

Discussion:
yannikschaelte: could you set the target to develop? we now have the "branching policy" of first committing to develop, and only from there fully functional releases to master.
codacy-bot: 
![Codacy](https://app.codacy.com/assets/images/favicon.png) Here is an overview of what got changed by this pull request:


```diff

Issues
======
- Added 7
           

Complexity increasing per file
==============================
- pypesto/visualize/profiles.py  1
- pypesto/visualize/parameters.py  1
         
```


See the complete overview on [Codacy](https://app.codacy.com/app/dweindl/pyPESTO/pullRequest?prid=3063240&bid=11238496)
        
codecov-io: # [Codecov](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114?src=pr&el=h1) Report
> :exclamation: No coverage uploaded for pull request base (`develop@2862636`). [Click here to learn what that means](https://docs.codecov.io/docs/error-reference#section-missing-base-commit).
> The diff coverage is `36.98%`.

[![Impacted file tree graph](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114/graphs/tree.svg?width=650&token=vaUsU0G1LY&height=150&src=pr)](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##             develop     #114   +/-   ##
==========================================
  Coverage           ?   82.34%           
==========================================
  Files              ?       37           
  Lines              ?     2016           
  Branches           ?        0           
==========================================
  Hits               ?     1660           
  Misses             ?      356           
  Partials           ?        0
```


| [Impacted Files](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [pypesto/visualize/\_\_init\_\_.py](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114/diff?src=pr&el=tree#diff-cHlwZXN0by92aXN1YWxpemUvX19pbml0X18ucHk=) | `100% <100%> (ø)` | |
| [pypesto/visualize/optimizer\_history.py](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114/diff?src=pr&el=tree#diff-cHlwZXN0by92aXN1YWxpemUvb3B0aW1pemVyX2hpc3RvcnkucHk=) | `12.79% <12.79%> (ø)` | |
| [pypesto/visualize/reference\_points.py](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114/diff?src=pr&el=tree#diff-cHlwZXN0by92aXN1YWxpemUvcmVmZXJlbmNlX3BvaW50cy5weQ==) | `34.37% <34.37%> (ø)` | |
| [pypesto/visualize/profiles.py](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114/diff?src=pr&el=tree#diff-cHlwZXN0by92aXN1YWxpemUvcHJvZmlsZXMucHk=) | `70.66% <48.57%> (ø)` | |
| [pypesto/visualize/waterfall.py](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114/diff?src=pr&el=tree#diff-cHlwZXN0by92aXN1YWxpemUvd2F0ZXJmYWxsLnB5) | `80% <52.63%> (ø)` | |
| [pypesto/visualize/parameters.py](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114/diff?src=pr&el=tree#diff-cHlwZXN0by92aXN1YWxpemUvcGFyYW1ldGVycy5weQ==) | `80% <65.78%> (ø)` | |
| [pypesto/profile/profile.py](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114/diff?src=pr&el=tree#diff-cHlwZXN0by9wcm9maWxlL3Byb2ZpbGUucHk=) | `88.11% <66.66%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114?src=pr&el=footer). Last update [2862636...308a041](https://codecov.io/gh/ICB-DCM/pyPESTO/pull/114?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

paulstapor: Okay, now, only super little stuff is missing:
  I'd like to cover the list plotting functionality with more tests,
  I'd like to implement legends for different results,
  I'm afraid it will be necessary to fill some documentation things I may have missed while coding,
  I think we should implement at least a part of this functionality in at least one example...

Besides that: All major implementations I wanted to have in here for plotting is done.
dweindl: >  This pr will fix #108, #115, #111, #98 

Note that this will auto-close only the first one. If individual ones are included in other commit messages, not problem.
dweindl: > +1,565 −179

What happened to small pull requests? :thinking: 
paulstapor: > > +1,565 −179
> 
> What happened to small pull requests? 🤔

That's an exception to the rule ;)
No, kidding: The idea is that this pr involved some changes throughout all plotting routines. I had to change some stuff over and over again, since not everything worked out as I wanted giving the first restructuring. I didn't want to merge as long as I knew that I would have to change many things again.
FFroehlich: > > > +1,565 −179
> > 
> > 
> > What happened to small pull requests? 🤔
> 
> That's an exception to the rule ;)
> No, kidding: The idea is that this pr involved some changes throughout all plotting routines. I had to change some stuff over and over again, since not everything worked out as I wanted giving the first restructuring. I didn't want to merge as long as I knew that I would have to change many things again.
> 
> Edit: anyhow: More than 50% of these 1565 lines are comments anyway...

The objective function history plotting could have been put in another PR, but other than that I can appreciate long PR if they include a lot of comments + unittests ;)
paulstapor: > The objective function history plotting could have been put in another PR, but other than that I can appreciate long PR if they include a lot of comments + unittests ;)

Test coverage increases, so should be fine.

Yes, the objective history could have been taken out. Besides that, all the three (now four) other lowlevel plotting routines had to be adapted, which directly multiplies the affected lines by a factor of three, four. Handling of the introduced reference points is somewhat intertwined with the changes in the high- and the lowlevel routines as well, the same holds for the legends, the customizable colors, etc...

paulstapor: Perfect, thanks already for all the input, guys! 👍 :)
Looking forward to getting even more. :)
dweindl: Not sure when and where the pytest warnings came up. If here, please fix.
paulstapor: @dweindl , @yannikschaelte :
...it seems like the version of petab changed on pypi and now, some of the corresponding tests for pypesto don't work any more...
dweindl: Yes, that was to expect. For this PR you can fix the PEtab version to the previous one. But should be upgraded soon.
FFroehlich: > unnest a nested if in color clustering a little
what?
yannikschaelte: >Yes, that was to expect. For this PR you can fix the PEtab version to the previous one. But should be upgraded soon.

Just set ``petab==0.0.0a10`` in setup.py I guess
paulstapor: > > unnest a nested if in color clustering a little
> 
> what?

Damn, I confused that with a typo fix... 😂 Sorry for that, but it's gonna be a squashed merge anyway...^^

Patch:
diff --git a/pypesto/profile/profile.py b/pypesto/profile/profile.py
index 865772b..c6890ed 100644
--- a/pypesto/profile/profile.py
+++ b/pypesto/profile/profile.py
@@ -297,11 +297,17 @@ def walk_along_profile(current_profile,
         # IMPORTANT: This optimization will need a proper exception
         # handling (coming soon)
         optimizer_result = optimizer.minimize(problem, startpoint, 0)
+        if optimizer_result["grad"] is not None:
+            gradnorm = np.linalg.norm(optimizer_result["grad"][
+                                          problem.x_free_indices])
+        else:
+            gradnorm = None
+
         current_profile.append_profile_point(
             optimizer_result.x,
             optimizer_result.fval,
             np.exp(global_opt - optimizer_result.fval),
-            np.linalg.norm(optimizer_result.grad[problem.x_free_indices]),
+            gradnorm,
             optimizer_result.exitflag,
             optimizer_result.time,
             optimizer_result.n_fval,
@@ -404,12 +410,17 @@ def fill_profile_list(
             number of parameters in the unreduced problem
         """
 
+    if optimize_result["grad"] is not None:
+        gradnorm = np.linalg.norm(optimize_result["grad"])
+    else:
+        gradnorm = None
+
     # create blanko profile
     new_profile = ProfilerResult(
         optimize_result["x"],
         np.array([optimize_result["fval"]]),
         np.array([1.]),
-        np.linalg.norm(optimize_result["grad"]),
+        gradnorm,
         optimize_result["exitflag"],
         np.array([0.]),
         np.array([0.]),
diff --git a/pypesto/visualize/__init__.py b/pypesto/visualize/__init__.py
index ca15541..cad6423 100644
--- a/pypesto/visualize/__init__.py
+++ b/pypesto/visualize/__init__.py
@@ -7,22 +7,38 @@ import pypesto.visualize.
 
 """
 
+from .reference_points import (ReferencePoint,
+                               create_references)
+from .clust_color import (assign_clusters,
+                          assign_clustered_colors,
+                          assign_colors)
+from .misc import (process_result_list,
+                   process_offset_y,
+                   process_y_limits)
 from .waterfall import (waterfall,
                         waterfall_lowlevel)
-from .clust_color import (assign_clusters,
-                          assign_clustered_colors)
 from .parameters import (parameters,
                          parameters_lowlevel)
+from .optimizer_history import (optimizer_history,
+                                optimizer_history_lowlevel)
 from .profiles import (profiles,
                        profiles_lowlevel,
                        profile_lowlevel)
 
-__all__ = ["waterfall",
-           "waterfall_lowlevel",
+__all__ = ["ReferencePoint",
+           "create_references",
            "assign_clusters",
            "assign_clustered_colors",
+           "assign_colors",
+           "process_result_list",
+           "process_offset_y",
+           "process_y_limits",
+           "waterfall",
+           "waterfall_lowlevel",
            "parameters",
            "parameters_lowlevel",
+           "optimizer_history",
+           "optimizer_history_lowlevel",
            "profiles",
            "profiles_lowlevel",
            "profile_lowlevel"
diff --git a/pypesto/visualize/clust_color.py b/pypesto/visualize/clust_color.py
index e287419..84b39f1 100644
--- a/pypesto/visualize/clust_color.py
+++ b/pypesto/visualize/clust_color.py
@@ -1,5 +1,4 @@
 from scipy import cluster
-import matplotlib.colors as colors
 import matplotlib.cm as cm
 import numpy as np
 
@@ -31,6 +30,8 @@ def assign_clusters(vals):
     # sanity checks
     if vals is None or len(vals) == 0:
         return [], [], []
+    elif len(vals) == 1:
+        return np.array([1]), np.array([1.]), np.array([0])
 
     # linkage requires (n, 1) data array
     vals = np.reshape(vals, (-1, 1))
@@ -59,7 +60,7 @@ def assign_clustered_colors(vals):
     Returns
     -------
 
-    Col: list of RGB
+    colors: list of RGBA
         One for each element in 'vals'.
     """
 
@@ -68,44 +69,131 @@ def assign_clustered_colors(vals):
         return []
 
     # assign clusters
-    clust, clustsize, ind_clust = assign_clusters(vals)
-
-    # assign colors
-    vmax = max(clust) - sum(clustsize == 1)
-    cnorm = colors.Normalize(vmin=0, vmax=vmax)
-    scalarmap = cm.ScalarMappable(norm=cnorm)
-
-    # colors for each cluster and one-size clusters
-    cols_clust = scalarmap.to_rgba(range(vmax + 1))
-
-    # grey color for 1-size clusters as the last color in 'cols_clust'
-    cols_clust[vmax] = (0.7, 0.7, 0.7, 1)
-
-    # pre-array of indices for colors
-    ind_col = np.zeros(len(ind_clust))
-
-    # number of colors assigned to cluster
-    sum_col = 0
-
-    # assign color indices for each point
-    for iind, value_ind in enumerate(ind_clust):
-        # if cluster size > 1
-        if clustsize[value_ind] > 1:
-            # assign a color with is not grey
-            ind_col[iind] = sum_col
-            # if element is not the last one in 'ind_col'
-            if iind < len(ind_clust) - 1:
-                # if the next element does not belongs to the seem cluster
-                if value_ind != ind_clust[iind + 1]:
-                    # use the next color
-                    sum_col = sum_col + 1
-        # if cluster size = 1
-        else:
-            # use grey
-            ind_col[iind] = vmax
-
-    # indices for colors
-    ind_col = [int(ind) for ind in ind_col]
-    col = cols_clust[ind_col]
-
-    return col
+    clust, cluster_size, cluster_indices = assign_clusters(vals)
+
+    # create list of colors, which has the correct shape
+    n_clusters = max(clust) - sum(cluster_size == 1)
+
+    # fill color array from colormap
+    colormap = cm.ScalarMappable().to_rgba
+    color_list = colormap(np.linspace(0., 1., n_clusters))
+
+    # We have clustered the results. However, clusters may have size 1,
+    # so we need to rearrange cluster_indices in order to map all
+    # 1-element-clusters to the cluster with index n_clusters
+
+    # create a dummy variable with the colors
+    colors = np.zeros((clust.size, 4))
+
+    # run through the cluster_indices and collect which element
+    # is in 1-element-clusters, and which is in real clusters
+    no_clusters = np.where(cluster_size == 1)[0]
+    real_clusters = np.unique(np.where(cluster_size > 1)[0])
+
+    # assign colors to real clusters
+    for icol, iclust in enumerate(real_clusters):
+        ind_of_iclust = np.argwhere(cluster_indices == iclust).flatten()
+        colors[ind_of_iclust, :] = color_list[icol, :]
+
+    # assign color to non-clustered indices
+    for noclust in no_clusters:
+        ind_of_noclust = np.argwhere(noclust in no_clusters).flatten()
+        colors[ind_of_noclust, :] = [0.7, 0.7, 0.7, 1]
+
+    return colors
+
+
+def assign_colors(vals, colors=None):
+    """
+    Assign colors or format user specified colors.
+
+    Parameters
+    ----------
+
+    vals: numeric list or array
+        List to be clustered and assigned colors.
+
+    colors: list, or RGBA, optional
+        list of colors, or single color
+
+    Returns
+    -------
+
+    colors: list of RGBA
+        One for each element in 'vals'.
+    """
+
+    # sanity checks
+    if vals is None or len(vals) == 0:
+        return np.array([])
+
+    # if the user did not specify any colors:
+    if colors is None:
+        return assign_clustered_colors(vals)
+
+    # The user passed values and colors: parse them first!
+    # we want everything to be numpy arrays, to not check everytime wther a
+    # list was passed or an ndarray
+    colors = np.array(colors)
+
+    # Get number of elements and use user assigned colors
+    n_vals = len(vals) if isinstance(vals, list) else vals.size
+
+    # Two usages are possible: One color for the whole data set, or one
+    # color for each value:
+    if colors.size == 4:
+        # Only one color was passed: flatten in case and repeat n_vals times
+        if colors.ndim == 2:
+            colors = colors[0]
+        return np.array([colors] * n_vals)
+    else:
+        if colors.shape[1] == 4 and n_vals == colors.shape[0]:
+            return colors
+        elif colors.shape[0] == 4:
+            colors = np.transpose(colors)
+            if n_vals == colors.shape[0]:
+                return colors
+
+        # Shape of array did not match n_vals. Error due to size mismatch:
+        raise ('Incorrect color input. Colors must be specified either as '
+               'list of [r, g, b, alpha] with length equal to function '
+               'values Number of function (here: ' + str(n_vals) + '), or as '
+               'one single [r, g, b, alpha] color.')
+
+
+def assign_colors_for_result_list(num_results, colors=None):
+    """
+    Creates a list of colors for a list of pypesto.Result objects or checks
+    a user-provided list of colors and uses this if everything is ok
+
+    Parameters
+    ----------
+
+    num_results: int
+        number of results in list
+
+    colors: list, or RGBA, optional
+        list of colors, or single color
+
+    Returns
+    -------
+
+    colors: list of RGBA
+        One for each element in 'vals'.
+    """
+
+    # if the user did not specify any colors:
+    if colors is None:
+        dummy_clusters = np.array(list(range(num_results)) * 2)
+        colors = assign_colors(dummy_clusters)
+        real_indices = list(range(int(colors.shape[0] / 2)))
+        return colors[real_indices]
+
+    # if the user specified color lies does not match the number of results
+    if len(colors) != num_results:
+        raise ('Incorrect color input. Colors must be specified either as '
+               'list of [r, g, b, alpha] with length equal to function '
+               'values Number of function (here: ' + str(num_results) + '), '
+               'or as one single [r, g, b, alpha] color.')
+
+    return colors
diff --git a/pypesto/visualize/misc.py b/pypesto/visualize/misc.py
new file mode 100644
index 0000000..fb23511
--- /dev/null
+++ b/pypesto/visualize/misc.py
@@ -0,0 +1,162 @@
+import numpy as np
+import warnings
+from .clust_color import assign_colors
+from .clust_color import assign_colors_for_result_list
+
+
+def process_result_list(results, colors=None, legends=None):
+    """
+    assigns colors and legends to a list of results, chekc user provided lists
+
+    Parameters
+    ----------
+
+    results: list or pypesto.Result
+        list of pypesto.Result objects or a single pypesto.Result
+
+    colors: list, optional
+        list of RGBA colors
+
+    legends: str or list
+        labels for line plots
+
+    Returns
+    -------
+
+    results: list of pypesto.Result
+       list of pypesto.Result objects
+
+    colors: list of RGBA
+        One for each element in 'results'.
+
+    legends: list of str
+        labels for line plots
+    """
+
+    # check how many results were passed
+    single_result = False
+    legend_error = False
+    if isinstance(results, list):
+        if len(results) == 1:
+            single_result = True
+    else:
+        single_result = True
+        results = [results]
+
+    # handle results according to their number
+    if single_result:
+        # assign colors and create list for later handling
+        if colors is not None:
+            colors = assign_colors(results, colors)
+        colors = [colors]
+
+        # create list of legends for later handling
+        if not isinstance(legends, list):
+            legends = [legends]
+    else:
+        # if more than one result is passed, we use one color per result
+        colors = assign_colors_for_result_list(len(results), colors)
+
+        # check whether list of legends has the correct length
+        if legends is None:
+            # No legends were passed: create some custom legends
+            legends = []
+            for i_leg in range(len(results)):
+                legends.append('Result ' + str(i_leg))
+        else:
+            # legends were passed by user: check length
+            if isinstance(legends, list):
+                if len(legends) != len(results):
+                    legend_error = True
+            else:
+                legend_error = True
+
+    # size of legend list and size of results does not match
+    if legend_error:
+        raise ('List of results passed and list of labels do not have the'
+               ' same length but should. Stopping.')
+
+    return results, colors, legends
+
+
+def process_offset_y(offset_y, scale_y, min_val):
+    """
+    compute offset for y-axis, depend on user settings
+
+    Parameters
+    ----------
+
+    offset_y: float
+       value for offsetting the later plotted values, in order to ensure
+       positivity if a semilog-plot is used
+
+    min_val: float
+       Can be 'lin' or 'log10', specifying whether values should be plotted
+       on linear or on log10-scale
+
+    min_val: float
+        Smallest value to be plotted
+
+    Returns
+    -------
+
+    offset_y: float
+       value for offsetting the later plotted values, in order to ensure
+       positivity if a semilog-plot is used
+    """
+
+    # check whether the offset specified by the user is sufficient
+    if offset_y is not None:
+        if (scale_y == 'log10') and (min_val + offset_y <= 0.):
+            warnings.warn("Offset specified by user is insufficient. "
+                          "Ignoring specified offset and using " +
+                          str(np.abs(min_val) + 1.) + " instead.")
+            offset_y = 1. - min_val
+    else:
+        # check whether scaling is lin or log10
+        if scale_y == 'lin':
+            # linear scaling doesn't need any offset
+            offset_y = 0.
+        else:
+            # log10 scaling does need offset
+            offset_y = 1. - min_val
+
+    return offset_y
+
+
+def process_y_limits(ax, y_limits):
+    """
+    apply user specified limits of y-axis
+
+    Parameters
+    ----------
+
+    ax: matplotlib.Axes, optional
+        Axes object to use.
+
+    y_limits: ndarray
+       y_limits, minimum and maximum, for current axes object
+
+    min_val: float
+        Smallest value to be plotted
+
+    Returns
+    -------
+
+    ax: matplotlib.Axes, optional
+        Axes object to use.
+    """
+
+    # apply y-limits, if they were specified by the user
+    if y_limits is not None:
+        y_limits = np.array(y_limits)
+        if y_limits.size == 1:
+            tmp_y_limits = ax.get_ylim()
+            y_limits = [tmp_y_limits[0], y_limits]
+        else:
+            y_limits = [y_limits[0], y_limits[1]]
+
+        # set limits
+        ax.set_ylim(y_limits)
+
+    return ax
diff --git a/pypesto/visualize/optimizer_history.py b/pypesto/visualize/optimizer_history.py
new file mode 100644
index 0000000..10eb68f
--- /dev/null
+++ b/pypesto/visualize/optimizer_history.py
@@ -0,0 +1,387 @@
+import matplotlib.pyplot as plt
+from matplotlib.ticker import MaxNLocator
+import numpy as np
+from .reference_points import create_references
+from .clust_color import assign_colors
+from .misc import process_result_list
+from .misc import process_y_limits
+from .misc import process_offset_y
+
+
+def optimizer_history(results,
+                      ax=None,
+                      size=(18.5, 10.5),
+                      trace_x='steps',
+                      trace_y='fval',
+                      scale_y='log10',
+                      offset_y=None,
+                      colors=None,
+                      y_limits=None,
+                      start_indices=None,
+                      reference=None,
+                      legends=None):
+    """
+    Plot history of optimizer. Can plot either the history of the cost
+    function or of the gradient norm, over either the optimizer steps or
+    the computation time.
+
+    Parameters
+    ----------
+
+    results: pypesto.Result or list
+        Optimization result obtained by 'optimize.py' or list of those
+
+    ax: matplotlib.Axes, optional
+        Axes object to use.
+
+    size: tuple, optional
+        Figure size (width, height) in inches. Is only applied when no ax
+        object is specified
+
+    trace_x: str, optional
+        What should be plotted on the x-axis?
+        Possibilities: 'time', 'steps'
+        Default: 'steps'
+
+    trace_y: str, optional
+        What should be plotted on the y-axis?
+        Possibilities: 'fval', 'gradnorm', 'stepsize'
+        Default: 'fval'
+
+    scale_y: str, optional
+        May be logarithmic or linear ('log10' or 'lin')
+
+    offset_y: float, optional
+        Offset for the y-axis-values, as these are plotted on a log10-scale
+        Will be computed automatically if necessary
+
+    colors: list, or RGBA, optional
+        list of colors, or single color
+        color or list of colors for plotting. If not set, clustering is done
+        and colors are assigned automatically
+
+    y_limits: float or ndarray, optional
+        maximum value to be plotted on the y-axis, or y-limits
+
+    start_indices: list or int
+        list of integers specifying the multistart to be plotted or
+        int specifying up to which start index should be plotted
+
+    reference: list, optional
+        List of reference points for optimization results, containing et
+        least a function value fval
+
+    legends: list or str
+        Labels for line plots, one label per result object
+
+    Returns
+    -------
+
+    ax: matplotlib.Axes
+        The plot axes.
+    """
+
+    # parse input
+    (results, colors, legends) = process_result_list(results, colors, legends)
+
+    for j, result in enumerate(results):
+        # extract cost function values from result
+        (x_label, y_label, vals) = get_trace(result, trace_x, trace_y)
+
+        # compute the necessary offset for the y-axis
+        vals = get_vals(vals, scale_y, offset_y, start_indices)
+
+        # call lowlevel plot routine
+        ax = optimizer_history_lowlevel(vals, scale_y=scale_y, ax=ax,
+                                        colors=colors[j],
+                                        size=size, x_label=x_label,
+                                        y_label=y_label,
+                                        legend_text=legends[j])
+
+    # parse and apply plotting options
+    ref = create_references(references=reference)
+
+    # handle options
+    ax = handle_options(ax, vals, ref, y_limits)
+
+    return ax
+
+
+def optimizer_history_lowlevel(vals, scale_y='log10', colors=None, ax=None,
+                               size=(18.5, 10.5), x_label='Optimizer steps',
+                               y_label='Objective value', legend_text=None):
+    """
+    Plot optimizer history using list of numpy arrays.
+
+    Parameters
+    ----------
+
+    vals: list of numpy arrays
+        list of 2xn-arrays (x_values and y_values of the trace)
+
+    scale_y: str, optional
+        May be logarithmic or linear ('log10' or 'lin')
+
+    colors: list, or RGBA, optional
+        list of colors, or single color
+        color or list of colors for plotting. If not set, clustering is done
+        and colors are assigned automatically
+
+    ax: matplotlib.Axes, optional
+        Axes object to use.
+
+    size: tuple, optional
+        see waterfall
+
+    x_label: str
+        label for x-axis
+
+    y_label: str
+        label for y-axis
+
+    legend_text: str
+        Label for line plots
+
+    Returns
+    -------
+
+    ax: matplotlib.Axes
+        The plot axes.
+    """
+
+    # axes
+    if ax is None:
+        ax = plt.subplots()[1]
+        fig = plt.gcf()
+        fig.set_size_inches(*size)
+
+    # parse input
+    fvals = []
+    if isinstance(vals, list):
+        # convert entries to numpy arrays
+        for val in vals:
+            val = np.array(val)
+            fvals.append(val[1, -1])
+    else:
+        # convert to a list of numpy arrays
+        vals = np.array(vals)
+        if vals.shape[0] != 2 or vals.ndim != 2:
+            raise('If numpy array is passed directly to lowlevel routine of'
+                  'optimizer_history, shape needs to be 2 x n.')
+        fvals = [vals[1, -1]]
+        vals = [vals]
+    n_fvals = len(fvals)
+
+    # assign colors
+    # note: this has to happen before sorting
+    # to get the same colors in different plots
+    colors = assign_colors(fvals, colors)
+
+    # sort
+    indices = sorted(range(n_fvals), key=lambda j: fvals[j])
+
+    # plot
+    ax.xaxis.set_major_locator(MaxNLocator(integer=True))
+    for j, val in enumerate(vals):
+        # collect and parse data
+        j_fval = indices[j]
+        color = colors[j_fval]
+        if j == 0:
+            tmp_legend = legend_text
+        else:
+            tmp_legend = None
+
+        # line plots
+        if scale_y == 'log10':
+            ax.semilogy(val[0, :], val[1, :], color=color, label=tmp_legend)
+        else:
+            ax.plot(val[0, :], val[1, :], color=color, label=tmp_legend)
+
+    # set labels
+    ax.set_xlabel(x_label)
+    ax.set_ylabel(y_label)
+    ax.set_title('Optimizer history')
+
+    return ax
+
+
+def get_trace(result, trace_x, trace_y):
+    """
+    Get the values of the optimizer trace from the pypesto.Result object
+
+    Parameters
+    ----------
+
+    result: pypesto.Result
+        Optimization result obtained by 'optimize.py'.
+
+    trace_x: str, optional
+        What should be plotted on the x-axis?
+        Possibilities: 'time', 'steps'
+        Default: 'steps'
+
+    trace_y: str, optional
+        What should be plotted on the y-axis?
+        Possibilities: 'fval'(later also: 'gradnorm', 'stepsize')
+        Default: 'fval'
+
+    Returns
+    -------
+
+    ax: matplotlib.Axes
+        The plot axes.
+
+    x_label: str
+        label for x-axis to be plotted later
+
+    y_label: str
+        label for y-axis to be plotted later
+    """
+
+    # get data frames
+    traces = result.optimize_result.get_for_key('trace')
+    vals = []
+
+    for trace in traces:
+        # SciPy optimizers seem to save more steps than actually taken: prune
+        indices = np.argwhere(np.isfinite(trace['fval']))
+        indices = indices.flatten()
+        indices.astype(int)
+
+        # retrieve values from dataframe
+        if trace_x == 'time':
+            x_vals = np.array(trace['time'][indices])
+            x_label = 'Computation time [s]'
+
+        else:  # trace_x == 'steps':
+            x_vals = np.array(list(range(len(indices))))
+            x_label = 'Optimizer steps'
+
+        if trace_y == 'gradnorm':
+            # retrieve gradient trace, if saved
+            if trace['grad'] is None:
+                raise("No gradient norm trace can be visualized: "
+                      "The pypesto.result object does not contain "
+                      "a gradient trace")
+
+            # Get gradient trace, prune Nones, compute norm
+            tmp_grad_trace = list(trace['grad'].values)
+            y_vals = np.array([np.linalg.norm(grad) for grad in
+                               tmp_grad_trace if grad is not None])
+            y_label = 'Gradient norm'
+
+        else:  # trace_y == 'fval':
+            y_label = 'Objective value'
+            y_vals = np.array(trace['fval'][indices])
+
+        # write down values
+        vals.append(np.array([x_vals, y_vals]))
+
+    return x_label, y_label, vals
+
+
+def get_vals(vals, scale_y, offset_y, start_indices):
+    """
+    Postprocesses the values of the optimization history, depending on the
+    options set by the user (e.g. scale_y, offset_y, start_indices)
+
+    Parameters
+    ----------
+
+    vals: list
+        list of numpy arrays of dimension 2 x len(start_indices)
+
+    scale_y: str, optional
+        May be logarithmic or linear ('log10' or 'lin')
+
+    offset_y:
+        offset for the y-axis, as this is supposed to be in log10-scale
+
+    start_indices: list or int
+        list of integers specifying the multi start indices to be plotted or
+        int specifying up to which start index trajectories should be plotted
+
+    Returns
+    -------
+
+    vals: list
+        list of numpy arrays of size 2 x len(start_indices)
+    """
+
+    # get list of indices
+    if start_indices is None:
+        start_indices = np.array(range(len(vals)))
+    else:
+        # check whether list or maximum value
+        start_indices = np.array(start_indices)
+        if start_indices.size == 1:
+            start_indices = np.array(range(start_indices))
+
+        # check, whether index set is not too big
+        existing_indices = np.array(range(len(vals)))
+        start_indices = np.intersect1d(start_indices, existing_indices)
+
+    # reduce values to listed values
+    vals = [val for i, val in enumerate(vals) if i in start_indices]
+
+    # get the minimal value which should be plotted
+    min_val = np.inf
+    for val in vals:
+        tmp_min = np.min(val[1, :])
+        min_val = np.min([min_val, tmp_min])
+
+    # check, whether offset can be used with this data
+    offset_y = process_offset_y(offset_y, scale_y, min_val)
+
+    if offset_y != 0:
+        for val in vals:
+            val[1, :] += offset_y * np.ones(val[1].shape)
+
+    return vals
+
+
+def handle_options(ax, vals, ref, y_limits):
+    """
+    Get the limits for the y-axis, plots the reference points, will do
+    more at a later time point. This function is there to apply whatever
+    kind of post-plotting transformations to the axis object.
+
+    Parameters
+    ----------
+
+    ref: list, optional
+        List of reference points for optimization results, containing et
+        least a function value fval
+
+    vals: list
+        list of numpy arrays of size 2 x number of values
+
+    ax: matplotlib.Axes, optional
+        Axes object to use.
+
+    y_limits: float or ndarray, optional
+        maximum value to be plotted on the y-axis, or y-limits
+
+    Returns
+    -------
+
+    ax: matplotlib.Axes
+        The plot axes.
+    """
+
+    # handle y-limits
+    ax = process_y_limits(ax, y_limits)
+
+    # handle reference points
+    if len(ref) > 0:
+        # plot reference points
+        # get length of longest trajectory
+        max_len = 0
+        for val in vals:
+            max_len = np.max([max_len, val[0, -1]])
+
+        for i_ref in ref:
+            ax.semilogy([0, max_len], [i_ref.fval, i_ref.fval], '--',
+                        color=i_ref.color, label=i_ref.legend)
+
+    return ax
diff --git a/pypesto/visualize/parameters.py b/pypesto/visualize/parameters.py
index 29ba2b3..ae4a264 100644
--- a/pypesto/visualize/parameters.py
+++ b/pypesto/visualize/parameters.py
@@ -1,19 +1,21 @@
 import matplotlib.pyplot as plt
 from matplotlib.ticker import MaxNLocator
 import numpy as np
-from .clust_color import assign_clustered_colors
+from .reference_points import create_references
+from .clust_color import assign_colors
+from .misc import process_result_list
 
 
-def parameters(result, ax=None, free_indices_only=True, lb=None, ub=None,
-               size=None):
+def parameters(results, ax=None, free_indices_only=True, lb=None, ub=None,
+               size=None, reference=None, colors=None, legends=None):
     """
     Plot parameter values.
 
     Parameters
     ----------
 
-    result: pypesto.Result
-        Optimization result obtained by 'optimize.py'.
+    results: pypesto.Result or list
+        Optimization result obtained by 'optimize.py' or list of those
 
     ax: matplotlib.Axes, optional
         Axes object to use.
@@ -30,6 +32,18 @@ def parameters(result, ax=None, free_indices_only=True, lb=None, ub=None,
         Figure size (width, height) in inches. Is only applied when no ax
         object is specified
 
+    reference: list, optional
+        List of reference points for optimization results, containing et
+        least a function value fval
+
+    colors: list, or RGBA, optional
+        list of colors, or single color
+        color or list of colors for plotting. If not set, clustering is done
+        and colors are assigned automatically
+
+    legends: list or str
+        Labels for line plots, one label per result object
+
     Returns
     -------
 
@@ -37,32 +51,34 @@ def parameters(result, ax=None, free_indices_only=True, lb=None, ub=None,
         The plot axes.
     """
 
-    if lb is None:
-        lb = result.problem.lb
-    if ub is None:
-        ub = result.problem.ub
+    # parse input
+    (results, colors, legends) = process_result_list(results, colors, legends)
 
-    fvals = result.optimize_result.get_for_key('fval')
-    xs = result.optimize_result.get_for_key('x')
+    for j, result in enumerate(results):
+        # handle results and bounds
+        (lb, ub, x_labels, fvals, xs) = \
+            handle_inputs(result=result, lb=lb, ub=ub,
+                          free_indices_only=free_indices_only)
 
-    x_labels = result.problem.x_names
+        # call lowlevel routine
+        ax = parameters_lowlevel(xs=xs, fvals=fvals, lb=lb, ub=ub,
+                                 x_labels=x_labels, ax=ax, size=size,
+                                 colors=colors[j], legend_text=legends[j])
 
-    if free_indices_only:
-        for ix, x in enumerate(xs):
-            xs[ix] = result.problem.get_reduced_vector(x)
-        lb = result.problem.get_reduced_vector(lb)
-        ub = result.problem.get_reduced_vector(ub)
-        x_labels = [x_labels[int(i)] for i in result.problem.x_free_indices]
-    else:
-        lb = result.problem.get_full_vector(lb)
-        ub = result.problem.get_full_vector(ub)
+    # parse and apply plotting options
+    ref = create_references(references=reference)
 
-    return parameters_lowlevel(xs=xs, fvals=fvals, lb=lb, ub=ub,
-                               x_labels=x_labels, ax=ax, size=size)
+    # plot reference points
+    for i_ref in ref:
+        ax = parameters_lowlevel([i_ref['x']], [i_ref['fval']], ax=ax,
+                                 colors=i_ref['color'],
+                                 legend_text=i_ref.legend)
+
+    return ax
 
 
 def parameters_lowlevel(xs, fvals, lb=None, ub=None, x_labels=None, ax=None,
-                        size=None):
+                        size=None, colors=None, legend_text=None):
     """
     Plot parameters plot using list of parameters.
 
@@ -88,6 +104,12 @@ def parameters_lowlevel(xs, fvals, lb=None, ub=None, x_labels=None, ax=None,
     size: tuple, optional
         see parameters
 
+    colors: list of RGBA
+        One for each element in 'fvals'.
+
+    legend_text: str
+        Label for line plots
+
     Returns
     -------
 
@@ -101,23 +123,30 @@ def parameters_lowlevel(xs, fvals, lb=None, ub=None, x_labels=None, ax=None,
 
     if size is None:
         # 0.5 inch height per parameter
-        size = (18.5, xs.shape[1]/2)
+        size = (18.5, xs.shape[1] / 2)
 
     if ax is None:
         ax = plt.subplots()[1]
         fig = plt.gcf()
         fig.set_size_inches(*size)
 
-    # assign color
-    colors = assign_clustered_colors(fvals)
+    # assign colors
+    colors = assign_colors(vals=fvals, colors=colors)
 
     # parameter indices
-    parameters_ind = range(1, xs.shape[1] + 1)
+    parameters_ind = list(range(1, xs.shape[1] + 1))[::-1]
 
     # plot parameters
     ax.xaxis.set_major_locator(MaxNLocator(integer=True))
     for j_x, x in reversed(list(enumerate(xs))):
-        ax.plot(x, parameters_ind, color=colors[j_x], marker='o')
+        if j_x == 0:
+            tmp_legend = legend_text
+        else:
+            tmp_legend = None
+        ax.plot(x, parameters_ind,
+                color=colors[j_x],
+                marker='o',
+                label=tmp_legend)
 
     plt.yticks(parameters_ind, x_labels)
 
@@ -133,3 +162,65 @@ def parameters_lowlevel(xs, fvals, lb=None, ub=None, x_labels=None, ax=None,
     ax.set_title('Estimated parameters')
 
     return ax
+
+
+def handle_inputs(result, free_indices_only, lb=None, ub=None):
+    """
+    Computes the correct bounds for the parameter indices to be plotted and
+    outputs the corrsponding parameters and their labels
+
+    Parameters
+    ----------
+
+    result: pypesto.Result
+        Optimization result obtained by 'optimize.py'.
+
+    free_indices_only: bool, optional
+        If True, only free parameters are shown. If
+        False, also the fixed parameters are shown.
+
+    lb, ub: ndarray, optional
+        If not None, override result.problem.lb, problem.problem.ub.
+        Dimension either result.problem.dim or result.problem.dim_full.
+
+    Returns
+    -------
+
+    lb, ub: ndarray
+        Dimension either result.problem.dim or result.problem.dim_full.
+
+    x_labels: list of str
+        ytick labels to be applied later on
+
+    fvals: ndarray
+        objective function values which are needed for plotting later
+
+    xs: ndarray
+        parameter values which will be plotted later
+    """
+
+    # retrieve results
+    fvals = result.optimize_result.get_for_key('fval')
+    xs = result.optimize_result.get_for_key('x')
+
+    # get bounds
+    if lb is None:
+        lb = result.problem.lb
+    if ub is None:
+        ub = result.problem.ub
+
+    # get labels
+    x_labels = result.problem.x_names
+
+    # handle fixed and free indices
+    if free_indices_only:
+        for ix, x in enumerate(xs):
+            xs[ix] = result.problem.get_reduced_vector(x)
+        lb = result.problem.get_reduced_vector(lb)
+        ub = result.problem.get_reduced_vector(ub)
+        x_labels = [x_labels[int(i)] for i in result.problem.x_free_indices]
+    else:
+        lb = result.problem.get_full_vector(lb)
+        ub = result.problem.get_full_vector(ub)
+
+    return lb, ub, x_labels, fvals, xs
diff --git a/pypesto/visualize/profiles.py b/pypesto/visualize/profiles.py
index cd9db96..c9812b7 100644
--- a/pypesto/visualize/profiles.py
+++ b/pypesto/visualize/profiles.py
@@ -1,9 +1,13 @@
 import matplotlib.pyplot as plt
 from matplotlib.ticker import MaxNLocator
 import numpy as np
+from .reference_points import create_references
+from .clust_color import assign_colors
+from .misc import process_result_list
 
 
-def profiles(result, ax=None, profile_indices=None, size=(18.5, 6.5)):
+def profiles(results, fig=None, profile_indices=None, size=(18.5, 6.5),
+             reference=None, colors=None, legends=None):
     """
     Plot classical 1D profile plot (using the posterior, e.g. Gaussian like
     profile)
@@ -11,11 +15,11 @@ def profiles(result, ax=None, profile_indices=None, size=(18.5, 6.5)):
     Parameters
     ----------
 
-    result: pypesto.Result
-        Optimization result obtained by 'optimize.py'
+    results: list or pypesto.Result
+        list of pypesto.Result or single pypesto.Result
 
-    ax: matplotlib.Axes, optional
-        Axes object to use.
+    fig: matplotlib.Figure, optional
+        Figure object to use.
 
     profile_indices: list of integer values
         list of integer values specifying which profiles should be plotted
@@ -24,6 +28,18 @@ def profiles(result, ax=None, profile_indices=None, size=(18.5, 6.5)):
         Figure size (width, height) in inches. Is only applied when no ax
         object is specified
 
+    reference: list, optional
+        List of reference points for optimization results, containing et
+        least a function value fval
+
+    colors: list, or RGBA, optional
+        list of colors, or single color
+        color or list of colors for plotting. If not set, clustering is done
+        and colors are assigned automatically
+
+    legends: list or str
+        Labels for line plots, one label per result object
+
     Returns
     -------
 
@@ -31,25 +47,37 @@ def profiles(result, ax=None, profile_indices=None, size=(18.5, 6.5)):
         The plot axes.
     """
 
+    # parse input
+    (results, colors, legends) = process_result_list(results, colors, legends)
+
+    # get the correct number of parameter indices, even if not the same in
+    # all result obejcts
     if profile_indices is None:
-        profile_indices = \
-            [i for i in range(0, len(result.profile_result.list[0]))]
+        profile_indices = []
+        for result in results:
+            tmp_indices = [i for i in
+                           range(len(result.profile_result.list[0]))]
+            profile_indices = list(set().union(profile_indices, tmp_indices))
 
-    # extract ratio values values from result
-    fvals = []
-    for i_par in range(0, len(result.profile_result.list[0])):
-        if i_par in profile_indices:
-            tmp = np.array(
-                [result.profile_result.list[0][i_par].x_path[i_par, :],
-                 result.profile_result.list[0][i_par].ratio_path[:]])
-        else:
-            tmp = None
-        fvals.append(tmp)
+    # loop over results
+    for j, result in enumerate(results):
+        fvals = handle_inputs(result, profile_indices)
+
+        # call lowlevel routine
+        ax = profiles_lowlevel(fvals=fvals, ax=fig, size=size,
+                               color=colors[j], legend_text=legends[j])
+
+    # parse and apply plotting options
+    ref = create_references(references=reference)
 
-    return profiles_lowlevel(fvals, ax)
+    # plot reference points
+    ax = handle_reference_points(ref, ax, fvals)
+
+    return ax
 
 
-def profiles_lowlevel(fvals, ax=None, size=(18.5, 6.5)):
+def profiles_lowlevel(fvals, ax=None, size=(18.5, 6.5), color=None,
+                      legend_text=None):
     """
     Lowlevel routine for profile plotting, working with a list of arrays
     only, opening different axes objects in case
@@ -60,13 +88,23 @@ def profiles_lowlevel(fvals, ax=None, size=(18.5, 6.5)):
     fvals: numeric list or array
         Including values need to be plotted.
 
-    ax: matplotlib.Axes, optional
-        Axes object to use.
+    ax: list of matplotlib.Axes, optional
+        list of axes object to use.
+
+    size: tuple, optional
+        Figure size (width, height) in inches. Is only applied when no ax
+        object is specified
 
     size: tuple, optional
         Figure size (width, height) in inches. Is only applied when no ax
         object is specified
 
+    color: RGBA, optional
+        color for profiles in plot.
+
+    legend_text: str
+        Label for line plots
+
     Returns
     -------
 
@@ -76,9 +114,12 @@ def profiles_lowlevel(fvals, ax=None, size=(18.5, 6.5)):
 
     # axes
     if ax is None:
-        ax = plt.subplots()[1]
-        fig = plt.gcf()
+        ax = []
+        fig = plt.figure()
         fig.set_size_inches(*size)
+    else:
+        plt.axes(ax)
+        fig = plt.gcf()
 
     if isinstance(fvals, list):
         n_fvals = 0
@@ -95,23 +136,35 @@ def profiles_lowlevel(fvals, ax=None, size=(18.5, 6.5)):
     else:
         rows = columns - 1
 
-    counter = 1
+    counter = 0
     for i_plot, fval in enumerate(fvals):
+        # handle legend
+        if i_plot == 0:
+            tmp_legend = legend_text
+        else:
+            tmp_legend = None
 
+        # plot if data
         if fval is not None:
-            ax = plt.subplot(rows, columns, counter)
-            ax = profile_lowlevel(fval, ax)
+            ax.append(fig.add_subplot(rows, columns, counter + 1))
+            ax[counter] = profile_lowlevel(fval, ax[counter],
+                                           size=size, color=color,
+                                           legend_text=tmp_legend)
 
             # labels
-            ax.set_xlabel(f'Parameter {i_plot} value')
-            if (counter-1) % columns == 0:
-                ax.set_ylabel('Log-posterior ratio')
+            ax[counter].set_xlabel(f'Parameter {i_plot} value')
+            if counter % columns == 0:
+                ax[counter].set_ylabel('Log-posterior ratio')
             else:
-                ax.set_yticklabels([''])
+                ax[counter].set_yticklabels([''])
             counter += 1
+            tmp_legend = None
 
+    return ax
 
-def profile_lowlevel(fvals, ax=None, size=(18.5, 6.5)):
+
+def profile_lowlevel(fvals, ax=None, size=(18.5, 6.5), color=None,
+                     legend_text=None):
     """
     Lowlevel routine for plotting one profile, working with a numpy array only
 
@@ -128,6 +181,12 @@ def profile_lowlevel(fvals, ax=None, size=(18.5, 6.5)):
         Figure size (width, height) in inches. Is only applied when no ax
         object is specified
 
+    color: RGBA, optional
+        color for profiles in plot.
+
+    legend_text: str
+        Label for line plots
+
     Returns
     -------
 
@@ -138,6 +197,9 @@ def profile_lowlevel(fvals, ax=None, size=(18.5, 6.5)):
     # parse input
     fvals = np.array(fvals)
 
+    # get colors
+    color = assign_colors([1.], color)
+
     # axes
     if ax is None:
         ax = plt.subplots()[1]
@@ -149,6 +211,73 @@ def profile_lowlevel(fvals, ax=None, size=(18.5, 6.5)):
     # plot
     if fvals.size != 0:
         ax.xaxis.set_major_locator(MaxNLocator(integer=True))
-        ax.plot(fvals[0, :], fvals[1, :], color=[.9, .2, .2, 1.])
+        ax.plot(fvals[0, :], fvals[1, :], color=color[0], label=legend_text)
 
     return ax
+
+
+def handle_reference_points(ref, ax, fvals):
+    """
+    Handle reference points.
+
+    Parameters
+    ----------
+
+    ref: list, optional
+        List of reference points for optimization results, containing et
+        least a function value fval
+
+    ax: matplotlib.Axes, optional
+        Axes object to use.
+    """
+
+    if len(ref) > 0:
+        # get the parameters which have profiles plotted
+        par_indices = []
+        for i_plot, fval in enumerate(fvals):
+            if fval is not None:
+                par_indices.append(i_plot)
+
+        # loop over axes objects
+        for i_par, i_ax in enumerate(ax):
+            for i_ref in ref:
+                current_x = i_ref['x'][par_indices[i_par]]
+                i_ax.plot([current_x, current_x], [0., 1.],
+                          color=i_ref.color, label=i_ref.legend)
+
+    return ax
+
+
+def handle_inputs(result, profile_indices):
+    """
+    Retrieves the values of the profiles to be plotted later from a
+    pypesto.ProfileResult object
+
+    Parameters
+    ----------
+
+    result: pypesto.Result
+        Profile result obtained by 'profile.py'
+
+    profile_indices: list of integer values
+        list of integer values specifying which profiles should be plotted
+
+    Returns
+    -------
+
+    fvals: numeric list
+        Including values need to be plotted.
+    """
+
+    # extract ratio values values from result
+    fvals = []
+    for i_par in range(0, len(result.profile_result.list[0])):
+        if i_par in profile_indices:
+            tmp = np.array(
+                [result.profile_result.list[0][i_par].x_path[i_par, :],
+                 result.profile_result.list[0][i_par].ratio_path[:]])
+        else:
+            tmp = None
+        fvals.append(tmp)
+
+    return fvals
diff --git a/pypesto/visualize/reference_points.py b/pypesto/visualize/reference_points.py
new file mode 100644
index 0000000..38297e6
--- /dev/null
+++ b/pypesto/visualize/reference_points.py
@@ -0,0 +1,170 @@
+import numpy as np
+
+
+class ReferencePoint(dict):
+    """
+    Reference point for plotting. Should contain a parameter value and an
+    objective function value, may alos contain a color and a legend.
+
+    Can be used like a dict.
+
+    Attributes
+    ----------
+
+    x: ndarray
+        Reference parameters.
+
+    fval: float
+        Function value, fun(x), for reference parameters.
+
+    color: RGBA, optional
+        Color which should be used for reference point.
+
+    auto_color: boolean
+        flag indicating whether color for this reference point should be
+        assigned automatically or whether it was assigned by user
+
+    legend: str
+        legend text for reference point
+    """
+
+    def __init__(self,
+                 reference=None,
+                 x=None,
+                 fval=None,
+                 color=None,
+                 legend=None):
+        super().__init__()
+
+        if (reference is not None) and ((x is not None) or (fval is not None)):
+            raise ("Please specify either an argument for reference or for x "
+                   "and fval, but not both.")
+
+        # assign legend, my be None
+        self.legend = legend
+        if isinstance(reference, dict) or \
+                isinstance(reference, ReferencePoint):
+            # Handle case of dict or ReferencePoint
+            self.x = np.array(reference["x"])
+            self.fval = reference["fval"]
+            if "color" in reference.keys():
+                self.color = reference["color"]
+                if "auto_color" in reference.keys():
+                    self.auto_color = reference["auto_color"]
+                else:
+                    self.auto_color = True
+            else:
+                self.color = None
+                self.auto_color = True
+        elif isinstance(reference, tuple):
+            # Handle case of tuple
+            self.x = np.array(reference[0])
+            self.fval = reference[1]
+            if len(reference) > 2:
+                self.color = reference[2]
+                self.auto_color = False
+            else:
+                self.color = None
+                self.auto_color = True
+        if reference is None:
+            if x is not None:
+                self.x = np.array(x)
+            else:
+                raise ('Parameter vector x not passed, but is a manadatory '
+                       'input when creating a reference point. Stopping.')
+            if fval is not None:
+                self.fval = fval
+            else:
+                raise ('Objective value fval not passed, but is a manadatory '
+                       'input when creating a reference point. Stopping.')
+            if color is not None:
+                self.color = color
+                self.auto_color = False
+            else:
+                self.color = None
+                self.auto_color = True
+
+    def __getattr__(self, key):
+        try:
+            return self[key]
+        except KeyError:
+            raise AttributeError(key)
+
+    __setattr__ = dict.__setitem__
+    __delattr__ = dict.__delitem__
+
+
+def assign_colors(ref):
+    """
+    Assigns colors to reference points, depending on user settings
+
+    Parameters
+    ----------
+
+    ref: list of ReferencePoint
+        Reference points, which need to get their color property filled
+
+    Returns
+    -------
+
+    ref: list of ReferencePoint
+        Reference points, which got their color property filled
+    """
+
+    # loop over reference points
+    auto_color_count = 0
+    for i_ref in ref:
+        if i_ref['auto_color']:
+            auto_color_count += 1
+
+    auto_colors = [[0., 0.5 * (1. + i_auto / auto_color_count), 0., 0.9] for
+                   i_auto in range(auto_color_count)]
+
+    # loop over reference points and assign auto_colors
+    auto_color_count = 0
+    for i_num, i_ref in enumerate(ref):
+        if i_ref['auto_color']:
+            i_ref['color'] = auto_colors[i_num]
+            auto_color_count += 1
+
+    return ref
+
+
+def create_references(references=None, x=None, fval=None):
+    """
+    This function creates a list of reference point objects from user inputs
+
+    Parameters
+    ----------
+
+    references: ReferencePoint or dict or list, optional
+        Will be converted into a list of RefPoints
+
+    x: ndarray, optional
+        Parameter vector which should be used for reference point
+
+    fval: float, optional
+        Objective function value which should be used for reference point
+
+    Returns
+    -------
+
+    colors: list of RGBA
+        One for each element in 'vals'.
+    """
+
+    # parse input (reference)
+    ref = []
+    if references is not None:
+        if isinstance(references, list):
+            for reference in references:
+                ref.append(ReferencePoint(reference))
+        else:
+            ref.append(ReferencePoint(references))
+
+    # parse input (x and fval)
+    if (x is not None) and (fval is not None):
+        ref.append(ReferencePoint(x=x, fval=fval))
+
+    # assign colors for reference points which have no user-specified colors
+    return assign_colors(ref)
diff --git a/pypesto/visualize/waterfall.py b/pypesto/visualize/waterfall.py
index 550d03a..55de666 100644
--- a/pypesto/visualize/waterfall.py
+++ b/pypesto/visualize/waterfall.py
@@ -1,18 +1,31 @@
 import matplotlib.pyplot as plt
 from matplotlib.ticker import MaxNLocator
 import numpy as np
-from .clust_color import assign_clustered_colors
+from .reference_points import create_references
+from .clust_color import assign_colors
+from .misc import process_result_list
+from .misc import process_y_limits
+from .misc import process_offset_y
 
 
-def waterfall(result, ax=None, size=(18.5, 10.5)):
+def waterfall(results,
+              ax=None,
+              size=(18.5, 10.5),
+              y_limits=None,
+              scale_y='log10',
+              offset_y=None,
+              start_indices=None,
+              reference=None,
+              colors=None,
+              legends=None):
     """
     Plot waterfall plot.
 
     Parameters
     ----------
 
-    result: pypesto.Result
-        Optimization result obtained by 'optimize.py'
+    results: pypesto.Result or list
+        Optimization result obtained by 'optimize.py' or list of those
 
     ax: matplotlib.Axes, optional
         Axes object to use.
@@ -21,6 +34,31 @@ def waterfall(result, ax=None, size=(18.5, 10.5)):
         Figure size (width, height) in inches. Is only applied when no ax
         object is specified
 
+    y_limits: float or ndarray, optional
+        maximum value to be plotted on the y-axis, or y-limits
+
+    scale_y: str, optional
+        May be logarithmic or linear ('log10' or 'lin')
+
+    offset_y:
+        offset for the y-axis, if this is supposed to be in log10-scale
+
+    start_indices: list or int
+        list of integers specifying the multistart to be plotted or
+        int specifying up to which start index should be plotted
+
+    reference: list, optional
+        List of reference points for optimization results, containing et
+        least a function value fval
+
+    colors: list, or RGBA, optional
+        list of colors, or single color
+        color or list of colors for plotting. If not set, clustering is done
+        and colors are assigned automatically
+
+    legends: list or str
+        Labels for line plots, one label per result object
+
     Returns
     -------
 
@@ -28,13 +66,33 @@ def waterfall(result, ax=None, size=(18.5, 10.5)):
         The plot axes.
     """
 
-    # extract cost function values from result
-    fvals = result.optimize_result.get_for_key('fval')
+    # parse input
+    (results, colors, legends) = process_result_list(results, colors, legends)
+
+    # plotting routine needs the maximum number of multistarts
+    max_len_fvals = np.array([0])
+
+    # loop over results
+    for j, result in enumerate(results):
+        # extract specific cost function values from result
+        fvals = get_fvals(result, scale_y, offset_y, start_indices)
+        max_len_fvals = np.max([max_len_fvals, len(fvals)])
+
+        # call lowlevel plot routine
+        ax = waterfall_lowlevel(fvals=fvals, scale_y=scale_y, ax=ax, size=size,
+                                colors=colors[j], legend_text=legends[j])
 
-    return waterfall_lowlevel(fvals, ax, size)
+    # parse and apply plotting options
+    ref = create_references(references=reference)
 
+    # apply changes specified be the user to the axis object
+    ax = handle_options(ax, max_len_fvals, ref, y_limits)
 
-def waterfall_lowlevel(fvals, ax=None, size=(18.5, 10.5)):
+    return ax
+
+
+def waterfall_lowlevel(fvals, scale_y='log10', ax=None, size=(18.5, 10.5),
+                       colors=None, legend_text=None):
     """
     Plot waterfall plot using list of function values.
 
@@ -44,12 +102,23 @@ def waterfall_lowlevel(fvals, ax=None, size=(18.5, 10.5)):
     fvals: numeric list or array
         Including values need to be plotted.
 
+    scale_y: str, optional
+        May be logarithmic or linear ('log10' or 'lin')
+
     ax: matplotlib.Axes, optional
         Axes object to use.
 
     size: tuple, optional
         see waterfall
 
+    colors: list, or RGBA, optional
+        list of colors, or single color
+        color or list of colors for plotting. If not set, clustering is done
+        and colors are assigned automatically
+
+    legend_text: str
+        Label for line plots
+
     Returns
     -------
 
@@ -72,7 +141,7 @@ def waterfall_lowlevel(fvals, ax=None, size=(18.5, 10.5)):
     # assign colors
     # note: this has to happen before sorting
     # to get the same colors in different plots
-    colors = assign_clustered_colors(fvals)
+    colors = assign_colors(fvals, colors=colors)
 
     # sort
     indices = sorted(range(n_fvals),
@@ -80,16 +149,133 @@ def waterfall_lowlevel(fvals, ax=None, size=(18.5, 10.5)):
 
     # plot
     ax.xaxis.set_major_locator(MaxNLocator(integer=True))
-    ax.plot(start_ind, fvals)
+    # plot line
+    if scale_y == 'log10':
+        ax.semilogy(start_ind, fvals)
+    else:
+        ax.plot(start_ind, fvals)
+
+    # plot points
     for j in range(n_fvals):
+        # parse data for plotting
         j_fval = indices[j]
         color = colors[j_fval]
         fval = fvals[j_fval]
-        ax.plot(j, fval, color=color, marker='o')
+        if j == 0:
+            tmp_legend = legend_text
+        else:
+            tmp_legend = None
+
+        # line plot (linear or logarithmic)
+        if scale_y == 'log10':
+            ax.semilogy(j, fval, color=color, marker='o', label=tmp_legend)
+        else:
+            ax.plot(j, fval, color=color, marker='o', label=tmp_legend)
 
     # labels
     ax.set_xlabel('Ordered optimizer run')
     ax.set_ylabel('Function value')
     ax.set_title('Waterfall plot')
+    ax.legend()
+
+    return ax
+
+
+def get_fvals(result, scale_y, offset_y, start_indices):
+    """
+    Get function values to be plotted later from results.
+
+    Parameters
+    ----------
+
+    result: pypesto.Result
+        Optimization result obtained by 'optimize.py'
+
+    scale_y: str, optional
+        May be logarithmic or linear ('log10' or 'lin')
+
+    offset_y:
+        offset for the y-axis, if this is supposed to be in log10-scale
+
+    start_indices: list or int
+        list of integers specifying the multistart to be plotted or
+        int specifying up to which start index should be plotted
+
+    Returns
+    -------
+
+    fvals: ndarray
+        function values
+    """
+
+    # extract cost function values from result
+    fvals = np.array(result.optimize_result.get_for_key('fval'))
+
+    # get list of indices
+    if start_indices is None:
+        start_indices = np.array(range(len(fvals)))
+    else:
+        # check whether list or maximum value
+        start_indices = np.array(start_indices)
+        if start_indices.size == 1:
+            start_indices = np.array(range(start_indices))
+
+        # check, whether index set is not too big
+        existing_indices = np.array(range(len(fvals)))
+        start_indices = np.intersect1d(start_indices, existing_indices)
+
+    # reduce to indices for which the user asked
+    fvals = fvals[start_indices]
+
+    # get the minimal value which should be plotted
+    min_val = np.min(fvals)
+
+    # check, whether offset can be used with this data
+    offset_y = process_offset_y(offset_y, scale_y, min_val)
+
+    # apply offset
+    if offset_y != 0.:
+        fvals += offset_y * np.ones(fvals.shape)
+
+    # get only the indices which the user asked for
+    return fvals
+
+
+def handle_options(ax, max_len_fvals, ref, y_limits):
+    """
+    Get the limits for the y-axis, plots the reference points, will do
+    more at a later time point. This function is there to apply whatever
+    kind of post-plotting transformations to the axis object.
+
+    Parameters
+    ----------
+
+    ax: matplotlib.Axes, optional
+        Axes object to use.
+
+    max_len_fvals: int
+        maximum number of points
+
+    ref: list, optional
+        List of reference points for optimization results, containing at
+        least a function value fval
+
+    y_limits: float or ndarray, optional
+        maximum value to be plotted on the y-axis, or y-limits
+
+    Returns
+    -------
+
+    ax: matplotlib.Axes
+        The plot axes.
+    """
+
+    # handle y-limits
+    ax = process_y_limits(ax, y_limits)
+
+    # handle reference points
+    for i_ref in ref:
+        ax.plot([0, max_len_fvals - 1], [i_ref.fval, i_ref.fval], '--',
+                color=i_ref.color, label=i_ref.legend)
 
     return ax
diff --git a/setup.py b/setup.py
index ceeab12..588701b 100644
--- a/setup.py
+++ b/setup.py
@@ -28,7 +28,7 @@ setup(name='pypesto',
                         'pandas>=0.23.4',
                         'matplotlib>=2.2.3',
                         'cloudpickle>=0.7.0',
-                        'petab>=0.0.0a10'],
+                        'petab==0.0.0a10'],
       tests_require=['pytest', 'flake8>=3.7.1', 'gitpython'],
       extras_require={'amici': ['amici>=0.9.2']},
       python_requires='>=3.6',


Test patch:
diff --git a/test/test_visualize.py b/test/test_visualize.py
index d5f0159..1f11adc 100644
--- a/test/test_visualize.py
+++ b/test/test_visualize.py
@@ -5,48 +5,151 @@ import scipy as sp
 import unittest
 
 
-# define a pypesto problem
-objective = pypesto.Objective(fun=sp.optimize.rosen,
-                              grad=sp.optimize.rosen_der,
-                              hess=sp.optimize.rosen_hess)
-lb = -7 * np.ones((1, 2))
-ub = 7 * np.ones((1, 2))
-problem = pypesto.Problem(objective=objective, lb=lb, ub=ub)
-
-# write some dummy results for optimization
-result = pypesto.Result(problem=problem)
-for j in range(0, 3):
-    optimizer_result = pypesto.OptimizerResult(fval=j * 0.01,
-                                               x=[j + 0.1, j + 1])
-    result.optimize_result.append(optimizer_result=optimizer_result)
-for j in range(0, 4):
-    optimizer_result = pypesto.OptimizerResult(fval=10 + j * 0.01,
-                                               x=[2.5 + j + 0.1, 2 + j + 1])
-    result.optimize_result.append(optimizer_result=optimizer_result)
-
-# write some dummy results for profiling
-ratio_path_1 = [0.15, 0.25, 0.7, 1., 0.8, 0.35, 0.15]
-ratio_path_2 = [0.1, 0.2, 0.7, 1., 0.8, 0.3, 0.1]
-x_path_1 = np.array([[2., 2.1, 2.3, 2.5, 2.7, 2.9, 3.],
-                     [1., 1.2, 1.4, 1.5, 1.6, 1.8, 2.]])
-x_path_2 = np.array([[1., 1.1, 1.3, 1.5, 1.7, 1.9, 2.1],
-                     [2.1, 2.2, 2.4, 2.5, 2.8, 2.9, 3.1]])
-fval_path_1 = [4., 3., 1., 0., 1.5, 2.5, 5.]
-fval_path_2 = [4.5, 3.5, 1.5, 0., 1.3, 2.3, 4.3]
-tmp_result_1 = pypesto.ProfilerResult(x_path_1, fval_path_1, ratio_path_1)
-tmp_result_2 = pypesto.ProfilerResult(x_path_2, fval_path_2, ratio_path_2)
-
-# use pypesto function to write the numeric values into the results
-result.profile_result.create_new_profile_list()
-result.profile_result.create_new_profile(tmp_result_1)
-result.profile_result.create_new_profile(tmp_result_2)
+# Define some helper functions, to have the test code more readable
+def create_bounds():
+    # define bounds for a pypesto problem
+    lb = -7 * np.ones((1, 2))
+    ub = 7 * np.ones((1, 2))
+
+    return lb, ub
+
+
+def create_problem():
+    # define a pypesto objective (with tracing options)
+    objective_options = pypesto.ObjectiveOptions(trace_record=True,
+                                                 trace_save_iter=1)
+    objective = pypesto.Objective(fun=sp.optimize.rosen,
+                                  grad=sp.optimize.rosen_der,
+                                  hess=sp.optimize.rosen_hess,
+                                  options=objective_options)
+
+    # define a pypesto problem
+    (lb, ub) = create_bounds()
+    problem = pypesto.Problem(objective=objective, lb=lb, ub=ub)
+
+    return problem
+
+
+def create_optimization_result():
+    # create the pypesto problem
+    problem = create_problem()
+
+    # write some dummy results for optimization
+    result = pypesto.Result(problem=problem)
+    for j in range(0, 3):
+        optimizer_result = pypesto.OptimizerResult(fval=j * 0.01,
+                                                   x=[j + 0.1, j + 1])
+        result.optimize_result.append(optimizer_result=optimizer_result)
+    for j in range(0, 4):
+        optimizer_result = pypesto.OptimizerResult(fval=10 + j * 0.01,
+                                                   x=[2.5 + j + 0.1,
+                                                      2 + j + 1])
+        result.optimize_result.append(optimizer_result=optimizer_result)
+
+    return result
+
+
+def create_optimization_history():
+    # create the pypesto problem
+    problem = create_problem()
+
+    # create optimizer
+    optimizer_options = {'maxiter': 200}
+    optimizer = pypesto.ScipyOptimizer(method='TNC', options=optimizer_options)
+
+    # run optimization
+    optimize_options = pypesto.OptimizeOptions(allow_failed_starts=True)
+    result_with_trace = pypesto.minimize(
+        problem=problem,
+        optimizer=optimizer,
+        n_starts=5,
+        startpoint_method=pypesto.startpoint.uniform,
+        options=optimize_options
+    )
+
+    return result_with_trace
+
+
+def create_profile_result():
+    # create a pypesto result
+    result = create_optimization_result()
+
+    # write some dummy results for profiling
+    ratio_path_1 = [0.15, 0.25, 0.7, 1., 0.8, 0.35, 0.15]
+    ratio_path_2 = [0.1, 0.2, 0.7, 1., 0.8, 0.3, 0.1]
+    x_path_1 = np.array([[2., 2.1, 2.3, 2.5, 2.7, 2.9, 3.],
+                         [1., 1.2, 1.4, 1.5, 1.6, 1.8, 2.]])
+    x_path_2 = np.array([[1., 1.1, 1.3, 1.5, 1.7, 1.9, 2.1],
+                         [2.1, 2.2, 2.4, 2.5, 2.8, 2.9, 3.1]])
+    fval_path_1 = [4., 3., 1., 0., 1.5, 2.5, 5.]
+    fval_path_2 = [4.5, 3.5, 1.5, 0., 1.3, 2.3, 4.3]
+    tmp_result_1 = pypesto.ProfilerResult(x_path_1, fval_path_1, ratio_path_1)
+    tmp_result_2 = pypesto.ProfilerResult(x_path_2, fval_path_2, ratio_path_2)
+
+    # use pypesto function to write the numeric values into the results
+    result.profile_result.create_new_profile_list()
+    result.profile_result.create_new_profile(tmp_result_1)
+    result.profile_result.create_new_profile(tmp_result_2)
+
+    return result
+
+
+def create_plotting_options():
+    # create sets of reference points (from tuple, dict and from list)
+    ref1 = ([1., 1.5], 0.2)
+    ref2 = ([1.8, 1.9], 0.6)
+    ref3 = {'x': np.array([1.4, 1.7]), 'fval': 0.4}
+    ref4 = [ref1, ref2]
+    ref_point = pypesto.visualize.create_references(ref4)
+
+    return ref1, ref2, ref3, ref4, ref_point
 
 
 class TestVisualize(unittest.TestCase):
 
     @staticmethod
     def test_waterfall():
-        pypesto.visualize.waterfall(result)
+        # create the necessary results
+        result_1 = create_optimization_result()
+        result_2 = create_optimization_result()
+
+        # test a standard call
+        pypesto.visualize.waterfall(result_1)
+
+        # test plotting of lists
+        pypesto.visualize.waterfall([result_1, result_2])
+
+    @staticmethod
+    def test_waterfall_with_options():
+        # create the necessary results
+        result_1 = create_optimization_result()
+        result_2 = create_optimization_result()
+
+        # alternative figure size and plotting options
+        (_, _, ref3, _, ref_point) = create_plotting_options()
+        alt_fig_size = (9.0, 8.0)
+
+        # Test with y-limits as vector
+        pypesto.visualize.waterfall(result_1,
+                                    reference=ref_point,
+                                    y_limits=[-0.5, 2.5],
+                                    start_indices=[0, 1, 4, 11],
+                                    size=alt_fig_size,
+                                    colors=[1., .3, .3, 0.5])
+
+        # Test with linear scale
+        pypesto.visualize.waterfall([result_1, result_2],
+                                    reference=ref3,
+                                    offset_y=-2.5,
+                                    start_indices=3,
+                                    y_limits=5.)
+
+        # Test with y-limits as float
+        pypesto.visualize.waterfall(result_1,
+                                    reference=ref3,
+                                    scale_y='lin',
+                                    offset_y=0.2,
+                                    y_limits=5.)
 
     @staticmethod
     def test_waterfall_lowlevel():
@@ -59,46 +162,51 @@ class TestVisualize(unittest.TestCase):
         fvals = np.array(fvals)
         pypesto.visualize.waterfall_lowlevel(fvals)
 
-    def test_assign_clusters(self):
-        # test empty input
-        pypesto.visualize.assign_clusters([])
+    @staticmethod
+    def test_parameters():
+        # create the necessary results
+        result_1 = create_optimization_result()
+        result_2 = create_optimization_result()
 
-        # test if it runs at all
-        fvals = [0.01, 0.02, 1.01, 2.02, 2.03, 2.04, 3, 4, 4.1, 4.11, 10]
-        pypesto.visualize.assign_clusters(fvals)
-        fvals = np.array(fvals)
-        clust, clustsize, ind_clust = pypesto.visualize.assign_clusters(fvals)
+        # test a standard call
+        pypesto.visualize.parameters(result_1)
 
-        # test if clustering works as intended
-        fvals = [0., 0.00001, 1., 2., 2.001]
-        clust, clustsize, ind_clust = pypesto.visualize.assign_clusters(fvals)
-        self.assertEqual(len(clust), 5)
-        self.assertEqual(len(clustsize), 3)
-        self.assertEqual(len(ind_clust), 5)
+        # test plotting of lists
+        pypesto.visualize.parameters([result_1, result_2])
 
     @staticmethod
-    def test_assign_clustered_colors():
-        # test empty input
-        pypesto.visualize.assign_clustered_colors([])
+    def test_parameters_with_options():
+        # create the necessary results
+        result_1 = create_optimization_result()
+        result_2 = create_optimization_result()
 
-        # test if it runs at all
-        fvals = [0.01, 0.02, 1.01, 2.02, 2.03, 2.04, 3, 4, 4.1, 4.11]
-        pypesto.visualize.assign_clustered_colors(fvals)
-        fvals = np.array(fvals)
-        pypesto.visualize.assign_clustered_colors(fvals)
+        # alternative figure size and plotting options
+        (_, _, _, _, ref_point) = create_plotting_options()
+        alt_fig_size = (9.0, 8.0)
 
-    @staticmethod
-    def test_parameters():
-        pypesto.visualize.parameters(result)
+        # test calls with specific options
+        pypesto.visualize.parameters(result_1,
+                                     free_indices_only=False,
+                                     reference=ref_point,
+                                     size=alt_fig_size,
+                                     colors=[1., .3, .3, 0.5])
+
+        pypesto.visualize.parameters([result_1, result_2],
+                                     free_indices_only=False,
+                                     reference=ref_point)
 
     @staticmethod
     def test_parameters_lowlevel():
         # test empty input
         xs = np.array([])
+
+        # test partial input
         xs.shape = (0, 0)  # we can assume in input that xs.ndim == 2
         fvals = np.array([])
         pypesto.visualize.parameters_lowlevel(xs, fvals)
 
+        # create some dummy results
+        (lb, ub) = create_bounds()
         fvals = [0.01, 0.02, 1.01, 2.02, 2.03, 2.04, 3, 4, 4.1, 4.11]
         xs = [[0.1, 1], [1.2, 3], [2, 4], [1.2, 4.1], [1.1, 3.5],
               [4.2, 3.5], [1, 4], [6.2, 5], [4.3, 3], [3, 2]]
@@ -116,14 +224,37 @@ class TestVisualize(unittest.TestCase):
 
     @staticmethod
     def test_profiles():
-        pypesto.visualize.profiles(result)
+        # create the necessary results
+        result_1 = create_profile_result()
+        result_2 = create_profile_result()
+
+        # test a standard call
+        pypesto.visualize.profiles(result_1)
+
+        # test plotting of lists
+        pypesto.visualize.profiles([result_1, result_2])
+
+    @staticmethod
+    def test_profiles_with_options():
+        # create the necessary results
+        result = create_profile_result()
+
+        # alternative figure size and plotting options
+        (_, _, _, _, ref_point) = create_plotting_options()
+        alt_fig_size = (9.0, 8.0)
+
+        # test a call with some specific options
+        pypesto.visualize.profiles(result,
+                                   reference=ref_point,
+                                   size=alt_fig_size,
+                                   colors=[1., .3, .3, 0.5])
 
     @staticmethod
     def test_profiles_lowlevel():
         # test empty input
         pypesto.visualize.profiles_lowlevel([])
 
-        # test if it runs at all
+        # test if it runs at all using dummy results
         p1 = np.array([[2., 2.1, 2.3, 2.5, 2.7, 2.9, 3.],
                        [0.15, 0.25, 0.7, 1., 0.8, 0.35, 0.15]])
         p2 = np.array([[1., 1.2, 1.4, 1.5, 1.6, 1.8, 2.],
@@ -136,11 +267,164 @@ class TestVisualize(unittest.TestCase):
         # test empty input
         pypesto.visualize.profile_lowlevel(fvals=[])
 
-        # test if it runs at all
+        # test if it runs at all using dummy results
         fvals = np.array([[2., 2.1, 2.3, 2.5, 2.7, 2.9, 3.],
                           [0.15, 0.25, 0.7, 1., 0.8, 0.35, 0.15]])
         pypesto.visualize.profile_lowlevel(fvals=fvals)
 
+    @staticmethod
+    def test_optimizer_history():
+        # create the necessary results
+        result_1 = create_optimization_history()
+        result_2 = create_optimization_history()
+
+        # test a standard call
+        pypesto.visualize.optimizer_history(result_1)
+
+        # test plotting of lists
+        pypesto.visualize.optimizer_history([result_1,
+                                             result_2])
+
+    @staticmethod
+    def test_optimizer_history_with_options():
+        # create the necessary results
+        result_1 = create_optimization_history()
+        result_2 = create_optimization_history()
+
+        # alternative figure size and plotting options
+        (_, _, ref3, _, ref_point) = create_plotting_options()
+        alt_fig_size = (9.0, 8.0)
+
+        # Test with y-limits as vector
+        pypesto.visualize.optimizer_history(result_1,
+                                            y_limits=[-0.5, 2.5],
+                                            start_indices=[0, 1, 4, 11],
+                                            reference=ref_point,
+                                            size=alt_fig_size,
+                                            trace_x='steps',
+                                            trace_y='fval',
+                                            offset_y=-10.,
+                                            colors=[1., .3, .3, 0.5])
+
+        # Test with linear scale
+        pypesto.visualize.optimizer_history([result_1,
+                                             result_2],
+                                            y_limits=[-0.5, 2.5],
+                                            start_indices=[0, 1, 4, 11],
+                                            reference=ref_point,
+                                            size=alt_fig_size,
+                                            scale_y='lin')
+
+        # Test with y-limits as float
+        pypesto.visualize.optimizer_history(result_1,
+                                            y_limits=5.,
+                                            start_indices=3,
+                                            reference=ref3,
+                                            trace_x='time',
+                                            trace_y='gradnorm',
+                                            offset_y=10.)
+
+    @staticmethod
+    def test_optimizer_history_lowlevel():
+        # test empty input
+        pypesto.visualize.optimizer_history_lowlevel([])
+
+        # pass numpy array
+        x_vals = np.array(list(range(10)))
+        y_vals = 11. * np.ones(10) - x_vals
+        vals1 = np.array([x_vals, y_vals])
+        vals2 = np.array([0.1 * np.ones(10) + x_vals, np.ones(10) + y_vals])
+        vals = [vals1, vals2]
+
+        # test with numpy arrays
+        pypesto.visualize.optimizer_history_lowlevel(vals1)
+        pypesto.visualize.optimizer_history_lowlevel(vals2)
+
+        # test with a list of arrays
+        pypesto.visualize.optimizer_history_lowlevel(vals)
+
+    def test_assign_clusters(self):
+        # test empty input
+        pypesto.visualize.assign_clusters([])
+
+        # test if it runs at all
+        fvals = [0.01, 0.02, 1.01, 2.02, 2.03, 2.04, 3, 4, 4.1, 4.11, 10]
+        pypesto.visualize.assign_clusters(fvals)
+        fvals = np.array(fvals)
+        clust, clustsize, ind_clust = pypesto.visualize.assign_clusters(fvals)
+
+        # test if clustering works as intended
+        fvals = [0., 0.00001, 1., 2., 2.001]
+        clust, clustsize, ind_clust = pypesto.visualize.assign_clusters(fvals)
+        self.assertEqual(len(clust), 5)
+        self.assertEqual(len(clustsize), 3)
+        self.assertEqual(len(ind_clust), 5)
+
+    @staticmethod
+    def test_assign_clustered_colors():
+        # test empty input
+        pypesto.visualize.assign_clustered_colors([])
+
+        # test if it runs at all
+        fvals = [0.01, 0.02, 1.01, 2.02, 2.03, 2.04, 3, 4, 4.1, 4.11]
+        pypesto.visualize.assign_clustered_colors(fvals)
+        fvals = np.array(fvals)
+        pypesto.visualize.assign_clustered_colors(fvals)
+
+    @staticmethod
+    def test_assign_colors():
+        # test empty input
+        pypesto.visualize.assign_colors([])
+
+        # test if it runs at all
+        fvals = [0.01, 0.02, 1.01, 2.02, 2.03, 2.04, 3, 4, 4.1, 4.11]
+        pypesto.visualize.assign_colors(fvals)
+        fvals = np.array(fvals)
+        pypesto.visualize.assign_colors(fvals)
+        fvals = [0.01, 0.02, 1.]
+        pypesto.visualize.assign_colors(fvals, colors=[.5, .9, .9, .3])
+        pypesto.visualize.assign_colors(fvals, colors=[[.5, .9, .9, .3],
+                                                       [.5, .8, .8, .5],
+                                                       [.9, .1, .1, .1]])
+
+    @staticmethod
+    def test_reference_points():
+        # create the necessary results
+        (ref1, ref2, ref3, ref4, _) = create_plotting_options()
+
+        # Try conversion from different inputs
+        pypesto.visualize.create_references(ref1)
+        pypesto.visualize.create_references(references=ref2)
+        ref_list_1 = pypesto.visualize.create_references(ref3)
+        ref_list_2 = pypesto.visualize.create_references(ref4)
+
+        # Try to append to a list
+        ref_list_2.append(ref_list_1[0])
+        ref_list_2 = pypesto.visualize.create_references(ref_list_2)
+
+        # Try to append one point to a list via the interface
+        pypesto.visualize.create_references(references=ref_list_2,
+                                            x=ref2[0], fval=ref2[1])
+
+    @staticmethod
+    def test_process_result_list():
+        # create the necessary results
+        result_1 = create_optimization_result()
+        result_2 = create_optimization_result()
+
+        # Test empty arguments
+        pypesto.visualize.process_result_list([])
+
+        # Test single argument
+        # Test single argument
+        pypesto.visualize.process_result_list(result_1)
+
+        # Test handling of a real list
+        res_list = [result_1]
+        pypesto.visualize.process_result_list(res_list)
+        res_list.append(result_2)
+        pypesto.visualize.process_result_list(res_list)
+
 
 if __name__ == '__main__':
     unittest.main()