diff --git a/ext/bottle.c b/ext/bottle.c
deleted file mode 100644
index aef6512b..00000000
--- a/ext/bottle.c
+++ /dev/null
@@ -1 +0,0 @@
-#include "./ma/ma.c"
diff --git a/ext/bottle.o b/ext/bottle.o
deleted file mode 100644
index a1eb9c87..00000000
Binary files a/ext/bottle.o and /dev/null differ
diff --git a/ext/bottle.so b/ext/bottle.so
deleted file mode 100755
index b4396aad..00000000
Binary files a/ext/bottle.so and /dev/null differ
diff --git a/ext/ma/ma.c b/ext/ma/ma.c
deleted file mode 100644
index 9ec33ea0..00000000
--- a/ext/ma/ma.c
+++ /dev/null
@@ -1,36 +0,0 @@
-#include <stdio.h>
-#include <gsl/gsl_vector.h>
-
-const int MULTIPLICITY = 1;
-
-gsl_vector_int *
-gsl_vector_ma_equal (gsl_vector * u, gsl_vector * v)
-{
-  const size_t n = v->size;
-  const size_t stride_u = u->stride;
-  const size_t stride_v = v->stride;
-
-  gsl_vector_int * out = gsl_vector_int_calloc(n);
-
-  size_t j;
-
-  if (u->size != v->size)
-    {
-      GSL_ERROR_VAL ("vectors must have same length", GSL_EBADLEN, 0);
-    }
-
-  for (j = 0; j < n; j++)
-    {
-      size_t k;
-
-      for (k = 0; k < MULTIPLICITY; k++)
-        {
-          if (u->data[MULTIPLICITY * stride_u * j + k] == v->data[MULTIPLICITY * stride_v * j + k])
-            {
-              out->data[j] = 1;
-            }
-        }
-    }
-
-  return out;
-}
diff --git a/src/bottle.cr b/src/bottle.cr
index 1cb40e60..9ab99492 100644
--- a/src/bottle.cr
+++ b/src/bottle.cr
@@ -1,3 +1,5 @@
+require "./vector/*"
+
 module Bottle
   VERSION = "0.1.0"
 end
diff --git a/src/util/exceptions.cr b/src/core/exceptions.cr
similarity index 71%
rename from src/util/exceptions.cr
rename to src/core/exceptions.cr
index bb9241a4..6b874270 100644
--- a/src/util/exceptions.cr
+++ b/src/core/exceptions.cr
@@ -1,4 +1,4 @@
-module Bottle::Util::Exceptions
+module Bottle::Core::Exceptions
   class RangeError < Exception
   end
 
diff --git a/src/core/matrix/index.cr b/src/core/matrix/index.cr
new file mode 100644
index 00000000..51926d69
--- /dev/null
+++ b/src/core/matrix/index.cr
@@ -0,0 +1,19 @@
+macro matrix_indexing_abstract(type_, dtype, prefix)
+  module Bottle::Util::Indexing
+    include Bottle::Util::Exceptions
+    extend self
+
+    def get_matrix_row_at_index(matrix : Pointer({{ type_ }}), index : Int32)
+      vv = LibGsl.{{ prefix }}_row(matrix, index)
+      return Vector.new vv.vector, vv.vector.data
+    end
+
+    def get_matrix_col_at_index(matrix : Pointer({{ type_ }}), column : Int32)
+      vv = LibGsl.{{ prefix }}_column(matrix, column)
+      return Vector.new vv.vector, vv.vector.data
+    end
+  end
+end
+
+matrix_indexing_abstract LibGsl::GslMatrix, Float64, gsl_matrix
+matrix_indexing_abstract LibGsl::GslMatrixInt, Int32, gsl_matrix_int
diff --git a/src/core/vector.cr b/src/core/vector.cr
deleted file mode 100644
index 380a2c81..00000000
--- a/src/core/vector.cr
+++ /dev/null
@@ -1,336 +0,0 @@
-require "../llib/gsl"
-require "../util/indexing"
-require "../util/arithmetic"
-require "../util/statistics"
-require "../util/generic"
-require "../util/blas_arithmetic"
-require "../util/ma.cr"
-# require "../ma/comparison"
-
-class Vector(T, D)
-  @ptr : Pointer(T)
-  @obj : T
-  @owner : Int32
-  @size : UInt64
-  @stride : UInt64
-  @data : Pointer(D)
-  @dtype : D.class
-
-  getter ptr
-  getter owner
-  getter size
-  getter stride
-  getter data
-  getter dtype
-
-  def self.new(data : Array(A)) forall A
-    new(*fetch_struct(data))
-  end
-
-  private def self.fetch_struct(items : Array(Int32))
-    vector = LibGsl.gsl_vector_int_alloc(items.size)
-    items.each_with_index { |e, i| LibGsl.gsl_vector_int_set(vector, i, e) }
-    return vector, vector.value.data
-  end
-
-  private def self.fetch_struct(items : Array(Float64))
-    vector = LibGsl.gsl_vector_alloc(items.size)
-    items.each_with_index { |e, i| LibGsl.gsl_vector_set(vector, i, e) }
-    return vector, vector.value.data
-  end
-
-  private def self.fetch_struct(items : Array(Float64 | Int32))
-    vector = LibGsl.gsl_vector_alloc(items.size)
-    items.each_with_index { |e, i| LibGsl.gsl_vector_set(vector, i, e) }
-    return vector, vector.value.data
-  end
-
-  def initialize(@ptr : Pointer(T), @data : Pointer(D))
-    @obj = @ptr.value
-    @owner = @obj.owner
-    @size = @obj.size
-    @stride = @obj.stride
-    @dtype = D
-  end
-
-  def initialize(@obj : T, @data : Pointer(D))
-    @ptr = pointerof(@obj)
-    @owner = @obj.owner
-    @size = @obj.size
-    @stride = @obj.stride
-    @dtype = D
-  end
-
-  def self.zeros(n : Int32, dtype : Float64.class)
-    vector = LibGsl.gsl_vector_calloc(n)
-    return Vector.new vector, vector.value.data
-  end
-
-  def self.zeros(n : Int32, dtype : Int32.class)
-    vector = LibGsl.gsl_vector_int_calloc(n)
-    return Vector.new vector, vector.value.data
-  end
-
-  def copy
-    Bottle::Util::Indexing.copy_vector(self)
-  end
-
-  # Gets a single element from a vector at a given index, the core
-  # indexing operation of a vector
-  #
-  # ```
-  # vec = Vector.new [1, 2, 3, 4, 5]
-  # vec[0] # => 1
-  # ```
-  def [](index : Int32)
-    Bottle::Util::Indexing.get_vector_element_at_index(@ptr, index)
-  end
-
-  # Gets multiple elements from a vector at given indexes.  This returns
-  # a `copy` since there is no way to create a contiguous slice of memory
-  #
-  # ```
-  # vec = Vector.new [1, 2, 3]
-  # vec[[1, 2]] # => [2, 3]
-  # ```
-  def [](indexes : Array(Int32))
-    Bottle::Util::Indexing.get_vector_elements_at_indexes(@ptr, indexes)
-  end
-
-  # Returns a view of a vector defined by a given range.  Currently only
-  # supports single strided ranges due to limitations of Crystal
-  #
-  # ```
-  # vec = Vector.new [1, 2, 3, 4, 5]
-  # vec[2...4] # => [3, 4]
-  # ```
-  def [](range : Range(Int32 | Nil, Int32 | Nil))
-    Bottle::Util::Indexing.get_vector_elements_at_range(@ptr, range, @size)
-  end
-
-  # Sets a single element from a vector at a given index
-  #
-  # ```
-  # vec = Vector.new [1, 2, 3]
-  # vec[0] = 10
-  # vec # => [10, 2, 3]
-  # ```
-  def []=(index : Int32, value : Number)
-    Bottle::Util::Indexing.set_vector_element_at_index(@ptr, index, value)
-  end
-
-  # Sets multiple elements of a vector by the given indexes.
-  #
-  # ```
-  # vec = Vector.new [1, 2, 3]
-  # vec[[0, 1]] = [10, 9]
-  # vec # => [10, 9, 3]
-  # ```
-  def []=(indexes : Array(Int32), values : Array(Number))
-    Bottle::Util::Indexing.set_vector_elements_at_indexes(@ptr, indexes, values)
-  end
-
-  # Sets elements of a vector to given values based on the given range
-  #
-  # ```
-  # vec = Vector.new [1, 2, 3, 4, 5]
-  # vec[1...] = [10, 9, 8, 7]
-  # vec # => [1, 10, 9, 8, 7]
-  # ```
-  def []=(range : Range(Int32 | Nil, Int32 | Nil), values : Array(Number))
-    Bottle::Util::Indexing.set_vector_elements_at_range(@ptr, range, @size, values)
-  end
-
-  # Elementwise addition of a vector to another equally sized vector
-  #
-  # ```
-  # v1 = Vector.new [1.0, 2.0, 3.0]
-  # v2 = Vector.new [2.0, 4.0, 6.0]
-  # v1 + v2 # => [3.0, 6.0, 9.0]
-  # ```
-  def +(other : Vector(T, D))
-    Bottle::Util::VectorMath.add(self.copy, other)
-  end
-
-  # Elementwise addition of a vector to a scalar
-  #
-  # ```
-  # v1 = Vector.new [1.0, 2.0, 3.0]
-  # v2 = 2
-  # v1 + v2 # => [3.0, 4.0, 5.0]
-  # ```
-  def +(other : Number)
-    Bottle::Util::VectorMath.add_constant(self.copy, other)
-  end
-
-  # Elementwise subtraction of a vector to another equally sized vector
-  #
-  # ```
-  # v1 = Vector.new [1.0, 2.0, 3.0]
-  # v2 = Vector.new [2.0, 4.0, 6.0]
-  # v1 - v2 # => [-1.0, -2.0, -3.0]
-  # ```
-  def -(other : Vector(T, D))
-    Bottle::Util::VectorMath.sub(self.copy, other)
-  end
-
-  # Elementwise subtraction of a vector with a scalar
-  #
-  # ```
-  # v1 = Vector.new [1.0, 2.0, 3.0]
-  # v2 = 2
-  # v1 - v2 # => [-1.0, 0.0, 1.0]
-  # ```
-  def -(other : Number)
-    Bottle::Util::VectorMath.sub_constant(self.copy, other)
-  end
-
-  # Elementwise multiplication of a vector to another equally sized vector
-  #
-  # ```
-  # v1 = Vector.new [1.0, 2.0, 3.0]
-  # v2 = Vector.new [2.0, 4.0, 6.0]
-  # v1 * v2 # => [3.0, 8.0, 18.0]
-  # ```
-  def *(other : Vector(T, D))
-    Bottle::Util::VectorMath.mul(self.copy, other)
-  end
-
-  # Elementwise multiplication of a vector to a scalar
-  #
-  # ```
-  # v1 = Vector.new [1.0, 2.0, 3.0]
-  # v2 = 2
-  # v1 + v2 # => [2.0, 4.0, 6.0]
-  # ```
-  def *(other : Number)
-    Bottle::Util::VectorMath.mul_constant(self.copy, other)
-  end
-
-  # Elementwise division of a vector to another equally sized vector
-  #
-  # ```
-  # v1 = Vector.new [1.0, 2.0, 3.0]
-  # v2 = Vector.new [2.0, 4.0, 6.0]
-  # v1 / v2 # => [0.5, 0.5, 0.5]
-  # ```
-  def /(other : Vector(T, D))
-    Bottle::Util::VectorMath.div(self.copy, other)
-  end
-
-  # Elementwise division of a vector to a scalar
-  #
-  # ```
-  # v1 = Vector.new [1.0, 2.0, 3.0]
-  # v2 = 2
-  # v1 / v2 # => [0.5, 1, 1.5]
-  # ```
-  def /(other : Number)
-    Bottle::Util::VectorMath.div_constant(self.copy, other)
-  end
-
-  def ==(other : Vector)
-    Bottle::Ma::Comparison.vector_equal(self, other)
-  end
-
-  # Computes the dot product of two vectors
-  #
-  # ```
-  # v1 = Vector.new [1.0, 2.0, 3.0]
-  # v2 = Vector.new [4.0, 5.0, 6.0]
-  # v1.dot(v2) # => 32
-  # ```
-  def dot(other : Vector)
-    Bottle::Util::VectorMath.vector_dot(self, other)
-  end
-
-  # Computes the euclidean norm of the vector
-  #
-  # ```
-  # vec = Vector.new [1.0, 2.0, 3.0]
-  # vec.norm # => 3.741657386773941
-  # ```
-  def norm
-    Bottle::Util::VectorMath.vector_norm(self)
-  end
-
-  # Computes the maximum value of a vector
-  #
-  # ```
-  # v = Vector.new [1, 2, 3, 4]
-  # v.max # => 4
-  # ```
-  def max
-    Bottle::Util::VectorStats.vector_max(@ptr)
-  end
-
-  # Computes the minimum value of a vector
-  #
-  # ```
-  # v = Vector.new [1, 2, 3, 4]
-  # v.min # => 1
-  # ```
-  def min
-    Bottle::Util::VectorStats.vector_min(@ptr)
-  end
-
-  # Computes the min and max values of a vector
-  #
-  # ```
-  # v = Vector.new [1, 2, 3, 4]
-  # v.ptpv # => {1, 4}
-  # ```
-  def ptpv
-    Bottle::Util::VectorStats.vector_ptpv(@ptr)
-  end
-
-  # Computes the "peak to peak" of a vector (max - min)
-  #
-  # ```
-  # v = Vector.new [1, 2, 3, 4]
-  # v.ptp # => 3
-  # ```
-  def ptp
-    Bottle::Util::VectorStats.vector_ptp(@ptr)
-  end
-
-  # Computes the index of the maximum value of a vector
-  #
-  # ```
-  # v = Vector.new [1, 2, 3, 4]
-  # v.idxmax # => 3
-  # ```
-  def idxmax
-    Bottle::Util::VectorStats.vector_idxmax(@ptr)
-  end
-
-  # Computes the index of the minimum value of a vector
-  #
-  # ```
-  # v = Vector.new [1, 2, 3, 4]
-  # v.idxmin # => 0
-  # ```
-  def idxmin
-    Bottle::Util::VectorStats.vector_idxmin(@ptr)
-  end
-
-  # Computes the indexes of the minimum and maximum values of a vector
-  #
-  # ```
-  # v = Vector.new [1, 2, 3, 4]
-  # v.ptpidx # => {0, 3}
-  # ```
-  def ptpidx
-    Bottle::Util::VectorStats.vector_ptpidx(@ptr)
-  end
-
-  def to_s(io)
-    vals = (0...@size).map { |i| Bottle::Util::Indexing.get_vector_element_at_index(@ptr, i) }
-    io << "[" << vals.join(", ") << "]"
-  end
-end
-
-v1 = Vector.new [1.0, 2.0, 3.0]
-v2 = Vector.new [1.0, 3.0, 3.0]
-
-puts v2 * v2
diff --git a/src/core/vector/blas.cr b/src/core/vector/blas.cr
new file mode 100644
index 00000000..4f654c48
--- /dev/null
+++ b/src/core/vector/blas.cr
@@ -0,0 +1,15 @@
+# require "../core/vector"
+# require "../llib/gsl"
+# require "../llib/cblas"
+#
+# module Bottle::Core::VectorBlas
+#   extend self
+#
+#   def vector_dot(vector : Vector(LibGsl::GslVector), other : Vector(LibGsl::GslVector))
+#     LibCblas.ddot(vector.size, vector.@obj.data, vector.stride, other.@obj.data, other.stride)
+#   end
+#
+#   def vector_norm(vector : Vector(LibGsl::GslVector))
+#     LibCblas.snrm2(vector.size, vector.@obj.data, vector.stride)
+#   end
+# end
diff --git a/src/util/indexing.cr b/src/core/vector/index.cr
similarity index 76%
rename from src/util/indexing.cr
rename to src/core/vector/index.cr
index 6fbcc08a..b7b21404 100644
--- a/src/util/indexing.cr
+++ b/src/core/vector/index.cr
@@ -1,15 +1,13 @@
-require "../llib/gsl"
-require "../core/vector"
-require "../core/matrix"
-require "./exceptions"
-
-# Creates indexing methods for each support data type of gsl.
-# Since all the calls are virtually identical, and since
-# only the pointers are passed around, this is an easy way
-# to avoid duplicate code for no reason.
-macro vector_indexing_abstract(type_, prefix)
-  module Bottle::Util::Indexing
-    include Bottle::Util::Exceptions
+require "../../libs/gsl"
+require "../../vector/*"
+
+# Creates indexing methods for each support data type of gsl.
+# Since all the calls are virtually identical, and since
+# only the pointers are passed around, this is an easy way
+# to avoid duplicate code for no reason.
+macro vector_indexing_abstract(type_, dtype, prefix)
+  module Bottle::Core::VectorIndex
+    include Bottle::Core::Exceptions
     extend self
 
     # Gets a single element from a vector at a given index, the core
@@ -107,34 +105,14 @@ macro vector_indexing_abstract(type_, prefix)
       end
     end
 
-    def copy_vector(vector : Vector({{ type_ }}))
+    def copy_vector(vector : Vector({{ type_ }}, {{ dtype }}))
       ptv = LibGsl.{{ prefix }}_alloc(vector.size)
       LibGsl.{{ prefix }}_memcpy(ptv, vector.ptr)
       return Vector.new ptv, ptv.value.data
     end
 
   end
-end
-
-vector_indexing_abstract LibGsl::GslVector, gsl_vector
-vector_indexing_abstract LibGsl::GslVectorInt, gsl_vector_int
-
-macro matrix_indexing_abstract(type_, prefix)
-  module Bottle::Util::Indexing
-    include Bottle::Util::Exceptions
-    extend self
-
-    def get_matrix_row_at_index(matrix : Pointer({{ type_ }}), index : Int32)
-      vv = LibGsl.{{ prefix }}_row(matrix, index)
-      return Vector.new vv.vector, vv.vector.data
-    end
+end
 
-    def get_matrix_col_at_index(matrix : Pointer({{ type_ }}), column : Int32)
-      vv = LibGsl.{{ prefix }}_column(matrix, column)
-      return Vector.new vv.vector, vv.vector.data
-    end
-  end
-end
-
-matrix_indexing_abstract LibGsl::GslMatrix, gsl_matrix
-matrix_indexing_abstract LibGsl::GslMatrixInt, gsl_matrix_int
+vector_indexing_abstract LibGsl::GslVector, Float64, gsl_vector
+vector_indexing_abstract LibGsl::GslVectorInt, Int32, gsl_vector_int
diff --git a/src/util/arithmetic.cr b/src/core/vector/math.cr
similarity index 68%
rename from src/util/arithmetic.cr
rename to src/core/vector/math.cr
index 35852ace..ebf3603f 100644
--- a/src/util/arithmetic.cr
+++ b/src/core/vector/math.cr
@@ -1,10 +1,10 @@
-require "../llib/gsl"
-require "../core/vector"
-require "./exceptions"
-
-macro vector_arithmetic_abstract(type_, prefix)
-  module Bottle::Util::VectorMath
-    include Bottle::Util::Exceptions
+require "../../libs/gsl"
+require "../../vector/*"
+require "../exceptions"
+
+macro vector_arithmetic_abstract(type_, dtype, prefix)
+  module Bottle::Core::VectorMath
+    include Bottle::Core::Exceptions
     extend self
 
     # Elementwise addition of a vector to another equally sized vector
@@ -14,7 +14,7 @@ macro vector_arithmetic_abstract(type_, prefix)
     # v2 = Vector.new [2.0, 4.0, 6.0]
     # v1 + v2 # => [3.0, 6.0, 9.0]
     # ```
-    def add(a : Vector({{ type_ }}), b : Vector)
+    def add(a : Vector({{ type_ }}, {{ dtype }}), b : Vector)
       LibGsl.{{ prefix }}_add(a.ptr, b.ptr)
       return a
     end
@@ -26,7 +26,7 @@ macro vector_arithmetic_abstract(type_, prefix)
     # v2 = 2
     # v1 + v2 # => [3.0, 4.0, 5.0]
     # ```
-    def add_constant(a : Vector({{ type_ }}), x : Number)
+    def add_constant(a : Vector({{ type_ }}, {{ dtype }}), x : Number)
       LibGsl.{{ prefix }}_add_constant(a.ptr, x)
       return a
     end
@@ -38,7 +38,7 @@ macro vector_arithmetic_abstract(type_, prefix)
     # v2 = Vector.new [2.0, 4.0, 6.0]
     # v1 - v2 # => [-1.0, -2.0, -3.0]
     # ```
-    def sub(a : Vector({{ type_ }}), b : Vector)
+    def sub(a : Vector({{ type_ }}, {{ dtype }}), b : Vector)
       LibGsl.{{ prefix }}_sub(a.ptr, b.ptr)
       return a
     end
@@ -50,7 +50,7 @@ macro vector_arithmetic_abstract(type_, prefix)
     # v2 = 2
     # v1 - v2 # => [-1.0, 0.0, 1.0]
     # ```
-    def sub_constant(a : Vector({{ type_ }}), x : Number)
+    def sub_constant(a : Vector({{ type_ }}, {{ dtype }}), x : Number)
       LibGsl.{{ prefix }}_add_constant(a.ptr, -x)
       return a
     end
@@ -62,7 +62,7 @@ macro vector_arithmetic_abstract(type_, prefix)
     # v2 = Vector.new [2.0, 4.0, 6.0]
     # v1 * v2 # => [3.0, 8.0, 18.0]
     # ```
-    def mul(a : Vector({{ type_ }}), b : Vector)
+    def mul(a : Vector({{ type_ }}, {{ dtype }}), b : Vector)
       LibGsl.{{ prefix }}_mul(a.ptr, b.ptr)
       return a
     end
@@ -74,7 +74,7 @@ macro vector_arithmetic_abstract(type_, prefix)
     # v2 = 2
     # v1 + v2 # => [2.0, 4.0, 6.0]
     # ```
-    def mul_constant(a : Vector({{ type_ }}), x : Number)
+    def mul_constant(a : Vector({{ type_ }}, {{ dtype }}), x : Number)
       LibGsl.{{ prefix }}_scale(a.ptr, x)
       return a
     end
@@ -86,7 +86,7 @@ macro vector_arithmetic_abstract(type_, prefix)
     # v2 = Vector.new [2.0, 4.0, 6.0]
     # v1 / v2 # => [0.5, 0.5, 0.5]
     # ```
-    def div(a : Vector({{ type_ }}), b : Vector)
+    def div(a : Vector({{ type_ }}, {{ dtype }}), b : Vector)
       LibGsl.{{ prefix }}_div(a.ptr, b.ptr)
       return a
     end
@@ -99,13 +99,13 @@ macro vector_arithmetic_abstract(type_, prefix)
     # v2 = 2
     # v1 / v2 # => [0.5, 1, 1.5]
     # ```
-    def div_constant(a : Vector({{ type_ }}), x : Number)
+    def div_constant(a : Vector({{ type_ }}, {{ dtype }}), x : Number)
       LibGsl.{{ prefix }}_scale(a.ptr, 1 / x)
       return a
     end
 
   end
-end
-
-vector_arithmetic_abstract LibGsl::GslVector, gsl_vector
-vector_arithmetic_abstract LibGsl::GslVectorInt, gsl_vector_int
+end
+
+vector_arithmetic_abstract LibGsl::GslVector, Float64, gsl_vector
+vector_arithmetic_abstract LibGsl::GslVectorInt, Int32, gsl_vector_int
diff --git a/src/util/statistics.cr b/src/core/vector/stats.cr
similarity index 87%
rename from src/util/statistics.cr
rename to src/core/vector/stats.cr
index 7349c6a4..6c3dfc1e 100644
--- a/src/util/statistics.cr
+++ b/src/core/vector/stats.cr
@@ -1,9 +1,8 @@
-require "../llib/gsl"
-require "../core/vector"
-
-macro statistics_abstract(type_, prefix)
-  module Bottle::Util::VectorStats
-    include Bottle::Util::Exceptions
+require "../../libs/gsl"
+require "../../vector/*"
+
+macro statistics_abstract(type_, prefix)
+  module Bottle::Core::VectorStats
     extend self
 
     def vector_max(a : Pointer({{ type_ }}))
@@ -41,7 +40,7 @@ macro statistics_abstract(type_, prefix)
       return imin, imax
     end
   end
-end
-
-statistics_abstract LibGsl::GslVector, gsl_vector
-statistics_abstract LibGsl::GslVectorInt, gsl_vector_int
+end
+
+statistics_abstract LibGsl::GslVector, gsl_vector
+statistics_abstract LibGsl::GslVectorInt, gsl_vector_int
diff --git a/src/core/vector/util.cr b/src/core/vector/util.cr
new file mode 100644
index 00000000..e69de29b
diff --git a/src/llib/cblas.cr b/src/libs/cblas.cr
similarity index 100%
rename from src/llib/cblas.cr
rename to src/libs/cblas.cr
diff --git a/src/llib/gsl.cr b/src/libs/gsl.cr
similarity index 98%
rename from src/llib/gsl.cr
rename to src/libs/gsl.cr
index 63a8484c..0ff543f1 100644
--- a/src/llib/gsl.cr
+++ b/src/libs/gsl.cr
@@ -230,8 +230,8 @@ lib LibGsl
   fun gsl_matrix_int_sub(a : GslMatrixInt*, b : GslMatrixInt*)
   fun gsl_matrix_int_mul_elements(a : GslMatrixInt*, b : GslMatrixInt*) : Integer
   fun gsl_matrix_int_div_elements(a : GslMatrixInt*, b : GslMatrixInt*) : Integer
-  fun gsl_matrix_int_scale(a : GslMatrixInt*, x : Integer) : Integer
-  fun gsl_matrix_int_add_constant(a : GslMatrixInt*, x : Integer)
+  fun gsl_matrix_int_scale(a : GslMatrixInt*, x : Double) : Integer
+  fun gsl_matrix_int_add_constant(a : GslMatrixInt*, x : Double)
   fun gsl_matrix_int_max(m : GslMatrixInt*) : Integer
   fun gsl_matrix_int_min(m : GslMatrixInt*) : Integer
   fun gsl_matrix_int_minmax(m : GslMatrixInt*, min_out : Integer*, max_out : Integer*)
diff --git a/src/llib/bottle.cr b/src/llib/bottle.cr
deleted file mode 100644
index 3a11118d..00000000
--- a/src/llib/bottle.cr
+++ /dev/null
@@ -1,6 +0,0 @@
-require "./gsl"
-
-@[Link(ldflags: "#{__DIR__}/../../ext/bottle.so")]
-lib LibBottle
-  fun gsl_vector_ma_equal(u : LibGsl::GslVector*, v : LibGsl::GslVector*) : LibGsl::GslVectorInt
-end
diff --git a/src/ma/comparison.cr b/src/ma/comparison.cr
deleted file mode 100644
index e0153e4c..00000000
--- a/src/ma/comparison.cr
+++ /dev/null
@@ -1,18 +0,0 @@
-require "../llib/gsl"
-require "../util/exceptions"
-require "../core/vector"
-require "../llib/bottle"
-
-macro mask_comparison_abstract(type_, dtype, prefix)
-  module Bottle::Ma::Comparison
-    include Bottle::Util::Exceptions
-    extend self
-
-    def vector_equal(a : Vector({{ type_ }}), b : Vector)
-      ptv = LibBottle.{{ type_ }}_ma_equal(a.ptr, b.ptr)
-      return Vector.new ptv, ptv.value.data
-    end
-  end
-end
-
-mask_comparison_abstract LibGsl::GslVector, Float64, gsl_vector
diff --git a/src/core/matrix.cr b/src/matrix/init.cr
similarity index 95%
rename from src/core/matrix.cr
rename to src/matrix/init.cr
index 442c2545..6dfb6e20 100644
--- a/src/core/matrix.cr
+++ b/src/matrix/init.cr
@@ -1,65 +1,65 @@
-require "../llib/gsl"
-require "../util/indexing"
-require "../util/arithmetic"
-require "../util/statistics"
-require "../util/generic"
-require "../util/blas_arithmetic"
-
-class Matrix(T)
-  @ptr : Pointer(T)
-  @obj : T
-  @owner : Int32
-  @tda : UInt64
-  @nrows : UInt64
-  @ncols : UInt64
-
-  getter ptr
-  getter owner
-  getter tda
-  getter nrows
-  getter ncols
-
-  def self.new(data : Array(Array(A))) forall A
-    new(fetch_struct(data))
-  end
-
-  def self.fetch_struct(data : Array(Array(Float64)))
-    nrows = data.size
-    ncols = data[0].size
-    ptm = LibGsl.gsl_matrix_alloc(nrows, ncols)
-    data.each_with_index do |row, i|
-      row.each_with_index do |col, j|
-        LibGsl.gsl_matrix_set(ptm, i, j, col)
-      end
-    end
-    return ptm
-  end
-
-  def self.fetch_struct(data : Array(Array(Int32)))
-    nrows = data.size
-    ncols = data[0].size
-    ptm = LibGsl.gsl_matrix_int_alloc(nrows, ncols)
-    data.each_with_index do |row, i|
-      row.each_with_index do |col, j|
-        LibGsl.gsl_matrix_int_set(ptm, i, j, col)
-      end
-    end
-    return ptm
-  end
-
-  private def initialize(@ptr : Pointer(T))
-    @obj = @ptr.value
-    @owner = @obj.owner
-    @tda = @obj.tda
-    @nrows = @obj.size1
-    @ncols = @obj.size2
-  end
-
-  def [](row : Int32)
-    Bottle::Util::Indexing.get_matrix_row_at_index(@ptr, row)
-  end
-
-  def [](range : Range(Nil, Nil), column : Int32)
-    Bottle::Util::Indexing.get_matrix_col_at_index(@ptr, column)
-  end
-end
+require "../llib/gsl"
+require "../util/indexing"
+require "../util/arithmetic"
+require "../util/statistics"
+require "../util/generic"
+require "../util/blas_arithmetic"
+
+class Matrix(T)
+  @ptr : Pointer(T)
+  @obj : T
+  @owner : Int32
+  @tda : UInt64
+  @nrows : UInt64
+  @ncols : UInt64
+
+  getter ptr
+  getter owner
+  getter tda
+  getter nrows
+  getter ncols
+
+  def self.new(data : Array(Array(A))) forall A
+    new(fetch_struct(data))
+  end
+
+  def self.fetch_struct(data : Array(Array(Float64)))
+    nrows = data.size
+    ncols = data[0].size
+    ptm = LibGsl.gsl_matrix_alloc(nrows, ncols)
+    data.each_with_index do |row, i|
+      row.each_with_index do |col, j|
+        LibGsl.gsl_matrix_set(ptm, i, j, col)
+      end
+    end
+    return ptm
+  end
+
+  def self.fetch_struct(data : Array(Array(Int32)))
+    nrows = data.size
+    ncols = data[0].size
+    ptm = LibGsl.gsl_matrix_int_alloc(nrows, ncols)
+    data.each_with_index do |row, i|
+      row.each_with_index do |col, j|
+        LibGsl.gsl_matrix_int_set(ptm, i, j, col)
+      end
+    end
+    return ptm
+  end
+
+  private def initialize(@ptr : Pointer(T))
+    @obj = @ptr.value
+    @owner = @obj.owner
+    @tda = @obj.tda
+    @nrows = @obj.size1
+    @ncols = @obj.size2
+  end
+
+  def [](row : Int32)
+    Bottle::Util::Indexing.get_matrix_row_at_index(@ptr, row)
+  end
+
+  def [](range : Range(Nil, Nil), column : Int32)
+    Bottle::Util::Indexing.get_matrix_col_at_index(@ptr, column)
+  end
+end
diff --git a/src/util/blas_arithmetic.cr b/src/util/blas_arithmetic.cr
deleted file mode 100644
index 491688e9..00000000
--- a/src/util/blas_arithmetic.cr
+++ /dev/null
@@ -1,15 +0,0 @@
-require "../core/vector"
-require "../llib/gsl"
-require "../llib/cblas"
-
-module Bottle::Util::VectorMath
-  extend self
-
-  def vector_dot(vector : Vector(LibGsl::GslVector), other : Vector(LibGsl::GslVector))
-    LibCblas.ddot(vector.size, vector.@obj.data, vector.stride, other.@obj.data, other.stride)
-  end
-
-  def vector_norm(vector : Vector(LibGsl::GslVector))
-    LibCblas.snrm2(vector.size, vector.@obj.data, vector.stride)
-  end
-end
diff --git a/src/util/generic.cr b/src/util/generic.cr
deleted file mode 100644
index c5a9a8df..00000000
--- a/src/util/generic.cr
+++ /dev/null
@@ -1,37 +0,0 @@
-# require "../core/vector"
-# require "../llib/gsl"
-#
-# DOUBLE_TYPE = LibGsl::GslVector
-# INTEGER_TYPE = LibGsl::GslVectorInt
-#
-# macro abstract_cast(type_, prefix)
-#   module Bottle::Util::Cast
-#     extend self
-#
-#     def vector_to_double(vector : Vector({{ type_ }}))
-#       if vector.is_a?(DOUBLE_TYPE)
-#         return vector.as(Vector(LibGsl::GslVector))
-#       end
-#       ptv = LibGsl.gsl_vector_alloc(vector.size)
-#       (0...vector.size).each do |i|
-#         LibGsl.gsl_vector_set(ptv, i, LibGsl.{{ prefix }}_get(vector.ptr, i))
-#       end
-#       vector = Vector.new(ptv)
-#       return vector.as(Vector(LibGsl::GslVector))
-#     end
-#
-#     def vector_to_int(vector : Vector({{ type_ }}))
-#       if vector.is_a?(INTEGER_TYPE)
-#         return vector
-#       end
-#       ptv = LibGsl.gsl_vector_int_alloc(vector.size)
-#       (0...vector.size).each do |i|
-#         LibGsl.gsl_vector_int_set(ptv, i, LibGsl.{{ prefix }}_get(vector.ptr, i))
-#       end
-#       return Vector.new ptv
-#     end
-#   end
-# end
-#
-# abstract_cast LibGsl::GslVectorInt, gsl_vector_int
-# abstract_cast LibGsl::GslVector, gsl_vector
diff --git a/src/util/ma.cr b/src/util/ma.cr
deleted file mode 100644
index 25de03a4..00000000
--- a/src/util/ma.cr
+++ /dev/null
@@ -1,23 +0,0 @@
-# require "../llib/gsl"
-# require "../core/vector"
-#
-# macro mask_abstract(type_, prefix)
-#   module Bottle::Util::VectorMask
-#     include Bottle::Util::Exceptions
-#     extend self
-#
-#     def vector_nonzero(a : Pointer({{ type_ }}))
-#       ptv = LibMa.{{ prefix }}_ma_nonzero(a)
-#       return Vector.new ptv
-#     end
-#
-#     def vector_gt(a : Vector({{ type_ }}), b : Vector({{ type_ }}))
-#       ptv = LibMa.{{ prefix }}_ma_gt(a.ptr, b.ptr)
-#       return Vector.new ptv
-#     end
-#
-#   end
-# end
-#
-# mask_abstract LibGsl::GslVector, gsl_vector
-# mask_abstract LibGsl::GslVectorInt, gsl_vector_int
diff --git a/src/vector/indexing.cr b/src/vector/indexing.cr
new file mode 100644
index 00000000..46cd7cb8
--- /dev/null
+++ b/src/vector/indexing.cr
@@ -0,0 +1,75 @@
+require "../libs/gsl"
+require "../core/vector/*"
+require "./*"
+
+class Vector(T, D)
+  def copy
+    Bottle::Core::VectorIndex.copy_vector(self)
+  end
+
+  # Gets a single element from a vector at a given index, the core
+  # indexing operation of a vector
+  #
+  # ```
+  # vec = Vector.new [1, 2, 3, 4, 5]
+  # vec[0] # => 1
+  # ```
+  def [](index : Int32)
+    Bottle::Core::VectorIndex.get_vector_element_at_index(@ptr, index)
+  end
+
+  # Gets multiple elements from a vector at given indexes.  This returns
+  # a `copy` since there is no way to create a contiguous slice of memory
+  #
+  # ```
+  # vec = Vector.new [1, 2, 3]
+  # vec[[1, 2]] # => [2, 3]
+  # ```
+  def [](indexes : Array(Int32))
+    Bottle::Core::VectorIndex.get_vector_elements_at_indexes(@ptr, indexes)
+  end
+
+  # Returns a view of a vector defined by a given range.  Currently only
+  # supports single strided ranges due to limitations of Crystal
+  #
+  # ```
+  # vec = Vector.new [1, 2, 3, 4, 5]
+  # vec[2...4] # => [3, 4]
+  # ```
+  def [](range : Range(Int32 | Nil, Int32 | Nil))
+    Bottle::Core::VectorIndex.get_vector_elements_at_range(@ptr, range, @size)
+  end
+
+  # Sets a single element from a vector at a given index
+  #
+  # ```
+  # vec = Vector.new [1, 2, 3]
+  # vec[0] = 10
+  # vec # => [10, 2, 3]
+  # ```
+  def []=(index : Int32, value : Number)
+    Bottle::Core::VectorIndex.set_vector_element_at_index(@ptr, index, value)
+  end
+
+  # Sets multiple elements of a vector by the given indexes.
+  #
+  # ```
+  # vec = Vector.new [1, 2, 3]
+  # vec[[0, 1]] = [10, 9]
+  # vec # => [10, 9, 3]
+  # ```
+  def []=(indexes : Array(Int32), values : Array(Number))
+    Bottle::Core::VectorIndex.set_vector_elements_at_indexes(@ptr, indexes, values)
+  end
+
+  # Sets elements of a vector to given values based on the given range
+  #
+  # ```
+  # vec = Vector.new [1, 2, 3, 4, 5]
+  # vec[1...] = [10, 9, 8, 7]
+  # vec # => [1, 10, 9, 8, 7]
+  # ```
+  def []=(range : Range(Int32 | Nil, Int32 | Nil), values : Array(Number))
+    Bottle::Core::VectorIndex.set_vector_elements_at_range(@ptr, range, @size, values)
+  end
+end
diff --git a/src/vector/init.cr b/src/vector/init.cr
new file mode 100644
index 00000000..8d463951
--- /dev/null
+++ b/src/vector/init.cr
@@ -0,0 +1,68 @@
+require "../libs/gsl"
+require "../core/vector/*"
+require "./*"
+
+class Vector(T, D)
+  @ptr : Pointer(T)
+  @obj : T
+  @owner : Int32
+  @size : UInt64
+  @stride : UInt64
+  @data : Pointer(D)
+  @dtype : D.class
+
+  getter ptr
+  getter owner
+  getter size
+  getter stride
+  getter data
+  getter dtype
+
+  def self.new(data : Array(A)) forall A
+    new(*fetch_struct(data))
+  end
+
+  private def self.fetch_struct(items : Array(Int32))
+    vector = LibGsl.gsl_vector_int_alloc(items.size)
+    items.each_with_index { |e, i| LibGsl.gsl_vector_int_set(vector, i, e) }
+    return vector, vector.value.data
+  end
+
+  private def self.fetch_struct(items : Array(Float64))
+    vector = LibGsl.gsl_vector_alloc(items.size)
+    items.each_with_index { |e, i| LibGsl.gsl_vector_set(vector, i, e) }
+    return vector, vector.value.data
+  end
+
+  private def self.fetch_struct(items : Array(Float64 | Int32))
+    vector = LibGsl.gsl_vector_alloc(items.size)
+    items.each_with_index { |e, i| LibGsl.gsl_vector_set(vector, i, e) }
+    return vector, vector.value.data
+  end
+
+  def initialize(@ptr : Pointer(T), @data : Pointer(D))
+    @obj = @ptr.value
+    @owner = @obj.owner
+    @size = @obj.size
+    @stride = @obj.stride
+    @dtype = D
+  end
+
+  def initialize(@obj : T, @data : Pointer(D))
+    @ptr = pointerof(@obj)
+    @owner = @obj.owner
+    @size = @obj.size
+    @stride = @obj.stride
+    @dtype = D
+  end
+
+  def self.zeros(n : Int32, dtype : Float64.class)
+    vector = LibGsl.gsl_vector_calloc(n)
+    return Vector.new vector, vector.value.data
+  end
+
+  def self.zeros(n : Int32, dtype : Int32.class)
+    vector = LibGsl.gsl_vector_int_calloc(n)
+    return Vector.new vector, vector.value.data
+  end
+end
diff --git a/src/vector/math.cr b/src/vector/math.cr
new file mode 100644
index 00000000..96394789
--- /dev/null
+++ b/src/vector/math.cr
@@ -0,0 +1,93 @@
+require "../libs/gsl"
+require "../core/vector/*"
+require "./*"
+
+class Vector(T, D)
+  # Elementwise addition of a vector to another equally sized vector
+  #
+  # ```
+  # v1 = Vector.new [1.0, 2.0, 3.0]
+  # v2 = Vector.new [2.0, 4.0, 6.0]
+  # v1 + v2 # => [3.0, 6.0, 9.0]
+  # ```
+  def +(other : Vector(T, D))
+    Bottle::Core::VectorMath.add(self.copy, other)
+  end
+
+  # Elementwise addition of a vector to a scalar
+  #
+  # ```
+  # v1 = Vector.new [1.0, 2.0, 3.0]
+  # v2 = 2
+  # v1 + v2 # => [3.0, 4.0, 5.0]
+  # ```
+  def +(other : Number)
+    Bottle::Core::VectorMath.add_constant(self.copy, other.to_f)
+  end
+
+  # Elementwise subtraction of a vector to another equally sized vector
+  #
+  # ```
+  # v1 = Vector.new [1.0, 2.0, 3.0]
+  # v2 = Vector.new [2.0, 4.0, 6.0]
+  # v1 - v2 # => [-1.0, -2.0, -3.0]
+  # ```
+  def -(other : Vector(T, D))
+    Bottle::Core::VectorMath.sub(self.copy, other)
+  end
+
+  # Elementwise subtraction of a vector with a scalar
+  #
+  # ```
+  # v1 = Vector.new [1.0, 2.0, 3.0]
+  # v2 = 2
+  # v1 - v2 # => [-1.0, 0.0, 1.0]
+  # ```
+  def -(other : Number)
+    Bottle::Core::VectorMath.sub_constant(self.copy, other.to_f)
+  end
+
+  # Elementwise multiplication of a vector to another equally sized vector
+  #
+  # ```
+  # v1 = Vector.new [1.0, 2.0, 3.0]
+  # v2 = Vector.new [2.0, 4.0, 6.0]
+  # v1 * v2 # => [3.0, 8.0, 18.0]
+  # ```
+  def *(other : Vector(T, D))
+    Bottle::Core::VectorMath.mul(self.copy, other)
+  end
+
+  # Elementwise multiplication of a vector to a scalar
+  #
+  # ```
+  # v1 = Vector.new [1.0, 2.0, 3.0]
+  # v2 = 2
+  # v1 + v2 # => [2.0, 4.0, 6.0]
+  # ```
+  def *(other : Number)
+    Bottle::Core::VectorMath.mul_constant(self.copy, other.to_f)
+  end
+
+  # Elementwise division of a vector to another equally sized vector
+  #
+  # ```
+  # v1 = Vector.new [1.0, 2.0, 3.0]
+  # v2 = Vector.new [2.0, 4.0, 6.0]
+  # v1 / v2 # => [0.5, 0.5, 0.5]
+  # ```
+  def /(other : Vector(T, D))
+    Bottle::Core::VectorMath.div(self.copy, other)
+  end
+
+  # Elementwise division of a vector to a scalar
+  #
+  # ```
+  # v1 = Vector.new [1.0, 2.0, 3.0]
+  # v2 = 2
+  # v1 / v2 # => [0.5, 1, 1.5]
+  # ```
+  def /(other : Number)
+    Bottle::Core::VectorMath.div_constant(self.copy, other.to_f)
+  end
+end
diff --git a/src/vector/stats.cr b/src/vector/stats.cr
new file mode 100644
index 00000000..78a482bb
--- /dev/null
+++ b/src/vector/stats.cr
@@ -0,0 +1,100 @@
+require "../core/vector/*"
+require "./*"
+
+class Vector(T, D)
+  # Computes the dot product of two vectors
+  #
+  # ```
+  # v1 = Vector.new [1.0, 2.0, 3.0]
+  # v2 = Vector.new [4.0, 5.0, 6.0]
+  # v1.dot(v2) # => 32
+  # ```
+  def dot(other : Vector)
+    Bottle::Core::Vector::Stats.vector_dot(self, other)
+  end
+
+  # Computes the euclidean norm of the vector
+  #
+  # ```
+  # vec = Vector.new [1.0, 2.0, 3.0]
+  # vec.norm # => 3.741657386773941
+  # ```
+  def norm
+    Bottle::Core::Vector::Stats.vector_norm(self)
+  end
+
+  # Computes the maximum value of a vector
+  #
+  # ```
+  # v = Vector.new [1, 2, 3, 4]
+  # v.max # => 4
+  # ```
+  def max
+    Bottle::Core::Vector::Stats.vector_max(@ptr)
+  end
+
+  # Computes the minimum value of a vector
+  #
+  # ```
+  # v = Vector.new [1, 2, 3, 4]
+  # v.min # => 1
+  # ```
+  def min
+    Bottle::Core::Vector::Stats.vector_min(@ptr)
+  end
+
+  # Computes the min and max values of a vector
+  #
+  # ```
+  # v = Vector.new [1, 2, 3, 4]
+  # v.ptpv # => {1, 4}
+  # ```
+  def ptpv
+    Bottle::Core::Vector::Stats.vector_ptpv(@ptr)
+  end
+
+  # Computes the "peak to peak" of a vector (max - min)
+  #
+  # ```
+  # v = Vector.new [1, 2, 3, 4]
+  # v.ptp # => 3
+  # ```
+  def ptp
+    Bottle::Core::Vector::Stats.vector_ptp(@ptr)
+  end
+
+  # Computes the index of the maximum value of a vector
+  #
+  # ```
+  # v = Vector.new [1, 2, 3, 4]
+  # v.idxmax # => 3
+  # ```
+  def idxmax
+    Bottle::Core::Vector::Stats.vector_idxmax(@ptr)
+  end
+
+  # Computes the index of the minimum value of a vector
+  #
+  # ```
+  # v = Vector.new [1, 2, 3, 4]
+  # v.idxmin # => 0
+  # ```
+  def idxmin
+    Bottle::Core::Vector::Stats.vector_idxmin(@ptr)
+  end
+
+  # Computes the indexes of the minimum and maximum values of a vector
+  #
+  # ```
+  # v = Vector.new [1, 2, 3, 4]
+  # v.ptpidx # => {0, 3}
+  # ```
+  def ptpidx
+    Bottle::Core::Vector::Stats.vector_ptpidx(@ptr)
+  end
+
+  def to_s(io)
+    vals = (0...@size).map { |i| Bottle::Core::VectorIndex.get_vector_element_at_index(@ptr, i) }
+    io << "[" << vals.join(", ") << "]"
+  end
+end