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

/**
 * @file kernel/pdevs/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/pdevs/qss/Data.hpp>

#include <cmath>

namespace artis {
    namespace pdevs {
        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:
                enum inputs {
                    IN = 1
                };

                enum outputs {
                    OUT = 1
                };

                Quantifier(const std::string& name,
                        const Context<Time, Quantifier<Time>, QuantifierParameters>& context)
                        :
                        artis::pdevs::Dynamics<Time, Quantifier<Time>, QuantifierParameters>(name,
                                context)
                {
                    this->input_port({IN, "in"});
                    this->output_port({OUT, "out"});
                    this->observables({{UP,    "up"},
                                       {DOWN,  "down"},
                                       {VALUE, "value"}});

                    _adaptive = context.parameters().allow_offsets;
                    _adaptive_state = _adaptive ? POSSIBLE : 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(typename Time::type t, typename Time::type e,
                        const common::Bag<Time>& bag)
                {
                    dint(t);
                    dext(t, e, bag);
                }

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

                virtual void dext(typename Time::type t, typename Time::type e,
                        const common::Bag<Time>& bag)
                {
                    std::for_each(bag.begin(), bag.end(),
                            [this, t, e](const common::ExternalEvent<Time>& event) {
                                IntegratorData data;
                                double shifting_factor;
                                double value;
                                int cnt;

                                event.data()(data);
                                value = data.value;
                                if (_state == INIT) {
                                    init_step_number_and_offset(value);
                                    update_thresholds();
                                    _state = 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 IMPOSSIBLE:
                                            update_thresholds();
                                            break;
                                        case 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_DOWN);
                                                } else {
                                                    update_thresholds(shifting_factor,
                                                            DIRECTION_UP);
                                                }
                                                _adaptive_state = DONE;
                                            } else {
                                                update_thresholds();
                                            }
                                            break;
                                        case DONE:
                                            init_step_number_and_offset(value);
                                            _adaptive_state = POSSIBLE;
                                            update_thresholds();
                                            break;
                                        }
                                    }
                                }
                            });
                    _state = RESPONSE;
                }

                virtual typename Time::type start(typename Time::type /* time */)
                {
                    _offset = 0;
                    _state = INIT;
                    return Time::infinity;
                }

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

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

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

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

            private:
                typedef enum {
                    DIRECTION_UP, DIRECTION_DOWN
                } Direction;

                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));
                }

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

                    if (d == 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));
                    }
                }

                typedef enum vars {
                    UP, DOWN, VALUE
                } Observable;

                typedef enum {
                    INIT, IDLE, RESPONSE
                } State;

                typedef enum {
                    IMPOSSIBLE, POSSIBLE, DONE
                } AdaptiveState;

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

                State _state;
                AdaptiveState _adaptive_state;
                bool _adaptive;
                bool _zero_init_offset;

                double _offset;
                double _step_size;

                long int _step_number;
                double _up_threshold;
                double _down_threshold;

                std::deque<record_t> _archive;

                unsigned int _past_length;
            };

        }
    }
}

#endif