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 0b146be91..3e8f43a19 100644
--- a/discopop_library/discopop_optimizer/OptimizationGraph.py
+++ b/discopop_library/discopop_optimizer/OptimizationGraph.py
@@ -138,20 +138,32 @@ def __init__(
                     FreeSymbolDistribution, symbol_distribution
                 )
 
-        # add function symbols to list of substitutions
-        # TODO NOTE THIS IS JUST A DUMMY!!!
-        #  collect substitutions
+        # 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[
+                function
+            ] = sequential_complete_performance_models[function][0]
+
+        # add function symbols to list of substitutions
+        # collect substitutions
+        for function in experiment.selected_paths_per_function:
             # register substitution
             substitutions[
                 cast(Symbol, function.sequential_costs)
-            ] = sequential_complete_performance_models[function][0][0].sequential_costs
+            ] = experiment.selected_paths_per_function[function][0].sequential_costs
             substitutions[
                 cast(Symbol, function.parallelizable_costs)
-            ] = sequential_complete_performance_models[function][0][0].parallelizable_costs
+            ] = experiment.selected_paths_per_function[function][0].parallelizable_costs
 
         print("FUNCTION SUBSTITUTIONS", substitutions)
-        # todo MAY NOT STAY AS IS! MUST BE DYNAMIC
 
         # perform iterative substitutions
         modification_found = True
@@ -159,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:
@@ -188,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:
@@ -238,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:
@@ -263,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/Variables/Experiment.py b/discopop_library/discopop_optimizer/Variables/Experiment.py
index 0097d3ab2..5a5e68e30 100644
--- a/discopop_library/discopop_optimizer/Variables/Experiment.py
+++ b/discopop_library/discopop_optimizer/Variables/Experiment.py
@@ -63,6 +63,7 @@ class Experiment(object):
     detection_result: DetectionResult
 
     function_models: Dict[FunctionRoot, List[Tuple[CostModel, ContextObject, str]]]
+    selected_paths_per_function: Dict[FunctionRoot, Tuple[CostModel, ContextObject]]
 
     optimization_graph: nx.DiGraph
 
@@ -100,6 +101,7 @@ def __init__(
             self.register_free_symbol(free_symbol, value_suggestion)
 
         self.function_models = dict()
+        self.selected_paths_per_function = dict()
 
         self.compile_check_command = arguments["--compile-command"]
 
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 3c941085e..46373dbb7 100644
--- a/discopop_library/discopop_optimizer/classes/nodes/Workload.py
+++ b/discopop_library/discopop_optimizer/classes/nodes/Workload.py
@@ -97,6 +97,16 @@ 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)
+
         return cm
 
     def __get_costs_of_function_call(self, experiment, all_function_nodes) -> CostModel:
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 41cd6b4f8..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,
 )
@@ -39,7 +42,7 @@ def show_options(
     free_symbol_distributions: Dict[Symbol, FreeSymbolDistribution],
     function_root: FunctionRoot,
     parent_frame: tkinter.Frame,
-    spawned_windows: List[tkinter.Toplevel] = [],
+    spawned_windows: List[tkinter.Toplevel],
     window_title=None,
 ) -> List[Tuple[CostModel, ContextObject, str]]:
     """Shows a tkinter table to browse and plot models"""
@@ -64,6 +67,31 @@ def show_options(
         header_cols.append(e)
     rows.append(header_cols)
 
+    label1 = Entry(root, relief=RIDGE)
+    label1.grid(row=1, column=0, sticky=NSEW)
+    label1.insert(END, "Current selection:")
+    label1.configure(state=DISABLED, disabledforeground="black")
+
+    label2 = Entry(root, relief=RIDGE)
+    label2.grid(row=1, column=1, sticky=NSEW)
+    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",
@@ -154,6 +182,26 @@ 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)
+            label2.delete(0, END)
+            label2.insert(0, str(cm.path_decisions))
+            label2.configure(state=DISABLED)
+
+        update_selection_button = Button(
+            options_field,
+            text="Update selection",
+            command=lambda opt=option, opt_name=option_name, ctx=context: __update_selection(  # type: ignore
+                opt, ctx
+            ),
+        )
+        update_selection_button.grid(row=0, column=3)
+
     root.mainloop()
 
     return options
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: