diff --git a/README.md b/README.md
index c38d3d4..14089bc 100644
--- a/README.md
+++ b/README.md
@@ -26,9 +26,10 @@ git clone https://github.com/schoonhovenrichard/BlooPy.git
 ```
 
 ### Dependencies
-**BlooPy** install the following dependencies:
+**BlooPy** installs the following dependencies:
 - bitarray
 - pyswarms
+- networkx
 
 ## Implemented algorithms
 ### Discrete local search algorithms
diff --git a/bloopy/analysis/FFG.py b/bloopy/analysis/FFG.py
new file mode 100644
index 0000000..044b7e9
--- /dev/null
+++ b/bloopy/analysis/FFG.py
@@ -0,0 +1,60 @@
+import networkx as nx
+import itertools
+
+import bloopy.analysis.analysis_utils as anutil
+
+def build_FFG(nidxs_dict, bound_list, method='bounded'):
+    if method not in ["circular", "bounded", "Hamming"]:
+        raise Exception("Unknown neighbour generator given!")
+    var_ranges = anutil.get_variable_ranges(bound_list)
+
+    G = nx.DiGraph()
+    for key, value in nidxs_dict.items():
+        nidx = value[0]
+        if nidx not in G:
+            G.add_node(nidx)
+
+    it = 0
+    for x in itertools.product(*var_ranges):
+        it += 1
+        print('Adding point to CPG, node {}...\r'.format(it), end="")
+        xidx, xfit = nidxs_dict[x]
+        if method == "circular":
+            nbours = anutil.generate_circular_neighbours(x, var_ranges)
+        elif method == "Hamming":
+            nbours = anutil.generate_Hamming_neighbours(x, var_ranges)
+        else:
+            nbours = anutil.generate_bounded_neighbours(x, var_ranges)
+        for nbour in nbours:
+             nidx, nfit = nidxs_dict[nbour]
+             if nfit <= xfit:
+             #if nfit < xfit:
+                 G.add_edge(xidx, nidx)
+    return G
+
+def average_centrality_nodes(centrality_dict, nodelist, space_dict, nidxs_pnts):
+    nodelist_centrality = 0.0
+    for node in nodelist:
+        nodelist_centrality += abs(centrality_dict[node])
+    minima_centralities = []
+    for idx, centr in centrality_dict.items():
+        pnt = nidxs_pnts[idx]
+        #[point, fitness, ptype]
+        if pnt[2] != 1:
+            continue
+        minima_centralities.append(abs(centr))
+    total_centrality = sum([abs(n) for i, n in centrality_dict.items()])
+    minima_centrality = sum(minima_centralities)
+    return nodelist_centrality, total_centrality, minima_centrality
+
+def average_centrality_crits(centrality_dict, nodelist, space_dict, nidxs_pnts):
+    nodelist_centrality = 0.0
+    for node in nodelist:
+        nodelist_centrality += abs(centrality_dict[node])
+    critical_centralities = []
+    for idx, centr in centrality_dict.items():
+        #[point, fitness, ptype]
+        pnt = nidxs_pnts[idx]
+        critical_centralities.append(abs(centr))
+    total_centrality = sum([abs(n) for i, n in centrality_dict.items()])
+
diff --git a/bloopy/analysis/__init__.py b/bloopy/analysis/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/bloopy/analysis/analysis_utils.py b/bloopy/analysis/analysis_utils.py
new file mode 100644
index 0000000..8268c64
--- /dev/null
+++ b/bloopy/analysis/analysis_utils.py
@@ -0,0 +1,139 @@
+import itertools
+
+#from bloopy.individual import individual, continuous individual
+from bloopy.individual import individual
+import bloopy.utils as utils
+
+def get_variable_ranges(bound_list):
+    var_ranges = []
+    for var in range(len(bound_list)):
+        rang = []
+        for i in range(bound_list[var][0], bound_list[var][1]+1):
+            rang.append(i)
+        var_ranges.append(rang)
+    return var_ranges
+
+def generate_Hamming_neighbours(x0, varranges):
+    neighbours = []
+    for var in range(len(x0)):
+        left = x0[:var]
+        right = x0[var+1:]
+        for y in varranges[var]:
+            if y == x0[var]:
+                continue
+            n = left + (y,) + right
+            neighbours.append(n)
+    return neighbours
+
+def generate_circular_neighbours(x0, varranges):
+    neighbours = []
+    for var in range(len(x0)):
+        left = x0[:var]
+        right = x0[var+1:]
+        if x0[var] - 1 >= varranges[var][0]:
+            y = x0[var] - 1
+        else:
+            y = varranges[var][-1]
+        n1 = left + (y,) + right
+        if x0[var] + 1 <= varranges[var][-1]:
+            y = x0[var] + 1
+        else:
+            y = varranges[var][0]
+        n2 = left + (y,) + right
+        neighbours.append(n1)
+        neighbours.append(n2)
+    return neighbours
+
+def generate_bounded_neighbours(x0, varranges):
+    neighbours = []
+    for var in range(len(x0)):
+        n1 = None
+        n2 = None
+        left = x0[:var]
+        right = x0[var+1:]
+        if x0[var] - 1 >= varranges[var][0]:
+            y = x0[var] - 1
+            n1 = left + (y,) + right
+        if x0[var] + 1 <= varranges[var][-1]:
+            y = x0[var] + 1
+            n2 = left + (y,) + right
+        if n1 is not None:
+            neighbours.append(n1)
+        if n2 is not None:
+            neighbours.append(n2)
+    return neighbours
+
+def generate_circular_neighbours_perdim(x0, varranges):
+    neighbours = []
+    for var in range(len(x0)):
+        left = x0[:var]
+        right = x0[var+1:]
+        if x0[var] - 1 >= varranges[var][0]:
+            y = x0[var] - 1
+        else:
+            y = varranges[var][-1]
+        n1 = left + (y,) + right
+        if x0[var] + 1 <= varranges[var][-1]:
+            y = x0[var] + 1
+        else:
+            y = varranges[var][0]
+        n2 = left + (y,) + right
+        neighbours.append([n1, n2])
+    return neighbours
+
+def generate_bounded_neighbours_perdim(x0, varranges):
+    neighbours = []
+    for var in range(len(x0)):
+        n1 = None
+        n2 = None
+        left = x0[:var]
+        right = x0[var+1:]
+        if x0[var] - 1 >= varranges[var][0]:
+            y = x0[var] - 1
+            n1 = left + (y,) + right
+        if x0[var] + 1 <= varranges[var][-1]:
+            y = x0[var] + 1
+            n2 = left + (y,) + right
+        lst = []
+        if n1 is not None:
+            lst.append(n1)
+        if n2 is not None:
+            lst.append(n2)
+        neighbours.append(lst)
+    return neighbours
+
+def generate_Hamming_neighbours_perdim(x0, varranges):
+    neighbours = []
+    for var in range(len(x0)):
+        left = x0[:var]
+        right = x0[var+1:]
+        lst = []
+        for y in varranges[var]:
+            if y == x0[var]:
+                continue
+            n = left + (y,) + right
+            lst.append(n)
+        neighbours.append(lst)
+    return neighbours
+
+def build_nodeidxs_dict(bound_list, fitfunc, bsize):
+    #Give all point node indices for the graph, so an entry
+    # is [point_type, fitness, node_index]
+    indiv = individual(bsize, boundary_list=bound_list)
+    var_ranges = get_variable_ranges(bound_list)
+    node_dict = dict()
+    node_idx = 0
+    for x in itertools.product(*var_ranges):
+        utils.set_bitstring(indiv, list(x))
+        xfit = fitfunc(indiv.bitstring)
+        node_dict[tuple(x)] = [node_idx, xfit]
+        node_idx += 1
+    return node_dict
+
+def indices_to_points(node_dict):
+    # We need a mapping from node-indices to points
+    idxs_to_point = dict()
+    for point, value in node_dict.items():
+        # Save [point, fitness, ptype]
+        idxs_to_point[value[2]] = [point, value[1], value[0]]
+    return idxs_to_point
diff --git a/bloopy/analysis/critical_points.py b/bloopy/analysis/critical_points.py
new file mode 100644
index 0000000..57b606e
--- /dev/null
+++ b/bloopy/analysis/critical_points.py
@@ -0,0 +1,171 @@
+import itertools
+import warnings
+
+import bloopy.analysis.analysis_utils as anutil
+
+
+def classify_points(bsize, bound_list, nidxs_dict, method='bounded'):
+    if method not in ["circular", "bounded", "Hamming"]:
+        raise Exception("Unknown neighbour generator given!")
+    var_ranges = anutil.get_variable_ranges(bound_list)
+
+    count_saddles = 0
+    count_minima = 0
+    count_maxima = 0
+    regular_points = 0
+    total = 0
+
+    #A point is classified as type 0,1,2,3
+    # 0 is regular, 1 is local minimum,
+    # 2 is local maximum, and 3 is saddlepoint
+    # we save points as tuple (type, fit, node_index) for convenience
+    point_class = dict()
+    for x in itertools.product(*var_ranges):
+        nr_dims_is_minimal = 0
+        isRegular = False
+        xidx, xfit = nidxs_dict[x]
+        if method == "circular":
+            nbours = anutil.generate_circular_neighbours_perdim(x, var_ranges)
+        elif method == "Hamming":
+            nbours = anutil.generate_Hamming_neighbours_perdim(x, var_ranges)
+        else:
+            nbours = anutil.generate_bounded_neighbours_perdim(x, var_ranges)
+        for dim in range(len(nbours)):
+            nbours_larger = 0
+            if len(nbours[dim]) == 0:
+                raise Exception("Should already have been removed when processing space.")
+            elif len(nbours[dim]) == 1:# Only one neighbour
+                n1 = nbours[dim][0]
+                n1_idxs, n1fit = nidxs_dict[n1]
+                if n1fit == xfit:#This is for the failfit
+                    if xfit > 10000:
+                        isRegular = True
+                        break
+                    else:
+                        print(xfit, "1")
+                        warnings.warn("Equal fitness found")
+                        #raise Exception("Equal fitness found")
+                if n1fit >= xfit:
+                    nbours_larger += 1
+                    nr_dims_is_minimal += 1
+            elif len(nbours[dim]) == 2:# Only one neighbour
+                n1, n2 = nbours[dim]
+                n1_idxs, n1fit = nidxs_dict[n1]
+                n2_idxs, n2fit = nidxs_dict[n2]
+                if n1fit == xfit or n2fit == xfit:#This is for the failfit
+                    if xfit > 10000:
+                        isRegular = True
+                        break
+                    else:
+                        print(xfit, "2")
+                        warnings.warn("Equal fitness found")
+                        #raise Exception("Equal fitness found")
+                if n1fit >= xfit:
+                    nbours_larger += 1
+                if n2fit >= xfit:
+                    nbours_larger += 1
+                if nbours_larger == 1:#Its not a critical point as it has
+                    # a non-zero gradient in this dimension
+                    isRegular = True
+                    break
+                elif nbours_larger == 2:
+                    nr_dims_is_minimal += 1
+            else:# More neighbours
+                nbours_larger = 0
+                for ni in nbours[dim]:
+                    ni_idxs, nifit = nidxs_dict[ni]
+                    if nifit == xfit:#This is for the failfit
+                        if xfit > 10000:
+                            isRegular = True
+                            break
+                        else:
+                            warnings.warn("Equal fitness found")
+                            print(xfit, "N")
+                            #raise Exception("Equal fitness found")
+                    if nifit >= xfit:
+                        nbours_larger += 1
+                if nbours_larger == len(nbours[dim]):#minimal in this dim
+                    nr_dims_is_minimal += 1
+                elif nbours_larger == 0:#maximal in this dim
+                    pass
+                else:#Its not a critical point as it has
+                    # a non-zero gradient in this dimension
+                    isRegular = True
+                    break
+
+        total += 1
+        if x in point_class.keys():
+            raise Exception("SOMETHING WRONG")
+
+        if isRegular:
+            regular_points += 1
+            point_class[x] = [0, xfit, xidx]
+            continue
+
+        if nr_dims_is_minimal == 0:
+            count_maxima += 1
+            point_class[x] = [2, xfit, xidx]
+        elif nr_dims_is_minimal == len(bound_list):
+            count_minima += 1
+            point_class[x] = [1, xfit, xidx]
+        else:
+            count_saddles += 1
+            point_class[x] = [3, xfit, xidx]
+    return total, count_minima, count_maxima, count_saddles, regular_points, point_class
+
+def strong_local_minima(p, globfit, space_dict):
+    loc_min_idxs = []
+    for pt, value in space_dict.items():
+        ptype, pfit, pidx = value
+        if ptype != 1:
+            continue
+        if pfit <= (1+p)*globfit:
+            loc_min_idxs.append(pidx)
+    return loc_min_idxs
+
+def strong_critical_points(p, globfit, space_dict):
+    loc_min_idxs = []
+    for pt, value in space_dict.items():
+        ptype, pfit, pidx = value
+        if ptype == 0:
+            continue
+        if pfit <= (1+p)*globfit:
+            loc_min_idxs.append(pidx)
+    return loc_min_idxs
+
+def sizes_minima(sspace, bsize, bound_list, fitfunc, method="bounded"):
+    if method not in ["circular", "bounded", "Hamming"]:
+        raise Exception("Unknown neighbour generator given!")
+    indiv = individual(bsize, boundary_list=bound_list)
+    var_ranges = anutil.get_variable_ranges(bound_list)
+    count_minima = 0
+    total = 0
+    hole_depths = []
+    for x in itertools.product(*var_ranges):
+        anutil.set_bitstring(indiv, list(x))
+        xfit = fitfunc(indiv.bitstring)
+
+        if method == "circular":
+            nbours = anutil.generate_circular_neighbours(x, var_ranges)
+        elif method == "Hamming":
+            nbours = anutil.generate_Hamming_neighbours(x, var_ranges)
+        else:
+            nbours = anutil.generate_bounded_neighbours(x, var_ranges)
+
+        if xfit > 10000:
+            is_minimal = False
+        else:
+            is_minimal = True
+            nfits = []
+            for nbour in nbours:
+                anutil.set_bitstring(indiv, nbour)
+                nfit = fitfunc(indiv.bitstring)
+                nfits.append(nfit)
+                if nfit < xfit:# Neighbour has lower fitness
+                    is_minimal = False
+            if is_minimal:
+                hole_depth = min(nfits) - xfit
+                hole_depths.append([xfit, hole_depth])
+                count_minima += 1
+        total += 1
+    return count_minima, total, np.array(hole_depths)
diff --git a/bloopy/utils.py b/bloopy/utils.py
index 24f75c3..f5f3c47 100644
--- a/bloopy/utils.py
+++ b/bloopy/utils.py
@@ -127,6 +127,16 @@ def calculate_bitstring_length(sspace):
         bsize += len(sspace[key])
     return bsize
 
+def set_bitstring(sol, arr):
+    if sol.boundary_list is None:
+        raise Exception("Method not implemented for standard bitstrings")
+    if len(arr) != len(sol.boundary_list):
+        raise Exception("Invalid settings for bitstring given!")
+    indices = [i for i, x in enumerate(list(sol.bitstring)) if x]
+    for k in range(len(arr)):
+        sol.bitstring[indices[k]] = 0 # Set old one to 0
+        sol.bitstring[arr[k]] = 1 # set new one to 1
+
 def generate_population(population_size, ffunc, sspace):
     splits = generate_boundary_list(sspace)
     input_pop = []
diff --git a/setup.py b/setup.py
index be5318b..3743bb5 100644
--- a/setup.py
+++ b/setup.py
@@ -3,7 +3,7 @@
 setup(
     name = 'bloopy',
     packages=find_packages(),
-    version = '0.2',  # Ideally should be same as your GitHub release tag varsion
+    version = '0.3',  # Ideally should be same as your GitHub release tag varsion
     description = 'BlooPy: Black-box optimization Python for bitstring, categorical, and numerical discrete problems with local, and population-based algorithms.',
     author = 'Richard Schoonhoven',
     author_email = 'r.a.schoonhoven@hotmail.com',
@@ -15,5 +15,6 @@
         'numpy>=1.19.0',
         'scipy>=1.6.0',
         'bitarray',
-        'pyswarms'],
+        'pyswarms',
+        'networkx'],
 )