diff --git a/discopop_library/discopop_optimizer/CostModels/CostModel.py b/discopop_library/discopop_optimizer/CostModels/CostModel.py
index 08d85e1c7..784880749 100644
--- a/discopop_library/discopop_optimizer/CostModels/CostModel.py
+++ b/discopop_library/discopop_optimizer/CostModels/CostModel.py
@@ -22,6 +22,8 @@ class CostModel(object):
     identifier: str
     parallelizable_costs: Expr
     sequential_costs: Expr
+    raw_parallelizable_costs: Optional[Expr]
+    raw_sequential_costs: Optional[Expr]
     free_symbol_ranges: Dict[Symbol, Tuple[float, float]]
     free_symbol_distributions: Dict[Symbol, FreeSymbolDistribution]
     symbol_value_suggestions: Dict[Symbol, Expr]
@@ -51,6 +53,8 @@ def __init__(
         self.identifier = identifier
         self.parallelizable_costs = parallelizable_costs
         self.sequential_costs = sequential_costs
+        self.raw_parallelizable_costs = None
+        self.raw_sequential_costs = None
 
     def __str__(self):
         return str(self.parallelizable_costs) + "\n" + str(self.sequential_costs)
diff --git a/discopop_library/discopop_optimizer/OptimizationGraph.py b/discopop_library/discopop_optimizer/OptimizationGraph.py
index f2f1b56ef..3e8f43a19 100644
--- a/discopop_library/discopop_optimizer/OptimizationGraph.py
+++ b/discopop_library/discopop_optimizer/OptimizationGraph.py
@@ -138,6 +138,14 @@ def __init__(
                     FreeSymbolDistribution, symbol_distribution
                 )
 
+        # save raw cost models for sequential functions
+        for function in sequential_complete_performance_models:
+            for model, ctx in sequential_complete_performance_models[function]:
+                if model.raw_sequential_costs is None:
+                    model.raw_sequential_costs = model.sequential_costs
+                if model.raw_parallelizable_costs is None:
+                    model.raw_parallelizable_costs = model.parallelizable_costs
+
         # by default, select the sequential version of each function for substitution
         for function in sequential_complete_performance_models:
             experiment.selected_paths_per_function[
@@ -163,7 +171,12 @@ def __init__(
             print("SUBSTITUTION LOOP")
             modification_found = False
             for idx, function in enumerate(sequential_complete_performance_models):
-                for midx, model in enumerate(sequential_complete_performance_models[function]):
+                for model, ctx in sequential_complete_performance_models[function]:
+                    # save raw cost models
+                    if model.raw_sequential_costs is None:
+                        model.raw_sequential_costs = model.sequential_costs
+                    if model.raw_parallelizable_costs is None:
+                        model.raw_parallelizable_costs = model.parallelizable_costs
                     # apply substitution to parallelizable costs
                     tmp_model = model.parallelizable_costs.subs(substitutions)
                     if tmp_model != model.parallelizable_costs:
@@ -192,6 +205,11 @@ def __init__(
             for idx, function in enumerate(sequential_complete_performance_models):
                 for midx, pair in enumerate(sequential_complete_performance_models[function]):
                     model, context = pair
+                    # save raw cost models
+                    if model.raw_sequential_costs is None:
+                        model.raw_sequential_costs = model.sequential_costs
+                    if model.raw_parallelizable_costs is None:
+                        model.raw_parallelizable_costs = model.parallelizable_costs
                     # apply substitution to parallelizable costs
                     tmp_model = model.parallelizable_costs.subs(substitutions)
                     if tmp_model != model.parallelizable_costs:
@@ -242,6 +260,11 @@ def __init__(
             for idx, function in enumerate(locally_optimized_models):
                 for midx, pair in enumerate(locally_optimized_models[function]):
                     model, context = pair
+                    # save raw cost models
+                    if model.raw_sequential_costs is None:
+                        model.raw_sequential_costs = model.sequential_costs
+                    if model.raw_parallelizable_costs is None:
+                        model.raw_parallelizable_costs = model.parallelizable_costs
                     # apply substitution to parallelizable costs
                     tmp_model = model.parallelizable_costs.subs(substitutions)
                     if tmp_model != model.parallelizable_costs:
@@ -267,6 +290,11 @@ def __init__(
             for idx, function in enumerate(exhaustive_performance_models):
                 for midx, pair in enumerate(exhaustive_performance_models[function]):
                     model, context = pair
+                    # save raw cost models
+                    if model.raw_sequential_costs is None:
+                        model.raw_sequential_costs = model.sequential_costs
+                    if model.raw_parallelizable_costs is None:
+                        model.raw_parallelizable_costs = model.parallelizable_costs
                     # apply substitution to parallelizable costs
                     tmp_model = model.parallelizable_costs.subs(substitutions)
                     if tmp_model != model.parallelizable_costs:
diff --git a/discopop_library/discopop_optimizer/classes/nodes/FunctionRoot.py b/discopop_library/discopop_optimizer/classes/nodes/FunctionRoot.py
index 71ec35ced..ddbf1f352 100644
--- a/discopop_library/discopop_optimizer/classes/nodes/FunctionRoot.py
+++ b/discopop_library/discopop_optimizer/classes/nodes/FunctionRoot.py
@@ -60,4 +60,11 @@ def get_cost_model(self, experiment, all_function_nodes) -> CostModel:
         cm.parallelizable_costs = cm.parallelizable_costs.subs({Expr(Integer(0)): Integer(0)})
         cm.sequential_costs = cm.sequential_costs.subs({Expr(Integer(0)): Integer(0)})
 
+        if cm.raw_parallelizable_costs is not None:
+            cm.raw_parallelizable_costs = cm.raw_parallelizable_costs.subs(
+                {Expr(Integer(0)): Integer(0)}
+            )
+        if cm.raw_sequential_costs is not None:
+            cm.raw_sequential_costs = cm.raw_sequential_costs.subs({Expr(Integer(0)): Integer(0)})
+
         return cm
diff --git a/discopop_library/discopop_optimizer/classes/nodes/Loop.py b/discopop_library/discopop_optimizer/classes/nodes/Loop.py
index b68878e75..72ded9bc7 100644
--- a/discopop_library/discopop_optimizer/classes/nodes/Loop.py
+++ b/discopop_library/discopop_optimizer/classes/nodes/Loop.py
@@ -105,6 +105,13 @@ def get_cost_model(self, experiment, all_function_nodes) -> CostModel:
         cm.parallelizable_costs = cm.parallelizable_costs.subs({Expr(Integer(0)): Integer(0)})
         cm.sequential_costs = cm.sequential_costs.subs({Expr(Integer(0)): Integer(0)})
 
+        if cm.raw_sequential_costs is not None:
+            cm.raw_sequential_costs = cm.raw_sequential_costs.subs({Expr(Integer(0)): Integer(0)})
+        if cm.raw_parallelizable_costs is not None:
+            cm.raw_parallelizable_costs = cm.raw_parallelizable_costs.subs(
+                {Expr(Integer(0)): Integer(0)}
+            )
+
         return cm
 
     def register_child(self, other, experiment, all_function_nodes):
diff --git a/discopop_library/discopop_optimizer/classes/nodes/Workload.py b/discopop_library/discopop_optimizer/classes/nodes/Workload.py
index 702a07888..46373dbb7 100644
--- a/discopop_library/discopop_optimizer/classes/nodes/Workload.py
+++ b/discopop_library/discopop_optimizer/classes/nodes/Workload.py
@@ -97,6 +97,13 @@ def get_cost_model(self, experiment, all_function_nodes) -> CostModel:
         cm.parallelizable_costs = cm.parallelizable_costs.subs({Expr(Integer(0)): Integer(0)})
         cm.sequential_costs = cm.sequential_costs.subs({Expr(Integer(0)): Integer(0)})
 
+        if cm.raw_sequential_costs is not None:
+            cm.raw_sequential_costs = cm.raw_sequential_costs.subs({Expr(Integer(0)): Integer(0)})
+        if cm.raw_parallelizable_costs is not None:
+            cm.raw_parallelizable_costs = cm.raw_parallelizable_costs.subs(
+                {Expr(Integer(0)): Integer(0)}
+            )
+
         print("CM: ")
         print(cm)
 
diff --git a/discopop_library/discopop_optimizer/gui/plotting/CostModels.py b/discopop_library/discopop_optimizer/gui/plotting/CostModels.py
index 916974fee..67973d8e1 100644
--- a/discopop_library/discopop_optimizer/gui/plotting/CostModels.py
+++ b/discopop_library/discopop_optimizer/gui/plotting/CostModels.py
@@ -5,16 +5,18 @@
 # This software may be modified and distributed under the terms of
 # the 3-Clause BSD License.  See the LICENSE file in the package base
 # directory for details.
-from typing import List, Dict, Tuple, Optional
+import copy
+from typing import List, Dict, Tuple, Optional, cast
 
 import numpy as np
 from matplotlib import pyplot as plt  # type: ignore
 import matplotlib
 from spb import plot3d, MB, plot  # type: ignore
-from sympy import Symbol
+from sympy import Symbol, Expr
 import sympy
 
 from discopop_library.discopop_optimizer.CostModels.CostModel import CostModel
+from discopop_library.discopop_optimizer.Variables.Experiment import Experiment
 
 
 def plot_CostModels(
@@ -54,6 +56,93 @@ def plot_CostModels(
         print("Plotiting not supported for", len(sorted_free_symbols), "free symbols!")
 
 
+def plot_CostModels_using_function_path_selections(
+    experiment: Experiment,
+    models: List[CostModel],
+    sorted_free_symbols: List[Symbol],
+    free_symbol_ranges: Dict[Symbol, Tuple[float, float]],
+    labels: Optional[List[str]] = None,
+    title: Optional[str] = None,
+    super_title: Optional[str] = None,
+):
+    print("PLOTTING: ")
+    for m in models:
+        print(m.raw_sequential_costs)
+        print(m.raw_parallelizable_costs)
+    print()
+
+    # apply selected substitutions
+    # collect substitutions
+    local_substitutions = copy.deepcopy(experiment.substitutions)
+    for function in experiment.selected_paths_per_function:
+        # register substitution
+        local_substitutions[
+            cast(Symbol, function.sequential_costs)
+        ] = experiment.selected_paths_per_function[function][0].sequential_costs
+        local_substitutions[
+            cast(Symbol, function.parallelizable_costs)
+        ] = experiment.selected_paths_per_function[function][0].parallelizable_costs
+
+    print("LOCAL FUNCTION SUBSTITUTIONS", local_substitutions)
+
+    # prepare models by loading raw costs
+    for model in models:
+        model.sequential_costs = cast(Expr, model.raw_sequential_costs)
+        model.parallelizable_costs = cast(Expr, model.raw_parallelizable_costs)
+
+    # perform iterative substitutions
+    modification_found = True
+    while modification_found:
+        print("LOCAL SUBSTITUTION LOOP")
+        modification_found = False
+        for model in models:
+            # apply substitution to parallelizable costs
+            tmp_model = model.parallelizable_costs.subs(local_substitutions)
+            if tmp_model != model.parallelizable_costs:
+                modification_found = True
+            model.parallelizable_costs = tmp_model
+
+            # apply substitutions to sequential costs
+            tmp_model = model.sequential_costs.subs(local_substitutions)
+            if tmp_model != model.sequential_costs:
+                modification_found = True
+            model.sequential_costs = model.sequential_costs.subs(local_substitutions)
+
+    print("PLOTTING AFTER SUBSTITUTION: ")
+    for m in models:
+        print(m.sequential_costs)
+        print(m.parallelizable_costs)
+    print()
+
+    if len(sorted_free_symbols) == 2:
+        __3d_plot(
+            models,
+            sorted_free_symbols,
+            free_symbol_ranges,
+            labels=labels,
+            title=str(title) + str(super_title) if super_title is not None else title,
+        )
+    elif len(sorted_free_symbols) == 1:
+        __2d_plot(
+            models,
+            sorted_free_symbols,
+            free_symbol_ranges,
+            labels=labels,
+            title=str(title) + str(super_title) if super_title is not None else title,
+        )
+    elif len(sorted_free_symbols) == 0:
+        __1d_plot(
+            models,
+            sorted_free_symbols,
+            free_symbol_ranges,
+            labels=labels,
+            title=title,
+            super_title=super_title,
+        )
+    else:
+        print("Plotiting not supported for", len(sorted_free_symbols), "free symbols!")
+
+
 __unique_plot_id = 0
 
 
diff --git a/discopop_library/discopop_optimizer/gui/presentation/OptionTable.py b/discopop_library/discopop_optimizer/gui/presentation/OptionTable.py
index 7857130d9..c03a39c10 100644
--- a/discopop_library/discopop_optimizer/gui/presentation/OptionTable.py
+++ b/discopop_library/discopop_optimizer/gui/presentation/OptionTable.py
@@ -19,7 +19,10 @@
 from discopop_library.discopop_optimizer.classes.context.ContextObject import ContextObject
 from discopop_library.discopop_optimizer.classes.enums.Distributions import FreeSymbolDistribution
 from discopop_library.discopop_optimizer.classes.nodes.FunctionRoot import FunctionRoot
-from discopop_library.discopop_optimizer.gui.plotting.CostModels import plot_CostModels
+from discopop_library.discopop_optimizer.gui.plotting.CostModels import (
+    plot_CostModels,
+    plot_CostModels_using_function_path_selections,
+)
 from discopop_library.discopop_optimizer.gui.presentation.ChoiceDetails import (
     display_choices_for_model,
 )
@@ -74,6 +77,21 @@ def show_options(
     label2.insert(END, str(experiment.selected_paths_per_function[function_root][0].path_decisions))
     label2.configure(state=DISABLED, disabledforeground="black")
 
+    plot_button = Button(
+        root,
+        text="Plot using selections",
+        command=lambda: plot_CostModels_using_function_path_selections(  # type: ignore
+            experiment,
+            [experiment.selected_paths_per_function[function_root][0]],  # type: ignore
+            sorted_free_symbols,
+            free_symbol_ranges,
+            [function_root.name],
+            title=function_root.name,
+            super_title=function_root.name,
+        ),
+    )
+    plot_button.grid(row=1, column=2, sticky=NSEW)
+
     Button(
         root,
         text="Plot All",
@@ -165,6 +183,9 @@ def show_options(
         export_code_button.grid(row=0, column=2)
 
         def __update_selection(cm, ctx):
+            # use raw models for selection updates
+            cm.parallelizable_costs = cm.raw_parallelizable_costs
+            cm.sequential_costs = cm.raw_sequential_costs
             experiment.selected_paths_per_function[function_root] = (cm, ctx)
             # update displayed value
             label2.configure(state=NORMAL)
diff --git a/discopop_library/discopop_optimizer/utilities/optimization/GlobalOptimization/RandomSamples.py b/discopop_library/discopop_optimizer/utilities/optimization/GlobalOptimization/RandomSamples.py
index 552c587c8..2c70c16f8 100644
--- a/discopop_library/discopop_optimizer/utilities/optimization/GlobalOptimization/RandomSamples.py
+++ b/discopop_library/discopop_optimizer/utilities/optimization/GlobalOptimization/RandomSamples.py
@@ -59,6 +59,11 @@ def find_quasi_optimal_using_random_samples(
         print("\tApplying substitutions...")
         print("\t" + str(substitutions))
     for model, context in random_paths:
+        # save raw cost models
+        if model.raw_sequential_costs is None:
+            model.raw_sequential_costs = model.sequential_costs
+        if model.raw_parallelizable_costs is None:
+            model.raw_parallelizable_costs = model.parallelizable_costs
         # apply substitutions iteratively
         modification_found = True
         while modification_found:
diff --git a/discopop_library/discopop_optimizer/utilities/optimization/LocalOptimization/TopDown.py b/discopop_library/discopop_optimizer/utilities/optimization/LocalOptimization/TopDown.py
index a57172ae8..424e2d0bf 100644
--- a/discopop_library/discopop_optimizer/utilities/optimization/LocalOptimization/TopDown.py
+++ b/discopop_library/discopop_optimizer/utilities/optimization/LocalOptimization/TopDown.py
@@ -90,6 +90,12 @@ def get_locally_optimized_models(
                 modification_found = False
                 for decision, pair in decision_models:
                     model, context = pair
+                    # save raw cost models
+                    if model.raw_sequential_costs is None:
+                        model.raw_sequential_costs = model.sequential_costs
+                    if model.raw_parallelizable_costs is None:
+                        model.raw_parallelizable_costs = model.parallelizable_costs
+
                     # apply substitutions to parallelizable costs
                     tmp_model = model.parallelizable_costs.subs(substitutions)
                     if tmp_model != model.parallelizable_costs: