_lcs.ext.1.c 2.83 KB
Newer Older
Celine Mercier committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
#include "_lcs.h"

#include <string.h>
#include <stdlib.h>
#include <limits.h>

#include <stdio.h>



// Allocate a band allowing to align sequences of length : 'length'

column_t* allocateColumn(int length,column_t *column, bool mode8bits)
{
	int size;
	bool newc = false;

			// The band length should be equal to the length
			// of the sequence + 7 for taking into account its
			// shape

	size = (length+1) * ((mode8bits) ? sizeof(int8_t):sizeof(int16_t));


			// If the pointer to the old column is NULL we allocate
			// a new column

	if (column==NULL)
	{

		column = malloc(sizeof(column_t));
		if (!column)
			return NULL;

		column->size = 0;
		column->data.shrt=NULL;
		column->score.shrt=NULL;
		newc = true;
	}

			// Otherwise we check if its size is sufficient
			// or if it should be extended

	if (size > column->size)
	{
		int16_t *old = column->data.shrt;
		int16_t *olds= column->score.shrt;

		column->data.shrt = malloc(size);
		column->score.shrt= malloc(size);

		if (column->data.shrt==NULL || column->score.shrt==NULL)
		{
			fprintf(stderr,"Allocation Error on column for a size of %d\n" , size);
			column->data.shrt = old;
			column->score.shrt= olds;

			if (newc)
			{
				free(column);
				column=NULL;
				return NULL;
			}
			return NULL;
		}
		else
			column->size = size;
	}

	return column;
}

void freeColumn(column_p column)
{
	if (column)
	{
		if (column->data.shrt)
			free(column->data.shrt);

		if (column->score.shrt)
			free(column->score.shrt);

		free(column);
	}
}

int fastLCSScore(const char* seq1, const char* seq2,column_pp column,int32_t* lpath)
{
	return fastLCSScore16(seq1,seq2,column,lpath);
}

int simpleLCS(const char* seq1, const char* seq2,column_pp ppcolumn,int32_t* lpath)
{
	int lseq1,lseq2;		// length of the both sequences
	int lcs;
	int itmp;				// tmp variables for swap
	const char* stmp;		//
	int32_t *score;
	int32_t *path;
	column_t *column;
	int32_t i,j;
	int32_t sl,su,sd;
	int32_t pl,pu,pd;

		// 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;
	}

	lseq1++;
	lseq2++;

							// a band sized to the smallest sequence is allocated

	if (ppcolumn)
		column = *ppcolumn;
	else
		column=NULL;

	column = allocateColumn(lseq1*2,column,0);
	score = (int32_t*) column->score.shrt;
	path = (int32_t*) column->data.shrt;

	memset(score,0,lseq1 * sizeof(int32_t));

	for (j=0; j < lseq1; j++)
		path[j]=j;

	for (i=1; i< lseq2; i++)
	{
		sl=0;
		pl=i;
		for (j=1; j < lseq1; j++)
		{
			sd=score[j-1] + (seq2[i-1]==seq1[j-1] ? 1:0);
			pd=path[j-1]  + 1;

			su=score[j];
			pu=path[j] + 1;

			score[j-1]=sl;

			if (su > sl) sl=su, pl=pu;
			if (sd > sl) sl=sd, pl=pd;
		}
	}

	lcs = sl;
	if(lpath) *lpath=pl;

	if (ppcolumn)
		*ppcolumn=column;
	else
		freeColumn(column);

	return lcs;
}