First version of alignment functions (imported from suma* programs)

python/obitools3/obidms/_obidms.cfiles

@@ -10,6 +10,8 @@
 ../../../src/encode.c
 ../../../src/murmurhash2.h
 ../../../src/murmurhash2.c
+../../../src/obi_align.h
+../../../src/obi_align.c
 ../../../src/obiavl.h
 ../../../src/obiavl.c
 ../../../src/obiblob_indexer.h
@@ -37,6 +39,8 @@
 ../../../src/obiview.c
 ../../../src/utils.h
 ../../../src/utils.c
+../../../src/sse_banded_LCS_alignment.h
+../../../src/sse_banded_LCS_alignment.c
 ../../../src/obidmscolumn_bool.c
 ../../../src/obidmscolumn_bool.h
 ../../../src/obidmscolumn_char.c

python/obitools3/obidms/_obidms.pxd

@@ -6,11 +6,12 @@ from .capi.obiview      cimport Obiview_p
 from .capi.obitypes     cimport obiversion_t, OBIType_t, index_t
 from ._obitaxo          cimport OBI_Taxonomy
 
+
 cdef class OBIDMS_column:
 
     cdef  OBIDMS_column_p*   pointer
     cdef  OBIDMS             dms
-    cdef  Obiview_p          view
+    cdef  OBIView            view
     cdef  str                data_type
     cdef  str                dms_name
     cdef  str                column_name

python/obitools3/obidms/_obidms.pyx

@@ -63,6 +63,7 @@ from .capi.obiview cimport Obiview_p, \
                            obi_view_delete_column, \
                            obi_view_add_column, \
                            obi_view_get_column, \
+                           obi_view_get_column, \
                            obi_view_get_pointer_on_column_in_view, \
                            obi_select_line, \
                            obi_select_lines, \
@@ -90,7 +91,7 @@ cdef class OBIDMS_column :
         # Fill structure
         self.pointer = column_pp
         self.dms = view.dms
-        self.view = view.pointer
+        self.view = view
         self.data_type = bytes2str(name_data_type((column_p.header).returned_data_type))
         self.column_name = bytes2str((column_p.header).name)
         self.nb_elements_per_line = (column_p.header).nb_elements_per_line
@@ -120,7 +121,7 @@ cdef class OBIDMS_column :
             yield self.get_line(line_nb)
     
     cpdef update_pointer(self):
-        self.pointer = <OBIDMS_column_p*> obi_view_get_pointer_on_column_in_view(self.view, str2bytes(self.column_name))
+        self.pointer = <OBIDMS_column_p*> obi_view_get_pointer_on_column_in_view(self.view.pointer, str2bytes(self.column_name))
 
     cpdef list get_elements_names(self):
         return self.elements_names
@@ -336,7 +337,6 @@ cdef class OBIView :
         cdef bytes column_name_b
         cdef bytes elements_names_b
         cdef object subclass
-        cdef OBIDMS_column_p* column_pp
         cdef OBIDMS_column_p column_p
         
         column_name_b = str2bytes(column_name)
@@ -369,10 +369,9 @@ cdef class OBIView :
             raise Exception("Problem adding a column in a view")
         
         # Get the column pointer
-        column_pp = obi_view_get_pointer_on_column_in_view(self.pointer, column_name_b)
+        column_p = obi_view_get_column(self.pointer, column_name_b)
         
         # Open and store the subclass
-        column_p = column_pp[0]     # TODO ugly cython dereferencing
         subclass = OBIDMS_column.get_subclass_type(column_p)
         (self.columns)[column_name] = subclass(self, column_name)
                 
@@ -495,18 +494,15 @@ cdef class OBIView_NUC_SEQS(OBIView):
     cpdef delete_column(self, str column_name) :
         
         cdef int i
-        cdef Obiview_p view
+        cdef Obiview_p view_p
         cdef OBIDMS_column column
         cdef OBIDMS_column_p column_p
         cdef OBIDMS_column_header_p header
         cdef str column_n
  
-        if ((column_name == bytes2str(ID_COLUMN)) or (column_name == bytes2str(NUC_SEQUENCE_COLUMN)) or (column_name == bytes2str(DEFINITION_COLUMN))) :
+        view_p = self.pointer
-            raise Exception("Can't delete an obligatory column from a NUC_SEQS view")
-        
-        view = self.pointer
  
-        if obi_view_delete_column(view, str2bytes(column_name)) < 0 :
+        if obi_view_delete_column(view_p, str2bytes(column_name)) < 0 :
             raise Exception("Problem deleting a column from a view")
         
         # Remove instance from the dictionary

python/obitools3/obidms/_obidmscolumn_bool.cfiles

@@ -10,6 +10,8 @@
 ../../../src/encode.c
 ../../../src/murmurhash2.h
 ../../../src/murmurhash2.c
+../../../src/obi_align.h
+../../../src/obi_align.c
 ../../../src/obiavl.h
 ../../../src/obiavl.c
 ../../../src/obiblob_indexer.h
@@ -37,6 +39,8 @@
 ../../../src/obiview.c
 ../../../src/utils.h
 ../../../src/utils.c
+../../../src/sse_banded_LCS_alignment.h
+../../../src/sse_banded_LCS_alignment.c
 ../../../src/obidmscolumn_bool.c
 ../../../src/obidmscolumn_bool.h
 ../../../src/obidmscolumn_char.c

python/obitools3/obidms/_obidmscolumn_bool.pyx

@@ -17,7 +17,7 @@ cdef class OBIDMS_column_bool(OBIDMS_column):
     cpdef object get_line(self, index_t line_nb):
         cdef obibool_t value
         cdef object result
-        value = obi_column_get_obibool_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, 0)
+        value = obi_column_get_obibool_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, 0)
         if obi_errno > 0 :
             raise IndexError(line_nb)
         if value == OBIBool_NA :
@@ -29,7 +29,7 @@ cdef class OBIDMS_column_bool(OBIDMS_column):
     cpdef set_line(self, index_t line_nb, object value):
         if value is None :
             value = OBIBool_NA
-        if obi_column_set_obibool_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, 0, <obibool_t> value) < 0:
+        if obi_column_set_obibool_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, 0, <obibool_t> value) < 0:
             raise Exception("Problem setting a value in a column")
 
 
@@ -38,7 +38,7 @@ cdef class OBIDMS_column_multi_elts_bool(OBIDMS_column_multi_elts):
     cpdef object get_item(self, index_t line_nb, str element_name):
         cdef obibool_t value
         cdef object result
-        value = obi_column_get_obibool_with_elt_name_in_view(self.view, (self.pointer)[0], line_nb, str2bytes(element_name))
+        value = obi_column_get_obibool_with_elt_name_in_view(self.view.pointer, (self.pointer)[0], line_nb, str2bytes(element_name))
         if obi_errno > 0 :
             raise IndexError(line_nb, element_name)
         if value == OBIBool_NA :
@@ -56,7 +56,7 @@ cdef class OBIDMS_column_multi_elts_bool(OBIDMS_column_multi_elts):
         result = {}
         all_NA = True
         for i in range(self.nb_elements_per_line) :
-            value = obi_column_get_obibool_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, i)
+            value = obi_column_get_obibool_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, i)
             if obi_errno > 0 :
                 raise IndexError(line_nb)
             if value == OBIBool_NA :
@@ -73,5 +73,5 @@ cdef class OBIDMS_column_multi_elts_bool(OBIDMS_column_multi_elts):
     cpdef set_item(self, index_t line_nb, str element_name, object value):
         if value is None :
             value = OBIBool_NA
-        if obi_column_set_obibool_with_elt_name_in_view(self.view, (self.pointer)[0], line_nb, str2bytes(element_name), <obibool_t> value) < 0:
+        if obi_column_set_obibool_with_elt_name_in_view(self.view.pointer, (self.pointer)[0], line_nb, str2bytes(element_name), <obibool_t> value) < 0:
             raise Exception("Problem setting a value in a column")

python/obitools3/obidms/_obidmscolumn_char.cfiles

@@ -10,6 +10,8 @@
 ../../../src/encode.c
 ../../../src/murmurhash2.h
 ../../../src/murmurhash2.c
+../../../src/obi_align.h
+../../../src/obi_align.c
 ../../../src/obiavl.h
 ../../../src/obiavl.c
 ../../../src/obiblob_indexer.h
@@ -37,6 +39,8 @@
 ../../../src/obiview.c
 ../../../src/utils.h
 ../../../src/utils.c
+../../../src/sse_banded_LCS_alignment.h
+../../../src/sse_banded_LCS_alignment.c
 ../../../src/obidmscolumn_bool.c
 ../../../src/obidmscolumn_bool.h
 ../../../src/obidmscolumn_char.c

python/obitools3/obidms/_obidmscolumn_char.pyx

@@ -15,7 +15,7 @@ cdef class OBIDMS_column_char(OBIDMS_column):
     cpdef object get_line(self, index_t line_nb):
         cdef obichar_t value
         cdef object result
-        value = obi_column_get_obichar_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, 0)
+        value = obi_column_get_obichar_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, 0)
         if obi_errno > 0 :
             raise IndexError(line_nb)
         if value == OBIChar_NA :
@@ -27,7 +27,7 @@ cdef class OBIDMS_column_char(OBIDMS_column):
     cpdef set_line(self, index_t line_nb, object value):
         if value is None :
             value = OBIChar_NA
-        if obi_column_set_obichar_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, 0, str2bytes(value)[0]) < 0:
+        if obi_column_set_obichar_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, 0, str2bytes(value)[0]) < 0:
             raise Exception("Problem setting a value in a column")
 
 
@@ -36,7 +36,7 @@ cdef class OBIDMS_column_multi_elts_char(OBIDMS_column_multi_elts):
     cpdef object get_item(self, index_t line_nb, str element_name):
         cdef obichar_t value
         cdef object result
-        value = obi_column_get_obichar_with_elt_name_in_view(self.view, (self.pointer)[0], line_nb, str2bytes(element_name))
+        value = obi_column_get_obichar_with_elt_name_in_view(self.view.pointer, (self.pointer)[0], line_nb, str2bytes(element_name))
         if obi_errno > 0 :
             raise IndexError(line_nb, element_name)
         if value == OBIChar_NA :
@@ -54,7 +54,7 @@ cdef class OBIDMS_column_multi_elts_char(OBIDMS_column_multi_elts):
         result = {}
         all_NA = True
         for i in range(self.nb_elements_per_line) :
-            value = obi_column_get_obichar_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, i)
+            value = obi_column_get_obichar_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, i)
             if obi_errno > 0 :
                 raise IndexError(line_nb)
             if value == OBIChar_NA :
@@ -71,6 +71,6 @@ cdef class OBIDMS_column_multi_elts_char(OBIDMS_column_multi_elts):
     cpdef set_item(self, index_t line_nb, str element_name, object value):
         if value is None :
             value = OBIChar_NA
-        if obi_column_set_obichar_with_elt_name_in_view(self.view, (self.pointer)[0], line_nb, str2bytes(element_name), str2bytes(value)[0]) < 0:
+        if obi_column_set_obichar_with_elt_name_in_view(self.view.pointer, (self.pointer)[0], line_nb, str2bytes(element_name), str2bytes(value)[0]) < 0:
             raise Exception("Problem setting a value in a column")
 

python/obitools3/obidms/_obidmscolumn_float.cfiles

@@ -10,6 +10,8 @@
 ../../../src/encode.c
 ../../../src/murmurhash2.h
 ../../../src/murmurhash2.c
+../../../src/obi_align.h
+../../../src/obi_align.c
 ../../../src/obiavl.h
 ../../../src/obiavl.c
 ../../../src/obiblob_indexer.h
@@ -37,6 +39,8 @@
 ../../../src/obiview.c
 ../../../src/utils.h
 ../../../src/utils.c
+../../../src/sse_banded_LCS_alignment.h
+../../../src/sse_banded_LCS_alignment.c
 ../../../src/obidmscolumn_bool.c
 ../../../src/obidmscolumn_bool.h
 ../../../src/obidmscolumn_char.c

python/obitools3/obidms/_obidmscolumn_float.pyx

@@ -15,7 +15,7 @@ cdef class OBIDMS_column_float(OBIDMS_column):
     cpdef object get_line(self, index_t line_nb):
         cdef obifloat_t value
         cdef object result
-        value = obi_column_get_obifloat_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, 0)
+        value = obi_column_get_obifloat_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, 0)
         if obi_errno > 0 :
             raise IndexError(line_nb)
         if value == OBIFloat_NA :
@@ -27,7 +27,7 @@ cdef class OBIDMS_column_float(OBIDMS_column):
     cpdef set_line(self, index_t line_nb, object value):
         if value is None :
             value = OBIFloat_NA
-        if obi_column_set_obifloat_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, 0, <obifloat_t> value) < 0:
+        if obi_column_set_obifloat_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, 0, <obifloat_t> value) < 0:
             raise Exception("Problem setting a value in a column")
 
 
@@ -36,7 +36,7 @@ cdef class OBIDMS_column_multi_elts_float(OBIDMS_column_multi_elts):
     cpdef object get_item(self, index_t line_nb, str element_name):
         cdef obifloat_t value
         cdef object result
-        value = obi_column_get_obifloat_with_elt_name_in_view(self.view, (self.pointer)[0], line_nb, str2bytes(element_name))
+        value = obi_column_get_obifloat_with_elt_name_in_view(self.view.pointer, (self.pointer)[0], line_nb, str2bytes(element_name))
         if obi_errno > 0 :
             raise IndexError(line_nb, element_name)
         if value == OBIFloat_NA :
@@ -54,7 +54,7 @@ cdef class OBIDMS_column_multi_elts_float(OBIDMS_column_multi_elts):
         result = {}
         all_NA = True
         for i in range(self.nb_elements_per_line) :
-            value = obi_column_get_obifloat_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, i)
+            value = obi_column_get_obifloat_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, i)
             if obi_errno > 0 :
                 raise IndexError(line_nb)
             if value == OBIFloat_NA :
@@ -71,6 +71,6 @@ cdef class OBIDMS_column_multi_elts_float(OBIDMS_column_multi_elts):
     cpdef set_item(self, index_t line_nb, str element_name, object value):
         if value is None :
             value = OBIFloat_NA
-        if obi_column_set_obifloat_with_elt_name_in_view(self.view, (self.pointer)[0], line_nb, str2bytes(element_name), <obifloat_t> value) < 0:
+        if obi_column_set_obifloat_with_elt_name_in_view(self.view.pointer, (self.pointer)[0], line_nb, str2bytes(element_name), <obifloat_t> value) < 0:
             raise Exception("Problem setting a value in a column")
  

python/obitools3/obidms/_obidmscolumn_int.cfiles

@@ -10,6 +10,8 @@
 ../../../src/encode.c
 ../../../src/murmurhash2.h
 ../../../src/murmurhash2.c
+../../../src/obi_align.h
+../../../src/obi_align.c
 ../../../src/obiavl.h
 ../../../src/obiavl.c
 ../../../src/obiblob_indexer.h
@@ -37,6 +39,8 @@
 ../../../src/obiview.c
 ../../../src/utils.h
 ../../../src/utils.c
+../../../src/sse_banded_LCS_alignment.h
+../../../src/sse_banded_LCS_alignment.c
 ../../../src/obidmscolumn_bool.c
 ../../../src/obidmscolumn_bool.h
 ../../../src/obidmscolumn_char.c

python/obitools3/obidms/_obidmscolumn_int.pyx

@@ -17,7 +17,7 @@ cdef class OBIDMS_column_int(OBIDMS_column):
     cpdef object get_line(self, index_t line_nb):
         cdef obiint_t value
         cdef object result
-        value = obi_column_get_obiint_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, 0)
+        value = obi_column_get_obiint_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, 0)
         if obi_errno > 0 :
             raise IndexError(line_nb)
         if value == OBIInt_NA :
@@ -29,7 +29,7 @@ cdef class OBIDMS_column_int(OBIDMS_column):
     cpdef set_line(self, index_t line_nb, object value):
         if value is None :
             value = OBIInt_NA
-        if obi_column_set_obiint_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, 0, <obiint_t> value) < 0:
+        if obi_column_set_obiint_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, 0, <obiint_t> value) < 0:
             raise Exception("Problem setting a value in a column")
 
 
@@ -38,7 +38,7 @@ cdef class OBIDMS_column_multi_elts_int(OBIDMS_column_multi_elts):
     cpdef object get_item(self, index_t line_nb, str element_name):
         cdef obiint_t value
         cdef object result
-        value = obi_column_get_obiint_with_elt_name_in_view(self.view, (self.pointer)[0], line_nb, str2bytes(element_name))
+        value = obi_column_get_obiint_with_elt_name_in_view(self.view.pointer, (self.pointer)[0], line_nb, str2bytes(element_name))
         if obi_errno > 0 :
             raise IndexError(line_nb, element_name)
         if value == OBIInt_NA :
@@ -56,7 +56,7 @@ cdef class OBIDMS_column_multi_elts_int(OBIDMS_column_multi_elts):
         result = {}
         all_NA = True
         for i in range(self.nb_elements_per_line) :
-            value = obi_column_get_obiint_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, i)
+            value = obi_column_get_obiint_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, i)
             if obi_errno > 0 :
                 raise IndexError(line_nb)
             if value == OBIInt_NA :
@@ -73,6 +73,6 @@ cdef class OBIDMS_column_multi_elts_int(OBIDMS_column_multi_elts):
     cpdef set_item(self, index_t line_nb, str element_name, object value):
         if value is None :
             value = OBIInt_NA
-        if obi_column_set_obiint_with_elt_name_in_view(self.view, (self.pointer)[0], line_nb, str2bytes(element_name), <obiint_t> value) < 0:
+        if obi_column_set_obiint_with_elt_name_in_view(self.view.pointer, (self.pointer)[0], line_nb, str2bytes(element_name), <obiint_t> value) < 0:
             raise Exception("Problem setting a value in a column")
 

python/obitools3/obidms/_obidmscolumn_seq.cfiles

@@ -10,6 +10,8 @@
 ../../../src/encode.c
 ../../../src/murmurhash2.h
 ../../../src/murmurhash2.c
+../../../src/obi_align.h
+../../../src/obi_align.c
 ../../../src/obiavl.h
 ../../../src/obiavl.c
 ../../../src/obiblob_indexer.h
@@ -37,6 +39,8 @@
 ../../../src/obiview.c
 ../../../src/utils.h
 ../../../src/utils.c
+../../../src/sse_banded_LCS_alignment.h
+../../../src/sse_banded_LCS_alignment.c
 ../../../src/obidmscolumn_bool.c
 ../../../src/obidmscolumn_bool.h
 ../../../src/obidmscolumn_char.c

python/obitools3/obidms/_obidmscolumn_seq.pxd

@@ -1,12 +1,20 @@
 #cython: language_level=3
 
 from .capi.obitypes cimport index_t
-from ._obidms       cimport OBIDMS_column, OBIDMS_column_multi_elts
+from ._obidms       cimport OBIView, OBIDMS_column, OBIDMS_column_multi_elts
 
 
 cdef class OBIDMS_column_seq(OBIDMS_column):
     cpdef object get_line(self, index_t line_nb)
     cpdef set_line(self, index_t line_nb, object value)
+    cpdef align(self, 
+                OBIView score_view,
+                OBIDMS_column score_column, 
+                double threshold = *, 
+                bint normalize = *, 
+                int reference = *,
+                bint similarity_mode = *)
+
 
 cdef class OBIDMS_column_multi_elts_seq(OBIDMS_column_multi_elts):
     cpdef object get_item(self, index_t line_nb, str element_name)

python/obitools3/obidms/_obidmscolumn_seq.pyx

@@ -4,9 +4,12 @@ from .capi.obiview      cimport obi_column_get_obiseq_with_elt_name_in_view, \
                                 obi_column_get_obiseq_with_elt_idx_in_view, \
                                 obi_column_set_obiseq_with_elt_name_in_view, \
                                 obi_column_set_obiseq_with_elt_idx_in_view
+from .capi.obialign     cimport obi_align_one_column
 from .capi.obierrno     cimport obi_errno
 from .capi.obitypes     cimport OBISeq_NA, const_char_p
 
+from ._obidms           cimport OBIView
+
 from obitools3.utils cimport str2bytes, bytes2str
 
 from libc.stdlib cimport free
@@ -18,7 +21,7 @@ cdef class OBIDMS_column_seq(OBIDMS_column):
     cpdef object get_line(self, index_t line_nb):
         cdef char*  value
         cdef object result
-        value = obi_column_get_obiseq_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, 0)
+        value = obi_column_get_obiseq_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, 0)
         if obi_errno > 0 :
             raise IndexError(line_nb)
         if strcmp(value, OBISeq_NA) == 0 :
@@ -34,16 +37,28 @@ cdef class OBIDMS_column_seq(OBIDMS_column):
             value_b = OBISeq_NA
         else :
             value_b = str2bytes(value)
-        if obi_column_set_obiseq_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, 0, value_b) < 0:
+        if obi_column_set_obiseq_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, 0, value_b) < 0:
             raise Exception("Problem setting a value in a column")
+    
+    # TODO choose alignment type (lcs or other) with supplementary argument
+    cpdef align(self, 
+                OBIView score_view,
+                OBIDMS_column score_column, 
+                double threshold = 0.0, 
+                bint normalize = True, 
+                int reference = 0,      # TODO
+                bint similarity_mode = True):
+        if (obi_align_one_column(self.view.pointer, (self.pointer)[0], score_view.pointer, (score_column.pointer)[0], threshold, normalize, reference, similarity_mode) < 0) :
+            raise Exception("An error occurred while aligning sequences")
 
+        
 
 cdef class OBIDMS_column_multi_elts_seq(OBIDMS_column_multi_elts):
 
     cpdef object get_item(self, index_t line_nb, str element_name):
         cdef char*  value
         cdef object result
-        value = obi_column_get_obiseq_with_elt_name_in_view(self.view, (self.pointer)[0], line_nb, str2bytes(element_name))
+        value = obi_column_get_obiseq_with_elt_name_in_view(self.view.pointer, (self.pointer)[0], line_nb, str2bytes(element_name))
         if obi_errno > 0 :
             raise IndexError(line_nb, element_name)
         if strcmp(value, OBISeq_NA) == 0 :
@@ -62,7 +77,7 @@ cdef class OBIDMS_column_multi_elts_seq(OBIDMS_column_multi_elts):
         result = {}
         all_NA = True
         for i in range(self.nb_elements_per_line) :
-            value = obi_column_get_obiseq_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, i)
+            value = obi_column_get_obiseq_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, i)
             if obi_errno > 0 :
                 raise IndexError(line_nb)
             if strcmp(value, OBISeq_NA) == 0 :
@@ -83,6 +98,9 @@ cdef class OBIDMS_column_multi_elts_seq(OBIDMS_column_multi_elts):
             value_b = OBISeq_NA
         else :
             value_b = str2bytes(value)
-        if obi_column_set_obiseq_with_elt_name_in_view(self.view, (self.pointer)[0], line_nb, str2bytes(element_name), value_b) < 0:
+        if obi_column_set_obiseq_with_elt_name_in_view(self.view.pointer, (self.pointer)[0], line_nb, str2bytes(element_name), value_b) < 0:
             raise Exception("Problem setting a value in a column")
     
+#    cpdef align(self, ):        # TODO
+#        raise Exception("Columns with multiple sequences per line can't be aligned")    # TODO discuss
+

python/obitools3/obidms/_obidmscolumn_str.cfiles

@@ -10,6 +10,8 @@
 ../../../src/encode.c
 ../../../src/murmurhash2.h
 ../../../src/murmurhash2.c
+../../../src/obi_align.h
+../../../src/obi_align.c
 ../../../src/obiavl.h
 ../../../src/obiavl.c
 ../../../src/obiblob_indexer.h
@@ -37,6 +39,8 @@
 ../../../src/obiview.c
 ../../../src/utils.h
 ../../../src/utils.c
+../../../src/sse_banded_LCS_alignment.h
+../../../src/sse_banded_LCS_alignment.c
 ../../../src/obidmscolumn_bool.c
 ../../../src/obidmscolumn_bool.h
 ../../../src/obidmscolumn_char.c

python/obitools3/obidms/_obidmscolumn_str.pyx

@@ -17,7 +17,7 @@ cdef class OBIDMS_column_str(OBIDMS_column):
     cpdef object get_line(self, index_t line_nb):
         cdef const_char_p value
         cdef object       result
-        value = obi_column_get_obistr_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, 0)
+        value = obi_column_get_obistr_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, 0)
         if obi_errno > 0 :
             raise IndexError(line_nb)
         if strcmp(value, OBIStr_NA) == 0 :
@@ -33,7 +33,7 @@ cdef class OBIDMS_column_str(OBIDMS_column):
             value_b = OBIStr_NA
         else :
             value_b = str2bytes(value)
-        if obi_column_set_obistr_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, 0, value_b) < 0:
+        if obi_column_set_obistr_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, 0, value_b) < 0:
             raise Exception("Problem setting a value in a column")
 
 
@@ -42,7 +42,7 @@ cdef class OBIDMS_column_multi_elts_str(OBIDMS_column_multi_elts):
     cpdef object get_item(self, index_t line_nb, str element_name):
         cdef const_char_p value
         cdef object       result
-        value = obi_column_get_obistr_with_elt_name_in_view(self.view, (self.pointer)[0], line_nb, str2bytes(element_name))
+        value = obi_column_get_obistr_with_elt_name_in_view(self.view.pointer, (self.pointer)[0], line_nb, str2bytes(element_name))
         if obi_errno > 0 :
             raise IndexError(line_nb, element_name)
         if strcmp(value, OBIStr_NA) == 0 :
@@ -61,7 +61,7 @@ cdef class OBIDMS_column_multi_elts_str(OBIDMS_column_multi_elts):
         result = {}
         all_NA = True
         for i in range(self.nb_elements_per_line) :
-            value = obi_column_get_obistr_with_elt_idx_in_view(self.view, (self.pointer)[0], line_nb, i)
+            value = obi_column_get_obistr_with_elt_idx_in_view(self.view.pointer, (self.pointer)[0], line_nb, i)
             if obi_errno > 0 :
                 raise IndexError(line_nb)
             if strcmp(value, OBIStr_NA) == 0 :
@@ -82,6 +82,6 @@ cdef class OBIDMS_column_multi_elts_str(OBIDMS_column_multi_elts):
             value_b = OBIStr_NA
         else :
             value_b = str2bytes(value)
-        if obi_column_set_obistr_with_elt_name_in_view(self.view, (self.pointer)[0], line_nb, str2bytes(element_name), value_b) < 0:
+        if obi_column_set_obistr_with_elt_name_in_view(self.view.pointer, (self.pointer)[0], line_nb, str2bytes(element_name), value_b) < 0:
             raise Exception("Problem setting a value in a column")
     

python/obitools3/obidms/_obiseq.cfiles

@@ -10,6 +10,8 @@
 ../../../src/encode.c
 ../../../src/murmurhash2.h
 ../../../src/murmurhash2.c
+../../../src/obi_align.h
+../../../src/obi_align.c
 ../../../src/obiavl.h
 ../../../src/obiavl.c
 ../../../src/obiblob_indexer.h
@@ -37,6 +39,8 @@
 ../../../src/obiview.c
 ../../../src/utils.h
 ../../../src/utils.c
+../../../src/sse_banded_LCS_alignment.h
+../../../src/sse_banded_LCS_alignment.c
 ../../../src/obidmscolumn_bool.c
 ../../../src/obidmscolumn_bool.h
 ../../../src/obidmscolumn_char.c

python/obitools3/obidms/_obitaxo.cfiles

@@ -10,6 +10,8 @@
 ../../../src/encode.c
 ../../../src/murmurhash2.h
 ../../../src/murmurhash2.c
+../../../src/obi_align.h
+../../../src/obi_align.c
 ../../../src/obiavl.h
 ../../../src/obiavl.c
 ../../../src/obiblob_indexer.h
@@ -37,6 +39,8 @@
 ../../../src/obiview.c
 ../../../src/utils.h
 ../../../src/utils.c
+../../../src/sse_banded_LCS_alignment.h
+../../../src/sse_banded_LCS_alignment.c
 ../../../src/obidmscolumn_bool.c
 ../../../src/obidmscolumn_bool.h
 ../../../src/obidmscolumn_char.c

python/obitools3/obidms/capi/obialign.pxd

@@ -0,0 +1,10 @@
+#cython: language_level=3
+
+from ..capi.obiview      cimport Obiview_p
+from ..capi.obidmscolumn cimport OBIDMS_column_p
+
+
+cdef extern from "obi_align.h" nogil:
+
+    int obi_align_one_column(Obiview_p seq_view, OBIDMS_column_p seq_column, Obiview_p score_view, OBIDMS_column_p score_column, double threshold, bint normalize, int reference, bint similarity_mode)
+

src/_sse.h

@@ -0,0 +1,961 @@
+#ifndef _SSE_H_
+#define _SSE_H_
+
+#include <string.h>
+
+#include <inttypes.h>
+#ifdef __SSE2__
+#include <xmmintrin.h>
+#else
+typedef long long __m128i __attribute__ ((__vector_size__ (16), __may_alias__));
+#endif /* __SSE2__ */
+
+#ifndef MAX
+#define MAX(x,y) (((x)>(y)) ? (x):(y))
+#define MIN(x,y) (((x)<(y)) ? (x):(y))
+#endif
+
+#define ALIGN __attribute__((aligned(16)))
+typedef __m128i vUInt8;
+typedef __m128i vInt8;
+
+typedef __m128i vUInt16;
+typedef __m128i vInt16;
+
+typedef __m128i vUInt64;
+
+typedef union
+{
+     __m128i i;
+     int64_t  s64[ 2];
+     int16_t  s16[ 8];
+     int8_t   s8 [16];
+     uint8_t  u8 [16];
+     uint16_t u16[8 ];
+	 uint32_t u32[4 ];
+     uint64_t u64[2 ];
+} um128;
+
+typedef union
+    {
+	    vUInt8 m;
+	    uint8_t c[16];
+    } uchar_v;
+
+typedef union
+	{
+		vUInt16 m;
+		uint16_t c[8];
+	} ushort_v;
+
+typedef union
+	{
+		vUInt64 m;
+		uint64_t c[2];
+	} uint64_v;
+
+
+#ifdef __SSE2__
+
+static inline int8_t _s2_extract_epi8(__m128i r, const int p)
+{
+#define ACTIONP(r,x) return _mm_extract_epi16(r,x) & 0xFF
+#define ACTIONI(r,x) return _mm_extract_epi16(r,x) >> 8
+	switch (p) {
+	case  0:  ACTIONP(r,0);
+	case  1:  ACTIONI(r,0);
+	case  2:  ACTIONP(r,1);
+	case  3:  ACTIONI(r,1);
+	case  4:  ACTIONP(r,2);
+	case  5:  ACTIONI(r,2);
+	case  6:  ACTIONP(r,3);
+	case  7:  ACTIONI(r,3);
+	case  8:  ACTIONP(r,4);
+	case  9:  ACTIONI(r,4);
+	case 10:  ACTIONP(r,5);
+	case 11:  ACTIONI(r,5);
+	case 12:  ACTIONP(r,6);
+	case 13:  ACTIONI(r,6);
+	case 14:  ACTIONP(r,7);
+	case 15:  ACTIONI(r,7);
+	}
+#undef ACTIONP
+#undef ACTIONI
+
+	return 0;
+}
+
+static inline __m128i _s2_max_epi8(__m128i a, __m128i b)
+{
+	__m128i mask  = _mm_cmpgt_epi8( a, b );
+	a = _mm_and_si128   (a,mask );
+	b = _mm_andnot_si128(mask,b);
+	return _mm_or_si128(a,b);
+}
+
+static inline __m128i _s2_min_epi8(__m128i a, __m128i b)
+{
+	__m128i mask  = _mm_cmplt_epi8( a, b );
+	a = _mm_and_si128   (a,mask );
+	b = _mm_andnot_si128(mask,b);
+	return _mm_or_si128(a,b);
+}
+
+static inline __m128i _s2_insert_epi8(__m128i r, int b, const int p)
+{
+#define ACTIONP(r,x) return _mm_insert_epi16(r,(_mm_extract_epi16(r,x) & 0xFF00) | (b & 0x00FF),x)
+#define ACTIONI(r,x) return _mm_insert_epi16(r,(_mm_extract_epi16(r,x) & 0x00FF) | ((b << 8)& 0xFF00),x)
+	switch (p) {
+	case  0:  ACTIONP(r,0);
+	case  1:  ACTIONI(r,0);
+	case  2:  ACTIONP(r,1);
+	case  3:  ACTIONI(r,1);
+	case  4:  ACTIONP(r,2);
+	case  5:  ACTIONI(r,2);
+	case  6:  ACTIONP(r,3);
+	case  7:  ACTIONI(r,3);
+	case  8:  ACTIONP(r,4);
+	case  9:  ACTIONI(r,4);
+	case 10:  ACTIONP(r,5);
+	case 11:  ACTIONI(r,5);
+	case 12:  ACTIONP(r,6);
+	case 13:  ACTIONI(r,6);
+	case 14:  ACTIONP(r,7);
+	case 15:  ACTIONI(r,7);
+	}
+#undef ACTIONP
+#undef ACTIONI
+
+	return _mm_setzero_si128();
+}
+
+// Fill a SSE Register with 16 time the same 8bits integer value
+#define _MM_SET1_EPI8(x)        _mm_set1_epi8(x)
+#define _MM_INSERT_EPI8(r,x,i)	_s2_insert_epi8((r),(x),(i))
+#define _MM_CMPEQ_EPI8(x,y)     _mm_cmpeq_epi8((x),(y))
+#define _MM_CMPGT_EPI8(x,y)     _mm_cmpgt_epi8((x),(y))
+#define _MM_CMPLT_EPI8(x,y)     _mm_cmplt_epi8((x),(y))
+#define _MM_MAX_EPI8(x,y)       _s2_max_epi8((x),(y))
+#define _MM_MIN_EPI8(x,y)       _s2_min_epi8((x),(y))
+#define _MM_ADD_EPI8(x,y)       _mm_add_epi8((x),(y))
+#define _MM_SUB_EPI8(x,y)       _mm_sub_epi8((x),(y))
+#define _MM_EXTRACT_EPI8(r,p)   _s2_extract_epi8((r),(p))
+
+#define _MM_MIN_EPU8(x,y)       _mm_min_epu8((x),(y))
+
+// Fill a SSE Register with 8 time the same 16bits integer value
+#define _MM_SET1_EPI16(x)       _mm_set1_epi16(x)
+
+#define _MM_INSERT_EPI16(r,x,i)	_mm_insert_epi16((r),(x),(i))
+#define _MM_CMPEQ_EPI16(x,y)    _mm_cmpeq_epi16((x),(y))
+#define _MM_CMPGT_EPI16(x,y)    _mm_cmpgt_epi16((x),(y))
+#define _MM_CMPGT_EPU16(x,y)    _mm_cmpgt_epu16((x),(y))    // n'existe pas ??
+#define _MM_CMPLT_EPI16(x,y)    _mm_cmplt_epi16((x),(y))
+#define _MM_MAX_EPI16(x,y)      _mm_max_epi16((x),(y))
+#define _MM_MIN_EPI16(x,y)      _mm_min_epi16((x),(y))
+#define _MM_ADD_EPI16(x,y)      _mm_add_epi16((x),(y))
+#define _MM_SUB_EPI16(x,y)      _mm_sub_epi16((x),(y))
+#define _MM_EXTRACT_EPI16(r,p)  _mm_extract_epi16((r),(p))
+#define _MM_UNPACKLO_EPI8(a,b)  _mm_unpacklo_epi8((a),(b))
+#define _MM_UNPACKHI_EPI8(a,b)  _mm_unpackhi_epi8((a),(b))
+#define _MM_ADDS_EPU16(x,y)     _mm_adds_epu16((x),(y))
+
+// Multiplication
+#define _MM_MULLO_EPI16(x,y)    _mm_mullo_epi16((x), (y))
+
+#define _MM_SRLI_EPI64(r,x)     _mm_srli_epi64((r),(x))
+#define _MM_SLLI_EPI64(r,x)     _mm_slli_epi64((r),(x))
+
+// Set a SSE Register to 0
+#define _MM_SETZERO_SI128() _mm_setzero_si128()
+
+#define _MM_AND_SI128(x,y)     _mm_and_si128((x),(y))
+#define _MM_ANDNOT_SI128(x,y)  _mm_andnot_si128((x),(y))
+#define _MM_OR_SI128(x,y)      _mm_or_si128((x),(y))
+#define _MM_XOR_SI128(x,y)     _mm_xor_si128((x),(y))
+#define _MM_SLLI_SI128(r,s)	   _mm_slli_si128((r),(s))
+#define _MM_SRLI_SI128(r,s)    _mm_srli_si128((r),(s))
+
+// Load a SSE register from an unaligned address
+#define _MM_LOADU_SI128(x) _mm_loadu_si128(x)
+
+// Load a SSE register from an aligned address (/!\ not defined when SSE not available)
+#define _MM_LOAD_SI128(x)  _mm_load_si128(x)
+
+// #define	_MM_UNPACKLO_EPI8(x,y) _mm_unpacklo_epi8((x),(y))
+
+#else /* __SSE2__ Not defined */
+
+static inline __m128i _em_set1_epi8(int x)
+{
+	um128 a;
+
+	x&=0xFF;
+	a.s8[0]=x;
+	a.s8[1]=x;
+	a.u16[1]=a.u16[0];
+	a.u32[1]=a.u32[0];
+	a.u64[1]=a.u64[0];
+
+	return a.i;
+}
+
+static inline __m128i _em_insert_epi8(__m128i r, int x, const int i)
+{
+	um128 a;
+	a.i=r;
+	a.s8[i]=x & 0xFF;
+	return a.i;
+}
+
+static inline __m128i _em_cmpeq_epi8(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s8[z]=(x.s8[z]==y.s8[z]) ? 0xFF:0
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+	R(8);
+	R(9);
+	R(10);
+	R(11);
+	R(12);
+	R(13);
+	R(14);
+	R(15);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_cmpgt_epi8(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s8[z]=(x.s8[z]>y.s8[z]) ? 0xFF:0
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+	R(8);
+	R(9);
+	R(10);
+	R(11);
+	R(12);
+	R(13);
+	R(14);
+	R(15);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_cmplt_epi8(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s8[z]=(x.s8[z]<y.s8[z]) ? 0xFF:0
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+	R(8);
+	R(9);
+	R(10);
+	R(11);
+	R(12);
+	R(13);
+	R(14);
+	R(15);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_max_epi8(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s8[z]=MAX(x.s8[z],y.s8[z])
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+	R(8);
+	R(9);
+	R(10);
+	R(11);
+	R(12);
+	R(13);
+	R(14);
+	R(15);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_min_epi8(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s8[z]=MIN(x.s8[z],y.s8[z])
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+	R(8);
+	R(9);
+	R(10);
+	R(11);
+	R(12);
+	R(13);
+	R(14);
+	R(15);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_add_epi8(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s8[z]=x.s8[z]+y.s8[z]
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+	R(8);
+	R(9);
+	R(10);
+	R(11);
+	R(12);
+	R(13);
+	R(14);
+	R(15);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_sub_epi8(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s8[z]=x.s8[z]+y.s8[z]
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+	R(8);
+	R(9);
+	R(10);
+	R(11);
+	R(12);
+	R(13);
+	R(14);
+	R(15);
+#undef R
+
+	return r.i;
+}
+
+
+static inline int _em_extract_epi8(__m128i r, const int i)
+{
+	um128 a;
+
+	a.i=r;
+
+	return a.s8[i] & 0xFF;
+}
+
+static inline __m128i _em_min_epu8(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.u8[z]=MIN(x.u8[z],y.u8[z])
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+	R(8);
+	R(9);
+	R(10);
+	R(11);
+	R(12);
+	R(13);
+	R(14);
+	R(15);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_set1_epi16(int x)
+{
+	um128 a;
+
+	x&=0xFFFF;
+	a.s16[0]=x;
+	a.s16[1]=x;
+	a.u32[1]=a.u32[0];
+	a.u64[1]=a.u64[0];
+
+	return a.i;
+}
+
+static inline __m128i _em_insert_epi16(__m128i r, int x, const int i)
+{
+	um128 a;
+	a.i=r;
+	a.s16[i]=x & 0xFFFF;
+	return a.i;
+}
+
+static inline __m128i _em_cmpeq_epi16(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s16[z]=(x.s16[z]==y.s16[z]) ? 0xFFFF:0
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_cmpgt_epi16(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s16[z]=(x.s16[z]>y.s16[z]) ? 0xFFFF:0
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_cmplt_epi16(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s16[z]=(x.s16[z]<y.s16[z]) ? 0xFFFF:0
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_max_epi16(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+
+#define R(z) 	r.s16[z]=MAX(x.s16[z],y.s16[z])
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_min_epi16(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+
+#define R(z) 	r.s16[z]=MIN(x.s16[z],y.s16[z])
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_add_epi16(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s16[z]=x.s16[z]+y.s16[z]
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_sub_epi16(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s16[z]=x.s16[z]+y.s16[z]
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+#undef R
+
+	return r.i;
+}
+
+static inline int _em_extract_epi16(__m128i r, const int i)
+{
+	um128 a;
+	a.i=r;
+	return a.s16[i] & 0xFFFF;
+}
+
+static inline __m128i _em_unpacklo_epi8(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s16[z]=(((int16_t)(y.s8[z])) << 8) | (int16_t)(x.s8[z])
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_unpackhi_epi8(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.s16[z]=(((int16_t)(y.s8[z+8])) << 8) | (int16_t)(x.s8[z+8])
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_adds_epu16(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+
+#define R(z) 	r.u16[z]=x.u16[z]+y.u16[z]
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_srli_epi64(__m128i a, int b)
+{
+	um128 x;
+
+	x.i=a;
+
+	x.s64[0]>>=b;
+	x.s64[1]>>=b;
+
+	return x.i;
+}
+
+static inline __m128i _em_slli_epi64(__m128i a, int b)
+{
+	um128 x;
+
+	x.i=a;
+
+	x.s64[0]<<=b;
+	x.s64[1]<<=b;
+
+	return x.i;
+}
+
+static inline __m128i _em_setzero_si128()
+{
+	um128 x;
+
+	x.s64[0]=x.s64[1]=0;
+
+	return x.i;
+}
+
+static inline __m128i _em_and_si128(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+
+#define R(z) 	r.u64[z]=x.u64[z] & y.u64[z]
+	R(0);
+	R(1);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_andnot_si128(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+
+#define R(z) 	r.u64[z]=(~x.u64[z]) & y.u64[z]
+	R(0);
+	R(1);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_or_si128(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.u64[z]=x.u64[z] | y.u64[z]
+	R(0);
+	R(1);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_xor_si128(__m128i a, __m128i b)
+{
+	um128 x;
+	um128 y;
+	um128 r;
+
+	x.i=a;
+	y.i=b;
+
+#define R(z) 	r.u64[z]=x.u64[z] ^ y.u64[z]
+	R(0);
+	R(1);
+#undef R
+
+	return r.i;
+}
+
+static inline __m128i _em_slli_si128(__m128i a, int b)
+{
+	um128 x;
+
+	x.i=a;
+
+#define R(z) 	x.u8[z]=(z>=b) ? x.u8[z-b]:0
+	R(15);
+	R(14);
+	R(13);
+	R(12);
+	R(11);
+	R(10);
+	R(9);
+	R(8);
+	R(7);
+	R(6);
+	R(5);
+	R(4);
+	R(3);
+	R(2);
+	R(1);
+	R(0);
+#undef R
+
+	return x.i;
+}
+
+static inline __m128i _em_srli_si128(__m128i a, int b)
+{
+	um128 x;
+
+	x.i=a;
+
+#define R(z) 	x.u8[z]=((b+z) > 15) ? 0:x.u8[z+b]
+	R(0);
+	R(1);
+	R(2);
+	R(3);
+	R(4);
+	R(5);
+	R(6);
+	R(7);
+	R(8);
+	R(9);
+	R(10);
+	R(11);
+	R(12);
+	R(13);
+	R(14);
+	R(15);
+#undef R
+
+	return x.i;
+}
+
+inline static __m128i _em_loadu_si128(__m128i const *P)
+{
+	um128 tmp;
+	um128 *pp=(um128*)P;
+
+	tmp.u8[0]=(*pp).u8[0];
+    tmp.u8[1]=(*pp).u8[1];
+	tmp.u8[2]=(*pp).u8[2];
+    tmp.u8[3]=(*pp).u8[3];
+	tmp.u8[4]=(*pp).u8[4];
+    tmp.u8[5]=(*pp).u8[5];
+	tmp.u8[6]=(*pp).u8[6];
+    tmp.u8[7]=(*pp).u8[7];
+	tmp.u8[8]=(*pp).u8[8];
+    tmp.u8[9]=(*pp).u8[9];
+	tmp.u8[10]=(*pp).u8[10];
+    tmp.u8[11]=(*pp).u8[11];
+	tmp.u8[12]=(*pp).u8[12];
+    tmp.u8[13]=(*pp).u8[13];
+	tmp.u8[14]=(*pp).u8[14];
+    tmp.u8[15]=(*pp).u8[15];
+	return tmp.i;
+}
+
+
+#define _MM_SET1_EPI8(x)        _em_set1_epi8(x)
+#define _MM_INSERT_EPI8(r,x,i)	_em_insert_epi8((r),(x),(i))
+#define _MM_CMPEQ_EPI8(x,y)     _em_cmpeq_epi8((x),(y))
+#define _MM_CMPGT_EPI8(x,y)     _em_cmpgt_epi8((x),(y))
+#define _MM_CMPLT_EPI8(x,y)     _em_cmplt_epi8((x),(y))
+#define _MM_MAX_EPI8(x,y)       _em_max_epi8((x),(y))
+#define _MM_MIN_EPI8(x,y)       _em_min_epi8((x),(y))
+#define _MM_ADD_EPI8(x,y)       _em_add_epi8((x),(y))
+#define _MM_SUB_EPI8(x,y)       _em_sub_epi8((x),(y))
+#define _MM_EXTRACT_EPI8(r,p)   _em_extract_epi8((r),(p))
+
+#define _MM_MIN_EPU8(x,y)       _em_min_epu8((x),(y))
+
+#define _MM_SET1_EPI16(x)       _em_set1_epi16(x)
+#define _MM_INSERT_EPI16(r,x,i)	_em_insert_epi16((r),(x),(i))
+#define _MM_CMPEQ_EPI16(x,y)    _em_cmpeq_epi16((x),(y))
+#define _MM_CMPGT_EPI16(x,y)    _em_cmpgt_epi16((x),(y))
+#define _MM_CMPLT_EPI16(x,y)    _em_cmplt_epi16((x),(y))
+#define _MM_MAX_EPI16(x,y)      _em_max_epi16((x),(y))
+#define _MM_MIN_EPI16(x,y)      _em_min_epi16((x),(y))
+#define _MM_ADD_EPI16(x,y)      _em_add_epi16((x),(y))
+#define _MM_SUB_EPI16(x,y)      _em_sub_epi16((x),(y))
+#define _MM_EXTRACT_EPI16(r,p)  _em_extract_epi16((r),(p))
+#define _MM_UNPACKLO_EPI8(a,b)  _em_unpacklo_epi8((a),(b))
+#define _MM_UNPACKHI_EPI8(a,b)  _em_unpackhi_epi8((a),(b))
+#define _MM_ADDS_EPU16(x,y)     _em_adds_epu16((x),(y))
+
+#define _MM_SRLI_EPI64(r,x)     _em_srli_epi64((r),(x))
+#define _MM_SLLI_EPI64(r,x)     _em_slli_epi64((r),(x))
+
+#define _MM_SETZERO_SI128()     _em_setzero_si128()
+
+#define _MM_AND_SI128(x,y)      _em_and_si128((x),(y))
+#define _MM_ANDNOT_SI128(x,y)   _em_andnot_si128((x),(y))
+#define _MM_OR_SI128(x,y)       _em_or_si128((x),(y))
+#define _MM_XOR_SI128(x,y)      _em_xor_si128((x),(y))
+#define _MM_SLLI_SI128(r,s)	    _em_slli_si128((r),(s))
+#define _MM_SRLI_SI128(r,s)	    _em_srli_si128((r),(s))
+
+#define _MM_LOADU_SI128(x)      _em_loadu_si128(x)
+#define _MM_LOAD_SI128(x)       _em_loadu_si128(x)
+
+
+#endif /* __SSE2__ */
+
+#define _MM_NOT_SI128(x)     _MM_XOR_SI128((x),(_MM_SET1_EPI8(0xFFFF)))
+
+#endif

src/obi_align.c

@@ -0,0 +1,109 @@
+/****************************************************************************
+ * Sequence alignment functions				                                *
+ ****************************************************************************/
+
+/**
+ * @file obi_align.c
+ * @author Celine Mercier
+ * @date May 4th 2016
+ * @brief Functions handling sequence alignments.
+ */
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdbool.h>
+
+#include "obidebug.h"
+#include "obierrno.h"
+#include "obitypes.h"
+#include "obiview.h"
+#include "sse_banded_LCS_alignment.h"
+
+
+#define DEBUG_LEVEL 0	// TODO has to be defined somewhere else (cython compil flag?)
+
+
+// TODO
+// use openMP pragmas : garder scores en memoire et ecrire a la fin?
+// tout ecrire en stdint?
+// check NUC_SEQS and score type (int or float if normalize)
+// what's with multiple sequence/line columns?
+// make function that put blobs in int16
+
+
+int obi_align_one_column(Obiview_p seq_view, OBIDMS_column_p seq_column, Obiview_p score_view, OBIDMS_column_p score_column, double threshold, bool normalize, int reference, bool similarity_mode)
+{
+	index_t i, j, k;
+	index_t seq_count;
+	char* seq1;
+	char* seq2;
+	double score;
+
+	k = 0;
+
+	if ((seq_column->header)->returned_data_type != OBI_SEQ)
+	{
+		obi_set_errno(OBI_ALIGN_ERROR);
+		obidebug(1, "\nTrying to align a column of a different type than OBI_SEQ");
+		return -1;
+	}
+
+	if ((normalize && ((score_column->header)->returned_data_type != OBI_FLOAT)) ||
+			(!normalize && ((score_column->header)->returned_data_type != OBI_INT)))
+	{
+		obi_set_errno(OBI_ALIGN_ERROR);
+		obidebug(1, "\nTrying to store alignment scores in a column of an inappropriate type");
+		return -1;
+	}
+
+	seq_count = (seq_column->header)->lines_used;
+
+	for (i=0; i < (seq_count - 1); i++)
+	{
+		for (j=i+1; j < seq_count; j++)
+		{
+			//fprintf(stderr, "\ni=%lld, j=%lld, k=%lld", i, j, k);
+
+			seq1 = obi_column_get_obiseq_with_elt_idx_in_view(seq_view, seq_column, i, 0);
+			seq2 = obi_column_get_obiseq_with_elt_idx_in_view(seq_view, seq_column, j, 0);
+
+			if ((seq1 == NULL) || (seq2 == NULL))
+			{
+				obidebug(1, "\nError retrieving sequences to align");
+				return -1;
+			}
+
+			// TODO filter
+			score = generic_sse_banded_lcs_align(seq1, seq2, threshold, normalize, reference, similarity_mode);
+
+			if (normalize)
+			{
+				if (obi_column_set_obifloat_with_elt_idx_in_view(score_view, score_column, k, 0, (obifloat_t) score) < 0)
+				{
+					obidebug(1, "\nError writing alignment score in a column");
+					return -1;
+				}
+			}
+			else
+			{
+				if (obi_column_set_obiint_with_elt_idx_in_view(score_view, score_column, k, 0, (obiint_t) score) < 0)
+				{
+					obidebug(1, "\nError writing alignment score in a column");
+					return -1;
+				}
+			}
+
+			k++;
+		}
+	}
+
+	return 0;
+}
+
+
+int obi_align_two_columns(OBIDMS_column_p seq_column_1, OBIDMS_column_p seq_column_2, OBIDMS_column_p score_column, double threshold, bool normalize, int reference, bool similarity_mode)
+{
+	// TODO
+	return 0;
+}

src/obi_align.h

@@ -0,0 +1,37 @@
+/****************************************************************************
+ * Sequence alignment functions header file	                                *
+ ****************************************************************************/
+
+/**
+ * @file obi_align.h
+ * @author Celine Mercier
+ * @date May 11th 2016
+ * @brief Header file for the functions handling the alignment of DNA sequences.
+ */
+
+
+#ifndef OBI_ALIGN_H_
+#define OBI_ALIGN_H_
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdbool.h>
+
+#include "obidms.h"
+#include "obidmscolumn.h"
+#include "obitypes.h"
+
+
+
+/**
+ * @brief
+ *
+ * TODO
+ *
+ */
+int obi_align_one_column(Obiview_p seq_view, OBIDMS_column_p seq_column, Obiview_p score_view, OBIDMS_column_p score_column, double threshold, bool normalize, int reference, bool similarity_mode);
+
+
+#endif /* OBI_ALIGN_H_ */
+

src/obierrno.h

@@ -114,6 +114,8 @@ extern int obi_errno;
 	 	 	 	 	 	 	 	 	 	 	 */
 #define OBI_INDEXER_ERROR			(27)    /** Error handling a blob indexer
 											 */
+#define OBI_ALIGN_ERROR				(28)    /** Error while aligning sequences
+ 	 	 	 	 	 	 	 	 	 	 	 */
 /**@}*/
 
 #endif /* OBIERRNO_H_ */

src/sse_banded_LCS_alignment.c

@@ -0,0 +1,696 @@
+/*
+ * sse_banded_LCS_alignment.c
+ *
+ *  Created on: 7 nov. 2012
+ *      Author: celine mercier
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <stdint.h>
+#include <stdbool.h>
+
+#include "obierrno.h"
+#include "obidebug.h"
+#include "utils.h"
+#include "_sse.h"
+#include "sse_banded_LCS_alignment.h"
+
+
+#define DEBUG_LEVEL 0	// TODO has to be defined somewhere else (cython compil flag?)
+
+
+static void  printreg(__m128i r)
+{
+	int16_t a0,a1,a2,a3,a4,a5,a6,a7;
+
+	a0= _MM_EXTRACT_EPI16(r,0);
+	a1= _MM_EXTRACT_EPI16(r,1);
+	a2= _MM_EXTRACT_EPI16(r,2);
+	a3= _MM_EXTRACT_EPI16(r,3);
+	a4= _MM_EXTRACT_EPI16(r,4);
+	a5= _MM_EXTRACT_EPI16(r,5);
+	a6= _MM_EXTRACT_EPI16(r,6);
+	a7= _MM_EXTRACT_EPI16(r,7);
+
+fprintf(stderr, "a00 :-> %7d  %7d  %7d  %7d "
+		" %7d  %7d  %7d  %7d "
+		"\n"
+		, a0,a1,a2,a3,a4,a5,a6,a7
+		);
+}
+
+
+static inline int extract_reg(__m128i r, int p)
+{
+	switch (p) {
+	case 0: return(_MM_EXTRACT_EPI16(r,0));
+	case 1: return(_MM_EXTRACT_EPI16(r,1));
+	case 2: return(_MM_EXTRACT_EPI16(r,2));
+	case 3: return(_MM_EXTRACT_EPI16(r,3));
+	case 4: return(_MM_EXTRACT_EPI16(r,4));
+	case 5: return(_MM_EXTRACT_EPI16(r,5));
+	case 6: return(_MM_EXTRACT_EPI16(r,6));
+	case 7: return(_MM_EXTRACT_EPI16(r,7));
+	}
+	return(0);
+}
+
+
+// TODO warning on length order
+void sse_banded_align_lcs_and_ali_len(int16_t* seq1, int16_t* seq2, int l1, int l2, int bandLengthLeft, int bandLengthTotal, int16_t* address, double* lcs_length, int* ali_length)
+{
+	register int j;
+	int k1, k2;
+	int max, diff;
+	int l_reg, l_loc;
+	int line;
+	int numberOfRegistersPerLine;
+	int numberOfRegistersFor3Lines;
+
+	bool even_line;
+	bool odd_line;
+	bool even_BLL;
+	bool odd_BLL;
+
+	um128*  SSEregisters;
+	um128*  p_diag;
+	um128*  p_gap1;
+	um128*  p_gap2;
+	um128*  p_diag_j;
+	um128*  p_gap1_j;
+	um128*  p_gap2_j;
+	um128   current;
+
+	um128*  l_ali_SSEregisters;
+	um128*  p_l_ali_diag;
+	um128*  p_l_ali_gap1;
+	um128*  p_l_ali_gap2;
+	um128*  p_l_ali_diag_j;
+	um128*  p_l_ali_gap1_j;
+	um128*  p_l_ali_gap2_j;
+	um128   l_ali_current;
+
+	um128  nucs1;
+	um128  nucs2;
+	um128  scores;
+
+	um128 boolean_reg;
+
+	// Initialisations
+
+	odd_BLL = bandLengthLeft & 1;
+	even_BLL  = !odd_BLL;
+
+	max = INT16_MAX - l1;
+
+	numberOfRegistersPerLine = bandLengthTotal / 8;
+	numberOfRegistersFor3Lines   = 3 * numberOfRegistersPerLine;
+
+	SSEregisters = (um128*) calloc(numberOfRegistersFor3Lines * 2, sizeof(um128));
+	l_ali_SSEregisters = SSEregisters + numberOfRegistersFor3Lines;
+
+	// preparer registres SSE
+
+	for (j=0; j<numberOfRegistersFor3Lines; j++)
+		l_ali_SSEregisters[j].i = _MM_LOAD_SI128(address+j*8);
+
+	p_diag    = SSEregisters;
+	p_gap1    = SSEregisters+numberOfRegistersPerLine;
+	p_gap2    = SSEregisters+2*numberOfRegistersPerLine;
+
+	p_l_ali_diag    = l_ali_SSEregisters;
+	p_l_ali_gap1    = l_ali_SSEregisters+numberOfRegistersPerLine;
+	p_l_ali_gap2    = l_ali_SSEregisters+2*numberOfRegistersPerLine;
+
+	// Loop on diagonals = 'lines' :
+	for (line=2; line <= l1+l2; line++)
+	{
+		odd_line = line & 1;
+		even_line  = !odd_line;
+
+		// loop on the registers of a line :
+		for (j=0; j < numberOfRegistersPerLine; j++)
+		{
+			p_diag_j       = p_diag+j;
+			p_gap1_j       = p_gap1+j;
+			p_gap2_j       = p_gap2+j;
+			p_l_ali_diag_j = p_l_ali_diag+j;
+			p_l_ali_gap1_j = p_l_ali_gap1+j;
+			p_l_ali_gap2_j = p_l_ali_gap2+j;
+
+		// comparing nucleotides for diagonal scores :
+
+			// k1 = position of the 1st nucleotide to align for seq1 and k2 = position of the 1st nucleotide to align for seq2
+			if (odd_line && odd_BLL)
+				k1 = (line / 2) + ((bandLengthLeft+1) / 2) - j*8;
+			else
+				k1 = (line / 2) + (bandLengthLeft/2) - j*8;
+
+			k2 = line - k1 - 1;
+
+			nucs1.i = _MM_LOADU_SI128(seq1+l1-k1);
+			nucs2.i = _MM_LOADU_SI128(seq2+k2);
+
+/*			fprintf(stderr, "\nnucs, r %d, k1 = %d, k2 = %d\n", j, k1, k2);
+			printreg(nucs1.i);
+			printreg(nucs2.i);
+*/
+
+		// computing diagonal score :
+			scores.i = _MM_AND_SI128(_MM_CMPEQ_EPI16(nucs1.i, nucs2.i), _MM_SET1_EPI16(1));
+			current.i = _MM_ADDS_EPU16(p_diag_j->i, scores.i);
+
+		// Computing alignment length
+
+			l_ali_current.i = p_l_ali_diag_j->i;
+			boolean_reg.i = _MM_CMPGT_EPI16(p_gap1_j->i, current.i);
+			l_ali_current.i = _MM_OR_SI128(
+								_MM_AND_SI128(p_l_ali_gap1_j->i, boolean_reg.i),
+								_MM_ANDNOT_SI128(boolean_reg.i, l_ali_current.i));
+			current.i = _MM_OR_SI128(
+							_MM_AND_SI128(p_gap1_j->i, boolean_reg.i),
+							_MM_ANDNOT_SI128(boolean_reg.i, current.i));
+			boolean_reg.i = _MM_AND_SI128(
+								_MM_CMPEQ_EPI16(p_gap1_j->i, current.i),
+								_MM_CMPLT_EPI16(p_l_ali_gap1_j->i, l_ali_current.i));
+			l_ali_current.i = _MM_OR_SI128(
+								_MM_AND_SI128(p_l_ali_gap1_j->i, boolean_reg.i),
+								_MM_ANDNOT_SI128(boolean_reg.i, l_ali_current.i));
+			current.i = _MM_OR_SI128(
+							_MM_AND_SI128(p_gap1_j->i, boolean_reg.i),
+							_MM_ANDNOT_SI128(boolean_reg.i, current.i));
+			boolean_reg.i = _MM_CMPGT_EPI16(p_gap2_j->i, current.i);
+			l_ali_current.i = _MM_OR_SI128(
+								_MM_AND_SI128(p_l_ali_gap2_j->i, boolean_reg.i),
+								_MM_ANDNOT_SI128(boolean_reg.i, l_ali_current.i));
+			current.i = _MM_OR_SI128(
+							_MM_AND_SI128(p_gap2_j->i, boolean_reg.i),
+							_MM_ANDNOT_SI128(boolean_reg.i, current.i));
+			boolean_reg.i = _MM_AND_SI128(
+								_MM_CMPEQ_EPI16(p_gap2_j->i, current.i),
+								_MM_CMPLT_EPI16(p_l_ali_gap2_j->i, l_ali_current.i));
+			l_ali_current.i = _MM_OR_SI128(
+								_MM_AND_SI128(p_l_ali_gap2_j->i, boolean_reg.i),
+								_MM_ANDNOT_SI128(boolean_reg.i, l_ali_current.i));
+			current.i = _MM_OR_SI128(
+							_MM_AND_SI128(p_gap2_j->i, boolean_reg.i),
+							_MM_ANDNOT_SI128(boolean_reg.i, current.i));
+
+/*
+			fprintf(stderr, "\nline = %d", line);
+			fprintf(stderr, "\nDiag, r %d : ", j);
+			printreg((*(p_diag_j)).i);
+			fprintf(stderr, "Gap1      : ");
+			printreg((*(p_gap1_j)).i);
+			fprintf(stderr, "Gap2      : ");
+			printreg((*(p_gap2_j)).i);
+			fprintf(stderr, "current   : ");
+			printreg(current.i);
+			fprintf(stderr, "L ALI\nDiag  r %d : ", j);
+			printreg((*(p_l_ali_diag_j)).i);
+			fprintf(stderr, "Gap1      : ");
+			printreg((*(p_l_ali_gap1_j)).i);
+			fprintf(stderr, "Gap2      : ");
+			printreg((*(p_l_ali_gap2_j)).i);
+			fprintf(stderr, "current   : ");
+			printreg(l_ali_current.i);
+*/
+
+		// diag = gap1 and gap1 = current
+			p_diag_j->i = p_gap1_j->i;
+			p_gap1_j->i = current.i;
+
+		// l_ali_diag = l_ali_gap1 and l_ali_gap1 = l_ali_current+1
+			p_l_ali_diag_j->i = p_l_ali_gap1_j->i;
+			p_l_ali_gap1_j->i = _MM_ADD_EPI16(l_ali_current.i, _MM_SET1_EPI16(1));
+		}
+
+		// shifts for gap2, to do only once all the registers of a line have been computed     Copier gap2 puis le charger depuis la copie?
+
+		for (j=0; j < numberOfRegistersPerLine; j++)
+		{
+			if ((odd_line && even_BLL) || (even_line && odd_BLL))
+			{
+				p_gap2[j].i       = _MM_LOADU_SI128((p_gap1[j].s16)-1);
+				p_l_ali_gap2[j].i = _MM_LOADU_SI128((p_l_ali_gap1[j].s16)-1);
+				if (j == 0)
+				{
+					p_gap2[j].i = _MM_INSERT_EPI16(p_gap2[j].i, 0, 0);
+					p_l_ali_gap2[j].i = _MM_INSERT_EPI16(p_l_ali_gap2[j].i, max, 0);
+				}
+			}
+			else
+			{
+				p_gap2[j].i       = _MM_LOADU_SI128(p_gap1[j].s16+1);
+				p_l_ali_gap2[j].i = _MM_LOADU_SI128(p_l_ali_gap1[j].s16+1);
+				if (j == numberOfRegistersPerLine - 1)
+				{
+					p_gap2[j].i = _MM_INSERT_EPI16(p_gap2[j].i, 0, 7);
+					p_l_ali_gap2[j].i = _MM_INSERT_EPI16(p_l_ali_gap2[j].i, max, 7);
+				}
+			}
+		}
+		// end shifts for gap2
+
+	}
+
+/*  /// Recovering LCS and alignment lengths  \\\  */
+
+	// finding the location of the results in the registers :
+	diff = l1-l2;
+	if ((diff & 1) && odd_BLL)
+		l_loc = (int) floor((double)(bandLengthLeft) / (double)2) - floor((double)(diff) / (double)2);
+	else
+		l_loc = (int) floor((double)(bandLengthLeft) / (double)2) - ceil((double)(diff) / (double)2);
+
+	l_reg = (int)floor((double)l_loc/(double)8.0);
+	//fprintf(stderr, "\nl_reg = %d, l_loc = %d\n", l_reg, l_loc);
+	l_loc = l_loc - l_reg*8;
+
+	// extracting the results from the registers :
+	*lcs_length = extract_reg(p_gap1[l_reg].i, l_loc);
+	*ali_length = extract_reg(p_l_ali_gap1[l_reg].i, l_loc) - 1;
+
+	// freeing the registers
+	free(SSEregisters);
+}
+
+
+// TODO warning on length order
+double sse_banded_align_just_lcs(int16_t* seq1, int16_t* seq2, int l1, int l2, int bandLengthLeft, int bandLengthTotal)
+{
+	register int j;
+	int k1, k2;
+	int diff;
+	int l_reg, l_loc;
+	int16_t l_lcs;
+	int line;
+	int numberOfRegistersPerLine;
+	int numberOfRegistersFor3Lines;
+
+	bool even_line;
+	bool odd_line;
+	bool even_BLL;
+	bool odd_BLL;
+
+	um128*  SSEregisters;
+	um128*  p_diag;
+	um128*  p_gap1;
+	um128*  p_gap2;
+	um128*  p_diag_j;
+	um128*  p_gap1_j;
+	um128*  p_gap2_j;
+	um128   current;
+
+	um128  nucs1;
+	um128  nucs2;
+	um128  scores;
+
+	// Initialisations
+
+	odd_BLL = bandLengthLeft & 1;
+	even_BLL  = !odd_BLL;
+
+	numberOfRegistersPerLine = bandLengthTotal / 8;
+	numberOfRegistersFor3Lines   = 3 * numberOfRegistersPerLine;
+
+	SSEregisters = malloc(numberOfRegistersFor3Lines * sizeof(um128));
+	if (SSEregisters == NULL)
+	{
+		obi_set_errno(OBI_MALLOC_ERROR);
+		obidebug(1, "\nError allocating memory for SSE registers for LCS alignment");
+		return 0;	// TODO DOUBLE_MIN?
+	}
+
+	// preparer registres SSE
+
+	for (j=0; j<numberOfRegistersFor3Lines; j++)
+		(*(SSEregisters+j)).i       = _MM_SETZERO_SI128();
+
+	p_diag    = SSEregisters;
+	p_gap1    = SSEregisters+numberOfRegistersPerLine;
+	p_gap2    = SSEregisters+2*numberOfRegistersPerLine;
+
+	// Loop on diagonals = 'lines' :
+	for (line=2; line <= l1+l2; line++)
+	{
+		odd_line = line & 1;
+		even_line  = !odd_line;
+
+		// loop on the registers of a line :
+		for (j=0; j < numberOfRegistersPerLine; j++)
+		{
+			p_diag_j = p_diag+j;
+			p_gap1_j = p_gap1+j;
+			p_gap2_j = p_gap2+j;
+
+		// comparing nucleotides for diagonal scores :
+
+			// k1 = position of the 1st nucleotide to align for seq1 and k2 = position of the 1st nucleotide to align for seq2
+			if (odd_line && odd_BLL)
+				k1 = (line / 2) + ((bandLengthLeft+1) / 2) - j*8;
+			else
+				k1 = (line / 2) + (bandLengthLeft/2) - j*8;
+
+			k2 = line - k1 - 1;
+
+			nucs1.i = _MM_LOADU_SI128(seq1+l1-k1);
+			nucs2.i = _MM_LOADU_SI128(seq2+k2);
+
+		// computing diagonal score :
+			scores.i = _MM_AND_SI128(_MM_CMPEQ_EPI16(nucs1.i, nucs2.i), _MM_SET1_EPI16(1));
+			current.i = _MM_ADDS_EPU16((*(p_diag_j)).i, scores.i);
+
+		// current = max(gap1, current)
+			current.i = _MM_MAX_EPI16((*(p_gap1_j)).i, current.i);
+
+		// current  = max(gap2, current)
+			current.i = _MM_MAX_EPI16((*(p_gap2_j)).i, current.i);
+
+		// diag = gap1 and gap1 = current
+			(*(p_diag_j)).i = (*(p_gap1_j)).i;
+			(*(p_gap1_j)).i = current.i;
+		}
+
+		// shifts for gap2, to do only once all the registers of a line have been computed
+
+			for (j=0; j < numberOfRegistersPerLine; j++)
+			{
+				if ((odd_line && even_BLL) || (even_line && odd_BLL))
+				{
+					(*(p_gap2+j)).i = _MM_LOADU_SI128(((*(p_gap1+j)).s16)-1);
+					if (j == 0)
+					{
+						(*(p_gap2+j)).i = _MM_INSERT_EPI16((*(p_gap2+j)).i, 0, 0);
+					}
+				}
+				else
+				{
+					(*(p_gap2+j)).i = _MM_LOADU_SI128(((*(p_gap1+j)).s16)+1);
+					if (j == numberOfRegistersPerLine - 1)
+					{
+						(*(p_gap2+j)).i = _MM_INSERT_EPI16((*(p_gap2+j)).i, 0, 7);
+					}
+				}
+			}
+		// end shifts for gap2
+	}
+
+/*  /// Recovering LCS and alignment lengths  \\\  */
+
+	// finding the location of the results in the registers :
+	diff = l1-l2;
+	if ((diff & 1) && odd_BLL)
+		l_loc = (int) floor((double)(bandLengthLeft) / (double)2) - floor((double)(diff) / (double)2);
+	else
+		l_loc = (int) floor((double)(bandLengthLeft) / (double)2) - ceil((double)(diff) / (double)2);
+
+	l_reg = (int)floor((double)l_loc/(double)8.0);
+	//fprintf(stderr, "\nl_reg = %d, l_loc = %d\n", l_reg, l_loc);
+	l_loc = l_loc - l_reg*8;
+
+	// extracting LCS from the registers :
+	l_lcs = extract_reg((*(p_gap1+l_reg)).i, l_loc);
+
+	// freeing the registers
+	free(SSEregisters);
+
+	return((double) l_lcs);
+}
+
+
+int calculateLeftBandLength(int lmax, int LCSmin)
+{
+	return (lmax - LCSmin);
+}
+
+
+int calculateRightBandLength(int lmin, int LCSmin)
+{
+	return (lmin - LCSmin);
+}
+
+
+int calculateSSEBandLength(int bandLengthRight, int bandLengthLeft)
+{
+	int bandLengthTotal= (double)(bandLengthRight + bandLengthLeft) / 2.0 + 1.0;
+
+	return (bandLengthTotal & (~ (int)7)) + (( bandLengthTotal & (int)7) ? 8:0); // Calcule le multiple de 8 superieur
+}
+
+
+// TODO that's gonna be fun to doc
+int calculateSizeToAllocate(int maxLen, int minLen, int LCSmin)
+{
+	int size;
+
+	size = calculateLeftBandLength(maxLen, LCSmin);
+
+	size *=  2;
+	size  =  (size & (~ (int)7)) + (( size & (int)7) ? 8:0); // Closest greater 8 multiple
+	size *=  3;
+	size +=  16;
+
+	return(size*sizeof(int16_t));
+}
+
+
+void iniSeq(int16_t* seq, int size, int16_t iniValue)
+{
+	int16_t* target = seq;
+	int16_t* end = target + (size_t)size;
+
+	for (; target < end; target++)
+		*target = iniValue;
+}
+
+
+void putSeqInSeq(int16_t* seq, char* s, int l, bool reverse)
+{
+	int16_t *target=seq;
+	int16_t *end = target + (size_t)l;
+	char    *source=s;
+
+	if (reverse)
+		for (source=s + (size_t)l-1; target < end; target++, source--)
+			*target=*source;
+	else
+		for (; target < end; source++,target++)
+			*target=*source;
+}
+
+
+void initializeAddressWithGaps(int16_t* address, int bandLengthTotal, int bandLengthLeft, int l1)
+{
+	int i;
+	int address_00, x_address_10, address_01, address_01_shifted;
+	int numberOfRegistersPerLine;
+	int bm;
+	int value=INT16_MAX-l1;
+
+	numberOfRegistersPerLine   = bandLengthTotal / 8;
+	bm = bandLengthLeft%2;
+
+	for (i=0; i < (3*numberOfRegistersPerLine*8); i++)
+		address[i] = value;
+
+
+	// 0,0 set to 1 and 0,1 and 1,0 set to 2
+
+	address_00   = bandLengthLeft / 2;
+
+	x_address_10 = address_00 + bm - 1;
+	address_01   = numberOfRegistersPerLine*8 + x_address_10;
+
+	address_01_shifted = numberOfRegistersPerLine*16 + address_00 - bm;
+
+	// fill address_00, address_01,+1, address_01_shifted,+1
+
+	address[address_00] = 1;
+	address[address_01] = 2;
+	address[address_01+1] = 2;
+	address[address_01_shifted] = 2;
+	address[address_01_shifted+1] = 2;
+}
+
+
+// TODO warning on length order
+double sse_banded_lcs_align(int16_t* seq1, int16_t* seq2, int l1, int l2, bool normalize, int reference, bool similarity_mode, int16_t* address, int LCSmin)
+{
+	double id;
+	int bandLengthRight, bandLengthLeft, bandLengthTotal;
+	int ali_length;
+
+	bandLengthLeft = calculateLeftBandLength(l1, LCSmin);
+	bandLengthRight = calculateRightBandLength(l2, LCSmin);
+
+//	fprintf(stderr, "\nBLL = %d, BLR = %d, LCSmin = %d\n", bandLengthLeft, bandLengthRight, LCSmin);
+
+	bandLengthTotal = calculateSSEBandLength(bandLengthRight, bandLengthLeft);
+
+//	fprintf(stderr, "\nBLT = %d\n", bandLengthTotal);
+
+	if ((reference == ALILEN) && (normalize || !similarity_mode))
+	{
+		initializeAddressWithGaps(address, bandLengthTotal, bandLengthLeft, l1);
+		sse_banded_align_lcs_and_ali_len(seq1, seq2, l1, l2, bandLengthLeft, bandLengthTotal, address, &id, &ali_length);
+	}
+	else
+		id = sse_banded_align_just_lcs(seq1, seq2, l1, l2, bandLengthLeft, bandLengthTotal);
+
+//	fprintf(stderr, "\nid before normalizations = %f", id);
+
+//	fprintf(stderr, "\nlcs = %f, ali = %d\n", id, ali_length);
+
+	if (!similarity_mode && !normalize)
+		switch(reference) {
+			case ALILEN: id = ali_length - id;
+						break;
+			case MAXLEN: id = l1 - id;
+						break;
+			case MINLEN: id = l2 - id;
+		}
+
+//	fprintf(stderr, "\n2>>> %f, %d\n", id, ali_length);
+	if (normalize)
+		switch(reference) {
+			case ALILEN: id = id / (double) ali_length;
+						break;
+			case MAXLEN: id = id / (double) l1;
+						break;
+			case MINLEN: id = id / (double) l2;
+		}
+
+//	fprintf(stderr, "\nid = %f\n", id);
+	return(id);
+}
+
+
+
+// PUBLIC FUNCTIONS
+
+
+int calculateLCSmin(int l1, int l2, double threshold, bool normalize, int reference, bool similarity_mode)
+{
+	int LCSmin;
+
+	if (threshold > 0)
+	{
+		if (normalize)
+		{
+			if (reference == MINLEN)
+				LCSmin = threshold*l2;
+			else 		// ref = maxlen or alilen
+				LCSmin = threshold*l1;
+		}
+		else if (similarity_mode)
+			LCSmin = threshold;
+		else if (reference == MINLEN) // not similarity_mode
+			LCSmin = l2 - threshold;
+		else	// not similarity_mode and ref = maxlen or alilen
+			LCSmin = l1 - threshold;
+	}
+	else
+		LCSmin = 0;
+
+	return(LCSmin);
+}
+
+
+double generic_sse_banded_lcs_align(char* seq1, char* seq2, double threshold, bool normalize, int reference, bool similarity_mode)
+{
+	double   id;
+	int      l1, l2;
+	int      lmax, lmin;
+	int      sizeToAllocateForBand, sizeToAllocateForSeqs;
+	int	     maxBLL;
+	int      LCSmin;
+	int      shift;
+	int16_t* address;
+	int16_t* iseq1;
+	int16_t* iseq2;
+
+	address = NULL;
+
+	l1 = strlen(seq1);
+	l2 = strlen(seq2);
+
+	if (l1 > l2)
+	{
+		lmax = l1;
+		lmin = l2;
+	}
+	else
+	{
+		lmax = l2;
+		lmin = l1;
+	}
+
+	// If the score is expressed as a normalized distance, get the corresponding identity
+	if (!similarity_mode && normalize)
+		threshold = 1.0 - threshold;
+
+	// Calculate the minimum LCS length corresponding to the threshold
+	LCSmin = calculateLCSmin(lmax, lmin, threshold, normalize, reference, similarity_mode);
+
+	// Allocate space for matrix band if the alignment length must be computed
+	if ((reference == ALILEN) && (normalize || !similarity_mode)) // cases in which alignment length must be computed
+	{
+		sizeToAllocateForBand = calculateSizeToAllocate(lmax, lmin, LCSmin);
+		address = obi_get_memory_aligned_on_16(sizeToAllocateForBand, &shift);
+		if (address == NULL)
+		{
+			obi_set_errno(OBI_MALLOC_ERROR);
+			obidebug(1, "\nError getting a memory address aligned on 16 bytes boundary");
+			return 0;	// TODO DOUBLE_MIN
+		}
+	}
+
+	// Allocate space for the int16_t arrays representing the sequences
+	maxBLL = calculateLeftBandLength(lmax, LCSmin);
+	sizeToAllocateForSeqs = 2*maxBLL+lmax;
+	iseq1 = (int16_t*) malloc(sizeToAllocateForSeqs*sizeof(int16_t));
+	iseq2 = (int16_t*) malloc(sizeToAllocateForSeqs*sizeof(int16_t));
+	if ((iseq1 == NULL) || (iseq2 == NULL))
+	{
+		obi_set_errno(OBI_MALLOC_ERROR);
+		obidebug(1, "\nError allocating memory for integer arrays to use in LCS alignment");
+		return 0; 	// TODO DOUBLE_MIN
+	}
+
+	// Initialize the int arrays
+	iniSeq(iseq1, (2*maxBLL)+lmax, 0);
+	iniSeq(iseq2, (2*maxBLL)+lmax, 255);
+
+	// Shift addresses to where the sequences have to be put
+	iseq1 = iseq1+maxBLL;
+	iseq2 = iseq2+maxBLL;
+
+	// Put the DNA sequences in the int arrays. Longest sequence must be first argument of sse_align function
+	if (l2 > l1)
+	{
+		putSeqInSeq(iseq1, seq2, l2, TRUE);
+		putSeqInSeq(iseq2, seq1, l1, FALSE);
+		// Compute alignment
+		id = sse_banded_lcs_align(iseq1, iseq2, l2, l1, normalize, reference, similarity_mode, address, LCSmin);
+	}
+	else
+	{
+		putSeqInSeq(iseq1, seq1, l1, TRUE);
+		putSeqInSeq(iseq2, seq2, l2, FALSE);
+		// Compute alignment
+		id = sse_banded_lcs_align(iseq1, iseq2, l1, l2, normalize, reference, similarity_mode, address, LCSmin);
+	}
+
+	// Free allocated elements
+	if (address != NULL)
+		free(address-shift);
+	free(iseq1-maxBLL);
+	free(iseq2-maxBLL);
+
+	return(id);
+}
+

src/sse_banded_LCS_alignment.h

@@ -0,0 +1,23 @@
+/*
+ * sse_banded_LCS_alignment.h
+ *
+ *  Created on: november 29, 2012
+ *      Author: mercier
+ */
+
+#ifndef SSE_BANDED_LCS_ALIGNMENT_H_
+#define SSE_BANDED_LCS_ALIGNMENT_H_
+
+
+#include <stdint.h>
+#include <stdbool.h>
+
+#define ALILEN (0)	// TODO enum
+#define MAXLEN (1)
+#define MINLEN (2)
+
+// TODO doc
+int calculateLCSmin(int l1, int l2, double threshold, bool normalize, int reference, bool lcsmode);
+double generic_sse_banded_lcs_align(char* seq1, char* seq2, double threshold, bool normalize, int reference, bool similarity_mode);
+
+#endif

src/utils.c

@@ -65,6 +65,12 @@ char* obi_format_date(time_t date)
 	struct tm*  tmp;
 
 	formatted_time = (char*) malloc(FORMATTED_TIME_LENGTH*sizeof(char));
+	if (formatted_time == NULL)
+	{
+		obi_set_errno(OBI_MALLOC_ERROR);
+		obidebug(1, "\nError allocating memory to format a date");
+		return NULL;
+	}
 
 	tmp = localtime(&date);
 	if (tmp == NULL)
@@ -84,3 +90,27 @@ char* obi_format_date(time_t date)
 	return formatted_time;
 }
 
+
+void* obi_get_memory_aligned_on_16(int size, int* shift)
+{
+	void* memory;
+
+	*shift = 0;
+
+	memory = (void*) malloc(size);
+	if (memory == NULL)
+	{
+		obi_set_errno(OBI_MALLOC_ERROR);
+		obidebug(1, "\nError allocating memory");
+		return NULL;
+	}
+
+	while ((((long long unsigned int) (memory))%16) != 0)
+	{
+		memory++;
+		(*shift)++;
+	}
+
+	return (memory);
+}
+

src/utils.h

@@ -53,4 +53,25 @@ int count_dir(char* dir_path);
 char* obi_format_date(time_t date);
 
 
+/**
+ * @brief Allocates a chunk of memory aligned on 16 bytes boundary.
+ *
+ * @warning The pointer returned must be freed by the caller.
+ * @warning The memory chunk pointed at by the returned pointer can be
+ *          smaller than the size asked for, since the address is shifted
+ *          to be aligned.
+ *
+ * @param size The size in bytes of the memory chunk to be allocated.
+ * 			   Ideally the closest multiple of 16 greater than the wanted size.
+ * @param shift A pointer on an integer corresponding to the shift made to align
+ * 				the address. Used to free the memory chunk.
+ *
+ * @returns The pointer on the aligned memory.
+ *
+ * @since May 2016
+ * @author Celine Mercier (celine.mercier@metabarcoding.org)
+ */
+void* obi_get_memory_aligned_on_16(int size, int* shift);
+
+
 #endif /* UTILS_H_ */