work in progress for mapping read extracts by chunk on GPU

_modules/GPU_chunkMapper.py

@@ -35,6 +35,11 @@ class GPU_chunkMapper(gpuBasicMapper):
         max_nb_occ = getMaxNbOccPerSeq(o_nb_occ_foreach_extract_in_seqs, chunk_size, 1)
         ficname1 = "read_occ_positions" + str(idx_seq) + "_" + str(chunk_idx)
         ficname2 = "reduced" + str(idx_seq) + "_" + str(chunk_idx)
+        ficname3 = "max_occ" + str(idx_seq) + "_" + str(chunk_idx)
+        fname=os.path.join(self.ficdir, ficname3)
+        of = open(fname,"w")
+        of.write(str(max_nb_occ))
+        of.close()
         if (max_nb_occ>0):
             host_read_occ_positions = self.getOccPositions(max_nb_occ, 1, chunk_size)
             # Gets reduced version of output_reads_nb_occ_in_seqs that enables me to iterate only on useful elements
@@ -66,9 +71,6 @@ class GPU_chunkMapper(gpuBasicMapper):
         for seq in self.seq_data.refseq_list:
             h = nb_extracts // self.nb_chunks # TODO: change code here to take nb_chunks into account.
             extracts_chunk = extracts_l[0:h+1] # Problem here: element at index h is not copied
-            # print extracts_chunk
-            # print len(extracts_chunk)
-            # print "h=",h," extcats_chunk[h]=",extracts_chunk[h]
             self.mapOneChunk(extracts_chunk,h+1,0, seq,seq_sizes_l[idx_seq], idx_seq)
             extracts_chunk=extracts_l[h+1:]
             self.mapOneChunk(extracts_chunk,nb_extracts-h-1, 1, seq, seq_sizes_l[idx_seq], idx_seq)

_modules/GPU_mapper.py

@@ -36,6 +36,18 @@ class gpuBasicMappingResultsForCov:
             self.rev_cpl_PE2_start_match=l_matches[3]
             self.rev_cpl_PE2_stop_match=l_matches[4]
 
+def getGPUMappingResults(reduced_fic,occ_pos_fic,max_nb_occ_fic,nb_reads,nb_extracts,ref_data,single=True):
+    reduced=np.loadtxt(reduced_fic, dtype=np.int64, delimiter='|')
+    occ_pos=np.loadtxt(occ_pos_fic, dtype=np.int32, delimiter='|')
+    f_nb_occ = open(max_nb_occ_fic)
+    str_nb_occ = f_nb_occ.read()
+    f_nb_occ.close()
+    max_nb_occ = int(str_nb_occ)
+    mres=gpuBasicMappingResults(reduced, occ_pos, max_nb_occ, nb_reads,
+                           nb_extracts, ref_data,single)
+    return mres
+
+
 ## Class to encapsulate the results of the mapping of all read extracts against all sequences.
 class gpuBasicMappingResults:
     ## Constructor
@@ -142,7 +154,7 @@ class gpuBasicMapper:
     def __init__(self,read_extracts,ref_data):
         self.read_extracts=read_extracts
         self.seq_data=ref_data
-        self.starting_from=0
+        #self.starting_from=0
 
     ## Prepare input parameters on device side, only available for parameters that are used by several kernels.
     def prepInputOnDevice(self,input_rextracts,nb_sequences,nb_extracts):

_modules/functions_PhageTerm_gpu.py

@@ -8,6 +8,10 @@ import numpy as np
 from IData_handling import ReadGetter
 from utilities import reverseComplement,hybridCoverage,correctEdge
 
+##
+# This function aims at processing the results of the mapping performed on GPU but partilly (a part of the reads against a part of a sequence).
+def readsCoverageGPU_partial(id_seq,id_chunk):
+    pass
 
 ## Equivalent to the readsCoverage function
 #

tests/test_GPUMapper.py

@@ -10,6 +10,7 @@ import sys
 sys.path.append("..")
 
 from _modules.GPU_chunkMapper import GPU_chunkMapper
+from _modules.GPU_mapper import getGPUMappingResults
 from _modules.IData_handling import refData,readsExtractsS
 
 from data_refseq import refseq_list
@@ -41,20 +42,7 @@ class TestGPUMapper (unittest.TestCase):
         re.nb_extracts=nb_extracts
         re.nb_reads=16895
         mapper=GPU_chunkMapper(re,rd,ficdir,2)
-        # print "nb extracts=",re.nb_extracts
-        # print "re.nb_extracts//2=",re.nb_extracts//2
-        print all_extracts_l[33790]
-        # print all_extracts_l[33791]
-        #print all_extracts_l
-        # print all_extracts_l[11104] # appears once in seq 0
-        # print all_extracts_l[11105] # appears 1 in seq 0
-        # print all_extracts_l[11106] # does not appear
-        # print all_extracts_l[11107] # does not appear
-        #
-        # print all_extracts_l[60000] # does not appear
-        # print all_extracts_l[60001] # does not appear
-        # print all_extracts_l[60002] # appears in seq3
-        #print all_extracts_l[
+        #print all_extracts_l[33790]
         mapper.doMapping()
         self.assertTrue(os.path.isfile(os.path.join(ficdir, "read_occ_positions0_0")))
         self.assertTrue(os.path.isfile(os.path.join(ficdir, "read_occ_positions0_1")))
@@ -78,6 +66,17 @@ class TestGPUMapper (unittest.TestCase):
         self.assertTrue(os.path.isfile(os.path.join(ficdir, "reduced4_0")))
         self.assertTrue(os.path.isfile(os.path.join(ficdir, "reduced4_1")))
 
+        self.assertTrue(os.path.isfile(os.path.join(ficdir, "max_occ0_0")))
+        self.assertTrue(os.path.isfile(os.path.join(ficdir, "max_occ0_1")))
+        self.assertTrue(os.path.isfile(os.path.join(ficdir, "max_occ1_0")))
+        self.assertTrue(os.path.isfile(os.path.join(ficdir, "max_occ1_1")))
+        self.assertTrue(os.path.isfile(os.path.join(ficdir, "max_occ2_0")))
+        self.assertTrue(os.path.isfile(os.path.join(ficdir, "max_occ2_1")))
+        self.assertTrue(os.path.isfile(os.path.join(ficdir, "max_occ3_0")))
+        self.assertTrue(os.path.isfile(os.path.join(ficdir, "max_occ3_1")))
+        self.assertTrue(os.path.isfile(os.path.join(ficdir, "max_occ4_0")))
+        self.assertTrue(os.path.isfile(os.path.join(ficdir, "max_occ4_1")))
+
         i_occ_pos = np.loadtxt(os.path.join(ficdir, "read_occ_positions0_0"), dtype=np.int32, delimiter='|')
         self.assertEqual(i_occ_pos[7677],1864)
         self.assertEqual(i_occ_pos[7678],1954)
@@ -106,6 +105,13 @@ class TestGPUMapper (unittest.TestCase):
         self.assertTrue(found1)
         self.assertTrue(found2)
 
+        f_nb_occ=open(os.path.join(ficdir,"max_occ0_0"))
+        str_nb_occ=f_nb_occ.read()
+        f_nb_occ.close()
+        max_nb_occ=int(str_nb_occ)
+        self.assertEqual(max_nb_occ,15690)
+
+
         reduced = np.loadtxt(os.path.join(ficdir, "reduced2_0"), dtype=np.int64, delimiter='|')
         found = False
         i = 0
@@ -116,8 +122,8 @@ class TestGPUMapper (unittest.TestCase):
                     found=True
                     nb_occ=reduced[i+1]
             i+=1
-        #self.assertTrue(found)
-        #self.assertEqual(nb_occ,2)
+        self.assertTrue(found)
+        self.assertEqual(nb_occ,2)
 
         reduced = np.loadtxt(os.path.join(ficdir, "reduced2_1"), dtype=np.int64, delimiter='|')
         found = False
@@ -134,6 +140,20 @@ class TestGPUMapper (unittest.TestCase):
 
 
     def testIterationOverMappingResults(self):
+        rd = refData(refseq_list, 20, "")
+        re = readsExtractsS(20)
+        re.r_extracts_list = all_extracts_l
+        re.nb_extracts = nb_extracts
+        re.nb_reads = 16895
+        mapper = GPU_chunkMapper(re, rd, ficdir, 2)
+        mapper.doMapping()
+        mres=getGPUMappingResults(os.path.join(ficdir, "reduced3_1"), os.path.join(ficdir, "read_occ_positions3_1"),
+                                  os.path.join(ficdir,"max_occ3_1"), re.nb_reads, re.nb_extracts, rd, single=True)
+        for r, r_hostseq, idx_read in mres.mappingResultsForSeq(refseq_list[3]):
+            if (idx_read==6691):
+                self.assertEqual(r.)
+
+
         pass
 
 if __name__ == "__main__":