/container/ClassContainerBase.cxx
C++ | 2189 lines | 1722 code | 399 blank | 68 comment | 348 complexity | f811424b32b6a6b5187de0553fdbb5ee MD5 | raw file
Possible License(s): GPL-3.0
Large files files are truncated, but you can click here to view the full file
- // 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…
Large files files are truncated, but you can click here to view the full file