diff --git a/kadmos/graph/graph_data.py b/kadmos/graph/graph_data.py
index 0e53adccc554c3df2e603c32569fa743dfac5fc7..bb68df26fdc3ed3a501c39bb4506938deed344a1 100644
--- a/kadmos/graph/graph_data.py
+++ b/kadmos/graph/graph_data.py
@@ -478,6 +478,12 @@ class DataGraph(KadmosGraph):
             assert not self.find_all_nodes(subcategory='all problematic variables'), \
                 'Graph still has problematic variables.'
 
+        # Check for partitions
+        if 'problem_formulation' in self.graph and 'partitions' in self.graph['problem_formulation']:
+            partitions = self.graph['problem_formulation']['partitions']
+        else:
+            partitions = None
+
         # Get function graph
         function_graph = self.get_function_graph()
         function_graph.remove_edges_from(function_graph.selfloop_edges())
@@ -486,22 +492,39 @@ class DataGraph(KadmosGraph):
         function_graph.add_node('super_node', category='function')
         coupled_functions = []
 
-        # As long as not all coupled functions are merged into the super node:
-        while not nx.is_directed_acyclic_graph(function_graph):
-            # Find a cycle
-            cycle = nx.find_cycle(function_graph)
-
-            # Find the functions in the cycle
-            functions_in_cycle = set()
-            functions_in_cycle.update(function_id for edges in cycle for function_id in edges)
-            functions_in_cycle = list(functions_in_cycle)
-
-            # Merge the coupled functions in the super node
-            for function_id in functions_in_cycle:
-                if function_id != 'super_node':
-                    coupled_functions.append(function_id)
-                    function_graph = nx.contracted_nodes(function_graph, 'super_node', function_id)
-                    function_graph.remove_edges_from(function_graph.selfloop_edges())
+        # If partitions check if all coupled nodes are assigned to a partition
+        if partitions:
+            functions_in_partitions = []
+            for partition in partitions:
+                nodes = list(partitions[partition])
+                functions_in_partitions.extend(nodes)
+            for function_id in functions_in_partitions:
+                function_graph = nx.contracted_nodes(function_graph, 'super_node', function_id)
+            function_graph.remove_edges_from(function_graph.selfloop_edges())
+            if not nx.is_directed_acyclic_graph(function_graph):
+                cycle = nx.find_cycle(function_graph)
+                functions_in_cycle = set()
+                functions_in_cycle.update(function_id for edges in cycle for function_id in edges)
+                functions_in_cycle.remove('super_node')
+                assert nx.is_directed_acyclic_graph(function_graph), 'Coupled functions {} should be added to a ' \
+                                                                     'partition'.format(list(functions_in_cycle))
+        else:
+            # As long as not all coupled functions are merged into the super node:
+            while not nx.is_directed_acyclic_graph(function_graph):
+                # Find a cycle
+                cycle = nx.find_cycle(function_graph)
+
+                # Find the functions in the cycle
+                functions_in_cycle = set()
+                functions_in_cycle.update(function_id for edges in cycle for function_id in edges)
+                functions_in_cycle = list(functions_in_cycle)
+
+                # Merge the coupled functions in the super node
+                for function_id in functions_in_cycle:
+                    if function_id != 'super_node':
+                        coupled_functions.append(function_id)
+                        function_graph = nx.contracted_nodes(function_graph, 'super_node', function_id)
+                        function_graph.remove_edges_from(function_graph.selfloop_edges())
 
         # Find a topological function order
         function_order = list(nx.topological_sort(function_graph))
@@ -511,8 +534,17 @@ class DataGraph(KadmosGraph):
         pre_coupling_functions_order = self.sort_nodes_for_process(pre_coupling_functions)
 
         # Sort coupled functions to minimize feedback
-        coupled_functions_order = self.minimize_feedback(coupled_functions, method, multi_start=multi_start)
-        coupled_functions_order = self.sort_nodes_for_process(coupled_functions_order)
+        if partitions:
+            coupled_functions_order = []
+            for partition in partitions:
+                nodes = self.minimize_feedback(list(partitions[partition]), method, multi_start=multi_start)
+                nodes = self.sort_nodes_for_process(nodes)
+                partitions[partition] = nodes
+                coupled_functions_order.extend(nodes)
+            self.graph['problem_formulation']['partitions'] = partitions
+        else:
+            coupled_functions_order = self.minimize_feedback(coupled_functions, method, multi_start=multi_start)
+            coupled_functions_order = self.sort_nodes_for_process(coupled_functions_order)
 
         # Get post-coupling functions and sort
         post_coupling_functions = function_order[function_order.index('super_node') + 1:]
@@ -2388,16 +2420,193 @@ class FundamentalProblemGraph(DataGraph, KeChainMixin):
             self.graph['problem_formulation']['doe_settings']['doe_seed'] = doe_settings['doe_seed']
             self.graph['problem_formulation']['doe_settings']['doe_runs'] = doe_settings['doe_runs']
 
-    def select_number_of_partitions(self, partition_range, edge_weight='coupling_strength', node_weights=None,
-                                    plot_pareto_front=False):
+    def partition_graph(self, n_parts, include_run_time=False, tpwgts=None, recursive=True, contig=False):
+        """Partition a graph using the Metis algorithm (http://glaros.dtc.umn.edu/gkhome/metis/metis/overview).
+
+        Note that partitioning can only be performed on undirected graphs. Therefore every graph input is translated
+        into an undirected graph.
+
+        :param n_parts: number of partitions requested (algorithm might provide less)
+        :type n_parts: int
+        :param include_run_time:
+        :type include_run_time: bool
+        :param tpwgts: list of target partition weights
+        :type tpwgts: list
+        :param recursive: option to use the recursive bisection or k-way partitioning algorithm
+        :type recursive: bool
+        :param contig: Metis option
+        :type contig: bool
+        """
+
+        import metis
+
+        # Input assertions
+        assert 'function_ordering' in self.graph['problem_formulation'], 'Function ordering is missing'
+        coupled_nodes = self.graph['problem_formulation']['function_ordering'][self.FUNCTION_ROLES[1]]
+        if include_run_time:
+            for node in coupled_nodes:
+                assert 'run_time' in self.nodes[node]['performance_info'], 'Run time missing for function ' \
+                                                                           '{}'.format(node)
+
+        # Get subgraph
+        subgraph = self.get_subgraph_by_function_nodes(coupled_nodes)
+        graph = subgraph.deepcopy()
+        coupling_dict = graph.get_coupling_dictionary()
+
+        calculated_partitions = dict()
+        iteration = 0
+
+        while iteration < 5:
+            iteration += 1
+
+            # Get undirected graph
+            g_und = graph.to_undirected()
+
+            # Add run time to the nodes of the undirected graph
+            for node in g_und.nodes():
+                if g_und.nodes[node]['category'] == 'variable':
+                    g_und.nodes[node]['run_time'] = 0
+                elif g_und.nodes[node]['category'] == 'function':
+                    g_und.nodes[node]['run_time'] = g_und.nodes[node]['performance_info']['run_time'] if \
+                        include_run_time else 1
+            g_und.graph['node_weight_attr'] = 'run_time'
+
+            # Partition graph
+            (edgecuts, parts) = metis.part_graph(g_und, n_parts, tpwgts=tpwgts, recursive=recursive, contig=contig)
+
+            # Get partition dict
+            partitions = dict()
+            for part in range(n_parts):
+                nodes = []
+                # Get function nodes in partition
+                for idx, node in enumerate(g_und.nodes):
+                    if parts[idx] == part and graph.nodes[node]['category'] == 'function':
+                        nodes.extend(node.split('--') if '--' in node else [node])
+                # Minimize feedback within the partition
+                nodes = self.minimize_feedback(nodes, 'single-swap')
+                nodes = self.sort_nodes_for_process(nodes)
+                # Add nodes to the partition dict
+                partitions[part+1] = nodes
+
+            # Reset graph
+            graph = subgraph.deepcopy()
+
+            # Evaluate the properties of the partitioning
+            partition_variables, system_variables, runtime = graph.get_partition_info(partitions,
+                                                                                      include_run_time=include_run_time)
+
+            # Merge nodes that can be merged based on process and calculate runtime of each partition
+            for partition in partitions:
+                nodes = list(partitions[partition])
+                while nodes:
+                    merge_nodes, run_times = [], []
+                    for idx, node in enumerate(nodes):
+                        if not set(nodes[:idx]).intersection(coupling_dict[node]):
+                            merge_nodes.append(node)
+                            run_times.append(self.nodes[node]['performance_info']['run_time'] if include_run_time else
+                                             1)
+                        else:
+                            break
+                    if len(merge_nodes) > 1:
+                        new_node = '--'.join(merge_nodes)
+                        try:
+                            graph = graph.merge_parallel_functions(merge_nodes, new_label=new_node)
+                            graph.nodes[new_node]['performance_info'] = {'run_time': max(run_times)}
+                        except AssertionError:
+                            print 'Could not merge nodes {}'.format(merge_nodes)
+                    for node in merge_nodes:
+                        nodes.pop(nodes.index(node))
+
+            # Remember current partition
+            calculated_partitions[iteration] = dict()
+            calculated_partitions[iteration]['partitions'] = partitions
+            calculated_partitions[iteration]['variables'] = len(system_variables) + sum([len(variables) for variables
+                                                                                         in partition_variables])
+            calculated_partitions[iteration]['run_time'] = runtime
+
+        # Determine the best solution todo
+        best_solution = 1
+
+        # Return best option
+        partitions = calculated_partitions[best_solution]['partitions']
+
+        return partitions
+
+    def get_partition_info(self, partitions, include_run_time=False):
+        """ Function to get the number of feedback variables in the partitions and the number system variables (feedback
+        and feedforward variables between partitions)
+
+        :param partitions: dictionary which indicates which nodes are in which partition
+        :type partitions: dict
+        :param include_run_time:
+        :type include_run_time: bool
+        :return partition_variables:
+        :rtype partition_variables: list of sets
+        :return system_variables:
+        :rtype system_variables: set
+        """
+
+        # Get complete function order of nodes in the partitions
+        function_order = []
+        for partition in partitions:
+            function_order.extend(partitions[partition])
+
+        # Input assertions
+        if include_run_time:
+            for node in function_order:
+                assert 'run_time' in self.nodes[node]['performance_info'], 'Run time missing for function ' \
+                                                                           '{}'.format(node)
+
+        # Get coupling dictionary
+        coupling_dict = self.get_coupling_dictionary()
+
+        # For each node in the partitions check whether its input comes from the same partition, another partition or
+        # neither
+        partition_variables = [set() for _ in partitions]
+        system_variables = set()
+        for partition in partitions:
+            for idx, target in enumerate(partitions[partition]):
+                for source in coupling_dict[target]:
+                    if source in partitions[partition][idx+1:]:
+                        paths = nx.all_shortest_paths(self, source, target)
+                        for path in paths:
+                            partition_variables[partition-1].update([path[1].split('/')[-1]])
+                    elif source in function_order and source not in partitions[partition]:
+                        paths = nx.all_shortest_paths(self, source, target)
+                        for path in paths:
+                            system_variables.update([path[1].split('/')[-1]])
+
+        # Calculate run time of each partition
+        run_time_partitions = []
+        for partition in partitions:
+            nodes = list(partitions[partition])
+            run_time_partition = 0
+            while nodes:
+                parallel_nodes = []
+                run_times = []
+                for idx, node in enumerate(nodes):
+                    if not set(nodes[:idx]).intersection(coupling_dict[node]):
+                        parallel_nodes.append(node)
+                        if include_run_time:
+                            run_times.append(self.nodes[node]['performance_info']['run_time'])
+                        else:
+                            run_times.append(1)
+                    else:
+                        break
+                run_time_partition += max(run_times)
+                for node in parallel_nodes:
+                    nodes.pop(nodes.index(node))
+            run_time_partitions.append(run_time_partition)
+
+        return partition_variables, system_variables, run_time_partitions
+
+    def select_number_of_partitions(self, partition_range, include_run_time=False, plot_pareto_front=False):
         """ Function to evaluate the properties of different number of partitions and to select the best one.
 
         :param partition_range: range of number of partitions that need to be evaluated
         :type partition_range: list
-        :param edge_weight: Edge attribute that will be used as edge weight in the partitioning algorithm
-        :type edge_weight: str
-        :param node_weights: Node attributes that will be used as node weights in the partitioning algorithm
-        :type node_weights: list
+        :param include_run_time:
+        :type include_run_time:
         :param plot_pareto_front: Option to plot the characteristics of different number of partitions
         :type plot_pareto_front: bool
         :return:
@@ -2405,91 +2614,29 @@ class FundamentalProblemGraph(DataGraph, KeChainMixin):
 
         # Input assertions
         assert 'function_ordering' in self.graph['problem_formulation'], 'Function ordering is missing'
-
-        # Get node functions
-        function_ordering = self.graph['problem_formulation']['function_ordering']
-        precoupled_nodes = function_ordering[self.FUNCTION_ROLES[0]]
-        coupled_nodes = function_ordering[self.FUNCTION_ROLES[1]]
-        postcoupled_nodes = function_ordering[self.FUNCTION_ROLES[2]]
-
-        # Get subgraph of nodes that needs to be partitioned
-        subgraph = self.get_subgraph_by_function_nodes(coupled_nodes)
+        coupled_nodes = self.graph['problem_formulation']['function_ordering'][self.FUNCTION_ROLES[1]]
+        if include_run_time:
+            for node in coupled_nodes:
+                assert 'run_time' in self.nodes[node]['performance_info'], 'Run time missing for function ' \
+                                                                           '{}'.format(node)
 
         partition_info = []
         partition_results = dict()
-        coupling_dict = self.get_coupling_dictionary()
 
         for idx, n_partitions in enumerate(partition_range):
-
             # Partition graph
-            print 'Number of partitions: ', n_partitions
-            partitioned_graph = subgraph.deepcopy()
-            partitioned_graph.partition_graph_incl_process(n_partitions, recursive=True, edge_weight=edge_weight,
-                                              node_weights=node_weights)
-            partitions = partitioned_graph.graph['problem_formulation']['partitions']
-
-            function_order = []
-            max_runtime = 0
-
-            for partition in partitions:
+            partitioned_graph = self.deepcopy()
+            partitions = partitioned_graph.partition_graph(n_partitions, include_run_time=include_run_time)
 
-                # Optimize function order within the partitions
-                function_order_partition = self.minimize_feedback(partitions[partition], 'single-swap')
-                function_order_partition = self.sort_nodes_for_process(function_order_partition)
-                partitions[partition] = function_order_partition
-                function_order.extend(function_order_partition)
-                print 'Partition {}: '.format(partition), function_order_partition
-
-                # Calculate runtime for each partition based on the process and data connections
-                runtime = 0
-                nodes = list(partitions[partition])
-                if all('runtime' in self.nodes[node] for node in nodes):
-                    while nodes:
-                        executable_nodes = []
-                        runtimes = []
-                        for i, node in enumerate(nodes):
-                            if not set(nodes[:i]).intersection(coupling_dict[node]):
-                                executable_nodes.append(node)
-                                runtimes.append(self.nodes[node]['runtime'])
-                                print 'runtimes', runtimes
-                            else:
-                                break
-                        for node in executable_nodes:
-                            nodes.pop(nodes.index(node))
-                        runtime += max(runtimes)
-                        print 'runtime', runtime
-                if runtime > max_runtime:
-                    max_runtime = runtime
-
-            # Save results
-            partition_results[n_partitions] = dict()
-            partition_results[n_partitions]['partitions'] = partitions
-            partition_results[n_partitions]['function_order'] = function_order
-
-            # Get all couplings
-            coupling_matrix = self.get_coupling_matrix(node_selection=function_order)
-            couplings = np.transpose(np.nonzero(coupling_matrix))
-
-            # Sort couplings in partition and system variables
-            partition_variables = [set() for _ in range(n_partitions)]
-            system_variables = set()
-
-            for coupling in couplings:
-                source = function_order[coupling[0]]
-                target = function_order[coupling[1]]
-                paths = nx.all_shortest_paths(self, source, target)
-                for path in paths:
-                    if partitioned_graph.nodes[source]['part_id'] == partitioned_graph.nodes[target]['part_id']:
-                        if coupling[0] > coupling[1]:
-                            partition_variables[partitioned_graph.nodes[source]['part_id']-1].update([path[1].split('/')
-                                                                                                      [-1]])
-                    else:
-                        system_variables.update([path[1].split('/')[-1]])
+            # Evaluate graph
+            partition_variables, system_variables, runtime = partitioned_graph.get_partition_info(
+                partitions, include_run_time=include_run_time)
             total_var = len(system_variables) + sum([len(variables) for variables in partition_variables])
 
             # Save partition information
             partition_info.append([idx, n_partitions, [len(variables) for variables in partition_variables],
-                                   len(system_variables), total_var, max_runtime])
+                                   len(system_variables), total_var, max(runtime)])
+            partition_results[n_partitions] = partitions
 
         # Print partition information in table
         header = ['Option', '# partitions', '# feedback in partitions', '# system variables', 'Total # variables',
@@ -2521,193 +2668,102 @@ class FundamentalProblemGraph(DataGraph, KeChainMixin):
         idx = partition_range.index(int(sel[0]))
 
         # Add partition ids to graph
-        partitions = partition_results[partition_range[idx]]['partitions']
-        for node in self.nodes:
-            for partition in partitions:
-                if node in partitions[partition]:
-                    self.nodes[node]['part_id'] = partition
-
-        # Get function order
-        precoupled_nodes_order = self.sort_nodes_for_process(precoupled_nodes)
-        coupled_nodes_order = partition_results[partition_range[idx]]['function_order']
-        post_coupled_nodes_order = self.sort_nodes_for_process(postcoupled_nodes)
-        function_order = precoupled_nodes_order + coupled_nodes_order + post_coupled_nodes_order
+        partitions = partition_results[partition_range[idx]]
 
-        self.graph['problem_formulation']['partitions'] = partitions
-        self.graph['problem_formulation']['function_order'] = function_order
-
-        return
-
-    def partition_graph_incl_process(self, n_parts, edge_weight='coupling_strength', node_weights=None, tpwgts=None,
-                                     recursive=False, contig=False):
-
-        import metis
-
-        # Make a copy of the graph
-        graph = self.deepcopy()
-        new_iteration = True
-        merged_nodes = dict()
-
-        while new_iteration:
-            new_iteration = False
-            previous_graph = graph.deepcopy()
-            # Get function graph
-            function_graph = graph.get_function_graph()
-
-            # Merge feedforward and feedback between two nodes to get a correct undirected graph
-            if edge_weight:
-                remove_edges = []
-                for edge in function_graph.edges(data=True):
-                    if (edge[0], edge[1]) in remove_edges:
-                        continue
-                    elif (edge[1], edge[0]) in function_graph.edges():
-                        function_graph.edges[edge[0], edge[1]][edge_weight] += function_graph.edges[edge[1], edge[0]][
-                            edge_weight]
-                        remove_edges.append((edge[1], edge[0]))
-                for edge in remove_edges:
-                    function_graph.remove_edge(edge[0], edge[1])
-
-            # Get undirected graph
-            g_und = function_graph.to_undirected()
-
-            # Add node and edge weight attributes to undirected graph
-            if edge_weight:
-                g_und.graph['edge_weight_attr'] = edge_weight
-            if node_weights:
-                g_und.graph['node_weight_attr'] = node_weights
-
-            # Partition graph
-            (edgecuts, parts) = metis.part_graph(g_und, n_parts, tpwgts=tpwgts, recursive=recursive, contig=contig)
-
-            # Get coupling dict of current graph
-            coupling_dict = graph.get_coupling_dictionary()
-
-            partitions = dict()
-            function_order = []
-
-            for part in range(n_parts):
-                runtime = 0
-                # Get nodes in partition
-                nodes = [node for idx, node in enumerate(g_und.nodes) if parts[idx] == part]
-                # Minimize feedback
-                nodes = graph.minimize_feedback(nodes, 'single-swap')
-                # Check which nodes can be executed in parallel and thus be merged
-                partition_function_order = []
-                while nodes:
-                    merge_nodes = []
-                    runtimes = []
-                    for idx, node in enumerate(nodes):
-                        if not set(nodes[:idx]).intersection(coupling_dict[node]):
-                            merge_nodes.append(node)
-                            runtimes.append(graph.nodes[node]['runtime'])
-                    partition_function_order.extend(merge_nodes)
-                    if len(merge_nodes) > 1:
-                        new_node = '--'.join(merge_nodes)
-                        graph = graph.merge_parallel_functions(merge_nodes, new_label=new_node)
-                        graph.nodes[new_node]['runtime'] = max(runtimes)
-                        new_iteration = True
-                    for node in merge_nodes:
-                        nodes.pop(nodes.index(node))
-                    runtime += max(runtimes)
-                partitions[part+1] = partition_function_order
-                function_order.extend(partition_function_order)
-                print 'Partition {} runtime:'.format(part+1), runtime
-
-            # Add partition information to graph
-            i = 0
-            for node, data in g_und.nodes(data=True):
-                previous_graph.nodes[node]['part_id'] = parts[i] + 1
-                i += 1
-
-            # Evaluate properties of current graph:
-            coupling_matrix = previous_graph.get_coupling_matrix(node_selection=function_order)
-            couplings = np.transpose(np.nonzero(coupling_matrix))
-            partition_variables = [set() for _ in range(n_parts)]
-            system_variables = set()
-            for coupling in couplings:
-                source = function_order[coupling[0]]
-                target = function_order[coupling[1]]
-                paths = nx.all_shortest_paths(previous_graph, source, target)
-                for path in paths:
-                    if previous_graph.nodes[source]['part_id'] == previous_graph.nodes[target]['part_id']:
-                        if coupling[0] > coupling[1]:
-                            partition_variables[previous_graph.nodes[source]['part_id']-1].update([path[1].split('/')[-1]])
-                    else:
-                        system_variables.update([path[1].split('/')[-1]])
-            print 'total_variables', len(system_variables) + sum([len(variables) for variables in partition_variables])
-
-        # Add partition information to graph
-        i = 0
-        for node1, data in g_und.nodes(data=True):
-            if '--' in node1:
-                nodes = node1.split('--')
-                print nodes
-                for node2 in nodes:
-                    self.nodes[node2]['part_id'] = parts[i] + 1
-            else:
-                self.nodes[node1]['part_id'] = parts[i] + 1
-            i += 1
-
-        function_order = []
-        # Split nodes in the partition dictionary
-        for partition in partitions:
-            order = []
-            for node in partitions[partition]:
-                order.extend(node.split('--') if '--' in node else [node])
-            partitions[partition] = order
-            function_order.extend(order)
-
-        self.graph['problem_formulation']['partitions'] = partitions
-
-        return
+        return partitions
 
     def select_distributed_architecture(self):
+        """ Function for easy selection of a distributed architecture for a partitioned graph.
+
+        :return Extended problem formulation
+        """
 
+        # Check if graph is partitioned
         assert 'partitions' in self.graph['problem_formulation'], 'Graph is not partitioned. ' \
                                                                   'Distributed architecture is not possible'
         assert len(self.graph['problem_formulation']['partitions']) > 1, 'Graph must have at least two partitions ' \
                                                                          'to select a distributed architecture'
 
+        # Get graph info
         partitions = self.graph['problem_formulation']['partitions']
+        coupling_dict = self.get_coupling_dictionary()
 
-        # For each partition decide whether the partition gets its own converger or not
-        selmsg = 'Please select which partitions must have a local converger:'
-        idx_list = [partition for partition in partitions]
-        local_convergers = prompting.user_prompt_select_options(*idx_list, allow_empty=True, allow_multi=True,
-                                                                message=selmsg)
-        # TODO: add check if converger can be added
-        # For each partition decide whether the partition must be solved using Jacobi or Gauss-Seidel convergence
-        selmsg = 'Please select which partitions must be solved using a Jacobi convergence instead of Gauss-Seidel'
-        idx_list = [partition for partition in partitions]
-        jacobi_convergence = prompting.user_prompt_select_options(*idx_list, allow_empty=True, allow_multi=True,
-                                                                message=selmsg)
-
-        # Choose system architecture
-        selmsg = 'Please select system architecture:'
+        # Select system architecture
         system_architectures = ['unconverged-MDA', 'converged-MDA', 'MDF', 'IDF', 'unconverged-OPT']
         options = []
         for idx, arch in enumerate(system_architectures):
             options.append([idx, arch])
-        header = ['Option', 'System architecture']
-        printing.print_in_table(options, headers=header)
-        idx_list = range(len(system_architectures))
-        system_architecture = prompting.user_prompt_select_options(*idx_list, allow_empty=False,
-                                                                   allow_multi=False, message=selmsg)
-        system_architecture = system_architectures[int(system_architecture[0])]
-
-        # Choose which partitions must be run in sequence # TODO: beter tuples ipv lists?
-        # TODO: add checks whether a valid input is given
-        selmsg = 'Please select which partitions must be run in sequence (e.g. [[1,2], [3,4]])'
-        sequence_partitions = prompting.user_prompt_string(allow_empty=True, message=selmsg)
-        if sequence_partitions:
-            sequence_partitions = eval(sequence_partitions)
-        else:
-            sequence_partitions = []
+        printing.print_in_table(options, headers=['Option', 'System architecture'])
+        msg = 'Please select system architecture:'
+        system_architecture = str(prompting.user_prompt_select_options(*system_architectures, allow_empty=False,
+                                                                       allow_multi=False, message=msg)[0])
+
+        # Select which partitions need a local converger
+        msg = 'Please select which partitions need a local converger:'
+        idx_list = range(len(partitions)+1)
+        while True:
+            local_convergers = prompting.user_prompt_select_options(*idx_list, allow_empty=True, allow_multi=True,
+                                                                    message=msg)
+            # Partitions start counting at one, so zero is not allowed
+            if 0 in local_convergers:
+                print 'Partition numbering starts from 1.'
+                continue
+            # Check if feedback exists in the chosen partitions
+            valid_input = True
+            for converger in local_convergers:
+                feedback = False
+                for idx, node in enumerate(partitions[converger]):
+                    if set(partitions[converger][idx+1:]).intersection(coupling_dict[node]):
+                        feedback = True
+                if not feedback:
+                    print 'Partition {} does not contain feedback and therefore it cannot have a local ' \
+                          'converger'.format(converger)
+                    valid_input = False
+                    continue
+            if not valid_input:
+                continue
+            break
+
+        # Select which partitions needs to be solved using the Jacobi method instead of Gauss-Seidel
+        msg = 'Please select which partitions must be solved using a Jacobi convergence instead of Gauss-Seidel'
+        idx_list = range(len(partitions) + 1)
+        while True:
+            jacobi_convergence = prompting.user_prompt_select_options(*idx_list, allow_empty=True, allow_multi=True,
+                                                                      message=msg)
+            # Partitions start counting at one, so zero is not allowed
+            if 0 in jacobi_convergence:
+                print 'Partition numbering starts from 1.'
+                continue
+            break
+
+        # Select which partitions must be executed in sequence # TODO: add more checks
+        msg = 'Please select which partitions must be run in sequence (e.g. [[1, 2], [3, 4]])'
+        while True:
+            valid_input = True
+            sequence_partitions = prompting.user_prompt_string(allow_empty=True, message=msg)
+            sequence_partitions = eval(sequence_partitions) if sequence_partitions else []
+            sequence_partitions = [sequence_partitions] if not any(isinstance(el, list) for el in sequence_partitions) \
+                else sequence_partitions
+            if not isinstance(sequence_partitions, list) or not any(isinstance(el, list) for el in sequence_partitions):
+                print 'Input should be a list or list of lists'
+                valid_input = False
+            print 'second alteration', sequence_partitions
+            for sequence in sequence_partitions:
+                if 0 in sequence:
+                    print 'Partition numbering starts from 1.'
+                    valid_input = False
+                for element in sequence:
+                    if not isinstance(element, int):
+                        print 'Invalid input given'
+                        valid_input = False
+            if not valid_input:
+                continue
+            break
 
         print 'local converger:', local_convergers
         print 'system architecture:', system_architecture
-        print 'jacobi convergence', jacobi_convergence
-        print 'sequence partitions', sequence_partitions
+        print 'jacobi convergence:', jacobi_convergence
+        print 'sequence partitions:', sequence_partitions
 
         self.graph['problem_formulation']['system_architecture'] = system_architecture
         self.graph['problem_formulation']['local_convergers'] = local_convergers
diff --git a/kadmos/graph/graph_kadmos.py b/kadmos/graph/graph_kadmos.py
index 8c85c5ff50f2990ee0105fcf2bc91a0e36bbb7eb..ad1d3743d203631209e30d50ef4514f170839ebd 100644
--- a/kadmos/graph/graph_kadmos.py
+++ b/kadmos/graph/graph_kadmos.py
@@ -2375,91 +2375,6 @@ class KadmosGraph(nx.DiGraph, EquationMixin, VistomsMixin):
         return {'dict of dicts': nx.convert.to_dict_of_dicts(self, edge_data=1),
                 'SciPy sparse matrix': nx.adjacency_matrix(self)}
 
-    def partition_graph(self, n_parts, edge_weight='coupling_strength', node_weights=None, tpwgts=None, recursive=False,
-                        contig=False):
-        """Partition a graph using the Metis algorithm (http://glaros.dtc.umn.edu/gkhome/metis/metis/overview).
-
-        Note that partitioning can only be performed on undirected graphs. Therefore every graph input is translated
-        into an undirected graph.
-
-        :param n_parts: number of partitions requested (algorithm might provide less)
-        :type n_parts: int
-        :param edge_weight: edge attribute that is used for the edge weight in the partitioning algorithm, only one edge
-        weight is allowed
-        :type edge_weight: str
-        :param node_weights: node attributes that are used for the node weights in the partitioning algorithm, multiple
-        node weights are allowed
-        :type node_weights: list
-        :param tpwgts: list of target partition weights
-        :type tpwgts: list
-        :param recursive: Metis option
-        :type recursive: bool
-        :param contig: Metis option
-        :type contig: bool
-        """
-
-        # TODO: check metis settings
-
-        import metis
-
-        # Get function graph
-        function_graph = self.get_function_graph()
-
-        # Check input
-        if edge_weight:
-            for edge in function_graph.edges.data():
-                assert edge_weight in edge[2], 'Attribute {0} not defined for edge {1}'.format(edge_weight, edge)
-                assert isinstance(edge[2][edge_weight], (int, long)), 'Edge weights must be integers'
-        if node_weights:
-            assert isinstance(node_weights, list)
-            for node_weight in node_weights:
-                for node, data in function_graph.nodes(data=True):
-                    assert node_weight in data, 'Attribute {0} not defined for node {1}'.format(node_weight, node)
-                    assert isinstance(data[node_weight], (int, long)), 'Node weights must be integers'
-
-        # Find nodes that are directly coupled with both feedforward and feedback
-        if edge_weight:
-            remove_edges = []
-            for edge in function_graph.edges(data=True):
-                if (edge[0], edge[1]) in remove_edges:
-                    continue
-                elif (edge[1], edge[0]) in function_graph.edges():
-                    # Combine edge weights of both edges
-                    function_graph.edges[edge[0], edge[1]][edge_weight] += function_graph.edges[edge[1], edge[0]][
-                        edge_weight]
-                    # Remove second edge
-                    remove_edges.append((edge[1], edge[0]))
-            # Remove edges
-            for edge in remove_edges:
-                function_graph.remove_edge(edge[0], edge[1])
-
-        # Get undirected graph
-        g_und = function_graph.to_undirected()
-
-        # Add node and edge weight attributes to undirected graph
-        if edge_weight:
-            g_und.graph['edge_weight_attr'] = edge_weight
-        if node_weights:
-            g_und.graph['node_weight_attr'] = node_weights
-
-        # Get partitions
-        (edgecuts, parts) = metis.part_graph(g_und, n_parts, tpwgts=tpwgts, recursive=recursive, contig=contig)
-
-        # Add partition information to graph
-        i = 0
-        for node, data in g_und.nodes(data=True):
-            self.nodes[node]['part_id'] = parts[i]+1
-            i += 1
-
-        # Create a partition dictionary and add to graph
-        partition_dict = dict()
-        for partition in range(n_parts):
-            nodes_in_partition = [node for idx, node in enumerate(g_und.nodes) if parts[idx] == partition]
-            partition_dict[partition+1] = nodes_in_partition
-        self.graph['problem_formulation']['partitions'] = partition_dict
-
-        return
-
     def add_objective_function_by_nodes(self, *args, **kwargs):
         """This function adds objective functions to the graph using lists of variable nodes.
 
diff --git a/kadmos/graph/graph_process.py b/kadmos/graph/graph_process.py
index 82cf4fe1b4e7d2c911aa30d3edcd68f8b4eac3f7..fbaa8a350e02c6bacd4720c216ed1c5b37587bfc 100644
--- a/kadmos/graph/graph_process.py
+++ b/kadmos/graph/graph_process.py
@@ -395,8 +395,9 @@ class MdaoProcessGraph(ProcessGraph):
             # Check if the partition is the first of a sequence
             partition_sequence = [part]
             for sequence in sequence_partitions:
-                if part == sequence[0]:
-                    partition_sequence = list(sequence)
+                if sequence:
+                    if part == sequence[0]:
+                        partition_sequence = list(sequence)
 
             # Connect the sequence of partitions
             for partition in partition_sequence: