KdJametPython/neuron_kd_switch3.c

#include <stdio.h>
#include <math.h>
#include <errno.h>
#include <stdlib.h>
#include <string.h>
#include <memory.h>
#include "neuron_switch3.h"
#include "iop_Rrs_neuron.h"

/* ------------------------------------------------------------------------------- 
 Lecture des differrentes Look Up Tables
 -------------------------------------------------------------------------------  */
void neuron_lect_LUTswitch(int idSensor, char *lutpath)
{
  FILE *fic;
  int i,j,poub;
  char *ligne=malloc(sizeof(char)*150), nomfic[1000];
  float fpoub;

  if( (poub = strcmp(lutpath, "None")) == 0){
    if( (lutpath = getenv("IOP_LUTS_PATH")) == NULL) {perror("IOP_LUTS_PATH"); exit(-1);}
  }
  /* ===================================== SeaWiFS ===================================== */
  if( idSensor == idSEAWIFS || idSensor == idMERIS){
    /* ----- LUTs for Rrs490/Rrs555 >= .85 ------ */
    sprintf(nomfic,"%s/%s", lutpath, rsupSW_LUT_POIDS);
    if( (fic=fopen(nomfic,"r")) == NULL) {perror(nomfic); exit(-1);}
    fgets(ligne,150,fic);
    for(i=0; i<rsupSW_NC1; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rsupSW_b1[i]);
    for(i=0; i<rsupSW_NC2; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rsupSW_b2[i]);
    fscanf(fic,"%d %d %f",&poub,&poub,&rsupSW_b3);
    for(j=0; j<rsupSW_NE; j++){
    for(i=0; i<rsupSW_NC1; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rsupSW_w1[j][i]);
    }
    for(j=0; j<rsupSW_NC1; j++){
      for(i=0; i<rsupSW_NC2; i++)
        fscanf(fic,"%d %d %f",&poub,&poub,&rsupSW_w2[j][i]);
    }
    for(i=0; i<rsupSW_NC2; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rsupSW_w3[i]);
    fclose(fic);

    sprintf(nomfic,"%s/%s", lutpath, rsupSW_LUT_MOY);
    if( (fic=fopen(nomfic,"r")) == NULL) {perror(nomfic); exit(-1);}
    fscanf(fic,"%f",&fpoub);
    for(i=0; i<rsupSW_NE-1; i++)
      fscanf(fic,"%f",&rsupSW_moy[i]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&rsupSW_moy[rsupSW_NE-1]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&rsupSW_moy[rsupSW_NES-1]);
    fclose(fic);
    sprintf(nomfic,"%s/%s", lutpath, rsupSW_LUT_ECART);
    if( (fic=fopen(nomfic,"r")) == NULL) {perror(nomfic); exit(-1);}
    fscanf(fic,"%f",&fpoub);
    for(i=0; i<rsupSW_NE-1; i++)
      fscanf(fic,"%f",&rsupSW_ecart[i]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&rsupSW_ecart[rsupSW_NE-1]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&rsupSW_ecart[rsupSW_NES-1]);
    fclose(fic);

  /* ----- LUTs for Rrs490/Rrs555 < .85 ------ */
    sprintf(nomfic,"%s/%s", lutpath, rinfSW_LUT_POIDS);
    if( (fic=fopen(nomfic,"r")) == NULL) {perror(nomfic); exit(-1);}
    fgets(ligne,150,fic);
    for(i=0; i<rinfSW_NC1; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rinfSW_b1[i]);
    for(i=0; i<rinfSW_NC2; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rinfSW_b2[i]);
    fscanf(fic,"%d %d %f",&poub,&poub,&rinfSW_b3);
    for(j=0; j<rinfSW_NE; j++){
    for(i=0; i<rinfSW_NC1; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rinfSW_w1[j][i]);
    }
    for(j=0; j<rinfSW_NC1; j++){
      for(i=0; i<rinfSW_NC2; i++)
        fscanf(fic,"%d %d %f",&poub,&poub,&rinfSW_w2[j][i]);
    }
    for(i=0; i<rinfSW_NC2; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rinfSW_w3[i]);
    fclose(fic);

    sprintf(nomfic,"%s/%s", lutpath, rinfSW_LUT_MOY);
    if( (fic=fopen(nomfic,"r")) == NULL) {perror(nomfic); exit(-1);}
    fscanf(fic,"%f",&fpoub);
    for(i=0; i<rinfSW_NE; i++)
      fscanf(fic,"%f",&rinfSW_moy[i]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&rinfSW_moy[rinfSW_NES-1]);
    fclose(fic);
    sprintf(nomfic,"%s/%s", lutpath, rinfSW_LUT_ECART);
    if( (fic=fopen(nomfic,"r")) == NULL) {perror(nomfic); exit(-1);}
    fscanf(fic,"%f",&fpoub);
    for(i=0; i<rinfSW_NE; i++)
      fscanf(fic,"%f",&rinfSW_ecart[i]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&rinfSW_ecart[rinfSW_NES-1]);
    fclose(fic);
  }
  
  /* =====================================  MODIS  ===================================== */
  if( idSensor == idMODIS){
    /* ----- LUTs for Rrs490/Rrs555 >= .85 ------ */
    sprintf(nomfic,"%s/%s", lutpath, rsupMO_LUT_POIDS);
    if( (fic=fopen(nomfic,"r")) == NULL) {perror(nomfic); exit(-1);}
    fgets(ligne,150,fic);
    for(i=0; i<rsupMO_NC1; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rsupMO_b1[i]);
    for(i=0; i<rsupMO_NC2; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rsupMO_b2[i]);
    fscanf(fic,"%d %d %f",&poub,&poub,&rsupMO_b3);
    for(j=0; j<rsupMO_NE; j++){
    for(i=0; i<rsupMO_NC1; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rsupMO_w1[j][i]);
    }
    for(j=0; j<rsupMO_NC1; j++){
      for(i=0; i<rsupMO_NC2; i++)
        fscanf(fic,"%d %d %f",&poub,&poub,&rsupMO_w2[j][i]);
    }
    for(i=0; i<rsupMO_NC2; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rsupMO_w3[i]);
    fclose(fic);

    sprintf(nomfic,"%s/%s", lutpath, rsupMO_LUT_MOY);
    if( (fic=fopen(nomfic,"r")) == NULL) {perror(nomfic); exit(-1);}
    fscanf(fic,"%f",&fpoub);/*412*/
    fscanf(fic,"%f",&rsupMO_moy[0]); /*443*/
    fscanf(fic,"%f",&rsupMO_moy[1]); /*488*/
    fscanf(fic,"%f",&rsupMO_moy[2]); /*531*/
    fscanf(fic,"%f",&rsupMO_moy[3]); /*547*/
    fscanf(fic,"%f",&fpoub); /*620*/
    fscanf(fic,"%f",&fpoub); /*667*/
    fscanf(fic,"%f",&rsupMO_moy[rsupMO_NE-1]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&rsupMO_moy[rsupMO_NES-1]);
    fclose(fic);
    sprintf(nomfic,"%s/%s", lutpath, rsupMO_LUT_ECART);
    if( (fic=fopen(nomfic,"r")) == NULL) {perror(nomfic); exit(-1);}
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&rsupMO_ecart[0]); /*443*/
    fscanf(fic,"%f",&rsupMO_ecart[1]); /*488*/
    fscanf(fic,"%f",&rsupMO_ecart[2]); /*531*/
    fscanf(fic,"%f",&rsupMO_ecart[3]); /*547*/
    fscanf(fic,"%f",&fpoub); /*620*/
    fscanf(fic,"%f",&fpoub); /*667*/
    fscanf(fic,"%f",&rsupMO_ecart[rsupMO_NE-1]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&rsupMO_ecart[rsupMO_NES-1]);
    fclose(fic);

  /* ----- LUTs for Rrs490/Rrs555 < .85 ------ */
    sprintf(nomfic,"%s/%s", lutpath, rinfMO_LUT_POIDS);
    if( (fic=fopen(nomfic,"r")) == NULL) {perror(nomfic); exit(-1);}
    fgets(ligne,150,fic);
    for(i=0; i<rinfMO_NC1; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rinfMO_b1[i]);
    for(i=0; i<rinfMO_NC2; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rinfMO_b2[i]);
    fscanf(fic,"%d %d %f",&poub,&poub,&rinfMO_b3);
    for(j=0; j<rinfMO_NE; j++){
    for(i=0; i<rinfMO_NC1; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rinfMO_w1[j][i]);
    }
    for(j=0; j<rinfMO_NC1; j++){
      for(i=0; i<rinfMO_NC2; i++)
        fscanf(fic,"%d %d %f",&poub,&poub,&rinfMO_w2[j][i]);
    }
    for(i=0; i<rinfMO_NC2; i++)
      fscanf(fic,"%d %d %f",&poub,&poub,&rinfMO_w3[i]);
    fclose(fic);

    sprintf(nomfic,"%s/%s", lutpath, rinfMO_LUT_MOY);
    if( (fic=fopen(nomfic,"r")) == NULL) {perror(nomfic); exit(-1);}
    fscanf(fic,"%f",&fpoub);/*412*/
    fscanf(fic,"%f",&rinfMO_moy[0]); /*443*/
    fscanf(fic,"%f",&rinfMO_moy[1]); /*488*/
    fscanf(fic,"%f",&rinfMO_moy[2]); /*531*/
    fscanf(fic,"%f",&rinfMO_moy[3]); /*547*/
    fscanf(fic,"%f",&fpoub); /*620*/
    fscanf(fic,"%f",&rinfMO_moy[4]); /*667*/
    fscanf(fic,"%f",&rinfMO_moy[rinfMO_NE-1]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&rinfMO_moy[rinfMO_NES-1]);
    fclose(fic);
    sprintf(nomfic,"%s/%s", lutpath, rinfMO_LUT_ECART);
    if( (fic=fopen(nomfic,"r")) == NULL) {perror(nomfic); exit(-1);}
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&rinfMO_ecart[0]); /*443*/
    fscanf(fic,"%f",&rinfMO_ecart[1]); /*488*/
    fscanf(fic,"%f",&rinfMO_ecart[2]); /*531*/
    fscanf(fic,"%f",&rinfMO_ecart[3]); /*547*/
    fscanf(fic,"%f",&fpoub); /*620*/
    fscanf(fic,"%f",&rinfMO_ecart[4]); /*667*/
    fscanf(fic,"%f",&rinfMO_ecart[rinfMO_NE-1]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&rinfMO_ecart[rinfMO_NES-1]);
    fclose(fic);
  }    
}

/* ===============================================  SEAWIFS  =============================================== */
/* -------------------------------------------------------------------------------
 Calcul du Kd a partir des poids
 - Input:
  input[NE] = Rrs 443 490 510 555 Lambda
 ------------------------------------------------------------------------------- */
float rsupSW_neuron_passe_avant(float input[rsupSW_NE+1])
{
  float a[rsupSW_NC1], b[rsupSW_NC2], y=0.0, x[rsupSW_NE];
  int i,j;

  /* Normalisation */
  for(i=0; i<rsupSW_NE; i++){
    x[i] = ((2./3.)*(input[i]-rsupSW_moy[i]))/rsupSW_ecart[i];
  }

  for(i=0;i<rsupSW_NC1;i++){
    a[i] = 0.0;
    for(j=0;j<rsupSW_NE;j++){
      a[i] += (x[j]*rsupSW_w1[j][i]);
    }
    a[i] = 1.715905*(float)tanh((2./3.)*(double)(a[i] + rsupSW_b1[i]));
  }
  for(i=0;i<rsupSW_NC2;i++){
    b[i] = 0.0;
    for(j=0;j<rsupSW_NC1;j++){
      b[i] += (a[j]*rsupSW_w2[j][i]);
    }
    b[i] = 1.715905*(float)tanh((2./3.)*(double)(b[i] + rsupSW_b2[i]));
  }
  for(j=0;j<rsupSW_NC2;j++){
    y += (b[j]*rsupSW_w3[j]);
  }

  /* Denormalisation */
  y = 1.5*(y + rsupSW_b3)*rsupSW_ecart[rsupSW_NES-1] + rsupSW_moy[rsupSW_NES-1];
  y = (float)pow(10.,y);
  return(y);
}

/* -------------------------------------------------------------------------------
 Calcul du Kd a partir des poids
 - Input:
  input[NE] = Rrs 443 490 510 555 670 Lambda
 ------------------------------------------------------------------------------- */
float rinfSW_neuron_passe_avant(float input[rinfSW_NE])
{
  float a[rinfSW_NC1], b[rinfSW_NC2], y=0.0, x[rinfSW_NE];
  int i,j;

  /* Normalisation */
  for(i=0; i<rinfSW_NE; i++){
    x[i] = ((2./3.)*(input[i]-rinfSW_moy[i]))/rinfSW_ecart[i];
  }

  for(i=0;i<rinfSW_NC1;i++){
    a[i] = 0.0;
    for(j=0;j<rinfSW_NE;j++){
      a[i] += (x[j]*rinfSW_w1[j][i]);
    }
    a[i] = 1.715905*(float)tanh((2./3.)*(double)(a[i] + rinfSW_b1[i]));
  }
  for(i=0;i<rinfSW_NC2;i++){
    b[i] = 0.0;
    for(j=0;j<rinfSW_NC1;j++){
      b[i] += (a[j]*rinfSW_w2[j][i]);
    }
    b[i] = 1.715905*(float)tanh((2./3.)*(double)(b[i] + rinfSW_b2[i]));
  }
  for(j=0;j<rinfSW_NC2;j++){
    y += (b[j]*rinfSW_w3[j]);
  }

  /* Denormalisation */
  y = 1.5*(y + rinfSW_b3)*rinfSW_ecart[rinfSW_NES-1] + rinfSW_moy[rinfSW_NES-1];
  y = (float)pow(10.,y);
  return(y);
}


/* ===============================================  MODIS  ================================================ */
/* -------------------------------------------------------------------------------
 Calcul du Kd a partir des poids
 - Input:
  input[NE] = Rrs 443 490 510 555 Lambda
 ------------------------------------------------------------------------------- */
float rsupMO_neuron_passe_avant(float input[rsupMO_NE+1])
{
  float a[rsupMO_NC1], b[rsupMO_NC2], y=0.0, x[rsupMO_NE];
  int i,j;

  /* Normalisation */
  for(i=0; i<rsupMO_NE; i++){
    x[i] = ((2./3.)*(input[i]-rsupMO_moy[i]))/rsupMO_ecart[i];
  }

  for(i=0;i<rsupMO_NC1;i++){
    a[i] = 0.0;
    for(j=0;j<rsupMO_NE;j++){
      a[i] += (x[j]*rsupMO_w1[j][i]);
    }
    a[i] = 1.715905*(float)tanh((2./3.)*(double)(a[i] + rsupMO_b1[i]));
  }
  for(i=0;i<rsupMO_NC2;i++){
    b[i] = 0.0;
    for(j=0;j<rsupMO_NC1;j++){
      b[i] += (a[j]*rsupMO_w2[j][i]);
    }
    b[i] = 1.715905*(float)tanh((2./3.)*(double)(b[i] + rsupMO_b2[i]));
  }
  for(j=0;j<rsupMO_NC2;j++){
    y += (b[j]*rsupMO_w3[j]);
  }

  /* Denormalisation */
  y = 1.5*(y + rsupMO_b3)*rsupMO_ecart[rsupMO_NES-1] + rsupMO_moy[rsupMO_NES-1];
  y = (float)pow(10.,y);
  return(y);
}

/* -------------------------------------------------------------------------------
 Calcul du Kd a partir des poids
 - Input:
  input[NE] = Rrs 443 490 510 555 670 Lambda
 ------------------------------------------------------------------------------- */
float rinfMO_neuron_passe_avant(float input[rinfMO_NE])
{
  float a[rinfMO_NC1], b[rinfMO_NC2], y=0.0, x[rinfMO_NE];
  int i,j;

  /* Normalisation */
  for(i=0; i<rinfMO_NE; i++){
    x[i] = ((2./3.)*(input[i]-rinfMO_moy[i]))/rinfMO_ecart[i];
  }

  for(i=0;i<rinfMO_NC1;i++){
    a[i] = 0.0;
    for(j=0;j<rinfMO_NE;j++){
      a[i] += (x[j]*rinfMO_w1[j][i]);
    }
    a[i] = 1.715905*(float)tanh((2./3.)*(double)(a[i] + rinfMO_b1[i]));
  }
  for(i=0;i<rinfMO_NC2;i++){
    b[i] = 0.0;
    for(j=0;j<rinfMO_NC1;j++){
      b[i] += (a[j]*rinfMO_w2[j][i]);
    }
    b[i] = 1.715905*(float)tanh((2./3.)*(double)(b[i] + rinfMO_b2[i]));
  }
  for(j=0;j<rinfMO_NC2;j++){
    y += (b[j]*rinfMO_w3[j]);
  }

  /* Denormalisation */
  y = 1.5*(y + rinfMO_b3)*rinfMO_ecart[rinfMO_NES-1] + rinfMO_moy[rinfMO_NES-1];
  y = (float)pow(10.,y);
  return(y);
}


/*int main (int argc, char *argv[])
{
  FILE *fic;
  float data_in[rinfSW_NE],result_in,result_out,res_norm, sumInf_rms=0., sumInf_rel=0., sumSup_rms=0., sumSup_rel=0.;
  int i,nb,lu;

  if( (fic=fopen("../LUTS/base_KdSeaWiFS_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_670_log_Kd_lambda_merge_seuil_15_publi_091213_test.dat","r")) == NULL) {perror("input"); exit(-1);}
  neuron_lect_LUTs();

  nb = 0;
  while((lu=fscanf(fic,"%f",&data_in[0])) == 1){
    for(i=1; i<rinfSW_NE; i++)
      fscanf(fic,"%f",&data_in[i]);
    fscanf(fic,"%f",&result_in);
    result_out = rinfSW_neuron_passe_avant(data_in);
    result_in = (float)pow(10.,result_in);
    sumInf_rms += (float)pow((double)(result_in-result_out),2.);
    sumInf_rel +=  (float)sqrt((double)(((result_in-result_out)/result_in) * ((result_in-result_out)/result_in)));
    printf("%f %f",result_in,result_out);

    data_in[5] = data_in[6];
    result_out = rsupSW_neuron_passe_avant(data_in);
    sumSup_rms += (float)pow((double)(result_in-result_out),2.);
    sumSup_rel +=  (float)sqrt((double)(((result_in-result_out)/result_in) * ((result_in-result_out)/result_in)));
    printf(" %f\n",result_out);
    nb++;
  }
  printf("rmseInf = %f, rel errInf= %f,  rmseSup = %f, rel errSup= %f  nb = %d\n",(float)sqrt((double)sumInf_rms/nb),sumInf_rel/nb, (float)sqrt((double)sumSup_rms/nb),sumSup_rel/nb, nb);



  fclose(fic);
  exit(1);
}*/

/*int main (int argc, char *argv[])
{
  FILE *fic;
  float fpoub, data_in[rinfSW_NE],result_in,result_out,res_norm, sumInf_rms=0., sumInf_rel=0., sumSup_rms=0., sumSup_rel=0.;
  int i,nb,lu;

  printf(" -- SUP --\n");
  if( (fic=fopen("../LUTS/base_KdMODIS_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_555_log_Kd_lambda_merge_seuil_15_ss_645_publi_test_sup_CSIRO.dat","r")) == NULL) {perror("input"); exit(-1);}
  neuron_lect_LUTswitch(idMODIS);

  nb = 0;
  while((lu=fscanf(fic,"%f",&fpoub)) == 1){
    for(i=0; i<4; i++)
      fscanf(fic,"%f",&data_in[i]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&data_in[4]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&result_in);
    
    result_out = rsupMO_neuron_passe_avant(data_in);
    result_in = (float)pow(10.,result_in);
    if (data_in[4] == 488 ){
    sumSup_rms += (float)pow((double)(result_in-result_out),2.);
    sumSup_rel +=  (float)sqrt((double)(((result_in-result_out)/result_in) * ((result_in-result_out)/result_in)));
    nb++;
    }
  }
  printf("rmseSup = %f, rel errSup= %f  nb = %d\n",(float)sqrt((double)sumSup_rms/nb),sumSup_rel/nb, nb);
  fclose(fic);
  
  printf(" -- INF --\n");
  if( (fic=fopen("../LUTS/base_KdMODIS_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_670_log_Kd_lambda_merge_seuil_15_ss_645_publi_test_inf_CSIRO.dat","r")) == NULL) {perror("input"); exit(-1);}
  nb = 0;
  while((lu=fscanf(fic,"%f",&fpoub)) == 1){
    for(i=0; i<4; i++)
      fscanf(fic,"%f",&data_in[i]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&data_in[4]);
    fscanf(fic,"%f",&data_in[5]);
    fscanf(fic,"%f",&fpoub);
    fscanf(fic,"%f",&result_in);
    
    result_out = rinfMO_neuron_passe_avant(data_in);
    result_in = (float)pow(10.,result_in);
    if (data_in[5] == 488 ){
    sumInf_rms += (float)pow((double)(result_in-result_out),2.);
    sumInf_rel +=  (float)sqrt((double)(((result_in-result_out)/result_in) * ((result_in-result_out)/result_in)));
    nb++;
    }
  }
  printf("rmseInf = %f, rel errInf= %f  nb = %d\n",(float)sqrt((double)sumInf_rms/nb),sumInf_rel/nb, nb);
  
  
  exit(1);
}*/