artis-star/src/artis-star/kernel/qss/Quantifier.hpp

/**
 * @file kernel/qss/Quantifier.hpp
 * @author The ARTIS Development Team
 * See the AUTHORS or Authors.txt file
 */

/*
 * ARTIS - the multimodeling and simulation environment
 * This file is a part of the ARTIS environment
 *
 * Copyright (C) 2013-2019 ULCO http://www.univ-littoral.fr
 *
 * This program 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.
 *
 * This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef QSS_QUANTIFIER
#define QSS_QUANTIFIER

#include <artis-star/kernel/pdevs/Dynamics.hpp>
#include <artis-star/kernel/qss/Data.hpp>

#include <cmath>

namespace artis {
namespace qss {

struct QuantifierParameters
{
  bool allow_offsets;
  bool zero_init_offset;
  double quantum;
  unsigned int archive_length;
};

template<class Time>
class Quantifier
    : public artis::pdevs::Dynamics<Time, Quantifier<Time>, QuantifierParameters>
{
public:
  struct input
  {
    enum values
    {
      IN, RESET
    };
  };

  struct output
  {
    enum values
    {
      OUT
    };
  };

  Quantifier(const std::string &name,
             const artis::pdevs::Context<Time, Quantifier<Time>, QuantifierParameters> &context)
      :
      artis::pdevs::Dynamics<Time, Quantifier<Time>, QuantifierParameters>(name, context)
  {
    DECLARE_STATES(int,
                   ((state::PHASE, &Quantifier<Time>::_phase),
                       (state::ADAPTIVE_STATE, &Quantifier<Time>::_adaptive_state)));
    DECLARE_STATES(unsigned int,
                   ((state::STEP_NUMBER, &Quantifier<Time>::_step_number)));
    DECLARE_STATES(double,
                   ((state::OFFSET, &Quantifier<Time>::_offset),
                       (state::UP_THRESHOLD, &Quantifier<Time>::_up_threshold),
                       (state::DOWN_THRESHOLD, &Quantifier<Time>::_down_threshold)));

    this->input_ports({{input::IN, "in"},
                       {input::RESET, "reset"}});
    this->output_port({output::OUT, "out"});
    this->observables({{var::UP, "up"},
                       {var::DOWN, "down"},
                       {var::VALUE, "value"}});

    _adaptive = context.parameters().allow_offsets;
    _adaptive_state = _adaptive ? adaptive_state::POSSIBLE : adaptive_state::IMPOSSIBLE;
    _zero_init_offset = context.parameters().zero_init_offset;
    _step_size = context.parameters().quantum;

    assert(_step_size > 0);

    _past_length = context.parameters().archive_length;

    assert(_past_length > 2);
  }

  virtual ~Quantifier()
  {}

  virtual void dconf(const typename Time::type &t, const typename Time::type &e,
                     const common::Bag<Time> &bag)
  {
    dint(t);
    dext(t, e, bag);
  }

  virtual void dint(const typename Time::type & /* t */)
  {
    switch (_phase) {
    case phase::INIT:break;
    case phase::IDLE:break;
    case phase::RESPONSE:_phase = phase::IDLE;
      break;
    }
  }

  virtual void dext(const typename Time::type &t, const typename Time::type &e,
                    const common::Bag<Time> &bag)
  {
    bool reset = false;

    std::for_each(bag.begin(), bag.end(),
                  [this, t, e, &reset](const common::ExternalEvent<Time> &event) {
                    if (event.on_port(input::IN)) {
                      IntegratorData data;
                      double shifting_factor;
                      double value;
                      int cnt;

                      event.data()(data);
                      value = data.value;
                      if (_phase == phase::INIT) {
                        init_step_number_and_offset(value);
                        update_thresholds();
                        _phase = phase::RESPONSE;
                      } else {
                        cnt = 0;
                        while (value >= _up_threshold or value <= _down_threshold) {
                          cnt++;
                          if (value >= _up_threshold) {
                            _step_number++;
                          } else {
                            _step_number--;
                          }
                          switch (_adaptive_state) {
                          case adaptive_state::IMPOSSIBLE:update_thresholds();
                            break;
                          case adaptive_state::POSSIBLE:
                            if (value >= _up_threshold) {
                              store_change(_step_size, t);
                            } else {
                              store_change(-_step_size, t);
                            }
                            shifting_factor = shift_quanta();

                            assert(shifting_factor >= 0
                                       and shifting_factor <= 1);

                            if (shifting_factor != 0 and shifting_factor != 1) {
                              if (value >= _up_threshold) {
                                update_thresholds(shifting_factor,
                                                  direction::DIRECTION_DOWN);
                              } else {
                                update_thresholds(shifting_factor,
                                                  direction::DIRECTION_UP);
                              }
                              _adaptive_state = adaptive_state::DONE;
                            } else {
                              update_thresholds();
                            }
                            break;
                          case adaptive_state::DONE:init_step_number_and_offset(value);
                            _adaptive_state = adaptive_state::POSSIBLE;
                            update_thresholds();
                            break;
                          }
                        }
                      }
                    } else if (event.on_port(input::RESET)) {
                      _offset = 0;
                      reset = true;
                      _archive.clear();
                    }
                  });
    if (reset) {
      _phase = phase::INIT;
    } else {
      _phase = phase::RESPONSE;
    }
  }

  virtual void start(const typename Time::type & /* time */)
  {
    _offset = 0;
    _phase = phase::INIT;
  }

  virtual typename Time::type ta(const typename Time::type & /* time */)
  {
    switch (_phase) {
    case phase::INIT:
    case phase::IDLE:return Time::infinity;
    case phase::RESPONSE:return 0.0;
    }
    return Time::infinity;
  }

  virtual common::Bag<Time> lambda(const typename Time::type & /* time */) const
  {
    common::Bag<Time> msgs;
    const QuantifierData data = {_up_threshold, _down_threshold};

    msgs.push_back(common::ExternalEvent<Time>(output::OUT, data));
    return msgs;
  }

  virtual common::Value observe(const typename Time::type & /* t */,
                                unsigned int index) const
  {
    switch (index) {
    case var::UP:return (double) _up_threshold;
    case var::DOWN:return (double) _down_threshold;
    case var::VALUE:return (double) (_up_threshold - _down_threshold);
    default:return common::Value();
    }
  }

private:
  void init_step_number_and_offset(double value)
  {
    _step_number = static_cast<long int>(std::floor(value / _step_size));
    if (_zero_init_offset) {
      _offset = 0;
    } else {
      _offset = value - static_cast<double>(_step_number) * _step_size;
    }
  }

  bool monotonous(unsigned int range)
  {
    if ((range + 1) > _archive.size()) {
      return false;
    }
    for (size_t i = 0; i < range; i++) {
      if (_archive[i].value * _archive[i + 1].value < 0) {
        return false;
      }
    }
    return true;
  }

  bool oscillating(unsigned int range)
  {
    if ((range + 1) > _archive.size()) {
      return false;
    }
    for (size_t i = _archive.size() - range; i < _archive.size() - 1; i++) {
      if (_archive[i].value * _archive[i + 1].value > 0) {
        return false;
      }
    }
    return true;
  }

  double shift_quanta()
  {
    double factor = 0;

    if (oscillating(_past_length - 1) and
        _archive.back().date - _archive.front().date != 0) {
      double acc;
      double local_estim;
      int cnt;

      acc = 0;
      cnt = 0;
      for (size_t i = 0; i < _archive.size() - 2; ++i) {
        if (0 != (_archive[i + 2].date - _archive[i].date)) {
          if ((_archive.back().value * _archive[i + 1].value) > 0) {
            local_estim =
                1 - (_archive[i + 1].date - _archive[i].date) /
                    (_archive[i + 2].date - _archive[i].date);
          } else {
            local_estim = (_archive[i + 1].date - _archive[i].date) /
                (_archive[i + 2].date - _archive[i].date);
          }
          acc += local_estim;
          cnt++;
        }
      }
      acc = acc / cnt;
      factor = acc;
      _archive.resize(0);
    }
    return factor;
  }

  void store_change(double val, const typename Time::type &time)
  {
    record_t record;

    record.date = time;
    record.value = val;
    _archive.push_back(record);
    while (_archive.size() > _past_length) {
      _archive.pop_front();
    }
  }

  void update_thresholds()
  {
    auto step_number = static_cast<double>(_step_number);

    _up_threshold = _offset + _step_size * (step_number + 1);
    _down_threshold = _offset + _step_size * (step_number - 1);
  }

  void update_thresholds(double factor)
  {
    auto step_number = static_cast<double>(_step_number);

    _up_threshold = _offset + _step_size * (step_number + (1 - factor));
    _down_threshold = _offset + _step_size * (step_number - (1 - factor));
  }

  struct direction
  {
    enum values
    {
      DIRECTION_UP, DIRECTION_DOWN
    };
  };

  void update_thresholds(double factor, const typename direction::values &d)
  {
    auto step_number = static_cast<double>(_step_number);

    if (d == direction::DIRECTION_UP) {
      _up_threshold = _offset + _step_size * (step_number + (1 - factor));
      _down_threshold = _offset + _step_size * (step_number - 1);
    } else {
      _up_threshold = _offset + _step_size * (step_number + 1);
      _down_threshold = _offset + _step_size * (step_number - (1 - factor));
    }
  }

  struct state
  {
    enum values
    {
      PHASE, ADAPTIVE_STATE, STEP_NUMBER, OFFSET, UP_THRESHOLD, DOWN_THRESHOLD
    };
  };

  struct var
  {
    enum values
    {
      UP, DOWN, VALUE
    };
  };

  struct phase
  {
    enum values
    {
      INIT, IDLE, RESPONSE
    };
  };

  struct adaptive_state
  {
    enum values
    {
      IMPOSSIBLE, POSSIBLE, DONE
    };
  };

  struct record_t
  {
    double value;
    typename Time::type date;
  };

  // parameters
  bool _adaptive;
  bool _zero_init_offset;
  unsigned int _past_length;
  double _step_size;

  // state
  int _phase;
  int _adaptive_state;

  unsigned int _step_number; // long int

  double _offset;
  double _up_threshold;
  double _down_threshold;

  std::deque<record_t> _archive;
};

}
}

#endif