diff --git a/pylinkage/__init__.py b/pylinkage/__init__.py
index 807145d..558a81a 100644
--- a/pylinkage/__init__.py
+++ b/pylinkage/__init__.py
@@ -20,12 +20,15 @@
     circle_intersect,
     intersection
 )
-from .interface import (
+from .exceptions import (
     UnbuildableError,
     HypostaticError,
     NotCompletelyDefinedError,
-    Pivot,
+)
+from .linkage import (
     Linkage,
+    kinematic_default_test,
+    bounding_box,
 )
 from .joints import (
     Crank,
@@ -34,16 +37,13 @@
     Revolute,
     Static,
 )
+from .joints.revolute import Pivot
 from .optimization import (
     generate_bounds,
     trials_and_errors_optimization,
-    particle_swarm_optimization
-)
-from .utility import (
-    kinematic_default_test,
+    particle_swarm_optimization,
     kinematic_maximization,
     kinematic_minimization,
-    bounding_box
 )
 from .visualizer import (
     plot_static_linkage,
diff --git a/pylinkage/collections/__init__.py b/pylinkage/collections/__init__.py
index 8c4be37..b8c9334 100644
--- a/pylinkage/collections/__init__.py
+++ b/pylinkage/collections/__init__.py
@@ -1,2 +1,5 @@
+"""
+Package for collection objects.
+"""
 from .agent import Agent
 from .mutable_agent import MutableAgent
diff --git a/pylinkage/interface/exceptions.py b/pylinkage/exceptions.py
similarity index 100%
rename from pylinkage/interface/exceptions.py
rename to pylinkage/exceptions.py
diff --git a/pylinkage/geometry/__init__.py b/pylinkage/geometry/__init__.py
new file mode 100644
index 0000000..fb216d0
--- /dev/null
+++ b/pylinkage/geometry/__init__.py
@@ -0,0 +1,16 @@
+"""
+Basic geometry package.
+"""
+
+from .core import (
+    dist,
+    sqr_dist,
+    norm,
+    cyl_to_cart,
+)
+from .secants import (
+    circle_intersect,
+    circle_line_intersection,
+    circle_line_from_points_intersection,
+    intersection,
+)
diff --git a/pylinkage/geometry/core.py b/pylinkage/geometry/core.py
new file mode 100644
index 0000000..29ca4d4
--- /dev/null
+++ b/pylinkage/geometry/core.py
@@ -0,0 +1,101 @@
+"""
+Basic geometry features.
+"""
+import sys
+import warnings
+import math
+
+
+def dist_builtin(point1, point2):
+    """Euclidian distance between two 2D points.
+
+    Legacy built-in unoptimized equivalent of `math.dist` in Python 3.8.
+
+    :param tuple[float, float] point1: First point
+    :param tuple[float, float] point2: Second point
+
+    """
+    return math.sqrt(
+        (point1[0] - point2[0]) ** 2 + (point1[1] - point2[1]) ** 2
+    )
+
+
+if sys.version_info >= (3, 8, 0):
+    dist = math.dist
+else:
+    warnings.warn('Unable to import dist from math. Using built-in function.')
+    dist = dist_builtin
+
+
+def sqr_dist(point1, point2):
+    """
+    Square of the distance between two points.
+
+    Faster than dist.
+
+    :param tuple[float, float] point1: First point to compare
+    :param tuple[float, float] point2: Second point
+
+    :return float: Computed distance
+    """
+    return (point1[0] - point2[0]) ** 2 + (point1[1] - point2[1]) ** 2
+
+
+def get_nearest_point(reference_point, first_point, second_point):
+    """
+    Return the point closer to the reference.
+
+    :param tuple[float, float] reference_point: Point to compare to
+    :param tuple[float, float] first_point: First point candidate
+    :param tuple[float, float] second_point: Second point candidate
+    :return tuple[float, float]: Either first point or second point
+    """
+    if reference_point == first_point or reference_point == second_point:
+        return reference_point
+    if sqr_dist(reference_point, first_point) < sqr_dist(reference_point, second_point):
+        return first_point
+    return second_point
+
+
+def norm(vec):
+    """
+    Return the norm of a 2-dimensional vector.
+
+    :param tuple[float, float] vec: Vector to get norm from
+    """
+    return math.sqrt(vec[0] ** 2 + vec[1] ** 2)
+
+
+def cyl_to_cart(radius, theta, ori=(0, 0)):
+    """Convert polar coordinates into cartesian.
+
+    :param radius: distance from ori
+    :param theta: angle is the angle starting from abscissa axis
+    :param ori: origin point (Default value = (0)).
+
+    """
+    return radius * math.cos(theta) + ori[0], radius * math.sin(theta) + ori[1]
+
+
+def line_from_points(first_point, second_point):
+    """
+    A cartesian equation of the line joining two points.
+
+    :param tuple[float, float] first_point: One point of the line.
+    :param tuple[float, float] second_point: Another point on the line.
+    :return tuple[float, float, float]: A cartesian equation of this line.
+    """
+    if first_point == second_point:
+        warnings.warn("Cannot choose a line, inputs points are the same!")
+        return 0, 0, 0
+    director = (
+        second_point[0] - first_point[0],
+        second_point[1] - first_point[1]
+    )
+    # The barycenter should give more precision
+    mean = (
+        (first_point[0] + second_point[0]) / 2,
+        (first_point[1] + second_point[1]) / 2
+    )
+    equilibrium = mean[0] * director[1] - mean[1] * director[0]
+    return -director[1], director[0], equilibrium
diff --git a/pylinkage/geometry.py b/pylinkage/geometry/secants.py
similarity index 71%
rename from pylinkage/geometry.py
rename to pylinkage/geometry/secants.py
index 79b5a07..fbbdf6b 100644
--- a/pylinkage/geometry.py
+++ b/pylinkage/geometry/secants.py
@@ -10,84 +10,13 @@
 @author: HugoFara
 """
 import math
-import warnings
 
-
-def dist_builtin(point1, point2):
-    """Euclidian distance between two 2D points.
-
-    Legacy built-in unoptimized equivalent of math.dist in Python 3.8.
-
-    :param point1:
-    :param point2:
-
-    """
-    return math.sqrt(
-        (point1[0] - point2[0]) ** 2 + (point1[1] - point2[1]) ** 2
-    )
-
-
-if hasattr(math, 'dist'):
-    dist = math.dist
-else:
-    print('Unable to import dist from math. Using built-in function.')
-    dist = dist_builtin
-
-
-def sqr_dist(point1, point2):
-    """
-    Square of the distance between two points.
-
-    Faster than dist.
-
-    :param tuple[float, float] point1: First point to compare
-    :param tuple[float, float] point2: Second point
-
-    :return float: Computed distance
-    """
-    return (point1[0] - point2[0]) ** 2 + (point1[1] - point2[1]) ** 2
-
-
-def get_nearest_point(reference_point, first_point, second_point):
-    """
-    Return the point closer to the reference.
-
-    :param tuple[float, float] reference_point: Point to compare to
-    :param tuple[float, float] first_point: First point candidate
-    :param tuple[float, float] second_point: Second point candidate
-    :return tuple[float, float]: Either first point or second point
-    """
-    if reference_point == first_point or reference_point == second_point:
-        return reference_point
-    if sqr_dist(reference_point, first_point) < sqr_dist(reference_point, second_point):
-        return first_point
-    return second_point
-
-
-def norm(vec):
-    """
-    Return the norm of a 2-dimensional vector.
-
-    :param vec:
-
-    """
-    return math.sqrt(vec[0] ** 2 + vec[1] ** 2)
-
-
-def cyl_to_cart(radius, theta, ori=(0, 0)):
-    """Convert polar coordinates into cartesian.
-
-    :param radius: distance from ori
-    :param theta: angle is the angle starting from abscissa axis
-    :param ori: origin point (Default value = (0)).
-
-    """
-    return radius * math.cos(theta) + ori[0], radius * math.sin(theta) + ori[1]
+from .core import dist
 
 
 def secant_circles_intersections(
-        distance, dist_x, dist_y, mid_dist, radius1, projected
-    ):
+    distance, dist_x, dist_y, mid_dist, radius1, projected
+):
     """Return the TWO intersections of secant circles."""
     # Distance between projected P and points
     # and the points of which P is projection
@@ -166,30 +95,6 @@ def circle_intersect(circle1, circle2, tol=0.0):
     return 1, projected
 
 
-def line_from_points(first_point, second_point):
-    """
-    A cartesian equation of the line joining two points.
-
-    :param tuple[float, float] first_point: One point of the line.
-    :param tuple[float, float] second_point: Another point on the line.
-    :return tuple[float, float, float]: A cartesian equation of this line.
-    """
-    if first_point == second_point:
-        warnings.warn("Cannot choose a line, inputs points are the same!")
-        return 0, 0, 0
-    director = (
-        second_point[0] - first_point[0],
-        second_point[1] - first_point[1]
-    )
-    # The barycenter should give more precision
-    mean = (
-        (first_point[0] + second_point[0]) / 2,
-        (first_point[1] + second_point[1]) / 2
-    )
-    equilibrium = mean[0] * director[1] - mean[1] * director[0]
-    return -director[1], director[0], equilibrium
-
-
 def circle_line_from_points_intersection(circle, first_point, second_point):
     """
     Intersection(s) of a circle and a line defined by two points.
diff --git a/pylinkage/interface/__init__.py b/pylinkage/interface/__init__.py
deleted file mode 100644
index 273730d..0000000
--- a/pylinkage/interface/__init__.py
+++ /dev/null
@@ -1,18 +0,0 @@
-"""
-Classes that represent the main interface with the library.
-"""
-from .exceptions import (
-    UnbuildableError,
-    HypostaticError,
-    NotCompletelyDefinedError,
-)
-from ..joints import (
-    Static,
-    Revolute,
-    Crank,
-    Linear,
-    Fixed,
-)
-# Will be deleted in next major release
-from ..joints.revolute import Pivot
-from .linkage import Linkage
diff --git a/pylinkage/joints/crank.py b/pylinkage/joints/crank.py
index e21c697..9f1d34c 100644
--- a/pylinkage/joints/crank.py
+++ b/pylinkage/joints/crank.py
@@ -5,7 +5,7 @@
 
 from . import joint as pl_joint
 from .. import geometry as pl_geom
-from ..interface import exceptions as pl_exceptions
+from .. import exceptions as pl_exceptions
 
 
 class Crank(pl_joint.Joint):
diff --git a/pylinkage/joints/fixed.py b/pylinkage/joints/fixed.py
index 5cb76cf..1c3e24c 100644
--- a/pylinkage/joints/fixed.py
+++ b/pylinkage/joints/fixed.py
@@ -4,7 +4,7 @@
 import math
 
 from .. import geometry as pl_geom
-from ..interface import exceptions as pl_exceptions
+from .. import exceptions as pl_exceptions
 from . import joint as pl_joint
 
 
diff --git a/pylinkage/joints/linear.py b/pylinkage/joints/linear.py
index e9be0b7..02db005 100644
--- a/pylinkage/joints/linear.py
+++ b/pylinkage/joints/linear.py
@@ -2,7 +2,7 @@
 Definition of a linear joint.
 """
 from . import joint as pl_joint
-from ..interface import exceptions as pl_exceptions
+from .. import exceptions as pl_exceptions
 from .. import geometry as geom
 
 
diff --git a/pylinkage/joints/revolute.py b/pylinkage/joints/revolute.py
index 5cb01ee..25f3842 100644
--- a/pylinkage/joints/revolute.py
+++ b/pylinkage/joints/revolute.py
@@ -8,7 +8,7 @@
 from math import atan2
 
 from .. import geometry as pl_geom
-from ..interface import exceptions as pl_exceptions
+from .. import exceptions as pl_exceptions
 from . import joint as pl_joint
 
 
@@ -113,7 +113,9 @@ def reload(self):
             if intersections[0] == 1:
                 self.x, self.y = intersections[1]
             elif intersections[0] == 2:
-                self.x, self.y = pl_geom.get_nearest_point(self.coord(), intersections[1], intersections[2])
+                self.x, self.y = pl_geom.core.get_nearest_point(
+                    self.coord(), intersections[1], intersections[2]
+                )
             elif intersections[0] == 3:
                 warnings.warn(
                     f"Joint {self.name} has an infinite number of"
diff --git a/pylinkage/linkage/__init__.py b/pylinkage/linkage/__init__.py
new file mode 100644
index 0000000..a382f33
--- /dev/null
+++ b/pylinkage/linkage/__init__.py
@@ -0,0 +1,8 @@
+"""
+Definition and analysis of a linkage as a dynamic set of joints.
+"""
+from .linkage import Linkage
+from .analysis import (
+    kinematic_default_test,
+    bounding_box,
+)
\ No newline at end of file
diff --git a/pylinkage/utility.py b/pylinkage/linkage/analysis.py
similarity index 68%
rename from pylinkage/utility.py
rename to pylinkage/linkage/analysis.py
index 3250107..7314b7b 100644
--- a/pylinkage/utility.py
+++ b/pylinkage/linkage/analysis.py
@@ -1,18 +1,12 @@
 """
-The utility module provides various useful functions.
-
-Created on Mon Jul 12 00:00:01 2021.
-
-@author: HugoFara
+Analysis tools for linkages.
 """
-import warnings
-
-from pylinkage.interface.exceptions import UnbuildableError
+from ..exceptions import UnbuildableError
 
 
 def kinematic_default_test(func, error_penalty):
     """Standard run for any linkage before a complete fitness evaluation.
-    
+
     This decorator makes a kinematic simulation, before passing the loci to the
     decorated function.
 
@@ -22,8 +16,9 @@ def kinematic_default_test(func, error_penalty):
             float('inf') and 0.
     :type error_penalty: float
 
-    
+
     """
+
     def wrapper(linkage, params, init_pos=None):
         """Decorated function.
 
@@ -35,7 +30,7 @@ def wrapper(linkage, params, init_pos=None):
             redefined at each successful iteration. (Default value = None)
         :type init_pos: tuple[tuple[float]]
 
-        
+        :return Callable: New optimization function wrapper
         """
         if init_pos is not None:
             linkage.set_coords(init_pos)
@@ -66,35 +61,8 @@ def wrapper(linkage, params, init_pos=None):
             return func(
                 linkage=linkage, params=params, init_pos=init_pos, loci=loci
             )
-    return wrapper
-
-
-def kinematic_maximization(func):
-    """Standard run for any linkage before a complete fitness evaluation.
-    
-    This decorator makes a kinematic simulation, before passing the loci to the
-    decorated function. In case of error, the penalty value is -float('inf')
-
-    :param func: Fitness function to be decorated.
-    :type func: Callable
 
-    
-    """
-    return kinematic_default_test(func, -float('inf'))
-
-
-def kinematic_minimization(func):
-    """Standard run for any linkage before a complete fitness evaluation.
-    
-    This decorator makes a kinematic simulation, before passing the loci to the
-    decorated function. In case of error, the penalty value is float('inf')
-
-    :param func: Fitness function to be decorated.
-    :type func: Callable
-
-    
-    """
-    return kinematic_default_test(func, float('inf'))
+    return wrapper
 
 
 def bounding_box(locus):
@@ -135,18 +103,3 @@ def movement_bounding_box(loci):
             max(new_bb[2], bb[2]), min(new_bb[3], bb[3])
         )
     return bb
-
-
-def movement_bounding_bow(loci):
-    """
-    Bounding box for a group of loci.
-
-    :param loci:
-
-    :rtype: tuple[float, float, float, float]
-
-    .. deprecated :: 0.6.0
-        Was replaced by :func:`movement_bounding_box`, will be removed in 0.7.0.
-    """
-    warnings.warn("movement_bounding_bow is deprecated, please use movement_bounding_box")
-    return movement_bounding_box(loci)
diff --git a/pylinkage/interface/linkage.py b/pylinkage/linkage/linkage.py
similarity index 99%
rename from pylinkage/interface/linkage.py
rename to pylinkage/linkage/linkage.py
index f2badbb..99f8046 100644
--- a/pylinkage/interface/linkage.py
+++ b/pylinkage/linkage/linkage.py
@@ -9,7 +9,7 @@
 """
 import warnings
 from math import gcd, tau
-from .exceptions import HypostaticError
+from ..exceptions import HypostaticError
 from ..joints import (Revolute, Fixed, Crank, Static)
 
 
diff --git a/pylinkage/optimization/__init__.py b/pylinkage/optimization/__init__.py
index aa084ae..fb92e6d 100644
--- a/pylinkage/optimization/__init__.py
+++ b/pylinkage/optimization/__init__.py
@@ -1,4 +1,8 @@
 """Optimization package."""
-from .utils import generate_bounds
 from .grid_search import trials_and_errors_optimization
 from .particle_swarm import particle_swarm_optimization
+from .utils import (
+    generate_bounds,
+    kinematic_maximization,
+    kinematic_minimization,
+)
\ No newline at end of file
diff --git a/pylinkage/optimization/utils.py b/pylinkage/optimization/utils.py
index cfe037d..e305e10 100644
--- a/pylinkage/optimization/utils.py
+++ b/pylinkage/optimization/utils.py
@@ -1,6 +1,14 @@
-"""Utility for optimization."""
+"""
+This utility module provides various useful functions for optimization.
+
+Created on Mon Jul 12 00:00:01 2021.
+
+@author: HugoFara
+"""
 import numpy as np
 
+from ..linkage.analysis import kinematic_default_test
+
 
 def generate_bounds(center, min_ratio=5, max_factor=5):
     """Simple function to generate bounds from a linkage.
@@ -20,3 +28,31 @@ def generate_bounds(center, min_ratio=5, max_factor=5):
     """
     np_center = np.array(center)
     return np_center / min_ratio, np_center * max_factor
+
+
+def kinematic_maximization(func):
+    """Standard run for any linkage before a complete fitness evaluation.
+
+    This decorator makes a kinematic simulation, before passing the loci to the
+    decorated function. In case of error, the penalty value is -float('inf')
+
+    :param func: Fitness function to be decorated.
+    :type func: Callable
+
+
+    """
+    return kinematic_default_test(func, -float('inf'))
+
+
+def kinematic_minimization(func):
+    """Standard run for any linkage before a complete fitness evaluation.
+
+    This decorator makes a kinematic simulation, before passing the loci to the
+    decorated function. In case of error, the penalty value is float('inf')
+
+    :param func: Fitness function to be decorated.
+    :type func: Callable
+
+
+    """
+    return kinematic_default_test(func, float('inf'))
diff --git a/pylinkage/visualizer/__init__.py b/pylinkage/visualizer/__init__.py
new file mode 100644
index 0000000..f47f1d8
--- /dev/null
+++ b/pylinkage/visualizer/__init__.py
@@ -0,0 +1,10 @@
+"""
+Linkage visualization features.
+"""
+from .core import COLOR_SWITCHER
+from .static import plot_static_linkage
+from .animated import (
+    plot_kinematic_linkage,
+    show_linkage,
+    swarm_tiled_repr,
+)
diff --git a/pylinkage/visualizer.py b/pylinkage/visualizer/animated.py
similarity index 70%
rename from pylinkage/visualizer.py
rename to pylinkage/visualizer/animated.py
index 8c823c8..8d3168f 100644
--- a/pylinkage/visualizer.py
+++ b/pylinkage/visualizer/animated.py
@@ -9,112 +9,16 @@
 """
 import matplotlib.pyplot as plt
 import matplotlib.animation as anim
-from .utility import movement_bounding_box
-from .interface.exceptions import UnbuildableError
-from .interface import Crank, Fixed, Static, Pivot
-from .joints import Revolute, Linear
+
+from ..linkage.analysis import movement_bounding_box
+from ..exceptions import UnbuildableError
+from ..joints import Crank, Static
+from .static import plot_static_linkage
+from .core import _get_color
 
 # List of animations
 ANIMATIONS = []
 
-# Colors to use for plotting
-COLOR_SWITCHER = {
-    Static: 'k',
-    Crank: 'g',
-    Fixed: 'r',
-    Pivot: 'b',
-    Revolute: 'b',
-    Linear: 'orange'
-}
-
-
-def _get_color(joint):
-    """Search in COLOR_SWITCHER for the corresponding color.
-
-    :param joint:
-
-    """
-    for joint_type, color in COLOR_SWITCHER.items():
-        if isinstance(joint, joint_type):
-            return color
-    return ''
-
-
-def plot_static_linkage(
-        linkage, axis, loci, locus_highlights=None,
-        show_legend=False
-):
-    """Plot a linkage without movement.
-
-    :param linkage: The linkage you want to see.
-    :type linkage: Linkage
-    :param axis: The graph we should draw on.
-    :type axis: Artist
-    :param loci: List of list of coordinates. They will be plotted.
-    :type loci: Iterable
-    :param locus_highlights: If a list, should be a list of list of coordinates you want to see
-        highlighted. The default is None.
-    :type locus_highlights: list
-    :param show_legend: To add an automatic legend to the graph. The default is False.
-    :type show_legend: bool
-
-
-    """
-    axis.set_aspect('equal')
-    axis.grid(True)
-    # Plot loci
-    for i, joint in enumerate(linkage.joints):
-        axis.plot(tuple(j[i][0] for j in loci), tuple(j[i][1] for j in loci))
-
-    # The plot linkage in initial positioning
-    # It as important to use separate loops, because we would have bad
-    # formatted legend otherwise
-    for i, joint in enumerate(linkage.joints):
-        # Then the linkage in initial position
-
-        # Draw a link to the first parent if it exists
-        if joint.joint0 is None:
-            continue
-        pos = joint.coord()
-        par_pos = joint.joint0.coord()
-        plot_kwargs = {
-            "c": _get_color(joint),
-            "linewidth": .3
-        }
-        axis.plot(
-            [par_pos[0], pos[0]], [par_pos[1], pos[1]],
-            **plot_kwargs
-        )
-        # Then second parent
-        if isinstance(joint, (Fixed, Pivot, Revolute)):
-            par_pos = joint.joint1.coord()
-            axis.plot(
-                [par_pos[0], pos[0]], [par_pos[1], pos[1]],
-                **plot_kwargs
-            )
-        elif isinstance(joint, Linear):
-            # Different ordering
-            par_pos = joint.joint2.coord()
-            other_pos = joint.joint1.coord()
-            axis.plot(
-                [par_pos[0], other_pos[0]], [par_pos[1], other_pos[1]],
-                **plot_kwargs
-            )
-
-    # Highlight for specific loci
-    if locus_highlights:
-        for locus in locus_highlights:
-            axis.scatter(
-                tuple(coord[0] for coord in locus),
-                tuple(coord[1] for coord in locus)
-            )
-
-    if show_legend:
-        axis.set_title("Static representation")
-        axis.set_xlabel("x")
-        axis.set_ylabel("y")
-        axis.legend(tuple(i.name for i in linkage.joints))
-
 
 def update_animated_plot(linkage, index, images, loci):
     """Modify im, instead of recreating it to make the animation run faster.
diff --git a/pylinkage/visualizer/core.py b/pylinkage/visualizer/core.py
new file mode 100644
index 0000000..a87f642
--- /dev/null
+++ b/pylinkage/visualizer/core.py
@@ -0,0 +1,34 @@
+"""
+Core features for visualization.
+"""
+from ..joints.revolute import Pivot
+from ..joints import (
+    Revolute,
+    Linear,
+    Crank,
+    Fixed,
+    Static,
+)
+
+# Colors to use for plotting
+COLOR_SWITCHER = {
+    Static: 'k',
+    Crank: 'g',
+    Fixed: 'r',
+    Pivot: 'b',
+    Revolute: 'b',
+    Linear: 'orange'
+}
+
+
+def _get_color(joint):
+    """Search in COLOR_SWITCHER for the corresponding color.
+
+    :param joint:
+
+    """
+    for joint_type, color in COLOR_SWITCHER.items():
+        if isinstance(joint, joint_type):
+            return color
+    return ''
+
diff --git a/pylinkage/visualizer/static.py b/pylinkage/visualizer/static.py
new file mode 100644
index 0000000..1f43d81
--- /dev/null
+++ b/pylinkage/visualizer/static.py
@@ -0,0 +1,83 @@
+"""
+Static (not animated) visualization.
+"""
+
+from .core import _get_color
+from ..joints import (Fixed, Revolute, Linear)
+from ..joints.revolute import Pivot
+
+
+def plot_static_linkage(
+        linkage, axis, loci, locus_highlights=None,
+        show_legend=False
+):
+    """Plot a linkage without movement.
+
+    :param linkage: The linkage you want to see.
+    :type linkage: Linkage
+    :param axis: The graph we should draw on.
+    :type axis: Artist
+    :param loci: List of list of coordinates. They will be plotted.
+    :type loci: Iterable
+    :param locus_highlights: If a list, should be a list of list of coordinates you want to see
+        highlighted. The default is None.
+    :type locus_highlights: list
+    :param show_legend: To add an automatic legend to the graph. The default is False.
+    :type show_legend: bool
+
+
+    """
+    axis.set_aspect('equal')
+    axis.grid(True)
+    # Plot loci
+    for i, joint in enumerate(linkage.joints):
+        axis.plot(tuple(j[i][0] for j in loci), tuple(j[i][1] for j in loci))
+
+    # The plot linkage in initial positioning
+    # It as important to use separate loops, because we would have bad
+    # formatted legend otherwise
+    for i, joint in enumerate(linkage.joints):
+        # Then the linkage in initial position
+
+        # Draw a link to the first parent if it exists
+        if joint.joint0 is None:
+            continue
+        pos = joint.coord()
+        par_pos = joint.joint0.coord()
+        plot_kwargs = {
+            "c": _get_color(joint),
+            "linewidth": .3
+        }
+        axis.plot(
+            [par_pos[0], pos[0]], [par_pos[1], pos[1]],
+            **plot_kwargs
+        )
+        # Then second parent
+        if isinstance(joint, (Fixed, Pivot, Revolute)):
+            par_pos = joint.joint1.coord()
+            axis.plot(
+                [par_pos[0], pos[0]], [par_pos[1], pos[1]],
+                **plot_kwargs
+            )
+        elif isinstance(joint, Linear):
+            # Different ordering
+            par_pos = joint.joint2.coord()
+            other_pos = joint.joint1.coord()
+            axis.plot(
+                [par_pos[0], other_pos[0]], [par_pos[1], other_pos[1]],
+                **plot_kwargs
+            )
+
+    # Highlight for specific loci
+    if locus_highlights:
+        for locus in locus_highlights:
+            axis.scatter(
+                tuple(coord[0] for coord in locus),
+                tuple(coord[1] for coord in locus)
+            )
+
+    if show_legend:
+        axis.set_title("Static representation")
+        axis.set_xlabel("x")
+        axis.set_ylabel("y")
+        axis.legend(tuple(i.name for i in linkage.joints))
diff --git a/tests/test_optimizer.py b/tests/test_optimizer.py
index d2a8f1c..87f62f2 100644
--- a/tests/test_optimizer.py
+++ b/tests/test_optimizer.py
@@ -4,7 +4,7 @@
 import pylinkage as pl
 from pylinkage import optimization
 from pylinkage.optimization.grid_search import fast_variator, sequential_variator
-from pylinkage.utility import kinematic_minimization
+from pylinkage.optimization.utils import kinematic_minimization
 
 
 def prepare_linkage():