/************************************************************************* * Last edited on 28/04/2005 by rcapua * * Last edited on 2005-05-27 18:35:27 by stolfi * * Last edited on 2005-06-09 15:59:00 by rcapua * * * * * * Universidade Federal Fluminense * * Mestrado em Computacao * * Aluna: Renatha Oliva Capua * * Orientadores: Helena Cristina da Gama Leitao * * Jorge Stolfi * * * * Pre.c - Programa com o preditor para o Reconhecimento de Exons * *************************************************************************/ #define PROG_NAME "pre" #define PROG_USAGE \ PROG_NAME " \\\n" \ " ARQ_EVENTOS \\\n" \ " ARQ_TABELA_1 ARQ_TABELA_2 ARQ_TABELA_3 \\\n" \ " TAM_MEDIO_CODING TAM_MEDIO_NONCODING \\\n" \ " ARQ_BASES [ ARQ_ROTULOS ] \\\n" \ " NOME_SAIDA \\\n" \ #define _GNU_SOURCE #include #include int main(int argc, char **argv){ int M; /* variavel que guarda o tamanho da tupla na sequencia alvo */ int NE; /* quantidade de eventos(e) */ tfreq *vfreq1; /* Freqs para cada par (1-evento,1-tupla) na base de treinamento */ tfreq *vfreq2; /* Freqs para cada par (2-evento,2-tupla) na base de treinamento */ tfreq *vfreq3; /* Freqs para cada par (3-evento,3-tupla) na base de treinamento */ int tam_vfreq1; /* Num linhas e num colunas da tabela vfreq1. */ int tam_vfreq2; /* Num linhas e num colunas da tabela vfreq2. */ int tam_vfreq3; /* Num linhas e num colunas da tabela vfreq3. */ char *bas; /* A seqüência de bases a analisar. */ int tam_bas; /* Tamanho da seqüência {bas}. */ char *gab; /* A classificação real das bases, ou NULL. */ int tam_gab; /* Tamanho da seqüência {gab}. */ double *probs; /* Vetor com probabilidades de E,F,G,I por base (4*tam_bas elementos) */ tevento *event; /* Eventos possíveis especificados pelo usuário. */ float tempoi_CPU, tempof_CPU; /* Variaveis para medir o tempo de CPU - tempo inicial e tempo final */ double tm_coding, tm_intron; /* tamanho medio de uma região codificadora e não-codificadora*/ /* para comecar a contar o tempo de execucao */ times(&tempoi); /* Pega o tempo de CPU */ tempoCPU(tempoi, &tempoi_CPU); /* para testar a quantidade de argumentos passados na chamada do programa */ if ((argc != 9) && (argc != 10)) { fprintf (stderr, "argc = %d\n", argc); fprintf (stderr, "use: %s\n", PROG_USAGE); return (-1); } char *evt_file = argv[1]; char *tb1_file = argv[2]; char *tb2_file = argv[3]; char *tb3_file = argv[4]; tm_coding = atof(argv[5]); tm_intron = atof(argv[6]); char *bas_file = argv[7]; char *rot_file = (argc == 10 ? argv[8] : NULL); char *sai_name = (argc == 10 ? argv[9] : argv[8]); /* Procedimento para ler o arquivo com os parametros de entrada do programa */ ler_param(evt_file, &NE, &M, &event, tm_coding, tm_intron); /* Para montar uma estrutura com as frequencias das triplas obtidas na fase de treinamento */ carregar_vfreq(tb1_file, 1, &vfreq1, &tam_vfreq1); carregar_vfreq(tb2_file, 2, &vfreq2, &tam_vfreq2); carregar_vfreq(tb3_file, 3, &vfreq3, &tam_vfreq3); /* Lê a sequência a analisar: */ fprintf(stderr, "Lendo o arquivo: %s\n", sai_name); ler_seq(bas_file, &tam_bas, &bas); if (rot_file != NULL) { ler_seq(rot_file, &tam_gab, &gab); if (tam_gab != tam_bas) { fprintf(stderr, "comprimentos bas = %d gab = %d são diferentes\n", tam_bas, tam_gab); exit (1); } } else { gab = NULL; tam_gab = 0; } probs = (double *) malloc(4*tam_bas*sizeof(double)); /* Calcula a probabilidade de cada base da entrada ser E,F,G ou I: */ analisar_bas( tam_bas, bas, vfreq1,tam_vfreq1, vfreq2,tam_vfreq2, vfreq3,tam_vfreq3, M,NE,event, probs ); /* Imprime cada base e as probabilidades estimadas de cada rótulo: */ char *est_file = NULL; asprintf(&est_file, "%s.est", sai_name); char *ert_file = NULL; asprintf(&ert_file, "%s.ert",sai_name); escrever_bas_probs( est_file, ert_file, tam_bas, bas, tam_gab, gab, probs ); if (rot_file != NULL) { int tipo; for (tipo = 0; tipo <= 1; tipo++) { double *escore_iefg = (double *) malloc (4*4 * sizeof(double)); /* I, E, F e G */ double *escore_ic = (double *) malloc (2*2 * sizeof(double)); /* I, C */ double *escore_efg = (double *) malloc (3*3 * sizeof(double)); /* E, F, G */ /* Calcula a quantidade de vezes que a predição feita para a sequencia obteve sucesso ou fracasso */ fprintf(stderr, "calculando escores... "); calcular_escore(&escore_iefg, &escore_ic, &escore_efg, probs, gab, tam_gab, tipo); /* Imprime o escore */ char *scr_file = NULL; char *tipo_ch = (tipo == 0 ? "int" : "flt"); asprintf(&scr_file, "%s-%s.scr", sai_name, tipo_ch); fprintf(stderr, "escrevendo escores em %s...\n", scr_file); escrever_escore(scr_file, escore_iefg, escore_ic, escore_efg, tipo); } } /* Finaliza o tempo de execucao */ times(&tempof); /* Tempo final */ tempoCPU(tempof, &tempof_CPU); fprintf(stderr, "\ntempo total de CPU(s): %.3f\n\n", ((tempof_CPU - tempoi_CPU)/60)); return(0); }//end main /**************************************************************************** * Procedimento para ler o arquivo com os parametros de entrada do programa. * * Os parametros de entrada sao NE(numero de eventos), M(o tamanho da tupla) * * e event(lista com os rotulos dos eventos e suas particoes). * ****************************************************************************/ void ler_param(char *file_param, int *NE, int *M, tevento **event, double tm_coding, double tm_intron){ int cnt_e; char letra; tevento *event_aux; /* lista para guardar os eventos possiveis */ FILE *pt_param; abrir_arq(file_param, &pt_param, "r"); fscanf(pt_param,"%c",&letra); consumir_cabecalho(&pt_param,&letra); ungetc(letra, pt_param); // Bota a letra de volta no arquivo. fscanf(pt_param,"%d", M); /* tamanho da tupla da sequencia alvo */ fscanf(pt_param,"%d", NE); /* quantidade de eventos */ /* Para incluir os rotulos dos eventos no vetor event */ event_aux = *event; event_aux = (tevento *)malloc((*NE) * sizeof(tevento)); /*cria o vetor do tipo tevento com NE elementos*/ char *temp = malloc ((2*(*M)) * sizeof(char)); for(cnt_e = 0; cnt_e < (*NE); cnt_e++){ tevento *ek = &(event_aux[cnt_e]); // int k; fscanf(pt_param, "%s", temp); extrai_evento_e_particao(temp, (*M), &ek->evento, &ek->num_particao, &ek->particao); ek->priori = get_prob_evento_a_priori(ek->evento, tm_coding, tm_intron); //fprintf(stderr, "%s %12.10f %3d", ek->evento, ek->priori, ek->num_particao); //for (k = 0; k < ek->num_particao; k++){ // fprintf(stderr, "%c%s", (k == 0 ? ' ' : '.'), ek->particao[k]); //} //fprintf(stderr, "\n"); }//fim do for *event = event_aux; fclose(pt_param); /*fecha o arquivo com os parametros*/ }//end ler_parametros /**************************************************************************** * Procedimento para criar o vetor que vai guardar as frequencias * ****************************************************************************/ int criar_vetor(int K, tfreq **vfreq, int tam_vfreq){ int cnt, cnt_prob; tfreq *vfreq_aux; vfreq_aux = *vfreq; vfreq_aux = (tfreq *) malloc(tam_vfreq * sizeof(tfreq)); if (vfreq_aux ==NULL){ printf("Erro na alocacao de memoria,vfreq nao pode ser criado!"); return(1); } for(cnt=0; cnt < tam_vfreq; cnt++){ vfreq_aux[cnt].tupla = (char *) malloc((K+1)* sizeof(char)); /* rotulo da tupla */ strcpy(vfreq_aux[cnt].tupla,"\0"); vfreq_aux[cnt].prob = (tprob *) malloc (tam_vfreq * sizeof(tprob));/*rotulo de cada evento e sua probabilidade*/ for (cnt_prob=0; cnt_prob < tam_vfreq; cnt_prob++){ vfreq_aux[cnt].prob[cnt_prob].evento = (char *) malloc((K+1) * sizeof(char)); /* rotulo do evento */ strcpy(vfreq_aux[cnt].prob[cnt_prob].evento,"\0"); /* limpa o valor do rotulo */ vfreq_aux[cnt].prob[cnt_prob].valor = 0.0; } }//fim do for i< vfreq_tam *vfreq = vfreq_aux; return(0); }//fim do criar_vetor /**************************************************************************** * Procedimento para liberar o espaco de memoria ocupado por vfreq * ****************************************************************************/ void desaloca_vfreq(tfreq **vfreq, int tam_vfreq){ int cnt, cnt_prob; tfreq *vfreq_aux; vfreq_aux = *vfreq; for (cnt=0; cnt < tam_vfreq;cnt++){ for (cnt_prob=0; cnt_prob < tam_vfreq; cnt_prob++) free(vfreq_aux[cnt].prob[cnt_prob].evento); free(vfreq_aux[cnt].prob); free(vfreq_aux[cnt].tupla); } free(vfreq_aux); } /**************************************************************************** * Procedimento que cria um vetor onde serao guardadas as frequencias com * * que as triplas aparecem nas estruturas de treinamento. Os valores sao * * lidos do arquivo e gravados no vetor vfreq para uma manipulacao mais * * rapida. * ****************************************************************************/ void carregar_vfreq(char *file_freq, int K, tfreq **vfreq, int *tam_vfreq){ char *tupla, /* rotulo da tupla */ *evento, /* rotulo do evento */ letra; int ja_cadastrada, /* testa se uma tupla ja' foi cadastrada */ ind, /* guarda o indice do ultimo valor de vfreq */ cnt, cnt_prob; int nocs; /* Número de ocorrências do par (rotulo,tupla) na amostra. */ double freq; /* Freqüência condicional F(tupla|rotulo). */ FILE *pt_freq; /* ponteiro para o arquivo com as frequencias */ tfreq *vfreq_aux; /* variavel auxiliar para facilitar a manipulacao com ponteiros */ ind = 0; freq = 0.0; tupla = (char*)malloc ((K + 1) * sizeof(char)); evento = (char*)malloc ((K + 1) * sizeof(char)); (*tam_vfreq) = (int)pow(4,K); /* tam = 4^k */ if (criar_vetor(K, vfreq, (*tam_vfreq))!= 0) exit(1); vfreq_aux = *vfreq; abrir_arq(file_freq,&pt_freq,"r"); /* Abre o arquivo com as frequencias */ fscanf(pt_freq,"%c",&letra); consumir_cabecalho(&pt_freq,&letra); ungetc(letra, pt_freq); // Bota a letra de volta no arquivo. fscanf(pt_freq, "%s", evento); fscanf(pt_freq, "%s", tupla); fscanf(pt_freq, "%d", &nocs); fscanf(pt_freq, "%lf", &freq); /* Para incluir os dados do arquivo no vetor criado */ while(!feof(pt_freq)){ /* Pesquisa nos elementos já exitentes em vfreq se a tupla foi cadastrada */ cnt = 0; /* contador para vfreq */ cnt_prob = 0; /* guardar a posicao vazia no vetor das prob, onde o evento deve ser inserido */ ja_cadastrada = 0; /* sinalizar se o elemento ja esta cadastrado */ while((cnt < ind) && (ja_cadastrada == 0)){ /*enquanto cnt_aux menor que quantidade de elem do vetor e nao achar*/ if(strcmp(vfreq_aux[cnt].tupla,tupla) == 0){ /* testa se eh igual*/ /* tenta encontrar uma posicao vazia no vetor*/ while((strcmp(vfreq_aux[cnt].prob[cnt_prob].evento,"\0") != 0) && (cnt_prob < (*tam_vfreq) - 1)) cnt_prob++; ja_cadastrada = 1; /*sinaliza que achou a tupla correspondente no vetor vfreq */ } cnt++; }//fim do while /* testa se a tupla jah foi realmente cadastrada(flag =1). Se foi, somente devemos incluir seu evento e sua probabilidade no vetor de eventos da tupla, caso contrario devemos incluir o rotulo da tupla tambem */ if (ja_cadastrada == 1){ strcpy(vfreq_aux[cnt - 1].prob[cnt_prob].evento,evento); vfreq_aux[cnt - 1].prob[cnt_prob].valor = freq; }else{ strcpy(vfreq_aux[ind].tupla,tupla); /*inclui o rotulo da tupla na ultima posicao(i) de vfreq*/ strcpy(vfreq_aux[ind].prob[cnt_prob].evento,evento);/*coloca rotulo do evento na primeira posicao de prob(0). Ja' que essa tupla nao era cadastrada(flag=0), prob nao continha nenhum elemento*/ vfreq_aux[ind].prob[cnt_prob].valor = freq; ind++; /* aumenta uma posicao no vetor */ }//fim if strcpy(tupla,"\0"); strcpy(evento,"\0"); freq = 0.0; fscanf(pt_freq,"%s",evento); fscanf(pt_freq,"%s",tupla); fscanf(pt_freq,"%d",&nocs); fscanf(pt_freq,"%lf",&freq); }//end while fclose(pt_freq); *vfreq = vfreq_aux; }//fim carregar_vfreq /**************************************************************************** * Procedimento que retorna a frequencia condicional F(tupla|evento), * * procurando na tabela {vfreq} de tamanho {tam_freq}. * * A tupla deve ter o mesmo tamanho {K} das entradas da tabela. * ****************************************************************************/ double get_freq_tupla(char *tupla, char *evento, int tam_vfreq, tfreq *vfreq){ double freq; /* variavel que guardara' a frequencia da tupla passada como parametro */ int cnt, ind; freq = 0.0; /* percorre o vfreq para pesquisar se a tupla esta' em vfreq */ for(cnt=0; cnt < tam_vfreq; cnt++){ if (strcmp(vfreq[cnt].tupla,tupla)==0){/* testa se e' igual a tupla passada como parametro */ for(ind=0; ind < tam_vfreq; ind++){ /* busca o evento correspondente ao passado por parametro */ if (strcmp(vfreq[cnt].prob[ind].evento,evento)==0) return(vfreq[cnt].prob[ind].valor); }//fim do for }//fim do if }//fim do for /* Se não achou, supõe que a freqüência é zero: */ return 0.0; } /**************************************************************************** * Calcula {ev_prob[i] = Pr(evento[i]|tupla), para cada evento valido * * {evento[i]} com i em 0..NE-1. Usa a fórmula de Bayes, obtendo * * Pr(tupla|evento[i]) a partir das tabelas e hipóteses de fatoração, e * * Pr(evento[i]) pelo modelo de Markov para transições intron-exon. * ****************************************************************************/ void calcular_ev_prob( int M, char *tupla, int tam_vfreq1, tfreq *vfreq1, int tam_vfreq2, tfreq *vfreq2, int tam_vfreq3, tfreq *vfreq3, int NE, tevento *event, double *ev_prob /* (Saida) probabilidades estimadas dos eventos. */ ){ int cnt, cnt_p, pos; double numerador; /* numerador para a formula da probabilidade */ double denom; char *tupla_p; /* variavel que guarda o valor da tupla correspondente a particao do evento*/ /*inicializacao das variaveis*/ tupla_p = (char*) malloc((M+1) * sizeof(char)); double freq_tupla; /* calcular a probabilidade para todos os eventos válidos */ denom = 0.0; for (cnt = 0; cnt < NE; cnt++){ pos = 0; numerador = 1.0; for (cnt_p = 0; cnt_p < event[cnt].num_particao; cnt_p++){ /* para cada evento pego as particoes */ /* copia o pedaco da tupla correspondente a particao do evento, comecando em ind e de tamanho strlen(event_aux[cnt].evento), e coloca na tupla_destino. */ strcpy(tupla_p,"\0"); str_incpy(tupla, tupla_p, pos, strlen(event[cnt].particao[cnt_p])); /* para o calculo do numerador devemos pegar a frequencia para cada particao */ switch(strlen(tupla_p)) { case 1: freq_tupla = get_freq_tupla(tupla_p, event[cnt].particao[cnt_p], tam_vfreq1, vfreq1); break; case 2: freq_tupla = get_freq_tupla(tupla_p, event[cnt].particao[cnt_p], tam_vfreq2, vfreq2); break; case 3: freq_tupla = get_freq_tupla(tupla_p, event[cnt].particao[cnt_p], tam_vfreq3, vfreq3); break; default: fprintf(stderr, "a parte \"%s\" tem comprimento inválido\n", event[cnt].particao[cnt_p]); exit(1); } numerador = numerador * freq_tupla; pos = pos + strlen(event[cnt].particao[cnt_p]); /* pos eh a posicao da proxima particao */ }//fim do for numerador = numerador * event[cnt].priori; ev_prob[cnt] = numerador; denom += numerador; }//fim do for for(cnt = 0; cnt < NE; cnt++) ev_prob[cnt] /= denom; }//fim do calcular_ev_prob /**************************************************************************** * Procedimento que retorna a probabilidade a priori de um evento. ****************************************************************************/ double get_prob_evento_a_priori(char *evento, double tm_coding, double tm_intron){ double prob; char *p = evento; /* Proxima letra do evento. */ int coding; /* 1 se letra anterior era EFG, 0 se era I. */ double probI = tm_intron/(tm_intron + tm_coding); /* Pr(janela inicia em I) */ double probC = tm_coding/(tm_intron + tm_coding); /* Pr(janela começa em E,F ou G. */ double probIC = 1.0/tm_intron; /* Pr(I ser seguido de E,F ou G). */ double probCI = 1.0/tm_coding; /* Pr(E,F ou G ser seguido de I). */ if ((*p) == 'I') { prob = probI; coding = 0; } else { prob = probC/3; coding = 1; } p++; while ((*p) != '\0'){ if (coding) { if ((*p) == 'I'){ prob *= probCI; coding = 0; }else{ prob *= (1.0 - probCI); coding = 1; } }else{ if ((*p) == 'I'){ prob *= (1.0 - probIC); coding = 0; }else{ prob *= probIC/3; coding = 1; } } p++; } return prob; }//fim do get_prob_evento_a_priori /**************************************************************************** * Procedimento que lê uma seqüência de bases ou rotulos do arquivo. * ****************************************************************************/ void ler_seq(char *nome_arq, int *tam, char **seq){ char letra; /* caracter(base) que e' lido do arquivo */ FILE *pt_seq; /* ponteiro para o arquivo com a sequencia alvo */ abrir_arq(nome_arq, &pt_seq, "r"); /* abrir arquivo com a sequencia alvo */ fscanf(pt_seq,"%c", &letra); printf("letra:%c\n",letra); consumir_cabecalho(&pt_seq, &letra); int na = 100; /* Tamanho alocado para a area (*seq) -- chute inicial. */ (*seq) = (char *)malloc(na*sizeof(char)); (*tam) = 0; /* Número de bases/rotulos já lidas. */ printf("Entrou ler_seq \n tam:%d\n ",(*tam)); /* Percorre a sequencia do arquivo */ while(!feof(pt_seq)){ str_chr_up(&letra); printf("B: %c - ", letra); if (((letra >= 'A') && (letra <='Z')) || (letra == '?') || (letra == '*')) { /* Testa se e' um caracter valido*/ if ((*tam) >= na) { /* Espaço esgotado, realocar sequência: */ int nr = 2*na; (*seq) = realloc((*seq), nr); if ((*seq) == NULL) { fprintf(stderr, "memoria esgotada"); } na = nr; } (*seq)[(*tam)] = letra; (*tam) += 1; } fscanf(pt_seq,"%c",&letra); }//fim while fclose(pt_seq); (*seq) = realloc((*seq), (*tam)+1); (*seq)[(*tam)] = '\0'; printf( "\ntam = %d\n", (*tam)); }//fim ler_seq /**************************************************************************** * Procedimento que percorre a seqüência de bases a analisar {bas} e coloca * * em {probs[4*i+k]} a probabilidade estimada de que a base {bas[i]} tenha * * rótulo {k} (0=I, 1=E, 2=F, 3=G). * ****************************************************************************/ void analisar_bas( int tam_bas, char *bas, tfreq *vfreq1, int tam_vfreq1, tfreq *vfreq2, int tam_vfreq2, tfreq *vfreq3, int tam_vfreq3, int M, int NE, tevento *event, double *probs ){ char letra, /* caracter(base) da sequência */ *tupla; int cnt, pos=0; /* posicao da base na seqüência */ int i,k; int meio; /* Indice da letra do meio da janela. */ int qual; /* Letra de {bas} que está sendo classificada. */ double *ev_prob; /* vetor que guardara' os eventos da tupla e a probabilidade */ tupla = (char*)malloc((M + 1) * sizeof(char)); /* inicializacao das variaveis */ strcpy(tupla,"\0"); letra = '\0'; cnt = 0; meio = M/2; // char letra, /* caracter(base) da sequência */ // *tupla; // int cnt, // meio = M/2; ev_prob = (double *) malloc (NE*sizeof(double)); /* Preenche início do vetor {probs} com zeros: */ for (i = 0; i < meio; i++) for (k = 0; k < 4; k++) { probs[4*i+k] = 0.0; } /* Percorre a sequencia do arquivo */ for (pos = 0; pos < tam_bas; pos++){ letra = bas[pos]; /* pega um caracter(uma base) do arquivo com a sequencia alvo */ str_ncat(tupla,letra,M); /* acrescenta o caracter ao rotulo da tupla */ /* Se o conjunto de caracteres dentro de tupla jah chegou no tamanho que a tupla deve ter, que eh M. Devemos calcular sua probabilidade para cada um dos eventos considerados. Essa condicao nao eh atendida nos primeiros caracteres lidos do arquivo */ if (strlen(tupla) == (M)){ /* Se tupla contem N ou outros caracteres nao devemos usa'-la */ if (caracter_valido(tupla)==0){ /* calcula a probabilidade condicional para cada evento da tupla */ calcular_ev_prob( M, tupla, tam_vfreq1,vfreq1, tam_vfreq2,vfreq2, tam_vfreq3,vfreq3, NE,event, ev_prob ); /* Condensa as probabilidades dos eventos segundo a letra do meio: */ qual = pos - (M - 1) + meio; condensar_prob(NE,event,ev_prob, meio, &(probs[4*qual])); } }//fim do if }//fim for /* Preenche final do vetor {probs} com zeros: */ for (i = tam_bas - (M-1-meio); i < tam_bas; i++) for (k = 0; k < 4; k++) { probs[4*i+k] = 0.0; } }//fim analisar_bas /**************************************************************************** * Coloca em {probs[k]} a soma das probabilidades {ev_prob[j]} de todos os * * eventos {event[j]} que possuem o mesmo rotulo {k} na posição {meio}, * * para cada {k} (0=I, 1=E, 2=F, 3=G). * ****************************************************************************/ void condensar_prob(int NE, tevento *event, double *ev_prob, int meio, double *probs){ int cnt; int k; for (k = 0; k < 4; k++) { probs[k] = 0; } /* Inicializa probs com zero */ for (cnt = 0; cnt < NE; cnt++) { char *evento = event[cnt].evento; int irot; switch(evento[meio]){ case 'I': irot = 0; break; case 'E': irot = 1; break; case 'F': irot = 2; break; case 'G': irot = 3; break; default: fprintf(stderr, "condensar_prob: evento = %s meio = %c\n", evento, evento[meio]); irot = 0; }//fim do switch probs[irot] += ev_prob[cnt]; }//fim do for }//fim do condensar_prob /**************************************************************************** * Devolve o número de caracteres '.' na cadeia {s}. * ****************************************************************************/ int conta_pontos(char *s){ int n = 0; while((*s) != '\0') { if ((*s) == '.') { n++; } s++; } return n; }//fim do conta_pontos /**************************************************************************** * Excreve um cabeçalho apropriado para o arquivo de saida. Se {hifens} é * * verdadeiro, imprime a linha de hifens correspondente. * ****************************************************************************/ void escrever_cabecalho(FILE *pt_arq, int com_gab, int hifens){ fprintf(pt_arq, ">%11s", (hifens ? "----------" : "posição")); fprintf(pt_arq, " %1s", (hifens ? "-" : "b")); fprintf(pt_arq, " "); fprintf(pt_arq, " %-9s", (hifens ? "---------" : " Pr(I)")); fprintf(pt_arq, " %-9s", (hifens ? "---------" : " Pr(E)")); fprintf(pt_arq, " %-9s", (hifens ? "---------" : " Pr(F)")); fprintf(pt_arq, " %-9s", (hifens ? "---------" : " Pr(G)")); if (com_gab) { fprintf(pt_arq, " %4s", (hifens ? "----" : "rots")); fprintf(pt_arq, " "); fprintf(pt_arq, " %-9s", (hifens ? "---------" : " Pr(I)")); fprintf(pt_arq, " %-9s", (hifens ? "---------" : " Pr(C)")); fprintf(pt_arq, " %4s", (hifens ? "----" : "clas")); fprintf(pt_arq, " "); fprintf(pt_arq, " %4s", (hifens ? "----" : "fase")); } else { fprintf(pt_arq, " %1s", (hifens ? "-" : "r")); fprintf(pt_arq, " "); fprintf(pt_arq, " %-9s", (hifens ? "---------" : " Pr(I)")); fprintf(pt_arq, " %-9s", (hifens ? "---------" : " Pr(C)")); fprintf(pt_arq, " %1s", (hifens ? "-" : "c")); fprintf(pt_arq, " "); fprintf(pt_arq, " %1s", (hifens ? "-" : "ø")); } fprintf(pt_arq,"\n"); }//fim do escrever_cabecalho /**************************************************************************** * Escreve o arquivo de saida, dadas as seqüências de bases {bas} e dos * * rótulos reais {gab}, e a probabilidade estimada {probs[4*i+k]} de * * {bas[i]} ter rótulo {k} (0=I, 1=E, 2=F, 3=G). * ****************************************************************************/ void escrever_bas_probs( char *est_file, char *ert_file, int tam_bas, char *bas, int tam_gab, char *gab, double *probs ){ FILE *est_pt; /* ponteiro para o arquivo da saida completa do preditor */ FILE *ert_pt; /* ponteiro para o arquivo de rótulos estimados pelo preditor */ int i; abrir_arq(est_file,&est_pt,"w+"); /* abrir arquivo com a saida do preditor */ abrir_arq(ert_file,&ert_pt,"w+"); /* abrir arquivo com rotulacao do preditor */ /* imprime cabeçalhos: */ escrever_cabecalho(est_pt, (gab != NULL), 0); escrever_cabecalho(est_pt, (gab != NULL), 1); int n_ert = 0; /* Numero de rótulos escritos em {ert_pt} */ for (i = 0; i < tam_bas; i++) { double *pr = &(probs[4*i]); /* Probabilidades estimadas pelo programa para {I,E,F,G}. */ int k; char letra = bas[i]; char rot_pre; /* Rotulo {IEFG} mais provável segundo programa, '?' se não sabe. */ int irot_pre; /* {rot_pre} convertido de {?,I,E,F,G} para {0,0,1,2,3}. */ double pr_tot, prI, prC; char classe_pre;/* Classe estimada pelo programa {I,C,?}. */ char fase_pre; /* Fase {E,F,G} estimada pelo programa, '?' se não sabe. */ int ifase_pre; /* Fase {?,E,F,G} convertida para {0,1,2,3}. */ pr_tot = pr[0] + pr[1] + pr[2] + pr[3]; fprintf(est_pt," %11d", i); fprintf(est_pt," %c", letra); fprintf(est_pt, " "); if (pr_tot == 0) { irot_pre = 0; rot_pre = '?'; } else { /* determina o rótulo mais provável: */ irot_pre = 0; for (k = 1; k < 4; k++) { if (pr[k] > pr[irot_pre]) { irot_pre = k; } } switch(irot_pre){ case 0: rot_pre = 'I'; break; case 1: rot_pre = 'E'; break; case 2: rot_pre = 'F'; break; case 3: rot_pre = 'G'; break; default: fprintf(stderr, "irot_pre inválido\n"); rot_pre = '?'; exit(1); }//fim do switch }//fim do else if for (k = 0; k < 4; k++) { /* Assinala a maior probabilidade: */ char marca = ((pr_tot > 0) && (k == irot_pre) ? '*' : ' '); fprintf(est_pt," %c%8.6lf", marca, pr[k]); }//fim do for fprintf(est_pt," %c", rot_pre); if (gab != NULL) { char rot_gab = gab[i]; /* Rotulo real. */ int irot_gab; /* {rot_gab} convertido de {I,E,F,G} para {0,1,2,3}. */ char rot_ok; /* '·' se rótulos batem, '!' se não batem, '?' se não sabe. */ /* Converte rotulo real para inteiro: */ switch(rot_gab){ case 'I': irot_gab = 0; break; case 'E': irot_gab = 1; break; case 'F': irot_gab = 2; break; case 'G': irot_gab = 3; break; default: fprintf(stderr, "rot_gab inválido\n"); irot_gab = -1; exit(1); } if ((rot_pre == '?') || (rot_gab == '?')) { rot_ok = '?'; } else { rot_ok = (rot_pre == rot_gab ? '·' : '*'); } fprintf(est_pt,":%c%c", rot_gab, rot_ok); }//fim do if fprintf(est_pt, " "); /* probabilidades das classes {I,C}: */ prI = pr[0]; prC = pr[1]+pr[2]+pr[3]; fprintf(est_pt," %c%8.6lf", ' ', prI); fprintf(est_pt," %c%8.6lf", ' ', prC); /* classe mais provável: */ if (prI + prC == 0.0) { classe_pre = '?'; } else { classe_pre = (prI + prC == 0 ? '?' : (prI > prC ? 'I' : 'C')); } fprintf(est_pt," %c", classe_pre); if (gab != NULL) { char rot_gab = gab[i]; /* Rotulo real. */ char classe_gab = (rot_gab == 'I' ? 'I' : 'C'); /* Classe real {I,C}. */ char classe_ok; /* '·' se classes batem, '!' se não batem, '?' se não sabe. */ if ((classe_pre == '?') || (classe_gab == '?')) { classe_ok = '?'; } else { classe_ok = (classe_pre == classe_gab ? '·' : '*'); } fprintf(est_pt,":%c%c", classe_gab, classe_ok); }//fim do if fprintf(est_pt, " "); /* fase mais provável: */ if ((pr[1] >= pr[2]) && (pr[1] >= pr[3])) { ifase_pre = 1; } else if ((pr[2] > pr[1]) && (pr[2] >= pr[3])) { ifase_pre = 2; } else { ifase_pre = 3; } if (pr[ifase_pre] == 0.0) { ifase_pre = 0; } switch (ifase_pre){ case 0: fase_pre = '?'; break; case 1: fase_pre = 'E'; break; case 2: fase_pre = 'F'; break; case 3: fase_pre = 'G'; break; default: fprintf(stderr, "ifase_pre inválida\n"); fase_pre = '?'; exit(1); }//fim do switch fprintf(est_pt," %c", fase_pre); if (gab != NULL) { char rot_gab = gab[i]; /* Rotulo real. */ char fase_gab; /* Fase real {E,F,G}ou '?' se {I}. */ int ifase_gab; /* Fase real {1,2,3} ou 0 se {I}. */ char delta; /* Diferenca de fase {·,+,-,?}. */ /* Extrai a fase real: */ switch(rot_gab){ case 'I': ifase_gab = 0; break; case 'E': ifase_gab = 1; break; case 'F': ifase_gab = 2; break; case 'G': ifase_gab = 3; break; default: fprintf(stderr, "rot_gab inválido\n"); ifase_gab = -1; exit(1); }//fim do switch(rot_gab) fase_gab = (rot_gab == 'I' ? '?' : rot_gab); /* Calcula a defasagem: */ if ((ifase_pre == 0) || (ifase_gab == 0)) { delta = '?'; } else { int idelta = (ifase_pre - ifase_gab + 3) % 3; switch (idelta){ case 0: delta = '·'; break; case 1: delta = '+'; break; case 2: delta = '-'; break; default: fprintf(stderr, "idelta inválido\n"); delta = '?'; exit(1); } //fim do switch } fprintf(est_pt,":%c%c", fase_gab, delta); }//fim do if(gab) fprintf(est_pt,"\n"); /* Escreve rótulo do preditor: */ if ((n_ert > 0) && (n_ert % 60 == 0)) { fprintf(ert_pt, "\n"); } if (pr_tot == 0) { fprintf(ert_pt, "?"); } else if (prI > 0.800) { fprintf(ert_pt, "I"); } else if (prC > 0.800) { fprintf(ert_pt, "%c", fase_pre); } else { fprintf(ert_pt, "?"); } n_ert++; if ((pr_tot != 0.0) && (fabs(pr_tot - 1) > 0.00000001)) { fprintf(stderr, "probabilidades não somam 0 ou 1\n"); } }//fim do for fclose(est_pt); fprintf(ert_pt, "\n"); fclose(ert_pt); }// fim do escrever_bas_probs /**************************************************************************** * Dado um evento pontuado, como "II.III.FG.EFG.E", extrai do mesmo o evento * * puro "IIIIIFGEFGE" e as partes "II", "III", "FG", "EFG", "E", guardadas * * em {parte[0..npartes-1]}. * ****************************************************************************/ void extrai_evento_e_particao(char *ev_pontuado, int M, char **evento, int *npartes, char ***parte) { (*evento) = (char *) malloc ((M + 1) * sizeof(char)); (*npartes) = conta_pontos(ev_pontuado) + 1; (*parte) = (char **)malloc((*npartes) * sizeof(char*)); int ip; char *p = ev_pontuado; char *r = (*evento); for(ip = 0; ip < ((*npartes)); ip++){ char *q = p; while (((*q) != '.') && ((*q) != '\0')) { (*r) = (*q); r++; q++; } int tam_parte = q - p; (*q) = '\0'; (*parte)[ip] = (char *)malloc((tam_parte + 1)*sizeof(char)); strcpy((*parte)[ip], p); p = q+1; } (*r) = '\0'; }//fim do extrai_evento_e_particao /************************************************************************* * Percorre o vetor com as probabilidades encontradas e conta quantas * * vezes o preditor errou ou acertou ao rotular uma base. Monta um vetor, * * cada posição é determinada por duas variaveis: linha e coluna. A linha * * indica o evento verdadeiro da base e a coluna o evento previsto pelo * * preditor. * * Se {tipo==0}, arredonda as probabilidades para um-1-e-o-resto-0 * * Se {tipo==1}, soma as probabilidades sem arredondar. * *************************************************************************/ void calcular_escore(double **escore_iefg, double **escore_ic, double **escore_efg, double *probs, char *gab, int tam_gab, int tipo){ int cnt, i; int lin_iefg = 0, lin_ic = 0, lin_efg = 0; int col_iefg = 0, col_ic = 0, col_efg = 0; /* Inicializa escores */ for(cnt = 0; cnt < 16; cnt++){ (*escore_iefg)[cnt] = 0; } for(cnt = 0; cnt < 4; cnt++){ (*escore_ic)[cnt] = 0; } for(cnt = 0; cnt < 9; cnt++){ (*escore_efg)[cnt] = 0; } if (gab != NULL){ for(i = 0; i < tam_gab; i++){ switch(gab[i]){ /* A linha indica o gabarito para a base */ case 'I': lin_iefg = 0; lin_ic = 0; break; case 'E': lin_iefg = 1; lin_ic = 1; lin_efg = 0; break; case 'F': lin_iefg = 2; lin_ic = 1; lin_efg = 1; break; case 'G': lin_iefg = 3; lin_ic = 1; lin_efg = 2; break; }//fim do switch double prob_i = probs[4 * i + 0]; double prob_c = probs[4 * i + 1] + probs[4 * i + 2] + probs[4 * i + 3]; /* soma = Pr(C) = Pr(E) + Pr(F) + Pr(G) */ if (prob_i + prob_c != 0.0){ /* Tira as primeiras bases que não foram rotuladas */ if (tipo == 0) { /* Determina a coluna para as tabelas {escore_iefg} e {escore_efg}: */ int col_iefg = 0; int col_efg = 1; for (cnt = 1; cnt < 4; cnt++){ if (probs[4 * i + cnt] > probs[4 * i + col_iefg]) { col_iefg = cnt; } if (probs[4 * i + cnt] > probs[4 * i + col_efg]) { col_efg = cnt; }} (*escore_iefg)[4 * lin_iefg + col_iefg] += 1; /* Determina a coluna para a tabela {escore_ic}: */ col_ic = (prob_i > prob_c ? 0 : 1); (*escore_ic)[2 * lin_ic + col_ic] += 1; /* Determina se o evento tem janela para a tabela {escore_efg}: */ if (lin_ic != 0) { (*escore_efg)[3 * lin_efg + (col_efg - 1)] += 1; } } else { for (col_iefg = 0; col_iefg < 4; col_iefg++){ (*escore_iefg)[4 * lin_iefg + col_iefg] += probs[4 * i + col_iefg]; } if ((lin_ic != 0) && (prob_c != 0)) { for(col_efg = 1; col_efg < 4; col_efg++){ /* Normalizar o valor das probabilidades de {escore_efg}: */ (*escore_efg)[3 * lin_efg + (col_efg - 1)] += (probs[4 * i + col_efg]/ prob_c); } } (*escore_ic)[2 * lin_ic + 0] += prob_i; (*escore_ic)[2 * lin_ic + 1] += prob_c; }//fim do if else } }//fim do for }//fim do if gab }//fim do procedimento /**************************************************************************** * Imprime o cabecalho para o arquivo de escore(avaliador) * ****************************************************************************/ void escrever_cabecalho_escore(FILE *pt_arq, int linha, int tipo_escore){ fprintf(pt_arq, " "); if (tipo_escore != 3){ /* Diferente de {escore_efg}: */ fprintf(pt_arq, " %-9s", (linha ? "---------" : " Est(I)")); } if (tipo_escore != 0){ /* Diferente de {escore_ic} */ fprintf(pt_arq, " %-9s", (linha ? "---------" : " Est(E)")); fprintf(pt_arq, " %-9s", (linha ? "---------" : " Est(F)")); fprintf(pt_arq, " %-9s", (linha ? "---------" : " Est(G)")); }else { fprintf(pt_arq, " %-9s", (linha ? "---------" : " Est(C)")); } fprintf(pt_arq, "\n"); }//fim do escrever_cabecalho_escore void escrever_escore(char *nome_arq, double *escore_iefg, double *escore_ic, double *escore_efg, int tipo) { FILE *pt_escore; abrir_arq(nome_arq, &pt_escore, "w+"); escrever_matriz_iefg(&pt_escore, escore_iefg, tipo); fprintf(pt_escore, "\n\t Pr(C) = Pr(E) + Pr(F) + Pr(G)\n\n"); escrever_matriz_ic(&pt_escore, escore_ic, tipo); escrever_matriz_efg(&pt_escore, escore_efg, tipo); fclose(pt_escore); }//fim do escrever_escore void escrever_matriz_efg(FILE **pt, double *matriz, int tipo){ FILE *pt_escore; int lin, col; char rot[6]; double soma = 0, soma_diag = 0; pt_escore = *pt; escrever_cabecalho_escore(pt_escore, 0, 3); escrever_cabecalho_escore(pt_escore, 1, 3); fflush(pt_escore); strcpy(rot, "\0"); for (lin = 0; lin < 3; lin++) { switch (lin){ case 0 : strcpy(rot, " Pr(E)"); break; case 1 : strcpy(rot, " Pr(F)"); break; case 2 : strcpy(rot, " Pr(G)"); break; } fprintf(pt_escore, " %-9s", rot); for ( col = 0; col < 3; col++) { double elem = matriz[3 * lin + col]; if (tipo == 1) { fprintf(pt_escore, " %9.2f", elem); } else { fprintf(pt_escore, " %9d", (int)elem); } soma += elem; if (col == lin) soma_diag += elem; } fprintf(pt_escore, "\n\n"); }//fim do for fprintf(pt_escore, "\n\n"); if (soma != 0) { fprintf(pt_escore, "acertos = %.2f/%.2f = %.3f\n", soma_diag, soma, soma_diag/soma); } else { fprintf(pt_escore, "acertos = - "); } fprintf(pt_escore, "\n\n"); *pt = pt_escore; } /**************************************************************************** * Imprime o arquivo com o resultado do escore(avaliador) * ****************************************************************************/ void escrever_matriz_iefg(FILE **pt, double *matriz, int tipo){ FILE *pt_escore; int lin, col; char rot[6]; double soma_total = 0, soma_diag_total = 0; double soma_cod = 0, soma_diag_cod = 0; pt_escore = *pt; /* (I) X (E,F,G) */ escrever_cabecalho_escore(pt_escore, 0, 1); escrever_cabecalho_escore(pt_escore, 1, 1); fflush(pt_escore); strcpy(rot, "\0"); for(lin = 0; lin < 4; lin++){ switch (lin){ case 0 : strcpy(rot, " Pr(I)"); break; case 1 : strcpy(rot, " Pr(E)"); break; case 2 : strcpy(rot, " Pr(F)"); break; case 3 : strcpy(rot, " Pr(G)"); break; } fprintf(pt_escore, " %-9s", rot ); for( col = 0; col < 4; col++){ double elem = matriz[4 * lin + col]; if (tipo == 1) { fprintf(pt_escore, " %9.2f", elem); } else { fprintf(pt_escore, " %9d", (int)elem); } soma_total += elem; if (col == lin) soma_diag_total += elem; if (lin > 0) { soma_cod += elem; if (col == lin) soma_diag_cod += elem; } } fprintf(pt_escore, "\n\n"); }//fim do for fprintf(pt_escore, "\n\n"); fprintf(pt_escore, "acertos - total = %.2f/%.2f = %.3f\n", soma_diag_total, soma_total, soma_diag_total/soma_total); if (soma_cod != 0) { fprintf(pt_escore, "acertos - codificadora = %.2f/%.2f = %3f \n", soma_diag_cod, soma_cod, soma_diag_cod/soma_cod); } else { fprintf(pt_escore, "acertos - codificadora = - "); } fprintf(pt_escore, "\n\n"); *pt = pt_escore; }//fim do escrever_matriz /**************************************************************************** * Imprime o arquivo com o resultado do escore(avaliador) * ****************************************************************************/ void escrever_matriz_ic(FILE **pt, double *matriz, int tipo){ FILE *pt_escore; int lin, col; char rot[6]; double soma_total = 0, soma_diag_total = 0; double soma_cod = 0, soma_diag_cod = 0; pt_escore = *pt; /* (C) X (I) */ escrever_cabecalho_escore(pt_escore, 0, 0); escrever_cabecalho_escore(pt_escore, 1, 0); strcpy(rot, "\0"); for(lin = 0; lin < 2; lin++){ switch (lin){ case 0 : strcpy(rot, " Pr(I)"); break; case 1 : strcpy(rot, " Pr(C)"); break; } fprintf(pt_escore, " %-9s", rot ); for( col = 0; col < 2; col++){ double elem = matriz[2 * lin + col]; if (tipo == 1) { fprintf(pt_escore, " %9.2f", elem); } else { fprintf(pt_escore, " %9d", (int)elem); } soma_total += elem; if (col == lin) soma_diag_total += elem; if (lin > 0 ) { soma_cod += elem; if (col == lin) soma_diag_cod += elem; } } fprintf(pt_escore, "\n\n"); }//fim do for fprintf(pt_escore, "\n\n"); fprintf(pt_escore, "acertos - total = %.2f/%.2f = %.3f\n", soma_diag_total, soma_total, soma_diag_total/soma_total); if (soma_cod != 0) { fprintf(pt_escore, "acertos - codificadora = %.2f/%.2f = %3f \n", soma_diag_cod, soma_cod, soma_diag_cod/soma_cod); } else { fprintf(pt_escore, "acertos - codificadora = - "); } fprintf(pt_escore, "\n\n"); *pt = pt_escore; }//fim do escrever_matriz