diff --git a/problems/n_dim/square.py b/problems/n_dim/square.py
index fad24577..022b5a27 100644
--- a/problems/n_dim/square.py
+++ b/problems/n_dim/square.py
@@ -16,7 +16,7 @@ def __init__(self, n):
         super().__init__()
         self.x_min = 1e-3 * np.ones(n)  # cuz a subproblem uses both, whereas a problem only has x_min
         self.x_max = np.ones(n)  # cuz a subproblem uses both, whereas a problem only has x_max
-        self.x0 = np.linspace(0.8, 0.9, n)
+        self.x0 = np.linspace(1, 2, n)
         self.n = n
         self.m = 1
         self.f = np.zeros(n)
diff --git a/sao/intervening_variables/mma.py b/sao/intervening_variables/mma.py
index 69b916ae..c7eb1e5b 100644
--- a/sao/intervening_variables/mma.py
+++ b/sao/intervening_variables/mma.py
@@ -1,5 +1,6 @@
 import numpy as np
 
+from .intervening import Intervening
 from .exponential import Exponential
 from .split import PositiveNegative
 
@@ -195,3 +196,75 @@ def get_asymptotes(self, x):
             low = x - self.dist * self.dx
             upp = x + self.dist * self.dx
         return low, upp
+
+
+class MMA07(Intervening):
+    """
+    MMA implementation as proposed in the note "MMA and GCMMA - two methods for nonlinear optimization"
+    Krister Svanberg (KTH Stockholm), 2007
+
+    Some practical considerations:
+    1. Scale the constraint such that 1 < fj(x) = gj(x) - gjmax < 100 and 1 < f0(x) < 100
+    2. Scale variables such that 0.1 < ximax - ximin < 100
+    3. When working with artificial variables, start with cj = 100, and increase by factors of 10 until all yj are zero
+    """
+
+    def __init__(self, x_min=0.0, x_max=1.0, sinit=0.5, sincr=1.2, sdecr=0.7, asybound=10.0, oscillation_tol=1e-10,
+                 p=-1, factor=0.01):
+        super().__init__(p=p, factor=factor)
+        self.x, self.xold1, self.xold2 = None, None, None
+        self.dx = x_max - x_min
+
+        self.asybound = asybound
+        self.sinit = sinit
+        self.sincr = sincr
+        self.sdecr = sdecr
+        self.oscillation_tol = oscillation_tol
+
+        self.dist = None
+        self.dist_min, self.dist_max = 1 / (self.asybound ** 2), self.asybound
+
+    def y(self, x):
+        return
+
+    def dydx(self, x):
+        return np.where(self.positive, self.right.dydx(x), self.left.dydx(x))
+
+    def ddyddx(self, x):
+        return np.where(self.positive, self.right.ddyddx(x), self.left.ddyddx(x))
+
+    def clip(self, x):
+        self.left.clip(x)
+        self.right.clip(x)
+        return x
+
+
+def get_asymptotes(self, x):
+    self.xold2, self.xold1, self.x = self.xold1, self.x, x.copy()
+    """Increases or decreases the asymptotes interval based on oscillations in the design vector"""
+    if self.dist is None:
+        self.dist = np.full_like(self.x, self.sinit)
+
+    if self.xold2 is None:
+        # Initial values of asymptotes
+        low = x - self.dist * self.dx
+        upp = x + self.dist * self.dx
+    else:
+        # Update scheme for asymptotes
+        # depending on if the signs of (x_k-xold) and (xold-xold2) are opposite, indicating an oscillation in xi
+        # if the signs are equal the asymptotes are slowing down the convergence and should be relaxed
+
+        # check for oscillations in variables (if > 0: no oscillations, if < 0: oscillations)
+        oscillation = ((x - self.xold1) * (self.xold1 - self.xold2)) / self.dx
+
+        # oscillating variables x_i are increase or decrease the factor
+        self.dist[oscillation > +self.oscillation_tol] *= self.sincr
+        self.dist[oscillation < -self.oscillation_tol] *= self.sdecr
+
+        # Clip the asymptote factor
+        np.clip(self.dist, self.dist_min, self.dist_max)
+
+        # update lower and upper asymptotes
+        low = x - self.dist * self.dx
+        upp = x + self.dist * self.dx
+    return low, upp
diff --git a/sao/mappings/approximations.py b/sao/mappings/approximations.py
new file mode 100644
index 00000000..c1cb9a0b
--- /dev/null
+++ b/sao/mappings/approximations.py
@@ -0,0 +1,120 @@
+import numpy as np
+from .mapping import Mapping, Linear
+
+
+class LinearApproximation(Mapping):
+    """
+    Linear Approximation (LA) f[x] of function g[x] at x0.
+
+    f[x]    = g[x0] + g'[x0]*(x - x0)
+            = (g[x0] - g'[x0]*x0) + g'[x0]*x
+
+    f'[x]   = g'[x0]
+
+    f''[x]  = 0
+    """
+
+    def __init__(self, mapping=Linear()):
+        """
+        Initialization of LA.
+
+        :param mapping: The dependent mapping.
+
+        If mapping is not provided, dependent mapping is Linear() which ends the chain.
+        """
+        super().__init__(mapping)  # Initialization of dependent mapping.
+        self.g0, self.dg0 = None, None
+
+    def _update(self, x0, dg0, ddg0=0):
+        self.g0 = -dg0 * x0
+        self.dg0 = dg0
+        return self._g(x0), self._dg(dg0), self._ddg(ddg0)  # Why not return self.g0, self.dg0, self.ddg0?
+
+    def _g(self, x):
+        """
+        Function value of LA function at x.
+
+        f[x]    = g[x0] + g'[x0]*(x - x0)
+                = (g[x0] - g'[x0]*x0) + g'[x0]*x
+                = g[x0] + a + b*x, with
+                    a = -g'[x0]*x0
+                    b = g'[x0]
+
+        :param x: Incoming variable (or function) value
+        :return: Function value at x
+        """
+        return self.g0 + self.dg0 * x  # Excludes g[x0] (on purpose)
+
+    def _dg(self, x): return self.dg0
+
+    def _ddg(self, x): return np.zeros_like(x)
+
+
+class DiagonalQuadraticApproximation(LinearApproximation):
+    """
+    Diagonal Quadratic Approximation (DQA) f[x] of function g[x] at x0.
+
+    DQA builds on top of LA.
+    Explanation of what meaning Diagonal Quadratic
+
+    f[x]    = g[x0] + g'[x0]*(x - x0) + 1/2*g''[x0]*(x - x0)^2
+            = (g[x0] - g'[x0]*x0 + 1/2*g''[x0]*x0^2) + (g'[x0] - g''[x0]*x0)*x + 1/2*g''[x0]*x^2
+            = LA[x] + 1/2*g''[x0]*(x0^2 - 2*x0*x + x^2)
+            = (LA[x] + 1/2*g''[x0]*x0^2) - g''[x0]*x0*x + 1/2*g''[x0]*x^2
+
+    f'[x]   = (LA'[x] - g''[x0]*x0) + g''[x0]*x
+
+    f''[x]  = LA''[x] + g''[x0], with LA''[x] = 0
+    """
+
+    def __init__(self, mapping=Linear()):
+        super().__init__(mapping)
+        self.ddg0 = None
+
+    def _update(self, x0, dg0, ddg0=0):
+
+        super()._update(x0, dg0)
+        self.g0 += 0.5 * ddg0 * x0 ** 2
+        self.dg0 -= ddg0 * x0
+        self.ddg0 = ddg0
+
+    def _g(self, x):
+        """
+        Function value of DQA function at x.
+
+        f[x]    = g[x0] + g'[x0]*(x - x0) + 1/2*g''[x0]*(x - x0)^2
+                = (g[x0] - g'[x0]*x0 + 1/2*g''[x0]*x0^2) + (g'[x0] - g''[x0]*x0)*x + 1/2*g''[x0]*x^2
+                = g[x0] + a + b*x + c*x^2, with
+                    a = -g'[x0]*x0 + 1/2*g''[x0]*x0^2
+                    b = g'[x0] - g''[x0]*x0
+                    c = 1/2*g''[x0]
+
+        :param x: Incoming variable (or function) value
+        :return: Function value at x
+        """
+
+        return self.g0 + self.dg0 * x + 0.5 * self.ddg0 * x ** 2
+
+    def _dg(self, x):
+        """
+        First derivative of DQA function at x.
+
+        f'[x] = (g'[x0] - g''[x0]*x0) + g''[x0]*x
+
+        :param x: Incoming variable (or function) value
+        :return: First derivative at x
+        """
+
+        return self.dg0 + self.ddg0 * x
+
+    def _ddg(self, x):
+        """
+        Second derivative of DQA function at x.
+
+        f''[x] = g''[x0]
+
+        :param x: Incoming variable (or function) value
+        :return: Second derivative at x
+        """
+
+        return self.ddg0
diff --git a/sao/mappings/change_of_variable.py b/sao/mappings/change_of_variable.py
new file mode 100644
index 00000000..347d9c23
--- /dev/null
+++ b/sao/mappings/change_of_variable.py
@@ -0,0 +1,52 @@
+from abc import abstractmethod, ABC
+from .mapping import Mapping, Linear, Conditional
+import numpy as np
+
+
+class Exponential(Mapping):
+    def __init__(self, p=1, mapping=Linear(), xlim=1e-10):
+        super().__init__(mapping)
+        assert p != 0, f"Invalid power x^{p}, will result in zero division."
+        self.p, self.xlim = p, xlim
+
+    def _clip(self, x): return np.maximum(x, self.xlim, out=x) if self.p < 0 else x
+
+    def _g(self, x): return x ** self.p
+
+    def _dg(self, x): return self.p * x ** (self.p - 1)
+
+    def _ddg(self, x): return self.p * (self.p - 1) * x ** (self.p - 2)
+
+
+class MMAp(Conditional, ABC):
+    def __init__(self, p=-1, factor=1e-3, low=-10.0, upp=10.0):
+        super().__init__(Exponential(p), Exponential(p))
+        self.low, self.upp = low, upp
+        self.factor = factor
+
+    def update(self, x0, dg0, ddg0=0):
+        super().update(x0, dg0, ddg0)
+        [self.low, self.upp] = self.get_asymptotes(x0)
+
+    @abstractmethod
+    def get_asymptotes(self, x): pass
+
+    def _g(self, x):
+        return super().g(np.where(self.condition, self.upp - x, x - self.low))
+
+    def _dg(self, x):
+        return super().dg(np.where(self.condition, self.upp - x, x - self.low)) * np.where(self.positive, -1, +1)
+
+    def _ddg(self, x):
+        return super().ddg(np.where(self.condition, self.upp - x, x - self.low))
+
+    def clip(self, x):
+        return np.clip(x, self.low + self.factor, self.upp - self.factor, out=x)
+
+
+class ConLin(Conditional):
+    def __init__(self): super().__init__(Exponential(-1), Exponential(1))
+
+
+class MMA(MMAp, ABC):
+    def __init__(self): super().__init__(-1)
diff --git a/sao/mappings/mapping.py b/sao/mappings/mapping.py
new file mode 100644
index 00000000..50a57f9e
--- /dev/null
+++ b/sao/mappings/mapping.py
@@ -0,0 +1,169 @@
+from abc import ABC
+import numpy as np
+
+
+class Mapping(ABC):
+    def __init__(self, mapping=None):
+        self.child = mapping if mapping is not None else Linear()
+
+    def update(self, x0, dg0, ddg0=0):
+        if self.child is not None:
+            x0, dg0, ddg0 = self.child.update(x0, dg0, ddg0)
+        self._update(x0, dg0, ddg0)
+        return self._g(x0), self._dg(dg0), self._ddg(ddg0)
+
+    @property
+    def name(self):
+        """
+        Set Mapping name.
+
+        :return: Mapping name
+        """
+        return self.__class__.name
+
+    def clip(self, x):
+        return self._clip(self.child.clip(x))
+
+    def g(self, x):
+        """
+        Chain rule.
+
+        f[x] = f[g[x]]
+
+        :param x:
+        :return:
+        """
+
+        return self._g(self.child.g(x))
+
+    def dg(self, x):
+        """
+        Chain rule of first derivative.
+
+        f'[x] = f'[g[x]]*g'[x]
+
+        :param x:
+        :return:
+        """
+
+        return self._dg(self.child.g(x)) * self.child.dg(x)
+
+    def ddg(self, x):
+        """
+        Chain rule of second derivative.
+
+        f''[x] = f''[g[x]]*(g'[x])^2 + f'[g[x]]*g''[x]
+
+        :param x:
+        :return:
+        """
+        return self._ddg(self.child.g(x)) * (self.child.dg(x)) ** 2 + \
+               self._dg(self.child.g(x)) * self.child.ddg(x)
+
+    def _update(self, x0, dg0, ddg0=0):
+        pass
+
+    def _clip(self, x):
+        return x
+
+    def _g(self, x):
+        return x
+
+    def _dg(self, x):
+        return np.ones_like(x, dtype=float)
+
+    def _ddg(self, x):
+        return np.zeros_like(x, dtype=float)
+
+
+class Linear(Mapping, ABC):
+    """
+    Linear Mapping; end of chain.
+    """
+
+    def __init__(self):
+        """
+        Initialization of the Linear Mapping.
+
+        Sets child to None.
+        No initialization of child.
+        """
+        self.child = None
+
+    def g(self, x): return x
+
+    def dg(self, x): return np.ones_like(x, dtype=float)
+
+    def ddg(self, x): return np.zeros_like(x, dtype=float)
+
+
+class MixedMapping(Mapping):
+    def __init__(self, n: int, m: int, default: Mapping = Linear()):
+        self.map = [(default, np.arange(0, m), np.arange(0, n))]
+
+    def __setitem__(self, key, inter: Mapping):
+        self.map.append((inter, key[0], key[1]))
+
+    def eval(self, x, fn: callable):
+        out = np.ones((len(self.map[0][1]), x.shape[0]), dtype=float)
+        for mp, responses, variables in self.map:
+            out[np.ix_(responses, variables)] = fn(mp, x[variables])
+        return out
+
+    def g(self, x):
+        return self.eval(x, lambda cls, y: cls.g(y))
+
+    def dg(self, x):
+        return self.eval(x, lambda cls, y: cls.dg(y))
+
+    def ddg(self, x):
+        return self.eval(x, lambda cls, y: cls.ddg(y))
+
+    def update(self, x0, dg0, ddg0=0):
+        for mp, resp, var in self.map:
+            mp.update(x0[var], dg0[np.ix_(resp, var)], ddg0=0)
+        return self
+
+    def clip(self, x):
+        for mp, _, var in self.map:
+            mp.clip(x[var])
+        return x
+
+
+class TwoMap(Mapping, ABC):
+    def __init__(self, left: Mapping, right: Mapping):
+        self.left = left
+        self.right = right
+
+    def update(self, x0, dg0, ddg0=0):
+        self.left.update(x0, dg0, ddg0)
+        self.right.update(x0, dg0, ddg0)
+
+    def clip(self, x):
+        self.left.clip(x)
+        self.right.clip(x)
+        return x
+
+
+class Sum(TwoMap):
+    def g(self, x): return self.left.g(x) + self.right.g(x)
+
+    def dg(self, x): return self.left.dg(x) + self.right.dg(x)
+
+    def ddg(self, x): return self.left.ddg(x) + self.right.ddg(x)
+
+
+class Conditional(TwoMap, ABC):
+    def __init__(self, left: Mapping, right: Mapping):
+        super().__init__(left, right)
+        self.condition = None
+
+    def update(self, x0, dg0, ddg0=0):
+        super().update(x0, dg0, ddg0)
+        self.condition = dg0 >= 0
+
+    def g(self, x): return np.where(self.condition, self.right.g(x), self.left.g(x))
+
+    def dg(self, x): return np.where(self.condition, self.right.dg(x), self.left.dg(x))
+
+    def ddg(self, x): return np.where(self.condition, self.right.ddg(x), self.left.ddg(x))
diff --git a/sao/solvers/primal_dual_interior_point.py b/sao/solvers/primal_dual_interior_point.py
index b58a835a..549e3109 100644
--- a/sao/solvers/primal_dual_interior_point.py
+++ b/sao/solvers/primal_dual_interior_point.py
@@ -1,7 +1,6 @@
 from abc import ABC
 from copy import deepcopy
 from dataclasses import dataclass, fields
-from numba import njit
 import numpy as np
 from scipy.sparse import diags
 
diff --git a/tests/mappings/test_approx.py b/tests/mappings/test_approx.py
new file mode 100644
index 00000000..c9e40594
--- /dev/null
+++ b/tests/mappings/test_approx.py
@@ -0,0 +1,178 @@
+from sao.mappings.approximations import LinearApproximation as LA
+from sao.mappings.approximations import DiagonalQuadraticApproximation as DQA
+from sao.mappings.change_of_variable import Exponential as Exp
+from sao import intervening_variables, approximations
+import numpy as np
+import pytest
+from sao.problems import Problem
+
+
+class Dummy(Problem):
+    def __init__(self, n):
+        super().__init__()
+        self.n = n
+        self.x0 = np.linspace(1.0, 2.0, self.n, dtype=float)
+
+    def g(self, x): return x ** 3
+
+    def dg(self, x): return 3 * x ** 2
+
+    def ddg(self, x): return 6 * x
+
+
+class TestLinearApproximation:
+    tol = 1e-4
+    dx = 1
+
+    def test_la(self):
+        problem = Dummy(10)
+        x0 = problem.x0
+        f = problem.g(x0)
+        df = problem.dg(x0)
+
+        lin_approx = LA()
+        lin_approx.update(x0, f, df)
+
+        # Check validity at x0
+        assert lin_approx.g(x0) + f == pytest.approx(f, self.tol)
+        assert lin_approx.dg(x0) == pytest.approx(df, self.tol)
+        assert lin_approx.ddg(x0) == pytest.approx(0, self.tol)
+
+# def test_ta(dx=1, tol=1e-4):
+#     prob = Dummy(4)
+#     x = prob.x0
+#     f = prob.g(x)
+#     df = prob.dg(x)
+#
+#     old = approximations.Taylor1()
+#     old.update(x, f, df)
+#
+#     new = LA()
+#     new.update(x, df)
+#
+#     assert f + np.sum(new.g(x), 1) == pytest.approx(old.g(x), tol)
+#     assert new.dg(x) == pytest.approx(old.dg(x), tol)
+#     assert new.ddg(x) == pytest.approx(old.ddg(x), tol)
+#
+#     y = x + dx
+#     assert f + np.sum(new.g(y), 1) == pytest.approx(old.g(y), tol)
+#     assert new.dg(y) == pytest.approx(old.dg(y), tol)
+#     assert new.ddg(y) == pytest.approx(old.ddg(y), tol)
+#
+#
+# def test_ta_lin(dx=1, tol=1e-4):
+#     prob = Dummy(4)
+#     x = prob.x0
+#     f = prob.g(x)
+#     df = prob.dg(x)
+#
+#     # Oldskool aka old
+#     old = approximations.Taylor1(intervening_variables.Exponential(1))
+#     old.update(x, f, df)
+#
+#     # Newskool aka new
+#     new = LA(Exp(1))
+#     new.update(x, df)
+#
+#     assert f + np.sum(new.g(x), 1) == pytest.approx(old.g(x), tol)
+#     assert new.dg(x) == pytest.approx(old.dg(x), tol)
+#     assert new.ddg(x) == pytest.approx(old.ddg(x), tol)
+#
+#     y = x + dx
+#     assert f + np.sum(new.g(y), 1) == pytest.approx(old.g(y), tol)
+#     assert new.dg(y) == pytest.approx(old.dg(y), tol)
+#     assert new.ddg(y) == pytest.approx(old.ddg(y), tol)
+#
+#
+# def test_ta_rec(dx=1, tol=1e-4):
+#     prob = Dummy(10)
+#     x = prob.x0
+#     f = prob.g(x)
+#     df = prob.dg(x)
+#
+#     # Oldskool aka old
+#     old = approximations.Taylor1(intervening_variables.Exponential(-1))
+#     old.update(x, f, df)
+#
+#     # Newskool aka new
+#     new = LA(Exp(-1))
+#     new.update(x, df)
+#
+#     assert f + np.sum(new.g(x), 1) == pytest.approx(old.g(x), tol)
+#     assert new.dg(x) == pytest.approx(old.dg(x), tol)
+#     assert new.ddg(x) == pytest.approx(old.ddg(x), tol)
+#
+#     y = x + dx
+#     assert f + np.sum(new.g(y), 1) == pytest.approx(old.g(y), tol)
+#     assert new.dg(y) == pytest.approx(old.dg(y), tol)
+#     assert new.ddg(y) == pytest.approx(old.ddg(y), tol)
+#
+#
+# def test_ta_ta_rec(dx=1, tol=1e-4):
+#     prob = Dummy(10)
+#     x = prob.x0
+#     f = prob.g(x)
+#     df = prob.dg(x)
+#
+#     # Oldskool aka old
+#     old = approximations.Taylor1(intervening_variables.Exponential(2))
+#     old.update(x, f, df)
+#
+#     # Newskool aka new
+#     new = LA(LA(Exp(2)))
+#     new.update(x, df)
+#
+#     assert f + np.sum(new.g(x), 1) == pytest.approx(old.g(x), tol)
+#     assert new.dg(x) == pytest.approx(old.dg(x), tol)
+#     assert new.ddg(x) == pytest.approx(old.ddg(x), tol)
+#
+#     y = x + dx
+#     assert f + np.sum(new.g(y), 1) == pytest.approx(old.g(y), tol)
+#     assert new.dg(y) == pytest.approx(old.dg(y), tol)
+#     assert new.ddg(y) == pytest.approx(old.ddg(y), tol)
+#
+#
+# def test_ta2(dx=1, tol=1e-4):
+#     prob = Dummy(4)
+#     x = prob.x0
+#     f = prob.g(x)
+#     df = prob.dg(x)
+#     ddf = prob.ddg(x)
+#
+#     old = approximations.Taylor2()
+#     old.update(x, f, df, ddf)
+#
+#     new = DQA()
+#     new.update(x, df, ddf)
+#
+#     assert f + np.sum(new.g(x), 1) == pytest.approx(old.g(x), tol)
+#     assert new.dg(x) == pytest.approx(old.dg(x), tol)
+#     assert new.ddg(x) == pytest.approx(old.ddg(x), tol)
+#
+#     y = x + dx
+#     assert f + np.sum(new.g(y), 1) == pytest.approx(old.g(y), tol)
+#     assert new.dg(y) == pytest.approx(old.dg(y), tol)
+#     assert new.ddg(y) == pytest.approx(old.ddg(y), tol)
+#
+#
+# def test_ta2_rec(dx=1, tol=1e-4):
+#     prob = Dummy(4)
+#     x = prob.x0
+#     f = prob.g(x)
+#     df = prob.dg(x)
+#     ddf = prob.ddg(x)
+#
+#     old = approximations.Taylor2(intervening_variables.Exponential(-1))
+#     old.update(x, f, df, ddf)
+#
+#     new = DQA(Exp(-1))
+#     new.update(x, df, ddf)
+#
+#     assert f + np.sum(new.g(x), 1) == pytest.approx(old.g(x), tol)
+#     assert new.dg(x) == pytest.approx(old.dg(x), tol)
+#     assert new.ddg(x) == pytest.approx(old.ddg(x), tol)
+#
+#     y = x + dx
+#     assert f + np.sum(new.g(y), 1) == pytest.approx(old.g(y), tol)
+#     assert new.dg(y) == pytest.approx(old.dg(y), tol)
+#     assert new.ddg(y) == pytest.approx(old.ddg(y), tol)
diff --git a/tests/mappings/test_change_of_variable.py b/tests/mappings/test_change_of_variable.py
new file mode 100644
index 00000000..9ccc1e66
--- /dev/null
+++ b/tests/mappings/test_change_of_variable.py
@@ -0,0 +1,71 @@
+from problems.n_dim.square import Square
+from sao.mappings.change_of_variable import Exponential as Exp
+from sao.mappings.change_of_variable import ConLin
+import numpy as np
+import pytest
+
+
+def mapping_test(mapping, x, g, dg, ddg, tol=1e-4):
+    assert mapping.g(x) == pytest.approx(g, tol)
+    assert mapping.dg(x) == pytest.approx(dg, tol)
+    assert mapping.ddg(x) == pytest.approx(ddg, tol)
+
+
+class TestExp:
+    """
+    Test simple exponential function give correct function values and derivatives.
+    """
+
+    x = np.array([1.0, 2.0, 3.0])  # Note this only tests for simple x; try some more
+
+    # Tests fails for a variable that takes zero value!
+    # Exp(0) should raise error.
+
+    def test_lin(self): mapping_test(Exp(1), self.x, self.x, 1, 0)
+
+    def test_exp2(self): mapping_test(Exp(2), self.x, self.x ** 2, 2 * self.x, 2)
+
+    def test_rec(self): mapping_test(Exp(-1), self.x, 1 / self.x, -1 / self.x ** 2, 2 / self.x ** 3)
+
+    @pytest.mark.parametrize('i', [-3, -2, -1, 1, 2, 3])
+    def test_exp(self, i): mapping_test(Exp(i), self.x, self.x ** i, i * self.x ** (i - 1),
+                                        i * (i - 1) * self.x ** (i - 2))
+
+
+class TestStackedExp:
+    """
+    Test stacked exponential functions give correct function values and derivatives.
+    """
+
+    x = np.array([1.0, 2.0, 3.0])  # Note this only tests for simple x; try some more
+
+    def test_lin_rec(self): mapping_test(Exp(1, Exp(-1)), self.x, 1 / self.x, -1 / self.x ** 2, 2 / self.x ** 3)
+
+    def test_rec_lin(self): mapping_test(Exp(-1, Exp(1)), self.x, 1 / self.x, -1 / self.x ** 2, 2 / self.x ** 3)
+
+    def test_rec_rec(self): mapping_test(Exp(-1, Exp(-1)), self.x, self.x, 1, 0)
+
+    def test_rec_exp2_rec(self):
+        map = Exp(-1, Exp(2, Exp(-1)))
+        map2 = Exp(2)
+        mapping_test(map, self.x, map2.g(self.x), map2.dg(self.x), map2.ddg(self.x))  # Requires some comments
+
+
+class TestTwoMap:
+    def test_conlin(self, tol=1e-4):
+        problem = Square(10)  # Also test different initial values
+        df = problem.dg(problem.x0)
+        conlin = ConLin()
+        conlin.update(problem.x0, df)
+
+        dx = 1
+        y = problem.x0 + dx
+
+        assert conlin.g(y)[0, :] == pytest.approx(y, tol)
+        assert conlin.g(y)[1, :] == pytest.approx(1 / y, tol)
+
+        assert conlin.dg(y)[0, :] == pytest.approx(1, tol)
+        assert conlin.dg(y)[1, :] == pytest.approx(-1 / y ** 2, tol)
+
+        assert conlin.ddg(y)[0, :] == pytest.approx(0, tol)
+        assert conlin.ddg(y)[1, :] == pytest.approx(2 / y ** 3, tol)
diff --git a/tests/mappings/test_dqa.py b/tests/mappings/test_dqa.py
new file mode 100644
index 00000000..8fbe9562
--- /dev/null
+++ b/tests/mappings/test_dqa.py
@@ -0,0 +1,94 @@
+from sao.problems import Problem
+from sao.mappings.approximations import LinearApproximation as LA
+from sao.mappings.approximations import DiagonalQuadraticApproximation as DQA
+from sao.mappings.change_of_variable import Exponential as Exp
+import numpy as np
+import pytest
+
+
+class Dummy(Problem):
+    def __init__(self, n):
+        super().__init__()
+        self.n = n
+        self.x0 = np.linspace(1.0, 2.0, self.n, dtype=float)
+
+    def g(self, x): return x ** 3
+
+    def dg(self, x): return 3 * x ** 2
+
+    def ddg(self, x): return 6 * x
+
+
+def test_aoa_rec(dx=1, tol=1e-4):
+    prob = Dummy(4)
+    x = prob.x0
+    f = prob.g(x)
+    df = prob.dg(x)
+
+    t1_rec = LA(Exp(-1))
+    t1_rec.update(x, df)
+
+    aoa = DQA()
+    aoa.update(x, df, ddg0=t1_rec.ddg(x))
+
+    assert aoa.g(x) == pytest.approx(t1_rec.g(x), tol)
+    assert (f + np.sum(aoa.g(x))) == pytest.approx(f, tol)
+    assert aoa.dg(x) == pytest.approx(t1_rec.dg(x), tol)
+    assert aoa.dg(x) == pytest.approx(df, tol)
+    assert aoa.ddg(x) == pytest.approx(df / t1_rec.map.dg(x) * t1_rec.map.ddg(x), tol)
+
+    y = x + dx
+
+    assert aoa.ddg(y) == pytest.approx(aoa.ddg(x), tol)
+    assert aoa.dg(y) == pytest.approx(df + t1_rec.ddg(x) * (y - x), tol)
+    assert aoa.g(y) == pytest.approx(df * (y - x) + 0.5 * t1_rec.ddg(x) * (y - x) ** 2)
+    assert (f + np.sum(aoa.g(y))) == pytest.approx(f + np.sum(df * (y - x) + 0.5 * t1_rec.ddg(x) * (y - x) ** 2), tol)
+
+
+def non_spherical(delta_dg, delta_x): return delta_dg / delta_x
+
+
+def test_non_spherical(dx=1, tol=1e-4):
+    prob = Dummy(4)
+    x0 = prob.x0
+    df0 = prob.dg(x0)
+
+    aoa = DQA()
+    aoa.update(x0, df0)
+
+    x1 = x0 + dx
+
+    df1 = prob.dg(x1)
+    aoa.update(x1, df1, ddg0=non_spherical(df0 - df1, x0 - x1))
+
+    assert aoa.ddg0 == pytest.approx((df0 - df1) / (x0 - x1), tol)
+    assert aoa.dg(x0) == pytest.approx(df0, tol)
+
+
+def spherical(delta_g, delta_x, dg): return 2 * (delta_g - dg @ delta_x) / np.sum(delta_x ** 2)
+
+
+def test_spherical(dx=1, tol=1e-4):
+    prob = Dummy(4)
+    x0 = prob.x0
+    f0 = prob.g(x0)
+    df0 = prob.dg(x0)
+
+    aoa = DQA()
+    aoa.update(x0, df0)
+
+    x1 = x0 + dx
+
+    f1 = prob.g(x1)
+    df1 = prob.dg(x1)
+
+    aoa.update(x1, df1, ddg0=spherical(f0 - f1, x0 - x1, df1))
+
+    assert aoa.ddg0 == pytest.approx(2 * (f0 - f1 - df1 @ (x0 - x1)) / np.sum((x0 - x1) ** 2), tol)
+    # assert (f1 + np.sum(aoa.g(x0))) == pytest.approx(f0, tol)
+
+
+if __name__ == "__main__":
+    test_aoa_rec()
+    test_non_spherical()
+    test_spherical()
diff --git a/tests/mappings/test_mapping.py b/tests/mappings/test_mapping.py
new file mode 100644
index 00000000..3c134ac1
--- /dev/null
+++ b/tests/mappings/test_mapping.py
@@ -0,0 +1,59 @@
+import numpy as np
+import pytest
+from sao.mappings.mapping import Mapping, Linear
+
+
+class Dummy(Mapping):
+    g0 = 0
+    dg0 = 1
+    ddg0 = 0
+
+    def _update(self, x0, dg0, ddg0=0):
+        self.g0 += x0
+        self.dg0 *= dg0
+        # I don't like this line, to the users adding this line does not add anything
+
+    def _g(self, x): return self.g0 + x
+
+    def _dg(self, x): return self.dg0 * x
+
+    def _ddg(self, x): return self.ddg0 + x
+
+
+class TestMapping:
+    def test_init(self):
+        """
+        Test empty initialization sets child to Linear
+        """
+
+        mymap = Dummy(Dummy())
+
+        assert isinstance(mymap, Dummy)
+        assert isinstance(mymap.child, Dummy)
+        assert isinstance(mymap.child.child, Linear)
+
+    @pytest.mark.parametrize('x', [-1, 0, 1, 2, 10, 1000])
+    def test_g_dg_ddg(self, x):
+        """
+        Test the responses of a chain of Mappings using Dummy.
+        """
+
+        mymap2 = Dummy(Dummy())
+
+        # f[x] = f[g[x]]
+        assert mymap2.g(x) == x
+
+        # f'[x] = f'[g[x]]*g'[x]
+        assert mymap2.dg(x) == x ** 2
+
+        # f''[x] = f''[g[x]]*(g'[x])^2 + f'[g[x]]*g''[x]
+        assert mymap2.ddg(x) == x ** 3 + x ** 2
+
+    @pytest.mark.parametrize('x', [-1, 0, 1, 2, 10, 1000])
+    def test_update(self, x):
+        mymap = Dummy(Dummy())
+        mymap.update(0, 1)
+
+        assert mymap.g(x) == x
+        assert mymap.dg(x) == x ** 2
+        assert mymap.ddg(x) == x ** 3 + x ** 2
diff --git a/tests/mappings/test_mixed_mapping.py b/tests/mappings/test_mixed_mapping.py
new file mode 100644
index 00000000..57bacd3c
--- /dev/null
+++ b/tests/mappings/test_mixed_mapping.py
@@ -0,0 +1,33 @@
+from problems.n_dim.square import Square
+from sao.mappings.change_of_variable import Exponential as Exp
+from sao.mappings.approximations import DiagonalQuadraticApproximation as DQA
+import pytest
+from sao.mappings.mapping import MixedMapping as MM
+import numpy as np
+
+
+def test_mmn(dx=1, tol=1e-4):
+    prob = Square(3)
+    x = prob.x0
+    f = prob.g(x)
+    df = prob.dg(x)
+    ddf = prob.ddg(x)
+
+    mm = MM(prob.n, prob.m + 1)
+    mm[np.array([0]), [0, 1, 2]] = Exp(-1)
+    mm[[1], [2]] = Exp(-2)
+
+    mm.update(prob.x0, df, ddf)
+
+    aoa = DQA()
+    aoa.update(x, df, ddg0=mm.ddg(x))
+
+    assert mm.g(x)[1, [0, 1]] == pytest.approx(x[0:2], tol)
+    assert mm.g(x)[0, [0, 1, 2]] == pytest.approx(1 / x, tol)
+    assert mm.g(x)[1, 2] == pytest.approx(1 / (x[2]) ** 2, tol)
+
+    assert aoa.ddg0 == pytest.approx(mm.ddg(x), tol)
+
+
+if __name__ == "__main__":
+    test_mmn()