/container/ClassContainerBase.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_CONTAINER_BASE_CXX

#include "ClassContainerBase.hxx"

namespace GPEC
{
  // Set number of time steps.
  template<class T>
  inline void ClassContainerBase<T>::set_Nt()
  {
    Nt_ = 0;
    if (date_begin_ < date_end_)
      Nt_ = get_time_index(date_end_.GetSecondsFrom(date_begin_), delta_t_);
  }

  // Set number of days.
  template<class T>
  inline void ClassContainerBase<T>::set_Nd()
  {
    Nd_ = 0;
    if (date_begin_ < date_end_)
      Nd_ = date_end_.GetDaysFrom(date_begin_);
    if (Nd_ * Nt_day_ < Nt_)
      Nd_++;
  }


  // Get location index.
  template<class T>
  inline int ClassContainerBase<T>::get_location_index(const int &t, const int &location) const
  {
    if (data_(t).GetSize() == Nl_)
      return location;

    int l;
    for (l = 0; l < data_(t).GetSize(); l++)
      if (location_index_(t)(l) > location)
        return data_(t).GetSize();
      else if (location_index_(t)(l) == location)
        break;
    return l;
  }

  
  // Construct.
  template<class T>
  void ClassContainerBase<T>::construct(const string &prefix_base)
  {
    // Read configuration.
    Ops::Ops ops(this->configuration_file_);
    ops.SetPrefix("general.");
    delta_t_ = ops.Get<T>("Delta_t");

    // Custom values for species.
    value_min_ = ops.Get<T>(species_ + ".value_min", "", T(CONTAINER_VALUE_MIN));
    value_max_ = ops.Get<T>(species_ + ".value_max", "", T(CONTAINER_VALUE_MAX));

    // Set type and name of container.
    ops.SetPrefix(prefix_base);

    if (type_.empty())
      type_ = ops.Get<string>(name_ + ".type", "", name_);

#ifdef GPEC_WITH_LOGGER
    if (! ops.Exists(type_))
      *(this->log_ptr_) << Warning() << "Container type \"" << type_
                        << "\" not found in configuration file \"" << ClassGPEC::configuration_file_ << "\"." << endl;
#endif

    // Set time step.
    ops.SetPrefix(prefix_base + type_ + ".");

    delta_t_ = ops.Get<T>("Delta_t", "", delta_t_);

    // If type appears to be already filtered.
    if (type_.substr(0, 4) == "peak" || type_.substr(0, 4) == "mean")
      delta_t_ = T(86400);

    if (delta_t_ > T(0))
      Nt_day_ = get_time_index<T>(T(86400), delta_t_);
    else
      throw GPEC::Error("Negative or null time step in configuration.");

    ops.Close();
  }


  // Init.
  template<class T>
  void ClassContainerBase<T>::Init()
  {
    Ops::Ops ops(ClassGPEC::configuration_file_);

#ifdef GPEC_WITH_LOGGER
    *(ClassGPEC::log_ptr_) << Bgreen() << Fblack() << Info() << Reset()
                           << "Init static data for base container class." << endl;
#endif

    ops.SetPrefix("general.forecast.");
    forecast_min_ = ops.Get<int>("min");
    forecast_max_ = ops.Get<int>("max");

#ifdef GPEC_WITH_LOGGER
    *(ClassGPEC::log_ptr_) << Fgreen() << Debug() << Reset() << "forecast_min = " << forecast_min_ << endl;
    *(ClassGPEC::log_ptr_) << Fgreen() << Debug() << Reset() << "forecast_max = " << forecast_max_ << endl;
    *(ClassGPEC::log_ptr_) << Bgreen() << Fblack() << User(0, "EXPLANATION") << Reset()
                            << "Minimal and maximal forecast days used in GPEC : 0 means today, 1 tomorrow,"
                            << " ..., and negative means an analysis, change with care." << endl;
#endif

    ops.Close();
  }


  // Constructors.
  template<class T>
  ClassContainerBase<T>::ClassContainerBase()
    : ClassGPEC(), type_(""), name_(""), comment_(""),
      species_(""), date_begin_(DATE_END_DEFAULT), date_end_(DATE_BEGIN_DEFAULT), Nt_(0),
      Nt_day_(0), Nd_(0), delta_t_(T(0)), Nl_(0), value_min_(T(0)), value_max_(T(1.e9))
  {
    return;
  }


  template<class T>
  ClassContainerBase<T>::ClassContainerBase(const string &type,
                                            const string &name,
                                            const string &species,
                                            const string &prefix_base)
    : ClassGPEC(), type_(type), name_(name), comment_(""),
      species_(species), date_begin_(DATE_END_DEFAULT), date_end_(DATE_BEGIN_DEFAULT), Nt_(0),
      Nt_day_(0), Nd_(0), delta_t_(T(0)), Nl_(0), value_min_(T(0)), value_max_(T(1.e9))
  {
    construct(prefix_base);

    return;
  }


  template<class T>
  ClassContainerBase<T>::ClassContainerBase(const string &name,
                                            const string &species,
                                            const string &prefix_base)
    : ClassGPEC(), type_(""), name_(name), comment_(""),
      species_(species), date_begin_(DATE_END_DEFAULT), date_end_(DATE_BEGIN_DEFAULT), Nt_(0),
      Nt_day_(0), Nd_(0), delta_t_(T(0)), Nl_(0), value_min_(T(0)), value_max_(T(1.e9))
  {
    construct(prefix_base);

    return;
  }


  // Destructor.
  template<class T>
  ClassContainerBase<T>::~ClassContainerBase()
  {
    return;
  }


  // Clear data.
  template<class T>
  void ClassContainerBase<T>::Clear()
  {
    Nt_ = 0;
    Nd_ = 0;
    date_begin_.SetDate(DATE_END_DEFAULT);
    date_end_.SetDate(DATE_BEGIN_DEFAULT);
    forecast_.Clear();
    location_index_.Clear();
    data_.Clear();
    string comment = "cleared container";
    AddComment(comment);
  }


  template<class T>
  void ClassContainerBase<T>::PartialClear()
  {
    for (int t = 0; t < Nt_; t++)
      {
        location_index_(t).Clear();
        if (data_(t).GetSize() > 0)
          data_(t).Resize(1);
      }
  }


  // Copy constructor.
  template<class T>
  void ClassContainerBase<T>::Copy(const ClassContainerBase<T> &container, const string name)
  {
    // Name of copy may be different.
    name_ = name;

    // Everything else is identical.
    type_ = container.type_;
    species_ = container.species_;
    date_begin_ = container.date_begin_;
    date_end_ = container.date_end_;
    Nt_ = container.Nt_;
    Nd_ = container.Nd_;
    Nt_day_ = container.Nt_day_;
    delta_t_ = container.delta_t_;
    Nl_ = container.Nl_;
    location_index_ = container.location_index_;
    data_ = container.data_;
    forecast_ = container.forecast_;

    AddComment("copied from container \"" + container.name_ + "/" + container.type_ + "\"");

    ClassGPEC::Copy(container);
  }


  template<class T>
  bool ClassContainerBase<T>::HasSameStructure(const ClassContainerBase<T> &container) const
  {
    return Nl_ == container.Nl_;
  }

  // Checks if same time structure.
  template<class T>
  bool ClassContainerBase<T>::HasSameTimeStructure(const ClassContainerBase<T> &container) const
  {
    if (date_begin_ != container.date_begin_ ||
        Nt_ != container.Nt_ ||
        Nd_ != container.Nd_ ||
        Nt_day_ != container.Nt_day_)
      return false;

    return true;
  }


  template<class T>
  void ClassContainerBase<T>::SameTimeStructureAs(const ClassContainerBase<T> &container)
  {
    if (delta_t_ != container.delta_t_)
      throw GPEC::Error("Cannot put same time structure between containers with different time steps.");

    int t_beg(0), t_end(Nt_);

    if (date_begin_ < container.date_begin_)
      t_beg = GetDateIndex(container.date_begin_);

    if (date_end_ > container.date_end_)
      t_end = GetDateIndex(container.date_end_);

    RestrictDate(t_beg, t_end);

    Date date_beg(date_begin_), date_end(date_end_);

    if (date_begin_ > container.date_begin_)
      date_beg = container.date_begin_;
    if (date_end_ < container.date_end_)
      date_end = container.date_end_;

    ExtendDate(get_int_from_date(date_beg),
               get_int_from_date(date_end));
  }


  // Set methods.
  template<class T>
  void ClassContainerBase<T>::SetName(const string &name)
  {
    name_ = name;
  }


  template<class T>
  void ClassContainerBase<T>::SetType(const string &type)
  {
    type_ = type;
  }


  template<class T>
  void ClassContainerBase<T>::SetSpecies(const string &species)
  {
    species_ = species;
  }


  template<class T>
  void ClassContainerBase<T>::Set(const ClassContainerBase<T> &container, const string &prefix_base)
  {
    Clear();

    type_ = container.type_;
    name_ = container.name_;
    species_ = container.species_;

    construct(prefix_base);
  }


  template<class T> template<class U>
  void ClassContainerBase<T>::SetData(const Date &date_begin,
                                      const Date &date_end,
                                      const Vector1I &time,
                                      const Vector1I &location,
                                      const Vector1I &forecast,
                                      Vector<U> &data,
                                      const bool mask_negative_value)
  {
    if (data.GetSize() != time.GetSize() ||
        data.GetSize() != location.GetSize())
      throw GPEC::Error("Data vectors given in argument are not of equal size.");

    Clear();

    date_begin_ = date_begin;
    date_end_ = date_end;
    set_Nt();
    set_Nd();

    if (forecast.GetSize() != Nt_)
      throw GPEC::Error("Forecast vector should be of size " + to_str(Nt_) + ".");

    for (int i = 0; i < time.GetSize(); i++)
      check_index(0, Nt_, time(i));

    for (int i = 0; i < location.GetSize(); i++)
      check_index(0, Nl_, location(i));

    forecast_.Reallocate(Nt_);
    location_index_.Reallocate(Nt_);
    data_.Reallocate(Nt_);

    Vector1I Nl(Nt_);
    Nl.Zero();

    Vector<int> mask = check_data<U>(U(value_min_), U(value_max_), data);

    // Treat negative values as low values.
    if (! mask_negative_value)
      for (int i = 0; i < data.GetSize(); i++)
        if (mask(i) == -3)
          mask(i) = 1;

    // Negative values may appear for some sequential aggregation algorithms
    // that do not garantee positiveness. However, it usually appears because
    // of very low simulation values, and it should not be considered as an error.

    for (int i = 0; i < data.GetSize(); i++)
      if (mask(i) == 0)
        Nl(time(i))++;
      else
        {
#ifdef GPEC_WITH_LOGGER
          *(this->log_ptr_) << User(0, "BAD-VALUE") << "At time " << time(i) << " and location "
                        << location(i) << " got value " << data(i) << " ." << endl;
#endif

          if (mask(i) > 0)
            {
#ifdef GPEC_WITH_LOGGER
              *(this->log_ptr_) << User(0, "LOW-HIGH-VALUE") << "At time " << time(i) << " and location "
                            << location(i) << " got value " << data(i) << " ." << endl;
#endif
              Nl(time(i))++;
              data(i) = min<U>(max<U>(value_min_, data(i)), value_max_);
            }
        }

    for (int t = 0; t < Nt_; t++)
      if (Nl(t) > 0)
        {
          forecast_(t).Reallocate(1);
          check_index(forecast_min_, forecast_max_, forecast(t));
          forecast_(t)(0) = forecast(t);
          location_index_(t).Reallocate(Nl(t));
          data_(t).Reallocate(Nl(t));
        }

    Vector1I l(Nt_);
    l.Zero();
    for (int i = 0; i < data.GetSize(); i++)
      if (mask(i) >= 0)
        {
          int t = time(i);
          location_index_(t)(l(t)) = location(i);
          data_(t)(l(t)) = T(data(i));
          l(t)++;
        }

    for (int t = 0; t < Nt_; t++)
      for (int l = 1; l < location_index_(t).GetSize(); l++)
        if (location_index_(t)(l - 1) > location_index_(t)(l))
          {
            Sort(location_index_(t), data_(t));
            break;
          }

    // Location index not necessary if all locations are present.
    for (int t = 0; t < Nt_; t++)
      if (location_index_(t).GetSize() == Nl_)
        location_index_(t).Clear();
  }


  template<class T>
  void ClassContainerBase<T>::SetValueMin(const T &value_min)
  {
    value_min_ = value_min;
  }


  template<class T>
  void ClassContainerBase<T>::SetValueMax(const T &value_max)
  {
    value_max_ = value_max;
  }


  template<class T>
  void ClassContainerBase<T>::SetValue(const int &t,
                                       const int &location,
                                       const T &value)
  {
    if (value < value_min_ || value >= value_max_)
      return;

    check_index(0, Nt_, t);
    check_index(0, Nl_, location);

    int l = get_location_index(t, location);

    if (l < data_(t).GetSize())
      data_(t)(l) = value;
  }


  template<class T>
  void ClassContainerBase<T>::AddComment(const string add_comment)
  {
    comment_ += "At " + get_current_time() + ": " + add_comment + ".\n";
  }


  template<class T>
  void ClassContainerBase<T>::SetComment(const string &comment)
  {
    comment_ = comment;
  }


  // Check container.
  template<class T>
  void ClassContainerBase<T>::CheckContainer() const
  {
    if (Nt_ != location_index_.GetSize() ||
        Nt_ != data_.GetSize() ||
        Nt_ != forecast_.GetSize())
      throw GPEC::Error(string("Location index, forecast or data vector length does")
                                   + string(" not correspond to the number of time steps."));

    for (int t = 0; t < Nt_; t++)
      if (data_(t).GetSize() < Nl_)
        {
          if (location_index_(t).GetSize() != data_(t).GetSize()) 
            throw GPEC::Error("At time " + to_str(t) + ", sizes of location index and data vectors differ.");
          for (int l = 1; l < location_index_(t).GetSize(); l++)
            if (location_index_(t)(l - 1) >= location_index_(t)(l))
              throw GPEC::Error("At time " + to_str(t) + ", location indexes are not stored in ascending order.");
        }
      else if (data_(t).GetSize() == Nl_)
        {
          if (location_index_(t).GetSize() > 0)
            throw GPEC::Error("At time " + to_str(t) + ", location index vector should be empty.");
        }
      else
        throw GPEC::Error("At time " + to_str(t) + ", data vector greater than number of locations.");
  }


  // Find missing time steps.
  template<class T>
  void ClassContainerBase<T>::FindMissingTimeStep(Vector1I &missing_time_step) const
  {
#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Info() << "Looking for missing time steps ..." << endl;
#endif

    for (int t = 0; t < Nt_; t++)
      if (data_(t).GetSize() == 0)
        missing_time_step.PushBack(t);

#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Info(1) << "... found " << missing_time_step.GetSize() << " of them." << endl;
#endif
  }


  // Operations.
  template<class T>
  void ClassContainerBase<T>::Add(const T &value)
  {
    for (int t = 0; t < Nt_; t++)
      for (int l = 0; l < data_(t).GetSize(); l++)
        data_(t)(l) += value;
  }


  template<class T>
  void ClassContainerBase<T>::Mlt(const T &value)
  {
    for (int t = 0; t < Nt_; t++)
      for (int l = 0; l < data_(t).GetSize(); l++)
        data_(t)(l) *= value;
  }


  template<class T>
  void ClassContainerBase<T>::Set(const T &value)
  {
    for (int t = 0; t < Nt_; t++)
      for (int l = 0; l < data_(t).GetSize(); l++)
        data_(t)(l) = value;
  }


  // Print.
  template<class T>
  void ClassContainerBase<T>::Print(string &str, const bool &with_data) const
  {
    ostringstream sout;

    if (with_data)
      for (int t = 0; t < Nt_; t++)
        {
          Date date = date_begin_;
          date.AddSeconds(t * delta_t_);

          sout << setw(15) << left << date.GetDate("%y-%m-%d_%h");
          sout << setw(6) << left << "(" + to_str(t) + ")";

          Vector1T data;
          if (data_(t).GetSize() == Nl_)
            data = data_(t);
          else
            {
              data.Resize(Nl_);
              data.Fill(this->bad_value_);

              int loc(-1), l(0);
              while (++loc < Nl_ && l < location_index_(t).GetSize())
                if (location_index_(t)(l) == loc)
                  data(loc) = data_(t)(l++);
            }

          for (int loc = 0; loc < Nl_; loc++)
            sout << setw(12) << setprecision(8) << left << data(loc);
          sout << endl;
        }

    for (int t = 0; t < Nt_; t++)
      {
        Date date = date_begin_;
        date.AddSeconds(t * delta_t_);

        sout << setw(15) << left << date.GetDate("%y-%m-%d_%h");
        sout << setw(8) << right << "(" + to_str(t) + ")";
        sout << setw(8) << right << data_(t).GetSize() << "/"  << Nl_;
        if (forecast_(t).GetSize() > 0)
          sout << setw(4) << right << "D" << showpos << forecast_(t)(0);
        if (data_(t).GetSize() > 0)
          sout << setw(8) << right << noshowpos << fixed << setprecision(2) << Norm1(data_(t)) / T(data_(t).GetSize())
               << setw(8) << right << noshowpos << fixed << setprecision(2) << data_(t).GetNormInf();
        sout << endl;
      }

    sout << "Container: " << name_ << "/" << type_ << endl
         << "Species: " << this->species_ << endl;
    sout << "Date: " << date_begin_.GetDate("%y-%m-%d_%h") << " ---> "
         << date_end_.GetDate("%y-%m-%d_%h") << " (excluded)" << endl
         << "Time step = " << setprecision(1) << fixed << delta_t_ << " seconds" << endl
         << "Number of time steps = " << Nt_ << endl
         << "Number of days = " << Nd_ << endl;
    if (! comment_.empty())
      {
        sout << "Comment:" << endl;
        vector<string> ls = split(comment_, "\n");
        for (int i = 0; i < int(ls.size()); i++)
          sout << "  " << ls[i] << endl;
      }
    str = sout.str();
  }


  // Get methods.
  template<class T>
  string ClassContainerBase<T>::GetType() const
  {
    return type_;
  }


  template<class T>
  string ClassContainerBase<T>::GetName() const
  {
    return name_;
  }


  template<class T>
  string ClassContainerBase<T>::GetComment() const
  {
    return comment_;
  }


  template<class T>
  string ClassContainerBase<T>::GetSpecies() const
  {
    return species_;
  }


  template<class T>
  int ClassContainerBase<T>::GetLocationIndex(const int &t, const int &location) const
  {
    return get_location_index(t, location);
  }


  template<class T>
  inline int ClassContainerBase<T>::GetNlTotal() const
  {
    int NlTotal(0);
    for (int t = 0; t < data_.GetSize(); t++)
      NlTotal += data_(t).GetSize();
    return NlTotal;
  }


  template<class T>
  int ClassContainerBase<T>::GetNl(const int &t) const
  {
    return data_(t).GetSize();
  }


  template<class T>
  int ClassContainerBase<T>::GetNlocation() const
  {
    return Nl_;
  }


  template<class T>
  int ClassContainerBase<T>::GetNt() const
  {
    return Nt_;
  }


  template<class T>
  int ClassContainerBase<T>::GetNtDay() const
  {
    return Nt_day_;
  }


  template<class T>
  int ClassContainerBase<T>::GetNd() const
  {
    return Nd_;
  }


  template<class T>
  T ClassContainerBase<T>::GetDeltat() const
  {
    return delta_t_;
  }


  template<class T>
  int ClassContainerBase<T>::GetDateIndex(const Date &date) const
  {
    return get_time_index<T>(date.GetSecondsFrom(date_begin_), delta_t_);
  }


  template<class T>
  Date ClassContainerBase<T>::GetDate(const int &t) const
  {
    Date date(date_begin_);
    date.AddSeconds(t * delta_t_);

    return date;
  }


  template<class T>
  Date ClassContainerBase<T>::GetDateBegin() const
  {
    return date_begin_;
  }


  template<class T>
  Date ClassContainerBase<T>::GetDateEnd() const
  {
    return date_end_;
  }


  template<class T>
  string ClassContainerBase<T>::GetDateStr(const int &t,
                                           const string format) const
  {
    return GetDate(t).GetDate(format);
  }


  template<class T>
  int ClassContainerBase<T>::GetDateIndexStr(const string &date_str) const
  {
    Date date(date_str);
    return GetDateIndex(date);
  }


  template<class T>
  string ClassContainerBase<T>::GetDateBeginStr(const string format) const
  {
    return date_begin_.GetDate(format);
  }


  template<class T>
  string ClassContainerBase<T>::GetDateEndStr(const string format) const
  {
    return date_end_.GetDate(format);
  }


  template<class T>
  T ClassContainerBase<T>::GetValue(const int &t, const int &location) const
  {
    check_index(0, Nt_, t);
    check_index(0, Nl_, location);

    int l = get_location_index(t, location);

    if (l < data_(t).GetSize())
      return data_(t)(l);
    else
      return this->bad_value_;
  }


  template<class T>
  inline T ClassContainerBase<T>::GetValueRaw(const int &t, const int &l) const
  {
    return data_(t)(l);
  }


  template<class T>
  int ClassContainerBase<T>::GetForecastDay(const int &t) const
  {
    check_index(0, Nt_, t);
    if (forecast_(t).GetSize() > 0)
      return forecast_(t)(0);
    else
      return int(this->bad_value_);
  }


  template<class T>
  T ClassContainerBase<T>::GetForecastDayMean(int t_begin, int t_end) const
  {
    if (t_end == 0)
      t_end = Nt_;

    int sum(0), count(0);
    for (int t = t_begin; t < t_end; t++)
      if (forecast_(t).GetSize() > 0)
        {
          sum += forecast_(t)(0);
          count++;
        }

    if (count > 0)
      return T(sum) / T(count);
    else
      return this->bad_value_;
  }


  template<class T>
  inline void ClassContainerBase<T>::CollectData(Vector1I &time,
                                                 Vector1I& location,
                                                 Vector1T &data) const
  {
    int Nsize(Nt_ * Nl_);

    time.Reallocate(Nsize);
    location.Reallocate(Nsize);
    data.Reallocate(Nsize);

    int i(0);
    for (int t = 0; t < Nt_; t++)
      if (data_(t).GetSize() == Nl_)
        for (int l = 0; l < Nl_; l++)
          {
            time(i) = t;
            location(i) = l;
            data(i++) = data_(t)(l);
          }
      else
        for (int l = 0; l < data_(t).GetSize(); l++)
          {
            time(i) = t;
            location(i) = location_index_(t)(l);
            data(i++) = data_(t)(l);
          }

    time.Resize(i);
    location.Resize(i);
    data.Resize(i);
  }


  template<class T>
  inline void ClassContainerBase<T>::CollectDataPeriod(const int &t_begin,
                                                       const int &t_end,
                                                       Vector1I &time,
                                                       Vector1I& location,
                                                       Vector1T &data) const
  {
    time.Clear();
    location.Clear();
    data.Clear();

    check_index(0, Nt_, t_begin);
    check_index(t_begin + 1, Nt_ + 1, t_end);

    int Nsize((t_end - t_begin) * Nl_);
    time.Reallocate(Nsize);
    location.Reallocate(Nsize);
    data.Reallocate(Nsize);

    int i(0);
    for (int t = t_begin; t < t_end; t++)
      if (data_(t).GetSize() == Nl_)
        for (int l = 0; l < data_(t).GetSize(); l++)
          {
            time(i) = t;
            location(i) = l;
            data(i++) = data_(t)(l);
          }
      else
        for (int l = 0; l < data_(t).GetSize(); l++)
          {
            time(i) = t;
            location(i) = location_index_(t)(l);
            data(i++) = data_(t)(l);
          }

    time.Resize(i);
    location.Resize(i);
    data.Resize(i);
  }


  template<class T>
  inline void ClassContainerBase<T>::CollectDataLocation(const int &location,
                                                         Vector1I &time,
                                                         Vector1T &data) const
  {
    time.Clear();
    data.Clear();

    check_index(0, Nl_, location);

    time.Reallocate(Nt_);
    data.Reallocate(Nt_);

    int i(0);
    for (int t = 0; t < Nt_; t++)
      {
        int l = get_location_index(t, location);

        if (l < data_(t).GetSize())
          {
            time(i) = t;
            data(i++) = data_(t)(l);
          }
      }

    time.Resize(i);
    data.Resize(i);    
  }


  template<class T>
  inline void ClassContainerBase<T>::CollectDataTime(const int &t,
                                                     Vector1I &location,
                                                     Vector1T &data) const
  {
    location.Clear();
    data.Clear();

    if (data_(t).GetSize() == Nl_)
      {
        location.Reallocate(Nl_);
        location.Fill();
      }
    else
      location = location_index_(t);

    data = data_(t);
  }


  template<class T>
  inline void ClassContainerBase<T>::CollectDataTimeAtLocation(const int &t,
                                                               const Vector1I &location,
                                                               Vector1I &location_index,
                                                               Vector1T &data) const
  {
    location_index.Reallocate(location.GetSize());
    data.Reallocate(location.GetSize());
    int j(0);

    if (data_(t).GetSize() == Nl_)
      for (int i = 0; i < location.GetSize(); i++)
        {
          location_index(j) = location(i);
          data(j++) = data_(t)(location(i));
        }
    else
      {
        int i(0), l(0);
        while (l < location_index_(t).GetSize() && i < location.GetSize())
          if (location_index_(t)(l) > location(i))
            i++;
          else if (location_index_(t)(l) < location(i))
            l++;
          else
            {
              location_index(j) = location(i);
              data(j++) = data_(t)(l);
              i++;
              l++;
            }
      }

    location_index.Resize(j);
    data.Resize(j);
  }


  template<class T>
  void ClassContainerBase<T>::CollectDataRaw(Vector1T &data,
                                             const bool with_bad_value) const
  {
    Vector1I time, location;
    Vector1T data_tmp;
    CollectData(time, location, data_tmp);

    if (with_bad_value)
      {
        data.Reallocate(Nt_ * Nl_);
        data.Fill(this->bad_value_);

        for (int i = 0; i < data_tmp.GetSize(); i++)
          data(Nl_ * time(i) + location(i)) = data_tmp(i);
      }
    else
      data = data_tmp;
  }


  template<class T>
  void ClassContainerBase<T>::CollectDataLocationRaw(const int &location,
                                                     Vector1T &data,
                                                     const bool with_bad_value) const
  {
    Vector1I time;
    Vector1T data_tmp;
    CollectDataLocation(location, time, data_tmp);

    if (with_bad_value)
      {
        data.Reallocate(Nt_);
        data.Fill(this->bad_value_);

        for (int i = 0; i < data_tmp.GetSize(); i++)
          data(time(i)) = data_tmp(i);
      }
    else
      data = data_tmp;
  }


  template<class T>
  void ClassContainerBase<T>::CollectDataTimeRaw(const int &t,
                                                 Vector1T &data,
                                                 const bool with_bad_value) const
  {
    Vector1I location;
    Vector1T data_tmp;
    CollectDataTime(t, location, data_tmp);

    if (with_bad_value)
      {
        data.Reallocate(Nl_);
        data.Fill(this->bad_value_);

        for (int i = 0; i < data_tmp.GetSize(); i++)
          data(location(i)) = data_tmp(i);
      }
    else
      data = data_tmp;
  }


  template<class T>
  void ClassContainerBase<T>::CollectDataPeriodRaw(const int &t_begin,
                                                   const int &t_end,
                                                   Vector1T &data,
                                                   const bool with_bad_value) const
  {
    Vector1I time, location;
    Vector1T data_tmp;
    CollectDataPeriod(t_begin, t_end, time, location, data_tmp);

    if (with_bad_value)
      {
        data.Reallocate((t_end - t_begin) * Nl_);
        data.Fill(this->bad_value_);

        for (int i = 0; i < data_tmp.GetSize(); i++)
          data(Nl_ * (time(i) - t_begin) + location(i)) = data_tmp(i);
      }
    else
      data = data_tmp;
  }


  // Remove or keep locations.
  template<class T>
  void ClassContainerBase<T>::RemoveTimeStep(const Vector1I &time_step)
  {
    // Abruptly remove some time steps.

    for (int i = 0; i < time_step.GetSize(); i++)
      {
        int t = time_step(i);
        forecast_(t).Clear();
        location_index_(t).Clear();
        data_(t).Clear();
      }
  }


  template<class T>
  void ClassContainerBase<T>::RemoveLocation(const Vector1I &location)
  {
    Vector1I location_ascending = location;
    Sort(location_ascending);

    Vector1I location_keep(Nl_);

    int l(-1), i(0), j(0);
    while (++l < Nl_)
      if (i < location_ascending.GetSize())
        {
          if (l == location_ascending(i))
            i++;
          else
            location_keep(j++) = l;
        }
      else
        location_keep(j++) = l;

    location_keep.Resize(j);

    KeepLocation(location_keep);
  }


  template<class T>
  void  ClassContainerBase<T>::KeepLocation(const Vector1I &location)
  {
    Vector1I location_ascending = location;
    Sort(location_ascending);

    for (int t = 0; t < Nt_; t++)
      {
        Vector1I location_index_new(data_(t).GetSize());
        Vector1T data_new(data_(t).GetSize());

        int i(0), count(0);
        if (data_(t).GetSize() == Nl_)
          {
            int l(-1);
            while (++l < Nl_ && i < location_ascending.GetSize())
              if (l == location_ascending(i))
                {
                  location_index_new(count) = l;
                  data_new(count++) = data_(t)(l);
                  i++;
                }
          }
        else
          {
            int l(0);
            while (l < location_index_(t).GetSize() && i < location_ascending.GetSize())
              if (location_index_(t)(l) > location_ascending(i))
                i++;
              else if (location_index_(t)(l) < location_ascending(i))
                l++;
              else
                {
                  location_index_new(count) = location_index_(t)(l);
                  data_new(count++) = data_(t)(l);
                  l++;
                  i++;
                }
          }

        location_index_new.Resize(count);
        data_new.Resize(count);

        if (count == 0)
          forecast_(t).Clear();

        location_index_(t) = location_index_new;
        data_(t) = data_new;
      }

    AddComment("kept " + to_str(location.GetSize()) + " locations");
  }


  // Restrict or extend date range.
  template<class T>
  void  ClassContainerBase<T>::RestrictDate(const int &t_begin, const int &t_end)
  {
    check_index(0, Nt_, t_begin);
    check_index(t_begin + 1, Nt_ + 1, t_end);

    Nt_ = t_end - t_begin;
    date_begin_.AddSeconds(t_begin * delta_t_);
    date_end_ = date_begin_;
    date_end_.AddSeconds(Nt_ * delta_t_);
    set_Nd();

    for (int t = 0; t < Nt_; t++)
      {
        int t_offset(t_begin + t);
        forecast_(t) = forecast_(t_offset);
        location_index_(t) = location_index_(t_offset);
        data_(t) = data_(t_offset);
      }

    forecast_.Resize(Nt_);
    location_index_.Resize(Nt_);
    data_.Resize(Nt_);

    AddComment("restricted container between time steps " + to_str(t_begin)
               + " and " + to_str(t_end) + " (excluded)");
  }


  template<class T>
  void ClassContainerBase<T>::ExtendDate(const int &date_begin_int, const int &date_end_int)
  {
    Date date_begin_ext = get_date_from_int(date_begin_int);
    Date date_end_ext = get_date_from_int(date_end_int);

    if (date_begin_ext > date_begin_ && date_end_ext <= date_end_)
      return;

    int t_old = get_time_index<T>(date_begin_.GetSecondsFrom(date_begin_ext), delta_t_);
    int Nt_old(Nt_);

    date_begin_ = date_begin_ext;
    date_end_ = date_end_ext;
    set_Nt();
    set_Nd();

    forecast_.Resize(Nt_);
    location_index_.Resize(Nt_);
    data_.Resize(Nt_);
    if (t_old != 0)
      for (int t = Nt_old - 1; t >= 0; t--)
        {
          int t_new(t + t_old);
          forecast_(t_new) = forecast_(t);
          location_index_(t_new) = location_index_(t);
          data_(t_new) = data_(t);
          forecast_(t).Clear();
          location_index_(t).Clear();
          data_(t).Clear();
        }

    AddComment("extended container between dates " + date_begin_ext.GetDate("%y-%m-%d_%h")
               + " and " + date_end_ext.GetDate("%y-%m-%d_%h") + " (excluded)");
  }


  // Merge.
  template<class T>
  void ClassContainerBase<T>::Merge(const ClassContainerBase<T> &container,
                                    const bool overwrite)
  {
    if (delta_t_ != container.delta_t_)
      throw GPEC::Error("Not equal time step between containers.");
    if (Nl_ != container.Nl_)
      throw GPEC::Error("Not equal number of locations between containers.");

    Date date_begin_merge = min<Date>(date_begin_, container.date_begin_);
    Date date_end_merge = max<Date>(date_end_, container.date_end_);

    // Extend current container to merge dates.
    ExtendDate(get_int_from_date(date_begin_merge),
               get_int_from_date(date_end_merge));

    // Merge with other container.
    int t_merge(0);
    Date date_running(date_begin_merge);
    while (date_running < date_end_merge)
      {
        if (date_running >= container.date_begin_ &&
            date_running < container.date_end_)
          {
            int t = container.GetDateIndex(date_running);

            if (data_(t_merge).GetSize() == 0)
              {
                forecast_(t_merge) = container.forecast_(t);
                location_index_(t_merge) = container.location_index_(t);
                data_(t_merge) = container.data_(t);
              }
            else
              {
                int l(0), l_merge(0), count(0);

                if (container.data_(t).GetSize() == Nl_ && data_(t_merge).GetSize() == Nl_ && overwrite)
                  {
                    data_(t_merge) = container.data_(t);
                    count = Nl_;
                  }
                else if (container.data_(t).GetSize() == Nl_)
                  {
                    while (l < Nl_ && l_merge < location_index_(t_merge).GetSize())
                      if (l < location_index_(t_merge)(l_merge))
                        {
                          location_index_(t_merge).PushBack(l);
                          data_(t_merge).PushBack(container.data_(t)(l));
                          l++;
                          count++;
                        }
                      else if (l == location_index_(t_merge)(l_merge))
                        {
                          if (overwrite)
                            {
                              data_(t_merge)(l_merge) = container.data_(t)(l);
                              count++;
                            }
                          l++;
                          l_merge++;
                        }
                  }
                else if (data_(t_merge).GetSize() == Nl_ && overwrite)
                  {
                    for (l = 0; l < container.data_(t).GetSize(); l++)
                      data_(t)(container.location_index_(t)(l)) = container.data_(t)(l);
                    count = container.data_(t).GetSize();
                  }
                else
                  {
                    while (l < container.location_index_(t).GetSize() &&
                           l_merge < location_index_(t_merge).GetSize())
                      if (container.location_index_(t)(l) > location_index_(t_merge)(l_merge))
                        l_merge++;
                      else if (container.location_index_(t)(l) < location_index_(t_merge)(l_merge))
                        {
                          location_index_(t_merge).PushBack(container.location_index_(t)(l));
                          data_(t_merge).PushBack(container.data_(t)(l));
                          l++;
                          count++;
                        }
                      else
                        {
                          if (overwrite)
                            {
                              data_(t_merge)(l_merge) = container.data_(t)(l);
                              count++;
                            }
                          l++;
                          l_merge++;
                        }
                  }

                if (count > 0)
                  forecast_(t_merge).PushBack(container.forecast_(t)(0));
              }

            if (data_(t_merge).GetSize() > 0)
              for (int l = 1; l < location_index_(t_merge).GetSize(); l++)
                if (location_index_(t_merge)(l - 1) > location_index_(t_merge)(l))
                  {
                    Sort(location_index_(t_merge), data_(t_merge));
                    break;
                  }
          }
        t_merge++;
        date_running.AddSeconds(delta_t_);
      }

    string add_comment = "merged with container \"" + container.name_ + "/" + container.type_ + "\"";
    if (overwrite)
      add_comment += ", with option overwrite";
    AddComment(add_comment);
  }


  // Filter data.
  template<class T> template<class F>
  void ClassContainerBase<T>::FilterTime(const int h_frequency, const int h_offset,
                                         string prefix_type,
                                         const int h_begin, int h_end, const T ratio)
  {
    if (h_frequency == 1)
      return;

    if (abs(h_offset) >= h_frequency)
      throw GPEC::Error("The absolute offset filter cannot equal or exceed the filter frequency.");

    // Only relevant if data frequency is less than daily.
    T delta_t_new = h_frequency * delta_t_;
    if (delta_t_new > T(86400))
      throw GPEC::Error("In GPEC, the time step cannot exceed one day.");

    if (h_end == 0)
      h_end = h_frequency;

    int count_min = int(ratio * h_frequency);

    int t_begin(0), t_end(Nt_);
    int Nt_new = int((Nt_ - h_offset) / h_frequency + 1.e-5);
    if ((Nt_ - h_offset) % h_frequency > count_min)
      Nt_new++;
    else
      t_end -= (Nt_ - h_offset) % h_frequency;

    date_begin_.AddSeconds((h_frequency / 2 + h_offset) * delta_t_);

    if (h_offset > 0)
      t_begin = h_offset;
    else if (h_offset < 0)
      if ((h_frequency + h_offset) < count_min)
        {
          Nt_new--;
          t_begin = h_frequency + h_offset;
          date_begin_.AddSeconds(h_frequency * delta_t_ );
        }

    Vector2I forecast_new(Nt_new);
    Vector2I location_index_new(Nt_new);
    Vector2T data_new(Nt_new);

    for (int t = 0; t < Nt_new; t++)
      {
        forecast_new(t).Reallocate(1);
        location_index_new(t).Reallocate(Nl_);
        data_new(t).Reallocate(Nl_);
      }

    Vector1I index(Nt_new);
    index.Zero();

    F filter;

    for (int location = 0; location < Nl_; location++)
      {
        Vector1I forecast_frequency(Nt_new);
        forecast_frequency.Zero();

        Vector1T data_frequency(Nt_new);
        data_frequency.Fill(filter.Begin());

        Vector1I count(Nt_new);
        count.Zero();

        int h(0), h_count(0);
        for (int t = t_begin; t < t_end; t++)
          {
            int l = get_location_index(t, location);

            if (l < data_(t).GetSize())
              if (h_count >= h_begin)
                if (h_count < h_end)
                  {
                    forecast_frequency(h) += forecast_(t)(0);
                    filter.Do(data_(t)(l), data_frequency(h));
                    count(h)++;
                  }

            h_count++;
            if (h_count == h_frequency)
              {
                h++;
                h_count = 0;
              }
          }

        for (int h = 0; h < Nt_new; h++)
          if (count(h) > count_min)
            {
              forecast_new(h)(0) = int(rint(T(forecast_frequency(h)) / T(count(h))));
              location_index_new(h)(index(h)) = location;
              data_new(h)(index(h)) = filter.End(count(h), data_frequency(h));
              index(h)++;
            }
      }

    for (int t = 0; t < Nt_new; t++)
      {
        if (index(t) == 0)
          forecast_new(t).Clear();

        if (index(t) == Nl_)
          location_index_new(t).Clear();
        else
          location_index_new(t).Resize(index(t));

        data_new(t).Resize(index(t));
      }

    // Update attributes.
    date_end_ = date_begin_;
    date_end_.AddSeconds(Nt_new * delta_t_new);
    Nt_ = Nt_new;
    delta_t_ = delta_t_new;
    set_Nd();
    Nt_day_ = int(T(86400) / delta_t_);

    forecast_ = forecast_new;
    location_index_ = location_index_new;
    data_ = data_new;

    AddComment("filtered " + filter.Name() + " with frequency " + to_str(h_frequency)
               + " and offset " + to_str(h_offset) + ", between " + to_str(h_begin) + " and " + to_str(h_end));

    // Prefix type with something to say we filtered time.
    if (prefix_type.empty())
      {
        prefix_type = filter.Name();
        if (h_frequency < 24)
          prefix_type += to_str(h_frequency);
      }

    if (prefix_type != "-")
      type_ = prefix_type + type_;
  }


  // Filter time steps with two few data.
  template<class T>
  void ClassContainerBase<T>::FilterTimeStep(const int location_min)
  {
#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Info() << "Filter time steps with"
                      << " number of locations < " << location_min << endl;
#endif

    int count_filtered(0);

    for (int t = 0; t < Nt_; t++)
      if (data_(t).GetSize() < location_min)
        {
          count_filtered++;

#ifdef GPEC_WITH_LOGGER
          *(this->log_ptr_) << Debug(2) << "Filter time step "
                            << t << " with number of locations = "
                            << data_(t).GetSize() << endl;
#endif

          location_index_(t).Clear();
          data_(t).Clear();
          forecast_(t).Clear();
        }

#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Debug(1) << "Filtered "
                      << count_filtered << " time steps." << endl;
#endif
    AddComment("filtered " + to_str(count_filtered) + " time steps");
  }


  // Filter value.
  template<class T>
  void ClassContainerBase<T>::FilterValue(T value_min, T value_max,
                                          const bool keep_threshold)
  {
    value_min = value_min == T(0) ? value_min_ : value_min;
    value_max = value_max == T(1e9) ? value_max_ : value_max;

    if (value_min < value_min_ || value_max > value_max_)
      return;

#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Debug() << "Filter values below "
                      << value_min << " and above " << value_max
                      << ", filtered values are "
                      << (keep_threshold ? "REPLACED by min or max" : "REMOVED")
                      << endl;
#endif
    int count_filtered(0);

    for (int t = 0; t < Nt_; t++)
      {
        int Nl = data_(t).GetSize();

        int i(0);

        if (Nl > 0)
          {
            if (location_index_(t).GetSize() == 0)
              {
                location_index_(t).Reallocate(Nl);
                location_index_(t).Fill();
              }

            for (int l = 0; l < Nl; l++)
              if (data_(t)(l) < value_min)
                {
                  count_filtered++;
#ifdef GPEC_WITH_LOGGER
                  *(this->log_ptr_) << Debug(2) << "Filter value "
                                    << setprecision(2) << fixed << data_(t)(l) << " at location "
                                    << l << " and time " << t << endl;
#endif
                  if (keep_threshold)
                    {
                      location_index_(t)(i) = location_index_(t)(l);
                      data_(t)(i++) = value_min;
                    }
                }
              else if (data_(t)(l) > value_max)
                {
                  count_filtered++;
#ifdef GPEC_WITH_LOGGER
                  *(this->log_ptr_) << Debug(2) << "Filter value "
                                    << setprecision(2) << fixed << data_(t)(l) << " at location "
                                    << l << " and time " << t << endl;
#endif
                  if (keep_threshold)
                    {
                      location_index_(t)(i) = location_index_(t)(l);
                      data_(t)(i++) = value_max;
                    }
                }
              else
                {
                  location_index_(t)(i) = location_index_(t)(l);
                  data_(t)(i++) = data_(t)(l);
                }

            location_index_(t).Resize(i);
            data_(t).Resize(i);

            if (data_(t).GetSize() == 0)
              forecast_(t).Clear();

            if (data_(t).GetSize() == Nl_)
              location_index_(t).Clear();
          }
      }

#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Debug(1) << "Filtered "
                      << count_filtered << " values." << endl;
#endif
    AddComment("filtered " + to_str(count_filtered) + " values below "
               + to_str(value_min) + " and above " + to_str(value_max));
  }


  // Filter values below or above given multiple of standard deviation.
  template<class T>
  void ClassContainerBase<T>::FilterStandardDeviation(const T multiple_std,
                                                      const int location_min,
                                                      int t_begin, int t_end,
                                                      const bool keep_threshold)
  {
    if (multiple_std <= T(0))
      throw GPEC::Error("Standard deviation multiple should be strictly positive");

    check_index(0, Nt_, t_begin);

    t_end = t_end == 0 ? Nt_ : t_end;

#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Info() << "Filter values with "
                      << multiple_std << " x standard deviation between timesteps ["
                      << t_begin << ", " << t_end << "[" << endl;
#endif

    int count_filtered(0);

    // At each time step, compute the mean and standard deviation value
    // and remove values which are superior or inferior to mean +- std * multiple_std.
    for (int t = t_begin; t < t_end; t++)
      if (data_(t).GetSize() > location_min)
        {
          int Nl = data_(t).GetSize();

          // Compute mean and std.
          T mean = Norm1(data_(t)) / T(Nl);
          T std(T(0));
          for (int l = 0; l < Nl; l++)
            {
              T tmp = data_(t)(l) - mean;
              std += tmp * tmp;
            }
          std = sqrt(std / T(Nl));

          T min = mean - std * multiple_std;
          min = min < T(0) ? T(0) : min; 
          T max = mean + std * multiple_std;

#ifdef GPEC_WITH_LOGGER
          *(this->log_ptr_) << Debug(2) << "At time t = "
                            << t << ", filter values below "
                            << setprecision(2) << fixed << min
                            << " and above " << setprecision(2)
                            << fixed << max << endl;
#endif

          // Now filter.
          int i(0);

          if (location_index_(t).GetSize() == 0)
            {
              location_index_(t).Reallocate(Nl);
              location_index_(t).Fill();
            }

          for (int l = 0; l < Nl; l++)
            if (data_(t)(l) < min)
              {
                count_filtered++;

                if (keep_threshold)
                  {
                    location_index_(t)(i) = location_index_(t)(l);
                    data_(t)(i++) = min;
                  }
#ifdef GPEC_WITH_LOGGER
                *(this->log_ptr_) << Debug(3) << "Filter value "
                                  << setprecision(2) << fixed << data_(t)(l) << " at location "
                                  << l << " and time " << t << endl;
#endif
              }
            else if (data_(t)(l) > max)
              {
                count_filtered++;

                if (keep_threshold)
                  {
                    location_index_(t)(i) = location_index_(t)(l);
                    data_(t)(i++) = max;
                  }
#ifdef GPEC_WITH_LOGGER
                *(this->log_ptr_) << Debug(3) << "Filter value "
                                  << setprecision(2) << fixed << data_(t)(l) << " at location "
                                  << l << " and time " << t << endl;
#endif
              }
            else
              {
                location_index_(t)(i) = location_index_(t)(l);
                data_(t)(i++) = data_(t)(l);
              }

          location_index_(t).Resize(i);
          data_(t).Resize(i);

          if (data_(t).GetSize() == 0)
            forecast_(t).Clear();

          if (data_(t).GetSize() == Nl_)
            location_index_(t).Clear();
        }

#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Debug(1) << "Filtered "
                      << count_filtered << " values." << endl;
#endif
    AddComment("filtered " + to_str(count_filtered) + " values with"
               + " standard deviation times " + to_str(multiple_std));
  }


  // Filter values below or above given percentile.
  template<class T>
  void ClassContainerBase<T>::FilterPercentile(const T &percentile, int t_begin, int t_end,
                                               const bool keep_threshold)
  {
    if (percentile <= T(0) || percentile >= T(1))
      throw GPEC::Error("Percentile should be between ]0, 1[, got " + to_str(percentile));

    check_index(0, Nt_, t_begin);

    t_end = t_end == 0 ? Nt_ : t_end;

#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Info() << "Filter values with percentile "
                      << percentile << " between timesteps ["
                      << t_begin << ", " << t_end << "[" << endl;
#endif

    int count_filtered(0);

    // At each time step, compute the percentile value and remove values which are
    // superior if percentile > 0.5, inferior otherwise.
    for (int t = t_begin; t < t_end; t++)
      if (data_(t).GetSize() > 0)
        {
          int Nl = data_(t).GetSize();

          Vector<T> data_ascending = data_(t);

          Sort(data_ascending);

          Vector<T> data_percentile_rank(Nl);

          for (int l = 0; l < Nl; l++)
            data_percentile_rank(l) = (T(l) - T(0.5)) / T(Nl);

          T percentile_value;

          if (percentile < data_percentile_rank(0))
            percentile_value = data_ascending(0);
          else if (percentile > data_percentile_rank(Nl - 1))
            percentile_value = data_ascending(Nl - 1);
          else
            {
              int k(0);
              for (k = 0; k < Nl; k++)
                if (percentile < data_percentile_rank(k))
                  {
                    k--;
                    break;
                  }
                  
              percentile_value = data_ascending(k)
                + (percentile - data_percentile_rank(k))
                / (data_percentile_rank(k + 1) - data_percentile_rank(k))
                * (data_ascending(k + 1) - data_ascending(k));
            }

#ifdef GPEC_WITH_LOGGER
          *(this->log_ptr_) << Debug(2) << "At time t = "
                            << t << ", computed percentile value = "
                            << percentile_value << endl;
#endif

          // Now filter.
          int i(0);

          if (location_index_(t).GetSize() == 0)
            {
              location_index_(t).Reallocate(Nl);
              location_index_(t).Fill();
            }

          for (int l = 0; l < Nl; l++)
            if (percentile < T(0.5) && data_(t)(l) < percentile_value)
              {
                count_filtered++;

                if (keep_threshold)
                  {
                    location_index_(t)(i) = location_index_(t)(l);
                    data_(t)(i++) = percentile_value;
                  }
#ifdef GPEC_WITH_LOGGER
                *(this->log_ptr_) << Debug(3) << "Filter value "
                                  << setprecision(2) << fixed << data_(t)(l) << " at location "
                                  << l << " and time " << t << endl;
#endif
              }
            else if (percentile >= T(0.5) && data_(t)(l) > percentile_value)
              {
                count_filtered++;

                if (keep_threshold)
                  {
                    location_index_(t)(i) = location_index_(t)(l);
                    data_(t)(i++) = percentile_value;
                  }
#ifdef GPEC_WITH_LOGGER
                *(this->log_ptr_) << Debug(3) << "Filter value "
                                  << setprecision(2) << fixed << data_(t)(l) << " at location "
                                  << l << " and time " << t << endl;
#endif
              }
            else
              {
                location_index_(t)(i) = location_index_(t)(l);
                data_(t)(i++) = data_(t)(l);
              }

          location_index_(t).Resize(i);
          data_(t).Resize(i);

          if (data_(t).GetSize() == 0)
            forecast_(t).Clear();

          if (data_(t).GetSize() == Nl_)
            location_index_(t).Clear();
        }

#ifdef GPEC_WITH_LOGGER
    *(this->log_ptr_) << Debug(1) << "Filtered "
                      << count_filtered << " values." << endl;
#endif
    AddComment("filtered " + to_str(count_filtered) + " values with"
               + " percentile " + to_str(percentile));
  }


  // Filter value below a threshold.
  template<class T>
  void ClassContainerBase<T>::FilterThresholdExceedance(const T &threshold,
                                                        const bool reduce_location)
  {
    if (threshold < value_min_)
      return;

    if (threshold >= value_max_)
      throw GPEC::Error("Threshold value is greater than the container maximum admissible value.");

    for (int t = 0; t < Nt_; t++)
      {
        int Nl = data_(t).GetSize();

        int i(0);

        if (Nl > 0)
          {
            if (location_index_(t).GetSize() == 0)
              {
                location_index_(t).Reallocate(Nl);
                location_index_(t).Fill();
              }

            if (reduce_location)
              {
                for (int l = 0; l < Nl; l++)
                  if (data_(t)(l) >= threshold)
                    {
                      location_index_(t)(i) = location_index_(t)(l);
                      data_(t)(i++) = data_(t)(l);
                    }

                location_index_(t).Resize(i);
                data_(t).Resize(i);
              }
            else
              for (int l = 0; l < Nl; l++)
                if (data_(t)(l) < threshold)
                  data_(t)(l) = T(0);

            if (data_(t).GetSize() == 0)
              forecast_(t).Clear();

            if (data_(t).GetSize() == Nl_)
              location_index_(t).Clear();
          }
      }
  }


  // MPI Communication.
#ifdef GPEC_WITH_MPI
  template<class T>
  void ClassContainerBase<T>::MPI_SendRecv(const int &send_rank, const int &recv_rank)
  {
    int isize, itag(0);
    if (this->rank_ == recv_rank)
      {
        MPI::COMM_WORLD.Recv(&isize, 1, MPI::INT, send_rank, itag++);
        char *str = new char[isize];
        MPI::COMM_WORLD.Recv(str, isize, MPI::CHAR, send_rank, itag++);
        type_ = str;

        MPI::COMM_WORLD.Recv(&isize, 1, MPI::INT, send_rank, itag++);
        str = new char[isize];
        MPI::COMM_WORLD.Recv(str, isize, MPI::CHAR, send_rank, itag++);
        name_ = str;

        MPI::COMM_WORLD.Recv(&isize, 1, MPI::INT, send_rank, itag++);
        str = new char[isize];
        MPI::COMM_WORLD.Recv(str, isize, MPI::CHAR, send_rank, itag++);
        species_ = str;

        MPI::COMM_WORLD.Recv(&isize, 1, MPI::INT, send_rank, itag++);
        str = new char[isize];
        MPI::COMM_WORLD.Recv(str, isize, MPI::CHAR, send_rank, itag++);
        comment_ = str;

        delete[] str;

        MPI::COMM_WORLD.Recv(&isize, 1, MPI::INT, send_rank, itag++);
        date_begin_ = get_date_from_int(isize);
        MPI::COMM_WORLD.Recv(&isize, 1, MPI::INT, send_rank, itag++);
        date_end_ = get_date_from_int(isize);

        MPI::COMM_WORLD.Recv(&Nt_day_, 1, MPI::INT, send_rank, itag++);
        MPI::COMM_WORLD.Recv(&Nt_, 1, MPI::INT, send_rank, itag++);
        MPI::COMM_WORLD.Recv(&Nd_, 1, MPI::INT, send_rank, itag++);
        MPI::COMM_WORLD.Recv(&delta_t_, 1, GPEC_MPI_REAL, send_rank, itag++);
        MPI::COMM_WORLD.Recv(&Nl_, 1, MPI::INT, send_rank, itag++);

        MPI::COMM_WORLD.Recv(&value_min_, 1, GPEC_MPI_REAL, send_rank, itag++);
        MPI::COMM_WORLD.Recv(&value_max_, 1, GPEC_MPI_REAL, send_rank, itag++);

        location_index_.Reallocate(Nt_);
        data_.Reallocate(Nt_);
        forecast_.Reallocate(Nt_);
        for (int t = 0; t < Nt_; t++)
          {
            MPI::COMM_WORLD.Recv(&isize, 1, MPI::INT, send_rank, itag++);

            if (isize > 0)
              {
                if (isize < Nl_)
                  {
                    location_index_(t).Reallocate(isize);
                    MPI::COMM_WORLD.Recv(location_index_(t).GetData(), isize, MPI::INT, send_rank, itag++);
                  }
                data_(t).Reallocate(isize);
                forecast_(t).Reallocate(1);
                MPI::COMM_WORLD.Recv(data_(t).GetData(), isize, GPEC_MPI_REAL, send_rank, itag++);
                MPI::COMM_WORLD.Recv(forecast_(t).GetData(), 1, MPI::INT, send_rank, itag++);
              }
          }
      }
    else if (this->rank_ == send_rank)
      {
        isize = int(type_.size());
        MPI::COMM_WORLD.Send(&isize, 1, MPI::INT, recv_rank, itag++);
        MPI::COMM_WORLD.Send(type_.data(), isize, MPI::CHAR, recv_rank, itag++);

        isize = int(name_.size());
        MPI::COMM_WORLD.Send(&isize, 1, MPI::INT, recv_rank, itag++);
        MPI::COMM_WORLD.Send(name_.data(), isize, MPI::CHAR, recv_rank, itag++);

        isize = int(species_.size());
        MPI::COMM_WORLD.Send(&isize, 1, MPI::INT, recv_rank, itag++);
        MPI::COMM_WORLD.Send(species_.data(), isize, MPI::CHAR, recv_rank, itag++);

        isize = int(comment_.size());
        MPI::COMM_WORLD.Send(&isize, 1, MPI::INT, recv_rank, itag++);
        MPI::COMM_WORLD.Send(comment_.data(), isize, MPI::CHAR, recv_rank, itag++);

        isize = get_int_from_date(date_begin_);
        MPI::COMM_WORLD.Send(&isize, 1, MPI::INT, recv_rank, itag++);
        isize = get_int_from_date(date_end_);
        MPI::COMM_WORLD.Send(&isize, 1, MPI::INT, recv_rank, itag++);

        MPI::COMM_WORLD.Send(&Nt_day_, 1, MPI::INT, recv_rank, itag++);
        MPI::COMM_WORLD.Send(&Nt_, 1, MPI::INT, recv_rank, itag++);
        MPI::COMM_WORLD.Send(&Nd_, 1, MPI::INT, recv_rank, itag++);
        MPI::COMM_WORLD.Send(&delta_t_, 1, GPEC_MPI_REAL, recv_rank, itag++);
        MPI::COMM_WORLD.Send(&Nl_, 1, MPI::INT, recv_rank, itag++);

        MPI::COMM_WORLD.Send(&value_min_, 1, GPEC_MPI_REAL, recv_rank, itag++);
        MPI::COMM_WORLD.Send(&value_max_, 1, GPEC_MPI_REAL, recv_rank, itag++);

        for (int t = 0; t < Nt_; t++)
          {
            isize = data_(t).GetSize();
            MPI::COMM_WORLD.Send(&isize, 1, MPI::INT, recv_rank, itag++);
            if (isize > 0)
              {
                if (isize < Nl_)
                  MPI::COMM_WORLD.Send(location_index_(t).GetData(), isize, MPI::INT, recv_rank, itag++);
                MPI::COMM_WORLD.Send(data_(t).GetData(), isize, GPEC_MPI_REAL, recv_rank, itag++);
                MPI::COMM_WORLD.Send(forecast_(t).GetData(), 1, MPI::INT, recv_rank, itag++);
              }
          }
      }
  }


  template<class T>
  void ClassContainerBase<T>::MPI_Gather(const int &recv_rank)
  {
    for (int t = 0; t < Nt_; t++)
      {
        Vector1I rcounts(this->Nrank_);
        rcounts.Zero();

        rcounts(this->rank_) = data_(t).GetSize();
        MPI::COMM_WORLD.Gather(rcounts.GetData() + this->rank_, 1, MPI::INT,
                               rcounts.GetData(), 1, MPI::INT, recv_rank);

        int rcount(0);
        Vector1I displs(this->Nrank_);
        displs.Zero();
        if (this->rank_ == recv_rank)
          {
            rcount = rcounts(0);
            for (int irank = 1; irank < this->Nrank_; irank++)
              {
                displs(irank) = rcount;
                rcount += rcounts(irank);
              }
          }

        MPI::COMM_WORLD.Bcast(&rcount, 1, MPI::INT, recv_rank);
        if (rcount == 0)
          continue;

        Vector1I recv_buf_int;
        Vector1T recv_buf_real;
        if (this->rank_ == recv_rank)
          {
            recv_buf_int.Reallocate(rcount);
            recv_buf_real.Reallocate(rcount);
          }

        int scount = data_(t).GetSize();
        MPI::COMM_WORLD.Gatherv(location_index_(t).GetData(), scount, MPI::INT,
                                recv_buf_int.GetData(), rcounts.GetData(),
                                displs.GetData(), MPI::INT, recv_rank);
        MPI::COMM_WORLD.Gatherv(data_(t).GetData(), scount, GPEC_MPI_REAL,
                                recv_buf_real.GetData(), rcounts.GetData(),
                                displs.GetData(), GPEC_MPI_REAL, recv_rank);

        if (this->rank_ == recv_rank)
          {
            location_index_(t) = recv_buf_int;
            data_(t) = recv_buf_real;

            for (int l = 1; l < location_index_(t).GetSize(); l++)
              if (location_index_(t)(l - 1) > location_index_(t)(l))
                {
                  Sort(location_index_(t), data_(t));
                  break;
                }
            
            if (rcount == Nl_)
              location_index_(t).Clear();
          }
      }
  }


  template<class T>
  void ClassContainerBase<T>::MPI_Scatter(const int &send_rank)
  {
    this->set_location_per_rank(this->Nl_);

    int idate_begin, idate_end;
    if (this->rank_ == send_rank)
      {
        idate_begin = get_int_from_date(date_begin_);
        idate_end = get_int_from_date(date_end_);
      }

    MPI::COMM_WORLD.Bcast(&idate_begin, 1, MPI::INT, send_rank);
    MPI::COMM_WORLD.Bcast(&idate_end, 1, MPI::INT, send_rank);

    if (this->rank_ != send_rank)
      {
        date_begin_ = get_date_from_int(idate_begin);
        date_end_ = get_date_from_int(idate_end);
        set_Nt();
        set_Nd();
        location_index_.Reallocate(Nt_);
        data_.Reallocate(Nt_);
        forecast_.Reallocate(Nt_);
      }

    for (int t = 0; t < Nt_; t++)
      {
        int Nl = data_(t).GetSize();
        MPI::COMM_WORLD.Bcast(&Nl, 1, MPI::INT, send_rank);

        if (Nl == 0)
          continue;

        int rcount;
        Vector1I scounts(this->Nrank_), displs(this->Nrank_);
        if (Nl < this->Nl_)
          {
            if (this->rank_ == send_rank)
              {
                int irank(0), count(0);
                int loc = this->Nlocation_rank_(0);
                for (int l = 0; l < location_index_(t).GetSize(); l++)
                  {
                    count++;
                    if (location_index_(t)(l) >= loc)
                      {
                        scounts(irank) = count;
                        displs(irank) = l;  
                        count = 0;
                        irank++;
                        loc = this->location_offset_rank_(irank) + this->Nlocation_rank_(irank);
                      }
                  }
              }
          }
        else
          {
            scounts = this->Nlocation_rank_;
            displs = this->location_offset_rank_;
          }

        MPI::COMM_WORLD.Bcast(scounts.GetData(), this->Nrank_, MPI::INT, send_rank);
        rcount = scounts(this->rank_);

        if (this->rank_ != send_rank)
          {

            data_(t).Reallocate(rcount);
            forecast_(t).Reallocate(1);
          }

        if (Nl < this->Nl_)
          {
            Vector1I send_buf;
            if (this->rank_ == send_rank)
              send_buf = location_index_(t);
            location_index_(t).Reallocate(rcount);
            MPI::COMM_WORLD.Scatterv(send_buf.GetData(), scounts.GetData(),
                                     displs.GetData(), MPI::INT, location_index_(t).GetData(),
                                     rcount, MPI::INT, send_rank);
          }
        else
          location_index_(t) = this->location_rank_;

        Vector1T send_buf;
        if (this->rank_ == send_rank)
          send_buf = data_(t);
        data_(t).Reallocate(rcount);
        MPI::COMM_WORLD.Scatterv(send_buf.GetData(), scounts.GetData(), displs.GetData(),
                                 GPEC_MPI_REAL, data_(t).GetData(),
                                 rcount, GPEC_MPI_REAL, send_rank);

        forecast_(t).Resize(1);
        MPI::COMM_WORLD.Bcast(forecast_(t).GetData(), 1, MPI::INT, send_rank);
      }
  }
#endif

}

#define GPEC_FILE_CLASS_CONTAINER_BASE_CXX
#endif