Commit 4ce163ae by Eric Coissac

Simplify computation of reserved memory area and avoid computation of LCS when…

Simplify computation of reserved memory area and avoid computation of LCS when LCSmin is > the shortest sequence.
parent 5117bf3a
......@@ -828,7 +828,11 @@ double sse_banded_lcs_align(int16_t* seq1, int16_t* seq2, int l1, int l2,
LCSmin = calculateLCSmin(l1, l2, threshold, normalize, reference, similarity_mode);
if (LCSmin > l2) return NA_REAL;
if (LCSmin > l2) {
*lcs_length = NA_REAL;
*ali_length = NA_REAL;
return NA_REAL;
}
bandLengthLeft = calculateLeftBandLength(l1, LCSmin);
bandLengthRight = calculateRightBandLength(l2, LCSmin);
......@@ -937,7 +941,6 @@ SEXP ROBI_lcs(SEXP sequences1,SEXP sequences2,SEXP length1,SEXP length2,
int c_lcs;
int c_alilength;
int sizeToAllocateForBand, sizeToAllocateForSeqs;
int maxBLL;
int shift;
int16_t* address;
......@@ -1028,18 +1031,10 @@ SEXP ROBI_lcs(SEXP sequences1,SEXP sequences2,SEXP length1,SEXP length2,
if (lmax > SHRT_MAX)
error("The longest sequence is too long (> 32000bp)");
// Calculate the minimum LCS length corresponding to the threshold
LCSmin = calculateLCSmin(lmin, lmin, c_threshold, c_normalize, c_reference, c_similarity_mode);
LCSmin = 0;
// printf("LCSmin=%d\n",LCSmin);
// printf("threshold=%f\n",(float)c_threshold);
SSEregisters = obi_aligned_on_16(R_alloc(6 * (lmax - LCSmin + 4),sizeof(um128)),
SSEregisters = obi_aligned_on_16(R_alloc(7 * lmax,sizeof(um128)),
&shift);
// printf("SSEregisters = %d bytes\n",6 * (int)((lmax - LCSmin + 4)*sizeof(um128)));
// Allocate space for matrix band if the alignment length must be computed
sizeToAllocateForBand = calculateSizeToAllocate(lmax, LCSmin);
......@@ -1048,8 +1043,7 @@ SEXP ROBI_lcs(SEXP sequences1,SEXP sequences2,SEXP length1,SEXP length2,
// Allocate space for the int16_t arrays representing the sequences
maxBLL = calculateLeftBandLength(lmax, LCSmin);
sizeToAllocateForSeqs = 2*maxBLL+lmax;
sizeToAllocateForSeqs = 3*lmax;
iseq1L = (int16_t*) R_alloc(sizeToAllocateForSeqs, sizeof(int16_t));
iseq1l = (int16_t*) R_alloc(sizeToAllocateForSeqs, sizeof(int16_t));
......@@ -1071,24 +1065,24 @@ SEXP ROBI_lcs(SEXP sequences1,SEXP sequences2,SEXP length1,SEXP length2,
if (csource1 != prev_source1) {
iniSeq(iseq1L,sizeToAllocateForSeqs,0);
putSeqInSeq(iseq1L+maxBLL, csource1, l1, TRUE);
putSeqInSeq(iseq1L+lmax, csource1, l1, TRUE);
iniSeq(iseq1l,sizeToAllocateForSeqs,255);
putSeqInSeq(iseq1l+maxBLL, csource1, l1, FALSE);
putSeqInSeq(iseq1l+lmax, csource1, l1, FALSE);
prev_source1 = csource1;
}
if (l1 > l2) {
Lseq = iseq1L + maxBLL;
Lseq = iseq1L + lmax;
iniSeq(iseq2,sizeToAllocateForSeqs,255);
lseq = iseq2 + maxBLL;
lseq = iseq2 + lmax;
putSeqInSeq(lseq, csource2, l2, FALSE);
Ll = l1;
ll = l2;
}
else {
lseq = iseq1l + maxBLL;
lseq = iseq1l + lmax;
iniSeq(iseq2,sizeToAllocateForSeqs,0);
Lseq = iseq2 + maxBLL;
Lseq = iseq2 + lmax;
putSeqInSeq(Lseq, csource2, l2, TRUE);
Ll = l2;
ll = l1;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment