/ensemble/ClassEnsembleForecast.cxx

// Copyright (C) 2013-2015 INERIS
// Author(s) : Edouard Debry
// 
// This file is part of GPEC.
// 
// GPEC is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// GPEC is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License
// along with GPEC.  If not, see <http://www.gnu.org/licenses/>.

#ifndef GPEC_FILE_CLASS_ENSEMBLE_FORECAST_CXX

#include "ClassEnsembleForecast.hxx"

namespace GPEC
{
  // Simulation index reverse at computation time.
  template<class C>
  inline void ClassEnsembleForecast<C>::compute_container_index_reverse()
  {
    int Nt = target_.GetNt();
    container_index_reverse_.Reallocate(Nt);

    for (int t = this->t_begin_; t < this->t_end_; t++)
      {
        // We could allocate only up to t_end, but this would break
        // ComputeWeight method for BCLC which runs on all time steps,
        // not only those before t.
        container_index_reverse_(t).Reallocate(Nt);

        if (this->container_index_(t).GetSize() == 0)
          continue;

        container_index_reverse_(t)(t) = this->container_index_(t);

        for (int h = 0; h < Nt; h++)
          {
            container_index_reverse_(t)(h).Reallocate(this->container_index_(t).GetSize());
            if (this->container_index_(h).GetSize() == this->Ncontainer_)
              container_index_reverse_(t)(h).Fill();
            else
              {
                int i(0), j(0);
                while (i < this->container_index_(h).GetSize() && j < this->container_index_(t).GetSize())
                  if (this->container_index_(h)(i) > this->container_index_(t)(j))
                    break;
                  else if (this->container_index_(h)(i) < this->container_index_(t)(j))
                    i++;
                  else
                    container_index_reverse_(t)(h)(j++) = i++;

                container_index_reverse_(t)(h).Resize(j);
              }
          } 
      }
  }


  // Load target.
  template<class C>
  inline void ClassEnsembleForecast<C>::load_target(const int t_begin, const int t_end,
                                                    const Vector1I &location_index,
                                                    Vector2I &target_index,
                                                    Vector2T &target_data) const
  {
#ifdef GPEC_WITH_TIME
    clock_t c = clock();
#endif

    int Nl = location_index.GetSize();
    int Nt = target_.GetNt();
    target_index.Reallocate(Nt);
    target_data.Reallocate(Nt);

    int n = target_time_step_.GetSize();
    for (int i = 0; i < n; i++)
      {
        int t = target_time_step_(i);
        target_index(t) = location_index;
        target_data(t).Reallocate(Nl);
      }

    if (this->in_memory_)
      {
        int location_offset = this->location_offset_rank_(this->rank_);

        for (int i = 0; i < n; i++)
          {
            int t = target_time_step_(i);
            if (target_index_(t).GetSize() > 0)
              {
                int i(0), j(0), l(0);
                while (l < target_index_(t).GetSize() && i < Nl)
                  if (target_index_(t)(l) > location_index(i))
                    i++;
                  else if (target_index_(t)(l) < location_index(i))
                    l++;
                  else
                    {
                      target_index(t)(j) = location_index(i);
                      target_data(t)(j++) = target_data_(t)(l);
                      i++;
                      l++;
                    }
                target_index(t).Resize(j);
                target_data(t).Resize(j);
              }
            else
              for (int l = 0; l < Nl; l++)
                target_data(t)(l) = target_data_(t)(location_index(l) - location_offset);
          }
      }
    else
      if (C::Type() == CONTAINER_SIM_TYPE)
        {
          Vector<real_io_container> data(n);

          ifstream fin(target_file_.c_str(), ifstream::binary);

          for (int l = 0; l < Nl; l++)
            {
              fin.seekg(n * location_index(l) * sizeof(real_io_container), ios_base::beg);
              fin.read(reinterpret_cast<char*>(data.GetData()), n * sizeof(real_io_container));

              for (int i = 0; i < n; i++)
                target_data(target_time_step_(i))(l) = real_container(data(i));
            }

          fin.close();
        }
      else
        {
          Vector1I loc_index;
          Vector1T data;

          int Nl;
          int Nl_container = target_.Nl_;

          ifstream fin(target_file_.c_str(), ifstream::binary);
          for (int t = 0; t < Nt; t++)
            {
              fin.seekg(sizeof(int), ios_base::cur);
              fin.read(reinterpret_cast<char*>(&Nl), sizeof(int));

              if (Nl > 0)
                {
                  fin.seekg(sizeof(int), ios_base::cur);

                  if (Nl < Nl_container)
                    {
                      loc_index.Reallocate(Nl);
                      fin.read(reinterpret_cast<char*>(loc_index.GetData()), Nl * sizeof(int));
                    }

                  data.Reallocate(Nl);
                  fin.read(reinterpret_cast<char*>(data.GetData()), Nl * sizeof(real_io_container));

                  if (Nl < Nl_container)
                    {
                      int l(0), loc(0);
                      while (l < Nl && loc < location_index.GetSize())
                        if (loc_index(l) > location_index(loc))
                          loc++;
                        else if (loc_index(l) < location_index(loc))
                          l++;
                        else
                          {
                            target_index(t)(loc) = loc_index(l);
                            target_data(t)(loc) = real_container(data(l));
                            l++;
                            loc++;
                          }

                      target_index(t).Resize(loc);
                      target_data(t).Resize(loc);
                    }
                  else
                    for (int l = 0; l < location_index.GetSize(); l++)
                      target_data(t)(l) = data(location_index(l));
                }
              else
                {
                  fin.seekg(sizeof(int), ios_base::cur);
                  fin.read(reinterpret_cast<char*>(&Nl), sizeof(int));
                  if (Nl > 0)
                    {
                      fin.seekg(sizeof(int), ios_base::cur);
                      if (Nl < Nl_container)
                        fin.seekg(Nl * sizeof(int), ios_base::cur);
                      fin.seekg(Nl * sizeof(real_io_container), ios_base::cur);
                    }
                }

              fin.close();
            }
        }

    for (int t = 0; t < Nt; t++)
      if (t < t_begin || t >= t_end)
        {
          target_index(t).Clear();
          target_data(t).Clear();
        }

#ifdef GPEC_WITH_TIME
    load_target_time += clock() - c;
#endif
  }


  // Constructors.
  template<class C>
  ClassEnsembleForecast<C>::ClassEnsembleForecast(const string &name)
    : ClassEnsembleBase<C>(name), target_file_("")
  {
#ifdef GPEC_WITH_LOGGER
    *(ClassGPEC::log_ptr_) << Bmagenta() << Fblack() << Info() << Reset()
                           << "Instantiate forecast ensemble \"" << name
                           << "\" (you will need one target : SetTarget() ...)." << endl;
#endif

#ifdef GPEC_WITH_TIME
    load_target_time = 0;
    compute_weight_time = 0;
#endif

    return;
  }


  // Init.
  template<class C>
  void ClassEnsembleForecast<C>::Init(const string &species, int date_reference_int, const bool operational)
  {
#ifdef GPEC_WITH_LOGGER
    *(ClassGPEC::log_ptr_) << Byellow() << Fblack() << Info() << Reset()
                           << "Master initialization for forecast ensemble \"" << this->name_
                           << "\" with species \"" << species << "\"." << Reset().Str() << endl;
    *(ClassGPEC::log_ptr_) << Fyellow() << Info(2) << Reset()
                           << "Operational ? " << (operational ? "YES" : "NO") << endl;
#endif

    // Read configuration.
    Ops::Ops ops(this->configuration_file_);

    string prefix = "ensemble." + this->name_ + ".init.";
    ops.SetPrefix(prefix);

    string frequency = ops.Get<string>("frequency", "", "hourly");
    string target_name = ops.Get<string>("target.name");
    string target_type = ops.Get<string>("target.type", "", "");
    int date_begin_int = ops.Get<int>("date.begin");
    int date_end_int = ops.Get<int>("date.end_excluded");

    vector<string> model_list = ops.GetEntryList("model");

    vector<string> model_name;
    vector<int> forecast_day;
    vector<string> method;

    for (int s = 0; s < int(model_list.size()); s++)
      {
        ops.SetPrefix(prefix + "model." + model_list[s] + ".");
        model_name.push_back(ops.Get<string>("name", "", model_list[s]));
        forecast_day.push_back(ops.Get<int>("forecast_day", "", 0));
        method.push_back(ops.Get<string>("method", "", "bilinear"));
      }

    ops.Close();

    if (date_reference_int == 0)
      date_reference_int = convert<int>(get_current_time("%Y%m%d"));

#ifdef GPEC_WITH_LOGGER
    *(ClassGPEC::log_ptr_) << Fyellow() << Info(2) << Reset()
                           << "Reference date is " << date_reference_int << endl;
#endif

    Date date_begin;
    try
      {
        date_begin = get_date_from_int(date_begin_int);
      }
    catch (...)
      {
        date_begin.SetDate(date_reference_int);
        date_begin.AddDays(date_begin_int);
        date_begin_int = get_int_from_date(date_begin);
      }

    Date date_end;
    try
      {
        date_end = get_date_from_int(date_end_int);
      }
    catch (...)
      {
        date_end.SetDate(date_reference_int);
        date_end.AddDays(date_end_int);
        date_end_int = get_int_from_date(date_end);
      }

#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Fcyan() << Info(1) << Reset()
                      << "Ensemble beginning date is \"" << date_begin_int
                      << "\" and ending (excluded) date is \"" << date_end_int << "\"." << endl;
#endif

    if (target_type.empty())
      target_ = C(target_name, species);
    else
      target_ = C(target_type, target_name, species);

    target_.LoadE(date_begin_int, date_end_int);

    if (frequency == "peak")
      target_.template FilterTime<ClassFilterTimePeak<real_container> >();
    else if (frequency == "mean")
      target_.template FilterTime<ClassFilterTimeMean<real_container> >();

    ostringstream sout;
    sout <<  this->container_directory_ << "/" << species
         << "_D" << showpos << int(rint(target_.GetForecastDayMean()))
         << "_" << target_name << "_" << target_type
         << C::GetFileExtension();

    target_file_ = sout.str();

    target_.Export(target_file_);

#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Bcyan() << Fblack() << Info() << Reset()
                      << "Target of ensemble \"" << this->name_
                      << "\" is \"" << target_.GetName() << "/" << target_.GetType() << "\"." << endl;
#endif

    for (int s = 0; s < int(model_name.size()); s++)
      {
        sim_type sim(model_name[s], species);
        // True means we will replace missing days as much as possible.
        sim.Load(date_begin_int, date_end_int, forecast_day[s], true, operational);

        if (frequency == "peak")
          sim.template FilterTime<ClassFilterTimePeak<real_container> >();
        else if (frequency == "mean")
          sim.template FilterTime<ClassFilterTimeMean<real_container> >();

        if (method[s] == "bilinear")
          AddSimulation<ClassMethodBilinear<real_container> >(sim);
        else if (method[s] == "max")
          AddSimulation<ClassMethodMax<real_container> >(sim);
        else if (method[s] == "min")
          AddSimulation<ClassMethodMin<real_container> >(sim);
        else if (method[s] == "upleft")
          AddSimulation<ClassMethodUpLeft<real_container> >(sim);
        else if (method[s] == "upright")
          AddSimulation<ClassMethodUpRight<real_container> >(sim);
        else if (method[s] == "downleft")
          AddSimulation<ClassMethodDownLeft<real_container> >(sim);
        else if (method[s] == "downright")
          AddSimulation<ClassMethodDownRight<real_container> >(sim);
        else if (method[s] == "nearest")
          AddSimulation<ClassMethodNearest<real_container> >(sim);
        else
          throw GPEC::Error("Unknown interpolation method \"" + method[s] + "\".");
      }
  }


  // Destructor.
  template<class C>
  ClassEnsembleForecast<C>::~ClassEnsembleForecast()
  {
    return;
  }


  // Add Target.
  template<class C>
  void ClassEnsembleForecast<C>::SetTarget(const string species,
                                           string target_file,
                                           string config_file)
  {
    // Read configuration or argument, if nothing in configuration.
    Ops::Ops ops(this->configuration_file_);
    target_file_ = ops.Get<string>("ensemble." + this->name_ + ".target", "", target_file);
    ops.Close();

    // If empty, then no way to go on.
    if (target_file_.empty())
      throw GPEC::Error("Empty target file.");

    // Put species name in generic target file name, if necessary.
    target_file_ = find_replace(target_file_, "&s", species);

    // Construct target.
    if (target_file_[0] != '/')
      target_file_ = this->container_directory_ + "/" + target_file_;

    check_file(target_file_);
    target_.LoadFile(target_file_, config_file);
  }


  template<class C>
  void ClassEnsembleForecast<C>::SetTarget(const C &target)
  {
    ostringstream sout;
    sout <<  this->container_directory_ << "/" << target.GetSpecies()
         << "_D" << showpos << int(rint(target.GetForecastDayMean()))
         << "_" << target.GetName() << "_" << target.GetType()
         << C::GetFileExtension();

    target_file_ = sout.str();

    target.Export(target_file_);

    target_.FullCopy(target);
  }


  // Add Simulations.
  template<class C>
  void ClassEnsembleForecast<C>::AddSimulation(Vector<int> &index)
  {
    ClassPolyContainer<C>::AddContainer(index, target_.species_);
  }


  template<class C> template<class M>
  void ClassEnsembleForecast<C>::AddSimulation(const ClassContainerSimulation<real, real_io> &simulation, string container_file)
  {
    M algorithm;

    if (container_file.empty())
      {
        ostringstream sout;
        sout <<  this->container_directory_ << "/" << simulation.GetSpecies()
             << "_D" << showpos << int(rint(simulation.GetForecastDayMean()))
             << "_" << simulation.GetName() << "_" << target_.type_
             << "_" << algorithm.Name() << C::GetFileExtension();
        //   << "_" << algorithm.Name() << "_" << noshowpos << this->Ncontainer_  << C::GetFileExtension();
        container_file = sout.str();
      }

    ClassPolyContainer<C>::template AddSimulation<M>(simulation, target_.type_, container_file);
  }


  template<class C>
  void ClassEnsembleForecast<C>::AddSimulationConstant(const real constant_value, string container_file)
  {
    if (this->Ncontainer_ == 0)
      throw GPEC::Error("You need to provide at least one container before adding a constant one.");

    C container;
    container.LoadFile(this->container_file_[0]);

    for (int t = 0; t < container.Nt_; t++)
      {
        container.forecast_(t).Reallocate(1);
        container.forecast_(t)(0) = 0;
        container.data_(t).Reallocate(container.Nl_);
        container.data_(t).Fill(real_container(constant_value));
        container.location_index_(t).Clear();
      }

    container.name_ = "Constant" + to_str(constant_value);

    if (container_file.empty())
      {
        ostringstream sout;
        sout <<  this->container_directory_ << "/" << container.species_
             << "_D" << showpos << int(rint(container.GetForecastDayMean()))
             << "_" << container.name_ << "_" << container.type_
             << C::GetFileExtension();
        //<< "_" << noshowpos << this->Ncontainer_  << C::GetFileExtension();
        container_file = sout.str();
      }

    ClassPolyContainer<C>::AddContainer(container, container_file);
  }


  template<class C>
  void ClassEnsembleForecast<C>::AddSimulationRandom(const real value_min, const real value_max,
                                                     string container_file, int seed)
  {
    if (this->Ncontainer_ == 0)
      throw GPEC::Error("You need to provide at least one container before adding a constant one.");

    srand48(seed);

    C container;
    container.LoadFile(this->container_file_[0]);

    real_container value_delta = value_max - value_min;

    for (int t = 0; t < container.Nt_; t++)
      {
        container.forecast_(t).Reallocate(1);
        container.forecast_(t)(0) = 0;
        container.data_(t).Reallocate(container.Nl_);
        for (int l = 0; l < container.Nl_; l++)
          container.data_(t)(l) = value_min + drand48() * value_delta;
        container.location_index_(t).Clear();
      }

    container.name_ = "Random" + to_str(seed);

    if (container_file.empty())
      {
        ostringstream sout;
        sout <<  this->container_directory_ << "/" << container.species_
             << "_D" << showpos << int(rint(container.GetForecastDayMean()))
             << "_" << container.name_ << "_" << container.type_
             << C::GetFileExtension();
        //<< "_" << noshowpos << this->Ncontainer_  << C::GetFileExtension();
        container_file = sout.str();
      }

    ClassPolyContainer<C>::AddContainer(container, container_file);
  }


  // Load ensemble.
  template<class C>
  void ClassEnsembleForecast<C>::Load(string config_file)
  {
    int Nt = target_.Nt_;
    if (Nt == 0 || target_file_ == "")
      throw GPEC::Error("Empty target in ensemble forecast \"" + this->name_ + "\", use SetTarget().");

    ClassPolyContainer<C>::Load(config_file);

    target_time_step_.Reallocate(Nt);
    int i(0);
    for (int t = 0; t < Nt; t++)
      if (target_.data_(t).GetSize() > 0)
        target_time_step_(i++) = t;
    target_time_step_.Resize(i);

    int Nl = target_.Nl_;
    this->set_location_per_rank(Nl);

    int location_offset(0);
    if (this->split_memory_)
      {
        Nl = this->Nlocation_rank_(this->rank_);
        location_offset = this->location_offset_rank_(this->rank_);
      }

    if (this->in_memory_)
      {
        int n = target_time_step_.GetSize();

        target_index_.Reallocate(Nt);
        target_data_.Reallocate(Nt);

        for (int i = 0; i < n; i++)
          {
            int t = target_time_step_(i);
            target_index_(t).Reallocate(Nl);
            target_data_(t).Reallocate(Nl);
            if (target_.location_index_(t).GetSize() > 0)
              {
                int i(location_offset), j(0), l(0), imax(Nl + location_offset);
                while (l < target_.location_index_(t).GetSize() && i < imax)
                  if (target_.location_index_(t)(l) > i)
                    i++;
                  else if (target_.location_index_(t)(l) < i)
                    l++;
                  else
                    {
                      target_index_(t)(j) = i;
                      target_data_(t)(j++) = target_.data_(t)(l);
                      i++;
                      l++;
                    }
                target_index_(t).Resize(j);
                target_data_(t).Resize(j);
              }
            else
              for (int j = 0; j < Nl; j++)
                target_data_(t)(j) = target_.data_(t)(j + location_offset);
          }
      }

    target_.PartialClear();
  }


  // Clear ensemble.
  template<class C>
  void ClassEnsembleForecast<C>::Clear()
  {
    target_.Clear();
    target_file_ = "";
    target_time_step_.Clear();
    target_index_.Clear();
    target_data_.Clear();
    container_index_reverse_.Clear();

    ClassEnsembleBase<C>::Clear();
  }


  // Dump configuration.
  template<class C>
  string ClassEnsembleForecast<C>::DumpConfiguration(string config_file) const
  {
    string output = ClassPolyContainer<C>::DumpConfiguration();

    // Write lua configuration.
    Ops::Ops ops;

    int pos = target_file_.find_last_of('/');
    string basedir = target_file_.substr(0, pos);

    string basename = target_file_.substr(pos + 1);

    if (basedir == this->container_directory_)
      ops.Push(this->name_ + ".target", basename);
    else
      ops.Push(this->name_ + ".target", target_file_);

    output += ops.LuaDefinition();
    ops.Close();

    if (! config_file.empty())
      {
        ofstream fout(config_file.c_str());
        fout << output;
        fout.close();
      }

    return output;
  }


  // Get methods.
  template<class C>
  C* ClassEnsembleForecast<C>::GetTargetPtr()
  {
    return &target_;
  }

  template<class C>
  string ClassEnsembleForecast<C>::GetTargetFile() const
  {
    return target_file_;
  }


  // Compute weight.
  template <class C> template<class M>
  void ClassEnsembleForecast<C>::ComputeWeight(const int &t_begin, const int &t_end,
                                               const Vector2I &location_index,
                                               const Vector2T &location_factor)
  {
    if (! this->is_loaded_)
      throw GPEC::Error("Ensemble not loaded yet.");

    // Count algorithm error.
    int count_algorithm_error(0);

    // Clear weights.
    this->ClearWeight();

    this->t_begin_ = t_begin;
    this->t_end_ = t_end;

#ifdef GPEC_WITH_LOGGER
    *(ClassGPEC::log_ptr_) << Byellow() << Fblack() << Info(2) << Reset()
                           << "Compute weights for ensemble \"" << this->name_ << "\" with species \""
                           << target_.species_ << "\" and target \"" << target_.name_ << "\", "
                           << "between time indexes [" << this->t_begin_ << ", " << this->t_end_ << "[." << endl;
#endif

    // Reallocate weights.
    int Nt = target_.GetNt();
    this->weight_.Reallocate(Nt);

    for (int t = this->t_begin_; t < this->t_end_; t++)
      this->weight_(t).Reallocate(location_index.GetSize());

    // Compute reverse index.
    compute_container_index_reverse();

    // Algorithm parameters.
    map<string, real_container> parameter = this->parameter_;
    M algorithm(parameter);

    algorithm.Init(target_, t_begin, t_end);

    this->algorithm_ = algorithm.GetName();

    // Each process rank has one part of weights.
    for (int j = 0; j < location_index.GetSize(); ++j)
      {
#ifdef GPEC_WITH_LOGGER
        *(ClassGPEC::log_ptr_) << Frequency(2) << Bred() << Fwhite() << Info(3) << Reset() << Logger::Date()
                               << "location " << j << "/" << location_index.GetSize() << endl;
#endif

        Vector2I target_index;
        Vector2T target_data;
        load_target(0, t_end, location_index(j), target_index, target_data);

        Vector3T simulation_data;
        this->load_container(location_index(j), target_index, simulation_data);

#ifdef GPEC_WITH_TIME
    clock_t c = clock();
#endif

        Vector2T location_factor2(t_end);
        for (int t = 0; t < t_end; ++t)
          if (target_index(t).GetSize() > 0)
            {
              if (target_index(t).GetSize() == location_index(j).GetSize())
                location_factor2(t) = location_factor(j);
              else
                {
                  location_factor2(t).Reallocate(target_index(t).GetSize());
                  int i(0), k(0);
                  while (i < target_index(t).GetSize())
                    if (location_index(j)(k) == target_index(t)(i))
                      location_factor2(t)(i++) = location_factor(j)(k++);
                    else
                      k++;
                }
            }

        for (int t = this->t_begin_; t < this->t_end_; ++t)
          if (this->container_index_(t).GetSize() > 0)
            if (! algorithm.ComputeWeight(t,
                                          target_,
                                          target_data,
                                          simulation_data,
                                          location_factor2,
                                          container_index_reverse_,
                                          this->weight_(t)(j)))
              count_algorithm_error++;

#ifdef GPEC_WITH_TIME
        compute_weight_time += clock() - c;
#endif
      }

#ifdef GPEC_WITH_TIME
#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Fyellow() << User(0, "TIME") << Reset()
                      << "load_container_time = " << real_container(load_container_time) / CLOCKS_PER_SEC << endl;
    *(this->log_ptr_) << Fyellow() << User(0, "TIME") << Reset()
                      << "load_target_time = " << real_container(load_target_time) / CLOCKS_PER_SEC << endl;
    *(this->log_ptr_) << Fyellow() << User(0, "TIME") << Reset()
                      << "compute_weight_time = " << real_container(compute_weight_time) / CLOCKS_PER_SEC << endl;
#else
    cout << "load_container_time = " << real_container(load_container_time) / CLOCKS_PER_SEC << endl;
    cout << "load_target_time = " << real_container(load_target_time) / CLOCKS_PER_SEC << endl;
    cout << "compute_weight_time = " << real_container(compute_weight_time) / CLOCKS_PER_SEC << endl;
#endif
    load_container_time = clock_t(0);
    load_target_time = clock_t(0);
    compute_weight_time = clock_t(0);
#endif

#ifdef GPEC_WITH_LOGGER
    if (count_algorithm_error > 0)
      {
        *(this->log_ptr_) << Bred() << Logger::Error() << Reset()
                          << "count_algorithm_error = " << count_algorithm_error << endl;
        *(this->log_ptr_) << Fred() << User(0, "EXPLANATION") << Reset()
                          << "Errors coming from the given algorithm are not fatal, they just mean "
                          << "computed weights are not optimal or just inexistant. Most of the time, "
                          << "it is a violated learning period." << endl;
      }
#endif
  }


  template<class C> template<class M>
  void ClassEnsembleForecast<C>::ComputeWeightByTime(const int &t_begin, const int &t_end)
  {
    // Each rank computes part of time steps.
    int Nt_global = t_end - t_begin;

    int Nt_tmp = Nt_global / this->Nrank_;
    int Nt_local = Nt_tmp;
    int Nt_remaining = Nt_global % this->Nrank_;
    if (this->rank_ < Nt_remaining)
      Nt_local++;

    int t_begin_local(t_begin);
    for (int irank = 0; irank < this->rank_; irank++)
      {
        t_begin_local += Nt_tmp;
        if (this->rank_ < Nt_remaining)
          t_begin_local++;
      }

    int t_end_local = t_begin_local + Nt_local;

    // All locations used for all ranks.
    Vector2I location_index(1);
    location_index(0).Reallocate(target_.GetNlocation());
    location_index(0).Fill();

    Vector2T location_factor(1);
    location_factor(0).Reallocate(target_.GetNlocation());
    location_factor(0).Fill(real_container(1));

#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Bgreen() << Info() << Reset() << "Compute weights BY TIME from time step "
                      << t_begin_local << " to " << t_end_local << " (excluded)." << endl;
#endif

    ComputeWeight<M>(t_begin_local, t_end_local, location_index, location_factor);

    // Each rank receives its part of locations.
    this->weight_location_ = location_index;
  }


  template<class C> template<class M>
  void ClassEnsembleForecast<C>::ComputeWeightByLocation(const int &t_begin, const int &t_end)
  {
    // All ranks receives their part of locations to compute.
    this->set_location_per_rank(target_.Nl_);

    Vector2I location_index(this->location_rank_.GetSize());
    Vector2T location_factor(this->location_rank_.GetSize());

    for (int j = 0; j < this->location_rank_.GetSize(); j++)
      {
        location_index(j).PushBack(this->location_rank_(j));
        location_factor(j).PushBack(real_container(1));
      }

#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Bred() << Info() << Reset() << "Compute weights BY LOCATION from location "
                      << this->location_offset_rank_(this->rank_) << " to "
                      << this->Nlocation_rank_(this->rank_) + this->location_offset_rank_(this->rank_) - 1
                      << "." << endl;
#endif

    ComputeWeight<M>(t_begin, t_end, location_index, location_factor);

    this->weight_location_ = location_index;
  }


  template<class C> template<class M>
  void ClassEnsembleForecast<C>::ComputeWeightByCell(const string &type,
                                                     const int &t_begin, const int &t_end)
  {
    real_container distance_cutoff = this->parameter_["distance_cutoff"];
    real_container distance_h = this->parameter_["distance_h"];
    real_container distance_beta = this->parameter_["distance_beta"];
    int location_min = int(this->parameter_["location_min"]);
    int location_max = int(this->parameter_["location_max"]);

    sim_type reference(type, target_.GetSpecies());

    this->set_location_per_rank(reference.GetNlocation());

    Vector2I location_index(this->location_rank_.GetSize());
    Vector2T location_factor(this->location_rank_.GetSize());

    for (int j = 0; j < this->location_rank_.GetSize(); j++)
      {
        int loc = this->location_rank_(j);

        real_container longitude = reference.GetLocationLongitude(loc);
        real_container latitude = reference.GetLocationLatitude(loc);

        int i(0);
        while (i < location_min && distance_cutoff < real_container(30))
          {
            location_index(j).Reallocate(target_.GetNlocation());
            location_factor(j).Reallocate(target_.GetNlocation());
            location_factor.Zero();

            for (int l = 0;  l < target_.GetNlocation(); l++)
              {
                real_container lon = longitude - target_.GetLocationLongitude(l);
                real_container lat = latitude - target_.GetLocationLatitude(l);
                real_container distance = sqrt(lon * lon + lat * lat);

                if (distance < distance_cutoff)
                  {
                    real_container ratio = distance / distance_h;
                    location_index(j)(i) = l;
                    location_factor(j)(i++) = distance_beta * exp(-ratio) * (real_container(1) + ratio) + real_container(1);
                    if (i == location_max)
                      break;
                  }
              }

            distance_cutoff *= real_container(2);
          }

        if (i < location_min)
          throw GPEC::Error("Unable to get enough locations for position " + to_str(loc) + " .");

        location_index(j).Resize(i);
        location_factor(j).Resize(i);
      }

#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Bmagenta() << Fblack() << Info() << Reset()
                      << "Compute weights BY CELL." << endl;
#endif

    ComputeWeight<M>(t_begin, t_end, location_index, location_factor);

    this->weight_location_.Reallocate(this->location_rank_.GetSize());
    for (int j = 0; j < this->location_rank_.GetSize(); j++)
      this->weight_location_(j).PushBack(this->location_rank_(j));
  }


  // Rank diagram.
  template<class C>
  void ClassEnsembleForecast<C>::RankDiagram(Vector1I &rank_diagram) const
  {
    rank_diagram.Reallocate(this->Ncontainer_ + 1);
    rank_diagram.Zero();

    int Nt = target_.GetNt();
    int Nlocation = target_.GetNlocation();

    for (int loc = 0; loc < Nlocation; loc++)
      {
        Vector1I target_location(1);
        target_location(0) = loc;

        Vector2I target_index;
        Vector2T target_data;
        load_target(0, Nt, target_location, target_index, target_data);

        Vector3T simulation_data;
        this->load_container(target_location, target_index, simulation_data);

        for (int t = 0; t < Nt; t++)
          if (target_index(t).GetSize() > 0)
            if (this->container_index_(t).GetSize() == this->Ncontainer_)
              {
                Sort(simulation_data(t)(0));

                if (target_data(t)(0) < simulation_data(t)(0)(0))
                  rank_diagram(0)++;
                else if (target_data(t)(0) > simulation_data(t)(0)(this->Ncontainer_ - 1))
                  rank_diagram(this->Ncontainer_)++;
                else
                  for (int s = 1; s < this->Ncontainer_; s++)
                    if (target_data(t)(0) > simulation_data(t)(0)(s))
                      {
                        rank_diagram(s)++;
                        break;
                      }
              }
      }
  }
}

#define GPEC_FILE_CLASS_ENSEMBLE_FORECAST_CXX
#endif