/*
 * \file Aeolus-SkillScore.cpp
 *
 * $Id$
 *
 * Cost function for parameter optimisation of Aeolus when coupled to FMS/POEM.
 *
 */

static HorizontalVariable** ua_agg, **va_agg;
static MeridionalVariable1D **u_za_agg;
static int *ntsteps;

static void
init_Aeolus_Score(//, const int ntsteps ///< select between monthly, seasonal, annual evaluation
                  ) {
  const char *monthabbr[12] = {
    "jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec"
  };

  // allocate arrays for temporal aggregation
  // 12 monthly means
  ua_agg = new HorizontalVariable *[12];
  va_agg = new HorizontalVariable *[12];
  u_za_agg = new MeridionalVariable1D *[12];
  ntsteps = new int[12];

  for (int i = 0; i < 12; i++) {
    ua_agg[i] = new HorizontalVariable(*grid, EBL, ((string)("ua_monthmean_")+monthabbr[i]).c_str(), "m/s",
                                       "",
                                       Variable::INTERP_OUT_REPLICATE|Variable::INTERP_IN_AVERAGE,
                                       0., -1. );
    va_agg[i] = new HorizontalVariable(*grid, EBL, ((string)("va_monthmean_")+monthabbr[i]).c_str(), "m/s",
                                       "",
                                       Variable::INTERP_OUT_REPLICATE|Variable::INTERP_IN_AVERAGE,
                                       0., -1. );
    u_za_agg[i] = new MeridionalVariable1D(*grid, EBL, ((string)("u_za_monthmean_")+monthabbr[i]).c_str(), "m/s",
                                           "",
                                       0., -1. );
    ua_agg[i]->UniformValue(0.0);
    va_agg[i]->UniformValue(0.0);
    u_za_agg[i]->UniformValue(0.0);
    ntsteps[i] = 0;
  }
#ifdef DEBUG_SKILLSCORE
  // make time units and calendar compatible with FMS diag manager output
  scorecdf = new NetCDFOutput(*grid, "debug_skillscore", false, "$Id$", "days since 0001-01-01 00:00:00", "NOLEAP");
  for (int i = 0; i < 12; i++) {
    scorecdf->AddVariable(*(ua_agg[i]));
    scorecdf->AddVariable(*(va_agg[i]));
    scorecdf->AddVariable(*(u_za_agg[i]));
  }
#endif
}

/// Aggregate to monthly mean values.
/// Start at 1-Jan of the 2nd year of simulation.
/// With fixed-surface-condition runs, we can assume that already in the 2nd year we have
/// reached a steady state.
static void
Aeolus_Score_Add_Tstep(const int year, const int month) {
  static int prevyear = -1;
  static int agg_year = 0;
  cout << __FUNCTION__ << " year " << year << " month " << month << endl;
  if (-1 != prevyear && year != prevyear) agg_year++; // we are at least in the 2nd year of simulation
  cout << " prevyear " << prevyear << " aggyear " << agg_year;
  if (1 == agg_year) {
    static HorizontalVariable tmp(*grid, EBL, "tmp", "", "");
    static MeridionalVariable1D tmp1d(*grid, EBL, "tmp1d", "", "");
    cout << __FUNCTION__ << " aggregating" << endl;
    tmp.ImportFrom(standing_wave->ua_());
    *(ua_agg[month-1]) += tmp;
    tmp.ImportFrom(standing_wave->va_());
    *(va_agg[month-1]) += tmp;
    tmp1d.ImportFrom(wind->u_za_());
    *(u_za_agg[month-1]) += tmp1d;
    ntsteps[month-1]++;
#ifdef DEBUG_SKILLSCORE
    scorecdf->OutputToFile(year*365+month*30); // approximate date
#endif
  }
  prevyear = year;
}

// from Sonja: 
///p/projects/climber3/molnos/Taylor_step1_20_05_2016_el_nino_la_nina_new/main_Aeolus_DynamicalCore_TaylorSkillScore.cpp::taylor_skillscore*()
static double
taylor_skillscore(const HorizontalVariable& xm, ///< model data
                  const HorizontalVariable& xo, ///< observation data
                  const HorizontalVariable& weights ///< weight function
                  ) {
  static HorizontalVariable tmp(*grid, EBL, "tmp", "", "");

  // weighted Mean values
  double weights_sum = weights.sum();
  tmp = weights;
  tmp *= xm;
  double mean_model=tmp.Mean() / weights_sum;
  tmp = weights;
  tmp *= xo;
  double mean_obs = tmp.Mean() / weights_sum;

  //Spatial statistical moments
  double std_model = 0., std_obs = 0., crossCorrelation = 0.;
  for (vector<Cell>::const_iterator it=grid->EBLCellsBegin(); it < grid->EBLCellsEnd(); it++) {
    unsigned int id = it->ID();
    std_model+=weights.Read(id)*(xm.Read(id)-mean_model)*(xm.Read(id)-mean_model);
    std_obs+=weights.Read(id)*(xo.Read(id)-mean_obs)*(xo.Read(id)-mean_obs);
    crossCorrelation+=weights.Read(id)*(xm.Read(id)-mean_model)*(xo.Read(id)-mean_obs);
  }
  double sigma_model=sqrt(std_model/weights_sum);
  double sigma_obs=sqrt(std_obs/weights_sum);

  //ratio between standard deviation of model and of obs -> 1 ==perfect
  //                                                       <1 == std of model smaller than std of obs
  //                                                       >1 == std of model greater than std of obs
  double amplitude=sigma_model/sigma_obs;
  //ratio of the spatial correlation -> mean between obs and model
  double ratio_spatialCorrelation=crossCorrelation/(weights_sum*sigma_model*sigma_obs);

  //finally calculate Taylor Skillscore
  return 0.25*pow(1.0+ratio_spatialCorrelation,4)/((amplitude+1.0/amplitude)*(amplitude+1.0/amplitude));
}

static double
taylor_skillscore(const MeridionalVariable1D& xm, ///< model data
                  const MeridionalVariable1D& xo, ///< observation data
                  const MeridionalVariable1D& weights ///< weight function
                  ) {
  static MeridionalVariable1D tmp(*grid, EBL, "tmp", "", "");

  // weighted Mean values
  double weights_sum = weights.sum();
  tmp = weights;
  tmp *= xm;
  double mean_model=tmp.Mean() / weights_sum;
  tmp = weights;
  tmp *= xo;
  double mean_obs = tmp.Mean() / weights_sum;

  //Spatial statistical moments
  double std_model = 0., std_obs = 0., crossCorrelation = 0.;
  for (unsigned int j = 0; j < grid->MeridionalCells(); j++) {
    std_model+=weights.Read(j)*(xm.Read(j)-mean_model)*(xm.Read(j)-mean_model);
    std_obs+=weights.Read(j)*(xo.Read(j)-mean_obs)*(xo.Read(j)-mean_obs);
    crossCorrelation+=weights.Read(j)*(xm.Read(j)-mean_model)*(xo.Read(j)-mean_obs);
  }
  double sigma_model=sqrt(std_model/weights_sum);
  double sigma_obs=sqrt(std_obs/weights_sum);

  //ratio between standard deviation of model and of obs -> 1 ==perfect
  //                                                       <1 == std of model smaller than std of obs
  //                                                       >1 == std of model greater than std of obs
  double amplitude=sigma_model/sigma_obs;
  //ratio of the spatial correlation -> mean between obs and model
  double ratio_spatialCorrelation=crossCorrelation/(weights_sum*sigma_model*sigma_obs);

  //finally calculate Taylor Skillscore
  return 0.25*pow(1.0+ratio_spatialCorrelation,4)/((amplitude+1.0/amplitude)*(amplitude+1.0/amplitude));
}


static double
Aeolus_Score_Calc() {
  /* first finish aggregation */
  for (int i = 0; i < 12; i++) {
    *(ua_agg[i]) /= ntsteps[i];
    *(va_agg[i]) /= ntsteps[i];
    *(u_za_agg[i]) /= ntsteps[i];
  }

  // Init weights

  // Needed for simple skillscore calculation --> weights for error calculation between model and observation (empirical formula)
  // Cells are weighted by area, normalized to [0 .. 1]
  MeridionalVariable1D weights1d(*grid, EBL, "weight1d"," ", "area weights for simple skillscore",  0.0, 1.0);
  
  double wsum=0;
  for (unsigned int lat = 0; lat < grid->MeridionalCells(); lat++) {
    const Cell *cell = grid->GetCell(0,lat,3);
    double latfromindex = cell->lat();
    double w = 0;
    if (latfromindex >-60) w = cell->FirstConnectedZInterface()->Area();
    weights1d.Write(lat,w);
    wsum+=w;
  }
  weights1d /= wsum;



  // weight is cell area, between 60S and 10S / 10N and 70N, else 0
  // dont use equator data because there were problems with wrong sign observed by Sonja
  HorizontalVariable weights2d(*grid, EBL, "weight2d"," ", "area weights for simple skillscore",  0.0, 1.0);

  wsum = 0;
  for (vector<Cell>::const_iterator it = grid->EBLCellsBegin(); it < grid->EBLCellsEnd(); it++){
    unsigned int id = it->ID();
    double lat = it->lat();

    double w = 0;
    if ((lat>-60 && lat<-10) or (lat>10 && lat<70)) {
      w = it->FirstConnectedZInterface()->Area();
    }
    weights2d.Write(id,w);//Level=3
    wsum+=w;
  }
  weights2d /= wsum;

  //
  // Prepare to read observation data (ERA interim reanalysis in this case)
  //
  MeridionalVariable1D u_zonal_mean1D(*grid, EBL, "u_za","m s**-1", "u_za observation data",  -100, 100);

  //MeridionalVariable1D EKE1D(*grid, EBL, "EKE","m**2 s**-2", "", 0.00000001, 500);
  //MeridionalVariable1D EKE1Dmodel(*grid, EBL, "EKE","m**2 s**-2", "", -300, 300);

  HorizontalVariable ua_obs(*grid, EBL, "ua_obs","m s**-1", "ua observation data",Variable::INTERP_OUT_REPLICATE|Variable::INTERP_IN_AVERAGE,  -100, 100);
  HorizontalVariable va_obs(*grid, EBL, "va_obs","m s**-1", "va observation data", Variable::INTERP_OUT_REPLICATE|Variable::INTERP_IN_AVERAGE, -100, 100);



  // Sonja does 3 series of experiments: all months, el nino months, la nina months.
  // Here we do only the all-months experiment
  char* Ekefile = "u_za_v_za_Eke_uava_1983_2009_3_75_3_75_-178_ymonmean.cdf";

  //load u_za and v_za
  NetCDFInput input_u_za_uava(*grid, Ekefile);
  input_u_za_uava.GetVariable(u_zonal_mean1D);
  // input_u_za_v_za_Eke_uava.GetVariable(v_zonal_mean1D);
  //load ua and va (pressure 500mb)
  input_u_za_uava.GetVariable(ua_obs);
  input_u_za_uava.GetVariable(va_obs);


  valarray<double> skillscore_vector_u_za(12), skillscore_vector_ua(12), skillscore_vector_va(12);

  // skill score is computed for each timestep, i.e. monthly mean, separately.
  for (int month = 0; month < 12; month++) {
    input_u_za_uava.InputFromFile(month);


    //Planetary waves -> ua va
    skillscore_vector_ua[month] = taylor_skillscore(*(ua_agg[month]), ua_obs, weights2d);
    skillscore_vector_va[month] = taylor_skillscore(*(va_agg[month]), va_obs, weights2d);

    //LargeScaleWind -> u_za & v_za
    skillscore_vector_u_za[month] = taylor_skillscore(*(u_za_agg[month]), u_zonal_mean1D, weights1d);
  }

  //skillscore is mean value of vector skillscore_vector values[t]

  //LargeScaleWind
  double skillscore_u_za = skillscore_vector_u_za.sum()/skillscore_vector_u_za.size();

  //Planetary waves
  double skillscore_ua = skillscore_vector_ua.sum()/skillscore_vector_ua.size();
  double skillscore_va = skillscore_vector_va.sum()/skillscore_vector_va.size();

  cout<<"skillscore_u_za; "<<skillscore_u_za<<endl;
  cout<<"skillscore_ua "<<skillscore_ua<<endl;
  cout<<"skillscore_va "<<skillscore_va<<endl;

  double skillscore = skillscore_u_za*skillscore_ua*skillscore_va;
  cout<<"skillscore "<<skillscore<<endl;
  if (skillscore>100000) skillscore = 100000;

  return skillscore;
}

/// Clean up.
/// \param days is number of days since begin of simulation
static void
Aeolus_Score_End(double days) {
#ifdef DEBUG_SKILLSCORE
    scorecdf->OutputToFile(days); // approximate date
    delete scorecdf;
#endif

  // deallocate storage for temporal aggregation
  for (int i = 0; i < 12; i++) {
    delete ua_agg[i];
    delete va_agg[i];
    delete u_za_agg[i];
  }

  delete[] ua_agg;
  delete[] va_agg;
  delete[] u_za_agg;
  delete[] ntsteps;
}