Partage

[FAIT][Défis] #9 : La recherche en génétique

23 décembre 2011 à 20:29:03

La recherche en génétique



Bonjour bonjour ! :)

Aujourd'hui, nous allons recoder la fonction strstr.

Mise en situation



La recherche en génétique progresse chaque jour. Elle a déjà donné d'excellents résultats, et continue d'être active. Seulement, la quantité d'informations contenue dans un brin d'ADN est simplement énorme... Il est donc logique de la traiter par ordinateur !

Les séquences ADN peuvent être représentées par un alphabet de quatre lettres : A, T, G ,C. Pour connaître la signification de ces lettres, je vous invite à consulter la page wikipédia concernant l'ADN. Une opération courante dans le domaine est la recherche d'un gène (qui sera ici une séquence de lettres) dans une séquence ADN bien plus grande.

Nous allons donc coder un programme capable de retrouver une sous-chaîne dans un texte !

Énoncé



Écrivez un programme qui lit un gène sur la première ligne, une séquence ADN sur la deuxième, puis qui donne l'indice ou démarre le gène dans la séquence.

Exemple :
ATGC
ATTGGCCTATGCGGA


En sortie :
9


Niveau 1

Écrivez le programme de recherche en utilisant l'algorithme le plus simple. Bien entendu, il est interdit de recourir à la fonction strstr, le but est que vous la fassiez vous-mêmes. :)

Niveau 2

Améliorez votre programme de recherche pour qu'il utilise l'algorithme de Rabin-Karp. L'idée est la suivante : nous allons représenter chaque lettre par un chiffre. Par exemple, A=0, T=1, G=2 et C=3. Nous pouvons alors écrire un gène, tel que ATGC, comme un nombre entier, qui vaut 123.

Pour chercher le gène, il suffit de prendre les 4 premiers caractères de la séquence, d'opérer la même conversion, puis de comparer les résultats. Il nous faut ensuite convertir les caractères 2 à 5 de la séquence, et comparer à nouveau... Mais, lorsque nous avons converti les caractères 1 à 4 tout à l'heure, nous avons en particulier converti les caractères 2 à 4 en entiers ! N'y a-t-il pas moyen d'en profiter ?

Reprenons l'exemple précédent pour nous fixer les idées. Au premier tour, de boucle, la situation était :
0123
0112GCCTATGCGGA


Comme 123 et 112 sont différents, le gène n'est pas en position 1. Nous nous décalons donc d'une case dans la séquence :
_0123
A1122CCTATGCGGA


Donc, sachant que la séquence actuelle vaut 0112, et que le prochain caractère sera associé au nombre « 2 », comment en déduire que le prochain entier à comparer sera 1122 ?

Niveau 3

Le code précédent constitue une bonne amélioration de la recherche naïve codée au niveau 1. Cependant, nous sommes limités par la taille des entiers que nous pouvons manipuler ! En effet, comme vous le savez sûrement, une variable de type int ou long comporte une limite de grandeur pour les entiers stockés. Cette limite maximale varie selon la machine ; vous pouvez l'afficher grâce au code suivant.
#include <limits.h>

int main(void)
{
    printf("%ul\n", ULONG_MAX);
    return 0;
}


Ainsi, nous ne pouvons pas rechercher des gènes arbitrairement longs ! Il faut que leur « version numérique » tienne dans un long (ou autre type que vous aurez choisi). Nous avons deux solutions pour pallier à ce problème.

La première est d'utiliser un type de données capables de manipuler des nombres arbitrairement grands. L'objet de l'un des exercices du forum est d'écrire un tel type, vous pouvez donc le réutiliser si ça vous amuse. :)
Cependant, les opérations sur de grands nombres sont relativement lentes, et nous perdrions le bonus en rapidité que nous fait gagner l'algorithme de Rabin-Karp.

La seconde solution est d'utiliser un hash. L'idée est la suivante : nous allons appliquer un modulo au gène et au bloc à tester. Notez que ceci ne modifie en rien l'astuce du niveau 2 : elle fonctionne toujours. :)

En revanche, lorsque nous trouverons deux nombres égaux, ne ne seront plus complètement sûrs que les gènes correspondants le sont bien aussi (du fait du modulo). Il devient donc nécessaire de faire une comparaison textuelle en cas d'égalité. Par contre, si les deux nombres trouvés sont différents, les deux séquences sont nécessairement différentes, et il n'y a plus besoin de test (ce qui est un peu l'intérêt quand même).

Conclusion



À l'issue de cet exercice, nous avons implémenté un algorithme efficace de recherche de sous-chaîne, qui évite « souvent » les comparaisons inutiles.

Il existe encore de nombreuses autres méthodes pour pour rechercher efficacement une sous-chaîne dans une chaîne. Vous pensiez que l'exercice de M@teo21 était fini, enterré, oublié ? Vous avez pourtant de quoi vous amuser encore longtemps... :diable:

Bon courage à tous ! :)
GuilOooo


______________________________________________________________________________________________

Si vous souhaitez plus de renseignements sur les défis, rendez vous sur le topic de recensement des défis, vous y trouverez également les règles principales.

Cet exercice a été écrit par GuilOooo, merci à lui !

Participants



Participants Code
redspikers Niveaux 1, 2 et 3
GurneyH Niveau 3
Lucas-84 Niveau 1
@che Niveau 1
'Housenka Niveau 1
paraze Niveau 1
yoch Niveaux 2 et 3
Magykdu13 Niveau 1
Taurre Niveau 2
23 décembre 2011 à 21:08:38

Niveau 1
#include <stdio.h>
#include <stdlib.h>
#include <string.h>


int recherche(char* motif, char* chaine){
    int i, k=0;
    int tailleMotif = strlen(motif);
    int tailleChaine = strlen(chaine);
    for(i=0; i<(tailleChaine-tailleMotif+1); i++){
        for(k=0; (k<tailleMotif)&&(chaine[i+k] == motif[k]); k++);
        if (k==tailleMotif)
            return (i+k)-tailleMotif+1;

    }
    return -1;

}


int main(int argc, char** argv)
{
    char* motif = malloc(sizeof(char)*50);
    char* chaine = malloc(sizeof(char)*50000);

    printf("Entrez le motif : ");
    scanf("%s", motif);

    printf("Entrez la chaine : ");
    scanf("%s", chaine);

    printf("Resultat : %d", recherche(motif, chaine));

    free(motif);
    free(chaine);

    return 0;
}


Niveau 2&3
Le code est assez moche, mais je me demande si on peut améliorer le calcul à la ligne 91.
Et je suis pas sur d'avoir mis le modulo au bon endroit :-°

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#define TRUE  1
#define FALSE 0
#define MODULO INT_MAX
int recherche(char* motif, char* chaine){
    int i, k=0;
    int tailleMotif = strlen(motif);
    int tailleChaine = strlen(chaine);
    for(i=0; i<(tailleChaine-tailleMotif+1); i++){
        for(k=0; (k<tailleMotif)&&(chaine[i+k] == motif[k]); k++);
        if (k==tailleMotif)
            return (i+k)-tailleMotif+1;

    }
    return -1;

}

int convert(char carac){
    //Transforme un caractere de nucléotide dans son entier correspondant
    int lettre;
    switch (carac){
            case 'A': lettre = 0; break;
            case 'T': lettre = 1; break;
            case 'G': lettre = 2; break;
            case 'C': lettre = 3; break;
            default: lettre = 0; break;
        }
    return lettre;
}

int hash(char* str){
    //Transforme une chaine de nucléotide en son entier correspondant
    int i;
    int tailleStr = strlen(str);
    int lettre, resultat = 0;
    for(i=0; i<tailleStr; i++){
        lettre=convert(str[i]);
        resultat += lettre*puiss10(tailleStr-i-1);
    }
    return resultat%MODULO;
}

int strCmp(char* str1, char* str2){
    //Histoire de pas trop utiliser string.h non plus ...
    int i = 0;
    if (strlen(str1) == strlen(str2)){
        for(i=0; i<strlen(str1); i++){
            if(str1[i] != str2[i])
                return FALSE;
        }
        return TRUE;
    }
    return FALSE;

}

int puiss10(int puissance){
    //Retourne 10^puissance
    int resultat = 1;
    for(; puissance>0;puissance--)
        resultat *= 10;
    return resultat;
}
int rechercheRabinKarp(char* motif, char* chaine){
    int i;
    int tailleMotif = strlen(motif);
    int tailleChaine = strlen(chaine);
    int hashMotif = hash(motif);

    //On extrait une sous-chaine du début de la séquence,
    //de la taille du motif
    char* sousChaine = malloc(tailleMotif+1);
    strncpy(sousChaine, chaine, tailleMotif);
    sousChaine[tailleMotif] = '\0';
    //Fin de l'extraction

    int hashSousChaine = hash(sousChaine);

    for(i=0; i<(tailleChaine-tailleMotif+1); i++){
        if(hashMotif == hashSousChaine){
            strncpy(sousChaine, chaine+i, tailleMotif);
            sousChaine[tailleMotif] = '\0';
            if(strCmp(motif, sousChaine))
                return i+1;
        }
        hashSousChaine = ((hashSousChaine%(puiss10(tailleMotif-1))*10
                          +
                          convert(chaine[i+tailleMotif])))%MODULO;

    }
    free(sousChaine);
    return -1;

}

int main(int argc, char** argv)
{
    char* motif = malloc(50000);
    char* chaine = malloc(500000);

    printf("Entrez le motif : ");
    scanf("%s", motif);

    printf("Entrez la chaine : ");
    scanf("%s", chaine);

    printf("Resultat : %d", rechercheRabinKarp(motif, chaine));

    free(motif);
    free(chaine);
    return 0;
}
24 décembre 2011 à 6:12:29

bonjour

C'est moi ou il n'y a que peu de différence entre les niveau 2 et 3(à un modulo et un strncmp près)? :-° Le modulo de plus s'effectue naturellement sur des unsigned.

Pour le niveau 3
#include <stdio.h>
#include <string.h>

struct t_hash
{
    unsigned long value;
    unsigned long mod;
};


static unsigned long powOf10(unsigned long n)
{
    unsigned long ret = 1;
    for(unsigned long i = 0; i < n; i++)
        ret *= 10;

    return ret;
}

static size_t getIndex(char c)
{
    switch(c)
    {
    case 'A':
        return 1;
    case 'C':
        return 2;
    case 'T':
        return 3;
    case 'G':
        return 4;
    default:
        return 0;
    }
}



struct t_hash hash_init(char const *s, size_t n)
{
    struct t_hash ret = {0, 0};
    ret.mod = powOf10(n - 1);
    for(size_t i = 0; s[i] && i < n; i++)
        ret.value = (ret.value * 10) + getIndex(s[i]);

    return ret;
}


void hash_update(struct t_hash *hash, char c)
{
    hash->value = ((hash->value % hash->mod) * 10) + getIndex(c);
}





char *find(char const *src, char const *needle)
{
    size_t len = strlen(needle);
    struct t_hash hash_needle = hash_init(needle, len);
    struct t_hash hash_crt = hash_init(src, len);

    for(char const *nxt_char = src + len - 1; *nxt_char; src++)
    {
        if(hash_needle.value == hash_crt.value && strncmp(src, needle, len) == 0)
            return (char *)src;
        hash_update(&hash_crt, *++nxt_char);
    }

    return NULL;
}


int main(void)
{
    char const *src[5] = {"", "ATGC", "AATGC", "AAATGC", "ATG"};
    char const *needle = "ATGC";
    for(size_t i = 0; i < 5; i++)
    {
         char *p = find(src[i], needle);
         if(p)
            printf("\"%s\" dans \"%s\" en position %zu\n", needle, src[i], p - src[i]);
        else
            printf("\"%s\" pas dans \"%s\"!\n", needle, src[i]);
    }

    return 0;
}


Ce serait pas mal de fournir un fichier de test, histoire qu'on fasse des tests communs.

PS: très sympa ces défis! :)
Zeste de Savoir, le site qui en a dans le citron !
24 décembre 2011 à 11:17:35

Pour le niveau 1.

#include <stdio.h>
#include <stdlib.h>

size_t search_gene(const char * , const char * );
char * get_text(void);

int main(void) {
  char * gene = get_text();
  char * sequence = get_text();  
  
  if (gene == NULL || sequence == NULL) {
     return EXIT_FAILURE;
  } else {
    size_t pos = search_gene(gene, sequence);
  
    if (pos == 0) {
      printf("/\n");
    } else {
      printf("%d\n", pos);
    }
  }
 
  free(gene);
  free(sequence);
  return EXIT_SUCCESS;
}

size_t search_gene(const char * gene, const char * sequence) {
  size_t i, j = 0, k = 0;

  for (i = 0; sequence[i] != '\0'; i++) { 
    if (sequence[i] == gene[j]) {
      k++;
      if (gene[++j] == '\0') {
        return i - k + 2;
      }
    } else {
      j = k = 0;
    }
  }
  return 0;
}

/* <IO> */
char * get_text(void) {
  size_t n = 8;
  char * s = malloc(n * sizeof *s);

  if(s != NULL) {
    size_t i = 0;
    unsigned int e = 1;

    do {
      if (i + 2 > n) {
        s = realloc(s, n + e);
        if (s == NULL) {
          perror("realloc");
          return NULL;
        }
        n += e;
        e *= 2;
      }
      s[i++] = getchar();
    } while(s[i - 1] != '\n');
    s[i - 1] = '\0'; 
    return s;
  }
  perror("malloc");
  return NULL;
}
/* </IO> */
Staff désormais retraité.
24 décembre 2011 à 14:45:22

Tiens intéressant, il y a du monde sur ce défis, allez un petit tour des codes présenté :D

Bon un petit tour des codes présenté ;)

@lucas-84: Humm :euh: :
ATT
TAATTF
/



@redspikers: Pas de problème niveau implémentation, l'exercice niveau 1 est correcte :D Je regarderais le niveau 3 plus tard ^^"
Après, je me permet quand même de faire 2-3 commentaires rapide :
  • L'allocation dynamique n'est pas franchement nécessaire ici ... Surtout si tu ne vérifies pas les retours de malloc ;)
  • sizeof(char) est toujours égal à 1 donc sizeof(char) * i est toujours égal à i
  • Tu pourrais peut-être éviter de chercher la taille de la chaine avec strlen, ce sont des parcours inutiles ;)
  • for(k=0; (k<tailleMotif)&&(chaine[i+k] == motif[k]); k++);// <3


@GurneyH: Rien à dire ^^"
Je vais voir ce que je peux faire pour le fichier de teste ;)
🍊 - Étudiant - Codeur en C | Zeste de Savoir apprenez avec une communauté | Articles  - ♡ Copying is an act of love.
24 décembre 2011 à 15:18:05

Citation : @che


@redspikers: Pas de problème niveau implémentation, l'exercice niveau 1 est correcte :D Je regarderais le niveau 3 plus tard ^^"
Après, je me permet quand même de faire 2-3 commentaires rapide :

  • L'allocation dynamique n'est pas franchement nécessaire ici ... Surtout si tu ne vérifies pas les retours de malloc ;)

Pour l'allocation dynamique, ça me permet d'avoir des tableaux plus grand. Exemples :
Chez moi :
</ul>

char chaine2[5000000]; //plante alors que : 
char* chaine2 = malloc(5000000); //ne plante pas.


Et pourtant je demande la même taille.
sizeof(char) est toujours égal à 1 donc sizeof(char) * i est toujours égal à i</li> Yep, je l'ai fais par réflexe, mais je l'ai pas fait pour le niveau 2&3 ;) Tu pourrais peut-être éviter de chercher la taille de la chaine avec strlen, ce sont des parcours inutiles ;)</li> La je vois pas pourquoi c'est inutile ? D'autant plus que strlen je le fais que au début de la recherche. for(k=0; (k<tailleMotif)&&(chaine[i+k] == motif[k]); k++);// <3 </li> Je pensais que j'allais me faire insulter pour avoir trop voulu compacter :-°


24 décembre 2011 à 15:27:50

Citation

@lucas-84: Humm :euh: :



Tiens donc, bizarre. :-°
Staff désormais retraité.
24 décembre 2011 à 16:08:03

Citation : redspikers

char chaine2[5000000]; //plante alors que : 
char* chaine2 = malloc(5000000); //ne plante pas.
Normal ^^ On peut accéder à plus de mémoire sur le tas que par la pile.
Mais après, on pas besoin de 50 000 cases ^^" 200 c'est largement suffisamment ...

Citation : redspikers

Je pensais que j'allais me faire insulter pour avoir trop voulu compacter :-°

Pas chez moi :-°

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _


Bon on va dire que je me suis amuser pour le niveau 1 ^^"

Pas optimiser du tout :
#include <stdio.h>
#include <string.h>

char* my_strstr(const char* s,const char *t)
{
	return (*t)?(s=strchr(s,*t))?strncmp(s,t, strlen(t) )?(*s)?strstr(s+1,t):NULL:(char*)s:NULL:(char*)s;
}
int main(int argc, char *argv[])
{
  if( argc != 3 )
    printf("Usage: %s STRING SUBSTRING\n", *argv);
  else
    printf("Position: %d\n", (argv[3] = my_strstr( argv[1], argv[2]))? argv[3] - argv[1]:-1);
  return 0;
}



Bonne chance :-°

Optimiser pour un minimum de tour :
#include <stdio.h>
#include <string.h>

int sstrcmp( const char* str, const char* sstr)
{
  for(;*sstr && *str && *str == *sstr; str++, sstr++);
  return !*sstr*2+!*str*3*!sstr-1;
}
int findSubstring( const char* str, const char* sstr)
{
  int i, pos = 0;
  while( ( (i = sstrcmp(str, sstr)) == -1) && *str++)++pos;
  return i>0?pos:-1;
}
int main(int argc, char *argv[])
{
  if( argc != 3 )
    printf("Usage: %s STRING SUBSTRING\n", *argv);
  else
    printf("Position: %d\n", findSubstring( argv[1], argv[2]) );
  return 0;
}
🍊 - Étudiant - Codeur en C | Zeste de Savoir apprenez avec une communauté | Articles  - ♡ Copying is an act of love.
24 décembre 2011 à 16:59:26

Salut à tous,

Je vous donne mon code (niveau 1) que je viens de finir, si vous avez des remarques à faire dessus n'hésites pas :p

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "lecture.h"

int main()
{
    char gene[20] = {0};
    char sequence[100] = {0};
    int vrai = 0;
    int i = 0, j = 0;

    printf("Rentrez votre gene a rechercher : ");
    lire(gene, 20);
    printf("Rentrez la sequence de recherche : ");
    lire(sequence, 1000);

    while (vrai != strlen(gene))
    {
        vrai = 0;

        for (i = 0; gene[i] != '\0'; i++)
        {
            if (gene[i] == sequence[j+i])
            {
                vrai+= 1;
            }
        }

        j++;
    }

    printf("\nIl y a une occurence a la position %d", j);

    return EXIT_SUCCESS;
}


(la fonction lire est celle étudiée dans le tutoriel de M@teo21, qui supprime le \n et vide le buffer).

En ce qui concerne le niveau 2 et 3 j'ai un peu de mal à comprendre le principe en fait : j'ai regardé un peu sur internet en quoi consistait une fonction de hachage mais j'ai pas vraiment compris :p Est-ce que dans ce cas-là on rentre "ATCG" par exemple et que la fonction de hachage va convertir ce ATCG en "0123" ou on rentre 0123 directement (je suppose que c'est la première solution) ?

Et sinon, là c'est juste à titre d'informations, à quoi cela sert de convertir cela en nombres pour comparer leurs valeurs numériques au lieu de directement comparer les caractères entre eux ? Si quelqu'un pourrait me répondre ce serait vraiment sympa :D

Merci d'avance et bonne journée.
24 décembre 2011 à 17:14:07

Citation : 'Housenka


Et sinon, là c'est juste à titre d'informations, à quoi cela sert de convertir cela en nombres pour comparer leurs valeurs numériques au lieu de directement comparer les caractères entre eux ? Si quelqu'un pourrait me répondre ce serait vraiment sympa :D


Plus rapide à comparer. En plus, ça réduis le nombre de comparaisons. Si tu veux plus de détails, demande ;)

Ton code ne donne pas une boucle infinie si le gène n'est pas trouvé ?


Aussi :
lire(sequence, 1000);
Avec :
char sequence[100] = {0};

:(

Dernier point :
while (vrai != strlen(gene))
Niveau performance, c'est horrible ! En effet, on recalcule la taille de gène à chaque tour de boucle ...
🍊 - Étudiant - Codeur en C | Zeste de Savoir apprenez avec une communauté | Articles  - ♡ Copying is an act of love.
24 décembre 2011 à 17:30:50

Ah ok, merci pour la précision :D

Pour la boucle infinie oui, avec cette solution je n'en trouvais pas d'autre pour arrêter le programme s'il n'en trouvait pas parce que la séquence n'a pas réellement de taille maximale prédéfinie :/
Pour le lire(sequence, 1000) c'est juste un faute de frappe, j'ai mis un 0 en trop :p
Pour le while d'accord, je vais créer une variable qui stockera strlen(gene) dans ce cas.

Merci pour tes corrections :D
24 décembre 2011 à 17:39:47

Mais derien, c'est le but du forum :-°

Pour éviter la boucle infinie, tu peux simplement faire :
while (sequence[j] != 0 && vrai != sizeGene) // oû sizeGene est égale à strlen(gene)

Attention ta boucle du for doit être :
for (i = 0; gene[i] != '\0' && sequence[j+i] != '\0'; i++)

Afin d'éviter tout segfaut ;)

Courage, ça m'a l'aire pas mal :D
🍊 - Étudiant - Codeur en C | Zeste de Savoir apprenez avec une communauté | Articles  - ♡ Copying is an act of love.
24 décembre 2011 à 19:00:54

Mon niveau 1 :
#include <stdio.h>
#include <string.h>

int niveau_1(char *s, char *sous_s)
{
    int taille_s = (int)strlen(s);
    int taille_sous_s = (int)strlen(sous_s);
    int indice = -1;
    int i, j;

    for(i = 0; i < taille_s; i++)
    {
        if(s[i] == sous_s[0])
        {
            for(j = 1; j < taille_sous_s; j++)
            {
                if(s[i + j] == sous_s[j])
                {
                    indice = i + 1;
                }

                else
                {
                    indice = -1;
                    break;
                }
            }
        }

        if(indice > -1)
            break;
    }

    return indice;
}

int main(void)
{
    printf("%d", niveau_1("ATTGGCCTATGCGGA", "ATGC"));
    return 0;
}


Sinon, on est obligé de faire l'algo de Rabin-Karp ou on à le droit de présenter un code implémentant un de ces algorithmes par exemple ?
24 décembre 2011 à 19:11:49

Après, tu fais ce que tu veux hein ;)
🍊 - Étudiant - Codeur en C | Zeste de Savoir apprenez avec une communauté | Articles  - ♡ Copying is an act of love.
25 décembre 2011 à 16:02:04

Bonjour (et bonnes fêtes à tous :) ),

Niveau 2 et 3, on peut optimiser un peu en utilisant la base 4 au lieu de la base 10.

Je ne sais pas si j'ai bien compris la technique du modulo, mais si oui, l’intérêt se retrouve limité dans le cas ou la séquence contient des gènes très ressemblants, et que l'on cherche un gène de longueur conséquente (donc demandant un modulo). En effet, le modulo permet seulement de savoir si la fin des 2 chaines (de longueur limitée, ici 31 caractères) sont équivalentes. Donc la fonction de hash n'en est pas véritablement une ici.

Mon code utilise les opérateurs bitwise à la place de l'addition et la multiplication. Le modulo est inutile, car le nombre est automatiquement tronqué par overflow.

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


long long unsigned val(char c)
{
    switch(c)
    {
        case 'A':
            return 0;
        case 'T':
            return 1;
        case 'G':
            return 2;
        case 'C':
            return 3;
    }

    fprintf(stderr, "Invalid sequence : %c found\n", c);
    exit(1);
}

long long unsigned hash (const char* str, size_t len)
{
    long long unsigned h = 0;
    size_t i = 0;

    for (i=0; i < len; ++i)
        h = (h << 2) | val(str[i]);

    return h;
}

long long unsigned secure_shift (long long unsigned n, size_t k)
{
    return (k >= sizeof(long long unsigned) * CHAR_BIT ? 0 : n << k);
}

int karp_rabin(const char* seq, const char* gene)
{
    size_t s_len = strlen(seq);
    size_t g_len = strlen(gene);

    if (s_len < g_len || g_len == 0)
        return -1;

    long long unsigned g_hash = hash(gene, g_len);
    long long unsigned s_hash = hash(seq, g_len);

    if (s_hash == g_hash && !strncmp(seq, gene, g_len))
            return 0;

    size_t i, j;
    for (i=1, j=g_len; j < s_len; ++i, ++j)
    {
        s_hash = ( ( s_hash - secure_shift(val(seq[i-1]), 2*(g_len-1)) ) << 2 ) | val(seq[j]);

        if (g_hash == s_hash && !strncmp(seq+i, gene, g_len))
            return i;
    }

    return -1;
}

int main (void)
{
    char* str = "ATTGGCCTATGCGGA";
    char* pattern = "ATGC";

    printf("%d\n", karp_rabin(str, pattern));

    return 0;
}


EDIT : Il semble que l'exemple donné sur la page française de Wikipédia soit mauvais (induit en erreur), la page anglaise indique d'utiliser un nombre premier pour la base, ce qui semble plus logique.

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


#define BASE    5

long long unsigned val(char c)
{
    switch(c)
    {
        case 'A':
            return 0;
        case 'T':
            return 1;
        case 'G':
            return 2;
        case 'C':
            return 3;
    }

    fprintf(stderr, "Invalid sequence : %c found\n", c);
    exit(1);
}

long long unsigned hash (const char* str, size_t len)
{
    long long unsigned h = 0;
    size_t i = 0;

    for (i=0; i < len; ++i)
        h = (h * BASE) + val(str[i]);

    return h;
}

long long unsigned secure_pow (long long unsigned n, unsigned exp)
{
    long long unsigned ret = 1;

    while (--exp)
    {
        if (ret > ULLONG_MAX / BASE)  // overflow
            return 0;
        ret *= n;
    }

    return ret;
}

int karp_rabin(const char* seq, const char* gene)
{
    size_t s_len = strlen(seq);
    size_t g_len = strlen(gene);

    if (s_len < g_len || g_len == 0)
        return -1;

    long long unsigned g_hash = hash(gene, g_len);
    long long unsigned s_hash = hash(seq, g_len);

    if (s_hash == g_hash && !strncmp(seq, gene, g_len))
            return 0;

    long long unsigned p = secure_pow(BASE, g_len);
    size_t i, j;
    for (i=1, j=g_len; j < s_len; ++i, ++j)
    {
        s_hash = ((s_hash - (val(seq[i-1]) * p)) * BASE) + val(seq[j]);

        if (g_hash == s_hash && !strncmp(seq+i, gene, g_len))
            return i;
    }

    return -1;
}

int main (void)
{
    char* str = "ATTGGCCTATGCGGA";
    char* pattern = "ATGC";

    printf("%d\n", karp_rabin(str, pattern));

    return 0;
}


EDIT 2 : correction du premier code.
25 décembre 2011 à 17:50:01

Toujours épaté pas la concision des codes de yoch. :)
Je ferais les mêmes si j'étais moins maladroit! :lol:
Zeste de Savoir, le site qui en a dans le citron !
25 décembre 2011 à 18:14:26

C'est gentil, GurneyH :) . Ton code à l'air sympa lui aussi (je suis assez minimaliste effectivement, mais ce n'est pas toujours une qualité...)

Au fait, j'ai fait des tests grandeur nature, mes deux codes sont faux. En fait, j'ai supposé à tort certains comportements d'overflow. J'ai corrigé le premier, reste à corriger le deuxième.
Ah, et au passage : dur dur de battre l'implémentation naïve dans les cas courants...
25 décembre 2011 à 18:31:20

Citation : yoch

En fait, j'ai supposé à tort certains comportements d'overflow.


Tu peux préciser le problème stp? j'ai la flemme de lire les corrections apportée, j'avoue! Sur des unsigned on a bien un modulo en cas d'overflow, non?

Citation : yoch

(je suis assez minimaliste effectivement, mais ce n'est pas toujours une qualité...)


Bah minimaliste sûrement :lol: , mais j'apprends toujours beaucoup avec tes codes! :)

edit: aie ortho
Zeste de Savoir, le site qui en a dans le citron !
25 décembre 2011 à 19:13:57

Citation : GurneyH

Tu peux préciser le problème stp? j'ai la flemme de lire les corrections apportée, j'avoue! Sur des unsigned on a bien un modulo en cas d'overflow, non?


En fait, j'ai travaillé sans utiliser de modulo, d’où mes soucis.

J'ai supposé que par exemple :
unsigned a = 1;
a << 100;
// ici, je pensais que a == 0, ce qui est faux...


Ensuite, pour mon second code (encore non corrigé), je pensais qu'un overflow de multiplication se comporte comme un modulo. Il me semble que ce n'est la cas qu'en base 2 (à confirmer que c'est bien le cas), mais pas dans les autres bases, donc on a bien besoin d'un modulo. (Au passage, même avec un modulo, il faut faire modulo une certaine puissance de la base pour avoir un résultat correct)
25 décembre 2011 à 19:33:34

Bonjours,
Je viens de voir le defis donc voici mon code pour le niveau 1
#include <stdio.h>
#include <stdlib.h>

int main()
{
    char gene[100] = {0}, sequenceADN[100] = {0};
    int i = 0, j = 0, verification = 0, taille[2] = {0}, indice = 0;

    printf("Entrez le gene voulu : ");
    scanf("%s", gene);
    printf("Entrez la sequence d'ADN etudier : ");
    scanf("%s", sequenceADN);

    for(i=0; gene[i]!='\0'; i++);
    taille[0] = i;
    for(i=0; sequenceADN[i]!='\0'; i++);
    taille[1] = 1;

    for(i=0; i<taille[1]; i++)
    {
        if(sequenceADN[i] == gene[0])
        {
            for(j=1; j<(taille[0] - 1); j++)
            {
                if(sequenceADN[i+j] == gene[j])
                    verification++;
            }
            if(verification == (taille[0] - 1) )
                indice = i + 1;
        }
    }

    if(indice != 0)
        printf("L'indice est : %d\n\n", indice);
    else
        printf("Le gene n'est pas dans la sequence d'ADN etudier\n\n");

    return 0;
}


Mais prise par une flemmingite aiguë très importante je ne ferais les deux autres niveaux que plus tard :)
26 décembre 2011 à 12:54:55

Salut,

Juste pour dire, il y a une petite faute dans l'énoncé:

printf("%ul\n", ULONG_MAX);


le "modificateur de longueur" comme "l" se positionne avant le "spécificateur de conversion". Autrement dit, c'est "%lu" et non "%ul" ;)

Sinon, voici ma participation pour le niveau 2 (pas super optimisé à mon avis :-° ):


#include <inttypes.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define GEN_MAX 18
#define SEQ_MAX 256
#define FMT_MAX 8


static void
delete_newline(char *s)
{
	for (; *s; s++)
		if (*s == '\n')
			break;
	*s = '\0';
}


static int
seqtonb(char *s)
{
	for (; *s; s++)
		switch(*s) {
		case 'A':
			*s = '0';
			break;

		case 'T':
			*s = '1';
			break;

		case 'G':
			*s = '2';
			break;

		case 'C' :
			*s = '3';
			break;

		default:
			return 0;
		}

	return 1;
}


static size_t
search_gen(char const *seq, char const *gen)
{
	static char fmt[FMT_MAX];
	uintmax_t nseq = 0, ngen = 0;
	size_t n;

	if (sscanf(gen, "%" SCNuMAX, &ngen) != 1)
		return 0;

	sprintf(fmt, "%%%zu" SCNuMAX, strlen(gen));
	for (n = 0; *seq && nseq != ngen; n++, seq++)
		if (sscanf(seq, fmt, &nseq) != 1)
			return 0;

	return (*seq) ? n : 0;
}


int
main(void)
{
	static char gen[GEN_MAX];
	static char seq[SEQ_MAX];
	size_t pos;

	if (fgets(gen, sizeof gen, stdin) == NULL ||
	    fgets(seq, sizeof seq, stdin) == NULL) {
	    	fputs("IO error\n", stderr);
	    	return EXIT_FAILURE;
	}

	delete_newline(gen);
	delete_newline(seq);
	if (!seqtonb(gen) || !seqtonb(seq)) {
		fputs("Invalid entry\n", stderr);
		return EXIT_FAILURE;
	}

	if ((pos = search_gen(seq, gen)) == 0) {
		fputs("No match\n", stderr);
		return EXIT_FAILURE;
	} else
		printf("%zu\n", pos);

	return EXIT_SUCCESS;
}



EDIT: petite modification.
26 décembre 2011 à 13:23:08

Oula, ça c'est vraiment original ! :lol:
26 décembre 2011 à 16:58:04

Citation : yoch


Oula, ça c'est vraiment original ! :lol:



Je ne sais pas comment je dois le prendre :p

Plus sérieusement, je voulais vraiment passé par sscanf pour la conversion chaine -> nombre, d'où le code obtenu (et la nécessité de crée un format de manière dynamique) :)
26 décembre 2011 à 17:04:26

Interresant je vais mis mettre

[FAIT][Défis] #9 : La recherche en génétique

× Après avoir cliqué sur "Répondre" vous serez invité à vous connecter pour que votre message soit publié.
× Attention, ce sujet est très ancien. Le déterrer n'est pas forcément approprié. Nous te conseillons de créer un nouveau sujet pour poser ta question.
  • Editeur
  • Markdown