diff --git a/ANN_Training.py b/ANN_Training.py
index 2c5316bb6742c112a6990e0b022a7208a26c8c46..41c71d66800872d5c98b2aa701461401bc9beade 100644
--- a/ANN_Training.py
+++ b/ANN_Training.py
@@ -10,8 +10,6 @@ TODO: Add README for ANN training
 """
 import numpy as np
 import time
-import matplotlib
-from matplotlib import pyplot as plt
 import os
 import torch
 import json
@@ -21,9 +19,6 @@ from sklearn.metrics import accuracy_score, precision_recall_fscore_support, \
     roc_auc_score
 
 import ANN_Model
-from Plotting import plot_classification_barplot, plot_classification_boxplot
-
-matplotlib.use('Agg')
 
 
 class ModelTrainer:
@@ -347,64 +342,3 @@ def evaluate_models(models: dict, directory: str, num_iterations: int = 100,
         json_file.write(json.dumps(classification_stats))
     toc = time.perf_counter()
     print(f'Total runtime: {toc - tic:0.4f}s')
-
-
-def plot_evaluation_results(evaluation_file: str, directory: str,
-                            colors: dict = None) -> None:
-    """Plots given evaluation results of model classifications.
-
-    Plots evaluation results for all measures for which a color is given. If
-    colors is set to None, all measures are plotted with a default color
-    scheme.
-
-    Parameters
-    ----------
-    evaluation_file: str
-        Path to file containing evaluation results.
-    directory : str
-        Path to directory for saving resulting plots.
-    colors : dict, optional
-        Dictionary containing plotting colors. If None, set to default colors.
-        Default: None.
-
-    """
-    tic = time.perf_counter()
-
-    # Set colors if not given
-    if colors is None:
-        colors = {'Accuracy': 'magenta', 'Precision_Smooth': 'red',
-                  'Precision_Troubled': '#8B0000', 'Recall_Smooth': 'blue',
-                  'Recall_Troubled': '#00008B', 'F-Score_Smooth': 'green',
-                  'F-Score_Troubled': '#006400', 'AUROC': 'yellow'}
-
-    # Read evaluation results
-    print('Reading evaluation results.')
-    with open(evaluation_file) as json_file:
-        classification_stats = json.load(json_file)
-
-    # Plot data
-    print('\nPlotting evaluation of trained models...')
-    print('Plotting data in boxplot.')
-    models = classification_stats[list(colors.keys())[0]].keys()
-    plot_classification_boxplot(classification_stats, colors)
-    print('Plotting averaged data in barplot.')
-    classification_stats = {measure: {model: np.array(
-        classification_stats[measure][model]).mean()
-                                      for model in models}
-                            for measure in colors}
-    plot_classification_barplot(classification_stats, colors)
-    print('Finished plotting evaluation of trained models!\n')
-
-    # Set paths for plot files if not existing already
-    plot_dir = directory + '/model evaluation'
-    if not os.path.exists(plot_dir):
-        os.makedirs(plot_dir)
-
-    # Save plots
-    print('Saving plots.')
-    file_name = evaluation_file.split('/')[-1].rstrip('.json')
-    for identifier in plt.get_figlabels():
-        plt.figure(identifier)
-        plt.savefig(plot_dir + '/' + file_name + '.' + identifier + '.pdf')
-    toc = time.perf_counter()
-    print(f'Total runtime: {toc - tic:0.4f}s')
diff --git a/DG_Approximation.py b/DG_Approximation.py
index 58e86da77ef389a07c8c984eda4ec1078ad3c198..d15fb5902ab94634f7f9506d8428345bdb48bbfd 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -5,6 +5,8 @@
 Discussion:
 
 Urgent:
+TODO: Rename Vector to Basis
+TODO: Extract ndarray encoding and decoding
 TODO: Hard-code simplification of cell average/reconstruction in basis
 TODO: Make basis variables public (if feasible)
 TODO: Contain polynomial degree in basis
@@ -13,7 +15,7 @@ TODO: Introduce Mesh class
 TODO: Check whether ghost cells are handled/set correctly
 TODO: Find error in centering for ANN training
 TODO: Investigate g-mesh(?)
-TODO: Contain all plotting in Plotting
+TODO: Contain all plotting in Plotting -> Done
 TODO: Remove use of DGScheme from ANN_Data_Generator
 TODO: Clean up docstrings
 TODO: Adapt TCD from Soraya
@@ -52,14 +54,11 @@ TODO: Discuss adding kwargs to attributes in documentation
 TODO: Add type annotations to function heads
 
 """
-import os
 import json
 import numpy as np
 from sympy import Symbol
 import math
 import seaborn as sns
-import matplotlib
-from matplotlib import pyplot as plt
 
 import Troubled_Cell_Detector
 import Initial_Condition
@@ -67,12 +66,8 @@ import Limiter
 import Quadrature
 import Update_Scheme
 from Basis_Function import OrthonormalLegendre
-from projection_utils import calculate_cell_average, \
-    calculate_exact_solution, calculate_approximate_solution
-from Plotting import plot_solution_and_approx, plot_semilog_error, \
-    plot_error, plot_shock_tube, plot_details
+from projection_utils import calculate_cell_average
 
-matplotlib.use('Agg')
 x = Symbol('x')
 sns.set()
 
@@ -406,195 +401,3 @@ def do_initial_projection(initial_condition, basis, quadrature,
 
     # print(np.array(output_matrix).shape)
     return np.transpose(np.array(output_matrix))
-
-
-def plot_approximation_results(data_file, directory, plot_name, quadrature,
-                               init_cond, basis):
-    """Plots given approximation results.
-
-    Generates plots based on given data, sets plot directory if not
-    already existing, and saves plots.
-
-    Parameters
-    ----------
-    data_file: str
-        Path to data file for plotting.
-    directory: str
-        Path to directory in which plots will be saved.
-    plot_name : str
-        Name of plot.
-    basis: Vector object
-        Basis used for calculation.
-    quadrature: Quadrature object
-        Quadrature used for evaluation.
-    init_cond : InitialCondition object
-        Initial condition used for calculation.
-
-    """
-    # Read approximation results
-    with open(data_file + '.json') as json_file:
-        approx_stats = json.load(json_file)
-
-    # Decode all ndarrays by converting lists
-    approx_stats = {key: decode_ndarray(approx_stats[key])
-                    for key in approx_stats.keys()}
-
-    # Plot exact/approximate results, errors, shock tubes,
-    # and any detector-dependant plots
-    plot_results(quadrature=quadrature, basis=basis,
-                 init_cond=init_cond, **approx_stats)
-
-    # Set paths for plot files if not existing already
-    if not os.path.exists(directory):
-        os.makedirs(directory)
-
-    # Save plots
-    for identifier in plt.get_figlabels():
-        # Set path for figure directory if not existing already
-        if not os.path.exists(directory + '/' + identifier):
-            os.makedirs(directory + '/' + identifier)
-
-        plt.figure(identifier)
-        plt.savefig(directory + '/' + identifier + '/' +
-                    plot_name + '.pdf')
-
-
-def plot_results(projection, troubled_cell_history, time_history, mesh,
-                 num_grid_cells, polynomial_degree, wave_speed, final_time,
-                 left_bound, right_bound, basis, quadrature, init_cond,
-                 colors=None, coarse_projection=None,
-                 multiwavelet_coeffs=None):
-    """Plots results and troubled cells of a projection.
-
-    Plots exact and approximate solution, errors, and troubled cells of a
-    projection given its evaluation history.
-
-    If coarse grid and projection are given, solutions are displayed for
-    both coarse and fine grid. Additionally, coefficient details are plotted.
-
-    Parameters
-    ----------
-    projection : ndarray
-        Matrix of projection for each polynomial degree.
-    troubled_cell_history : list
-        List of detected troubled cells for each time step.
-    time_history : list
-        List of value of each time step.
-    mesh : ndarray
-        List of mesh valuation points.
-    num_grid_cells : int
-        Number of cells in the mesh. Usually exponential of 2.
-    polynomial_degree : int
-        Polynomial degree.
-    wave_speed : float
-        Speed of wave in rightward direction.
-    final_time : float
-        Final time for which approximation is calculated.
-    left_bound : float
-        Left boundary of interval.
-    right_bound : float
-        Right boundary of interval.
-    basis: Vector object
-        Basis used for calculation.
-    quadrature: Quadrature object
-        Quadrature used for evaluation.
-    init_cond : InitialCondition object
-        Initial condition used for calculation.
-    colors: dict
-        Dictionary of colors used for plots.
-    coarse_projection: ndarray, optional
-        Matrix of projection on coarse grid for each polynomial degree.
-        Default: None.
-    multiwavelet_coeffs: ndarray, optional
-        Matrix of wavelet coefficients. Default: None.
-
-    """
-    # Set colors
-    if colors is None:
-        colors = {}
-    colors = _check_colors(colors)
-
-    # Calculate needed variables
-    interval_len = right_bound-left_bound
-    cell_len = interval_len/num_grid_cells
-
-    # Plot troubled cells
-    plot_shock_tube(num_grid_cells, troubled_cell_history, time_history)
-
-    # Determine exact and approximate solution
-    grid, exact = calculate_exact_solution(
-        mesh[2:-2], cell_len, wave_speed,
-        final_time, interval_len, quadrature,
-        init_cond)
-    approx = calculate_approximate_solution(
-        projection[:, 1:-1], quadrature.get_eval_points(),
-        polynomial_degree, basis.get_basis_vector())
-
-    # Plot multiwavelet solution (fine and coarse grid)
-    if coarse_projection is not None:
-        coarse_cell_len = 2*cell_len
-        coarse_mesh = np.arange(left_bound - (0.5*coarse_cell_len),
-                                right_bound + (1.5*coarse_cell_len),
-                                coarse_cell_len)
-
-        # Plot exact and approximate solutions for coarse mesh
-        coarse_grid, coarse_exact = calculate_exact_solution(
-            coarse_mesh[1:-1], coarse_cell_len, wave_speed,
-            final_time, interval_len, quadrature,
-            init_cond)
-        coarse_approx = calculate_approximate_solution(
-            coarse_projection, quadrature.get_eval_points(),
-            polynomial_degree, basis.get_basis_vector())
-        plot_solution_and_approx(
-            coarse_grid, coarse_exact, coarse_approx, colors['coarse_exact'],
-            colors['coarse_approx'])
-
-        # Plot multiwavelet details
-        num_coarse_grid_cells = num_grid_cells//2
-        plot_details(projection[:, 1:-1], mesh[2:-2], coarse_projection,
-                     basis.get_basis_vector(),
-                     basis.get_wavelet_vector(), multiwavelet_coeffs,
-                     num_coarse_grid_cells,
-                     polynomial_degree)
-
-        plot_solution_and_approx(grid, exact, approx,
-                                 colors['fine_exact'],
-                                 colors['fine_approx'])
-        plt.legend(['Exact (Coarse)', 'Approx (Coarse)', 'Exact (Fine)',
-                    'Approx (Fine)'])
-    # Plot regular solution (fine grid)
-    else:
-        plot_solution_and_approx(grid, exact, approx, colors['exact'],
-                                 colors['approx'])
-        plt.legend(['Exact', 'Approx'])
-
-    # Calculate errors
-    pointwise_error = np.abs(exact-approx)
-    max_error = np.max(pointwise_error)
-
-    # Plot errors
-    plot_semilog_error(grid, pointwise_error)
-    plot_error(grid, exact, approx)
-
-    print('p =', polynomial_degree)
-    print('N =', num_grid_cells)
-    print('maximum error =', max_error)
-
-
-def _check_colors(colors):
-    """Checks plot colors.
-
-    Checks whether colors for plots were given and sets them if required.
-
-    """
-    # Set colors for general plots
-    colors['exact'] = colors.get('exact', 'k-')
-    colors['approx'] = colors.get('approx', 'y')
-
-    # Set colors for multiwavelet plots
-    colors['fine_exact'] = colors.get('fine_exact', 'k-.')
-    colors['fine_approx'] = colors.get('fine_approx', 'b-.')
-    colors['coarse_exact'] = colors.get('coarse_exact', 'k-')
-    colors['coarse_approx'] = colors.get('coarse_approx', 'y')
-
-    return colors
diff --git a/Plotting.py b/Plotting.py
index 5c63e04dcce824b480b00d81ced9b58a2ca65940..43b854b6d6ffffac3805cdbe729291ef48c867e3 100644
--- a/Plotting.py
+++ b/Plotting.py
@@ -6,6 +6,9 @@ TODO: Give option to select plotting color
 
 """
 
+import os
+import time
+import json
 import numpy as np
 import matplotlib
 from matplotlib import pyplot as plt
@@ -13,6 +16,13 @@ import seaborn as sns
 from numpy import ndarray
 from sympy import Symbol
 
+from Quadrature import Quadrature
+from Initial_Condition import InitialCondition
+from Basis_Function import Vector
+from projection_utils import calculate_exact_solution,\
+    calculate_approximate_solution
+from DG_Approximation import decode_ndarray
+
 
 matplotlib.use('Agg')
 x = Symbol('x')
@@ -248,3 +258,270 @@ def plot_classification_boxplot(evaluation_dict: dict, colors: dict) -> None:
     ax.legend([bp["boxes"][0] for bp in boxplots], evaluation_dict.keys(),
               loc='center right', bbox_to_anchor=(1.4, 0.75), shadow=True,
               ncol=1, fancybox=True, fontsize=8)
+
+
+def plot_evaluation_results(evaluation_file: str, directory: str,
+                            colors: dict = None) -> None:
+    """Plots given evaluation results of model classifications.
+
+    Plots evaluation results for all measures for which a color is given. If
+    colors is set to None, all measures are plotted with a default color
+    scheme.
+
+    Parameters
+    ----------
+    evaluation_file: str
+        Path to file containing evaluation results.
+    directory : str
+        Path to directory for saving resulting plots.
+    colors : dict, optional
+        Dictionary containing plotting colors. If None, set to default colors.
+        Default: None.
+
+    """
+    tic = time.perf_counter()
+
+    # Set colors if not given
+    if colors is None:
+        colors = {'Accuracy': 'magenta', 'Precision_Smooth': 'red',
+                  'Precision_Troubled': '#8B0000', 'Recall_Smooth': 'blue',
+                  'Recall_Troubled': '#00008B', 'F-Score_Smooth': 'green',
+                  'F-Score_Troubled': '#006400', 'AUROC': 'yellow'}
+
+    # Read evaluation results
+    print('Reading evaluation results.')
+    with open(evaluation_file) as json_file:
+        classification_stats = json.load(json_file)
+
+    # Plot data
+    print('\nPlotting evaluation of trained models...')
+    print('Plotting data in boxplot.')
+    models = classification_stats[list(colors.keys())[0]].keys()
+    plot_classification_boxplot(classification_stats, colors)
+    print('Plotting averaged data in barplot.')
+    classification_stats = {measure: {model: np.array(
+        classification_stats[measure][model]).mean()
+                                      for model in models}
+                            for measure in colors}
+    plot_classification_barplot(classification_stats, colors)
+    print('Finished plotting evaluation of trained models!\n')
+
+    # Set paths for plot files if not existing already
+    plot_dir = directory + '/model evaluation'
+    if not os.path.exists(plot_dir):
+        os.makedirs(plot_dir)
+
+    # Save plots
+    print('Saving plots.')
+    file_name = evaluation_file.split('/')[-1].rstrip('.json')
+    for identifier in plt.get_figlabels():
+        plt.figure(identifier)
+        plt.savefig(plot_dir + '/' + file_name + '.' + identifier + '.pdf')
+    toc = time.perf_counter()
+    print(f'Total runtime: {toc - tic:0.4f}s')
+
+
+def plot_approximation_results(data_file: str, directory: str, plot_name: str,
+                               basis: Vector, quadrature: Quadrature,
+                               init_cond: InitialCondition) -> None:
+    """Plots given approximation results.
+
+    Generates plots based on given data, sets plot directory if not
+    already existing, and saves plots.
+
+    Parameters
+    ----------
+    data_file: str
+        Path to data file for plotting.
+    directory: str
+        Path to directory in which plots will be saved.
+    plot_name : str
+        Name of plot.
+    basis: Vector object
+        Basis used for calculation.
+    quadrature: Quadrature object
+        Quadrature used for evaluation.
+    init_cond : InitialCondition object
+        Initial condition used for calculation.
+
+    """
+    # Read approximation results
+    with open(data_file + '.json') as json_file:
+        approx_stats = json.load(json_file)
+
+    # Decode all ndarrays by converting lists
+    approx_stats = {key: decode_ndarray(approx_stats[key])
+                    for key in approx_stats.keys()}
+
+    # Plot exact/approximate results, errors, shock tubes,
+    # and any detector-dependant plots
+    plot_results(quadrature=quadrature, basis=basis,
+                 init_cond=init_cond, **approx_stats)
+
+    # Set paths for plot files if not existing already
+    if not os.path.exists(directory):
+        os.makedirs(directory)
+
+    # Save plots
+    for identifier in plt.get_figlabels():
+        # Set path for figure directory if not existing already
+        if not os.path.exists(directory + '/' + identifier):
+            os.makedirs(directory + '/' + identifier)
+
+        plt.figure(identifier)
+        plt.savefig(directory + '/' + identifier + '/' +
+                    plot_name + '.pdf')
+
+
+def plot_results(projection: ndarray, troubled_cell_history: list,
+                 time_history: list, mesh: ndarray, num_grid_cells: int,
+                 polynomial_degree: int, wave_speed: float, final_time: float,
+                 left_bound: float, right_bound: float, basis: Vector,
+                 quadrature: Quadrature, init_cond: InitialCondition,
+                 colors: dict = None, coarse_projection: ndarray = None,
+                 multiwavelet_coeffs: ndarray = None) -> None:
+    """Plots results and troubled cells of a projection.
+
+    Plots exact and approximate solution, errors, and troubled cells of a
+    projection given its evaluation history.
+
+    If coarse grid and projection are given, solutions are displayed for
+    both coarse and fine grid. Additionally, coefficient details are plotted.
+
+    Parameters
+    ----------
+    projection : ndarray
+        Matrix of projection for each polynomial degree.
+    troubled_cell_history : list
+        List of detected troubled cells for each time step.
+    time_history : list
+        List of value of each time step.
+    mesh : ndarray
+        List of mesh valuation points.
+    num_grid_cells : int
+        Number of cells in the mesh. Usually exponential of 2.
+    polynomial_degree : int
+        Polynomial degree.
+    wave_speed : float
+        Speed of wave in rightward direction.
+    final_time : float
+        Final time for which approximation is calculated.
+    left_bound : float
+        Left boundary of interval.
+    right_bound : float
+        Right boundary of interval.
+    basis: Vector object
+        Basis used for calculation.
+    quadrature: Quadrature object
+        Quadrature used for evaluation.
+    init_cond : InitialCondition object
+        Initial condition used for calculation.
+    colors: dict
+        Dictionary of colors used for plots.
+    coarse_projection: ndarray, optional
+        Matrix of projection on coarse grid for each polynomial degree.
+        Default: None.
+    multiwavelet_coeffs: ndarray, optional
+        Matrix of wavelet coefficients. Default: None.
+
+    """
+    # Set colors
+    if colors is None:
+        colors = {}
+    colors = _check_colors(colors)
+
+    # Calculate needed variables
+    interval_len = right_bound-left_bound
+    cell_len = interval_len/num_grid_cells
+
+    # Plot troubled cells
+    plot_shock_tube(num_grid_cells, troubled_cell_history, time_history)
+
+    # Determine exact and approximate solution
+    grid, exact = calculate_exact_solution(
+        mesh[2:-2], cell_len, wave_speed,
+        final_time, interval_len, quadrature,
+        init_cond)
+    approx = calculate_approximate_solution(
+        projection[:, 1:-1], quadrature.get_eval_points(),
+        polynomial_degree, basis.get_basis_vector())
+
+    # Plot multiwavelet solution (fine and coarse grid)
+    if coarse_projection is not None:
+        coarse_cell_len = 2*cell_len
+        coarse_mesh = np.arange(left_bound - (0.5*coarse_cell_len),
+                                right_bound + (1.5*coarse_cell_len),
+                                coarse_cell_len)
+
+        # Plot exact and approximate solutions for coarse mesh
+        coarse_grid, coarse_exact = calculate_exact_solution(
+            coarse_mesh[1:-1], coarse_cell_len, wave_speed,
+            final_time, interval_len, quadrature,
+            init_cond)
+        coarse_approx = calculate_approximate_solution(
+            coarse_projection, quadrature.get_eval_points(),
+            polynomial_degree, basis.get_basis_vector())
+        plot_solution_and_approx(
+            coarse_grid, coarse_exact, coarse_approx, colors['coarse_exact'],
+            colors['coarse_approx'])
+
+        # Plot multiwavelet details
+        num_coarse_grid_cells = num_grid_cells//2
+        plot_details(projection[:, 1:-1], mesh[2:-2], coarse_projection,
+                     basis.get_basis_vector(),
+                     basis.get_wavelet_vector(), multiwavelet_coeffs,
+                     num_coarse_grid_cells,
+                     polynomial_degree)
+
+        plot_solution_and_approx(grid, exact, approx,
+                                 colors['fine_exact'],
+                                 colors['fine_approx'])
+        plt.legend(['Exact (Coarse)', 'Approx (Coarse)', 'Exact (Fine)',
+                    'Approx (Fine)'])
+    # Plot regular solution (fine grid)
+    else:
+        plot_solution_and_approx(grid, exact, approx, colors['exact'],
+                                 colors['approx'])
+        plt.legend(['Exact', 'Approx'])
+
+    # Calculate errors
+    pointwise_error = np.abs(exact-approx)
+    max_error = np.max(pointwise_error)
+
+    # Plot errors
+    plot_semilog_error(grid, pointwise_error)
+    plot_error(grid, exact, approx)
+
+    print('p =', polynomial_degree)
+    print('N =', num_grid_cells)
+    print('maximum error =', max_error)
+
+
+def _check_colors(colors: dict) -> dict:
+    """Checks plot colors.
+
+    Checks whether colors for plots were given and sets them if required.
+
+    Parameters
+    ----------
+    colors: dict
+        Dictionary containing color strings for plotting.
+
+    Returns
+    -------
+    dict
+        Dictionary containing color strings for plotting.
+
+
+    """
+    # Set colors for general plots
+    colors['exact'] = colors.get('exact', 'k-')
+    colors['approx'] = colors.get('approx', 'y')
+
+    # Set colors for multiwavelet plots
+    colors['fine_exact'] = colors.get('fine_exact', 'k-.')
+    colors['fine_approx'] = colors.get('fine_approx', 'b-.')
+    colors['coarse_exact'] = colors.get('coarse_exact', 'k-')
+    colors['coarse_approx'] = colors.get('coarse_approx', 'y')
+
+    return colors
diff --git a/workflows/ANN_training.smk b/workflows/ANN_training.smk
index af8085efed9435e0011c0b1a8017a7430f7e8f8f..f4ccfc40753da3178ffecfdf141e03dd11c30384 100644
--- a/workflows/ANN_training.smk
+++ b/workflows/ANN_training.smk
@@ -3,6 +3,7 @@ import numpy as np
 
 import ANN_Training
 from ANN_Training import *
+from Plotting import plot_evaluation_results
 
 configfile: 'config.yaml'
 
diff --git a/workflows/approximation.smk b/workflows/approximation.smk
index f5cc877d4c56d7b9075d18de37f633610aaeff6e..eb8a84c17ad7f67c964a110afedfe0cb212f24c5 100644
--- a/workflows/approximation.smk
+++ b/workflows/approximation.smk
@@ -2,10 +2,10 @@ import sys
 import time
 import numpy as np
 
-import Troubled_Cell_Detector
 import Initial_Condition
 import Quadrature
-from DG_Approximation import DGScheme, plot_approximation_results
+from DG_Approximation import DGScheme
+from Plotting import plot_approximation_results
 from Basis_Function import OrthonormalLegendre
 
 configfile: 'config.yaml'