diff --git a/.githooks/pre-commit b/.githooks/pre-commit
new file mode 100755
index 0000000..2ee06ca
--- /dev/null
+++ b/.githooks/pre-commit
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+
+SCR=$(cat << EOF
+import re
+import subprocess
+
+repo = subprocess.check_output("git rev-parse --show-toplevel", shell=True).decode('utf8').strip()
+with open(f"{repo}/README.md", "r") as f:
+  md = f.read()
+with open(f"{repo}/setup.py", "r") as f:
+  content = f.read()
+with open(f"{repo}/setup.py", "w") as f:
+  f.write(re.sub(r"long_description=.*?    url", 'long_description="""' + md + '""",\n    url', content, flags=re.DOTALL))
+subprocess.check_output(f"git add {repo}/setup.py", shell=True)
+EOF
+)
+
+python3 -c "$SCR"
\ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5b236a3..79c4e65 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,19 +1,12 @@
 cmake_minimum_required(VERSION 3.11.0)
 project(elkai)
 
-# set(CMAKE_POSITION_INDEPENDENT_CODE ON)
-
 find_package(PythonExtensions REQUIRED)
 
 include_directories(${CMAKE_SOURCE_DIR}/LKH-3.0.8/SRC/INCLUDE)
 file(GLOB LKH_SRC "LKH-3.0.8/SRC/*.c")
 
-# add_library(LKH ${LKH_SRC})
-# target_compile_definitions(LKH PRIVATE TWO_LEVEL_TREE)
-# target_compile_options(LKH PRIVATE -flto -fcommon)
-
 add_library(_elkai MODULE elkai/_elkai.c ${LKH_SRC})
 python_extension_module(_elkai)
-# target_link_libraries(_elkai LKH)
 
 install(TARGETS _elkai LIBRARY DESTINATION elkai)
diff --git a/LKH-3.0.8/SRC/INCLUDE/LKH.h b/LKH-3.0.8/SRC/INCLUDE/LKH.h
index 793eb04..732e343 100755
--- a/LKH-3.0.8/SRC/INCLUDE/LKH.h
+++ b/LKH-3.0.8/SRC/INCLUDE/LKH.h
@@ -1,17 +1,6 @@
 #ifndef _LKH_H
 #define _LKH_H
 
-#define DEFAULT_PARAMETERS "RUNS = 10\nTRACE_LEVEL = 0\nPROBLEM_FILE = :empty:\n"
-#define DEFAULT_PROBLEM "TYPE : ATSP\n\
-DIMENSION : 3\n\
-EDGE_WEIGHT_TYPE: EXPLICIT\n\
-EDGE_WEIGHT_FORMAT: FULL_MATRIX\n\
-EDGE_WEIGHT_SECTION\n\
-0 4 0\n\
-0 0 5\n\
-0 0 0\n\
-"
-
 /*
  * This header is used by almost all functions of the program. It defines 
  * macros and specifies data structures and function prototypes.
diff --git a/LKH-3.0.8/SRC/LKHmain.c b/LKH-3.0.8/SRC/LKHmain.c
index 61d57b0..823c157 100755
--- a/LKH-3.0.8/SRC/LKHmain.c
+++ b/LKH-3.0.8/SRC/LKHmain.c
@@ -23,6 +23,7 @@ static void ReadableTourClear() {
     ReadableTourAllocated = 100;
     if(ReadableTour != 0) {
         free(ReadableTour);
+        ReadableTour = 0;
     }
 }
 
@@ -60,20 +61,160 @@ void TheTour(int *Tour, GainType Cost)
     }
 }
 
-extern char *ReadLineBuf;
-
 #define GB_STRING_IMPLEMENTATION
 #include "gb_string.h"
 
+extern void ClearLines();
+extern void WriteLine(gbString str);
+
+void ElkaiSolveProblem(gbString params, gbString problem, int *tourSize, int **tourPtr) {
+    ClearLines();
+    WriteLine(params);
+    ReadParameters();
+    GainType Cost, OldOptimum;
+    double Time, LastTime;
+    StartTime = LastTime = GetTime();
+    MaxMatrixDimension = 20000;
+    MergeWithTour =
+        Recombination == GPX2 ? MergeWithTourGPX2 :
+        Recombination == CLARIST ? MergeWithTourCLARIST :
+                                   MergeWithTourIPT;
+
+    WriteLine(problem);
+
+    ReadProblem();
+    AllocateStructures();
+    CreateCandidateSet();
+    InitializeStatistics();
+
+    if (Norm != 0 || Penalty) {
+        Norm = 9999;
+        BestCost = PLUS_INFINITY;
+        BestPenalty = CurrentPenalty = PLUS_INFINITY;
+    } else {
+        /* The ascent has solved the problem! */
+        Optimum = BestCost = (GainType) LowerBound;
+        UpdateStatistics(Optimum, GetTime() - LastTime);
+        RecordBetterTour();
+        RecordBestTour();
+        CurrentPenalty = PLUS_INFINITY;
+        BestPenalty = CurrentPenalty = Penalty ? Penalty() : 0;
+        TheTour(BestTour, BestCost);
+        Runs = 0;
+    }
+
+    /* Find a specified number (Runs) of local optima */
+
+    for (Run = 1; Run <= Runs; Run++) {
+        LastTime = GetTime();
+        if (LastTime - StartTime >= TotalTimeLimit) {
+            if (TraceLevel >= 1)
+                printff("*** Time limit exceeded ***\n");
+            Run--;
+            break;
+        }
+        Cost = FindTour();      /* using the Lin-Kernighan heuristic */
+        if (MaxPopulationSize > 1 && !TSPTW_Makespan) {
+            /* Genetic algorithm */
+            int i;
+            for (i = 0; i < PopulationSize; i++) {
+                Cost = MergeTourWithIndividual(i);
+            }
+            if (!HasFitness(CurrentPenalty, Cost)) {
+                if (PopulationSize < MaxPopulationSize) {
+                    AddToPopulation(CurrentPenalty, Cost);
+                    if (TraceLevel >= 1)
+                        PrintPopulation();
+                } else if (SmallerFitness(CurrentPenalty, Cost,
+                                          PopulationSize - 1)) {
+                    i = ReplacementIndividual(CurrentPenalty, Cost);
+                    ReplaceIndividualWithTour(i, CurrentPenalty, Cost);
+                    if (TraceLevel >= 1)
+                        PrintPopulation();
+                }
+            }
+        } else if (Run > 1 && !TSPTW_Makespan)
+            Cost = MergeTourWithBestTour();
+        if (CurrentPenalty < BestPenalty ||
+            (CurrentPenalty == BestPenalty && Cost < BestCost)) {
+            BestPenalty = CurrentPenalty;
+            BestCost = Cost;
+            RecordBetterTour();
+            RecordBestTour();
+            TheTour(BestTour, BestCost);
+        }
+        OldOptimum = Optimum;
+        if (!Penalty ||
+            (MTSPObjective != MINMAX && MTSPObjective != MINMAX_SIZE)) {
+            if (CurrentPenalty == 0 && Cost < Optimum)
+                Optimum = Cost;
+        } else if (CurrentPenalty < Optimum)
+            Optimum = CurrentPenalty;
+        if (Optimum < OldOptimum) {
+            if (FirstNode->InputSuc) {
+                Node *N = FirstNode;
+                while ((N = N->InputSuc = N->Suc) != FirstNode);
+            }
+        }
+        Time = fabs(GetTime() - LastTime);
+        UpdateStatistics(Cost, Time);
+        if (StopAtOptimum && MaxPopulationSize >= 1) {
+            if (ProblemType != CCVRP && ProblemType != TRP &&
+                ProblemType != MLP &&
+                MTSPObjective != MINMAX &&
+                MTSPObjective != MINMAX_SIZE ?
+                CurrentPenalty == 0 && Cost == Optimum :
+                CurrentPenalty == Optimum) {
+                Runs = Run;
+                break;
+            }
+        }
+        if (PopulationSize >= 2 &&
+            (PopulationSize == MaxPopulationSize ||
+             Run >= 2 * MaxPopulationSize) && Run < Runs) {
+            Node *N;
+            int Parent1, Parent2;
+            Parent1 = LinearSelection(PopulationSize, 1.25);
+            do
+                Parent2 = LinearSelection(PopulationSize, 1.25);
+            while (Parent2 == Parent1);
+            ApplyCrossover(Parent1, Parent2);
+            N = FirstNode;
+            do {
+                if (ProblemType != HCP && ProblemType != HPP) {
+                    int d = C(N, N->Suc);
+                    AddCandidate(N, N->Suc, d, INT_MAX);
+                    AddCandidate(N->Suc, N, d, INT_MAX);
+                }
+                N = N->InitialSuc = N->Suc;
+            } while (N != FirstNode);
+        }
+        SRandom(++Seed);
+    }
+
+    *tourSize = ReadableTourSize;
+
+    if(tourPtr != 0) {
+        *tourPtr = ReadableTour;
+    }
+
+//    if(tour != 0) {
+//        for(int M = 0; M < ReadableTourSize; M++) {
+//            tour[M] = ReadableTour[M];
+//        }
+//    }
+}
+
 int ElkaiSolveATSP(int dimension, float *weights, int *tour, int runs)
 {
     GainType Cost, OldOptimum;
     double Time, LastTime;
+    ClearLines();
+    WriteLine(gb_make_string("TRACE_LEVEL = 0\nPROBLEM_FILE = :stdin:\n"));
 
-    ReadLineBuf = gb_make_string("TRACE_LEVEL = 0\nPROBLEM_FILE = :empty:\n");
     char runsStr[100];
     sprintf(runsStr, "RUNS = %d\n", runs);
-    ReadLineBuf = gb_append_cstring(ReadLineBuf, runsStr);
+    WriteLine(gb_make_string(runsStr));
 
     ReadParameters();
 
@@ -83,42 +224,23 @@ int ElkaiSolveATSP(int dimension, float *weights, int *tour, int runs)
         Recombination == GPX2 ? MergeWithTourGPX2 :
         Recombination == CLARIST ? MergeWithTourCLARIST :
                                    MergeWithTourIPT;
-
-    ReadLineBuf = (char *) malloc(sizeof(char) * 1000);
     
     char dimensionStr[100];
     sprintf(dimensionStr, "DIMENSION : %d\n", dimension);
-    ReadLineBuf = gb_make_string("TYPE: ATSP\n");
-    ReadLineBuf = gb_append_cstring(ReadLineBuf, dimensionStr);
-    ReadLineBuf = gb_append_cstring(ReadLineBuf, "EDGE_WEIGHT_TYPE: EXPLICIT\nEDGE_WEIGHT_FORMAT: FULL_MATRIX\nEDGE_WEIGHT_SECTION\n");
+    WriteLine(gb_make_string("TYPE: ATSP\n"));
+    WriteLine(gb_make_string(dimensionStr));
+    WriteLine(gb_make_string("EDGE_WEIGHT_TYPE: EXPLICIT\nEDGE_WEIGHT_FORMAT: FULL_MATRIX\nEDGE_WEIGHT_SECTION\n"));
 
     for(int y = 0; y < dimension; y++) {
         for(int x = 0; x < dimension; x++) {
             char weightStr[64];
             sprintf(weightStr, "%f ", weights[y * dimension + x]);
-            ReadLineBuf = gb_append_cstring(ReadLineBuf, weightStr);
+            WriteLine(gb_make_string(weightStr));
         }
-        ReadLineBuf = gb_append_cstring(ReadLineBuf, "\n");
+        WriteLine(gb_make_string("\n"));
     }
 
     ReadProblem();
-    if (SubproblemSize > 0) {
-        if (DelaunayPartitioning)
-            SolveDelaunaySubproblems();
-        else if (KarpPartitioning)
-            SolveKarpSubproblems();
-        else if (KCenterPartitioning)
-            SolveKCenterSubproblems();
-        else if (KMeansPartitioning)
-            SolveKMeansSubproblems();
-        else if (RohePartitioning)
-            SolveRoheSubproblems();
-        else if (MoorePartitioning || SierpinskiPartitioning)
-            SolveSFCSubproblems();
-        else
-            SolveTourSegmentSubproblems();
-        return EXIT_SUCCESS;
-    }
     AllocateStructures();
     CreateCandidateSet();
     InitializeStatistics();
diff --git a/LKH-3.0.8/SRC/ReadLine.c b/LKH-3.0.8/SRC/ReadLine.c
index e04cfa0..7c64911 100755
--- a/LKH-3.0.8/SRC/ReadLine.c
+++ b/LKH-3.0.8/SRC/ReadLine.c
@@ -20,33 +20,62 @@ static int EndOfLine(FILE * InputFile, int c)
     return EOL;
 }
 
-char *ReadLineBuf = ":empty:";
+char *ReadLineBuf = 0;
+int ReadLinePtr = 0;
 #include "gb_string.h"
 
+void WriteLine(gbString str) {
+    if(ReadLineBuf == 0) {
+        ReadLineBuf = gb_make_string("");
+    }
+    ReadLineBuf = gb_append_string(ReadLineBuf, str);
+}
+
+void ClearLines() {
+    ReadLinePtr = 0;
+    if(ReadLineBuf != 0) {
+        gb_free_string(ReadLineBuf);
+        ReadLineBuf = 0;
+    }
+}
+
+double ReadNumber() {
+    if(ReadLinePtr == 0) return 0;
+    char *k = ReadLineBuf + ReadLinePtr;
+    double output = strtof(ReadLineBuf + ReadLinePtr, &k);
+    ReadLinePtr += k - (ReadLineBuf + ReadLinePtr);
+    return output;
+}
+
 char *ReadLine(FILE * InputFile)
 {
     if(InputFile == 0) {
-        if(ReadLineBuf[0] == '\0') {
+        if(ReadLineBuf[ReadLinePtr] == '\0') {
             return 0;
         }
 
         gbString currentLine = gb_make_string("");
 
-        while(ReadLineBuf[0] != '\0') {
+        while(ReadLineBuf[ReadLinePtr] != '\0') {
             char singleCh[2];
-            singleCh[0] = ReadLineBuf[0];
+            singleCh[0] = ReadLineBuf[ReadLinePtr];
             singleCh[1] = '\0';
 
             currentLine = gb_append_cstring(currentLine, singleCh);
 
-            ReadLineBuf++;
-            if(ReadLineBuf[0] == '\n') {
-                ReadLineBuf++;
+            ReadLinePtr++;
+            if(ReadLineBuf[ReadLinePtr] == '\n') {
+                ReadLinePtr++;
                 break;
             }
         }
 
-        return (char*)currentLine;
+        gbUsize lineSize = gb_string_length(currentLine);
+        char *L = malloc(lineSize + 1);
+        memcpy(L, currentLine, lineSize + 1);
+        gb_free_string(currentLine);
+
+        return L;
     }
 
     int i, c;
diff --git a/LKH-3.0.8/SRC/ReadParameters.c b/LKH-3.0.8/SRC/ReadParameters.c
index c81d9f0..18dec48 100755
--- a/LKH-3.0.8/SRC/ReadParameters.c
+++ b/LKH-3.0.8/SRC/ReadParameters.c
@@ -485,8 +485,6 @@ static char *ReadYesOrNo(int *V);
 #undef max
 static size_t max(size_t a, size_t b);
 
-extern char *ReadLineBuf;
-
 void _Reset1();
 void _Reset2();
 void _Reset3();
@@ -612,9 +610,6 @@ void ReadParameters()
         }
     } else {
         ParameterFile = 0;
-        if(!strcmp(ReadLineBuf, ":empty:")) {
-            ReadLineBuf = DEFAULT_PARAMETERS;
-        }
     }
     while ((Line = ReadLine(ParameterFile))) {
         if (!(Keyword = strtok(Line, Delimiters)))
@@ -1202,9 +1197,8 @@ void ReadParameters()
         DelaunayPure = 0;
     if(ParameterFile != 0)
     fclose(ParameterFile);
-    free(LastLine);
+//    free(LastLine);
     LastLine = 0;
-    ReadLineBuf = ":empty:";
 }
 
 static char *GetFileName(char *Line)
diff --git a/LKH-3.0.8/SRC/ReadProblem.c b/LKH-3.0.8/SRC/ReadProblem.c
index ba9f486..8a3c7d4 100755
--- a/LKH-3.0.8/SRC/ReadProblem.c
+++ b/LKH-3.0.8/SRC/ReadProblem.c
@@ -373,18 +373,15 @@ static int TwoDWeightType(void);
 static int ThreeDWeightType(void);
 static void Convert2FullMatrix(void);
 
-extern char *ReadLineBuf;
+extern double ReadNumber();
 
 void ReadProblem()
 {
     int i, j, K;
     char *Line, *Keyword;
 
-    if(!strcmp(ProblemFileName, ":empty:")) {
+    if(ProblemFileName == 0 || !strcmp(ProblemFileName, ":stdin:")) {
         ProblemFile = 0;
-        if(!strcmp(ReadLineBuf, ":empty:")) {
-            ReadLineBuf = DEFAULT_PROBLEM;
-        }
     } else {
         if (!(ProblemFile = fopen(ProblemFileName, "r")))
             eprintf("Cannot open PROBLEM_FILE: \"%s\"", ProblemFileName);
@@ -800,7 +797,7 @@ void ReadProblem()
         for (i = 0; i < MergeTourFiles; i++)
             ReadTour(MergeTourFileName[i], &MergeTourFile[i]);
     }
-    free(LastLine);
+//    free(LastLine);
     LastLine = 0;
 }
 
@@ -1341,6 +1338,7 @@ static void Read_EDGE_WEIGHT_SECTION()
     if(WeightFormat != FULL_MATRIX) {
         eprintf("EDGE_WEIGHT_SECTION: Weight format not FULL_MATRIX");
     }
+
     switch (WeightFormat) {
     case FULL_MATRIX:
         for (i = 1; i <= Dim; i++) {
@@ -1350,9 +1348,7 @@ static void Read_EDGE_WEIGHT_SECTION()
                 if (!fscanf(ProblemFile, "%lf", &w))
                     eprintf("EDGE_WEIGHT_SECTION: Missing weight");
                 } else {
-                    w = strtof(ReadLineBuf, &ReadLineBuf);
-                    // if(!sscanf(ReadLineBuf, "%lf", &w))
-                    // eprintf("EDGE_WEIGHT_SECTION: Missing weight");
+                    w = ReadNumber();
                 }
                 W = round(Scale * w);
                 if (Asymmetric) {
diff --git a/LKH-3.0.8/SRC/eprintf.c b/LKH-3.0.8/SRC/eprintf.c
index 3a9552b..9433a8a 100755
--- a/LKH-3.0.8/SRC/eprintf.c
+++ b/LKH-3.0.8/SRC/eprintf.c
@@ -1,5 +1,6 @@
 #include "LKH.h"
 #include <stdarg.h>
+#include <Python.h>
 
 /* 
  * The eprintf function prints an error message and exits.
@@ -9,12 +10,20 @@ void eprintf(const char *fmt, ...)
 {
     va_list args;
 
-    if (LastLine && *LastLine)
-        fprintf(stderr, "\n%s\n", LastLine);
-    fprintf(stderr, "\n*** Error ***\n");
+    char err[1024];
     va_start(args, fmt);
-    vfprintf(stderr, fmt, args);
+    int res = vsprintf(err, fmt, args);
     va_end(args);
-    fprintf(stderr, "\n");
-    exit(EXIT_FAILURE);
+
+    PyErr_SetString(PyExc_TypeError, err);
+//    exit(0);
+
+//    if (LastLine && *LastLine)
+//        fprintf(stderr, "\n%s\n", LastLine);
+//    fprintf(stderr, "\n*** Error ***\n");
+//    va_start(args, fmt);
+//    vfprintf(stderr, fmt, args);
+//    va_end(args);
+//    fprintf(stderr, "\n");
+//    exit(EXIT_FAILURE);
 }
diff --git a/README-wip.md b/README-wip.md
new file mode 100644
index 0000000..72b0a94
--- /dev/null
+++ b/README-wip.md
@@ -0,0 +1,58 @@
+# elkai - a Python library for solving TSP problems
+
+[elkai](https://pypi.org/project/elkai/) is a Python 3 library for solving [travelling salesman problems](https://en.wikipedia.org/wiki/Travelling_salesman_problem):
+
+* ⚡ running fast native code with prebuilt wheels for most platforms
+* 🗺️ based on [LKH](http://akira.ruc.dk/~keld/research/LKH/) by Keld Helsgaun, with proven optimal solutions up to N=315
+* 🛣️ supports asymmetric distances (ATSP)
+
+[![Python build](https://github.com/fikisipi/elkai/actions/workflows/python-app.yml/badge.svg)](https://github.com/fikisipi/elkai/actions/workflows/python-app.yml)
+[![image](https://img.shields.io/pypi/v/elkai.svg)](https://pypi.org/project/elkai/)
+
+## Example usage 
+
+```python
+import elkai
+
+cities = elkai.DistanceMatrix([
+    [0, 4, 0],
+    [0, 0, 5],
+    [0, 0, 0]
+])
+
+print(cities.atsp())
+```
+
+```python
+import elkai
+
+cities = elkai.PositionMatrix([
+ [20, 20],
+ [30, 30],
+ [5, -1]   
+])
+
+print(cities.atsp())
+```
+
+```mermaid
+graph TD;
+    0-->|4|1;
+    1-->|0|0;
+    0-->|0|2;
+    1-->|5|2;
+    2-->|0|1;
+```
+
+> **Note**
+> solve_int_matrix and solve_float_matrix are deprecated in v1.
+
+## Installation
+
+💾 **To install it** run `pip install elkai`
+
+## Notes
+
+⚠️ elkai takes the **global interpreter lock (GIL)** during the solving phase which means two threads cannot solve problems at the same time. If you want to run other workloads at the same time, you have to run another process - for example by using the `multiprocessing` module.
+
+The LKH native code by Helsgaun is released for non-commercial use only. Therefore the same restriction applies to elkai, which is explained in the `LICENSE` file. If there isn't a prebuilt wheel for your platform, you'll have to follow the `scikit-build` process.
\ No newline at end of file
diff --git a/README.md b/README.md
index c56befd..fbc03ef 100644
--- a/README.md
+++ b/README.md
@@ -1,10 +1,11 @@
-# elkai - A Python solver for the Travelling Salesman Problem / TSP
+#  elkai - a Python library for solving TSP problems
 
 [elkai](https://pypi.org/project/elkai/) is a Python 3 library for solving [travelling salesman problems](https://en.wikipedia.org/wiki/Travelling_salesman_problem):
 
-* ⚡ running fast native code with prebuilt wheels for most platforms
 * 🗺️ based on [LKH](http://akira.ruc.dk/~keld/research/LKH/) by Keld Helsgaun, with proven optimal solutions up to N=315
+* ⚡ running fast native code with prebuilt wheels for most platforms
 * 🛣️ supports asymmetric distances (ATSP)
+* cleaner API and more accurate results than [Google's OR tools](https://developers.google.com/optimization/routing/tsp)
 
 [![Python build](https://github.com/fikisipi/elkai/actions/workflows/python-app.yml/badge.svg)](https://github.com/fikisipi/elkai/actions/workflows/python-app.yml)
 [![image](https://img.shields.io/pypi/v/elkai.svg)](https://pypi.org/project/elkai/)
@@ -12,17 +13,18 @@
 ## Example usage 
 
 ```python
-import numpy as np
 import elkai
 
-M = np.zeros((3, 3), dtype=int)
-M[0, 1] = 4
-M[1, 2] = 5
+distance_matrix = [
+    [0, 4, 0],
+    [0, 0, 5],
+    [0, 0, 0]
+]
 
-solution = elkai.solve_int_matrix(M)
+cities = elkai.solve(distance_matrix, skip_end=False)
 
-print(solution)
-# Output: [0, 2, 1]
+print(cities)
+# Output: [0, 2, 1, 0]
 ```
 
 ```mermaid
@@ -34,41 +36,15 @@ graph TD;
     2-->|0|1;
 ```
 
+> **Note**
+> solve_int_matrix and solve_float_matrix are deprecated in v1.
+
 ## Installation
 
 💾 **To install it** run `pip install elkai`
 
-Documentation
--------------
-
-
-**elkai.solve_int_matrix(matrix: List[List[int]], runs=10) -> List**
-
-* `matrix` is a list of lists or **2D numpy array** containing the distances between cities
-* `runs` is the solver iteration count
-
-An example matrix with 3 cities would be:
-
-```python
-[                 # cities are zero indexed, d() is distance
-   [0, 4,  3],    # d(0, 0), d(0, 1), d(0, 2)
-   [4, 0, 10],    # d(1, 0), d(1, 1), ...
-   [2, 4,  0]     # ... and so on
-]
-```
-
-So, the output would be `[0, 2, 1]` because it's best to visit `0 => 2 => 1 => 0`.
-
-⚠️ The final return to the start is implied and **is NOT** part of the output list, i.e. `len(output) == N`
-
-----
-
-**elkai.solve_float_matrix(matrix: List[List[float]], runs=10) -> List**
-
-Same behaviour as above, with float distances supported. Note that there may be precision issues.
-
 ## Notes
 
 ⚠️ elkai takes the **global interpreter lock (GIL)** during the solving phase which means two threads cannot solve problems at the same time. If you want to run other workloads at the same time, you have to run another process - for example by using the `multiprocessing` module.
 
-The LKH native code by Helsgaun is released for non-commercial use only. Therefore the same restriction applies to elkai, which is explained in the `LICENSE` file. If there isn't a prebuilt wheel for your platform, you'll have to follow the `scikit-build` process.
+The LKH native code by Helsgaun is released for non-commercial use only. Therefore the same restriction applies to elkai, which is explained in the `LICENSE` file. If there isn't a prebuilt wheel for your platform, you'll have to follow the `scikit-build` process.
\ No newline at end of file
diff --git a/build_wheels.sh b/build_wheels.sh
deleted file mode 100644
index 551bcfc..0000000
--- a/build_wheels.sh
+++ /dev/null
@@ -1,23 +0,0 @@
-#!/bin/bash
-set -e -x
-
-# Install a system package required by our library
-yum install -y atlas-devel
-
-MPATH="$PATH"
-
-# Compile wheels
-for PYBIN in /opt/python/cp3*/bin; do
-    PATH="${PYBIN}/:${MPATH}"
-    "${PYBIN}/python" -m pip install cmake
-    "${PYBIN}/python" -m pip install scikit-build ninja
-    "${PYBIN}/python" -m pip wheel /io/ -w dist/
-    /bin/rm -rf /io/_skbuild
-    /bin/rm -rf /io/CMakeFiles || true
-done
-
-# Bundle external shared libraries into the wheels
-PLAT="manylinux1_x86_64"
-for whl in dist/*.whl; do
-    auditwheel repair "$whl" --plat $PLAT -w /io/dist/
-done
diff --git a/elkai/__init__.py b/elkai/__init__.py
index 503acfb..e60ac05 100644
--- a/elkai/__init__.py
+++ b/elkai/__init__.py
@@ -1,41 +1,29 @@
-from ._elkai import solve
-from . import utils
-import struct
-
-
-def float_to_int(x):
-    """Converts a float into some unknown integer such that
-    the order is preserved: x <= y <=> f(x) <= f(y)"""
-    x_sign = int(x < 0)
-
-    int_repr = struct.unpack("<i", struct.pack("<f", abs(x / 100)))[0]
-    SQUEEZE_FACTOR = 0.7
-    int_sq = int(abs(int_repr) ** SQUEEZE_FACTOR) * ((-1) ** x_sign)
-    return int_sq
-
-
-def flatten_matrix(mat, mat_type=int):
-    MAXINT = 2147483647
-    N = len(mat)
-    if N < 2:
-        raise TypeError("Argument must be a N*N matrix with N > 1")
-    flattened = []
-    for row in mat:
-        if len(row) != N:
-            raise TypeError("Argument must be a N*N matrix")
-        for column in row:
-            int_col = int(column)
-            if int_col > MAXINT:
-                raise TypeError("Distance must be < 2^31. Try the float version")
-            flattened.append(mat_type(column))
-    return flattened
-
-
-def solve_int_matrix(matrix, runs=10):
-    flat_mat = flatten_matrix(matrix)
-    return solve(flat_mat, runs)
-
-
-def solve_float_matrix(matrix, runs=10):
-    flat_mat = [float_to_int(float(x)) for x in flatten_matrix(matrix, mat_type=float)]
-    return solve(flat_mat, runs)
+from . import utils, _elkai
+from typing import List
+
+
+def solve(matrix: List[List[float]], runs=10, skip_end=True):
+    """Given a 2D matrix of distances between cities, returns the indices of the cities ordered optimally for TSP."""
+
+    flattened = utils.flatten_matrix(matrix)
+    cities: List[int] = _elkai.solve(flattened, runs)
+    if not skip_end:
+        # include the return to first city
+        cities.append(cities[0])
+    return cities
+
+
+def solve_int_matrix(matrix: List[List[int]], runs=10, skip_end=True):
+    """[Deprecated] Given a 2D matrix of distances between cities, returns the indices of the cities ordered optimally for TSP."""
+
+    # Note: this function (and the split between int and float solvers) is redundant.
+    # The C library accepts floats natively now!
+    return solve(matrix, runs, skip_end)
+
+
+def solve_float_matrix(matrix, runs=10, skip_end=True):
+    """[Deprecated] Given a 2D matrix of distances between cities, returns the indices of the cities ordered optimally for TSP."""
+
+    # Note: this function (and the split between int and float solvers) is redundant.
+    # The C library accepts floats natively now!
+    return solve(matrix, runs, skip_end)
\ No newline at end of file
diff --git a/elkai/_elkai.c b/elkai/_elkai.c
index 837b116..b9c79bc 100644
--- a/elkai/_elkai.c
+++ b/elkai/_elkai.c
@@ -1,5 +1,6 @@
 #include <Python.h>
 #include "math.h"
+#include "gb_string.h"
 
 // TODO:
 // 1. Change float matrix handling
@@ -10,6 +11,49 @@
 // 6. Add readthedocs.io and better README / graph images
 
 extern int ElkaiSolveATSP(int dimension, float *weights, int *tour, int runs);
+extern void ElkaiSolveProblem(gbString params, gbString problem, int *tourSize, int *tour);
+
+static PyObject *PySolveProblem(PyObject *self, PyObject *args)
+{
+    if(PyObject_Length(args) != 2) {
+        PyErr_SetString(PyExc_TypeError, "Expected two arguments");
+        return 0;
+    }
+
+    PyObject *arg1 = PyObject_GetItem(args, PyLong_FromLong(0));
+    PyObject *arg2 = PyObject_GetItem(args, PyLong_FromLong(1));
+
+    if(!PyUnicode_Check(arg1) || !PyUnicode_Check(arg2)) {
+        PyErr_SetString(PyExc_TypeError, "Arguments should be strings");
+        return 0;
+    }
+
+    const char *paramsStr = PyUnicode_AsUTF8(arg1);
+    const char *problemStr = PyUnicode_AsUTF8(arg2);
+
+    gbString params = gb_make_string(paramsStr);
+    gbString problem = gb_make_string(problemStr);
+
+    int tourSize = 0;
+    int *tourPtr;
+
+    ElkaiSolveProblem(params, problem, &tourSize, &tourPtr);
+
+    if(PyErr_Occurred() != 0) {
+        return 0;
+    }
+
+    PyObject *list = PyList_New(tourSize);
+    int i = 0;
+    for (i = 0; i < tourSize; i++)
+    {
+        PyObject *tourElement = PyLong_FromLong((long)(tourPtr[i]));
+        PyList_SetItem(list, i, tourElement);
+    }
+    gb_free_string(params);
+    gb_free_string(problem);
+    return list;
+}
 
 static PyObject *ElkSolve(PyObject *self, PyObject *args)
 {
@@ -82,14 +126,17 @@ static char elk_docs[] =
 static PyMethodDef funcs[] = {
     {"solve", (PyCFunction) ElkSolve,
      METH_VARARGS, elk_docs},
-    {NULL}};
+     {"solve_problem", (PyCFunction) PySolveProblem, METH_VARARGS, ""},
+    {NULL}
+};
 
 static struct PyModuleDef elkDef = {
     PyModuleDef_HEAD_INIT,
     "_elkai",
     "",
     -1,
-    funcs};
+    funcs
+};
 
 PyMODINIT_FUNC PyInit__elkai(void)
 {
diff --git a/elkai/utils.py b/elkai/utils.py
index 62eb159..8c96f1c 100644
--- a/elkai/utils.py
+++ b/elkai/utils.py
@@ -1,13 +1,27 @@
-def path_distance(path, matrix):
+def path_distance(path, distance_matrix):
+    """Returns the distance of a path given a distance matrix."""
     dist = 0
     last_idx = path[0]
     tour = path[1:]
 
-    if len(path) != len(matrix[0]) + 1:
+    if len(path) != len(distance_matrix[0]) + 1:
         tour.append(path[0])
 
     for idx in tour:
-        dist += matrix[last_idx][idx]
+        dist += distance_matrix[last_idx][idx]
         last_idx = idx
 
     return dist
+
+def flatten_matrix(matrix):
+    """Flattens a 2D matrix into a list."""
+    N = len(matrix)
+    if N <= 2:
+        raise TypeError("Argument must be a N*N matrix with N > 2")
+    flattened = []
+    for row in matrix:
+        if len(row) != N:
+            raise TypeError("Argument must be a N*N matrix")
+        for column in row:
+            flattened.append(column)
+    return flattened
\ No newline at end of file
diff --git a/pytest.ini b/pytest.ini
index 9d59c18..e960201 100644
--- a/pytest.ini
+++ b/pytest.ini
@@ -1,3 +1,4 @@
 [pytest]
 log_cli = 1
-log_cli_level = INFO
\ No newline at end of file
+log_cli_level = INFO
+addopts=-s
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 24ee1ae..3763601 100644
--- a/setup.py
+++ b/setup.py
@@ -1,16 +1,60 @@
 from skbuild import setup
-from os import path
-this_directory = path.abspath(path.dirname(__file__))
-with open(path.join(this_directory, 'README.md'), encoding='utf-8') as f:
-    long_description = f.read()
 
 setup(
     name="elkai",
     description="A Python solver for the travelling salesman problem - TSP based on native/fast LKH",
-    #long_description=long_description,
-    #long_description_content_type='text/markdown',
+    long_description="""#  elkai - a Python library for solving TSP problems
+
+[elkai](https://pypi.org/project/elkai/) is a Python 3 library for solving [travelling salesman problems](https://en.wikipedia.org/wiki/Travelling_salesman_problem):
+
+* 🗺️ based on [LKH](http://akira.ruc.dk/~keld/research/LKH/) by Keld Helsgaun, with proven optimal solutions up to N=315
+* ⚡ running fast native code with prebuilt wheels for most platforms
+* 🛣️ supports asymmetric distances (ATSP)
+* cleaner API and more accurate results than [Google's OR tools](https://developers.google.com/optimization/routing/tsp)
+
+[![Python build](https://github.com/fikisipi/elkai/actions/workflows/python-app.yml/badge.svg)](https://github.com/fikisipi/elkai/actions/workflows/python-app.yml)
+[![image](https://img.shields.io/pypi/v/elkai.svg)](https://pypi.org/project/elkai/)
+
+## Example usage 
+
+```python
+import elkai
+
+distance_matrix = [
+    [0, 4, 0],
+    [0, 0, 5],
+    [0, 0, 0]
+]
+
+cities = elkai.solve(distance_matrix, skip_end=False)
+
+print(cities)
+# Output: [0, 2, 1, 0]
+```
+
+```mermaid
+graph TD;
+    0-->|4|1;
+    1-->|0|0;
+    0-->|0|2;
+    1-->|5|2;
+    2-->|0|1;
+```
+
+> **Note**
+> solve_int_matrix and solve_float_matrix are deprecated in v1.
+
+## Installation
+
+💾 **To install it** run `pip install elkai`
+
+## Notes
+
+⚠️ elkai takes the **global interpreter lock (GIL)** during the solving phase which means two threads cannot solve problems at the same time. If you want to run other workloads at the same time, you have to run another process - for example by using the `multiprocessing` module.
+
+The LKH native code by Helsgaun is released for non-commercial use only. Therefore the same restriction applies to elkai, which is explained in the `LICENSE` file. If there isn't a prebuilt wheel for your platform, you'll have to follow the `scikit-build` process.""",
     url="https://github.com/fikisipi/elkai",
-    version="0.1.3",
+    version="1.0.2",
     packages=['elkai'],
     author="Filip Dimitrovski",
     license="MIT",
diff --git a/tests/base_test.py b/tests/base_test.py
index 0f26b29..b5bb998 100644
--- a/tests/base_test.py
+++ b/tests/base_test.py
@@ -1,15 +1,7 @@
-# import unittest
-#
-# class SimpleTest(unittest.TestCase):
-#     def runTest(self):
-#         import elkai
-#         res = elkai.solve_int_matrix([[0, 0], [0, 0]])
-#         print(res)
-#         self.assertEqual(1, 1)
-
 import sys
 sys.path.pop(0)
 
+import pytest
 import logging
 import elkai
 from .known_solutions import cases
@@ -17,10 +9,36 @@
 logging.basicConfig(level=logging.DEBUG)
 log = logging.getLogger()
 
+def test_invalid_inputs():
+    with pytest.raises(Exception):
+        elkai.solve_int_matrix([])
+
+    with pytest.raises(Exception):
+        elkai.solve_int_matrix([1, 2, 3])
+
+    with pytest.raises(Exception):
+        elkai.solve_int_matrix("")
 
-def test_elkai():
-    for case in cases[1:]:
-        log.info("%s %s %s", case.name, case.input, case.distance)
+    with pytest.raises(Exception):
+        elkai.solve_int_matrix([[1, 2, 3], [1, 2, 3]])
+
+def test_valid():
+    for case in cases:
         solution = elkai.solve_int_matrix(case.input)
+        log.info("%s %s %s %s", case.name, case.input, case.distance, solution)
         dist = elkai.utils.path_distance(solution, case.input)
         assert dist == case.distance
+
+def test_with_return():
+    assert len(elkai.solve_int_matrix(
+        [
+            [0, 4, 3],
+            [4, 0, 10],
+            [2, 4, 0]
+        ], skip_end=False
+    )) == 4
+
+def test_internal():
+    res = elkai._elkai.solve_problem("TRACE_LEVEL = 0\nPROBLEM_FILE = :stdin:\nRUNS = 10\n",
+                   "TYPE : ATSP\nDIMENSION : 3 \nEDGE_WEIGHT_TYPE: EXPLICIT\nEDGE_WEIGHT_FORMAT: FULL_MATRIX\nEDGE_WEIGHT_SECTION\n0 0 0\n0 0 0\n0 0 0\n")
+    assert isinstance(res, list)