ClassEnsembleBase.cxx

/ensemble/ClassEnsembleBase.cxx

https://gitlab.com/ineris/GPEC
C++ | 934 lines | 711 code | 186 blank | 37 comment | 142 complexity | 0f20400f36ba362b5638dc06d6b973cf MD5 | raw file
Possible License(s): GPL-3.0

// 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_BASE_CXX

#include "ClassEnsembleBase.hxx"

namespace GPEC
{
  // Construct.
  template<class C>
  void ClassEnsembleBase<C>::construct()
  {
    // Read configuration.
    Ops::Ops ops(this->configuration_file_);

    // Local default ensemble directory.
    if (this->container_directory_[0] != '-')
      {
        if (this->container_directory_[0] != '/')
          this->container_directory_ = ClassGPEC::ensemble_directory_
            + "/" + this->container_directory_;

#ifdef GPEC_WITH_LOGGER
        *(ClassGPEC::log_ptr_) << Fblue() << Info(3) << Reset()
                               << "Overwrites container directory with \""
                               << this->container_directory_ << "\"." << endl;
#endif
      }

    ops.SetPrefix("ensemble.");

    if (ops.Exists("parameter"))
      {
        vector<string> str_vector = ops.GetEntryList("parameter");
        for (vector<string>::iterator it = str_vector.begin(); it != str_vector.end(); it++)
          parameter_[*it] = ops.Get<real_container>("parameter." + *it);

#ifdef GPEC_WITH_LOGGER
        *(ClassGPEC::log_ptr_) << Fmagenta() << Info(3) << Reset()
                               << "Load global parameters :" << parameter_ << endl;
#endif
      }

    if (ops.Exists(this->name_))
      {
        ops.SetPrefix("ensemble." + this->name_ + ".");

        if (ops.Exists("parameter"))
          {
            vector<string> str_vector = ops.GetEntryList("parameter");
            for (vector<string>::iterator it = str_vector.begin(); it != str_vector.end(); it++)
              parameter_[*it] = ops.Get<real_container>("parameter." + *it);

#ifdef GPEC_WITH_LOGGER
            *(ClassGPEC::log_ptr_) << Fmagenta() << Info(3) << Reset()
                                   << "Load local parameters :" << parameter_ << endl;
#endif
          }
      }
#ifdef GPEC_WITH_LOGGER
    else
      *(this->log_ptr_) << Bred() << Warning() << Reset() << "Nothing known about ensemble \"" << this->name_
                        << "\" in ensemble section of configuration file \"" << ClassGPEC::configuration_file_
                        << "\"." << endl;
#endif

    ops.Close();

    if (this->container_directory_[0] != '-')
      {
        make_directory(this->container_directory_);

#ifdef GPEC_WITH_LOGGER
        *(ClassGPEC::log_ptr_) << Fyellow() << User(0, "EXPLANATION") << Reset()
                               << "Create container directory (\"" << this->container_directory_
                               << "\") for this ensemble if it does not exists (never overwritten)."
                               << " If you do not want it at all, set it to \"-\"." << endl;
#endif
      }
  }


  // Constructors.
  template<class C>
  ClassEnsembleBase<C>::ClassEnsembleBase(const string &name)
    : ClassPolyContainer<C>(name, "ensemble"), t_begin_(0), t_end_(0), algorithm_("")
  {
#ifdef GPEC_WITH_LOGGER
    *(ClassGPEC::log_ptr_) << Bcyan() << Fblack() << Info() << Reset() << "Instantiate base ensemble :" << endl;
#endif

    construct();

    return;
  }


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


  // Clear.
  template<class C>
  void ClassEnsembleBase<C>::Clear()
  {
    ClassPolyContainer<C>::Clear();

    parameter_.clear();

    ClearWeight();
  }


  // Add one container.
  template<class C>
  void ClassEnsembleBase<C>::AddSimulation(const C &simulation, string container_file)
  {
    if (container_file.empty())
      {
        ostringstream sout;
        sout <<  this->container_directory_ << "/" << simulation.species_
             << "_D" << showpos << int(rint(simulation.GetForecastDayMean()))
             << "_" << simulation.GetName() << "_" << simulation.type_
             << C::GetFileExtension();
        //   << "_" << algorithm.Name() << "_" << noshowpos << this->Ncontainer_  << C::GetFileExtension();
        container_file = sout.str();
      }

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


  // Get methods.
  template<class C>
  string ClassEnsembleBase<C>::GetAlgorithm() const
  {
    return algorithm_;
  }


  template<class C>
  int ClassEnsembleBase<C>::GetWeightTimeBegin() const
  {
    return t_begin_;
  }


  template<class C>
  int ClassEnsembleBase<C>::GetWeightTimeEnd() const
  {
    return t_end_;
  }


  template<class C>
  real ClassEnsembleBase<C>::GetParameter(const string &name) const
  {
    typename map<string, real_container>::const_iterator it = parameter_.begin();

    while (it != parameter_.end())
      if (it->first == name)
        break;
      else
        it++;

    return it->second;
  }


  template<class C>
  void ClassEnsembleBase<C>::CollectParameterList(vector<string> &list) const
  {
    list.clear();
    for (typename map<string, real_container>::const_iterator it = parameter_.begin();
         it != parameter_.end(); it++)
      list.push_back(it->first);
  }


  template<class C>
  void ClassEnsembleBase<C>::CollectWeight(Vector1I &shape, Vector<real> &weight) const
  {
    shape.Reallocate(3);
    shape(0) = t_end_ - t_begin_;
    shape(1) = weight_(t_begin_).GetSize();
    shape(2) = this->Ncontainer_;

    int Nsize = shape(0) * shape(1) * shape(2);

    if (Nsize == 0)
      return;

    weight.Reallocate(Nsize);
    weight.Fill(this->bad_value_);

    int i(0);
    for (int t = t_begin_; t < t_end_; t++)
      if (weight_(t).GetSize() > 0)
        for (int l = 0; l < weight_(t).GetSize(); l++)
          {
            if (weight_(t)(l).GetSize() > 0)
              for (int j = 0; j < this->container_index_(t).GetSize(); j++)
                weight(i + this->container_index_(t)(j)) = weight_(t)(l)(j);
            i += this->Ncontainer_;
          }
  }


  template<class C>
  void ClassEnsembleBase<C>::CollectWeightTime(const int &t, Vector1I &shape, Vector<real> &weight) const
  {
    shape.Reallocate(2);
    shape(0) = weight_(t).GetSize();
    shape(1) = this->Ncontainer_;

    int Nsize = shape(0) * shape(1);

    if (Nsize == 0)
      return;

    weight.Reallocate(Nsize);
    weight.Fill(this->bad_value_);

    int i(0);
    for (int l = 0; l < weight_(t).GetSize(); l++)
      {
        if (weight_(t)(l).GetSize() > 0)
          for (int j = 0; j < this->container_index_(t).GetSize(); j++)
            weight(i + this->container_index_(t)(j)) = weight_(t)(l)(j);
        i += this->Ncontainer_;
      }
  }


  template<class C>
  void ClassEnsembleBase<C>::CollectWeightLocation(const int &l, Vector1I &shape, Vector<real> &weight) const
  {
    check_index(0, weight_location_.GetSize(), l);

    shape.Reallocate(2);
    shape(0) = t_end_ - t_begin_;
    shape(1) = this->Ncontainer_;

    int Nsize = shape(0) * shape(1);

    if (Nsize == 0)
      return;

    weight.Reallocate(Nsize);
    weight.Fill(this->bad_value_);

    int i(0);
    for (int t = t_begin_; t < t_end_; t++)
      if (weight_(t).GetSize() > 0)
        {
          if (weight_(t)(l).GetSize() > 0)
            for (int j = 0; j < this->container_index_(t).GetSize(); j++)
              weight(i + this->container_index_(t)(j)) = weight_(t)(l)(j);

          i += this->Ncontainer_;
        }
  }


  template<class C>
  void ClassEnsembleBase<C>::CollectLocationWeight(const int &j, Vector<int> &location) const
  {
    check_index(0, weight_location_.GetSize(), j);
    location = weight_location_(j);
  }


  // Set methods.
  template<class C>
  void ClassEnsembleBase<C>::SetParameter(const string &name, const real &value)
  {
    parameter_[name] = value;
  }


  template<class C>
  void ClassEnsembleBase<C>::ReloadParameter()
  {
    // Read configuration.
    Ops::Ops ops(this->configuration_file_);

    ops.SetPrefix("ensemble.");
    if (ops.Exists("parameter"))
      {
        vector<string> str_vector = ops.GetEntryList("parameter");
        for (vector<string>::iterator it = str_vector.begin(); it != str_vector.end(); it++)
          parameter_[*it] = ops.Get<real_container>("parameter." + *it);
      }

    if (ops.Exists(this->name_))
      {
        ops.SetPrefix("ensemble." + this->name_ + ".");

        if (ops.Exists("parameter"))
          {
            vector<string> str_vector = ops.GetEntryList("parameter");
            for (vector<string>::iterator it = str_vector.begin(); it != str_vector.end(); it++)
              parameter_[*it] = ops.Get<real_container>("parameter." + *it);
          }
      }

    ops.Close();
  }


  template<class C>
  void ClassEnsembleBase<C>::SetLocationWeightRaw(const int &j, const Vector<int> &location)
  {
    check_index(0, weight_location_.GetSize(), j);
    weight_location_(j) = location;
  }


  template<class C>
  void ClassEnsembleBase<C>::SetDirectory(const string directory)
  {
    ClassPolyContainer<C>::SetDirectory(directory, "ensemble");

    if (this->container_directory_[0] != '/')
      this->container_directory_ = ClassGPEC::ensemble_directory_
        + "/" + this->container_directory_;

    make_directory(this->container_directory_);

#ifdef GPEC_WITH_LOGGER
    *(ClassGPEC::log_ptr_) << Fmagenta() << Info(1) << Reset()
                           << "Set container directory to \"" << this->container_directory_ << "\"." << endl;
#endif
  }


  // Read and write weights.
  template<class C>
  void ClassEnsembleBase<C>::ReadWeight(const string &algorithm, string weight_file, const int rank)
  {
    if (this->rank_ == rank)
      {
        if (weight_.GetSize() > 0)
          throw GPEC::Error("Weights seems already computed, clear them before.");

        algorithm_ = algorithm;

        if (weight_file.empty())
          weight_file = this->container_directory_ + "/weight_" + this->name_ + "_" + algorithm_ + "_"
            + this->container_(0).GetSpecies() + "_" + this->container_(0).GetType() + ".bin";
        else if (weight_file[0] != '/')
          weight_file = this->container_directory_ + "/" + weight_file;

        check_file(weight_file);
        string config_file = get_config_file(weight_file);

        // Read configuration.
        Ops::Ops ops(config_file);

        t_begin_ = ops.Get<int>("t_begin");
        t_end_ = ops.Get<int>("t_end");
        int Nt = ops.Get<int>("Nt");
        int Nl = ops.Get<int>("Nlocation");
        string algorithm = ops.Get<string>("algorithm");
        if (algorithm != algorithm_)
          throw GPEC::Error("Weights in file \"" + weight_file +
                            "\" seem computed with algorithm \"" + algorithm + "\".");
        ops.Close();

        weight_location_.Reallocate(Nl);
        weight_.Reallocate(Nt);

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

        for (int j = 0; j < weight_location_.GetSize(); j++)
          {
            fin.read(reinterpret_cast<char*>(&Nl), sizeof(int));
            weight_location_(j).Reallocate(Nl);
            fin.read(reinterpret_cast<char*>(weight_location_(j).GetData()), Nl * sizeof(int));
          }

        for (int t = t_begin_; t < t_end_; t++)
          {
            weight_(t).Reallocate(weight_location_.GetSize());
            for (int j = 0; j < weight_location_.GetSize(); j++)
              {
                int Ns(0);
                fin.read(reinterpret_cast<char*>(&Ns), sizeof(int));
                weight_(t)(j).Reallocate(Ns);
                fin.read(reinterpret_cast<char*>(weight_(t)(j).GetData()), Ns * sizeof(real_container));
              }
          }

        fin.close();
      }
  }


  template<class C>
  void ClassEnsembleBase<C>::WriteWeight(string weight_file, const int rank) const
  {
    if (this->rank_ == rank)
      {
        if (weight_.GetSize() == 0)
          throw GPEC::Error("Weights not yet computed.");

        if (weight_file.empty())
          weight_file = this->container_directory_ + "/weight_" + this->name_ + "_" + algorithm_ + "_"
            + this->container_(0).GetSpecies() + "_" + this->container_(0).GetType() + ".bin";
        else if (weight_file[0] != '/')
          weight_file = this->container_directory_ + "/" + weight_file;

        string config_file = get_config_file(weight_file);

#ifdef GPEC_WITH_LOGGER
        *(this->log_ptr_) << Fyellow() << Info() << Reset() << "Write weights of ensemble \""
                      << this->name_ << "\" to binary file \"" << weight_file
                      << "\" with lua configuration file \"" << config_file << "\"." << endl;
#endif

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

        ops.Push("name", this->name_);
        ops.Push("Ncontainer", this->Ncontainer_);
        ops.Push("t_begin", t_begin_);
        ops.Push("t_end", t_end_);
        ops.Push("Nt", weight_.GetSize());
        ops.Push("Nlocation", weight_location_.GetSize());
        ops.Push("algorithm", algorithm_);
        typename map<string, real_container>::const_iterator it;
        for (it = parameter_.begin(); it != parameter_.end(); it++)
          ops.Push("parameter." + it->first, it->second);

        ops.WriteLuaDefinition(config_file);
        ops.Close();

        ofstream fout(weight_file.c_str(), ofstream::binary);

        for (int j = 0; j < weight_location_.GetSize(); j++)
          {
            int Nl = weight_location_(j).GetSize();
            fout.write(reinterpret_cast<char*>(&Nl), sizeof(int));
            fout.write(reinterpret_cast<char*>(weight_location_(j).GetData()), Nl * sizeof(int));
          }

        for (int t = t_begin_; t < t_end_; t++)
          for (int j = 0; j < weight_(t).GetSize(); j++)
            {
              int Ns = weight_(t)(j).GetSize();
              fout.write(reinterpret_cast<char*>(&Ns), sizeof(int));
              fout.write(reinterpret_cast<char*>(weight_(t)(j).GetData()), Ns * sizeof(real_container));
            }

        fout.close();
      }
  }


#ifdef GPEC_WITH_MPI
  template<class C>
  void ClassEnsembleBase<C>::MPI_GatherWeightLocation(const int recv_rank)
  {
    if (this->rank_ == recv_rank)
      {
        for (int irank = 0; irank < this->Nrank_; irank++)
          if (irank != recv_rank)
            {
              int nbuf;
              MPI::COMM_WORLD.Recv(&nbuf, 1, MPI::INT, irank, 0);

              Vector1I ibuf(nbuf);
              MPI::COMM_WORLD.Recv(ibuf.GetData(), nbuf, MPI::INT, irank, 1);

              int n = weight_location_.GetSize();
              weight_location_.Resize(ibuf(0) + n);

              int i(1);
              for (int j = 0; j < ibuf(0); j++)
                {
                  Vector1I tmp(ibuf(i++));
                  for (int k = 0; k < tmp.GetSize(); k++)
                    tmp(k) = ibuf(i++);
                  weight_location_(j + n) = tmp;
                }

              MPI::COMM_WORLD.Recv(&nbuf, 1, MPI::INT, irank, 2);

              Vector1T rbuf(nbuf);
              MPI::COMM_WORLD.Recv(rbuf.GetData(), nbuf, GPEC_MPI_REAL, irank, 3);

              for (int t = t_begin_; t < t_end_; t++)
                weight_(t).Resize(ibuf(0) + n);

              int k(0);
              for (int j = 0; j < ibuf(0); j++)
                for (int t = t_begin_; t < t_end_; t++)
                  {
                    if (ibuf(i) > 0)
                      {
                        Vector1T tmp(ibuf(i));
                        for (int s = 0; s < tmp.GetSize(); s++)
                          tmp(s) = rbuf(k++);
                        weight_(t)(j + n) = tmp;
                      }
                    i++;
                  }
            }
      }
    else
      {
        int nbuf(1);
        for (int j = 0; j < weight_location_.GetSize(); j++)
          nbuf += (1 + weight_location_(j).GetSize());
        nbuf += (t_end_ - t_begin_) * weight_location_.GetSize();

        Vector1I ibuf(nbuf);
        int i(0);
        ibuf(i++) = weight_location_.GetSize();
        for (int j = 0; j < weight_location_.GetSize(); j++)
          {
            ibuf(i++) = weight_location_(j).GetSize();
            for (int k = 0; k < weight_location_(j).GetSize(); k++)
              ibuf(i++) = weight_location_(j)(k);
          }

        for (int j = 0; j < weight_location_.GetSize(); j++)
          for (int t = t_begin_; t < t_end_; t++)
            ibuf(i++) = weight_(t)(j).GetSize();

        MPI::COMM_WORLD.Send(&nbuf, 1, MPI::INT, recv_rank, 0);
        MPI::COMM_WORLD.Send(ibuf.GetData(), nbuf, MPI::INT, recv_rank, 1);

        nbuf = 0;
        for (int j = 0; j < weight_location_.GetSize(); j++)
          for (int t = t_begin_; t < t_end_; t++)
            nbuf += weight_(t)(j).GetSize();

        Vector1T rbuf(nbuf);
        i = 0;
        for (int j = 0; j < weight_location_.GetSize(); j++)
          for (int t = t_begin_; t < t_end_; t++)
            for (int s = 0; s < weight_(t)(j).GetSize(); s++)
              rbuf(i++) = weight_(t)(j)(s);

        MPI::COMM_WORLD.Send(&nbuf, 1, MPI::INT, recv_rank, 2);
        MPI::COMM_WORLD.Send(rbuf.GetData(), nbuf, GPEC_MPI_REAL, recv_rank, 3);
      }

    MPI::COMM_WORLD.Barrier();
  }


  template<class C>
  void ClassEnsembleBase<C>::MPI_GatherWeightTime(const int recv_rank)
  {
    if (this->rank_ == recv_rank)
      {
        for (int irank = 0; irank < this->Nrank_; irank++)
          if (irank != recv_rank)
            {
              int t_begin, t_end;
              MPI::COMM_WORLD.Recv(&t_begin, 1, MPI::INT, irank, 0);
              MPI::COMM_WORLD.Recv(&t_end, 1, MPI::INT, irank, 1);

              int n = weight_location_.GetSize();
              for (int t = t_begin; t < t_end; t++)
                weight_(t).Reallocate(n);

              int nbuf;
              MPI::COMM_WORLD.Recv(&nbuf, 1, MPI::INT, irank, 2);

              Vector1I ibuf(nbuf);
              MPI::COMM_WORLD.Recv(ibuf.GetData(), nbuf, MPI::INT, irank, 3);

              nbuf = Norm1(ibuf);
              Vector1T rbuf(nbuf);
              MPI::COMM_WORLD.Recv(rbuf.GetData(), nbuf, GPEC_MPI_REAL, irank, 5);

              int i(0), k(0);
              for (int j = 0; j < n; j++)
                for (int t = t_begin; t < t_end; t++)
                  {
                    if (ibuf(i) > 0)
                      {
                        Vector1T tmp(ibuf(i));
                        for (int s = 0; s < tmp.GetSize(); s++)
                          tmp(s) = rbuf(k++);
                        weight_(t)(j) = tmp;
                      }
                    i++;
                  }

              t_begin_ = min<int>(t_begin_, t_begin);
              t_end_ = max<int>(t_end_, t_end);
            }
      }
    else
      {
        MPI::COMM_WORLD.Send(&t_begin_, 1, MPI::INT, recv_rank, 0);
        MPI::COMM_WORLD.Send(&t_end_, 1, MPI::INT, recv_rank, 1);

        int nbuf = (t_end_ - t_begin_) * weight_location_.GetSize();

        Vector1I ibuf(nbuf);
        int i(0);
        for (int j = 0; j < weight_location_.GetSize(); j++)
          for (int t = t_begin_; t < t_end_; t++)
            ibuf(i++) = weight_(t)(j).GetSize();

        MPI::COMM_WORLD.Send(&nbuf, 1, MPI::INT, recv_rank, 2);
        MPI::COMM_WORLD.Send(ibuf.GetData(), nbuf, MPI::INT, recv_rank, 3);

        nbuf = Norm1(ibuf);
        Vector1T rbuf(nbuf);
        i = 0;
        for (int j = 0; j < weight_location_.GetSize(); j++)
          for (int t = t_begin_; t < t_end_; t++)
            for (int s = 0; s < weight_(t)(j).GetSize(); s++)
              rbuf(i++) = weight_(t)(j)(s);

        MPI::COMM_WORLD.Send(rbuf.GetData(), nbuf, GPEC_MPI_REAL, recv_rank, 5);
      }

    MPI::COMM_WORLD.Barrier();
  }
#endif


  // Clear weights.
  template<class C>
  void ClassEnsembleBase<C>::ClearWeight()
  {
    algorithm_ = "";
    t_begin_ = 0;
    t_end_ = 0;
    weight_location_.Clear();
    weight_.Clear();
  }


  template<class C>
  void ClassEnsembleBase<C>::Aggregate(C &aggregate, string name, string comment, const bool mask_negative_value) const
  {
    if (weight_.GetSize() == 0)
      throw GPEC::Error("Weights not yet computed.");

    int Nt = t_end_ - t_begin_;
    int Nlocation_local(0);
    for (int j = 0; j < weight_location_.GetSize(); j++) 
      Nlocation_local += weight_location_(j).GetSize();
    int Nsize = Nt * Nlocation_local;

    Vector1I time(Nsize), location(Nsize);
    Vector1T data(Nsize);
    data.Zero();

    int i(0);
    for (int j = 0; j < weight_location_.GetSize(); j++)
      {
        Vector2I simulation_location(this->container_(0).GetNt());
        for (int t = t_begin_; t < t_end_; t++)
          simulation_location(t) = weight_location_(j);

        Vector3T simulation_data;
        this->load_container(weight_location_(j), simulation_location, simulation_data);

        for (int l = 0; l < weight_location_(j).GetSize(); l++)
          for (int t = t_begin_; t < t_end_; t++)
            if (weight_(t).GetSize() > 0)
              if (weight_(t)(j).GetSize() > 0)
                {
                  time(i) = t - t_begin_;
                  location(i) = weight_location_(j)(l);
                  for (int s = 0; s < weight_(t)(j).GetSize(); s++)
                    data(i) += simulation_data(t)(l)(s) * weight_(t)(j)(s);
                  i++;
                }
      }

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

    aggregate.Set(this->container_(0));

    Vector1I forecast;
    this->CollectForecastDayMean(forecast, t_begin_, t_end_);

    aggregate.template SetData<real_container>(this->container_(0).GetDate(t_begin_),
                                              this->container_(0).GetDate(t_end_),
                                              time,
                                              location,
                                              forecast,
                                              data,
                                              mask_negative_value);
    if (name.empty())
      name = this->name_ + algorithm_;

    if (comment.empty())
      comment = "aggregated ensemble \"" + this->name_ + "\" with algorithm \"" + algorithm_ + "\"";

    aggregate.SetName(name);
    aggregate.AddComment(comment);

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


  // Best model.
  template<class C>
  void ClassEnsembleBase<C>::BestModel(C &best_model)
  {
    best_model.Clear();

    string name = this->name_ + "BestModel";
    string comment = "Best model of ensemble \"" + this->name_ + "\"";
    Aggregate(best_model, name, comment, true);
  }


  template<class C>
  void ClassEnsembleBase<C>::BestModelIndex(C &best_model_index)
  {
    int Nt = t_end_ - t_begin_;
    int Nlocation_local(0);
    for (int j = 0; j < weight_location_.GetSize(); j++) 
      Nlocation_local += weight_location_(j).GetSize();
    int Nsize = Nt * Nlocation_local;

    Vector1I time(Nsize), location(Nsize);
    Vector1T data(Nsize);
    data.Zero();

    int i(0);
    for (int j = 0; j < weight_location_.GetSize(); j++)
      {
        Vector2I simulation_location(this->container_(0).GetNt());
        for (int t = t_begin_; t < t_end_; t++)
          simulation_location(t) = weight_location_(j);

        Vector3T simulation_data;
        this->load_container(weight_location_(j), simulation_location, simulation_data);

        for (int l = 0; l < weight_location_(j).GetSize(); l++)
          for (int t = t_begin_; t < t_end_; t++)
            if (weight_(t).GetSize() > 0)
              if (weight_(t)(j).GetSize() > 0)
                {
                  time(i) = t - t_begin_;
                  location(i) = weight_location_(j)(l);
                  int s;
                  for (s = 0; s < weight_(t)(j).GetSize(); s++)
                    if (weight_(t)(j)(s) > real_container(0))
                      break;
                  data(i) = real_container(this->container_index_(t)(s));
                  i++;
                }
      }

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

    best_model_index.Set(this->container_(0));

    Vector1I forecast;
    this->CollectForecastDayMean(forecast, t_begin_, t_end_);

    best_model_index.template SetData<real_container>(this->container_(0).GetDate(t_begin_),
                                                     this->container_(0).GetDate(t_end_),
                                                     time,
                                                     location,
                                                     forecast,
                                                     data);

    best_model_index.SetName(this->name_ + "BestModelIndex");
    best_model_index.AddComment("Best model index of ensemble \"" + this->name_ + "\""); 
  }


  // Median.
  template<class C>
  void ClassEnsembleBase<C>::Median(C &median, int t_begin, int t_end)
  {
    if (! this->is_loaded_)
      throw GPEC::Error("Ensemble not loaded yet.");

    algorithm_ = "Median";

    median.Clear();
    int Nt = this->container_(0).GetNt();

    t_begin_ = t_begin;
    t_end_ = t_end;
    if (t_end_ == 0)
      t_end_ = Nt;

    weight_.Reallocate(Nt);
    for (int t = t_begin_; t < t_end_; t++)
      {
        weight_(t).Reallocate(1);
        int Nsim = this->container_index_(t).GetSize();
        if (Nsim > 0)
          {
            weight_(t)(0).Reallocate(Nsim);
            weight_(t)(0).Fill(1 / real_container(Nsim));
          }
      }

    this->set_location_per_rank(this->container_(0).GetNlocation());
    weight_location_.Reallocate(1);
    weight_location_(0) = this->location_rank_;

    string name = this->name_ + "Median";
    string comment = "Median of ensemble \"" + this->name_ + "\"";
    Aggregate(median, name, comment, true);
  }


  // Standard deviation.
  template<class C>
  void ClassEnsembleBase<C>::StandardDeviation(const C &median, C &std) const
  {
    if (median.GetNt() == 0)
      throw GPEC::Error("Median not computed yet.");

    std.Clear();
    int Nt = median.GetNt();
    int Nlocation_local = this->location_rank_.GetSize();
    int Nsize = Nt * Nlocation_local;

    Vector1I std_time(Nsize), std_location(Nsize);
    Vector1T std_data(Nsize);
    std_data.Zero();

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

        // Get median data.
        Vector1I median_time;
        Vector1T median_data;
        median.CollectDataLocation(loc, median_time, median_data);

        Vector2I simulation_location2(this->container_(0).GetNt());
        for (int t = t_begin_; t < t_end_; t++)
          simulation_location2(t).PushBack(loc);
        Vector1I &simulation_location1 = simulation_location2(t_begin_);

        Vector3T simulation_data;
        this->load_container(simulation_location1, simulation_location2, simulation_data);

        for (int j = 0; j < median_time.GetSize(); j++)
          {
            int t = median_time(j);

            std_time(i) = t;
            std_location(i) = loc;
            for (int s = 0; s < this->container_index_(t + t_begin_).GetSize(); s++)
              {
                real_container tmp = simulation_data(t + t_begin_)(0)(s) - median_data(j);
                std_data(i) += tmp * tmp;
              }

            // If median exist for this time, there should be at least one model.
            std_data(i) = sqrt(std_data(i) / real_container(this->container_index_(t + t_begin_).GetSize()));
            i++;
          }
      }

    std_time.Resize(i);
    std_location.Resize(i);
    std_data.Resize(i);

    std.Set(median);

    Vector1I std_forecast(Nt);
    for (int t = 0; t < Nt; t++)
      std_forecast(t) = median.GetForecastDay(t);

    std.template SetData<real_container>(median.GetDateBegin(),
                                         median.GetDateEnd(),
                                         std_time,
                                         std_location,
                                         std_forecast,
                                         std_data);

    std.SetName(this->name_ + "Std");
    std.AddComment("Standard deviation of ensemble \"" + this->name_ + "\""); 

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

#define GPEC_FILE_CLASS_ENSEMBLE_BASE_CXX
#endif