obitools/obitools
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Remove dependency on SSE2

Browse Source
This commit is contained in:
Eric Coissac
2010-09-04 05:40:41 +00:00
parent ea44edcaac
commit 9b8752aa15
6 changed files with 403 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/obitools/align/__init__.py
View File

@ -2,7 +2,7 @@
from _nws import NWS from _nws import NWS
from _upperbond import indexSequences from _upperbond import indexSequences
from _lcs import LCS from _lcs import LCS,lenlcs
from _assemble import DirectAssemble, ReverseAssemble from _assemble import DirectAssemble, ReverseAssemble
from _qsassemble import QSolexaDirectAssemble,QSolexaReverseAssemble from _qsassemble import QSolexaDirectAssemble,QSolexaReverseAssemble
from _rassemble import RightDirectAssemble as RightReverseAssemble from _rassemble import RightDirectAssemble as RightReverseAssemble

2
src/obitools/align/_lcs.c
View File

@ -1,4 +1,4 @@
/* Generated by Cython 0.13 on Tue Aug 31 15:06:12 2010 */ /* Generated by Cython 0.13 on Wed Sep 1 07:59:42 2010 */
#define PY_SSIZE_T_CLEAN #define PY_SSIZE_T_CLEAN
#include "Python.h" #include "Python.h"

26
src/obitools/align/_lcs.ext.2.c Normal file
View File

@ -0,0 +1,26 @@
#include "_lcs.h"
#include <string.h>
#include <stdlib.h>
#include <limits.h>
#include <stdio.h>
#define VSIZE (8)
#define VTYPE vInt16
#define STYPE int16_t
#define CMENB shrt
#define VMODE false
#define FASTLCSSCORE fastLCSScore16
#define INSERT_REG _MM_INSERT_EPI16
#define EXTRACT_REG _MM_EXTRACT_EPI16
#define EQUAL_REG _MM_CMPEQ_EPI16
#define ADD_REG _MM_ADD_EPI16
#define SET_CONST _MM_SET1_EPI16
#define GET_MAX _MM_MAX_EPI16
#define MIN_SCORE INT16_MIN
#include "_lcs_fast.h"

26
src/obitools/align/_lcs.ext.3.c Normal file
View File

@ -0,0 +1,26 @@
#include "_lcs.h"
#include <string.h>
#include <stdlib.h>
#include <limits.h>
#include <stdio.h>
#define VSIZE (16)
#define VTYPE vInt8
#define STYPE int8_t
#define CMENB byte
#define VMODE true
#define FASTLCSSCORE fastLCSScore8
#define INSERT_REG _MM_INSERT_EPI8
#define EXTRACT_REG _MM_EXTRACT_EPI8
#define EQUAL_REG _MM_CMPEQ_EPI8
#define ADD_REG _MM_ADD_EPI8
#define SET_CONST _MM_SET1_EPI8
#define GET_MAX _MM_MAX_EPI8
#define MIN_SCORE INT8_MIN
#include "_lcs_fast.h"

348
src/obitools/align/_lcs_fast.h Normal file
View File

@ -0,0 +1,348 @@
static void printreg(VTYPE r)
{
STYPE a0,a1,a2,a3,a4,a5,a6,a7;
#if VMODE
STYPE a8,a9,a10,a11,a12,a13,a14,a15;
#endif
a0= EXTRACT_REG(r,0);
a1= EXTRACT_REG(r,1);
a2= EXTRACT_REG(r,2);
a3= EXTRACT_REG(r,3);
a4= EXTRACT_REG(r,4);
a5= EXTRACT_REG(r,5);
a6= EXTRACT_REG(r,6);
a7= EXTRACT_REG(r,7);
#if VMODE
a8= EXTRACT_REG(r,8);
a9= EXTRACT_REG(r,9);
a10= EXTRACT_REG(r,10);
a11= EXTRACT_REG(r,11);
a12= EXTRACT_REG(r,12);
a13= EXTRACT_REG(r,13);
a14= EXTRACT_REG(r,14);
a15= EXTRACT_REG(r,15);
#endif
printf( "a00 : %4d a01 : %4d a02 : %4d a03 : %4d "
"a04 : %4d a05 : %4d a06 : %4d a07 : %4d "
#if VMODE
"a08 : %4d a09 : %4d a10 : %4d a11 : %4d "
"a12 : %4d a13 : %4d a14 : %4d a15 : %4d "
#endif
"\n"
, a0,a1,a2,a3,a4,a5,a6,a7
#if VMODE
, a8,a9,a10,a11,a12,a13,a14,a15
#endif
);
}
static inline VTYPE insert_reg(VTYPE r, STYPE v, int p)
{
switch (p) {
case 0: return INSERT_REG(r,v,0);
case 1: return INSERT_REG(r,v,1);
case 2: return INSERT_REG(r,v,2);
case 3: return INSERT_REG(r,v,3);
case 4: return INSERT_REG(r,v,4);
case 5: return INSERT_REG(r,v,5);
case 6: return INSERT_REG(r,v,6);
case 7: return INSERT_REG(r,v,7);
#if VMODE
case 8: return INSERT_REG(r,v,8);
case 9: return INSERT_REG(r,v,9);
case 10: return INSERT_REG(r,v,10);
case 11: return INSERT_REG(r,v,11);
case 12: return INSERT_REG(r,v,12);
case 13: return INSERT_REG(r,v,13);
case 14: return INSERT_REG(r,v,14);
case 15: return INSERT_REG(r,v,15);
#endif
}
return _MM_SETZERO_SI128();
}
static inline STYPE extract_reg(VTYPE r, int p)
{
switch (p) {
case 0: return EXTRACT_REG(r,0);
case 1: return EXTRACT_REG(r,1);
case 2: return EXTRACT_REG(r,2);
case 3: return EXTRACT_REG(r,3);
case 4: return EXTRACT_REG(r,4);
case 5: return EXTRACT_REG(r,5);
case 6: return EXTRACT_REG(r,6);
case 7: return EXTRACT_REG(r,7);
#if VMODE
case 8: return EXTRACT_REG(r,8);
case 9: return EXTRACT_REG(r,9);
case 10: return EXTRACT_REG(r,10);
case 11: return EXTRACT_REG(r,11);
case 12: return EXTRACT_REG(r,12);
case 13: return EXTRACT_REG(r,13);
case 14: return EXTRACT_REG(r,14);
case 15: return EXTRACT_REG(r,15);
#endif
}
return 0;
}
#define GET_H_SYMBOLE(s,p) ((p && p < lseq1) ? (s)[(p)-1]:255)
#define GET_V_SYMBOLE(s,p) ((p && p < lseq2) ? (s)[(p)-1]:0)
#define LSHIFT_SCORE(r) { r = _MM_SLLI_SI128((r),sizeof(STYPE)); }
#define SET_H_SYMBOLE(r,p,s) { r = insert_reg((r),(STYPE)GET_H_SYMBOLE(seq1,(s)),(p)); }
#define PUSH_V_SYMBOLE(r,s) { r = insert_reg(_MM_SLLI_SI128((r),sizeof(STYPE)),(STYPE)GET_V_SYMBOLE(seq2,(s)),0); }
#define EQUAL(f1,f2) _MM_AND_SI128(EQUAL_REG((f1),(f2)),SET_CONST(1))
int FASTLCSSCORE(const char* seq1, const char* seq2,column_pp ppcolumn)
{
int lseq1,lseq2; // length of the both sequences
int itmp; // tmp variables for swap
const char* stmp; //
int nbands; // Number of bands of width eight in the score matrix
int lastband; // width of the last band
VTYPE minus1;
VTYPE minus2;
VTYPE current;
VTYPE fhseq;
VTYPE fvseq;
VTYPE match;
int band;
int line;
int limit;
int lcs;
int h;
int i;
column_t *column;
// Made seq1 the longest sequences
lseq1=strlen(seq1);
lseq2=strlen(seq2);
if (lseq1 < lseq2)
{
itmp=lseq1;
lseq1=lseq2;
lseq2=itmp;
stmp=seq1;
seq1=seq2;
seq2=stmp;
}
// we add one to the both length for taking into
// account the extra line and column in the score
// matrix
lseq1++;
lseq2++;
// a band sized to the smallest sequence is allocated
if (ppcolumn)
column = *ppcolumn;
else
column=NULL;
column = allocateColumn(lseq2,column,VMODE);
// Check memory allocation
if (column == NULL)
return -1;
for (i=0; i<lseq2;i++)
column->data.CMENB[i]=MIN_SCORE;
nbands = lseq1 / VSIZE; // You have VSIZE element in one SSE register
// Alignment will be realized in nbands
lastband = lseq1 - (nbands * VSIZE); // plus one of width lastband except if
// lastband==0
if (lastband) nbands++;
else lastband=VSIZE;
lastband--;
// printf("seq1 : %s seq2 : %s\n",seq1,seq2);
minus2 = SET_CONST(MIN_SCORE);
minus1 = _MM_SETZERO_SI128();
h=0;
fhseq = _MM_SETZERO_SI128();
fvseq = _MM_SETZERO_SI128();
//
// Beginnig of the first band
//
for (line = 0; line < VSIZE-1; line++,h++)
{
// printf("line= %4d h= %4d\n",line,h);
SET_H_SYMBOLE(fhseq,line,h)
PUSH_V_SYMBOLE(fvseq,line)
minus2 = insert_reg(minus2,0,line);
minus1 = insert_reg(minus1,MIN_SCORE,line); // 0 avant
match = EQUAL(fhseq,fvseq);
// printreg(fvseq);
// printreg(fhseq);
// printreg(match);
// printf("================================\n");
current = minus1;
// printf("Vert = "); printreg(current);
LSHIFT_SCORE(minus1)
minus1=insert_reg(minus1,(column)->data.CMENB[line],0);
// printf("Horz = "); printreg(minus1);
current = GET_MAX(current,minus1);
// printf("BstHV= "); printreg(current);
//
// printf("Diag = "); printreg(ADD_REG(minus2,match));
current = GET_MAX(current,ADD_REG(minus2,match));
printf("line %d :Best = ",line); printreg(current);
//
// printf("================================\n");
// printf("================================\n");
minus2=minus1;
minus1=current;
// printf("min2 = "); printreg(minus2);
// printf("min1 = "); printreg(minus1);
// printf("================================\n");
}
(column)->data.CMENB[lseq2-VSIZE+line]=EXTRACT_REG(current,VSIZE-1);
for (band=0; band < nbands; band++)
{
printf("band : %d\n",band);
SET_H_SYMBOLE(fhseq,line,h)
minus2 = insert_reg(minus2,0,line);
minus1 = insert_reg(minus1,MIN_SCORE,line); // 0 avant
h++;
for (; line < lseq2; line++)
{
printf("Je tourne avec line= %d \n",line);
PUSH_V_SYMBOLE(fvseq,line)
match = EQUAL(fhseq,fvseq);
// printreg(fvseq);
// printreg(fhseq);
// printreg(match);
// printf("================================\n");
current = minus1;
// printf("Vert = "); printreg(current);
LSHIFT_SCORE(minus1)
// printf("line = %d --> get = %d\n",line,(column)->data.CMENB[line]);
minus1=insert_reg(minus1,(column)->data.CMENB[line],0);
// printf("Horz = "); printreg(minus1);
current = GET_MAX(current,minus1);
// printf("BstHV= "); printreg(current);
//
// printf("Diag = "); printreg(ADD_REG(minus2,match));
current = GET_MAX(current,ADD_REG(minus2,match));
// printf("Best = "); printreg(current);
//
// printf("================================\n");
// printf("================================\n");
printf("line %d :Best = ",line); printreg(current);
minus2=minus1;
minus1=current;
// printf("min2 = "); printreg(minus2);
// printf("min1 = "); printreg(minus1);
// printf("================================\n");
// Store the last current score in extra column
// if (line >= VSIZE)
(column)->data.CMENB[line-VSIZE+1]=EXTRACT_REG(current,VSIZE-1);
// printf("Je stocke line= %d la valeur %d\n",line-VSIZE,EXTRACT_REG(current,7));
}
// end of the band and beginnig of the next one
limit=(band==(nbands-1)) ? lastband:VSIZE-1;
printf("band : %d nbands : %d lastband %d limit %d\n", band,nbands,lastband,limit);
for (line = 0; line < limit; line++,h++)
{
// printf("Je fini avec line= %d \n",line);
SET_H_SYMBOLE(fhseq,line,h)
PUSH_V_SYMBOLE(fvseq,line)
minus2 = insert_reg(minus2,0,line);
minus1 = insert_reg(minus1,MIN_SCORE,line);
match = EQUAL(fhseq,fvseq);
printf("\n");
printf("fhseq = "); printreg(fhseq);
printf("fvseq = "); printreg(fvseq);
printf("----------------------------------------------------------------\n");
printf("match = "); printreg(match);
current = minus1;
LSHIFT_SCORE(minus1)
minus1=insert_reg(minus1,(column)->data.CMENB[line],0);
current = GET_MAX(current,minus1);
current = GET_MAX(current,ADD_REG(minus2,match));
printf("currt = "); printreg(current);
minus2=minus1;
minus1=current;
(column)->data.CMENB[lseq2-VSIZE+line]=EXTRACT_REG(current,VSIZE-1);
// printf("Je stocke line= %d la valeur %d\n",lseq2-VSIZE+line,(column)->data.CMENB[lseq2-VSIZE+line]);
}
}
printf("\n");
printf("line = %d, h= %d, lastband = %d\n",line,h,lastband);
printf("currt = "); printreg(current);
lcs = extract_reg(current,lastband);
// printf("lastband = %d (%d) lcs = %d\n",lastband,lseq2,lcs);
if (ppcolumn)
*ppcolumn=column;
else
freeColumn(column);
return lcs;
}

2
src/obitools/align/_upperbond.ext.1.c
View File

@ -63,7 +63,7 @@ int buildTable(const char* sequence, unsigned char *table, int *count)
memset(table,0,256*sizeof(unsigned char)); memset(table,0,256*sizeof(unsigned char));
// encode ascii sequence with A : 00 C : 01 T: 10 G : 11 // encode ascii sequence with A : 00 C : 01 T: 10 G : 11
_MM_LOADU_SI128((vUInt8*)s);
for(frag.m=_MM_LOADU_SI128((vUInt8*)s); for(frag.m=_MM_LOADU_SI128((vUInt8*)s);
! anyzerom128(frag.m); ! anyzerom128(frag.m);
s+=12,frag.m=_MM_LOADU_SI128((vUInt8*)s)) s+=12,frag.m=_MM_LOADU_SI128((vUInt8*)s))