From 7647ee8f606479655c6a5d9381cae811de1dba7e Mon Sep 17 00:00:00 2001
From: philippeitis <33013301+philippeitis@users.noreply.github.com>
Date: Sun, 15 Dec 2019 01:26:08 -0800
Subject: [PATCH 01/10] Create triangulation.c

This provides faster implementations for fast_norm, fast_2d_circumcircle, fast_3d_circumcircle, and fast_2d_point_in_simplex.
---
 learner/triangulation.c | 272 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 272 insertions(+)
 create mode 100644 learner/triangulation.c

diff --git a/learner/triangulation.c b/learner/triangulation.c
new file mode 100644
index 000000000..44a391557
--- /dev/null
+++ b/learner/triangulation.c
@@ -0,0 +1,272 @@
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include <math.h>
+#include <stdio.h>
+#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+
+static PyObject *
+fast_norm(PyObject *self, PyObject *args)
+{
+    PyObject* vect;
+	if (!PyArg_ParseTuple(args, "O", &vect))
+		return NULL;
+
+    if (PyList_Check(vect)) {
+        Py_ssize_t numElements = PyObject_Length(vect);
+        if (numElements <= 0) {
+            PyErr_SetString(PyExc_ValueError, "fast_norm requires a list of length 1 or greater.");
+            return NULL;
+        }
+        double sum = 0;
+        for (Py_ssize_t i = 0; i < numElements; ++i) {
+            double val = PyFloat_AsDouble(PyList_GetItem(vect, i));
+            if (PyErr_Occurred()) {
+                return NULL;
+            }
+            sum += val * val;
+        }
+        return Py_BuildValue("f", sqrt(sum));
+    }
+
+    if (PyTuple_Check(vect)) {
+        Py_ssize_t numElements = PyTuple_Size(vect);
+        if (numElements <= 0) {
+            PyErr_SetString(PyExc_ValueError, "fast_norm requires a list of length 1 or greater.");
+            return NULL;
+        }
+        double sum = 0;
+        for (Py_ssize_t i = 0; i < numElements; ++i) {
+            double val = PyFloat_AsDouble(PyTuple_GetItem(vect, i));
+            if (PyErr_Occurred()) {
+                return NULL;
+            }
+            sum += val * val;
+        }
+        return Py_BuildValue("f", sqrt(sum));
+    }
+
+    if (PyArray_Check(vect)) {
+        PyObject *npArray = PyArray_FROM_OTF(vect, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+        if (npArray == NULL)
+            return NULL;
+        int nd = PyArray_NDIM(npArray);
+        if (nd != 1)
+            return NULL;
+        npy_intp * shape = PyArray_DIMS(npArray);
+        Py_ssize_t numElements = PyArray_SIZE(npArray);
+        double* data = (double *) PyArray_DATA(npArray);
+        double sum = 0;
+        for (Py_ssize_t i = 0; i < numElements; ++i) {
+            double val = data[i];
+            sum += val * val;
+        }
+        return Py_BuildValue("f", sqrt(sum));
+    }
+
+}
+
+static int get_tuple_elems(PyObject * tuple, double * first, double * second) {
+    PyObject * elem = PyTuple_GetItem(tuple, 0);
+    if (elem == NULL)
+        return 0;
+    *first = PyFloat_AsDouble(elem);
+    elem = PyTuple_GetItem(tuple, 1);
+    if (elem == NULL)
+        return 0;
+    *second = PyFloat_AsDouble(elem);
+    return 1;
+}
+
+static int get_triple_tuple_elems(PyObject * tuple, double * first, double * second, double* third) {
+    PyObject * elem = PyTuple_GetItem(tuple, 0);
+    if (elem == NULL)
+        return 0;
+    *first = PyFloat_AsDouble(elem);
+    elem = PyTuple_GetItem(tuple, 1);
+    if (elem == NULL)
+        return 0;
+    *second = PyFloat_AsDouble(elem);
+    elem = PyTuple_GetItem(tuple, 2);
+    if (elem == NULL)
+        return 0;
+    *third = PyFloat_AsDouble(elem);
+    return 1;
+}
+
+static PyObject *
+fast_2d_point_in_simplex(PyObject *self, PyObject *args, PyObject *kwargs)
+{
+    PyObject* point, * simplex;
+    double eps = 0.00000001;
+
+    // SystemError: bad format char passed to Py_BuildValue ??
+	static char *kwlist[] = {"point", "simplex", "eps", NULL};
+
+	if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|d", kwlist, &point, &simplex, &eps))
+		return NULL;
+
+
+    if (!PyTuple_Check(point))
+        return NULL;
+
+    if (!PyList_Check(simplex))
+        return NULL;
+
+    double px, py, p0x, p0y, p1x, p1y, p2x, p2y;
+    if (!get_tuple_elems(point, &px, &py))
+        return NULL;
+
+    point = PyList_GetItem(simplex, 0);
+    if (point == NULL)
+        return NULL;
+
+    if (!get_tuple_elems(point, &p0x, &p0y))
+        return NULL;
+
+    point = PyList_GetItem(simplex, 1);
+    if (point == NULL)
+        return NULL;
+
+    if (!get_tuple_elems(point, &p1x, &p1y))
+        return NULL;
+
+    point = PyList_GetItem(simplex, 2);
+    if (point == NULL)
+        return NULL;
+
+    if (!get_tuple_elems(point, &p2x, &p2y))
+        return NULL;
+
+    double area = 0.5 * (-p1y * p2x + p0y * (p2x - p1x) + p1x * p2y + p0x * (p1y - p2y));
+    double s = 1 / (2 * area) * (p0y * p2x + (p2y - p0y) * px - p0x * p2y + (p0x - p2x) * py);
+    if ((s < -eps) || (s > 1 + eps))
+        return Py_BuildValue("O", Py_False);
+    double t = 1 / (2 * area) * (p0x * p1y + (p0y - p1y) * px - p0y * p1x + (p1x - p0x) * py);
+    return Py_BuildValue("O", (t >= -eps) && (s + t <= 1 + eps) ? Py_True : Py_False);
+}
+
+static PyObject *
+fast_2d_circumcircle(PyObject *self, PyObject *args)
+{
+    PyObject* points;
+	if (!PyArg_ParseTuple(args, "O", &points))
+		return NULL;
+
+    if (PyList_Check(points)) {
+        Py_ssize_t numElements = PyObject_Length(points);
+        if (numElements <= 0) {
+            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a list of length 1 or greater.");
+            return NULL;
+        }
+        double x0, y0, x1, y1, x2, y2;
+        if (!get_tuple_elems(PyList_GetItem(points, 0), &x0, &y0))
+            return NULL;
+        if (!get_tuple_elems(PyList_GetItem(points, 1), &x1, &y1))
+            return NULL;
+        if (!get_tuple_elems(PyList_GetItem(points, 2), &x2, &y2))
+            return NULL;
+
+        x1 -= x0;
+        y1 -= y0;
+
+        x2 -= x0;
+        y2 -= y0;
+
+        double l1 = x1 * x1 + y1 * y1;
+        double l2 = x2 * x2 + y2 * y2;
+        double dx = l1 * y2 - l2 * y1;
+        double dy = -l1 * x2 + l2 * x1;
+
+        double aa = 2 * (x1 * y2 - x2 * y1);
+        double x = dx / aa;
+        double y = dy / aa;
+        return Py_BuildValue("(ff)f", x + x0, y + y0, sqrt(x * x + y * y));
+    }
+    return NULL;
+}
+
+static PyObject *
+fast_3d_circumcircle(PyObject *self, PyObject *args)
+{
+    PyObject* points;
+	if (!PyArg_ParseTuple(args, "O", &points))
+		return NULL;
+
+    if (PyList_Check(points)) {
+        Py_ssize_t numElements = PyObject_Length(points);
+        if (numElements <= 0) {
+            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a list of length 1 or greater.");
+            return NULL;
+        }
+        double x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3;
+        if (!get_triple_tuple_elems(PyList_GetItem(points, 0), &x0, &y0, &z0))
+            return NULL;
+        if (!get_triple_tuple_elems(PyList_GetItem(points, 1), &x1, &y1, &z1))
+            return NULL;
+        if (!get_triple_tuple_elems(PyList_GetItem(points, 2), &x2, &y2, &z2))
+            return NULL;
+        if (!get_triple_tuple_elems(PyList_GetItem(points, 3), &x3, &y3, &z3))
+            return NULL;
+
+        x1 -= x0;
+        y1 -= y0;
+        z1 -= z0;
+
+        x2 -= x0;
+        y2 -= y0;
+        z2 -= z0;
+
+        x3 -= x0;
+        y3 -= y0;
+        z3 -= z0;
+
+        double l1 = x1 * x1 + y1 * y1 + z1 * z1;
+        double l2 = x2 * x2 + y2 * y2 + z2 * z2;
+        double l3 = x3 * x3 + y3 * y3 + z3 * z3;
+        double dx = l1 * (y2 * z3 - z2 * y3) - l2 * (y1 * z3 - z1 * y3) + l3 * (y1 * z2 - z1 * y2);
+        double dy = -l1 * (x2 * z3 - z2 * x3) - l2 * (x1 * z3 - z1 * x3) + l3 * (x1 * z2 - z1 * x2);
+        double dz = l1 * (x2 * y3 - y2 * x3) - l2 * (x1 * y3 - y1 * x3) + l3 * (x1 * y2 - y1 * x2);
+
+        double aa = 2 * (x1 * (y2 * z3 - z2 * y3) - x2 * (y1 * z3 - z1 * y3) + x3 * (y1 * z2 - z1 * y2));
+        double x = dx / aa;
+        double y = dy / aa;
+        double z = dz / aa;
+        return Py_BuildValue("(fff)f", x + x0, y + y0, z + z0, sqrt(x * x + y * y + z * z));
+    }
+    return NULL;
+}
+
+
+static PyMethodDef triangulation_functions[] = {
+    {"fast_norm", fast_norm, METH_VARARGS,
+		"Returns the norm of the given array."},
+    {"fast_2d_circumcircle", fast_2d_circumcircle, METH_VARARGS,
+		"Returns the norm of the given array."},
+    {"fast_3d_circumcircle", fast_3d_circumcircle, METH_VARARGS,
+		"Returns the norm of the given array."},
+    {"fast_2d_point_in_simplex", (PyCFunction) fast_2d_point_in_simplex, METH_VARARGS | METH_KEYWORDS,
+		"Refer to docs."},
+	{NULL}
+};
+
+PyDoc_STRVAR(module_doc, "Fast implementations of triangulation functions.");
+
+static struct PyModuleDef moduledef = {
+	PyModuleDef_HEAD_INIT,
+	"triangulation",
+	module_doc,
+	-1,                 /* m_size */
+	triangulation_functions,   /* m_methods */
+	NULL,               /* m_reload (unused) */
+	NULL,               /* m_traverse */
+	NULL,               /* m_clear */
+	NULL                /* m_free */
+};
+
+PyObject *
+PyInit_triangulation(void)
+{
+    if(PyArray_API == NULL)
+        import_array();
+    return PyModule_Create(&moduledef);
+}

From 2c4f18690e191e5cf69a7a8b9fe7fb9dbf1ba40d Mon Sep 17 00:00:00 2001
From: philippeitis <33013301+philippeitis@users.noreply.github.com>
Date: Sun, 15 Dec 2019 01:28:33 -0800
Subject: [PATCH 02/10] Create triangulation.c

triangulation.c provides fast implementations for some of the functions in learner/triangulation.py
---
 adaptive/learner/triangulation.c | 272 +++++++++++++++++++++++++++++++
 1 file changed, 272 insertions(+)
 create mode 100644 adaptive/learner/triangulation.c

diff --git a/adaptive/learner/triangulation.c b/adaptive/learner/triangulation.c
new file mode 100644
index 000000000..44a391557
--- /dev/null
+++ b/adaptive/learner/triangulation.c
@@ -0,0 +1,272 @@
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include <math.h>
+#include <stdio.h>
+#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+
+static PyObject *
+fast_norm(PyObject *self, PyObject *args)
+{
+    PyObject* vect;
+	if (!PyArg_ParseTuple(args, "O", &vect))
+		return NULL;
+
+    if (PyList_Check(vect)) {
+        Py_ssize_t numElements = PyObject_Length(vect);
+        if (numElements <= 0) {
+            PyErr_SetString(PyExc_ValueError, "fast_norm requires a list of length 1 or greater.");
+            return NULL;
+        }
+        double sum = 0;
+        for (Py_ssize_t i = 0; i < numElements; ++i) {
+            double val = PyFloat_AsDouble(PyList_GetItem(vect, i));
+            if (PyErr_Occurred()) {
+                return NULL;
+            }
+            sum += val * val;
+        }
+        return Py_BuildValue("f", sqrt(sum));
+    }
+
+    if (PyTuple_Check(vect)) {
+        Py_ssize_t numElements = PyTuple_Size(vect);
+        if (numElements <= 0) {
+            PyErr_SetString(PyExc_ValueError, "fast_norm requires a list of length 1 or greater.");
+            return NULL;
+        }
+        double sum = 0;
+        for (Py_ssize_t i = 0; i < numElements; ++i) {
+            double val = PyFloat_AsDouble(PyTuple_GetItem(vect, i));
+            if (PyErr_Occurred()) {
+                return NULL;
+            }
+            sum += val * val;
+        }
+        return Py_BuildValue("f", sqrt(sum));
+    }
+
+    if (PyArray_Check(vect)) {
+        PyObject *npArray = PyArray_FROM_OTF(vect, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+        if (npArray == NULL)
+            return NULL;
+        int nd = PyArray_NDIM(npArray);
+        if (nd != 1)
+            return NULL;
+        npy_intp * shape = PyArray_DIMS(npArray);
+        Py_ssize_t numElements = PyArray_SIZE(npArray);
+        double* data = (double *) PyArray_DATA(npArray);
+        double sum = 0;
+        for (Py_ssize_t i = 0; i < numElements; ++i) {
+            double val = data[i];
+            sum += val * val;
+        }
+        return Py_BuildValue("f", sqrt(sum));
+    }
+
+}
+
+static int get_tuple_elems(PyObject * tuple, double * first, double * second) {
+    PyObject * elem = PyTuple_GetItem(tuple, 0);
+    if (elem == NULL)
+        return 0;
+    *first = PyFloat_AsDouble(elem);
+    elem = PyTuple_GetItem(tuple, 1);
+    if (elem == NULL)
+        return 0;
+    *second = PyFloat_AsDouble(elem);
+    return 1;
+}
+
+static int get_triple_tuple_elems(PyObject * tuple, double * first, double * second, double* third) {
+    PyObject * elem = PyTuple_GetItem(tuple, 0);
+    if (elem == NULL)
+        return 0;
+    *first = PyFloat_AsDouble(elem);
+    elem = PyTuple_GetItem(tuple, 1);
+    if (elem == NULL)
+        return 0;
+    *second = PyFloat_AsDouble(elem);
+    elem = PyTuple_GetItem(tuple, 2);
+    if (elem == NULL)
+        return 0;
+    *third = PyFloat_AsDouble(elem);
+    return 1;
+}
+
+static PyObject *
+fast_2d_point_in_simplex(PyObject *self, PyObject *args, PyObject *kwargs)
+{
+    PyObject* point, * simplex;
+    double eps = 0.00000001;
+
+    // SystemError: bad format char passed to Py_BuildValue ??
+	static char *kwlist[] = {"point", "simplex", "eps", NULL};
+
+	if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|d", kwlist, &point, &simplex, &eps))
+		return NULL;
+
+
+    if (!PyTuple_Check(point))
+        return NULL;
+
+    if (!PyList_Check(simplex))
+        return NULL;
+
+    double px, py, p0x, p0y, p1x, p1y, p2x, p2y;
+    if (!get_tuple_elems(point, &px, &py))
+        return NULL;
+
+    point = PyList_GetItem(simplex, 0);
+    if (point == NULL)
+        return NULL;
+
+    if (!get_tuple_elems(point, &p0x, &p0y))
+        return NULL;
+
+    point = PyList_GetItem(simplex, 1);
+    if (point == NULL)
+        return NULL;
+
+    if (!get_tuple_elems(point, &p1x, &p1y))
+        return NULL;
+
+    point = PyList_GetItem(simplex, 2);
+    if (point == NULL)
+        return NULL;
+
+    if (!get_tuple_elems(point, &p2x, &p2y))
+        return NULL;
+
+    double area = 0.5 * (-p1y * p2x + p0y * (p2x - p1x) + p1x * p2y + p0x * (p1y - p2y));
+    double s = 1 / (2 * area) * (p0y * p2x + (p2y - p0y) * px - p0x * p2y + (p0x - p2x) * py);
+    if ((s < -eps) || (s > 1 + eps))
+        return Py_BuildValue("O", Py_False);
+    double t = 1 / (2 * area) * (p0x * p1y + (p0y - p1y) * px - p0y * p1x + (p1x - p0x) * py);
+    return Py_BuildValue("O", (t >= -eps) && (s + t <= 1 + eps) ? Py_True : Py_False);
+}
+
+static PyObject *
+fast_2d_circumcircle(PyObject *self, PyObject *args)
+{
+    PyObject* points;
+	if (!PyArg_ParseTuple(args, "O", &points))
+		return NULL;
+
+    if (PyList_Check(points)) {
+        Py_ssize_t numElements = PyObject_Length(points);
+        if (numElements <= 0) {
+            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a list of length 1 or greater.");
+            return NULL;
+        }
+        double x0, y0, x1, y1, x2, y2;
+        if (!get_tuple_elems(PyList_GetItem(points, 0), &x0, &y0))
+            return NULL;
+        if (!get_tuple_elems(PyList_GetItem(points, 1), &x1, &y1))
+            return NULL;
+        if (!get_tuple_elems(PyList_GetItem(points, 2), &x2, &y2))
+            return NULL;
+
+        x1 -= x0;
+        y1 -= y0;
+
+        x2 -= x0;
+        y2 -= y0;
+
+        double l1 = x1 * x1 + y1 * y1;
+        double l2 = x2 * x2 + y2 * y2;
+        double dx = l1 * y2 - l2 * y1;
+        double dy = -l1 * x2 + l2 * x1;
+
+        double aa = 2 * (x1 * y2 - x2 * y1);
+        double x = dx / aa;
+        double y = dy / aa;
+        return Py_BuildValue("(ff)f", x + x0, y + y0, sqrt(x * x + y * y));
+    }
+    return NULL;
+}
+
+static PyObject *
+fast_3d_circumcircle(PyObject *self, PyObject *args)
+{
+    PyObject* points;
+	if (!PyArg_ParseTuple(args, "O", &points))
+		return NULL;
+
+    if (PyList_Check(points)) {
+        Py_ssize_t numElements = PyObject_Length(points);
+        if (numElements <= 0) {
+            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a list of length 1 or greater.");
+            return NULL;
+        }
+        double x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3;
+        if (!get_triple_tuple_elems(PyList_GetItem(points, 0), &x0, &y0, &z0))
+            return NULL;
+        if (!get_triple_tuple_elems(PyList_GetItem(points, 1), &x1, &y1, &z1))
+            return NULL;
+        if (!get_triple_tuple_elems(PyList_GetItem(points, 2), &x2, &y2, &z2))
+            return NULL;
+        if (!get_triple_tuple_elems(PyList_GetItem(points, 3), &x3, &y3, &z3))
+            return NULL;
+
+        x1 -= x0;
+        y1 -= y0;
+        z1 -= z0;
+
+        x2 -= x0;
+        y2 -= y0;
+        z2 -= z0;
+
+        x3 -= x0;
+        y3 -= y0;
+        z3 -= z0;
+
+        double l1 = x1 * x1 + y1 * y1 + z1 * z1;
+        double l2 = x2 * x2 + y2 * y2 + z2 * z2;
+        double l3 = x3 * x3 + y3 * y3 + z3 * z3;
+        double dx = l1 * (y2 * z3 - z2 * y3) - l2 * (y1 * z3 - z1 * y3) + l3 * (y1 * z2 - z1 * y2);
+        double dy = -l1 * (x2 * z3 - z2 * x3) - l2 * (x1 * z3 - z1 * x3) + l3 * (x1 * z2 - z1 * x2);
+        double dz = l1 * (x2 * y3 - y2 * x3) - l2 * (x1 * y3 - y1 * x3) + l3 * (x1 * y2 - y1 * x2);
+
+        double aa = 2 * (x1 * (y2 * z3 - z2 * y3) - x2 * (y1 * z3 - z1 * y3) + x3 * (y1 * z2 - z1 * y2));
+        double x = dx / aa;
+        double y = dy / aa;
+        double z = dz / aa;
+        return Py_BuildValue("(fff)f", x + x0, y + y0, z + z0, sqrt(x * x + y * y + z * z));
+    }
+    return NULL;
+}
+
+
+static PyMethodDef triangulation_functions[] = {
+    {"fast_norm", fast_norm, METH_VARARGS,
+		"Returns the norm of the given array."},
+    {"fast_2d_circumcircle", fast_2d_circumcircle, METH_VARARGS,
+		"Returns the norm of the given array."},
+    {"fast_3d_circumcircle", fast_3d_circumcircle, METH_VARARGS,
+		"Returns the norm of the given array."},
+    {"fast_2d_point_in_simplex", (PyCFunction) fast_2d_point_in_simplex, METH_VARARGS | METH_KEYWORDS,
+		"Refer to docs."},
+	{NULL}
+};
+
+PyDoc_STRVAR(module_doc, "Fast implementations of triangulation functions.");
+
+static struct PyModuleDef moduledef = {
+	PyModuleDef_HEAD_INIT,
+	"triangulation",
+	module_doc,
+	-1,                 /* m_size */
+	triangulation_functions,   /* m_methods */
+	NULL,               /* m_reload (unused) */
+	NULL,               /* m_traverse */
+	NULL,               /* m_clear */
+	NULL                /* m_free */
+};
+
+PyObject *
+PyInit_triangulation(void)
+{
+    if(PyArray_API == NULL)
+        import_array();
+    return PyModule_Create(&moduledef);
+}

From d39e319abcf92ddea889b45288960757aefb196b Mon Sep 17 00:00:00 2001
From: philippeitis <33013301+philippeitis@users.noreply.github.com>
Date: Sun, 15 Dec 2019 01:29:19 -0800
Subject: [PATCH 03/10] Delete triangulation.c

Wrong directory.
---
 learner/triangulation.c | 272 ----------------------------------------
 1 file changed, 272 deletions(-)
 delete mode 100644 learner/triangulation.c

diff --git a/learner/triangulation.c b/learner/triangulation.c
deleted file mode 100644
index 44a391557..000000000
--- a/learner/triangulation.c
+++ /dev/null
@@ -1,272 +0,0 @@
-#include "Python.h"
-#include "numpy/arrayobject.h"
-#include <math.h>
-#include <stdio.h>
-#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
-
-static PyObject *
-fast_norm(PyObject *self, PyObject *args)
-{
-    PyObject* vect;
-	if (!PyArg_ParseTuple(args, "O", &vect))
-		return NULL;
-
-    if (PyList_Check(vect)) {
-        Py_ssize_t numElements = PyObject_Length(vect);
-        if (numElements <= 0) {
-            PyErr_SetString(PyExc_ValueError, "fast_norm requires a list of length 1 or greater.");
-            return NULL;
-        }
-        double sum = 0;
-        for (Py_ssize_t i = 0; i < numElements; ++i) {
-            double val = PyFloat_AsDouble(PyList_GetItem(vect, i));
-            if (PyErr_Occurred()) {
-                return NULL;
-            }
-            sum += val * val;
-        }
-        return Py_BuildValue("f", sqrt(sum));
-    }
-
-    if (PyTuple_Check(vect)) {
-        Py_ssize_t numElements = PyTuple_Size(vect);
-        if (numElements <= 0) {
-            PyErr_SetString(PyExc_ValueError, "fast_norm requires a list of length 1 or greater.");
-            return NULL;
-        }
-        double sum = 0;
-        for (Py_ssize_t i = 0; i < numElements; ++i) {
-            double val = PyFloat_AsDouble(PyTuple_GetItem(vect, i));
-            if (PyErr_Occurred()) {
-                return NULL;
-            }
-            sum += val * val;
-        }
-        return Py_BuildValue("f", sqrt(sum));
-    }
-
-    if (PyArray_Check(vect)) {
-        PyObject *npArray = PyArray_FROM_OTF(vect, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
-        if (npArray == NULL)
-            return NULL;
-        int nd = PyArray_NDIM(npArray);
-        if (nd != 1)
-            return NULL;
-        npy_intp * shape = PyArray_DIMS(npArray);
-        Py_ssize_t numElements = PyArray_SIZE(npArray);
-        double* data = (double *) PyArray_DATA(npArray);
-        double sum = 0;
-        for (Py_ssize_t i = 0; i < numElements; ++i) {
-            double val = data[i];
-            sum += val * val;
-        }
-        return Py_BuildValue("f", sqrt(sum));
-    }
-
-}
-
-static int get_tuple_elems(PyObject * tuple, double * first, double * second) {
-    PyObject * elem = PyTuple_GetItem(tuple, 0);
-    if (elem == NULL)
-        return 0;
-    *first = PyFloat_AsDouble(elem);
-    elem = PyTuple_GetItem(tuple, 1);
-    if (elem == NULL)
-        return 0;
-    *second = PyFloat_AsDouble(elem);
-    return 1;
-}
-
-static int get_triple_tuple_elems(PyObject * tuple, double * first, double * second, double* third) {
-    PyObject * elem = PyTuple_GetItem(tuple, 0);
-    if (elem == NULL)
-        return 0;
-    *first = PyFloat_AsDouble(elem);
-    elem = PyTuple_GetItem(tuple, 1);
-    if (elem == NULL)
-        return 0;
-    *second = PyFloat_AsDouble(elem);
-    elem = PyTuple_GetItem(tuple, 2);
-    if (elem == NULL)
-        return 0;
-    *third = PyFloat_AsDouble(elem);
-    return 1;
-}
-
-static PyObject *
-fast_2d_point_in_simplex(PyObject *self, PyObject *args, PyObject *kwargs)
-{
-    PyObject* point, * simplex;
-    double eps = 0.00000001;
-
-    // SystemError: bad format char passed to Py_BuildValue ??
-	static char *kwlist[] = {"point", "simplex", "eps", NULL};
-
-	if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|d", kwlist, &point, &simplex, &eps))
-		return NULL;
-
-
-    if (!PyTuple_Check(point))
-        return NULL;
-
-    if (!PyList_Check(simplex))
-        return NULL;
-
-    double px, py, p0x, p0y, p1x, p1y, p2x, p2y;
-    if (!get_tuple_elems(point, &px, &py))
-        return NULL;
-
-    point = PyList_GetItem(simplex, 0);
-    if (point == NULL)
-        return NULL;
-
-    if (!get_tuple_elems(point, &p0x, &p0y))
-        return NULL;
-
-    point = PyList_GetItem(simplex, 1);
-    if (point == NULL)
-        return NULL;
-
-    if (!get_tuple_elems(point, &p1x, &p1y))
-        return NULL;
-
-    point = PyList_GetItem(simplex, 2);
-    if (point == NULL)
-        return NULL;
-
-    if (!get_tuple_elems(point, &p2x, &p2y))
-        return NULL;
-
-    double area = 0.5 * (-p1y * p2x + p0y * (p2x - p1x) + p1x * p2y + p0x * (p1y - p2y));
-    double s = 1 / (2 * area) * (p0y * p2x + (p2y - p0y) * px - p0x * p2y + (p0x - p2x) * py);
-    if ((s < -eps) || (s > 1 + eps))
-        return Py_BuildValue("O", Py_False);
-    double t = 1 / (2 * area) * (p0x * p1y + (p0y - p1y) * px - p0y * p1x + (p1x - p0x) * py);
-    return Py_BuildValue("O", (t >= -eps) && (s + t <= 1 + eps) ? Py_True : Py_False);
-}
-
-static PyObject *
-fast_2d_circumcircle(PyObject *self, PyObject *args)
-{
-    PyObject* points;
-	if (!PyArg_ParseTuple(args, "O", &points))
-		return NULL;
-
-    if (PyList_Check(points)) {
-        Py_ssize_t numElements = PyObject_Length(points);
-        if (numElements <= 0) {
-            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a list of length 1 or greater.");
-            return NULL;
-        }
-        double x0, y0, x1, y1, x2, y2;
-        if (!get_tuple_elems(PyList_GetItem(points, 0), &x0, &y0))
-            return NULL;
-        if (!get_tuple_elems(PyList_GetItem(points, 1), &x1, &y1))
-            return NULL;
-        if (!get_tuple_elems(PyList_GetItem(points, 2), &x2, &y2))
-            return NULL;
-
-        x1 -= x0;
-        y1 -= y0;
-
-        x2 -= x0;
-        y2 -= y0;
-
-        double l1 = x1 * x1 + y1 * y1;
-        double l2 = x2 * x2 + y2 * y2;
-        double dx = l1 * y2 - l2 * y1;
-        double dy = -l1 * x2 + l2 * x1;
-
-        double aa = 2 * (x1 * y2 - x2 * y1);
-        double x = dx / aa;
-        double y = dy / aa;
-        return Py_BuildValue("(ff)f", x + x0, y + y0, sqrt(x * x + y * y));
-    }
-    return NULL;
-}
-
-static PyObject *
-fast_3d_circumcircle(PyObject *self, PyObject *args)
-{
-    PyObject* points;
-	if (!PyArg_ParseTuple(args, "O", &points))
-		return NULL;
-
-    if (PyList_Check(points)) {
-        Py_ssize_t numElements = PyObject_Length(points);
-        if (numElements <= 0) {
-            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a list of length 1 or greater.");
-            return NULL;
-        }
-        double x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3;
-        if (!get_triple_tuple_elems(PyList_GetItem(points, 0), &x0, &y0, &z0))
-            return NULL;
-        if (!get_triple_tuple_elems(PyList_GetItem(points, 1), &x1, &y1, &z1))
-            return NULL;
-        if (!get_triple_tuple_elems(PyList_GetItem(points, 2), &x2, &y2, &z2))
-            return NULL;
-        if (!get_triple_tuple_elems(PyList_GetItem(points, 3), &x3, &y3, &z3))
-            return NULL;
-
-        x1 -= x0;
-        y1 -= y0;
-        z1 -= z0;
-
-        x2 -= x0;
-        y2 -= y0;
-        z2 -= z0;
-
-        x3 -= x0;
-        y3 -= y0;
-        z3 -= z0;
-
-        double l1 = x1 * x1 + y1 * y1 + z1 * z1;
-        double l2 = x2 * x2 + y2 * y2 + z2 * z2;
-        double l3 = x3 * x3 + y3 * y3 + z3 * z3;
-        double dx = l1 * (y2 * z3 - z2 * y3) - l2 * (y1 * z3 - z1 * y3) + l3 * (y1 * z2 - z1 * y2);
-        double dy = -l1 * (x2 * z3 - z2 * x3) - l2 * (x1 * z3 - z1 * x3) + l3 * (x1 * z2 - z1 * x2);
-        double dz = l1 * (x2 * y3 - y2 * x3) - l2 * (x1 * y3 - y1 * x3) + l3 * (x1 * y2 - y1 * x2);
-
-        double aa = 2 * (x1 * (y2 * z3 - z2 * y3) - x2 * (y1 * z3 - z1 * y3) + x3 * (y1 * z2 - z1 * y2));
-        double x = dx / aa;
-        double y = dy / aa;
-        double z = dz / aa;
-        return Py_BuildValue("(fff)f", x + x0, y + y0, z + z0, sqrt(x * x + y * y + z * z));
-    }
-    return NULL;
-}
-
-
-static PyMethodDef triangulation_functions[] = {
-    {"fast_norm", fast_norm, METH_VARARGS,
-		"Returns the norm of the given array."},
-    {"fast_2d_circumcircle", fast_2d_circumcircle, METH_VARARGS,
-		"Returns the norm of the given array."},
-    {"fast_3d_circumcircle", fast_3d_circumcircle, METH_VARARGS,
-		"Returns the norm of the given array."},
-    {"fast_2d_point_in_simplex", (PyCFunction) fast_2d_point_in_simplex, METH_VARARGS | METH_KEYWORDS,
-		"Refer to docs."},
-	{NULL}
-};
-
-PyDoc_STRVAR(module_doc, "Fast implementations of triangulation functions.");
-
-static struct PyModuleDef moduledef = {
-	PyModuleDef_HEAD_INIT,
-	"triangulation",
-	module_doc,
-	-1,                 /* m_size */
-	triangulation_functions,   /* m_methods */
-	NULL,               /* m_reload (unused) */
-	NULL,               /* m_traverse */
-	NULL,               /* m_clear */
-	NULL                /* m_free */
-};
-
-PyObject *
-PyInit_triangulation(void)
-{
-    if(PyArray_API == NULL)
-        import_array();
-    return PyModule_Create(&moduledef);
-}

From ca14591b48efbd0e665d18f7df8f9c3326d467a7 Mon Sep 17 00:00:00 2001
From: philippeitis <33013301+philippeitis@users.noreply.github.com>
Date: Sun, 15 Dec 2019 01:30:53 -0800
Subject: [PATCH 04/10] Update setup.py

On running pip install adaptive, adaptive.triangulation will also be set up.
---
 setup.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/setup.py b/setup.py
index 0b9bd474d..d58442790 100644
--- a/setup.py
+++ b/setup.py
@@ -1,6 +1,7 @@
 #!/usr/bin/env python3
 
 import sys
+import numpy as np
 
 from setuptools import find_packages, setup
 
@@ -43,6 +44,8 @@ def get_version_and_cmdclass(package_name):
     ]
 }
 
+triangulation_module = Extension("adaptive.triangulation", sources=['./adaptive/learner/triangulation.c'],
+                                 include_dirs=[np.get_include()])
 
 setup(
     name="adaptive",
@@ -63,4 +66,5 @@ def get_version_and_cmdclass(package_name):
     install_requires=install_requires,
     extras_require=extras_require,
     cmdclass=cmdclass,
+    ext_modules=[triangulation_module])
 )

From 448bf2c216cc361867618f9d236bd3608fe2e4c0 Mon Sep 17 00:00:00 2001
From: philippeitis <33013301+philippeitis@users.noreply.github.com>
Date: Sun, 15 Dec 2019 01:43:26 -0800
Subject: [PATCH 05/10] Update setup.py

Missed a bracket.
---
 setup.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.py b/setup.py
index d58442790..7e52d421c 100644
--- a/setup.py
+++ b/setup.py
@@ -66,5 +66,5 @@ def get_version_and_cmdclass(package_name):
     install_requires=install_requires,
     extras_require=extras_require,
     cmdclass=cmdclass,
-    ext_modules=[triangulation_module])
+    ext_modules=[triangulation_module]
 )

From 4d9e6b2a43fd9b854b6989374545863bc8cf8784 Mon Sep 17 00:00:00 2001
From: philippeitis <33013301+philippeitis@users.noreply.github.com>
Date: Sun, 15 Dec 2019 01:44:39 -0800
Subject: [PATCH 06/10] Update setup.py

Missing import.
---
 setup.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.py b/setup.py
index 7e52d421c..ff00c96ca 100644
--- a/setup.py
+++ b/setup.py
@@ -3,7 +3,7 @@
 import sys
 import numpy as np
 
-from setuptools import find_packages, setup
+from setuptools import Extension, find_packages, setup
 
 if sys.version_info < (3, 6):
     print("adaptive requires Python 3.6 or above.")

From 1b6ef17a7dd817c5021f70f0987881cacdee8834 Mon Sep 17 00:00:00 2001
From: Bas Nijholt <basnijholt@gmail.com>
Date: Sun, 15 Dec 2019 11:19:27 +0100
Subject: [PATCH 07/10] fix black and isort issue in setup.py

---
 setup.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/setup.py b/setup.py
index ff00c96ca..a533a2c79 100644
--- a/setup.py
+++ b/setup.py
@@ -1,8 +1,8 @@
 #!/usr/bin/env python3
 
 import sys
-import numpy as np
 
+import numpy as np
 from setuptools import Extension, find_packages, setup
 
 if sys.version_info < (3, 6):
@@ -44,8 +44,11 @@ def get_version_and_cmdclass(package_name):
     ]
 }
 
-triangulation_module = Extension("adaptive.triangulation", sources=['./adaptive/learner/triangulation.c'],
-                                 include_dirs=[np.get_include()])
+triangulation_module = Extension(
+    "adaptive.triangulation",
+    sources=["./adaptive/learner/triangulation.c"],
+    include_dirs=[np.get_include()],
+)
 
 setup(
     name="adaptive",
@@ -66,5 +69,5 @@ def get_version_and_cmdclass(package_name):
     install_requires=install_requires,
     extras_require=extras_require,
     cmdclass=cmdclass,
-    ext_modules=[triangulation_module]
+    ext_modules=[triangulation_module],
 )

From 2f2813c7727aaa0940d066cf65c990c649ceb3de Mon Sep 17 00:00:00 2001
From: philippeitis <33013301+philippeitis@users.noreply.github.com>
Date: Sun, 15 Dec 2019 13:07:08 -0800
Subject: [PATCH 08/10] Update triangulation.c

Added more error strings and updated docstrings to include details about what the functions accept and when to use them.
Expanded fast_2d_circumcircle, fast_3d_circumcircle to accept tuples and appropriately sized numpy arrays.
---
 adaptive/learner/triangulation.c | 200 +++++++++++++++++++++++--------
 1 file changed, 147 insertions(+), 53 deletions(-)

diff --git a/adaptive/learner/triangulation.c b/adaptive/learner/triangulation.c
index 44a391557..1c81c5831 100644
--- a/adaptive/learner/triangulation.c
+++ b/adaptive/learner/triangulation.c
@@ -31,7 +31,7 @@ fast_norm(PyObject *self, PyObject *args)
     if (PyTuple_Check(vect)) {
         Py_ssize_t numElements = PyTuple_Size(vect);
         if (numElements <= 0) {
-            PyErr_SetString(PyExc_ValueError, "fast_norm requires a list of length 1 or greater.");
+            PyErr_SetString(PyExc_ValueError, "fast_norm requires a tuple of length 1 or greater.");
             return NULL;
         }
         double sum = 0;
@@ -105,17 +105,26 @@ fast_2d_point_in_simplex(PyObject *self, PyObject *args, PyObject *kwargs)
 	if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|d", kwlist, &point, &simplex, &eps))
 		return NULL;
 
+    if (!PyTuple_Check(point)) {
+        PyErr_SetString(PyExc_ValueError, "fast_2d_point_in_simplex requires a tuple as its first argument.");
+        return NULL;
+    }
 
-    if (!PyTuple_Check(point))
+    if (!PyList_Check(simplex)) {
+        PyErr_SetString(PyExc_ValueError, "fast_2d_point_in_simplex requires a python list as its argument.");
         return NULL;
+    }
 
-    if (!PyList_Check(simplex))
+    if (PyObject_Length(simplex) < 3) {
+        PyErr_SetString(PyExc_ValueError, "fast_2d_point_in_simplex requires simplex to contain at least 3 elements.");
         return NULL;
+    }
 
     double px, py, p0x, p0y, p1x, p1y, p2x, p2y;
     if (!get_tuple_elems(point, &px, &py))
         return NULL;
 
+
     point = PyList_GetItem(simplex, 0);
     if (point == NULL)
         return NULL;
@@ -152,37 +161,71 @@ fast_2d_circumcircle(PyObject *self, PyObject *args)
 	if (!PyArg_ParseTuple(args, "O", &points))
 		return NULL;
 
+    double x0, y0, x1, y1, x2, y2;
+
     if (PyList_Check(points)) {
         Py_ssize_t numElements = PyObject_Length(points);
-        if (numElements <= 0) {
-            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a list of length 1 or greater.");
+        if (numElements < 3) {
+            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a list of length 3 or greater.");
             return NULL;
         }
-        double x0, y0, x1, y1, x2, y2;
         if (!get_tuple_elems(PyList_GetItem(points, 0), &x0, &y0))
             return NULL;
         if (!get_tuple_elems(PyList_GetItem(points, 1), &x1, &y1))
             return NULL;
         if (!get_tuple_elems(PyList_GetItem(points, 2), &x2, &y2))
             return NULL;
+    } else if (PyTuple_Check(points)) {
+        Py_ssize_t numElements = PyTuple_Size(points);
+        if (numElements < 3) {
+            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a tuple of length 3 or greater.");
+            return NULL;
+        }
+
+        if (!get_tuple_elems(PyTuple_GetItem(points, 0), &x0, &y0))
+            return NULL;
+        if (!get_tuple_elems(PyTuple_GetItem(points, 1), &x1, &y1))
+            return NULL;
+        if (!get_tuple_elems(PyTuple_GetItem(points, 2), &x2, &y2))
+            return NULL;
+    } else if (PyArray_Check(points)) {
+        int dims = PyArray_NDIM(points);
+        if (dims != 2) {
+            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a two dimensional numpy array.");
+            return NULL;
+        }
+        npy_intp * shape = PyArray_DIMS(points);
+        if ((shape[0] < 3) && (shape[1] != 2)) {
+            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a numpy array of width 3, height 2.");
+            return NULL;
+        }
+        double* values = (double *) PyArray_DATA(points);
+        x0 = values[0];
+        y0 = values[1];
+        x1 = values[2];
+        y1 = values[3];
+        x2 = values[4];
+        y2 = values[5];
+    } else {
+        PyErr_SetString(PyExc_ValueError, "Points must be a list, tuple, or numpy array.");
+        return NULL;
+    }
 
-        x1 -= x0;
-        y1 -= y0;
+    x1 -= x0;
+    y1 -= y0;
 
-        x2 -= x0;
-        y2 -= y0;
+    x2 -= x0;
+    y2 -= y0;
 
-        double l1 = x1 * x1 + y1 * y1;
-        double l2 = x2 * x2 + y2 * y2;
-        double dx = l1 * y2 - l2 * y1;
-        double dy = -l1 * x2 + l2 * x1;
+    double l1 = x1 * x1 + y1 * y1;
+    double l2 = x2 * x2 + y2 * y2;
+    double dx = l1 * y2 - l2 * y1;
+    double dy = -l1 * x2 + l2 * x1;
 
-        double aa = 2 * (x1 * y2 - x2 * y1);
-        double x = dx / aa;
-        double y = dy / aa;
-        return Py_BuildValue("(ff)f", x + x0, y + y0, sqrt(x * x + y * y));
-    }
-    return NULL;
+    double aa = 2 * (x1 * y2 - x2 * y1);
+    double x = dx / aa;
+    double y = dy / aa;
+    return Py_BuildValue("(ff)f", x + x0, y + y0, sqrt(x * x + y * y));
 }
 
 static PyObject *
@@ -192,13 +235,14 @@ fast_3d_circumcircle(PyObject *self, PyObject *args)
 	if (!PyArg_ParseTuple(args, "O", &points))
 		return NULL;
 
+    double x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3;
+
     if (PyList_Check(points)) {
         Py_ssize_t numElements = PyObject_Length(points);
-        if (numElements <= 0) {
-            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a list of length 1 or greater.");
+        if (numElements < 4) {
+            PyErr_SetString(PyExc_ValueError, "fast_3d_circumcircle requires a list of length 4 or greater.");
             return NULL;
         }
-        double x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3;
         if (!get_triple_tuple_elems(PyList_GetItem(points, 0), &x0, &y0, &z0))
             return NULL;
         if (!get_triple_tuple_elems(PyList_GetItem(points, 1), &x1, &y1, &z1))
@@ -207,45 +251,95 @@ fast_3d_circumcircle(PyObject *self, PyObject *args)
             return NULL;
         if (!get_triple_tuple_elems(PyList_GetItem(points, 3), &x3, &y3, &z3))
             return NULL;
-
-        x1 -= x0;
-        y1 -= y0;
-        z1 -= z0;
-
-        x2 -= x0;
-        y2 -= y0;
-        z2 -= z0;
-
-        x3 -= x0;
-        y3 -= y0;
-        z3 -= z0;
-
-        double l1 = x1 * x1 + y1 * y1 + z1 * z1;
-        double l2 = x2 * x2 + y2 * y2 + z2 * z2;
-        double l3 = x3 * x3 + y3 * y3 + z3 * z3;
-        double dx = l1 * (y2 * z3 - z2 * y3) - l2 * (y1 * z3 - z1 * y3) + l3 * (y1 * z2 - z1 * y2);
-        double dy = -l1 * (x2 * z3 - z2 * x3) - l2 * (x1 * z3 - z1 * x3) + l3 * (x1 * z2 - z1 * x2);
-        double dz = l1 * (x2 * y3 - y2 * x3) - l2 * (x1 * y3 - y1 * x3) + l3 * (x1 * y2 - y1 * x2);
-
-        double aa = 2 * (x1 * (y2 * z3 - z2 * y3) - x2 * (y1 * z3 - z1 * y3) + x3 * (y1 * z2 - z1 * y2));
-        double x = dx / aa;
-        double y = dy / aa;
-        double z = dz / aa;
-        return Py_BuildValue("(fff)f", x + x0, y + y0, z + z0, sqrt(x * x + y * y + z * z));
+    } else if (PyTuple_Check(points)) {
+        Py_ssize_t numElements = PyTuple_Size(points);
+        if (numElements < 4) {
+            PyErr_SetString(PyExc_ValueError, "fast_3d_circumcircle requires a tuple of length 4 or greater.");
+            return NULL;
+        }
+        if (!get_triple_tuple_elems(PyTuple_GetItem(points, 0), &x0, &y0, &z0))
+            return NULL;
+        if (!get_triple_tuple_elems(PyTuple_GetItem(points, 1), &x1, &y1, &z1))
+            return NULL;
+        if (!get_triple_tuple_elems(PyTuple_GetItem(points, 2), &x2, &y2, &z2))
+            return NULL;
+        if (!get_triple_tuple_elems(PyTuple_GetItem(points, 3), &x3, &y3, &z3))
+            return NULL;
+    } else if (PyArray_Check(points)) {
+        int dims = PyArray_NDIM(points);
+        if (dims != 2) {
+            PyErr_SetString(PyExc_ValueError, "fast_3d_circumcircle requires a two dimensional numpy array.");
+            return NULL;
+        }
+        npy_intp * shape = PyArray_DIMS(points);
+        if ((shape[0] < 4) && (shape[1] != 3)) {
+            PyErr_SetString(PyExc_ValueError, "fast_3d_circumcircle requires a numpy array of width 4, height 3.");
+            return NULL;
+        }
+        double* values = (double *) PyArray_DATA(points);
+        x0 = values[0];
+        y0 = values[1];
+        z0 = values[2];
+        x1 = values[3];
+        y1 = values[4];
+        z1 = values[5];
+        x2 = values[6];
+        y2 = values[7];
+        z2 = values[8];
+        x3 = values[9];
+        y3 = values[10];
+        z3 = values[11];
+
+    } else {
+        PyErr_SetString(PyExc_ValueError, "Points must be a list, tuple, or numpy array.");
+        return NULL;
     }
-    return NULL;
+
+    x1 -= x0;
+    y1 -= y0;
+    z1 -= z0;
+
+    x2 -= x0;
+    y2 -= y0;
+    z2 -= z0;
+
+    x3 -= x0;
+    y3 -= y0;
+    z3 -= z0;
+
+    double l1 = x1 * x1 + y1 * y1 + z1 * z1;
+    double l2 = x2 * x2 + y2 * y2 + z2 * z2;
+    double l3 = x3 * x3 + y3 * y3 + z3 * z3;
+    double dx = l1 * (y2 * z3 - z2 * y3) - l2 * (y1 * z3 - z1 * y3) + l3 * (y1 * z2 - z1 * y2);
+    double dy = -l1 * (x2 * z3 - z2 * x3) - l2 * (x1 * z3 - z1 * x3) + l3 * (x1 * z2 - z1 * x2);
+    double dz = l1 * (x2 * y3 - y2 * x3) - l2 * (x1 * y3 - y1 * x3) + l3 * (x1 * y2 - y1 * x2);
+
+    double aa = 2 * (x1 * (y2 * z3 - z2 * y3) - x2 * (y1 * z3 - z1 * y3) + x3 * (y1 * z2 - z1 * y2));
+    double x = dx / aa;
+    double y = dy / aa;
+    double z = dz / aa;
+    return Py_BuildValue("(fff)f", x + x0, y + y0, z + z0, sqrt(x * x + y * y + z * z));
 }
 
 
 static PyMethodDef triangulation_functions[] = {
     {"fast_norm", fast_norm, METH_VARARGS,
-		"Returns the norm of the given array."},
+		"fast_norm(vec)"
+		"Returns the norm of the given array. Requires one dimensional tuple, list, or numpy array."
+		"For large matrices, it is better to use numpy's linear algebra functions. Additionally, if the values in the"
+		"list are particularly small, squaring them may cause them to become zero. If they're very large, squaring them"
+		"may result in numerical overflow (in the case where they're above 10^150). Can not handle complex numbers."},
     {"fast_2d_circumcircle", fast_2d_circumcircle, METH_VARARGS,
-		"Returns the norm of the given array."},
+     		"fast_2d_circumcircle(list)"
+		"Returns center and radius of the circle touching the first three points in the list. Requires a list, tuple,"
+		 "or numpy array that is 3x2 in shape."},
     {"fast_3d_circumcircle", fast_3d_circumcircle, METH_VARARGS,
-		"Returns the norm of the given array."},
+		"fast_3d_circumcircle(list)"
+		"Returns center and radius of the sphere touching the first four points in the list. Requires a list, tuple,"
+		 "or numpy array that is 4x3 in shape."},
     {"fast_2d_point_in_simplex", (PyCFunction) fast_2d_point_in_simplex, METH_VARARGS | METH_KEYWORDS,
-		"Refer to docs."},
+		"fast_2d_point_in_simplex(point, simplex, eps=1e-8)"
+		"Returns true if the given 2d point is in the simplex, minus some error eps."},
 	{NULL}
 };
 

From 107ad4dc17ba05f0bb37595a0ae6a2ebc0c92722 Mon Sep 17 00:00:00 2001
From: philippeitis <33013301+philippeitis@users.noreply.github.com>
Date: Sun, 15 Dec 2019 14:10:47 -0800
Subject: [PATCH 09/10] Update triangulation.c

Add newlines to docstring. Remove null checks that are unneeded after adding size checks, and add checks for successful PyFloat conversion.
---
 adaptive/learner/triangulation.c | 102 ++++++++++++++++---------------
 1 file changed, 52 insertions(+), 50 deletions(-)

diff --git a/adaptive/learner/triangulation.c b/adaptive/learner/triangulation.c
index 1c81c5831..7681b2f33 100644
--- a/adaptive/learner/triangulation.c
+++ b/adaptive/learner/triangulation.c
@@ -4,6 +4,46 @@
 #include <stdio.h>
 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
 
+static int
+get_tuple_elems(PyObject * tuple, double * first, double * second) {
+    if (PyTuple_Size(tuple) < 2) {
+        PyErr_SetString(PyExc_ValueError, "Tuple must be len 2 or more.");
+        return 0;
+    }
+
+    PyObject * elem = PyTuple_GetItem(tuple, 0);
+    *first = PyFloat_AsDouble(elem);
+    if (first == NULL)
+        return 0;
+    elem = PyTuple_GetItem(tuple, 1);
+    *second = PyFloat_AsDouble(elem);
+    if (second == NULL)
+        return 0;
+    return 1;
+}
+
+static int
+get_triple_tuple_elems(PyObject * tuple, double * first, double * second, double* third) {
+    if (PyTuple_Size(tuple) < 3) {
+        PyErr_SetString(PyExc_ValueError, "Tuple must be len 3 or more.");
+        return 0;
+    }
+
+    PyObject * elem = PyTuple_GetItem(tuple, 0);
+    *first = PyFloat_AsDouble(elem);
+    if (first == NULL)
+        return 0;
+    elem = PyTuple_GetItem(tuple, 1);
+    *second = PyFloat_AsDouble(elem);
+    if (second == NULL)
+        return 0;
+    elem = PyTuple_GetItem(tuple, 2);
+    *third = PyFloat_AsDouble(elem);
+    if (third == NULL)
+        return 0;
+    return 1;
+}
+
 static PyObject *
 fast_norm(PyObject *self, PyObject *args)
 {
@@ -65,41 +105,12 @@ fast_norm(PyObject *self, PyObject *args)
 
 }
 
-static int get_tuple_elems(PyObject * tuple, double * first, double * second) {
-    PyObject * elem = PyTuple_GetItem(tuple, 0);
-    if (elem == NULL)
-        return 0;
-    *first = PyFloat_AsDouble(elem);
-    elem = PyTuple_GetItem(tuple, 1);
-    if (elem == NULL)
-        return 0;
-    *second = PyFloat_AsDouble(elem);
-    return 1;
-}
-
-static int get_triple_tuple_elems(PyObject * tuple, double * first, double * second, double* third) {
-    PyObject * elem = PyTuple_GetItem(tuple, 0);
-    if (elem == NULL)
-        return 0;
-    *first = PyFloat_AsDouble(elem);
-    elem = PyTuple_GetItem(tuple, 1);
-    if (elem == NULL)
-        return 0;
-    *second = PyFloat_AsDouble(elem);
-    elem = PyTuple_GetItem(tuple, 2);
-    if (elem == NULL)
-        return 0;
-    *third = PyFloat_AsDouble(elem);
-    return 1;
-}
-
 static PyObject *
 fast_2d_point_in_simplex(PyObject *self, PyObject *args, PyObject *kwargs)
 {
-    PyObject* point, * simplex;
+    PyObject * point, * simplex;
     double eps = 0.00000001;
 
-    // SystemError: bad format char passed to Py_BuildValue ??
 	static char *kwlist[] = {"point", "simplex", "eps", NULL};
 
 	if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|d", kwlist, &point, &simplex, &eps))
@@ -126,23 +137,14 @@ fast_2d_point_in_simplex(PyObject *self, PyObject *args, PyObject *kwargs)
 
 
     point = PyList_GetItem(simplex, 0);
-    if (point == NULL)
-        return NULL;
-
     if (!get_tuple_elems(point, &p0x, &p0y))
         return NULL;
 
     point = PyList_GetItem(simplex, 1);
-    if (point == NULL)
-        return NULL;
-
     if (!get_tuple_elems(point, &p1x, &p1y))
         return NULL;
 
     point = PyList_GetItem(simplex, 2);
-    if (point == NULL)
-        return NULL;
-
     if (!get_tuple_elems(point, &p2x, &p2y))
         return NULL;
 
@@ -267,15 +269,19 @@ fast_3d_circumcircle(PyObject *self, PyObject *args)
             return NULL;
     } else if (PyArray_Check(points)) {
         int dims = PyArray_NDIM(points);
+
         if (dims != 2) {
             PyErr_SetString(PyExc_ValueError, "fast_3d_circumcircle requires a two dimensional numpy array.");
             return NULL;
         }
+
         npy_intp * shape = PyArray_DIMS(points);
+
         if ((shape[0] < 4) && (shape[1] != 3)) {
             PyErr_SetString(PyExc_ValueError, "fast_3d_circumcircle requires a numpy array of width 4, height 3.");
             return NULL;
         }
+
         double* values = (double *) PyArray_DATA(points);
         x0 = values[0];
         y0 = values[1];
@@ -289,7 +295,6 @@ fast_3d_circumcircle(PyObject *self, PyObject *args)
         x3 = values[9];
         y3 = values[10];
         z3 = values[11];
-
     } else {
         PyErr_SetString(PyExc_ValueError, "Points must be a list, tuple, or numpy array.");
         return NULL;
@@ -318,27 +323,24 @@ fast_3d_circumcircle(PyObject *self, PyObject *args)
     double x = dx / aa;
     double y = dy / aa;
     double z = dz / aa;
+
     return Py_BuildValue("(fff)f", x + x0, y + y0, z + z0, sqrt(x * x + y * y + z * z));
 }
 
 
 static PyMethodDef triangulation_functions[] = {
     {"fast_norm", fast_norm, METH_VARARGS,
-		"fast_norm(vec)"
-		"Returns the norm of the given array. Requires one dimensional tuple, list, or numpy array."
-		"For large matrices, it is better to use numpy's linear algebra functions. Additionally, if the values in the"
-		"list are particularly small, squaring them may cause them to become zero. If they're very large, squaring them"
-		"may result in numerical overflow (in the case where they're above 10^150). Can not handle complex numbers."},
+		"Returns the norm of the given array. Requires one dimensional tuple, list, or numpy array.\n"
+		"For large matrices, it is better to use numpy's linear algebra functions.\nAdditionally, if the values in the "
+		"list are particularly small, squaring them may cause them to become zero.\nIf they're very large, squaring them "
+		"may result in numerical overflow (in the case where they're above 10^150).\nCan not handle complex numbers."},
     {"fast_2d_circumcircle", fast_2d_circumcircle, METH_VARARGS,
-     		"fast_2d_circumcircle(list)"
-		"Returns center and radius of the circle touching the first three points in the list. Requires a list, tuple,"
+		"Returns center and radius of the circle touching the first three points in the list.\nRequires a list, tuple,"
 		 "or numpy array that is 3x2 in shape."},
     {"fast_3d_circumcircle", fast_3d_circumcircle, METH_VARARGS,
-		"fast_3d_circumcircle(list)"
-		"Returns center and radius of the sphere touching the first four points in the list. Requires a list, tuple,"
+		"Returns center and radius of the sphere touching the first four points in the list.\nRequires a list, tuple,"
 		 "or numpy array that is 4x3 in shape."},
     {"fast_2d_point_in_simplex", (PyCFunction) fast_2d_point_in_simplex, METH_VARARGS | METH_KEYWORDS,
-		"fast_2d_point_in_simplex(point, simplex, eps=1e-8)"
 		"Returns true if the given 2d point is in the simplex, minus some error eps."},
 	{NULL}
 };

From 0cfdd030e7edadc3e450c74f40e94882a820a819 Mon Sep 17 00:00:00 2001
From: philippeitis <33013301+philippeitis@users.noreply.github.com>
Date: Mon, 16 Dec 2019 12:09:41 -0800
Subject: [PATCH 10/10] Fix typo

Accidentally moved the minus around.
---
 adaptive/learner/triangulation.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/adaptive/learner/triangulation.c b/adaptive/learner/triangulation.c
index 7681b2f33..6927fba3b 100644
--- a/adaptive/learner/triangulation.c
+++ b/adaptive/learner/triangulation.c
@@ -316,12 +316,12 @@ fast_3d_circumcircle(PyObject *self, PyObject *args)
     double l2 = x2 * x2 + y2 * y2 + z2 * z2;
     double l3 = x3 * x3 + y3 * y3 + z3 * z3;
     double dx = l1 * (y2 * z3 - z2 * y3) - l2 * (y1 * z3 - z1 * y3) + l3 * (y1 * z2 - z1 * y2);
-    double dy = -l1 * (x2 * z3 - z2 * x3) - l2 * (x1 * z3 - z1 * x3) + l3 * (x1 * z2 - z1 * x2);
+    double dy = l1 * (x2 * z3 - z2 * x3) - l2 * (x1 * z3 - z1 * x3) + l3 * (x1 * z2 - z1 * x2);
     double dz = l1 * (x2 * y3 - y2 * x3) - l2 * (x1 * y3 - y1 * x3) + l3 * (x1 * y2 - y1 * x2);
 
     double aa = 2 * (x1 * (y2 * z3 - z2 * y3) - x2 * (y1 * z3 - z1 * y3) + x3 * (y1 * z2 - z1 * y2));
     double x = dx / aa;
-    double y = dy / aa;
+    double y = -dy / aa;
     double z = dz / aa;
 
     return Py_BuildValue("(fff)f", x + x0, y + y0, z + z0, sqrt(x * x + y * y + z * z));