Updated DEVStone for aDEVS (testing included)

.gitignore

@@ -3,6 +3,7 @@
 
 cmake-build-debug/
 build/
+bin/
 
 .DS_Store
 

devstone/adevs/CMakeLists.txt

@@ -1,29 +1,28 @@
 cmake_minimum_required(VERSION 3.10)
 project(devstone-adevs)
-
 set (CMAKE_CXX_STANDARD 17)
-set (CMAKE_CXX_COMPILER "g++")
-add_compile_options(-g)
+set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/bin)
 
-set(Boost_USE_STATIC_LIBS        ON)
-set(Boost_USE_MULTITHREADED     OFF)
+if(NOT CMAKE_BUILD_TYPE)
+  set(CMAKE_BUILD_TYPE Release)
+endif()
+set(CMAKE_CXX_FLAGS "-Wall -Wextra")
+set(CMAKE_CXX_FLAGS_DEBUG "-g")
+set(CMAKE_CXX_FLAGS_RELEASE "-O3")
 
-include_directories(${CMAKE_CURRENT_SOURCE_DIR}/../../simulators/adevs/include)
+include_directories(SYSTEM ${CMAKE_CURRENT_SOURCE_DIR}/../../simulators/adevs-3.3/include)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR})
 
-find_package(Boost COMPONENTS program_options thread unit_test_framework REQUIRED)
-
-include_directories(${Boost_INCLUDE_DIRS})
-
-
-set(CMAKE_CXX_STANDARD 17)
-
-add_executable(DEVStone src/DEVStone.cpp src/DummyAtomic.cpp src/DEVSWrapper.cpp
-        src/li.cpp src/hoMod.cpp src/hi.cpp src/hi.hpp src/ho.cpp src/seeder.cpp)
-target_link_libraries(DEVStone ${Boost_PROGRAM_OPTIONS_LIBRARY})
-
-#add_executable(DEVStoneParallel src/DEVStone.cpp src/DummyAtomic.cpp src/DEVSWrapper.cpp
-#        src/li.cpp src/hoMod.cpp src/hi.cpp src/hi.hpp src/ho.cpp src/seeder.cpp)
-
-#target_link_libraries(DEVStoneParallel ${Boost_PROGRAM_OPTIONS_LIBRARY})
-#set_target_properties(DEVStoneParallel PROPERTIES COMPILE_FLAGS "-fopenmp")
+add_executable(devstone src/devstone.cpp)
+
+# Testing stuff
+find_package(Boost COMPONENTS system filesystem unit_test_framework)
+if (Boost_FOUND)
+  add_definitions(-DBOOST_TEST_DYN_LINK)
+  enable_testing()
+  add_executable(test_devstone test/test_devstone.cpp)
+  target_link_libraries(test_devstone ${Boost_FILESYSTEM_LIBRARY} ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+  add_test(NAME test_devstone COMMAND test_devstone)
+else()
+  message(STATUS "Boost not found. You won't be able to run tests for adevs.")
+endif()

devstone/adevs/src/DEVSWrapper.cpp

@@ -1,31 +0,0 @@
-//
-// Created by Román Cárdenas Rodríguez on 09/03/2020.
-//
-
-#include "DEVSWrapper.hpp"
-
-const int DEVSWrapper::in = 0;
-const int DEVSWrapper::out = 1;
-
-DEVSWrapper::DEVSWrapper(int depthIn, int widthIn, int intDelayIn, int extDelayIn, double procTimeIn):
-        adevs::Digraph<int*>()
-{
-    depth = depthIn;
-    width = widthIn;
-    intDelay = intDelayIn;
-    extDelay = extDelayIn;
-    procTime = procTimeIn;
-
-    if (depth == 1) {
-        DummyAtomic *a = new DummyAtomic(intDelay, extDelay, procTime);
-        add(a);
-        couple(this, this->in, a, a->in);
-        couple(a, a->out, this, this->out);
-    } else {
-        for (int i = 0; i < width - 1; ++i) {
-            DummyAtomic *a = new DummyAtomic(intDelay, extDelay, procTime);
-            add(a);
-            atomics.push_back(a);
-        }
-    }
-}

devstone/adevs/src/DEVSWrapper.hpp

@@ -1,28 +0,0 @@
-//
-// Created by rcardenas on 3/9/20.
-//
-
-#ifndef DEVSTONE_ADEVS_DEVSWRAPPER_HPP
-#define DEVSTONE_ADEVS_DEVSWRAPPER_HPP
-
-#include "utils.hpp"
-#include "DummyAtomic.hpp"
-#include <vector>
-
-class DEVSWrapper: public adevs::Digraph<int*>
-{
-protected:
-    int depth, width;
-    int intDelay, extDelay;
-    double procTime;
-    std::vector<DummyAtomic*> atomics;
-public:
-    /// Model input port
-    static const int in;
-    /// Model output port
-    static const int out;
-    /// Constructor
-    DEVSWrapper(int depthIn, int widthIn, int intDelayIn, int extDelayIn, double procTimeIn);
-};
-
-#endif //DEVSTONE_ADEVS_DEVSWRAPPER_HPP

devstone/adevs/src/DEVStone.cpp

@@ -1,103 +0,0 @@
-//
-// Created by rcardenas on 3/9/20.
-//
-
-#include "DEVStone.hpp"
-
-
-DEVStone::DEVStone(const string& typeIn, int depthIn, int widthIn, int intDelayIn, int extDelayIn, double procTimeIn):
-        adevs::Digraph<int*>()
-{
-    Seeder *seeder = new Seeder(1);
-    add(seeder);
-
-    adevs::Digraph<int*> devstone;
-    if (typeIn == "LI") {
-        LI *li = new LI(depthIn, widthIn, intDelayIn, extDelayIn, procTimeIn);
-        add(li);
-        couple(seeder, seeder->out, li, li->in);
-    } else if (typeIn == "HI") {
-        HI *hi = new HI(depthIn, widthIn, intDelayIn, extDelayIn, procTimeIn);
-        add(hi);
-        couple(seeder, seeder->out, hi, hi->in);
-    } else if (typeIn == "HO") {
-        HO *ho = new HO(depthIn, widthIn, intDelayIn, extDelayIn, procTimeIn);
-        add(ho);
-        couple(seeder, seeder->out, ho, ho->in);
-    } else if (typeIn == "HOmod") {
-        HOmod *homod = new HOmod(depthIn, widthIn, intDelayIn, extDelayIn, procTimeIn);
-        add(homod);
-        couple(seeder, seeder->out, homod, homod->in1);
-        couple(seeder, seeder->out, homod, homod->in2);
-    } else {
-        throw runtime_error("Unable to determine DEVStone model type");
-    }
-}
-
-
-int main(int argc, char* argv[]) {
-    auto start = hclock::now();
-
-    // Declare the supported options.
-    po::options_description desc("Allowed options");
-    desc.add_options()
-            ("help", "produce help message")
-            ("kind", po::value<string>()->required(), "set kind of devstone: LI, HI, HO, or HOmod")
-            ("width", po::value<int>()->required(), "set width of the DEVStone: integer value")
-            ("depth", po::value<int>()->required(), "set depth of the DEVStone: integer value")
-            ("int-cycles", po::value<int>()->required(), "set the Dhrystone cycles to expend in internal transtions: integer value")
-            ("ext-cycles", po::value<int>()->required(), "set the Dhrystone cycles to expend in external transtions: integer value")
-            ("time-advance", po::value<int>()->default_value(1), "set the time expend in external transtions by the Dhrystone in miliseconds: integer value")
-            ;
-
-    po::variables_map vm;
-    try {
-        po::store(po::parse_command_line(argc, argv, desc), vm);
-        po::notify(vm);
-    } catch ( boost::program_options::required_option be ){
-        if (vm.count("help")) {
-            cout << desc << "\n";
-            return 0;
-        } else {
-            cout << be.what() << endl;
-            cout << endl;
-            cout << "for mode information run: " << argv[0] << " --help" << endl;
-            return 1;
-        }
-    }
-    string kind = vm["kind"].as<string>();
-    if (kind != "LI"  && kind != "HI" && kind != "HO" && kind != "HOmod") {
-        cout << "The kind needs to be LI, HI, HO, or HOmod and received value was: " << kind;
-        cout << endl;
-        cout << "for mode information run: " << argv[0] << " --help" << endl;
-        return 1;
-    }
-
-    int width = vm["width"].as<int>();
-    int depth = vm["depth"].as<int>();
-    int int_cycles = vm["int-cycles"].as<int>();
-    int ext_cycles = vm["ext-cycles"].as<int>();
-    int time_advance = vm["time-advance"].as<int>();
-
-    auto processed_parameters = hclock::now();
-
-    DEVStone *devstone = new DEVStone(kind, depth, width, int_cycles, ext_cycles, time_advance);
-
-    auto model_built = hclock::now();
-
-    cout << endl;
-    cout << "Model creation time: " << chrono::duration_cast<chrono::duration<double, ratio<1>>>( model_built - processed_parameters).count() << endl;
-    adevs::Simulator<IO_Type> sim(devstone);
-
-    auto model_init = hclock::now();
-
-    cout << "Engine setup time: " << chrono::duration_cast<chrono::duration<double, ratio<1>>>( model_init - model_built).count() << endl;
-    while (sim.nextEventTime() < DBL_MAX)
-    {
-        sim.execNextEvent();
-    }
-    /// Done!
-    auto finished_simulation = hclock::now();
-    cout << "Simulation time: " << chrono::duration_cast<chrono::duration<double, ratio<1>>>( finished_simulation - model_init).count() << endl;
-    return 0;
-}

devstone/adevs/src/DEVStone.hpp

@@ -1,31 +0,0 @@
-//
-// Created by Román Cárdenas Rodríguez on 09/03/2020.
-//
-
-#ifndef DEVSTONE_ADEVS_DEVSTONE_HPP
-#define DEVSTONE_ADEVS_DEVSTONE_HPP
-
-#include <boost/program_options.hpp>
-#include <iostream>
-#include <chrono>
-#include <string>
-
-#include "utils.hpp"
-#include "li.hpp"
-#include "hi.hpp"
-#include "ho.hpp"
-#include "hoMod.hpp"
-#include "seeder.hpp"
-
-using namespace std;
-namespace po=boost::program_options;
-using hclock=chrono::high_resolution_clock;
-
-
-class DEVStone: public adevs::Digraph<int*>
-{
-public:
-    DEVStone(const string& typeIn, int depthIn, int widthIn, int intDelayIn, int extDelayIn, double propTimeIn);
-};
-
-#endif //DEVSTONE_ADEVS_DEVSTONE_HPP

devstone/adevs/src/DummyAtomic.cpp

@@ -1,57 +0,0 @@
-//
-// Created by Román Cárdenas Rodríguez on 09/03/2020.
-//
-
-#include "DummyAtomic.hpp"
-#include "../../dhry/dhry_1.c"
-
-const int DummyAtomic::in = 0;
-const int DummyAtomic::out = 1;
-
-DummyAtomic::DummyAtomic(int intDelayIn, int extDelayIn, double procTimeIn): adevs::Atomic<IO_Type>() {
-    intDelay = intDelayIn;
-    extDelay = extDelayIn;
-    procTime = procTimeIn;
-    sigma = DBL_MAX;
-    t = 0.0;
-}
-
-/// External transition function
-void DummyAtomic::delta_ext(double e, const adevs::Bag<IO_Type>& x) {
-    DhryStone().dhrystoneRun(extDelay);
-    t += e;
-    if (sigma != DBL_MAX) {
-        sigma -= e;
-    } else {
-        sigma = procTime;
-    }
-}
-
-/// Internal transition function
-void DummyAtomic::delta_int() {
-    DhryStone().dhrystoneRun(intDelay);
-    t += sigma;
-    sigma = DBL_MAX;
-}
-
-/// Confluent transition function.
-void DummyAtomic::delta_conf(const adevs::Bag<IO_Type>& x) {
-    delta_int();
-    delta_ext(0.0, x);
-}
-
-/// Output function.
-void DummyAtomic::output_func(adevs::Bag<IO_Type >& y) {
-    y.insert(IO_Type(out, 0));
-}
-
-/// Time advance function.
-double DummyAtomic::ta() {
-    return sigma;
-}
-
-/// Garbage collection. No heap allocation in output_func, so do nothing.
-void DummyAtomic::gc_output(adevs::Bag<IO_Type>&) {}
-
-/// Destructor
-DummyAtomic::~DummyAtomic() = default;

devstone/adevs/src/DummyAtomic.hpp

@@ -1,45 +0,0 @@
-//
-// Created by rcardenas on 3/9/20.
-//
-
-#ifndef DEVSTONE_ADEVS_DUMMYATOMIC_HPP
-#define DEVSTONE_ADEVS_DUMMYATOMIC_HPP
-
-#include <adevs.h>
-#include <list>
-#include "utils.hpp"
-
-
-class DummyAtomic: public adevs::Atomic<IO_Type> {
-    private:
-        /// Model state variables
-        int intDelay;  // TODO include Dhrystone stuff
-        int extDelay;
-        double sigma;
-        double procTime;
-        double t;
-    public:
-        /// Model input port
-        static const int in;
-        /// Model output port
-        static const int out;
-
-        /// Constructor.  The processing time is provided as an argument.
-        DummyAtomic(int intDelayIn, int extDelayIn, double procTimeIn);
-        /// External transition function
-        void delta_ext(double e, const adevs::Bag<IO_Type>& x);
-        /// Internal transition function
-        void delta_int();
-        /// Confluent transition function.
-        void delta_conf(const adevs::Bag<IO_Type>& x);
-        /// Output function.
-        void output_func(adevs::Bag<IO_Type>& y);
-        /// Time advance function.
-        double ta();
-        /// Garbage collection.
-        void gc_output(adevs::Bag<IO_Type>&);
-        /// Destructor
-        ~DummyAtomic();
-};
-
-#endif //DEVSTONE_ADEVS_DUMMYATOMIC_HPP

devstone/adevs/src/devstone.cpp

@@ -0,0 +1,53 @@
+#include <chrono>
+#include <iostream>
+#include <string>
+#include "devstone.hpp"
+using namespace std;
+using hclock=chrono::high_resolution_clock;
+
+int main(int argc, char* argv[]) {
+    // First, we parse the arguments
+    if (argc < 4) {
+        std::cerr << "ERROR: not enough arguments" << std::endl;
+        std::cerr << "    Usage:" << std::endl;
+        std::cerr << "    > devstone MODEL_TYPE WIDTH DEPTH INTDELAY EXTDELAY" << std::endl;
+        std::cerr << "        (MODEL_TYPE must be either LI, HI, HO, or HOmod)" << std::endl;
+        std::cerr << "        (WIDTH and DEPTH must be greater than or equal to 1)" << std::endl;
+        std::cerr << "        (INTDELAY and EXTDELAY must be greater than or equal to 0 ms)" << std::endl;
+        std::cerr << "    Alternative usages:" << std::endl;
+        std::cerr << "    > main_devstone MODEL_TYPE WIDTH DEPTH DELAY" << std::endl;
+        std::cerr << "        (INTDELAY and EXTDELAY are set to DELAY ms)" << std::endl;
+        std::cerr << "    > main_devstone MODEL_TYPE WIDTH DEPTH" << std::endl;
+        std::cerr << "        (INTDELAY and EXTDELAY are set to 0 ms)" << std::endl;
+        return -1;
+    }
+    std::string type = argv[1];
+    int width = std::stoi(argv[2]);
+    int depth = std::stoi(argv[3]);
+    int intDelay = 0;
+    int extDelay = 0;
+    if (argc > 4) {
+        intDelay = std::stoi(argv[4]);
+        extDelay = (argc == 5) ? intDelay : std::stoi(argv[5]);
+    }
+    auto processed_parameters = hclock::now();
+
+    auto *devstone = new DEVStone(type, width, depth, intDelay, extDelay);
+
+    auto model_built = hclock::now();
+
+    cout << endl;
+    cout << "Model creation time: " << chrono::duration_cast<chrono::duration<double, ratio<1>>>( model_built - processed_parameters).count() << " seconds" << endl;
+    adevs::Simulator<IO_Type> sim(devstone);
+
+    auto model_init = hclock::now();
+
+    cout << "Engine setup time: " << chrono::duration_cast<chrono::duration<double, ratio<1>>>( model_init - model_built).count() << " seconds" << endl;
+    while (sim.nextEventTime() < DBL_MAX) {
+        sim.execNextEvent();
+    }
+    /// Done!
+    auto finished_simulation = hclock::now();
+    cout << "Simulation time: " << chrono::duration_cast<chrono::duration<double, ratio<1>>>( finished_simulation - model_init).count() << " seconds" << endl;
+    return 0;
+}

devstone/adevs/src/devstone.hpp

@@ -0,0 +1,47 @@
+#ifndef DEVSTONE_ADEVS_DEVSTONE_HPP
+#define DEVSTONE_ADEVS_DEVSTONE_HPP
+
+#include "utils.hpp"
+#include "li.hpp"
+#include "hi.hpp"
+#include "ho.hpp"
+#include "homod.hpp"
+#include "seeder.hpp"
+
+class DEVStone: public adevs::Digraph<int*> {
+public:
+    [[maybe_unused]] DEVStoneCoupled *model;
+    DEVStone(const string& type, int width, int depth, int intDelay, int extDelay): adevs::Digraph<int*>() {
+        auto *seeder = new Seeder();
+        add(seeder);
+
+        //adevs::Digraph<int*> devstone;
+        if (type == "LI") {
+            auto *li = new LI(width, depth, intDelay, extDelay);
+            add(li);
+            couple(seeder, Seeder::out, li, LI::in);
+            model = li;
+        } else if (type == "HI") {
+            auto *hi = new HI(width, depth, intDelay, extDelay);
+            add(hi);
+            couple(seeder, Seeder::out, hi, HI::in);
+            model = hi;
+        } else if (type == "HO") {
+            auto *ho = new HO(width, depth, intDelay, extDelay);
+            add(ho);
+            couple(seeder, Seeder::out, ho, HO::in);
+            couple(seeder, Seeder::out, ho, HOmod::in2);
+            model = ho;
+        } else if (type == "HOmod") {
+            auto *homod = new HOmod(width, depth, intDelay, extDelay);
+            add(homod);
+            couple(seeder, Seeder::out, homod, HOmod::in);
+            couple(seeder, Seeder::out, homod, HOmod::in2);
+            model = homod;
+        } else {
+            throw runtime_error("Unable to determine DEVStone model type");
+        }
+    }
+};
+
+#endif //DEVSTONE_ADEVS_DEVSTONE_HPP

devstone/adevs/src/devstone_atomic.hpp

@@ -0,0 +1,61 @@
+#ifndef DEVSTONE_ADEVS_DEVSTONEATOMIC_HPP
+#define DEVSTONE_ADEVS_DEVSTONEATOMIC_HPP
+
+#include <adevs.h>
+#include <list>
+#include "utils.hpp"
+
+class DEVStoneAtomic: public adevs::Atomic<IO_Type> {
+    private:
+        /// Model state variables
+        __attribute__ ((unused)) int intDelay, extDelay;  // TODO include Dhrystone stuff
+        double sigma;
+    public:
+        /// Model input port
+        static const int in = 0;
+        /// Model output port
+        static const int out = 1;
+        /// To count the number of events etc.
+        unsigned int nInternals, nExternals, nEvents;
+
+        /// Constructor
+        DEVStoneAtomic(int intDelay, int extDelay): adevs::Atomic<IO_Type>(), intDelay(intDelay), extDelay(extDelay), sigma(DBL_MAX) {
+            nInternals = 0;
+            nExternals = 0;
+            nEvents = 0;
+        }
+
+        /// External transition function
+        void delta_ext(__attribute__ ((unused)) double e, const adevs::Bag<IO_Type>& x) override {
+            // DhryStone().dhrystoneRun(extDelay); // TODO
+            nExternals += 1;
+            nEvents += x.size();
+            sigma = 0;
+        }
+
+        /// Internal transition function
+        void delta_int() override {
+            // DhryStone().dhrystoneRun(intDelay); // TODO
+            nInternals += 1;
+            sigma = DBL_MAX;
+        }
+
+        /// Confluent transition function.
+        void delta_conf(const adevs::Bag<IO_Type>& x) override {
+            delta_int();
+            delta_ext(0.0, x);
+        }
+        /// Output function.
+        void output_func(adevs::Bag<IO_Type>& y) override {
+            y.insert(IO_Type(out, 0));
+        }
+
+        /// Time advance function.
+        double ta() override {
+            return sigma;
+        }
+        /// Garbage collection.
+        void gc_output(adevs::Bag<IO_Type>&) override {}
+};
+
+#endif //DEVSTONE_ADEVS_DEVSTONEATOMIC_HPP

devstone/adevs/src/devstone_coupled.hpp

@@ -0,0 +1,98 @@
+#ifndef DEVSTONE_ADEVS_DEVSTONECOUPLED_HPP
+#define DEVSTONE_ADEVS_DEVSTONECOUPLED_HPP
+
+#include "utils.hpp"
+#include "devstone_atomic.hpp"
+#include <vector>
+
+class DEVStoneCoupled: public adevs::Digraph<int*> {
+protected:
+    DEVStoneCoupled* childCoupled;
+    std::vector<DEVStoneAtomic*> atomics;
+    int intDelay, extDelay;
+    std::size_t nEIC, nIC, nEOC;
+
+public:
+    /// Model input port
+    static const int in = 0;
+    /// Model output port
+    static const int out = 1;
+    /// Constructor
+    DEVStoneCoupled(int width, int depth, int intDelay, int extDelay): adevs::Digraph<int*>(), childCoupled(), atomics(), intDelay(intDelay), extDelay(extDelay), nEIC(), nIC(), nEOC() {
+        if (depth < 1 || width < 1) {
+            throw std::bad_exception();
+        }
+        if (depth == 1) {
+            auto *a = addAtomic();
+            couple(this, DEVStoneCoupled::in, a, DEVStoneAtomic::in);
+            nEIC += 1;
+            couple(a, DEVStoneAtomic::out, this, DEVStoneCoupled::out);
+            nEOC += 1;
+        }
+    }
+
+    DEVStoneAtomic* addAtomic() {
+        auto *a = new DEVStoneAtomic(intDelay, extDelay);
+        add(a);
+        atomics.push_back(a);
+        return a;
+    }
+
+    [[nodiscard]] unsigned long nEICs() const {
+        auto res = nEIC;
+        if (childCoupled != nullptr) {
+            res += childCoupled->nEICs();
+        }
+        return res;
+    }
+
+    [[nodiscard]] unsigned long nICs() const {
+        auto res = nIC;
+        if (childCoupled != nullptr) {
+            res += childCoupled->nICs();
+        }
+        return res;
+    }
+
+    [[nodiscard]] unsigned long nEOCs() const {
+        auto res = nEOC;
+        if (childCoupled != nullptr) {
+            res += childCoupled->nEOCs();
+        }
+        return res;
+    }
+
+    [[nodiscard]] unsigned long nAtomics() const {
+        auto res = atomics.size();
+        if (childCoupled != nullptr) {
+            res += childCoupled->nAtomics();
+        }
+        return res;
+    }
+
+    [[nodiscard]] unsigned long nInternals() const {
+        auto res = (childCoupled != nullptr) ? childCoupled->nInternals() : 0;
+        for (const auto& atomic: atomics) {
+            res += atomic->nInternals;
+        }
+        return res;
+    }
+
+    [[nodiscard]] unsigned long nExternals() const {
+        auto res = (childCoupled != nullptr) ? childCoupled->nExternals() : 0;
+        for (const auto& atomic: atomics) {
+            res += atomic->nExternals;
+        }
+        return res;
+    }
+
+    [[nodiscard]] unsigned long nEvents() const {
+        auto res = (childCoupled != nullptr) ? childCoupled->nEvents() : 0;
+        for (const auto& atomic: atomics) {
+            res += atomic->nEvents;
+        }
+        return res;
+    }
+};
+
+#endif //DEVSTONE_ADEVS_DEVSTONECOUPLED_HPP

devstone/adevs/src/hi.cpp

@@ -1,26 +0,0 @@
-//
-// Created by rcardenas on 3/10/20.
-//
-
-#include "hi.hpp"
-
-HI::HI(int depthIn, int widthIn, int intDelayIn, int extDelayIn, double procTimeIn):
-        DEVSWrapper(depthIn, widthIn, intDelayIn, extDelayIn, procTimeIn)
-{
-    if (depth > 1) {
-        HI *c =  new HI(depth - 1, width, intDelay, extDelay, procTime);
-        add(c);
-        couple(this, this->in, c, c->in);
-        couple(c, c->out, this, this->out);
-
-        for (auto i = atomics.begin(); i < atomics.end(); ++i) {
-            DummyAtomic *atomic = *i;
-            couple(this, this->in, atomic, atomic->in);
-
-            if(i != atomics.end() - 1) {
-                DummyAtomic *nextAtomic = *(i+1);
-                couple(atomic, atomic->in, nextAtomic, nextAtomic->in);
-            }
-        }
-    }
-}

devstone/adevs/src/hi.hpp

@@ -1,16 +1,32 @@
-//
-// Created by rcardenas on 3/10/20.
-//
-
 #ifndef DEVSTONE_ADEVS_HI_HPP
 #define DEVSTONE_ADEVS_HI_HPP
 
 #include "utils.hpp"
-#include "DEVSWrapper.hpp"
+#include "devstone_coupled.hpp"
 
-class HI: public DEVSWrapper {
+class HI: public DEVStoneCoupled {
 public:
-    HI(int depthIn, int widthIn, int intDelayIn, int extDelayIn, double procTimeIn);
+    HI(int width, int depth, int intDelay, int extDelay): DEVStoneCoupled(width, depth, intDelay, extDelay) {
+        if (depth > 1) {
+            childCoupled = new HI(width, depth - 1, intDelay, extDelay);
+            add(childCoupled);
+            couple(this, HI::in, childCoupled, HI::in);
+            nEIC += 1;
+            couple(childCoupled, HI::out, this, HI::out);
+            nEOC += 1;
+            DEVStoneAtomic* prevAtomic = nullptr;
+            for (auto i = 1; i < width; ++i) {
+                auto atomic = addAtomic();
+                couple(this, HI::in, atomic, DEVStoneAtomic::in);
+                nEIC += 1;
+                if (prevAtomic != nullptr) {
+                    couple(prevAtomic, DEVStoneAtomic::out, atomic, DEVStoneAtomic::in);
+                    nIC += 1;
+                }
+                prevAtomic = atomic;
+            }
+        }
+    }
 };
 
 

devstone/adevs/src/ho.cpp

@@ -1,27 +0,0 @@
-//
-// Created by rcardenas on 3/10/20.
-//
-
-#include "ho.hpp"
-
-HO::HO(int depthIn, int widthIn, int intDelayIn, int extDelayIn, double procTimeIn):
-        DEVSWrapper(depthIn, widthIn, intDelayIn, extDelayIn, procTimeIn)
-{
-    if (depth > 1) {
-        HO *c =  new HO(depth - 1, width, intDelay, extDelay, procTime);
-        add(c);
-        couple(this, this->in, c, c->in);
-        couple(c, c->out, this, this->out);
-
-        for (auto i = atomics.begin(); i < atomics.end(); ++i) {
-            DummyAtomic *atomic = *i;
-            couple(this, this->in, atomic, atomic->in);
-            couple(atomic, atomic->out, this, this->out);
-
-            if(i != atomics.end() - 1) {
-                DummyAtomic *nextAtomic = *(i+1);
-                couple(atomic, atomic->in, nextAtomic, nextAtomic->in);
-            }
-        }
-    }
-}

devstone/adevs/src/ho.hpp

@@ -1,18 +1,41 @@
-//
-// Created by rcardenas on 3/10/20.
-//
-
 #ifndef DEVSTONE_ADEVS_HO_HPP
 #define DEVSTONE_ADEVS_HO_HPP
 
 #include <adevs.h>
-#include "DEVStone.hpp"
-#include "DEVSWrapper.hpp"
+#include "devstone.hpp"
+#include "devstone_coupled.hpp"
 
-class HO: public DEVSWrapper {
+class HO: public DEVStoneCoupled {
 public:
-    HO(int depthIn, int widthIn, int intDelayIn, int extDelayIn, double procTimeIn);
+    /// Additional model input port
+    static const int in2 = 2;
+    /// Additional model output port
+    static const int out2 = 3;
+    HO(int width, int depth, int intDelay, int extDelay): DEVStoneCoupled(width, depth, intDelay, extDelay) {
+        if (depth > 1) {
+            childCoupled = new HO(width, depth - 1, intDelay, extDelay);
+            add(childCoupled);
+            couple(this, HO::in, childCoupled, HO::in);
+            nEIC += 1;
+            couple(this, HO::in2, childCoupled, HO::in2);
+            nEIC += 1;
+            couple(childCoupled, HO::out, this, HO::out);
+            nEOC += 1;
+            DEVStoneAtomic* prevAtomic = nullptr;
+            for (auto i = 1; i < width; ++i) {
+                auto atomic = addAtomic();
+                couple(this, HO::in2, atomic, DEVStoneAtomic::in);
+                nEIC += 1;
+                if (prevAtomic != nullptr) {
+                    couple(prevAtomic, DEVStoneAtomic::out, atomic, DEVStoneAtomic::in);
+                    nIC += 1;
+                }
+                couple(atomic, DEVStoneAtomic::out, this, HO::out2);
+                nEOC += 1;
+                prevAtomic = atomic;
+            }
+        }
+    }
 };
 
-
 #endif //DEVSTONE_ADEVS_HO_HPP

devstone/adevs/src/hoMod.cpp

@@ -1,71 +0,0 @@
-//
-// Created by rcardenas on 3/10/20.
-//
-
-#include "hoMod.hpp"
-#include <map>
-#include <list>
-
-const int HOmod::in1 = 0;
-const int HOmod::in2 = 1;
-const int HOmod::out = 2;
-
-HOmod::HOmod(int depthIn, int widthIn, int intDelayIn, int extDelayIn, double procTimeIn):
-adevs::Digraph<int*>()
-{
-    depth = depthIn;
-    width = widthIn;
-    intDelay = intDelayIn;
-    extDelay = extDelayIn;
-    procTime = procTimeIn;
-
-    if (depth == 1) {
-        DummyAtomic *a = new DummyAtomic(intDelay, extDelay, procTime);
-        add(a);
-        couple(this, this->in1, a, a->in);
-        couple(a, a->out, this, this->out);
-    } else {
-        HOmod *c = new HOmod(depth - 1, width, intDelay, extDelay, procTime);
-        couple(this, this->in1, c, c->in1);
-        couple(c, c->out, this, this->out);
-        if (width > 1) {
-            std::map<int, std::list<DummyAtomic*>> atomics;
-            for (int i = 0; i < width; i++) {
-                int minRowIdx = (i < 2)? 0 : (i - 1);
-                for (int j = minRowIdx; j < width - 1; j++) {
-                    DummyAtomic *a = new DummyAtomic(intDelay, extDelay, procTime);
-                    add(a);
-
-                    atomics.insert(std::pair<int, list<DummyAtomic*>>(i, std::list<DummyAtomic*>()));  // en teoria insert no machaca
-                    atomics.find(i)->second.push_back(a);
-                }
-            }
-            // Connect EIC
-            for (auto a: atomics.find(0)->second) {
-                couple(this, this->in2, a, a->in);
-            }
-            for (int i = 1; i < width; i++) {
-                DummyAtomic *a = atomics.find(i)->second.front();
-                couple(this, this->in2, a, a->in);
-            }
-            // Connect IC
-            for (auto a: atomics.find(0)->second) {
-                couple(a, a->out, c, c->in2);
-            }
-            for (int i = 0; i < atomics.find(0)->second.size(); i++) {
-                DummyAtomic *aTop = *std::next(atomics.find(0)->second.begin(), i);
-                for (int j = 0; j < atomics.find(1)->second.size(); j++) {
-                    DummyAtomic *aDown = *std::next(atomics.find(1)->second.begin(), j);
-                    couple(aDown, aDown->out, aTop, aTop->in);
-                }
-            }
-            for (int i = 2; i < width; i++) {
-                for (int j = 0; j < atomics.find(i)->second.size(); j++) {
-                    DummyAtomic *aFrom = *std::next(atomics.find(i)->second.begin(), j);
-                    DummyAtomic *aTo = *std::next(atomics.find(i - 1)->second.begin(), j + 1);
-                    couple(aFrom, aFrom->out, aTo, aTo->in);
-                }
-            }
-        }
-    }
-}

devstone/adevs/src/hoMod.hpp

@@ -1,27 +0,0 @@
-//
-// Created by rcardenas on 3/10/20.
-//
-
-#ifndef DEVSTONE_ADEVS_HOMOD_HPP
-#define DEVSTONE_ADEVS_HOMOD_HPP
-
-#include "utils.hpp"
-#include "DummyAtomic.hpp"
-
-
-class HOmod: public adevs::Digraph<int*> {
-protected:
-    int depth, width;
-    int intDelay, extDelay;
-    double procTime;
-public:
-    /// Model input ports
-    static const int in1, in2;
-    /// Model output port
-    static const int out;
-    // Constructor
-    HOmod(int depthIn, int widthIn, int intDelayIn, int extDelayIn, double procTimeIn);
-};
-
-
-#endif //DEVSTONE_ADEVS_HOMOD_HPP

devstone/adevs/src/homod.hpp

@@ -0,0 +1,64 @@
+#ifndef DEVSTONE_ADEVS_HOMOD_HPP
+#define DEVSTONE_ADEVS_HOMOD_HPP
+
+#include "utils.hpp"
+#include "devstone_coupled.hpp"
+
+class HOmod: public DEVStoneCoupled {
+ public:
+    /// Additional model input port
+    static const int in2 = 2;
+    HOmod(int width, int depth, int intDelay, int extDelay): DEVStoneCoupled(width, depth, intDelay, extDelay) {
+        if (depth > 1) {
+            childCoupled = new HOmod(width, depth - 1, intDelay, extDelay);
+            add(childCoupled);
+            couple(this, HOmod::in, childCoupled, HOmod::in);
+            nEIC += 1;
+            couple(childCoupled, HOmod::out, this, HOmod::out);
+            nEOC += 1;
+
+            std::vector<DEVStoneAtomic*> prevLayer, currentLayer;
+            // First layer of atomic models:
+            for (int i = 1; i < width; ++i) {
+                auto atomic = addAtomic();
+                couple(this, HOmod::in2, atomic, DEVStoneAtomic::in);
+                nEIC += 1;
+                couple(atomic, DEVStoneAtomic::out, childCoupled, HOmod::in2);
+                nIC += 1;
+                prevLayer.push_back(atomic);
+            }
+            // Second layer of atomic models:
+            for (int i = 1; i < width; ++i) {
+                auto atomic = addAtomic();
+                if (i == 1) {
+                    couple(this, HOmod::in2, atomic, DEVStoneAtomic::in);
+                    nEIC += 1;
+                }
+                for (const auto& prevAtomic: prevLayer) {
+                    couple(atomic, DEVStoneAtomic::out, prevAtomic, DEVStoneAtomic::in);
+                    nIC += 1;
+                }
+                currentLayer.push_back(atomic);
+            }
+            // Rest of the tree
+            prevLayer = currentLayer;
+            currentLayer = std::vector<DEVStoneAtomic*>();
+            for (int layer = 3; layer <= width; ++layer) {
+                for (unsigned long i = 1; i < prevLayer.size(); ++i) {
+                    auto atomic = addAtomic();
+                    if (i == 1) {
+                        couple(this, HOmod::in2, atomic, DEVStoneAtomic::in);
+                        nEIC += 1;
+                    }
+                    couple(atomic, DEVStoneAtomic::out, prevLayer.at(i), DEVStoneAtomic::in);
+                    nIC += 1;
+                    currentLayer.push_back(atomic);
+                }
+                prevLayer = currentLayer;
+                currentLayer = std::vector<DEVStoneAtomic*>();
+            }
+        }
+    }
+};
+
+#endif //DEVSTONE_ADEVS_HOMOD_HPP

devstone/adevs/src/li.cpp

@@ -1,20 +0,0 @@
-//
-// Created by rcardenas on 3/10/20.
-//
-
-#include "li.hpp"
-
-LI::LI(int depthIn, int widthIn, int intDelayIn, int extDelayIn, double procTimeIn):
-DEVSWrapper(depthIn, widthIn, intDelayIn, extDelayIn, procTimeIn)
-{
-    if (depth > 1) {
-        LI *c =  new LI(depth - 1, width, intDelay, extDelay, procTime);
-        add(c);
-        couple(this, this->in, c, c->in);
-        couple(c, c->out, this, this->out);
-
-        for(auto atomic: atomics) {
-            couple(this, this->in, atomic, atomic->in);
-        }
-    }
-}

devstone/adevs/src/li.hpp

@@ -1,16 +1,26 @@
-//
-// Created by rcardenas on 3/10/20.
-//
-
 #ifndef DEVSTONE_ADEVS_LI_HPP
 #define DEVSTONE_ADEVS_LI_HPP
 
 #include "utils.hpp"
-#include "DEVSWrapper.hpp"
+#include "devstone_coupled.hpp"
 
-class LI: public DEVSWrapper {
-public:
-    LI(int depthIn, int widthIn, int intDelayIn, int extDelayIn, double procTimeIn);
+class LI: public DEVStoneCoupled {
+  public:
+    LI(int width, int depth, int intDelay, int extDelay): DEVStoneCoupled(width, depth, intDelay, extDelay) {
+        if (depth > 1) {
+            childCoupled = new LI(width, depth - 1, intDelay, extDelay);
+            add(childCoupled);
+            couple(this, LI::in, childCoupled, LI::in);
+            nEIC += 1;
+            couple(childCoupled, LI::out, this, LI::out);
+            nEOC += 1;
+            for (auto i = 1; i < width; ++i) {
+                auto atomic = addAtomic();
+                couple(this, LI::in, atomic, DEVStoneAtomic::in);
+                nEIC += 1;
+            }
+        }
+    }
 };
 
 #endif //DEVSTONE_ADEVS_LI_HPP

devstone/adevs/src/seeder.cpp

@@ -1,49 +0,0 @@
-//
-// Created by rcardenas on 3/13/20.
-//
-
-#include "seeder.hpp"
-
-const int Seeder::out = 0;
-
-Seeder::Seeder(int nMessagesIn): adevs::Atomic<IO_Type>() {
-    nMessages = nMessagesIn;
-    sigma = 0;
-    t = 0.0;
-}
-
-/// External transition function
-void Seeder::delta_ext(double e, const adevs::Bag<IO_Type>& x) { }
-
-/// Internal transition function
-void Seeder::delta_int() {
-    t += sigma;
-    sigma = DBL_MAX;
-}
-
-/// Confluent transition function.
-void Seeder::delta_conf(const adevs::Bag<IO_Type>& x) {
-    // Discard the old job
-    delta_int();
-    // Process the incoming job
-    delta_ext(0.0, x);
-}
-
-/// Output function.
-void Seeder::output_func(adevs::Bag<IO_Type >& y) {
-    // Get the departing customer
-    for (int i = 0; i < nMessages; i++) {
-        y.insert(IO_Type(out, 0));
-    }
-}
-
-/// Time advance function.
-double Seeder::ta() {
-    return sigma;
-}
-
-/// Garbage collection. No heap allocation in output_func, so do nothing.
-void Seeder::gc_output(adevs::Bag<IO_Type>&) {}
-
-/// Destructor
-Seeder::~Seeder() = default;

devstone/adevs/src/seeder.hpp

@@ -1,7 +1,3 @@
-//
-// Created by rcardenas on 3/13/20.
-//
-
 #ifndef DEVSTONE_ADEVS_SEEDER_HPP
 #define DEVSTONE_ADEVS_SEEDER_HPP
 
@@ -9,33 +5,39 @@
 #include <list>
 #include "utils.hpp"
 
-
 class Seeder: public adevs::Atomic<IO_Type> {
 private:
     /// Model state variables
     double sigma;
-    int nMessages;
-    double t;
 public:
     /// Model output port
-    static const int out;
+    static const int out = 0;
+
+    /// Constructor.
+    explicit Seeder(): adevs::Atomic<IO_Type>(), sigma() {}
+
+    /// External transition function (does not apply)
+    void delta_ext(__attribute__ ((unused)) double e, __attribute__ ((unused)) const adevs::Bag<IO_Type>& x) override {}
 
-    /// Constructor.  The processing time is provided as an argument.
-    Seeder(int nMessagesIn);
-    /// External transition function
-    void delta_ext(double e, const adevs::Bag<IO_Type>& x);
-    /// Internal transition function
-    void delta_int();
-    /// Confluent transition function.
-    void delta_conf(const adevs::Bag<IO_Type>& x);
-    /// Output function.
-    void output_func(adevs::Bag<IO_Type>& y);
-    /// Time advance function.
-    double ta();
+    /// Internal transition function (just passivates)
+    void delta_int() override {
+            sigma = std::numeric_limits<double>::infinity();
+    }
+    /// Confluent transition function (standard one)
+    void delta_conf(const adevs::Bag<IO_Type>& x) override {
+        delta_int();
+        delta_ext(0.0, x);
+    }
+    /// Output function (sends a single message via the single port).
+    void output_func(adevs::Bag<IO_Type>& y) override {
+        y.insert(IO_Type(out, 0));
+    }
+    /// Time advance function (just returns sigma)
+    double ta() override {
+        return sigma;
+    }
     /// Garbage collection.
-    void gc_output(adevs::Bag<IO_Type>&);
-    /// Destructor
-    ~Seeder();
+    void gc_output(adevs::Bag<IO_Type>&) override {}
 };
 
 #endif //DEVSTONE_ADEVS_SEEDER_HPP

devstone/adevs/src/utils.hpp

@@ -1,7 +1,3 @@
-//
-// Created by rcardenas on 3/13/20.
-//
-
 #ifndef DEVSTONE_ADEVS_UTILS_HPP
 #define DEVSTONE_ADEVS_UTILS_HPP
 

devstone/adevs/test/test_devstone.cpp

@@ -0,0 +1,156 @@
+#define BOOST_TEST_MODULE DEVStoneTests
+#include <boost/test/unit_test.hpp>
+#include <string>
+#include "../src/devstone.hpp"
+
+#define STEP 5
+#define MAX_WIDTH 50
+#define MAX_DEPTH 50
+
+int expectedAtomics(const std::string& type, int width, int depth) {
+    auto wFactor = width - 1;
+    if (type == "HOmod") {
+        wFactor += (width - 1) * width / 2;
+    }
+    return wFactor * (depth - 1) + 1;
+}
+
+int expectedEICs(const std::string& type, int width, int depth) {
+    auto wFactor = width;
+    if (type == "HO") {
+        wFactor = width + 1;
+    } else if (type == "HOmod") {
+        wFactor = 2 * (width - 1) + 1;
+    }
+    return wFactor * (depth - 1) + 1;
+}
+
+int expectedICs(const std::string& type, int width, int depth) {
+    auto wFactor = 0;
+    if (width > 2 && (type == "HI" || type == "HO")) {
+        wFactor = (width - 2);
+    } else if (type == "HOmod") {
+        wFactor = (width - 1) * (width - 1) + (width - 1) * width / 2;
+    }
+    return wFactor * (depth - 1);
+}
+
+int expectedEOCs(const std::string& type, int width, int depth) {
+    return (type == "HO")? width * (depth - 1) + 1 : depth;
+}
+
+int expectedInternals(const std::string& type, int width, int depth) {
+    if (type == "LI") {
+        return (width - 1) * (depth - 1) + 1;
+    } else if (type == "HI" || type == "HO") {
+        return (width - 1) * width / 2 * (depth - 1) + 1;
+    }
+    auto n = 1;
+    for (int d = 1; d < depth; ++d) {
+        n += (1 + (d - 1) * (width - 1)) * (width - 1) * width / 2 + (width - 1) * (width + (d - 1) * (width - 1));
+    }
+    return n;
+}
+
+int expectedExternals(const std::string& type, int width, int depth) {
+    return expectedInternals(type, width, depth);
+}
+
+int expectedEvents(const std::string& type, int width, int depth) {
+    if (type != "HOmod") {
+        return expectedInternals(type, width, depth);
+    }
+    auto n = 1;
+    if (width > 1 && depth > 1) {
+        n += 2 * (width - 1);
+        auto aux = 0;
+        for (int i = 2; i < depth; ++i) {
+            aux += 1 + (i - 1) * (width - 1);
+        }
+        n += aux * 2 * (width - 1) * (width - 1);
+        n += (aux + 1) * ((width - 1) * (width - 1) + (width - 2) * (width - 1) / 2);
+    }
+    return n;
+}
+
+adevs::Simulator<IO_Type> createEngine(DEVStone* devstone) {
+    return {devstone};
+}
+
+[[maybe_unused]] void runSimulation(adevs::Simulator<IO_Type>& sim) {
+    while (sim.nextEventTime() < DBL_MAX) {
+        sim.execNextEvent();
+    }
+}
+
+BOOST_AUTO_TEST_CASE(DEVStoneLI)
+{
+    for (int w = 1; w <= MAX_WIDTH; w += STEP) {
+        for (int d = 1; d <= MAX_DEPTH; d += STEP) {
+            auto coupled = new DEVStone("LI", w, d, 0, 0);
+            BOOST_CHECK_EQUAL(coupled->model->nAtomics(), expectedAtomics("LI", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nEICs(), expectedEICs("LI", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nICs(), expectedICs("LI", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nEOCs(), expectedEOCs("LI", w, d));
+            auto coordinator = createEngine(coupled);
+            runSimulation(coordinator);
+            BOOST_CHECK_EQUAL(coupled->model->nInternals(), expectedInternals("LI", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nExternals(), expectedExternals("LI", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nEvents(), expectedEvents("LI", w, d));
+        }
+    }
+}
+
+BOOST_AUTO_TEST_CASE(DEVStoneHI)
+{
+    for (int w = 1; w <= MAX_WIDTH; w += STEP) {
+        for (int d = 1; d <= MAX_DEPTH; d += STEP) {
+            auto coupled = new DEVStone("HI", w, d, 0, 0);
+            BOOST_CHECK_EQUAL(coupled->model->nAtomics(), expectedAtomics("HI", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nEICs(), expectedEICs("HI", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nICs(), expectedICs("HI", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nEOCs(), expectedEOCs("HI", w, d));
+            auto coordinator = createEngine(coupled);
+            runSimulation(coordinator);
+            BOOST_CHECK_EQUAL(coupled->model->nInternals(), expectedInternals("HI", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nExternals(), expectedExternals("HI", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nEvents(), expectedEvents("HI", w, d));
+        }
+    }
+}
+
+BOOST_AUTO_TEST_CASE(DEVStoneHO)
+{
+    for (int w = 1; w <= MAX_WIDTH; w += STEP) {
+        for (int d = 1; d <= MAX_DEPTH; d += STEP) {
+            auto coupled = new DEVStone("HO", w, d, 0, 0);
+            BOOST_CHECK_EQUAL(coupled->model->nAtomics(), expectedAtomics("HO", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nEICs(), expectedEICs("HO", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nICs(), expectedICs("HO", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nEOCs(), expectedEOCs("HO", w, d));
+            auto coordinator = createEngine(coupled);
+            runSimulation(coordinator);
+            BOOST_CHECK_EQUAL(coupled->model->nInternals(), expectedInternals("HO", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nExternals(), expectedExternals("HO", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nEvents(), expectedEvents("HO", w, d));
+        }
+    }
+}
+
+BOOST_AUTO_TEST_CASE(DEVStoneHOmod)
+{
+    for (int w = 1; w <= MAX_WIDTH / 5; w += STEP / 5) {
+        for (int d = 1; d <= MAX_DEPTH / 5; d += STEP / 5) {
+            auto coupled = new DEVStone("HOmod", w, d, 0, 0);
+            BOOST_CHECK_EQUAL(coupled->model->nAtomics(), expectedAtomics("HOmod", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nEICs(), expectedEICs("HOmod", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nICs(), expectedICs("HOmod", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nEOCs(), expectedEOCs("HOmod", w, d));
+            auto coordinator = createEngine(coupled);
+            runSimulation(coordinator);
+            BOOST_CHECK_EQUAL(coupled->model->nInternals(), expectedInternals("HOmod", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nExternals(), expectedExternals("HOmod", w, d));
+            BOOST_CHECK_EQUAL(coupled->model->nEvents(), expectedEvents("HOmod", w, d));
+        }
+    }
+}

setup.sh

@@ -28,16 +28,12 @@ cd ../..
 cd devstone
 # Set up aDEVS project
 cd adevs
-mkdir build
-cd build
-cmake ..
-make
+cmake -S . -B build/ -D CMAKE_BUILD_TYPE=Release
+cmake --build build/ --target devstone
 # Set up Cadmium project
 cd ../../cadmium
-mkdir build
-cd build
-cmake ..
-make
+cmake -S . -B build/ -D CMAKE_BUILD_TYPE=Release
+cmake --build build/ --target devstone
 # Set up CDBoost project
 cd ../../cdboost
 mkdir build

simulators/adevs-3.3/README

@@ -0,0 +1,81 @@
+This file contains instructions for building the C++ library and
+the Java bindings. Development versions of this library are available 
+from SourceForge and github.
+
+---- C++ -----
+Adevs is implemented primarily as a set of C++ header files,
+and for many applications these are sufficient. The rest of
+the simulation engine is in the static library file 'libadevs.a'.
+To build this library file this, go to the directory called 'src'
+and run 'make'. The default options for the compiler should
+work if you are using the GNU g++ compiler. For other compilers,
+you make set options, build commands, etc. by editing the Makefile.
+
+If you are using Visual Studio and Windows, then to build the
+static library do the following: 1) open a command shell, 2) go the
+'src' directory, and 3) run the batch file 'build.bat'. This 
+creates a static library to link with your simulation applications.
+
+---- Java -----
+To build the java bindings for Linux, go the the 'src' directory and  
+edit the Makefile so that JAVA_HOME points to your jdk installation.
+Then run 'make java_adevs'. This builds the shared library libjava_adevs.so
+and the jar-file adevs.jar. To use these put the shared library into 
+your LD_LIBRARY_PATH and the jar file into your CLASSPATH.
+
+To build the java bindings For Windows, first make sure that 1) your
+JAVA_HOME environment variable points to your jdk installation and
+2) that the java compiler is in your executable PATH. Then execute
+the batch file src/adevs_jni/build.bat. This creates the jar-file
+adevs.jar and DLL-file java_adevs.dll. You can use these by putting
+the former into your CLASSPATH and latter into your PATH.
+
+---- FMI ----
+Adevs supports importing models that use the FMI standard for Model
+Exchange. Information about this standard can be found at
+http://www.fmi-standard.org. This feature is particularly useful
+for importing continuous system models built with the Modelica
+language into your discrete event simulation. The examples and
+test cases shipped with adevs assume that you have the 
+OpenModelica compiler installed. You can get this compiler from
+http://www.openmodelica.org. The bash script build-omc.sh for building
+the OpenModelica compiler from source code is located in the util
+directory. The fmi test cases are located in test/fmi, and you can
+execute these by going to that directory and running 'make all'.
+
+---- Useful scripts ----
+
+xml2cpp.py
+This script will convert an FMI model description file to a C++ class
+that imports the FMI described into an adevs simulation by using the
+adevs::FMI class. Executing the script will display rudimentary information
+about command line arguments.
+
+build-omc.sh
+This bash shell script will download and build the OpenModelica compiler,
+which can be used to translate Modelica models into FMI modules. These FMI
+modules can then be loaded into your adevs simulation.
+
+---- Testing -----
+There are four sets of test cases that can be run from the
+directory 'test'. The first set are for the C++ library without
+and these tests are run with 'make check_cpp'.
+The java bindings can be tested by
+executing 'make check_java'. The fmi import support can be
+tested by executing 'make check_fmi'. Executing 'make check'
+runs all of the test cases.
+
+The make.common file contains platform specific settings for
+the test scripts. It is setup by default for a GNU/Linux type of
+environment. Change as needed for your environment, but if you are
+using a more or less typical Linux distribution then the defaults
+should work fine.
+
+When testing the Java bindings, be sure your PATH, LD_LIBRARY_PATH,
+and CLASSPATH are setup properly before running the Java test
+cases. These need to include your DLL (for Windows) and Java binaries,
+shared object file (for Linux), and jar file and the '.' directory
+respectively.
+
+For the FMI test cases, you will need to have the OpenModelica compiler
+installed.

simulators/adevs-3.3/copyright.txt

@@ -0,0 +1,30 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */

simulators/adevs-3.3/include/adevs.h

@@ -0,0 +1,51 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+/**
+ * \mainpage
+ * Links to the user manual, tutorials on DEVS, and other such
+ * stuff are <a href="../index.html">here</a>.
+ */
+/**
+ * Include everything needed to compile user models.
+ */
+#include "adevs_exception.h"
+#include "adevs_models.h"
+#include "adevs_simulator.h"
+#include "adevs_digraph.h"
+#include "adevs_simpledigraph.h"
+#include "adevs_cellspace.h"
+#include "adevs_rand.h"
+#include "adevs_hybrid.h"
+#include "adevs_corrected_euler.h"
+#include "adevs_event_locators.h"
+#include "adevs_rk_45.h"
+#include "adevs_poly.h"
+#include "adevs_wrapper.h"

simulators/adevs-3.3/include/adevs_abstract_simulator.h

@@ -0,0 +1,119 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef __adevs_abstract_simulator_h_
+#define __adevs_abstract_simulator_h_
+#include "adevs_models.h"
+#include "adevs_event_listener.h"
+#include "adevs_bag.h"
+
+namespace adevs
+{
+
+/**
+ * This is the base class for all simulators. It defines an interface that is
+ * supported by all derived classes and provides some basic helper routines
+ * for those derived classes.
+ */
+template <class X, class T = double> class AbstractSimulator 
+{
+	public:
+		/// Constructor
+		AbstractSimulator(){}
+		/**
+		 * Add an event listener that will be notified of output events 
+		 * produced by all components within the model.
+		 * @param l The listener to be notified
+		 */
+		void addEventListener(EventListener<X,T>* l)
+		{
+			listeners.insert(l);
+		}
+		/**
+		 * Remove an event listener that was previously added.
+		 * @param l The listener to be removed
+		 */
+		void removeEventListener(EventListener<X,T>* l)
+		{
+			listeners.erase(l);
+		}
+		/// Get the model's next event time
+		virtual T nextEventTime() = 0;
+		/// Execute the simulator until the next event time is greater than tend
+		virtual T execUntil(T tend) = 0;
+		/// Destructor leaves the model intact.
+		virtual ~AbstractSimulator(){}
+		/// Notify listeners of an output event.
+		void notify_output_listeners(Devs<X,T>* model, const X& value, T t);
+		/// Notify listeners of an input event.
+		void notify_input_listeners(Devs<X,T>* model, const X& value, T t);
+		/// Notify listeners of a state change.
+		void notify_state_listeners(Atomic<X,T>* model, T t);
+	private:
+		/// Eternal event listeners
+		Bag<EventListener<X,T>*> listeners;
+
+};
+
+template <class X, class T>
+void AbstractSimulator<X,T>::notify_output_listeners(Devs<X,T>* model, const X& value, T t)
+{
+	Event<X,T> event(model,value);
+	typename Bag<EventListener<X,T>*>::iterator iter;
+	for (iter = listeners.begin(); iter != listeners.end(); iter++)
+	{
+		(*iter)->outputEvent(event,t);
+	}
+}
+
+template <class X, class T>
+void AbstractSimulator<X,T>::notify_input_listeners(Devs<X,T>* model, const X& value, T t)
+{
+	Event<X,T> event(model,value);
+	typename Bag<EventListener<X,T>*>::iterator iter;
+	for (iter = listeners.begin(); iter != listeners.end(); iter++)
+	{
+		(*iter)->inputEvent(event,t);
+	}
+}
+
+template <class X, class T>
+void AbstractSimulator<X,T>::notify_state_listeners(Atomic<X,T>* model, T t)
+{
+	typename Bag<EventListener<X,T>*>::iterator iter;
+	for (iter = listeners.begin(); iter != listeners.end(); iter++)
+	{
+		(*iter)->stateChange(model,t);
+	}
+}
+
+} // end of namespace
+
+#endif

simulators/adevs-3.3/include/adevs_bag.h

@@ -0,0 +1,193 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef _adevs_bag_h
+#define _adevs_bag_h
+#include <cstdlib>
+
+namespace adevs
+{
+
+/**
+ * The Bag is (almost) a model of a STL Multiple Associative Container.
+ * This implementation is optimized for the simulation engine, and it 
+ * does not satisfy the STL complexity requirements. Neither does it implement
+ * the full set of required methods, but those methods that are implemented
+ * conform to the standard (except for the time complexity requirement).
+ */
+template <class T> class Bag
+{
+	public:
+		/// A bidirectional iterator for the Bag
+		class iterator
+		{
+			public:
+				/// Default constructor
+				iterator(){}
+				/// Copy constructor
+				iterator(const iterator& src):i(src.i),b(src.b){}
+				/// Copy by assignment
+				const iterator& operator=(const iterator& src) 
+				{ 
+					i = src.i; 
+					b = src.b;
+					return *this;
+				}
+				/// True if two iterators point at the same position in the bag
+				bool operator==(const iterator& src) const { return i==src.i; }
+				/// True if two iterators point at different positions in the bag
+				bool operator!=(const iterator& src) const { return i!=src.i; }
+				/// Get the element at the current position
+				T& operator*() { return b[i]; }
+				/// Get the element at the current position
+				const T& operator*() const { return b[i]; }
+				/// Move to the next position and get the element there
+				iterator& operator++() { i++; return *this; }
+				/// Move to the next position and get the element there
+				iterator& operator++(int) { ++i; return *this; }
+				/// Move to the previous position and get the element there
+				iterator& operator--() { i--; return *this; }
+				/// Move to the previous position and get the element there
+				iterator& operator--(int) { --i; return *this; }
+			private:
+				friend class Bag<T>;	
+				unsigned int i;
+				T* b;
+				iterator(unsigned int start, T* b):i(start),b(b){}
+		};
+		/// A bidirection iterator for the Bag that should not change its contents 
+		typedef iterator const_iterator;
+		/// Create an empty bag with an initial capacity
+		Bag(unsigned int cap = 8):
+		cap_(cap),size_(0),b(new T[cap]){}
+		/// Copy constructor uses the = operator of T
+		Bag(const Bag<T>& src):
+		cap_(src.cap_),
+		size_(src.size_)
+		{
+			b = new T[src.cap_];
+			for (unsigned int i = 0; i < size_; i++) 
+				b[i] = src.b[i];
+		}
+		/// Assignment operator uses the = operator of T
+		const Bag<T>& operator=(const Bag<T>& src)
+		{
+			cap_ = src.cap_;
+			size_ = src.size_;
+			delete [] b;
+			b = new T[src.cap_];
+			for (unsigned int i = 0; i < size_; i++) 
+				b[i] = src.b[i];
+			return *this;
+		}
+		/// Swaps contents of this bag with the contents of the supplied bag. Returns this bag.
+		Bag<T>& swap(Bag<T>& src)
+		{
+			unsigned tmp_cap_, tmp_size_;
+			T* tmp_b;
+			tmp_cap_ = src.cap_;
+			tmp_size_ = src.size_;
+			tmp_b = src.b;
+			src.cap_ = cap_;
+			src.size_ = size_;
+			src.b = b;
+			cap_ = tmp_cap_;
+			size_ = tmp_size_;
+			b = tmp_b;
+			return *this;
+		}
+		/// Count the instances of a stored in the bag
+		unsigned count(const T& a) const
+		{
+			unsigned result = 0;
+			for (unsigned i = 0; i < size_; i++)
+				if (b[i] == a) result++;
+			return result;
+		}
+		/// Get the number of elements in the bag
+		unsigned size() const { return size_; }
+		/// Same as size()==0
+		bool empty() const { return size_ == 0; }
+		/// Get an iterator pointing to the first element in the bag
+		iterator begin() const { return iterator(0,b); }
+		/// Get an iterator to the end of the bag (i.e., just after the last element)
+		iterator end() const { return iterator(size_,b); }
+		/// Erase the first instance of k
+		void erase(const T& k) 
+		{
+			iterator p = find(k);
+			if (p != end()) erase(p);
+		}
+		/// Erase the element pointed to by p
+		void erase(iterator p)
+		{
+			size_--;
+			b[p.i] = b[size_];
+		}
+		/// Remove all of the elements from the bag
+		void clear() { size_ = 0; }
+		/// Find the first instance of k, or end() if no instance is found. Uses == for comparing T.
+		iterator find(const T& k) const
+		{
+			for (unsigned i = 0; i < size_; i++)
+				if (b[i] == k) return iterator(i,b);
+			return end();
+		}
+		/// Put t into the bag
+		void insert(const T& t)
+		{
+			if (cap_ == size_) enlarge(2*cap_);
+			b[size_] = t;
+			size_++;
+		}
+		/// Access the element at position k
+		const T& operator[](unsigned k) const { return b[k]; }
+		/// Access the element at position k
+		T& operator[](unsigned k) { return b[k]; }
+		/// Delete the bag 
+		~Bag() { delete [] b; }
+	private:	
+		unsigned cap_, size_;
+		T* b;
+		/// Adds the specified capacity to the bag.
+		void enlarge(unsigned adjustment)
+		{
+			cap_ = cap_ + adjustment;
+			T* rb = new T[cap_];
+			for (unsigned i = 0; i < size_; i++) 
+				rb[i] = b[i];
+			delete [] b;
+			b = rb;
+		}
+	};
+
+} // end of namespace
+
+#endif

simulators/adevs-3.3/include/adevs_cellspace.h

@@ -0,0 +1,216 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef __adevs_cellspace_h_
+#define __adevs_cellspace_h_
+#include "adevs.h"
+#include <cstdlib>
+
+namespace adevs
+{
+
+/**
+ * Input and output events produced by components of a CellSpace must
+ * be of the type CellEvent.  A CellEvent has an event value (i.e., the actual
+ * input/output value) and a target cell for the event.
+ */
+template <class X> class CellEvent
+{
+	public:
+		/// Default constructor. Sets x = y = z = 0.
+		CellEvent():value() { x = y = z = 0; }
+		/// Copy constructor
+		CellEvent(const CellEvent<X>& src):
+		x(src.x),y(src.y),z(src.z),value(src.value){}
+		/// Assignment operator
+		const CellEvent& operator=(const CellEvent<X>& src)
+		{
+			x = src.x; y = src.y; z = src.z; value = src.value;
+			return *this;
+		}
+		/// The x coordinate of the event target
+		long int x;
+		/// The y coordinate of the event target
+		long int y;
+		/// The z coordinate of the event target
+		long int z;
+		/// The event value
+		X value;
+};
+
+/**
+ * This class describes a 3D cell space whose components accept and produce CellEvent objects.
+ * This class is meant to be useful for solving PDEs, simulating
+ * next event cell spaces, and for building other types of models represented as a space of 
+ * discrete, interacting points. Output events produced by component models must be of type CellEvent, and
+ * the CellEvent (x,y,z) coordinate indicates the
+ * target cell for the event. The corresponding input event will have the same (x,y,z) value as
+ * the output event. Targets that are outside of the CellSpace will become external output
+ * events for the CellSpace model.  Similarly, CellEvent objects that are injected into the
+ * CellSpace (i.e., external input events) will be delivered to the targeted cell.
+ */
+template <class X, class T = double> class CellSpace: public Network<CellEvent<X>,T>
+{
+	public:
+		/// A component model in the CellSpace
+		typedef Devs<CellEvent<X>,T> Cell;
+		/// Create an Width x Height x Depth CellSpace with NULL entries in the cell locations.
+		CellSpace(long int width, long int height = 1, long int depth = 1);
+		/// Insert a model at the x,y,z position.
+		void add(Cell* model, long int x, long int y = 0, long int z = 0) 
+		{
+			space[x][y][z] = model;
+			model->setParent(this);
+		}
+		/// Get the model at location x,y,z.
+		const Cell* getModel(long int x, long int y = 0, long int z = 0) const
+		{
+			return space[x][y][z];
+		}
+		/// Get a mutable version of the model at x,y,z.
+		Cell* getModel(long int x, long int y = 0, long int z = 0)
+		{
+			return space[x][y][z];
+		}
+		/// Get the width of the CellSpace.
+		long int getWidth() const { return w; }
+		/// Get the height of the CellSpace.
+		long int getHeight() const { return h; }
+		/// Get the depth of the CellSpace.
+		long int getDepth() const { return d; }
+		/// Get the model's set of components
+		void getComponents(Set<Cell*>& c);
+		/// Route events within the Cellspace
+		void route(const CellEvent<X>& event, Cell* model, 
+		Bag<Event<CellEvent<X>,T> >& r);
+		/// Destructor; this destroys the components as well.
+		~CellSpace();
+	private:	
+		long int w, h, d;
+		Cell**** space;
+};
+
+// Implementation of constructor
+template <class X, class T>
+CellSpace<X,T>::CellSpace(long int width, long int height, long int depth):
+Network<CellEvent<X>,T>()
+{
+	w = width;
+	h = height;
+	d = depth;
+	// Allocate space for the cells and set the entries to NULL
+	space = new Cell***[w];
+	for (long int x = 0; x < w; x++)
+	{
+		space[x] = new Cell**[h];
+		for (long int y = 0; y < h; y++)
+		{
+			space[x][y] = new Cell*[h];
+			for (long int z = 0; z < d; z++)
+			{
+				space[x][y][z] = NULL;
+			}
+		}
+	}
+}
+
+// Implementation of destructor
+template <class X, class T>
+CellSpace<X,T>::~CellSpace()
+{
+	for (long int x = 0; x < w; x++)
+	{
+		for (long int y = 0; y < h; y++)
+		{
+			for (long int z = 0; z < d; z++)
+			{
+				if (space[x][y][z] != NULL)
+				{
+					delete space[x][y][z];
+				}
+			}
+			delete [] space[x][y];
+		}
+		delete [] space[x];
+	}
+	delete [] space;
+}
+
+// Implementation of the getComponents() method
+template <class X, class T>
+void CellSpace<X,T>::getComponents(Set<Cell*>& c)
+{
+	// Add all non-null entries to the set c
+	for (long int x = 0; x < w; x++)
+	{
+		for (long int y = 0; y < h; y++)
+		{
+			for (long int z = 0; z < d; z++)
+			{
+				if (space[x][y][z] != NULL)
+				{
+					c.insert(space[x][y][z]);
+				}
+			}
+		}
+	}
+}
+
+// Event routing function for the net_exec
+template <class X, class T>
+void CellSpace<X,T>::route(
+const CellEvent<X>& event, Cell* model, Bag<Event<CellEvent<X>,T> >& r)
+{
+	Cell* target = NULL;
+	// If the target cell is inside of the cellspace
+	if (event.x >= 0 && event.x < w &&  // check x dimension
+	event.y >= 0 && event.y < h && // check y dimension
+	event.z >= 0 && event.z < d) // check z dimension
+	{
+		// Get the interior target
+		target = space[event.x][event.y][event.z];
+	}
+	else
+	{
+		// Otherwise, the event becomes an external output from the cellspace
+		target = this;
+	}
+	// If the target exists
+	if (target != NULL)
+	{
+		// Add an appropriate event to the receiver bag
+		Event<CellEvent<X>,T> io(target,event);
+		r.insert(io);
+	}
+}
+
+} // end of namespace
+
+#endif

simulators/adevs-3.3/include/adevs_corrected_euler.h

@@ -0,0 +1,143 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef _adevs_corrected_euler_h_
+#define _adevs_corrected_euler_h_
+#include <cmath>
+#include <algorithm>
+#include "adevs_hybrid.h"
+
+namespace adevs
+{
+
+/**
+ * This is the second order accurate RK2 method with adaptive step sizing for
+ * error control.
+ */
+template <typename X> class corrected_euler:
+	public ode_solver<X>
+{
+	public:
+		/**
+		 * Create an integrator that will use the specified per step error
+		 * tolerance and maximum step size.
+		 */
+		corrected_euler(ode_system<X>* sys, double err_tol, double h_max);
+		/// Destructor
+		~corrected_euler();
+		double integrate(double* q, double h_lim);
+		void advance(double* q, double h);
+	private:
+		double *dq, // derivative
+			   *qq, // trial solution
+			   *t,  // temporary variable for computing k2
+			   *k[2]; // k1 and k2
+		const double err_tol; // Error tolerance
+		const double h_max; // Maximum time step
+		double h_cur; // Previous time step that satisfied error constraint
+		// Compute a step of size h, put it in qq, and return the error
+		double trial_step(double h);
+};
+
+template <typename X>
+corrected_euler<X>::corrected_euler(ode_system<X>* sys, double err_tol,
+		double h_max):
+	ode_solver<X>(sys),err_tol(err_tol),h_max(h_max),h_cur(h_max)
+{
+	for (int i = 0; i < 2; i++) k[i] = new double[sys->numVars()];
+	dq = new double[sys->numVars()];
+	qq = new double[sys->numVars()];
+	t = new double[sys->numVars()];
+}
+
+template <typename X>
+corrected_euler<X>::~corrected_euler()
+{
+	delete [] t; delete [] qq; delete [] dq;
+	for (int i = 0; i < 2; i++) delete [] k[i];
+}
+
+template <typename X>
+void corrected_euler<X>::advance(double* q, double h)
+{
+	double dt;
+	while ((dt = integrate(q,h)) < h) h -= dt;
+}
+
+template <typename X>
+double corrected_euler<X>::integrate(double* q, double h_lim)
+{
+	// Initial error estimate and step size
+	double err = DBL_MAX,
+		   h = std::min<double>(h_cur*1.1,std::min<double>(h_max,h_lim));
+	for (;;) {
+		// Copy q to the trial vector
+		for (int i = 0; i < this->sys->numVars(); i++) qq[i] = q[i];
+		// Make the trial step which will be stored in qq
+		err = trial_step(h);
+		// If the error is ok, then we have found the proper step size
+		if (err <= err_tol) { // Keep h if shrunk to control the error
+			if (h_lim >= h_cur) h_cur = h; 
+			break;
+		}
+		// Otherwise shrink the step size and try again
+		else {
+			double h_guess = 0.8*err_tol*h/fabs(err);
+			if (h < h_guess) h *= 0.8;
+			else h = h_guess;
+		}
+	}
+	// Put the trial solution in q and return the selected step size
+	for (int i = 0; i < this->sys->numVars(); i++) q[i] = qq[i];
+	return h;
+}
+
+template <typename X>
+double corrected_euler<X>::trial_step(double step)
+{
+	int j;
+	// Compute k1
+	this->sys->der_func(qq,dq); 
+	for (j = 0; j < this->sys->numVars(); j++) k[0][j] = step*dq[j];
+	// Compute k2
+	for (j = 0; j < this->sys->numVars(); j++) t[j] = qq[j] + 0.5*k[0][j];
+	this->sys->der_func(t,dq);
+	for (j = 0; j < this->sys->numVars(); j++) k[1][j] = step*dq[j];
+	// Compute next state and approximate error
+	double err = 0.0;
+	for (j = 0; j < this->sys->numVars(); j++) {
+		qq[j] += k[1][j]; // Next state
+		err = std::max<double>(err,fabs(k[0][j]-k[1][j])); // Maximum error
+	}
+	return err; // Return the error
+}
+
+} // end of namespace
+#endif

simulators/adevs-3.3/include/adevs_digraph.h

@@ -0,0 +1,212 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef __adevs_digraph_h_
+#define __adevs_digraph_h_
+#include "adevs.h"
+#include <cassert>
+#include <map>
+#include <set>
+#include <cstdlib>
+
+namespace adevs
+{
+
+/**
+ * The components of a digraph model must use PortValue objects
+ * as their basic I/O type: the port and value types are template
+ * arguments. The default port type is an integer.
+ */
+template <class VALUE, class PORT=int> class PortValue
+{
+	public:
+		/// Constructor
+		PortValue():
+		port(),
+		value()
+		{
+		}
+		/// Copy constructor
+		PortValue(const PortValue& src):
+		port(src.port),
+		value(src.value)
+		{
+		}
+		/// Create an object with the specified port and value
+		PortValue(PORT port, const VALUE& value):
+		port(port),
+		value(value)
+		{
+		}
+		/// Assignment operator
+		const PortValue<VALUE,PORT>& operator=(const PortValue<VALUE,PORT>& src)
+		{
+			port = src.port;
+			value = src.value;
+			return *this;
+		}
+		/// Destructor
+		~PortValue()
+		{
+		}
+		/// The port on which the value appears
+		PORT port;
+		/// The value appearing on the port
+		VALUE value;
+};
+
+/**
+ * The digraph model is used to build block-diagrams from network and atomic components.
+ * Its components must have PortValue objects as their input/output type.
+ */
+template <class VALUE, class PORT=int, class T = double> class Digraph: 
+public Network<PortValue<VALUE,PORT>,T>
+{
+	public:
+		/// An input or output to a component model
+		typedef PortValue<VALUE,PORT> IO_Type;
+		/// A component of the Digraph model
+		typedef Devs<IO_Type,T> Component;
+
+		/// Construct a network with no components.
+		Digraph():
+		Network<IO_Type,T>()
+		{
+		}
+		/// Add a model to the network.
+		void add(Component* model);
+		/// Couple the source model to the destination model.  
+		void couple(Component* src, PORT srcPort, 
+		Component* dst, PORT dstPort);
+		/// Puts the network's components into to c
+		void getComponents(Set<Component*>& c);
+		/// Route an event based on the coupling information.
+		void route(const IO_Type& x, Component* model, 
+		Bag<Event<IO_Type,T> >& r);
+		/// Destructor.  Destroys all of the component models.
+		~Digraph();
+
+	private:	
+		// A node in the coupling graph
+		struct node
+		{
+			node():
+			model(NULL),
+			port()
+			{
+			}
+			node(Component* model, PORT port):
+			model(model),
+			port(port)
+			{
+			}
+			const node& operator=(const node& src)
+			{
+				model = src.model;
+				port = src.port;
+				return *this;
+			}
+			Component* model;
+			PORT port;
+		
+			// Comparison for STL map
+			bool operator<(const node& other) const
+			{
+				if (model == other.model) return port < other.port;
+				return model < other.model;
+			}
+		};
+		// Component model set
+		Set<Component*> models;
+		// Coupling information
+		std::map<node,Bag<node> > graph;
+};
+
+template <class VALUE, class PORT, class T>
+void Digraph<VALUE,PORT,T>::add(Component* model)
+{
+	assert(model != this);
+	models.insert(model);
+	model->setParent(this);
+}
+
+template <class VALUE, class PORT, class T>
+void Digraph<VALUE,PORT,T>::couple(Component* src, PORT srcPort, 
+Component* dst, PORT dstPort)
+{
+	if (src != this) add(src);
+	if (dst != this) add(dst);
+	node src_node(src,srcPort);
+	node dst_node(dst,dstPort);
+	graph[src_node].insert(dst_node);
+}
+
+template <class VALUE, class PORT, class T>
+void Digraph<VALUE,PORT,T>::getComponents(Set<Component*>& c)
+{
+	c = models;
+}
+
+template <class VALUE, class PORT, class T>
+void Digraph<VALUE,PORT,T>::
+route(const IO_Type& x, Component* model, 
+Bag<Event<IO_Type,T> >& r)
+{
+	// Find the list of target models and ports
+	node src_node(model,x.port);
+	typename std::map<node,Bag<node> >::iterator graph_iter;
+	graph_iter = graph.find(src_node);
+	// If no target, just return
+	if (graph_iter == graph.end()) return;
+	// Otherwise, add the targets to the event bag
+	Event<IO_Type,T> event;
+	typename Bag<node>::iterator node_iter;
+	for (node_iter = (*graph_iter).second.begin();
+	node_iter != (*graph_iter).second.end(); node_iter++)
+	{
+		event.model = (*node_iter).model;
+		event.value.port = (*node_iter).port;
+		event.value.value = x.value;
+		r.insert(event);
+	}
+}
+template <class VALUE, class PORT, class T>
+Digraph<VALUE,PORT,T>::~Digraph()
+{ 
+	typename Set<Component*>::iterator i;
+	for (i = models.begin(); i != models.end(); i++)
+	{
+		delete *i;
+	}
+}
+
+} // end of namespace 
+
+#endif

simulators/adevs-3.3/include/adevs_event_listener.h

@@ -0,0 +1,74 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef __adevs_event_listener_h_
+#define __adevs_event_listener_h_
+#include "adevs_models.h"
+#include "adevs_bag.h"
+
+namespace adevs
+{
+
+/**
+ * The EventListener interface is used to receive output events produced
+ * by model and to be notified of state changes at Atomic models.
+ */
+template <class X, class T = double> class EventListener
+{
+	public:
+		/**
+		 * This callback is invoked when a model, network or atomic,
+		 * produces an output. The default implementation is empty.
+		 * @param x The model that produced the output and the output's value
+		 * @param t The absolute time at which the output occurred
+		 */
+		virtual void outputEvent(Event<X,T> x, T t){}
+		/**
+		 * This callback is invoked when a model, network or atomic,
+		 * recieves an input. The default implementation is empty.
+		 * @param x The model and input value
+		 * @param t The absolute time at which the output occurred
+		 */
+		virtual void inputEvent(Event<X,T> x, T t){}
+		/**
+		 * This callback is invoked by the simulator after an Atomic
+		 * model changes its state. This method has an empty default
+		 * implementation.
+		 * @param model The model that changed state
+		 * @param t The absolute time at which the state change occurred
+		 */
+		virtual void stateChange(Atomic<X,T>* model, T t){}
+		/// Destructor
+		virtual ~EventListener(){}
+};
+
+} // end of namespace
+
+#endif

simulators/adevs-3.3/include/adevs_event_locators.h

@@ -0,0 +1,233 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef _adevs_event_locators_h_
+#define _adevs_event_locators_h_
+#include "adevs_hybrid.h"
+#include <cmath>
+#include <iostream>
+
+namespace adevs
+{
+
+/**
+ * This is a state event locator that uses either bisection or
+ * linear interpolation to pinpoints events in time.
+ */
+template <typename X> class event_locator_impl:
+	public event_locator<X>
+{
+	public:
+		/**
+		 * The locator will attempt to pinpoint events within err_tol of
+		 * zero for each state event function; i.e., an event occurs
+		 * at the first instant t' >= t where z(t)*z(t') <= 0 and
+		 * |z(t')| < err_tol.
+		 */
+		enum Mode { INTERPOLATE, BISECTION, DISCONTINUOUS };
+		/**
+		 * Create an event locator that will act on the supplied ode system.
+		 * @param sys The system of equations to solve
+		 * @param err_tol Threshold for the zero crossing function that will
+		 * trigger an event detect.
+		 * @param mode The event localization method to be used
+		 */
+		event_locator_impl(ode_system<X>* sys, double err_tol, Mode mode);
+		/// Destructor leaves the ode_system object intact
+		~event_locator_impl();
+		/**
+		 * Location discrete events between now and now + h.
+		 * @param events This array contains one entry for each zero crossing
+		 * function. An entry is true of that function is within err_tol of
+		 * zero and false otherwise.
+		 * @param qstart The continuous variable state vector now.
+		 * @param qend The continuous value state vector when the event
+		 * is detected or at now+h if no event is detected.
+		 * @param solver The ode solver to use when computing continuous
+		 * state vector values.
+		 * @param h The interval of time over which to look for events.
+		 * @return true if at least one event was found and false otherwise.
+		 */
+		bool find_events(bool* events, const double* qstart, double* qend,
+				ode_solver<X>* solver, double& h);
+	private:
+		double *z[2]; // State events at the start and end of [0,h]
+		const double err_tol; // Error tolerance
+		Mode mode;
+
+		int sign(double x) const
+		{
+			if (x < 0.0) return -1;
+			else if (x > 0.0) return 1;
+			else return 0;
+		}
+};
+
+template <typename X>
+event_locator_impl<X>::event_locator_impl(ode_system<X>* sys,
+		double err_tol, Mode mode):
+	event_locator<X>(sys),
+	err_tol(err_tol),
+	mode(mode)
+{
+	z[0] = new double[sys->numEvents()];
+	z[1] = new double[sys->numEvents()];
+}
+
+template <typename X>
+event_locator_impl<X>::~event_locator_impl()
+{
+	delete [] z[0]; delete [] z[1];
+}
+
+template <typename X>
+bool event_locator_impl<X>::find_events(bool* events,
+	const double* qstart, double* qend, ode_solver<X>* solver, double& h)
+{
+	// Calculate the state event functions at the start 
+	// of the interval
+	this->sys->state_event_func(qstart,z[0]);
+	// Look for the first event inside of the interval [0,h]
+	while (this->sys->numEvents() > 0)
+	{
+		double tguess = h;
+		bool event_in_interval = false, found_event = false;
+		this->sys->state_event_func(qend,z[1]);
+		// Do any of the z functions change sign? Have we found an event?
+		for (int i = 0; i < this->sys->numEvents(); i++)
+		{
+			events[i] = false;
+			if (sign(z[1][i]) != sign(z[0][i]))
+			{
+				// Event at h > 0
+				if (
+						// End condition when z is continuous
+						((mode != DISCONTINUOUS) && (fabs(z[1][i]) <= err_tol)) ||
+						// End condition when z is discontinuous
+						((mode == DISCONTINUOUS) && (h <= err_tol))
+					)
+				{
+					events[i] = found_event = true; 
+				}
+				// There is an event in (0,h)
+				else 
+				{
+					if (mode == INTERPOLATE)
+					{
+						double tcandidate = z[0][i]*h/(z[0][i]-z[1][i]);
+						// Don't let the step size go to zero
+						if (tcandidate < h/4.0) tcandidate = h/4.0;
+						if (tcandidate < tguess) tguess = tcandidate;
+					}
+					event_in_interval = true;
+				}
+			}
+		}
+		// Guess at a new h and calculate qend for that time
+		if (event_in_interval)
+		{
+			if (mode == BISECTION || mode == DISCONTINUOUS) h /= 2.0;
+			else h = tguess;
+			for (int i = 0; i < this->sys->numVars(); i++)
+				qend[i] = qstart[i];
+			solver->advance(qend,h);
+		}
+		else return found_event;
+	}
+	// Will never reach this line
+	return false;
+}
+
+/**
+ * Locate events using the bisection method. Your z functions must be
+ * continuous for this to work.
+ */
+template <typename X>
+class bisection_event_locator:
+	public event_locator_impl<X>
+{
+	public:
+		/// Create an event locator that implements the bisection search mode
+		bisection_event_locator(ode_system<X>* sys, double err_tol):
+			event_locator_impl<X>(
+					sys,err_tol,event_locator_impl<X>::BISECTION){}
+};
+
+/**
+ * Locate events using linear interpolation. Your z functions must be
+ * continuous for this to work.
+ */
+template <typename X>
+class linear_event_locator:
+	public event_locator_impl<X>
+{
+	public:
+		/// Create an event locator that implements the linear interpolation search mode
+		linear_event_locator(ode_system<X>* sys, double err_tol):
+			event_locator_impl<X>(
+					sys,err_tol,event_locator_impl<X>::INTERPOLATE){}
+};
+
+/**
+ * Locate events using bisection assuming discontinuous z functions.
+ */
+template <typename X>
+class discontinuous_event_locator:
+	public event_locator_impl<X>
+{
+	public:
+		/**
+		 * Create an event locator that implements a bisection search while
+		 * assuming the z functions are not continuous.
+		 */
+		discontinuous_event_locator(ode_system<X>* sys, double err_tol):
+			event_locator_impl<X>(
+					sys,err_tol,event_locator_impl<X>::DISCONTINUOUS){}
+};
+
+/**
+ * This event locator is for models that have no state events. Its find_events
+ * method simply returns false.
+ */
+template <typename X> class null_event_locator:
+	public event_locator<X>
+{
+	public:
+		null_event_locator():event_locator<X>(NULL){}
+		~null_event_locator(){}
+		bool find_events(bool*, const double*, double*, ode_solver<X>*, double&)
+		{
+			return false;
+		}
+};
+
+} // end of namespace 
+
+#endif

simulators/adevs-3.3/include/adevs_exception.h

@@ -0,0 +1,98 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef _adevs_exception_h_
+#define _adevs_exception_h_
+#include <string>
+#include <exception>
+
+namespace adevs
+{
+
+/**
+ * The adevs::exception class is derived from the standard template
+ * library exception class.
+ */
+class exception: public std::exception
+{
+	public:
+		/**
+		 * Create an exception with an error message and, if appropriate,
+		 * a pointer to the model that created the error.  To avoid
+		 * templated exceptions, the model pointer is just a void*.
+		 */
+		exception(const char* msg, void* model = NULL):
+		std::exception(),
+		msg(msg),
+		model(model) 
+		{}
+		/// Copy constructor.
+		exception(const adevs::exception& src):
+		std::exception(src),
+		msg(src.msg),
+		model(src.model)
+		{}
+		/// Get the error message.
+		const char* what() const throw()
+		{
+			return msg.c_str();
+		}
+		/// Get a pointer to the model that created the error.
+		void* who() const { return model; }
+		/// Destructor.
+		~exception() throw(){}
+	private:
+		std::string msg;
+		void* model;
+};
+
+/**
+ * The unsupported method exception is raised if an optional virtual method
+ * is not supported by a model.
+ */
+class method_not_supported_exception:
+	public exception
+{
+	public:
+		/**
+		 * Constructor should be supplied with the model throwing
+		 * the exception and the name of the method that is not supported.
+		 */
+		method_not_supported_exception(const char* method, void* model):
+			exception((std::string("Unsupported method: ")+std::string(method)).c_str(),
+					model)
+		{
+		}
+};
+
+} // end of namespace
+
+#endif
+

simulators/adevs-3.3/include/adevs_fmi.h

@@ -0,0 +1,592 @@
+/**
+ * Copyright (c) 2014, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef _adevs_fmi_h_
+#define _adevs_fmi_h_
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+#include <cstdlib>
+#include <cstdio>
+#include <cassert>
+#include "adevs_hybrid.h"
+#include "fmi2Functions.h"
+#include "fmi2FunctionTypes.h"
+#include "fmi2TypesPlatform.h"
+
+// Functions for loading DLL and so files
+#ifdef _WIN32
+#include <windows.h>
+#define OPEN_LIB(name) LoadLibrary(name)
+#define GET_FUNC(hndl,name) GetProcAddress(hndl,name)
+#define CLOSE_LIB(hndl) FreeLibrary(hndl)
+#else
+#include <dlfcn.h>
+#define OPEN_LIB(name) dlopen(name,RTLD_LAZY)
+#define GET_FUNC(hndl,name) dlsym(hndl,name)
+#define CLOSE_LIB(hndl) dlclose(hndl)
+#endif
+
+namespace adevs
+{
+
+/**
+ * Load an FMI wrapped continuous system model for use in a
+ * discrete event simulation. The FMI can then be attached
+ * to any of the ODE solvers and event detectors in adevs
+ * for simulation with the Hybrid class. This FMI loader 
+ * does not automatically extract model information from the
+ * description XML, and so that information must be provided
+ * explicitly by the end-user, but you probably need to know
+ * this information regardless if you are using the FMI inside
+ * of a larger discrete event simulation.
+ */
+template <typename X> class FMI:
+	public ode_system<X>
+{
+	public:
+		/**
+		 * This constructs a wrapper around an FMI. The constructor
+		 * must be provided with the FMI's GUID, the number of state variables,
+		 * number of event indicators, and the path to the .so file
+		 * that contains the FMI functions for this model.
+		 */
+		FMI(const char* modelname,
+			const char* guid,
+			const char* resource_location,
+			int num_state_variables,
+			int num_event_indicators,
+			const char* shared_lib_name,
+			const double tolerance = 1E-8,
+			int num_extra_event_indicators = 0,
+			double start_time = 0.0);
+		/// Copy the initial state of the model to q
+		virtual void init(double* q);
+		/// Compute the derivative for state q and put it in dq
+		virtual void der_func(const double* q, double* dq);
+		/// Compute the state event functions for state q and put them in z
+		virtual void state_event_func(const double* q, double* z);
+		/// Compute the time event function using state q
+		virtual double time_event_func(const double* q);
+		/**
+		 * This method is invoked immediately following an update of the
+		 * continuous state variables and signal to the FMI the end
+		 * of an integration state.
+		 */
+		virtual void postStep(double* q);
+		/**
+		 * Like the postStep() method, but this is called at the end
+		 * of a trial step made by the numerical integrator. Trial steps
+		 * are used to locate state events and to estimate numerical
+		 * errors. The state vector q passed to this method may not be
+		 * the final vector assigned to the model at the end of
+		 * the current integration step. The get that value use
+		 * the postStep() method.
+		 */
+		virtual void postTrialStep(double* q);
+		/**
+		 * The internal transition function. This function will process all events
+		 * required by the FMI. Any derived class should call this method for the
+		 * parent class, then set or get any variables as appropriate, and then
+		 * call the base class method again to account for these changes. 
+		 */
+		virtual void internal_event(double* q,
+				const bool* state_event);
+		/**
+		 * The external transition See the notes on the internal_event function for
+		 * derived classes.
+		 */
+		virtual void external_event(double* q, double e,
+				const Bag<X>& xb);
+		/**
+		 * The confluent transition function. See the notes on the internal_event function for
+		 * derived classes.
+		 */
+		virtual void confluent_event(double *q, const bool* state_event,
+				const Bag<X>& xb);
+		/**
+		 * The output function. This can read variables for the FMI, but should
+		 * not make any modifications to those variables.
+		 */
+		virtual void output_func(const double *q, const bool* state_event,
+				Bag<X>& yb);
+		/**
+		 * Garbage collection function. This works just like the Atomic gc_output method.
+		 * The default implementation does nothing.
+		 */
+		virtual void gc_output(Bag<X>& gb);
+		/// Destructor
+		virtual ~FMI();
+		/// Get the current time
+		double get_time() const { return t_now; }
+		/// Get the value of a real variable
+		double get_real(int k);
+		/// Set the value of a real variable
+		void set_real(int k, double val);
+		/// Get the value of an integer variable
+		int get_int(int k);
+		/// Set the value of an integer variable
+		void set_int(int k, int val);
+		/// Get the value of a boolean variable
+		bool get_bool(int k);
+		/// Set the value of a boolean variable
+		void set_bool(int k, bool val);
+		/// Get the value of a string variable
+		std::string get_string(int k);
+		/// Set the value of a string variable
+		void set_string(int k, std::string& val);
+	private:
+		// Reference to the FMI
+		fmi2Component c;
+		// Pointer to the FMI interface
+		fmi2Component (*_fmi2Instantiate)(fmi2String, fmi2Type,
+				fmi2String, fmi2String, const fmi2CallbackFunctions*,
+				fmi2Boolean, fmi2Boolean);
+		void (*_fmi2FreeInstance)(fmi2Component);
+		fmi2Status (*_fmi2SetupExperiment)(fmi2Component, fmi2Boolean,
+				fmi2Real, fmi2Real, fmi2Boolean, fmi2Real);
+		fmi2Status (*_fmi2EnterInitializationMode)(fmi2Component);
+		fmi2Status (*_fmi2ExitInitializationMode)(fmi2Component);
+		fmi2Status (*_fmi2GetReal)(fmi2Component, const fmi2ValueReference*, size_t, fmi2Real*);
+		fmi2Status (*_fmi2GetInteger)(fmi2Component, const fmi2ValueReference*, size_t, fmi2Integer*);
+		fmi2Status (*_fmi2GetBoolean)(fmi2Component, const fmi2ValueReference*, size_t, fmi2Boolean*);
+		fmi2Status (*_fmi2GetString)(fmi2Component, const fmi2ValueReference*, size_t, fmi2String*);
+		fmi2Status (*_fmi2SetReal)(fmi2Component, const fmi2ValueReference*, size_t, const fmi2Real*);
+		fmi2Status (*_fmi2SetInteger)(fmi2Component, const fmi2ValueReference*, size_t, const fmi2Integer*);
+		fmi2Status (*_fmi2SetBoolean)(fmi2Component, const fmi2ValueReference*, size_t, const fmi2Boolean*);
+		fmi2Status (*_fmi2SetString)(fmi2Component, const fmi2ValueReference*, size_t, const fmi2String*);
+		fmi2Status (*_fmi2EnterEventMode)(fmi2Component);
+		fmi2Status (*_fmi2NewDiscreteStates)(fmi2Component,fmi2EventInfo*);
+		fmi2Status (*_fmi2EnterContinuousTimeMode)(fmi2Component);
+		fmi2Status (*_fmi2CompletedIntegratorStep)(fmi2Component, fmi2Boolean, fmi2Boolean*, fmi2Boolean*);
+		fmi2Status (*_fmi2SetTime)(fmi2Component, fmi2Real);
+		fmi2Status (*_fmi2SetContinuousStates)(fmi2Component, const fmi2Real*, size_t);
+		fmi2Status (*_fmi2GetDerivatives)(fmi2Component, fmi2Real*, size_t);
+		fmi2Status (*_fmi2GetEventIndicators)(fmi2Component, fmi2Real*, size_t);
+		fmi2Status (*_fmi2GetContinuousStates)(fmi2Component, fmi2Real*, size_t);	
+		
+		// Instant of the next time event
+		double next_time_event;
+		// Current time
+		double t_now;
+		// so library handle
+		#ifdef _WIN32
+		HINSTANCE so_hndl;
+		#else
+		void* so_hndl;
+		#endif
+		// Are we in continuous time mode?
+		bool cont_time_mode;
+		// Number of event indicators that are not governed by the FMI
+		int num_extra_event_indicators;
+		// Start time of the simulation
+		double start_time;
+
+		static void fmilogger(
+			fmi2ComponentEnvironment componentEnvironment,
+			fmi2String instanceName,
+			fmi2Status status,
+			fmi2String category,
+			fmi2String message,...)
+		{
+			if (message != NULL){			
+		
+				fprintf(stderr, message,"\n");
+			}
+		}
+
+		fmi2CallbackFunctions* callbackFuncs;
+
+		void iterate_events();
+};
+
+template <typename X>
+FMI<X>::FMI(const char* modelname,
+			const char* guid,
+			const char* resource_location,
+			int num_state_variables,
+			int num_event_indicators,
+			const char* so_file_name,
+			const double tolerance,
+			int num_extra_event_indicators,
+			double start_time):
+	// One extra variable at the end for time
+	ode_system<X>(num_state_variables+1,num_event_indicators+num_extra_event_indicators),
+	next_time_event(adevs_inf<double>()),
+	t_now(start_time),
+	so_hndl(NULL),
+	cont_time_mode(false),
+	num_extra_event_indicators(num_extra_event_indicators)	
+{
+	fmi2CallbackFunctions tmp = {adevs::FMI<X>::fmilogger,calloc,free,NULL,NULL};
+	callbackFuncs = new fmi2CallbackFunctions(tmp);
+	so_hndl = OPEN_LIB(so_file_name);
+	if (so_hndl == NULL)
+	{
+		throw adevs::exception("Could not load so file",this);
+    }
+	// This only works with a POSIX compliant compiler/system
+	_fmi2Instantiate = (fmi2Component (*)(fmi2String, fmi2Type,
+		fmi2String, fmi2String, const fmi2CallbackFunctions*,
+		fmi2Boolean, fmi2Boolean))GET_FUNC(so_hndl,"fmi2Instantiate");
+	assert(_fmi2Instantiate != NULL);
+	_fmi2FreeInstance = (void (*)(fmi2Component))GET_FUNC(so_hndl,"fmi2FreeInstance");
+	assert(_fmi2FreeInstance != NULL);
+	_fmi2SetupExperiment = (fmi2Status (*)(fmi2Component, fmi2Boolean,
+		fmi2Real, fmi2Real, fmi2Boolean, fmi2Real))GET_FUNC(so_hndl,"fmi2SetupExperiment");
+	assert(_fmi2SetupExperiment != NULL);
+	_fmi2EnterInitializationMode = (fmi2Status (*)(fmi2Component))GET_FUNC(so_hndl,"fmi2EnterInitializationMode");
+	assert(_fmi2EnterInitializationMode != NULL);
+	_fmi2ExitInitializationMode = (fmi2Status (*)(fmi2Component))GET_FUNC(so_hndl,"fmi2ExitInitializationMode");
+	assert(_fmi2ExitInitializationMode != NULL);
+	_fmi2GetReal = (fmi2Status (*)(fmi2Component, const fmi2ValueReference*, size_t, fmi2Real*))
+		GET_FUNC(so_hndl,"fmi2GetReal");
+	assert(_fmi2GetReal != NULL);
+	_fmi2GetInteger = (fmi2Status (*)(fmi2Component, const fmi2ValueReference*, size_t, fmi2Integer*)) 
+		GET_FUNC(so_hndl,"fmi2GetInteger");
+	assert(_fmi2GetInteger != NULL);
+	_fmi2GetBoolean = (fmi2Status (*)(fmi2Component, const fmi2ValueReference*, size_t, fmi2Boolean*))
+		GET_FUNC(so_hndl,"fmi2GetBoolean");
+	assert(_fmi2GetBoolean != NULL);
+	_fmi2GetString = (fmi2Status (*)(fmi2Component, const fmi2ValueReference*, size_t, fmi2String*))
+		GET_FUNC(so_hndl,"fmi2GetString");
+	assert(_fmi2GetString != NULL);
+	_fmi2SetReal = (fmi2Status (*)(fmi2Component, const fmi2ValueReference*, size_t, const fmi2Real*))
+		GET_FUNC(so_hndl,"fmi2SetReal");
+	assert(_fmi2SetReal != NULL);
+	_fmi2SetInteger = (fmi2Status (*)(fmi2Component, const fmi2ValueReference*, size_t, const fmi2Integer*))
+		GET_FUNC(so_hndl,"fmi2SetInteger");
+	assert(_fmi2SetInteger != NULL);
+	_fmi2SetBoolean = (fmi2Status (*)(fmi2Component, const fmi2ValueReference*, size_t, const fmi2Boolean*))
+		GET_FUNC(so_hndl,"fmi2SetBoolean");
+	assert(_fmi2SetBoolean != NULL);
+	_fmi2SetString = (fmi2Status (*)(fmi2Component, const fmi2ValueReference*, size_t, const fmi2String*))
+		GET_FUNC(so_hndl,"fmi2SetString");
+	assert(_fmi2SetString != NULL);
+	_fmi2EnterEventMode = (fmi2Status (*)(fmi2Component))GET_FUNC(so_hndl,"fmi2EnterEventMode");
+	assert(_fmi2EnterEventMode != NULL);
+	_fmi2NewDiscreteStates = (fmi2Status (*)(fmi2Component,fmi2EventInfo*))GET_FUNC(so_hndl,"fmi2NewDiscreteStates");
+	assert(_fmi2NewDiscreteStates != NULL);
+	_fmi2EnterContinuousTimeMode = (fmi2Status (*)(fmi2Component))GET_FUNC(so_hndl,"fmi2EnterContinuousTimeMode");
+	assert(_fmi2EnterContinuousTimeMode != NULL);
+	_fmi2CompletedIntegratorStep = (fmi2Status (*)(fmi2Component, fmi2Boolean, fmi2Boolean*, fmi2Boolean*))
+		GET_FUNC(so_hndl,"fmi2CompletedIntegratorStep");
+	assert(_fmi2CompletedIntegratorStep != NULL);
+	_fmi2SetTime = (fmi2Status (*)(fmi2Component, fmi2Real))GET_FUNC(so_hndl,"fmi2SetTime");
+	assert(_fmi2SetTime != NULL);
+	_fmi2SetContinuousStates = (fmi2Status (*)(fmi2Component, const fmi2Real*, size_t))
+		GET_FUNC(so_hndl,"fmi2SetContinuousStates");
+	assert(_fmi2SetContinuousStates != NULL);
+	_fmi2GetDerivatives = (fmi2Status (*)(fmi2Component, fmi2Real*, size_t))GET_FUNC(so_hndl,"fmi2GetDerivatives");
+	assert(_fmi2GetDerivatives != NULL);
+	_fmi2GetEventIndicators = (fmi2Status (*)(fmi2Component, fmi2Real*, size_t))GET_FUNC(so_hndl,"fmi2GetEventIndicators");
+	assert(_fmi2GetEventIndicators != NULL);
+	_fmi2GetContinuousStates = (fmi2Status (*)(fmi2Component, fmi2Real*, size_t))GET_FUNC(so_hndl,"fmi2GetContinuousStates");
+	assert(_fmi2GetContinuousStates != NULL);
+	// Create the FMI component
+	c = _fmi2Instantiate(modelname,fmi2ModelExchange,guid,resource_location,callbackFuncs,fmi2False,fmi2False);
+	assert(c != NULL);
+	_fmi2SetupExperiment(c,fmi2True,tolerance,-1.0,fmi2False,-1.0);
+}
+
+template <typename X>
+void FMI<X>::iterate_events()
+{
+	fmi2Status status;
+	// Put into consistent initial state
+	fmi2EventInfo eventInfo;
+	do
+	{
+		status = _fmi2NewDiscreteStates(c,&eventInfo);
+		assert(status == fmi2OK);
+	}
+	while (eventInfo.newDiscreteStatesNeeded == fmi2True);
+	if (eventInfo.nextEventTimeDefined == fmi2True)
+		next_time_event = eventInfo.nextEventTime;
+	else
+		next_time_event = adevs_inf<double>();
+	assert(status == fmi2OK);
+}
+
+template <typename X>
+void FMI<X>::init(double* q)
+{
+	fmi2Status status;
+	// Initialize all variables
+	status = _fmi2EnterInitializationMode(c);
+	assert(status == fmi2OK);
+	// Done with initialization
+	status = _fmi2ExitInitializationMode(c);
+	assert(status == fmi2OK);
+	// Put into consistent initial state
+	iterate_events();
+	// Enter continuous time mode to start integration
+	status = _fmi2EnterContinuousTimeMode(c);
+	assert(status == fmi2OK);
+	// Set initial value for time
+	status = _fmi2SetTime(c,t_now);
+	assert(status == fmi2OK);
+	// Get starting state variable values
+	status = _fmi2GetContinuousStates(c,q,this->numVars()-1);
+	assert(status == fmi2OK);
+	q[this->numVars()-1] = t_now;
+	cont_time_mode = true;
+}
+
+template <typename X>
+void FMI<X>::der_func(const double* q, double* dq)
+{
+	fmi2Status status;
+	if (!cont_time_mode)
+	{
+		status = _fmi2EnterContinuousTimeMode(c);
+		assert(status == fmi2OK);
+		cont_time_mode = true;
+	}
+	status =_fmi2SetTime(c,q[this->numVars()-1]);
+	assert(status == fmi2OK);
+	status = _fmi2SetContinuousStates(c,q,this->numVars()-1);
+	assert(status == fmi2OK);
+	status = _fmi2GetDerivatives(c,dq,this->numVars()-1);
+	assert(status == fmi2OK);
+	dq[this->numVars()-1] = 1.0;
+}
+
+template <typename X>
+void FMI<X>::state_event_func(const double* q, double* z)
+{
+	fmi2Status status;
+	if (!cont_time_mode)
+	{
+		status = _fmi2EnterContinuousTimeMode(c);
+		assert(status == fmi2OK);
+		cont_time_mode = true;
+	}
+	status = _fmi2SetTime(c,q[this->numVars()-1]);
+	assert(status == fmi2OK);
+	status = _fmi2SetContinuousStates(c,q,this->numVars()-1);
+	assert(status == fmi2OK);
+	status = _fmi2GetEventIndicators(c,z,this->numEvents()-num_extra_event_indicators);
+	assert(status == fmi2OK);
+}
+
+template <typename X>
+double FMI<X>::time_event_func(const double* q)
+{
+	return next_time_event-q[this->numVars()-1];
+}
+
+template <typename X>
+void FMI<X>::postStep(double* q)
+{
+	assert(cont_time_mode);
+	// Don't advance the FMI state by zero units of time
+	// when in continuous mode
+	if (q[this->numVars()-1] <= t_now)
+		return;
+	// Try to complete the integration step
+	fmi2Status status;
+	fmi2Boolean enterEventMode;
+	fmi2Boolean terminateSimulation;
+	t_now = q[this->numVars()-1];
+	status = _fmi2SetTime(c,t_now);
+	assert(status == fmi2OK);
+	status = _fmi2SetContinuousStates(c,q,this->numVars()-1);
+	assert(status == fmi2OK);
+	status = _fmi2CompletedIntegratorStep(c,fmi2True,&enterEventMode,&terminateSimulation);
+	assert(status == fmi2OK);
+	// Force an event if one is indicated
+	if (enterEventMode == fmi2True)
+		next_time_event = t_now;
+}
+
+template <typename X>
+void FMI<X>::postTrialStep(double* q)
+{
+	assert(cont_time_mode);
+	// Restore values changed by der_func and state_event_func
+	fmi2Status status;
+	status = _fmi2SetTime(c,q[this->numVars()-1]);
+	assert(status == fmi2OK);
+	status = _fmi2SetContinuousStates(c,q,this->numVars()-1);
+	assert(status == fmi2OK);
+}
+
+template <typename X>
+void FMI<X>::internal_event(double* q, const bool* state_event)
+{
+	fmi2Status status;
+	// postStep will have updated the continuous variables, so 
+	// we just process discrete events here.
+	if (cont_time_mode)
+	{
+		status = _fmi2EnterEventMode(c);
+		assert(status == fmi2OK);
+		cont_time_mode = false;
+	}
+	// Process events
+	iterate_events();
+	// Update the state variable array
+	status = _fmi2GetContinuousStates(c,q,this->numVars()-1);
+	assert(status == fmi2OK);
+}
+				
+template <typename X>
+void FMI<X>::external_event(double* q, double e, const Bag<X>& xb)
+{
+	fmi2Status status;
+	// Go to event mode if we have not yet done so
+	if (cont_time_mode)
+	{
+		status = _fmi2EnterEventMode(c);
+		assert(status == fmi2OK);
+		cont_time_mode = false;
+	}
+	// process any events that need processing
+	iterate_events();
+	status = _fmi2GetContinuousStates(c,q,this->numVars()-1);
+	assert(status == fmi2OK);
+}
+				
+template <typename X>
+void FMI<X>::confluent_event(double *q, const bool* state_event, const Bag<X>& xb)
+{
+	fmi2Status status;
+	// postStep will have updated the continuous variables, so 
+	// we just process discrete events here.
+	if (cont_time_mode)
+	{
+		status = _fmi2EnterEventMode(c);
+		assert(status == fmi2OK);
+		cont_time_mode = false;
+	}
+	iterate_events();
+	status = _fmi2GetContinuousStates(c,q,this->numVars()-1);
+	assert(status == fmi2OK);
+}
+				
+template <typename X>
+void FMI<X>::output_func(const double *q, const bool* state_event, Bag<X>& yb)
+{
+}
+
+template <typename X>
+void FMI<X>::gc_output(Bag<X>& gb)
+{
+}
+
+template <typename X>
+FMI<X>::~FMI()
+{
+	_fmi2FreeInstance(c);
+	delete callbackFuncs;
+	CLOSE_LIB(so_hndl);
+}
+
+template <typename X>
+double FMI<X>::get_real(int k)
+{
+	const fmi2ValueReference ref = k;
+	fmi2Real val;
+	fmi2Status status = _fmi2GetReal(c,&ref,1,&val);
+	assert(status == fmi2OK);
+	return val;
+}
+
+template <typename X>
+void FMI<X>::set_real(int k, double val)
+{
+	const fmi2ValueReference ref = k;
+	fmi2Real fmi_val = val;
+	fmi2Status status = _fmi2SetReal(c,&ref,1,&fmi_val);
+	assert(status == fmi2OK);
+}
+
+template <typename X>
+int FMI<X>::get_int(int k)
+{
+	const fmi2ValueReference ref = k;
+	fmi2Integer val;
+	fmi2Status status = _fmi2GetInteger(c,&ref,1,&val);
+	assert(status == fmi2OK);
+	return val;
+}
+
+template <typename X>
+void FMI<X>::set_int(int k, int val)
+{
+	const fmi2ValueReference ref = k;
+	fmi2Integer fmi_val = val;
+	fmi2Status status = _fmi2SetInteger(c,&ref,1,&fmi_val);
+	assert(status == fmi2OK);
+}
+
+template <typename X>
+bool FMI<X>::get_bool(int k)
+{
+	const fmi2ValueReference ref = k;
+	fmi2Boolean val;
+	fmi2Status status = _fmi2GetBoolean(c,&ref,1,&val);
+	assert(status == fmi2OK);
+	return (val == fmi2True);
+}
+
+template <typename X>
+void FMI<X>::set_bool(int k, bool val)
+{
+	const fmi2ValueReference ref = k;
+	fmi2Boolean fmi_val = fmi2False;
+	if (val) fmi_val = fmi2True;
+	fmi2Status status = _fmi2SetBoolean(c,&ref,1,&fmi_val);
+	assert(status == fmi2OK);
+}
+
+template <typename X>
+std::string FMI<X>::get_string(int k)
+{
+	const fmi2ValueReference ref = k;
+	fmi2String val;
+	fmi2Status status = _fmi2GetString(c,&ref,1,&val);
+	assert(status == fmi2OK);
+	return val;
+}
+
+template <typename X>
+void FMI<X>::set_string(int k, std::string& val)
+{
+	const fmi2ValueReference ref = k;
+	fmi2String fmi_val = fmi2False;
+	fmi2Status status = _fmi2SetString(c,&ref,1,&fmi_val);
+	assert(status == fmi2OK);
+}
+
+} // end of namespace
+
+#endif

simulators/adevs-3.3/include/adevs_hybrid.h

@@ -0,0 +1,617 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef _adevs_hybrid_h_
+#define _adevs_hybrid_h_
+#include <algorithm>
+#include <cmath>
+#include "adevs_models.h"
+
+namespace adevs
+{
+
+/**
+ * This is the base class for all hybrid systems. The init and der_func methods
+ * are used to implement the model's continuous dynamics. The other functions
+ * are for discrete event dynamics.
+ */
+template <typename X> class ode_system
+{
+	public:
+		/// Make a system with N state variables and M state event functions
+		ode_system(int N_vars, int M_event_funcs):
+			N(N_vars),M(M_event_funcs){}
+		/// Get the number of state variables
+		int numVars() const { return N; }
+		/// Get the number of state events
+		int numEvents() const { return M; }
+		/// Copy the initial state of the model to q
+		virtual void init(double* q) = 0;
+		/// Compute the derivative for state q and put it in dq
+		virtual void der_func(const double* q, double* dq) = 0;
+		/// Compute the state event functions for state q and put them in z
+		virtual void state_event_func(const double* q, double* z) = 0;
+		/// Compute the time event function using state q
+		virtual double time_event_func(const double* q) = 0;
+		/**
+		 * This method is invoked immediately following an update of the
+		 * continuous state variables. The main use of this callback is to
+		 * update algberaic variables. The default implementation does nothing.
+		 */
+		virtual void postStep(double* q){};
+		/**
+		 * This is called after a trial step. It can be used to restore the values
+		 * of any variables that might have been changed by der_func or state_event_func
+		 * while calculating the trial step. By default this method does nothing.
+		 */
+		virtual void postTrialStep(double* q){};
+		/// The internal transition function
+		virtual void internal_event(double* q,
+				const bool* state_event) = 0;
+		/// The external transition function
+		virtual void external_event(double* q, double e,
+				const Bag<X>& xb) = 0;
+		/// The confluent transition function
+		virtual void confluent_event(double *q, const bool* state_event,
+				const Bag<X>& xb) = 0;
+		/// The output function
+		virtual void output_func(const double *q, const bool* state_event,
+				Bag<X>& yb) = 0;
+		/// Garbage collection function. This works just like the Atomic gc_output method.
+		virtual void gc_output(Bag<X>& gb) = 0;
+		/// Destructor
+		virtual ~ode_system(){}
+	private:
+		const int N, M;
+};
+
+/**
+ * <p>This extension of the ode_system provides for modeling some semi-explicit
+ * DAEs of index 1, specifically those in the form dx/dt = f(x,y), y = g(x,y).
+ * The solution to y=g(x,y) is found by iteration on y.
+ * See "The Numerical Solution of Differential-Algebraic Systems by Runge-Kutta
+ * Methods" by Ernst Hairer, Michel Roche and Christian Lubich, published
+ * by Springer as Lecture Notes in Mathematics, Volum 1409, c. 1989.
+ * The section on Half-explicit methods (starting pg. 20 of my copy)
+ * describes the procedure.</p>
+ * <p>Only the methods that include the algebraic variables should be overriden.
+ * Any explicit, single step ODE solver can be used to generate trajectories for
+ * this object (e.g., the Runge-Kutta methods included with adevs will work).</p>
+ */
+template <typename X> class dae_se1_system:
+	public ode_system<X>
+{
+	public:
+		/**
+		 * Make a system with N state variables, M state event functions
+		 * and A algebraic variables. The error tolerance determines
+		 * the acceptable accuracy for the algebraic solution to y=g(x,y).
+		 * The max_iters argument determines how many iterations will
+		 * be used to try and generate a solution.
+		 */
+		dae_se1_system(int N_vars, int M_event_funcs, int A_alg_vars,
+				double err_tol = 1E-10, int max_iters = 30, double alpha = -1.0):
+			ode_system<X>(N_vars,M_event_funcs),
+			A(A_alg_vars),
+			max_iters(max_iters),
+			err_tol(err_tol),
+			alpha(alpha)
+			{
+				failed = 0;
+				max_err = 0.0;
+				a = new double[A];
+				atmp = new double[A];
+				d = new double[A];
+				f[1] = new double[A];
+				f[0] = new double[A];
+			}
+		/// Get the number of algebraic variables
+		int numAlgVars() const { return A; }
+		/// Get the algebraic variables
+		double getAlgVar(int i) const { return a[i]; }
+		/**
+		 * Write an intial solution for the state variables q
+		 * and algebraic variables a.
+		 */
+		virtual void init(double* q, double* a) = 0;
+		/**
+		 * Calculate the algebraic function for the state vector
+		 * q and algebraic variables a and store the result to af.
+		 * A solution to alg_func(a,q) = a will be found by 
+		 * iteration on this function.
+		 */
+		virtual void alg_func(const double* q, const double* a, double* af) = 0;
+		/**
+		 * Calculate the derivative of the state variables and store
+		 * the result in dq.
+		 */
+		virtual void der_func(const double* q, const double* a, double* dq) = 0;
+		/**
+		 * Calculate the state event functions and store the result in z.
+		 */
+		virtual void state_event_func(const double* q, const double* a, double* z) = 0;
+		/// Compute the time event function using state q and algebraic variables a
+		virtual double time_event_func(const double* q, const double* a) = 0;
+		/**
+		 * Update any variables that need updating at then end of a simulation step.
+		 * This is called both when postStep(q) would be called and also immediately
+		 * after the execution a discrete state transition.
+		 */
+		virtual void postStep(double* q, double* a) = 0;
+		/**
+		 * Default implementation calls postTrialStep(q)
+		 */
+		virtual void postTrialStep(double *q, double* a) { ode_system<X>::postTrialStep(q); }
+		/// The internal transition function
+		virtual void internal_event(double* q, double* a,
+				const bool* state_event) = 0;
+		/// The external transition function
+		virtual void external_event(double* q, double* a,
+				double e, const Bag<X>& xb) = 0;
+		/// The confluent transition function
+		virtual void confluent_event(double *q, double* a, 
+				const bool* state_event, const Bag<X>& xb) = 0;
+		/// The output function
+		virtual void output_func(const double *q, const double* a,
+				const bool* state_event, Bag<X>& yb) = 0;
+		/// Destructor
+		virtual ~dae_se1_system()
+		{
+			delete [] d;
+			delete [] a;
+			delete [] atmp;
+			delete [] f[1];
+			delete [] f[0];
+		}
+		/**
+		 * Get the number of times that the error tolerance was not satisfied
+		 * before the iteration limit was reached.
+		 */
+		int getIterFailCount() const { return failed; }
+		/**
+		 * Get the worst error in a case where the algebraic solver did
+		 * not satisfy the error tolerance. This will be zero if
+		 * there were no failures of the algebraic solver.
+		 */
+		double getWorseError() const { return max_err; }
+		/// Do not override
+		void init(double* q)
+		{
+			init(q,a);
+		}
+		/// Do not override
+		void der_func(const double* q, double* dq)
+		{
+			solve(q);
+			der_func(q,a,dq);
+		}
+		/// Override only if you have no state event functions.
+		void state_event_func(const double* q, double* z)
+		{
+			solve(q);
+			state_event_func(q,a,z);
+		}
+		/// Override only if you have no time events.
+		double time_event_func(const double* q)
+		{
+			solve(q);
+			return time_event_func(q,a);
+		}
+		/// Do not override
+		void postStep(double* q)
+		{
+			solve(q);
+			postStep(q,a);
+		}
+		/// Do not override
+		void internal_event(double* q, const bool* state_event)
+		{
+			// The variable a was solved for in the post step
+			internal_event(q,a,state_event);
+			// Make sure the algebraic variables are consistent with q
+			solve(q);
+			postStep(q,a);
+		}
+		/// Do not override
+		void external_event(double* q, double e, const Bag<X>& xb)
+		{
+			// The variable a was solved for in the post step
+			external_event(q,a,e,xb);
+			// Make sure the algebraic variables are consistent with q
+			solve(q);
+			postStep(q,a);
+		}
+		/// Do not override
+		void confluent_event(double *q, const bool* state_event, const Bag<X>& xb)
+		{
+			// The variable a was solved for in the post step
+			confluent_event(q,a,state_event,xb);
+			// Make sure the algebraic variables are consistent with q
+			solve(q);
+			postStep(q,a);
+		}
+		/// Do not override
+		void output_func(const double *q, const bool* state_event, Bag<X>& yb)
+		{
+			// The variable a was solved for in the post step
+			output_func(q,a,state_event,yb);
+		}
+	protected:
+		/**
+		 * Solve the algebraic equations. This should
+		 * not usually need to be called by the derived class. An exception might
+		 * be where updated values for the algebraic variables are needed from
+		 * within an event function due to some discrete change in q or the
+		 * structure of the systems.
+		 */
+		void solve(const double* q);
+	private:
+		const int A, max_iters;
+		const double err_tol, alpha;
+		// Guess at the algebraic solution
+		double *a, *atmp;
+		// Guesses at g(y)-y
+		double* f[2];
+		// Direction
+		double* d;
+		// Maximum error in the wake of a failure
+		double max_err;
+		// Number of failures
+		int failed;
+};
+
+template <typename X>
+void dae_se1_system<X>::solve(const double* q)
+{
+	int iter_count = 0, alt, good;
+	double prev_err, err = 0.0, ee, beta, g2, alpha_tmp = alpha;
+	/**
+	 * Solve y=g(x,y) by the conjugate gradient method.
+	 * Method iterates on f(x,y)=g(x,y)-y to find
+	 * f(x,y)=0.
+	 */
+	_adevs_dae_se_1_system_solve_try_it_again:
+	alt = 0;
+	good = 1;
+	prev_err = DBL_MAX;
+	// First step by steepest descent
+	alg_func(q,a,f[alt]);
+	for (int i = 0; i < A; i++)
+	{
+		// Calculate f(x,y)
+		f[alt][i] -= a[i];
+		// First direction
+		d[i] = -f[alt][i];
+		// Make the move
+		atmp[i] = a[i];
+		a[i] += alpha_tmp*d[i];
+	}
+	// Otherwise, first guess by steepest descent
+	// Finish search by conjugate gradiant
+	while (iter_count < max_iters)
+	{
+		iter_count++;
+		err = 0.0;
+		// Calculate y = g(x,y) 
+		alg_func(q,a,f[good]);
+		// Check the quality of the solution
+		for (int i = 0; i < A; i++)
+		{
+			// Calculate f(x,y) 
+			f[good][i] -= a[i];
+			// Get the largest error 
+			ee = fabs(f[good][i]);
+			if (ee > err) err = ee;
+		}
+		// If the solution is good enough then return
+		if (err < err_tol) return;
+		// If the solution is not converging...
+		if (err > prev_err)
+		{
+			// Restore previous solution
+			for (int i = 0; i < A; i++)
+				a[i] = atmp[i];
+			// Restart with a new value for alpha
+			if (alpha_tmp < 0.0) alpha_tmp = -alpha_tmp;
+			else alpha_tmp *= -0.5;
+			goto _adevs_dae_se_1_system_solve_try_it_again;
+		}
+		prev_err = err;
+		// Calculate beta. See Strang's "Intro. to Applied Mathematics",
+		// pg. 379.
+		beta = g2 = 0.0;
+		for (int i = 0; i < A; i++)
+			g2 += f[alt][i]*f[alt][i];
+		for (int i = 0; i < A; i++)
+			beta += f[good][i]*(f[good][i]-f[alt][i]);
+		beta /= g2;
+		// Calculate a new guess at the solution
+		for (int i = 0; i < A; i++)
+		{
+			d[i] = beta*d[i]-f[good][i];
+			atmp[i] = a[i];
+			a[i] += alpha_tmp*d[i];
+		}
+		// Swap buffers
+		good = alt;
+		alt = (good+1)%2;
+	}
+	// Increment the failed count and worse case error if an 
+	// acceptible solution was not found.
+	failed++;
+	if (err > max_err)
+		max_err = err;
+}
+
+/**
+ * This is the interface for numerical integrators that are to be used with the
+ * Hybrid class.
+ */
+template <typename X> class ode_solver
+{
+	public:
+		/**
+		 * Create and ode_solver that will integrate the der_func method of the
+		 * supplied ode_system.
+		 */
+		ode_solver(ode_system<X>* sys):sys(sys){}
+		/**
+		 * Take an integration step from state q of at most size h_lim and
+		 * return the step size that was actually used. Copy the result of
+		 * the integration step to q.
+		 */
+		virtual double integrate(double* q, double h_lim) = 0;
+		/**
+		 * Advance the system through exactly h units of time.
+		 */
+		virtual void advance(double* q, double h) = 0;
+		/// Destructor
+		virtual ~ode_solver(){}
+	protected:
+		/// The system of odes to be acted upon by the solver
+		ode_system<X>* sys;
+};
+
+/**
+ * This is the interface for algorithms that detect state events in the trajectory
+ * of an ode_system. The ode_solver provided to this class is used to compute
+ * intermediate states during the detection process.
+ */
+template <typename X> class event_locator
+{
+	public:
+		/**
+		 * The locator will use the der_func and state_event_func of the supplied
+		 * ode_system object.
+		 */
+		event_locator(ode_system<X>* sys):sys(sys){}
+		/**
+		 * Find the first state event in the interval [0,h] starting from
+		 * state qstart. The method returns true if an event is found,
+		 * setting the events flags to true if the corresponding z entry in
+		 * the state_event_func above triggered the event. The value of
+		 * h is overwritten with the event time, and the state of the model
+		 * at that time is copied to qend. The event finding method should
+		 * select an instant of time when the zero crossing function is zero or
+		 * has changed sign to trigger an event.
+		 */
+		virtual bool find_events(bool* events, const double *qstart, 
+				double* qend, ode_solver<X>* solver, double& h) = 0;
+		/// Destructor
+		virtual ~event_locator(){}
+	protected:
+		/// The system of odes to be acted upon by the event locator
+		ode_system<X>* sys;
+};
+
+/**
+ * This Atomic model encapsulates an ode_system and numerical solvers for it.
+ * Output from the Hybrid model is produced by the output_func method of the
+ * ode_system whenever a state event or time event occurs. Internal, external,
+ * and confluent events for the Hybrid model are computed with the corresponding
+ * methods of the ode_system. The time advance of the Hybrid class ensures that
+ * its internal events coincide with state and time events in the ode_system.
+ */
+template <typename X, class T = double> class Hybrid:
+	public Atomic<X,T>
+{
+	public:
+		/**
+		 * Create and initialize a simulator for the system. All objects 
+		 * are adopted by the Hybrid object and are deleted when it is.
+		 */
+		Hybrid(ode_system<X>* sys, ode_solver<X>* solver,
+				event_locator<X>* event_finder):
+			sys(sys),solver(solver),event_finder(event_finder),
+			e_accum(0.0)
+		{
+			q = new double[sys->numVars()];
+			q_trial = new double[sys->numVars()];
+			event = new bool[sys->numEvents()+1];
+			event_exists = false;
+			sys->init(q_trial); // Get the initial state of the model
+			for (int i = 0; i < sys->numVars(); i++) q[i] = q_trial[i];
+			tentative_step(); // Take the first tentative step
+		}
+		/// Get the value of the kth continuous state variable
+		double getState(int k) const { return q[k]; }
+		/// Get the array of state variables
+		const double* getState() const { return q; }
+		/// Get the system that this solver is operating on
+		ode_system<X>* getSystem() { return sys; }
+		/// Did a discrete event occur at the last state transition?
+		bool eventHappened() const { return event_happened; }
+		/**
+		 * Do not override this method. It performs numerical integration and
+		 * invokes the ode_system method for internal events as needed.
+		 */
+		void delta_int()
+		{
+			if (!missedOutput.empty())
+			{
+				missedOutput.clear();
+				return;
+			}
+			e_accum += ta();
+			// Execute any discrete events
+			event_happened = event_exists;
+			if (event_exists) // Execute the internal event
+			{
+				sys->internal_event(q_trial,event); 
+				e_accum = 0.0;
+			}
+			// Copy the new state vector to q
+			for (int i = 0; i < sys->numVars(); i++) q[i] = q_trial[i];
+			tentative_step(); // Take a tentative step
+		}
+		/**
+		 * Do not override this method. It performs numerical integration and
+		 * invokes the ode_system for external events as needed.
+		 */
+		void delta_ext(T e, const Bag<X>& xb)
+		{
+			bool state_event_exists = false;
+			event_happened = true;
+			// Check that we have not missed a state event
+			if (event_exists)
+			{
+				for (int i = 0; i < sys->numVars(); i++)
+					q_trial[i] = q[i];
+				solver->advance(q_trial,e);
+				state_event_exists =
+					event_finder->find_events(event,q,q_trial,solver,e);
+				// We missed an event
+				if (state_event_exists)
+				{
+					output_func(missedOutput);
+					sys->confluent_event(q_trial,event,xb); 
+					for (int i = 0; i < sys->numVars(); i++)
+						q[i] = q_trial[i];
+				}
+			}
+			if (!state_event_exists)// We didn't miss an event
+			{
+				solver->advance(q,e); // Advance the state q by e
+				// Let the model adjust algebraic variables, etc. for the new state
+				sys->postStep(q);
+				// Process the discrete input
+				sys->external_event(q,e+e_accum,xb);
+			}
+			e_accum = 0.0;
+			// Copy the new state to the trial solution 
+			for (int i = 0; i < sys->numVars(); i++) q_trial[i] = q[i];
+			tentative_step(); // Take a tentative step
+		}
+		/**
+		 * Do not override. This method invokes the ode_system method
+		 * for confluent events as needed.
+		 */
+		void delta_conf(const Bag<X>& xb)
+		{
+			if (!missedOutput.empty())
+			{
+				missedOutput.clear();
+				if (sigma > 0.0) event_exists = false;
+			}
+			// Execute any discrete events
+			event_happened = true;
+			if (event_exists) 
+				sys->confluent_event(q_trial,event,xb); 
+			else sys->external_event(q_trial,e_accum+ta(),xb);
+			e_accum = 0.0;
+			// Copy the new state vector to q
+			for (int i = 0; i < sys->numVars(); i++) q[i] = q_trial[i];
+			tentative_step(); // Take a tentative step 
+		}
+		/// Do not override.
+		T ta()
+		{
+			if (missedOutput.empty()) return sigma;
+			else return 0.0;
+		}
+		/// Do not override. Invokes the ode_system output function as needed.
+		void output_func(Bag<X>& yb)
+		{
+			if (!missedOutput.empty())
+			{
+				typename Bag<X>::iterator iter = missedOutput.begin();
+				for (; iter != missedOutput.end(); iter++)
+					yb.insert(*iter);
+				if (sigma == 0.0) // Confluent event
+					sys->output_func(q_trial,event,yb);
+			}
+			else
+			{
+				// Let the model adjust algebraic variables, etc. for the new state
+				sys->postStep(q_trial);
+				if (event_exists)
+					sys->output_func(q_trial,event,yb);
+			}
+		}
+		/// Do not override. Invokes the ode_system gc_output method as needed.
+		void gc_output(Bag<X>& gb) { sys->gc_output(gb); }
+		/// Destructor deletes everything.
+		virtual ~Hybrid()
+		{
+			delete [] q; delete [] q_trial; delete [] event;
+			delete event_finder; delete solver; delete sys;
+		}
+	private:
+		ode_system<X>* sys; // The ODE system
+		ode_solver<X>* solver; // Integrator for the ode set
+		event_locator<X>* event_finder; // Event locator
+		double sigma; // Time to the next internal event
+		double *q, *q_trial; // Current and tentative states
+		bool* event; // Flags indicating the encountered event surfaces
+		bool event_exists; // True if there is at least one event
+		bool event_happened; // True if a discrete event in the ode_system took place
+		double e_accum; // Accumlated time between discrete events
+		Bag<X> missedOutput; // Output missed at an external event
+		// Execute a tentative step and calculate the time advance function
+		void tentative_step()
+		{
+			// Check for a time event
+			double time_event = sys->time_event_func(q);
+			// Integrate up to that time at most
+			double step_size = solver->integrate(q_trial,time_event);
+			// Look for state events inside of the interval [0,step_size]
+			bool state_event_exists =
+				event_finder->find_events(event,q,q_trial,solver,step_size);
+			// Find the time advance and set the time event flag
+			sigma = std::min<double>(step_size,time_event);
+			event[sys->numEvents()] = time_event <= sigma;
+			event_exists = event[sys->numEvents()] || state_event_exists;
+			sys->postTrialStep(q);
+		}
+};
+
+} // end of namespace
+
+#endif

simulators/adevs-3.3/include/adevs_models.h

@@ -0,0 +1,311 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef __adevs_models_h_
+#define __adevs_models_h_
+#include "adevs_time.h"
+#include "adevs_bag.h"
+#include "adevs_set.h"
+#include "adevs_exception.h"
+#include <cstdlib>
+
+namespace adevs
+{
+
+/*
+ * Declare network and atomic model so types can be used as the type of
+ * parent in the basic Devs model and for type ID functions.  
+ */
+template <typename X, typename T> class Network;
+template <typename X, typename T> class Atomic;
+template <typename X, typename T> class MealyAtomic;
+template <typename X, typename T> class Schedule;
+template <typename X, typename T> class Simulator;
+
+/**
+ * The Devs class provides basic operations for all devs models.
+ * The model I/O type is set by the template argument X. The 
+ * type to be used for time is set with the template argument
+ * T. The default type for time is double.
+ */
+template <typename X, typename T = double> class Devs
+{
+	public:
+		/// Default constructor.
+		Devs():
+		parent(NULL)
+		{
+		}
+		/// Destructor.
+		virtual ~Devs()
+		{
+		}
+		/**
+		 * Returns NULL if this is not a network model; returns a pointer to
+		 * itself otherwise. This method is used to avoid a relatively expensive
+		 * dynamic cast.
+		 */
+		virtual Network<X,T>* typeIsNetwork() { return NULL; }
+		/// Returns NULL if this is not an atomic model; returns itself otherwise.
+		virtual Atomic<X,T>* typeIsAtomic() { return NULL; }
+		/// Returns NULL if this is not a mealy atomic model; returns itself otherwise.
+		virtual MealyAtomic<X,T>* typeIsMealyAtomic() { return NULL; }
+		/**
+		 * Get the model that contains this model as a component.  Returns
+		 * NULL if this model is at the top of the hierarchy.
+		 */
+		const Network<X,T>* getParent() const { return parent; }
+		/// Get the model that contains this model as a component.
+		Network<X,T>* getParent() { return parent; }
+		/**
+		 * Assign a new parent to this model. Network model's should always
+		 * call this method to make themselves the parent of their components.
+		 * If the parent is not set correctly, then the event routing algorithm
+		 * in the simulator will fail.
+		 */
+		void setParent(Network<X,T>* parent) { this->parent = parent; }
+		/**
+		 * This is the structure transition function, which is evaluated following
+		 * every change of the model's state. It should return true
+		 * if a structure change is to occur, and false otherwise. False is the
+		 * default return value.
+		 * This method is used by the simulator to limit the execution
+		 * of potentially expensive structure changes. 
+		 * If the return value is true, then the parent's model_transition()
+		 * will also be evaluated. For network models, the model_transition() function is
+		 * preceded and anteceded by a call to getComponents(). The difference
+		 * of these two sets is used to determine if any models were added or removed
+		 * as part of the model transition.
+		 */
+		virtual bool model_transition() { return false; }
+	private:
+		Network<X,T>* parent;
+};
+
+/**
+ * Event objects are used for routing within a network model,
+ * for notifying event listeners of output events, and for injecting
+ * input into a running simulation.
+ */
+template <typename X, typename T = double> class Event
+{
+	public:
+		/// Constructor.  Sets the model to NULL.
+		Event():
+		model(NULL),
+		value()
+		{
+		}
+		/**
+		 * Constructor sets the model and value. The input into a
+		 * Simulator and in a network's routing method,
+		 * the model is the target of the input value.
+		 * In a callback to an event listener, the model is the
+		 * source of the output value. 
+		 */
+		Event(Devs<X,T>* model, const X& value):
+		model(model),
+		value(value)
+		{
+		}
+		/// Copy constructor.
+		Event(const Event<X,T>& src):
+		model(src.model),
+		value(src.value)
+		{
+		}
+		/// Assignment operator.
+		const Event<X,T>& operator=(const Event<X,T>& src)
+		{
+			model = src.model;
+			value = src.value;
+			return *this;
+		}
+		/// The model associated with the event.
+		Devs<X,T>* model;
+		/// The value associated with the event.
+		X value;
+		/// Destructor
+		~Event()
+		{
+		}
+};
+
+/**
+ * Base type for all atomic DEVS models.
+ */
+template <typename X, typename T = double> class Atomic: public Devs<X,T>
+{
+	public:
+		/// The constructor should place the model into its initial state.
+		Atomic():
+			Devs<X,T>(),
+			tL(adevs_zero<T>()),
+			q_index(0), // The Schedule requires this to be zero
+			proc(-1),
+			x(NULL),
+			y(NULL)
+		{
+		}
+		/// Internal transition function.
+		virtual void delta_int() = 0;
+		/**
+		 * External transition function.
+		 * @param e Time elapsed since the last change of state
+		 * @param xb Input for the model.
+		 */  
+		virtual void delta_ext(T e, const Bag<X>& xb) = 0;
+		/**
+		 * Confluent transition function.
+		 * @param xb Input for the model.
+		 */
+		virtual void delta_conf(const Bag<X>& xb) = 0;
+		/**
+		 * Output function.  Output values should be added to the bag yb.
+		 * @param yb Empty bag to be filled with the model's output
+		 */
+		virtual void output_func(Bag<X>& yb) = 0;
+		/**
+		 * Time advance function. adevs_inf<T>() is used for infinity.
+		 * @return The time to the next internal event
+		 */
+		virtual T ta() = 0;
+		/**
+		 * Garbage collection function.  The objects in g are
+		 * no longer in use by the simulation engine and should be disposed of. 
+`		 * Note that the elements in g are only those objects produced as
+		 * output by this model.
+		 */
+		virtual void gc_output(Bag<X>& g) = 0;
+		/// Destructor.
+		virtual ~Atomic(){}
+		/// Returns a pointer to this model.
+		Atomic<X,T>* typeIsAtomic() { return this; }
+	protected:
+		/**
+		 * Get the last event time for this model. This is 
+		 * provided primarily for use with the backwards compatibility
+		 * module and should not be relied on. It is likely to be
+		 * removed in later versions of the code.
+		 */
+		T getLastEventTime() const { return tL; }
+
+	private:
+
+		friend class Simulator<X,T>;
+		friend class Schedule<X,T>;
+
+		// Time of last event
+		T tL;
+		// Index in the priority queue
+		unsigned int q_index;
+		// Thread assigned to this model
+		int proc;
+		// Input and output event bags
+		Bag<X> *x, *y;
+};
+
+/**
+ * This is a Mealy type atomic model where its output
+ * may depend on its input. Mealy machines cannot be
+ * connected to other Mealy machines. An exception
+ * will be generated by the Simulator if you attempt
+ * to do so.
+ */
+
+template <typename X, typename T = double> class MealyAtomic:
+	public Atomic<X,T>
+{
+	public:
+		MealyAtomic<X,T>():Atomic<X,T>(),imm(false){}
+		MealyAtomic<X,T>* typeIsMealyAtomic() { return this; }
+		/**
+		 * Produce output at e < ta(q) in response to xb.
+		 * This is output preceding an external event.
+		 */
+		virtual void output_func(T e, const Bag<X>& xb, Bag<X>& yb) = 0;
+		/**
+		 * Produce output at e = ta(q) in response to xb.
+		 * This is output preceding a confluent event.
+		 */
+		virtual void output_func(const Bag<X>& xb, Bag<X>& yb) = 0;
+		virtual ~MealyAtomic(){}
+
+	private:
+
+		friend class Simulator<X,T>;
+
+		bool imm;
+};
+
+/**
+ * Base class for DEVS network models.
+ */
+template <typename X, typename T = double> class Network: public Devs<X,T>
+{
+	public:
+		/// Constructor.
+		Network():
+		Devs<X,T>()
+		{
+		}
+		/**
+		 * This method should fill the
+		 * set c with all the Network's components, excluding the 
+		 * Network model itself.
+		 * @param c An empty set to the filled with the Network's components.
+		 */
+		virtual void getComponents(Set<Devs<X,T>*>& c) = 0;
+		/**
+		 * This method is called by the Simulator to route an output value
+		 * produced by a model. This method should fill the bag r
+		 * with Events that point to the target model and carry the value
+		 * to be delivered to the target. The target may be a component 
+		 * of the Network or the Network itself, the latter causing the
+		 * Network to produce an output.
+		 * @param model The model that produced the output value
+		 * @param value The output value produced by the model
+		 * @param r A bag to be filled with (target,value) pairs
+		 */
+		virtual void route(const X& value, Devs<X,T>* model, Bag<Event<X,T> >& r) = 0;
+		/**
+		 * Destructor.  This destructor does not delete any component models.
+		 * Any necessary cleanup should be done by the derived class.
+		 */
+		virtual ~Network()
+		{
+		}
+		/// Returns a pointer to this model.
+		Network<X,T>* typeIsNetwork() { return this; }
+};
+
+} // end of namespace
+
+#endif

simulators/adevs-3.3/include/adevs_poly.h

@@ -0,0 +1,95 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef _adevs_poly_h_
+#define _adevs_poly_h_
+#include <cstdlib>
+
+namespace adevs
+{
+    /**
+     * This class implements Lagrange interpolating
+     * polynomials of arbitrary order for functions of
+     * a single variable. This can be 
+     * useful in many discrete event simulations where
+     * tabular values, possible produced by discrete time
+     * calculations, need to be interpolated for use
+     * in a discrete event system. GDEVS is one particular
+     * example of this.
+     */
+    class InterPoly
+    {
+	public:
+	    /**
+	     * Construct a polynomial to interpolate
+	     * a function u(t). The u- values
+	     * are the dependent variable and t- the independent
+	     * variables.
+	     */
+	    InterPoly(const double* u, const double* t, unsigned int n);
+		/**
+		 * Construct a polynomial to interpolate u(t) with data points
+		 * that are regularly spaced in time from an offset t0
+		 */
+		InterPoly(const double* u, double dt, unsigned int n, double t0 = 0.0);
+		/**
+		 * Assign new values to the data set. If t is NULL, then
+		 * only new u values will be assigned and the old t data is
+		 * kept.
+		 */
+		void setData(const double* u, const double* t = NULL);
+	    /**
+	     * Get the interpolated value at t
+	     */
+	    double interpolate(double t) const;
+	    /**
+	     * Overloaded operator for the interpolate method
+	     */
+	    double operator()(double t) const;
+	    /**
+	     * Approximate the function derivative at t
+	     */
+	    double derivative(double t) const;
+	    /**
+	     * Destructor
+	     */
+	    ~InterPoly();
+	private:
+	    InterPoly(){}
+	    InterPoly(const InterPoly&){}
+	    void operator=(const InterPoly&){}
+
+	    double* tdat;
+	    double* udat;
+	    unsigned int n;
+    };
+}
+
+#endif

simulators/adevs-3.3/include/adevs_rand.h

@@ -0,0 +1,233 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef __adevs_rand_h_
+#define __adevs_rand_h_
+#include "adevs.h"
+#include <cstdlib>
+
+namespace adevs
+{
+
+/**
+ * The random_seq class is an abstract interface to a random sequence
+ * generator.
+ */
+class random_seq 
+{
+	public:
+		/// Set the seed for the random number generator
+		virtual void set_seed(unsigned long seed) = 0;
+		/// Get the next double uniformly distributed in [0, 1]
+		virtual double next_dbl() = 0;
+		/// Copy the random number generator
+		virtual random_seq* copy() const = 0;
+		/// Get the next unsigned long
+		virtual unsigned long next_long() = 0;
+		/// Destructor
+		virtual ~random_seq(){}
+};
+
+/**
+ * The crand class provides random number sequences using the standard
+ * C rand_r() function. Each instance of crand generates its own random
+ * number sequence, and the clone method saves the state of the random
+ * number generator. This class can be used in parallel simulations. 
+ */
+class crand: public random_seq 
+{
+	public:
+		/// Create a generator with the default seed
+		crand():seedp(0){}
+		/// Copy constructor
+		crand(const crand& src):seedp(src.seedp){}
+		/// Create a generator with the given seed
+		crand(unsigned long seed):seedp((unsigned int)seed){}
+		/// Set the seed for the random number generator
+		void set_seed(unsigned long seed) { seedp = (unsigned int)seed; }
+		/// Get the next double uniformly distributed in [0, 1]
+		double next_dbl()
+		{
+			#ifdef _WIN32
+			return ((double)next_long()/(double)UINT_MAX);
+			#else
+			return ((double)next_long()/(double)RAND_MAX);
+			#endif
+		}
+		/// Get the next unsigned long
+		unsigned long next_long();
+		/// Copy the random number generator
+		random_seq* copy() const { return new crand(*this); }
+		/// Destructor
+		~crand(){}
+	private:
+		unsigned int seedp;
+};
+
+/**
+ * <p>The rv class provides a random variable based on a selectable
+ * implementation.  By default, this implementation is crand.
+ * I recommend that you find a modern random variate generator,
+ * such as the ones included in the new C++ standards,
+ * rather than use the one here.</p>
+ * <p>The assortment of random variable types was contributed by
+ * Alex Cave (who, at the time, was with the Intelligent 
+ * Systems Automation Group in the School of Engineering at
+ * Deakin University).</p>
+ */
+class rv 
+{
+	public:
+		/// Create a random variable with the default implementation.
+		rv (unsigned long seed = 1);
+		/**
+		 * Create a random variable with the desired implementation.  The
+		 * implementation class is adopted by the rv.
+		 */
+		rv(random_seq* rand);
+		/// Copy constructor relies on copy method of underlying stream.
+		rv(const rv& src);
+		/// Assignment operator relies on copy method of underlying stream.
+		const rv& operator=(const rv& src);
+		/// See the random number generator implementation
+		void set_seed(unsigned long seed);
+		/// Get a raw value from the underlying random number generator
+		unsigned long next_long();
+		/// Sample a triangular distribution with minimum a, maximum b,
+		/// and mode c.
+		double triangular(double a, double b, double c);
+		/// Sample a uniform distribution in the range [a, b]
+		double uniform(double a, double b);
+		/**
+		 * Sample a normally distributed random variable with mean m and 
+		 * standard deviation s.
+		 */
+		double normal(double m, double s);
+		/**
+		 * return a negative exponentially distributed random number with
+		 * the mean as parameter
+		 */
+		double exponential(double a);
+		/**
+		 * return a hyperexponentially distributed random number with
+		 * the means as parameters to two exponentially distributed variates,
+		 * the first with a chance of p of being generated, the second with a
+		 * chance of 1-p.
+		 */
+		double hyperexponential(double p,double a,double b);
+		/// Return a laplace number with the given parameter
+		double laplace(double a);
+		/// Sample a chisquare random variable with n degrees of freedom
+		double chisquare(unsigned int n);
+		/// Sample a student-t random variable
+		double student(unsigned int n);
+		/// Sample a lognormal random variable
+		double lognormal(double a,double b);
+		/// Sample an erlang distribution
+		double erlang(unsigned int n,double a);
+		/// Sample a gamma random variable
+		double gamma(double a,double b);
+		/// Sample a beta random variable
+		double beta(double a,double b);
+		/// Sample a F-distributed random variable
+		double fdistribution(unsigned int n,unsigned int m);
+		/// Sample a Poisson random variable
+		double poisson(double a);
+		/// Sample a geometric random variable with event probability p
+		double geometric(double p);
+		/**
+		 * return a variate from the hypergeometric distribution with m the
+		 * population, p the chance on success and n the number of items drawn
+		 */
+		double hypergeometric(unsigned int m,unsigned int n,double p);
+		/// Sample a weibull random variable
+		double weibull(double a,double b);
+		/**
+		 *  An event count for a binomial distribution with event
+		 *  probability p and n the number of trials
+		 */ 
+		double binomial(double p,unsigned int n);
+		/**
+		 * return a random variable with probabilty of success equal to p
+		 * and n as the number of successes
+		 */
+		double negativebinomial(double p,unsigned int n);
+		/// Sample a triangular random variable
+		double triangular(double a);
+		/// return TRUE with probability p and FALSE with probability 1-p
+		bool probability(double p);
+		/// Sample a lognormal gamma random variable
+		double lngamma(double xx);
+		/// Destructor
+		~rv();
+	private:	
+		random_seq* _impl;
+		typedef enum 
+		{
+			NDURATION,
+			NCOUNT,
+			EMPTY,
+			NDELAY,
+			ADDTOQ,
+			EMPTYQ,
+			HNCLMS,
+			HMODE,
+			PROBVAL,
+			ERRUNIFORM,
+			ERRNORMAL,
+			ERRLOGNORM,
+			ERRTRIANG,
+			ERRGAMMA,
+			ERRBETA,
+			ERREXPONENT,
+			ERRERLANG,
+			ERRHYPGEO,
+			NULLEV,
+			NOHISTO,
+			INITERR,
+			AMODE,
+			HFORM,
+			ERRFILE,
+			SAMPLE,
+			FRACTION,
+			LEVEL,
+			SCAN,
+			SUPPRESS,
+			SEED
+		} errorType;
+
+		void err(errorType n);
+};
+
+} // end of namespace
+
+#endif
+
+

simulators/adevs-3.3/include/adevs_rk_45.h

@@ -0,0 +1,171 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef _adevs_rk_45_h_
+#define _adevs_rk_45_h_
+#include "adevs_hybrid.h"
+#include <cmath>
+
+namespace adevs
+{
+
+/**
+ * This ode_solver implements a 4th/5th order integrator that adjust
+ * its step size to control error.
+ */
+template <typename X> class rk_45:
+	public ode_solver<X>
+{
+	public:
+		/**
+		 * The integrator will adjust its step size to maintain a per
+		 * step error less than err_tol, and will use a step size
+		 * strictly less than h_max.
+		 */
+		rk_45(ode_system<X>* sys, double err_tol, double h_max);
+		/// Destructor
+		~rk_45();
+		double integrate(double* q, double h_lim);
+		void advance(double* q, double h);
+	private:
+		double *dq, // derivative
+			   *qq, // trial solution
+			   *t,  // temporary variables for computing stages
+			   *k[6]; // the six RK stages
+		const double err_tol; // Error tolerance
+		const double h_max; // Maximum time step
+		double h_cur; // Previous successful step size
+		// Compute a trial step of size h, store the result in qq, and return the error
+		double trial_step(double h);
+};
+
+template <typename X>
+rk_45<X>::rk_45(ode_system<X>* sys, double err_tol, double h_max):
+	ode_solver<X>(sys),err_tol(err_tol),h_max(h_max),h_cur(h_max)
+{
+	for (int i = 0; i < 6; i++)
+		k[i] = new double[sys->numVars()];
+	dq = new double[sys->numVars()];
+	qq = new double[sys->numVars()];
+	t = new double[sys->numVars()];
+}
+
+template <typename X>
+rk_45<X>::~rk_45()
+{
+	delete [] dq;
+	delete [] t;
+	for (int i = 0; i < 6; i++)
+		delete [] k[i];
+}
+
+template <typename X>
+void rk_45<X>::advance(double* q, double h)
+{
+	double dt;
+	while ((dt = integrate(q,h)) < h) h -= dt;
+}
+
+template <typename X>
+double rk_45<X>::integrate(double* q, double h_lim)
+{
+	// Initial error estimate and step size
+	double err = DBL_MAX,
+		   h = std::min<double>(h_cur*1.1,std::min<double>(h_max,h_lim));
+	for (;;) {
+		// Copy q to the trial vector
+		for (int i = 0; i < this->sys->numVars(); i++) qq[i] = q[i];
+		// Make the trial step which will be stored in qq
+		err = trial_step(h);
+		// If the error is ok, then we have found the proper step size
+		if (err <= err_tol) {
+			if (h_cur <= h_lim) h_cur = h;
+			break;
+		}
+		// Otherwise shrink the step size and try again
+		else {
+			double h_guess = 0.8*pow(err_tol*pow(h,4.0)/fabs(err),0.25);
+			if (h < h_guess) h *= 0.8;
+			else h = h_guess;
+		}
+	}
+	// Copy the trial solution to q and return the step size that was selected
+	for (int i = 0; i < this->sys->numVars(); i++) q[i] = qq[i];
+	return h;
+}
+
+template <typename X>
+double rk_45<X>::trial_step(double step)
+{
+	// Compute k1
+	this->sys->der_func(qq,dq); 
+	for (int j = 0; j < this->sys->numVars(); j++) k[0][j] = step*dq[j];
+	// Compute k2
+	for (int j = 0; j < this->sys->numVars(); j++) t[j] = qq[j] + 0.5*k[0][j];
+	this->sys->der_func(t,dq);
+	for (int j = 0; j < this->sys->numVars(); j++) k[1][j] = step*dq[j];
+	// Compute k3
+	for (int j = 0; j < this->sys->numVars(); j++) t[j] = qq[j] + 0.25*(k[0][j]+k[1][j]);
+	this->sys->der_func(t,dq);
+	for (int j = 0; j < this->sys->numVars(); j++) k[2][j] = step*dq[j];
+	// Compute k4
+	for (int j = 0; j < this->sys->numVars(); j++) t[j] = qq[j] - k[1][j] + 2.0*k[2][j];
+	this->sys->der_func(t,dq);
+	for (int j = 0; j < this->sys->numVars(); j++) k[3][j] = step*dq[j];
+	// Compute k5
+	for (int j = 0; j < this->sys->numVars(); j++)
+		t[j] = qq[j] + (7.0/27.0)*k[0][j] + (10.0/27.0)*k[1][j] + (1.0/27.0)*k[3][j];
+	this->sys->der_func(t,dq);
+	for (int j = 0; j < this->sys->numVars(); j++) k[4][j] = step*dq[j];
+	// Compute k6
+	for (int j = 0; j < this->sys->numVars(); j++)
+		t[j] = qq[j] + (28.0/625.0)*k[0][j] - 0.2*k[1][j] + (546.0/625.0)*k[2][j]
+			+ (54.0/625.0)*k[3][j] - (378.0/625.0)*k[4][j];
+	this->sys->der_func(t,dq);
+	for (int j = 0 ; j < this->sys->numVars(); j++) k[5][j] = step*dq[j];
+	// Compute next state and the approximate error
+	double err = 0.0;
+	for (int j = 0; j < this->sys->numVars(); j++)
+	{
+		// Next state
+		qq[j] += (1.0/24.0)*k[0][j] + (5.0/48.0)*k[3][j] + 
+			(27.0/56.0)*k[4][j] + (125.0/336.0)*k[5][j];
+		// Componennt wise maximum of the approximate error
+		err = std::max(err,
+				fabs(k[0][j]/8.0+2.0*k[2][j]/3.0+k[3][j]/16.0-27.0*k[4][j]/56.0
+					-125.0*k[5][j]/336.0));
+	}
+	// Return the error
+	return err;
+}
+
+} // end of namespace
+#endif
+

simulators/adevs-3.3/include/adevs_sched.h

@@ -0,0 +1,240 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef __adevs_schedule_h_
+#define __adevs_schedule_h_
+#include "adevs_time.h"
+#include "adevs_models.h"
+#include <cfloat>
+#include <cstdlib>
+#include <list>
+
+using namespace std;
+
+namespace adevs
+{
+
+/**
+ * This is a binary heap for scheduling Atomic models. The Schedule uses 
+ * the q_index attribute of the Atomic model to keep track of the object's
+ * position in the heap. Therefore, a model can be put into only one Schedule.
+ * Please observe that the q_index value for a model must be initialized
+ * to zero before it is placed into the heap for the first time.
+ */
+template <class X, class T = double> class Schedule
+{
+	public:
+		/**
+		 * An interface for objects that want to visit the imminent models
+		 * in the schedule.
+		 */
+		class ImminentVisitor
+		{
+			public:
+				virtual void visit(Atomic<X,T>* model) = 0;
+				virtual ~ImminentVisitor(){}
+		};
+		/// Creates a scheduler with the default or specified initial capacity.
+		Schedule(unsigned int capacity = 100):
+		capacity(capacity),size(0),heap(new heap_element[capacity])
+		{
+			heap[0].priority = adevs_sentinel<T>(); // This is a sentinel value
+		}
+		/// Get the model at the front of the queue.
+		Atomic<X,T>* getMinimum() const { return heap[1].item; }
+		/// Get the time of the next event.
+		T minPriority() const { return heap[1].priority; }
+		/// Visit the imminent models.
+		void visitImminent(ImminentVisitor* visitor) const {
+			visitImminent(visitor,1);
+		}
+		/// Remove the model at the front of the queue.
+		void removeMinimum();
+		/// Add, remove, or move a model as required by its priority.
+		void schedule(Atomic<X,T>* model, T priority);
+		/// Returns true if the queue is empty, and false otherwise.
+		bool empty() const { return size == 0; }
+		/// Get the number of elements in the heap.
+		unsigned int getSize() const { return size; }
+		/// Destructor.
+		~Schedule() { delete [] heap; }
+	private:
+		// Definition of an element in the heap.
+		struct heap_element 
+		{
+			Atomic<X,T>* item;
+			T priority;
+			// Constructor initializes the item and priority
+			heap_element():
+				item(NULL),priority(adevs_inf<T>()){}
+		};
+		unsigned int capacity, size;
+		heap_element* heap;
+		/// Double the schedule capacity
+		void enlarge();
+		/// Move the item at index down and return its new position
+		unsigned int percolate_down(unsigned int index, T priority);
+		/// Move the item at index up and return its new position
+		unsigned int percolate_up(unsigned int index, T priority);
+		/// Visit the imminent set recursively
+		void visitImminent(ImminentVisitor* visitor, unsigned int root) const;
+};
+
+template <class X, class T>
+void Schedule<X,T>::visitImminent(typename Schedule<X,T>::ImminentVisitor* visitor,
+		unsigned int root) const
+{
+	// Stop if the bottom is reached or the next priority is not equal to the minimum
+	if (root > size || heap[1].priority < heap[root].priority)
+		return;
+	visitor->visit(heap[root].item);
+	// Look for more imminent models in the left sub-tree
+	visitImminent(visitor,root*2);
+	// Look in the right sub-tree
+	visitImminent(visitor,root*2+1);
+}
+
+template <class X, class T>
+void Schedule<X,T>::removeMinimum()
+{
+	// Don't do anything if the heap is empty
+	if (size == 0) return;
+	size--;
+	// Set index to 0 to show that this model is not in the schedule
+	heap[1].item->q_index = 0;
+	// If the schedule is empty, set the priority of the last element to adevs_inf
+	if (size == 0)
+	{
+		heap[1].priority = adevs_inf<T>();
+		heap[1].item = NULL;
+	}
+	// Otherwise fill the hole left by the deleted model
+	else
+	{
+		unsigned int i = percolate_down(1,heap[size+1].priority);
+		heap[i] = heap[size+1];
+		heap[i].item->q_index = i;
+		heap[size+1].item = NULL;
+	}
+}
+
+template <class X, class T>
+void Schedule<X,T>::schedule(Atomic<X,T>* model, T priority)
+{
+	// If the model is in the schedule
+	if (model->q_index != 0)
+	{
+		// Remove the model if the next event time is infinite
+		if (!(priority < adevs_inf<T>())) 
+		{
+			// Move the item to the top of the heap
+			T min_priority = minPriority();
+			model->q_index = percolate_up(model->q_index,min_priority);
+			heap[model->q_index].priority = min_priority;
+			heap[model->q_index].item = model;
+			// Remove it and return
+			removeMinimum();
+			return;
+		}
+		// Decrease the time to next event
+		else if (priority < heap[model->q_index].priority)
+			model->q_index = percolate_up(model->q_index,priority);
+		// Increase the time to next event
+		else if (heap[model->q_index].priority < priority)
+			model->q_index = percolate_down(model->q_index,priority);
+		// Don't do anything if the priority is unchanged
+		else return;
+		heap[model->q_index].priority = priority;
+		heap[model->q_index].item = model;
+	}
+	// If it is not in the schedule and the next event time is
+	// not at infinity, then add it to the schedule
+	else if (priority < adevs_inf<T>())
+	{
+		// Enlarge the heap to hold the new model
+		size++;
+		if (size == capacity) enlarge();
+		// Find a slot and put the item into it
+		model->q_index = percolate_up(size,priority);
+		heap[model->q_index].priority = priority;
+		heap[model->q_index].item = model;
+	}
+	// Otherwise, the model is not enqueued and has no next event
+}
+
+template <class X, class T>
+unsigned int Schedule<X,T>::percolate_down(unsigned int index, T priority)
+{
+	unsigned int child;
+	for (; index*2 <= size; index = child) 
+	{
+		child = index*2;
+		if (child != size && heap[child+1].priority < heap[child].priority) 
+		{
+			child++;
+		}
+		if (heap[child].priority < priority) 
+		{
+			heap[index] = heap[child];
+			heap[index].item->q_index = index;
+		}
+		else break;
+	}
+	return index;
+}
+
+template <class X, class T>
+unsigned int Schedule<X,T>::percolate_up(unsigned int index, T priority)
+{
+	// Position 0 has priority -1 and this method is always called
+	// with priority >= 0 and index > 0. 
+	while (priority <= heap[index/2].priority) 
+	{
+		heap[index] = heap[index/2];
+		heap[index].item->q_index = index;
+		index /= 2;
+	}
+	return index;
+}
+
+template <class X, class T>
+void Schedule<X,T>::enlarge()
+{
+	heap_element* rheap = new heap_element[capacity*2];
+	for (unsigned int i = 0; i < capacity; i++)
+		rheap[i] = heap[i];
+	capacity *= 2;
+	delete [] heap;
+	heap = rheap;
+}
+
+} // end of namespace
+
+#endif

simulators/adevs-3.3/include/adevs_set.h

@@ -0,0 +1,59 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef _adevs_set_h
+#define _adevs_set_h
+#include <set>
+#include <algorithm>
+
+namespace adevs
+{
+
+/**
+ * This Set is just an STL set. 
+ */
+template <class T> class Set: public std::set<T>
+{
+};
+
+/// Set difference operator. Returns the set A-B.
+template <class T> 
+void set_assign_diff(Bag<T>& result, const Set<T>& A, const Set<T>& B)
+{
+	typename Set<T>::const_iterator iter = A.begin();
+	for (; iter != A.end(); iter++)
+	{
+		if (B.find(*iter) == B.end()) result.insert(*iter);
+	}
+}
+
+} // end of namespace
+
+#endif

simulators/adevs-3.3/include/adevs_simpledigraph.h

@@ -0,0 +1,133 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef __adevs_simpledigraph_h_
+#define __adevs_simpledigraph_h_
+#include "adevs.h"
+#include <map>
+#include <set>
+#include <cstdlib>
+
+namespace adevs
+{
+
+/**
+ * This is a very simple digraph model for connecting single input/single
+ * output systems. Output generated by a component model is sent to all
+ * components connected to it.
+ */
+template <class VALUE, class T = double> class SimpleDigraph: 
+public Network<VALUE,T>
+{
+	public:
+		/// A component of the SimpleDigraph model
+		typedef Devs<VALUE,T> Component;
+
+		/// Construct a network without components.
+		SimpleDigraph():
+		Network<VALUE,T>()
+		{
+		}
+		/// Add a model to the network.
+		void add(Component* model);
+		/// Couple the source model to the destination model.  
+		void couple(Component* src, Component* dst);
+		/// Puts the network's set of components into c
+		void getComponents(Set<Component*>& c);
+		/// Route an event according to the network's couplings
+		void route(const VALUE& x, Component* model, 
+		Bag<Event<VALUE,T> >& r);
+		/// Destructor.  Destroys all of the component models.
+		~SimpleDigraph();
+
+	private:	
+		// The set of components
+		Set<Component*> models;
+		// Coupling information
+		std::map<Component*,Bag<Component*> > graph;
+};
+
+template <class VALUE, class T>
+void SimpleDigraph<VALUE,T>::add(Component* model)
+{
+	assert(model != this);
+	models.insert(model);
+	model->setParent(this);
+}
+
+template <class VALUE, class T>
+void SimpleDigraph<VALUE,T>::couple(Component* src, Component* dst) 
+{
+	if (src != this) add(src);
+	if (dst != this) add(dst);
+	graph[src].insert(dst);
+}
+
+template <class VALUE, class T>
+void SimpleDigraph<VALUE,T>::getComponents(Set<Component*>& c)
+{
+	c = models;
+}
+
+template <class VALUE, class T>
+void SimpleDigraph<VALUE,T>::
+route(const VALUE& x, Component* model, 
+Bag<Event<VALUE,T> >& r)
+{
+	// Find the list of target models and ports
+	typename std::map<Component*,Bag<Component*> >::iterator graph_iter;
+	graph_iter = graph.find(model);
+	// If no target, just return
+	if (graph_iter == graph.end()) return;
+	// Otherwise, add the targets to the event bag
+	Event<VALUE,T> event;
+	typename Bag<Component*>::iterator node_iter;
+	for (node_iter = (*graph_iter).second.begin();
+	node_iter != (*graph_iter).second.end(); node_iter++)
+	{
+		event.model = *node_iter;
+		event.value = x;
+		r.insert(event);
+	}
+}
+
+template <class VALUE, class T>
+SimpleDigraph<VALUE,T>::~SimpleDigraph()
+{ 
+	typename Set<Component*>::iterator i;
+	for (i = models.begin(); i != models.end(); i++)
+	{
+		delete *i;
+	}
+}
+
+} // end of namespace 
+
+#endif

simulators/adevs-3.3/include/adevs_simulator.h

@@ -0,0 +1,718 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef __adevs_simulator_h_
+#define __adevs_simulator_h_
+#include "adevs_abstract_simulator.h"
+#include "adevs_models.h"
+#include "adevs_event_listener.h"
+#include "adevs_sched.h"
+#include "adevs_bag.h"
+#include "adevs_set.h"
+#include "object_pool.h"
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+
+namespace adevs
+{
+
+/**
+ * This Simulator class implements the DEVS simulation algorithm.
+ * Its methods throw adevs::exception objects if any of the DEVS model
+ * constraints are violated (i.e., a negative time advance, a model
+ * attempting to send an input directly to itself, or coupled Mealy
+ * type systems).
+ */
+template <class X, class T = double> class Simulator:
+	public AbstractSimulator<X,T>,
+	private Schedule<X,T>::ImminentVisitor
+{
+	public:
+		/**
+		 * Create a simulator for a model. The simulator
+		 * constructor will fail and throw an adevs::exception if the
+		 * time advance of any component atomic model is less than zero.
+		 * @param model The model to simulate
+		 */
+		Simulator(Devs<X,T>* model):
+			AbstractSimulator<X,T>(),
+			Schedule<X,T>::ImminentVisitor(),
+			io_up_to_date(false)
+		{
+			schedule(model,adevs_zero<T>());
+		}
+		/**
+		 * Get the model's next event time
+		 * @return The absolute time of the next event
+		 */
+		T nextEventTime()
+		{
+			return sched.minPriority();
+		}
+		/**
+		 * Execute the simulation cycle at time nextEventTime()
+		 * @return The updated simulation time
+		 */
+		T execNextEvent()
+		{
+			return computeNextState();
+		}
+		/**
+		 * Execute until nextEventTime() > tend.
+		 * @return The updated simulation time
+		 */
+		T execUntil(T tend)
+		{
+			T t = tend+adevs_epsilon<T>();
+			while (nextEventTime() <= tend 
+                   && nextEventTime() < adevs_inf<T>()) {
+				t = execNextEvent();
+            }
+			return t;
+		}
+		/**
+		 * Compute the output values of the imminent component models 
+		 * if these values have not already been computed.  This will
+		 * notify registered EventListeners as the outputs are produced. 
+		 */
+		void computeNextOutput();
+		/**
+		 * Compute the output value of the model in response to an input
+		 * at some time in lastEventTime() <= t <= nextEventTime().
+		 * This will notify registered EventListeners as the outputs
+		 * are produced. If this is the first call since the prior
+		 * state change with the given t, then the new output is computed.
+		 * Subsequent calls for the same time t simply
+		 * append to the input already supplied at time t.
+		 * @param input A bag of (input target,value) pairs
+		 * @param t The time at which the input takes effect
+		 */
+		void computeNextOutput(Bag<Event<X,T> >& input, T t);
+		/**
+		 * Apply the bag of inputs at time t and then compute the next model
+		 * states. Requires that lastEventTime() <= t <= nextEventTime().
+		 * This, in effect, implements the state transition function of 
+		 * the resultant model. If the output has already been computed
+		 * at time t, then the new input at t is simply appended to the
+		 * prior input. Otherwise, the old results are discarded and input
+		 * is calculated at the given time.
+		 * @param input A bag of (input target,value) pairs
+		 * @param t The time at which the input takes effect
+		 * @return The new, current simulation time
+		 */
+		T computeNextState(Bag<Event<X,T> >& input, T t);
+		/**
+		 * Compute the next state at the time at the time t and with
+		 * input supplied at the prior call to computeNextOutput
+		 * assuming no computeNextState has intervened. Assumes
+		 * t = nextEventTime() and input an empty bag if there was
+		 * no prior call to computeNextOutput.
+		 * @return The new, current simulation time
+		 */
+		T computeNextState();
+		/**
+		 * Deletes the simulator, but leaves the model intact. The model must
+		 * exist when the simulator is deleted.  Delete the model only after
+		 * the Simulator is deleted.
+		 */
+		~Simulator();
+		/**
+		 * Assign a model to the simulator. This has the same effect as passing
+		 * the model to the constructor.
+		 */
+		void addModel(Atomic<X,T>* model) 
+		{
+			schedule(model,adevs_zero<T>());
+		}
+	private:
+		// Bogus input bag for execNextEvent() method
+		Bag<Event<X,T> > bogus_input;
+		// The event schedule
+		Schedule<X,T> sched;
+		// List of models that are imminent or activated by input
+		Bag<Atomic<X,T>*> activated;
+		// Mealy systems that we need to process
+		bool allow_mealy_input, io_up_to_date;
+		T io_time;
+		Bag<MealyAtomic<X,T>*> mealy;
+		// Pools of preallocated, commonly used objects
+		object_pool<Bag<X> > io_pool;
+		object_pool<Bag<Event<X,T> > > recv_pool;
+		// Sets for computing structure changes.
+		Bag<Devs<X,T>*> added;
+		Bag<Devs<X,T>*> removed;
+		Set<Devs<X,T>*> next;
+		Set<Devs<X,T>*> prev;
+		// Model transition functions are evaluated from the bottom up!
+		struct bottom_to_top_depth_compare
+		{
+			bool operator()(const Network<X,T>* m1, const Network<X,T>* m2) const
+			{
+				unsigned long int d1 = 0, d2 = 0;
+				// Compute depth of m1
+				const Network<X,T>* m = m1->getParent();
+				while (m != NULL) 
+				{
+					d1++;
+					m = m->getParent();
+				}
+				// Compute depth of m2
+				m = m2->getParent();
+				while (m != NULL)
+				{
+					d2++;
+					m = m->getParent();
+				}
+				// Models at the same depth are sorted by name
+				if (d1 == d2) return m1 < m2;
+				// Otherwise, sort by depth
+				return d1 > d2;
+			}
+		};
+		struct top_to_bottom_depth_compare
+		{
+			bool operator()(const Devs<X,T>* m1, const Devs<X,T>* m2) const
+			{
+				unsigned long int d1 = 0, d2 = 0;
+				// Compute depth of m1
+				const Network<X,T>* m = m1->getParent();
+				while (m != NULL) 
+				{
+					d1++;
+					m = m->getParent();
+				}
+				// Compute depth of m2
+				m = m2->getParent();
+				while (m != NULL)
+				{
+					d2++;
+					m = m->getParent();
+				}
+				// Models at the same depth are sorted by name
+				if (d1 == d2) return m1 < m2;
+				// Otherwise, sort by depth
+				return d1 < d2;
+			}
+		};
+		std::set<Network<X,T>*,bottom_to_top_depth_compare> model_func_eval_set;
+		std::set<Devs<X,T>*,top_to_bottom_depth_compare> sorted_removed;
+		/**
+		 * Recursively add the model and its elements to the schedule 
+		 * using t as the time of last event.
+		 */
+		void schedule(Devs<X,T>* model, T t);
+		/// Route an event generated by the source model contained in the parent model.
+		void route(Network<X,T>* parent, Devs<X,T>* src, X& x);
+		/**	
+		 * Add an input to the input bag of an an atomic model. If the 
+		 * model is not already active , then this method adds the model to
+		 * the activated bag.
+		 */
+		void inject_event(Atomic<X,T>* model, X& value);
+		/**
+		 * Recursively remove a model and its components from the schedule 
+		 * and the imminent/activated bags
+		 */
+		void unschedule_model(Devs<X,T>* model);
+		/**
+		 * Delete any thing in the output bag, and return the input
+		 * and output bags to the pools.
+		 * Recursively clean up network model components.
+		 */
+		void clean_up(Devs<X,T>* model);
+		/**
+		 * Construct the complete descendant set of a network model and store it in s.
+		 */
+		void getAllChildren(Network<X,T>* model, Set<Devs<X,T>*>& s);
+		/**
+		 * Visit method inhereted from ImminentVisitor
+		 */
+		void visit(Atomic<X,T>* model);
+};
+
+template <class X, class T>
+void Simulator<X,T>::visit(Atomic<X,T>* model)
+{
+	assert(model->y == NULL);
+	// Mealy models are processed after the Moore models
+	if (model->typeIsMealyAtomic() != NULL)
+	{
+		model->typeIsMealyAtomic()->imm = true;
+		assert(model->y == NULL);
+		// May be in the mealy list because of a route call
+		if (model->x == NULL)
+			mealy.insert(model->typeIsMealyAtomic());
+		return;
+	}
+	model->y = io_pool.make_obj();
+	// Put it in the active list if it is not already there
+	if (model->x == NULL)
+		activated.insert(model);
+	// Compute output functions and route the events. The bags of output
+	// are held for garbage collection at a later time.
+	model->output_func(*(model->y));
+	// Route each event in y
+	for (typename Bag<X>::iterator y_iter = model->y->begin(); 
+		y_iter != model->y->end(); y_iter++)
+	{
+		route(model->getParent(),model,*y_iter);
+	}
+}
+
+template <class X, class T>
+void Simulator<X,T>::computeNextOutput(Bag<Event<X,T> >& input, T t)
+{
+	// Undo any prior output calculation at another time
+	if (io_up_to_date && !(io_time == t))
+	{
+		typename Bag<Atomic<X,T>*>::iterator iter;
+		for (iter = activated.begin(); iter != activated.end(); iter++)
+		{
+			clean_up(*iter);
+		}
+		activated.clear();
+	}
+	// Input and output happen at the current time.
+	io_time = t;
+	// Get the imminent Moore models from the schedule if we have not
+	// already done so.
+	allow_mealy_input = true;
+	if (t == sched.minPriority() && !io_up_to_date)
+		sched.visitImminent(this);
+	// Apply the injected inputs. 
+	for (typename Bag<Event<X,T> >::iterator iter = input.begin(); 
+		iter != input.end(); iter++)
+	{
+		Atomic<X,T>* amodel = (*iter).model->typeIsAtomic();
+		if (amodel != NULL)
+		{
+			inject_event(amodel,(*iter).value);
+		}
+		else
+		{
+			route((*iter).model->typeIsNetwork(),(*iter).model,(*iter).value);
+		}
+	}
+	// Only Moore models can influence Mealy models. 
+	allow_mealy_input = false;
+	// Iterate over activated Mealy models to calculate their output
+	for (typename Bag<MealyAtomic<X,T>*>::iterator m_iter = mealy.begin();
+		m_iter != mealy.end(); m_iter++)
+	{
+		MealyAtomic<X,T> *model = *m_iter;
+		assert(model->y == NULL);
+		model->y = io_pool.make_obj();
+		// Put it in the active set if it is not already there
+		if (model->x == NULL)
+			activated.insert(model);
+		// Compute output functions and route the events. The bags of output
+		// are held for garbage collection at a later time.
+		if (model->imm) // These are the imminent Mealy models
+		{
+			if (model->x == NULL)
+				model->typeIsAtomic()->output_func(*(model->y));
+			else
+				model->output_func(*(model->x),*(model->y));
+		}
+		else
+		{
+			assert(model->x != NULL);
+			// These are the Mealy models activated by input
+			model->output_func(sched.minPriority()-model->tL,*(model->x),*(model->y));
+		}
+	}
+	// Translate Mealy output to inputs for Moore models. The route method
+	// will throw an exception if an event is sent to a Mealy model.
+	for (typename Bag<MealyAtomic<X,T>*>::iterator m_iter = mealy.begin();
+		m_iter != mealy.end(); m_iter++)
+	{
+		MealyAtomic<X,T> *model = *m_iter;
+		// Route each event in y
+		for (typename Bag<X>::iterator y_iter = model->y->begin(); 
+			y_iter != model->y->end(); y_iter++)
+		{
+			route(model->getParent(),model,*y_iter);
+		}
+	}
+	mealy.clear();
+	io_up_to_date = true;
+}
+
+template<class X, class T>
+void Simulator<X,T>::computeNextOutput()
+{
+	computeNextOutput(bogus_input,sched.minPriority());
+}
+
+template <class X, class T>
+T Simulator<X,T>::computeNextState(Bag<Event<X,T> >& input, T t)
+{
+	computeNextOutput(input,t);
+	assert(io_time == t && io_up_to_date);
+	return computeNextState();
+}
+
+template <class X, class T>
+T Simulator<X,T>::computeNextState()
+{
+	if (!io_up_to_date)
+		computeNextOutput();
+	io_up_to_date = false;
+	T t = io_time, tQ = io_time + adevs_epsilon<T>();
+	/*
+	 * Compute the states of atomic models.  Store Network models that 
+	 * need to have their model transition function evaluated in a
+	 * special container that will be used when the structure changes are
+	 * computed.
+	 */
+	for (unsigned i = 0; i < activated.size(); i++)
+	{
+		Atomic<X,T>* model = activated[i];
+		// Internal event if no input
+		if (model->x == NULL)
+			model->delta_int();
+		// Confluent event if model is imminent and has input
+		else if (
+				(model->typeIsMealyAtomic() == NULL && model->y != NULL)
+				|| (model->typeIsMealyAtomic() != NULL && model->typeIsMealyAtomic()->imm)
+			)
+			model->delta_conf(*(model->x));
+		// External event if model is not imminent and has input
+		else
+			model->delta_ext(t-model->tL,*(model->x));
+		// Notify listeners 
+		this->notify_state_listeners(model,tQ);
+		// Check for a model transition
+		if (model->model_transition() && model->getParent() != NULL)
+		{
+			model_func_eval_set.insert(model->getParent());
+		}
+		// Adjust position in the schedule
+		schedule(model,tQ);
+	}
+	/**
+	 * The new states are in effect at t + eps so advance t
+	 */
+	t = tQ;
+	/**
+	 * Compute model transitions and build up the prev (pre-transition)
+	 * and next (post-transition) component sets. These sets are built
+	 * up from only the models that have the model_transition function
+	 * evaluated.
+	 */
+	if (model_func_eval_set.empty() == false)
+	{
+		while (!model_func_eval_set.empty())
+		{
+			Network<X,T>* network_model = *(model_func_eval_set.begin());
+			model_func_eval_set.erase(model_func_eval_set.begin());
+			getAllChildren(network_model,prev);
+			if (network_model->model_transition() &&
+					network_model->getParent() != NULL)
+			{
+				model_func_eval_set.insert(network_model->getParent());
+			}
+			getAllChildren(network_model,next);
+		}
+		// Find the set of models that were added.
+		set_assign_diff(added,next,prev);
+		// Find the set of models that were removed
+		set_assign_diff(removed,prev,next);
+		// Intersection of added and removed is always empty, so no need to look
+		// for models in both (an earlier version of the code did this).
+		next.clear();
+		prev.clear();
+		/** 
+		 * The model adds are processed first.  This is done so that, if any
+		 * of the added models are components something that was removed at
+		 * a higher level, then the models will not have been deleted when
+		 * trying to schedule them.
+		 */
+		for (typename Bag<Devs<X,T>*>::iterator iter = added.begin(); 
+			iter != added.end(); iter++)
+		{
+			schedule(*iter,t);
+		}
+		// Done with the additions
+		added.clear();
+		// Remove the models that are in the removed set.
+		for (typename Bag<Devs<X,T>*>::iterator iter = removed.begin(); 
+			iter != removed.end(); iter++)
+		{
+			clean_up(*iter);
+			unschedule_model(*iter);
+			// Add to a sorted remove set for deletion
+			sorted_removed.insert(*iter); 
+		}
+		// Done with the unsorted remove set
+		removed.clear();
+		// Delete the sorted removed models
+		while (!sorted_removed.empty())
+		{
+			// Get the model to erase
+			Devs<X,T>* model_to_remove = *(sorted_removed.begin());
+			// Remove the model
+			sorted_removed.erase(sorted_removed.begin());
+			/**
+			 * If this model has children, then remove them from the 
+			 * deletion set. This will avoid double delete problems.
+			 */
+			if (model_to_remove->typeIsNetwork() != NULL)
+			{
+				getAllChildren(model_to_remove->typeIsNetwork(),prev);
+				typename Set<Devs<X,T>*>::iterator iter = prev.begin();
+				for (; iter != prev.end(); iter++)
+					sorted_removed.erase(*iter);
+				prev.clear();
+			}
+			// Delete the model and its children
+			delete model_to_remove;
+		}
+		// Removed sets should be empty now
+		assert(prev.empty());
+		assert(sorted_removed.empty());
+	} // End of the structure change
+	// Cleanup and reschedule models that changed state in this iteration
+	// and survived the structure change phase.
+	for (typename Bag<Atomic<X,T>*>::iterator iter = activated.begin(); 
+		iter != activated.end(); iter++)
+	{
+		clean_up(*iter);
+	}
+	// Empty the bags
+	activated.clear();
+	// Return the current simulation time
+	return t;
+}
+
+template <class X, class T>
+void Simulator<X,T>::clean_up(Devs<X,T>* model)
+{
+	Atomic<X,T>* amodel = model->typeIsAtomic();
+	if (amodel != NULL)
+	{
+		if (amodel->x != NULL)
+		{
+			amodel->x->clear();
+			io_pool.destroy_obj(amodel->x);
+			amodel->x = NULL;
+		}
+		if (amodel->y != NULL)
+		{
+			amodel->gc_output(*(amodel->y));
+			amodel->y->clear();
+			io_pool.destroy_obj(amodel->y);
+			amodel->y = NULL;
+		}
+		if (amodel->typeIsMealyAtomic() != NULL)
+		{
+			amodel->typeIsMealyAtomic()->imm = false;
+		}
+	}
+	else
+	{
+		Set<Devs<X,T>*> components;
+		model->typeIsNetwork()->getComponents(components);
+		for (typename Set<Devs<X,T>*>::iterator iter = components.begin();
+		iter != components.end(); iter++)
+		{
+			clean_up(*iter);
+		}
+	}
+}
+
+template <class X, class T>
+void Simulator<X,T>::unschedule_model(Devs<X,T>* model)
+{
+	if (model->typeIsAtomic() != NULL)
+	{
+		sched.schedule(model->typeIsAtomic(),adevs_inf<T>());
+		activated.erase(model->typeIsAtomic());
+	}
+	else
+	{
+		Set<Devs<X,T>*> components;
+		model->typeIsNetwork()->getComponents(components);
+		for (typename Set<Devs<X,T>*>::iterator iter = components.begin();
+		iter != components.end(); iter++)
+		{
+			unschedule_model(*iter);
+		}
+	}
+}
+
+template <class X, class T>
+void Simulator<X,T>::schedule(Devs<X,T>* model, T t)
+{
+	Atomic<X,T>* a = model->typeIsAtomic();
+	if (a != NULL)
+	{
+		a->tL = t;
+		T dt = a->ta();
+		if (dt == adevs_inf<T>())
+		{
+			sched.schedule(a,adevs_inf<T>());
+		}
+		else
+		{
+			T tNext = a->tL+dt;
+			if (tNext < a->tL)
+			{
+				exception err("Negative time advance",a);
+				throw err;
+			}
+			sched.schedule(a,tNext);
+		}
+	}
+	else
+	{
+		Set<Devs<X,T>*> components;
+		model->typeIsNetwork()->getComponents(components);
+		typename Set<Devs<X,T>*>::iterator iter = components.begin();
+		for (; iter != components.end(); iter++)
+		{
+			schedule(*iter,t);
+		}
+	}
+}
+
+template <class X, class T>
+void Simulator<X,T>::inject_event(Atomic<X,T>* model, X& value)
+{
+	if (io_time < model->tL)
+	{
+		exception err("Attempt to apply input in the past",model);
+		throw err;
+	}
+	// If this is a Mealy model, add it to the list of models that
+	// will need their input calculated
+	if (model->typeIsMealyAtomic())
+	{
+		if (allow_mealy_input)
+		{
+			assert(model->y == NULL);
+			// Add it to the list of its not already there
+			if (model->x == NULL && !model->typeIsMealyAtomic()->imm)
+				mealy.insert(model->typeIsMealyAtomic());
+		}
+		else
+		{
+			exception err("Mealy model coupled to a Mealy model",model);
+			throw err;
+		}
+	}
+	// Add the output to the model's bag of output to be processed
+	if (model->x == NULL)
+	{
+		if (model->y == NULL)
+			activated.insert(model);
+		model->x = io_pool.make_obj();
+	}
+	this->notify_input_listeners(model,value,io_time);
+	model->x->insert(value);
+}
+
+template <class X, class T>
+void Simulator<X,T>::route(Network<X,T>* parent, Devs<X,T>* src, X& x)
+{
+	// Notify event listeners if this is an output event
+	if (parent != src)
+		this->notify_output_listeners(src,x,io_time);
+	// No one to do the routing, so return
+	if (parent == NULL) return;
+	// Compute the set of receivers for this value
+	Bag<Event<X,T> >* recvs = recv_pool.make_obj();
+	parent->route(x,src,*recvs);
+	// Deliver the event to each of its targets
+	Atomic<X,T>* amodel = NULL;
+	typename Bag<Event<X,T> >::iterator recv_iter = recvs->begin();
+	for (; recv_iter != recvs->end(); recv_iter++)
+	{
+		/**
+		 * If the destination is an atomic model, add the event to the IO bag
+		 * for that model and add model to the list of activated models
+		 */
+		amodel = (*recv_iter).model->typeIsAtomic();
+		if (amodel != NULL)
+		{
+			inject_event(amodel,(*recv_iter).value);
+		}
+		// if this is an external output from the parent model
+		else if ((*recv_iter).model == parent)
+		{
+			route(parent->getParent(),parent,(*recv_iter).value);
+		}
+		// otherwise it is an input to a coupled model
+		else
+		{
+			this->notify_input_listeners((*recv_iter).model,(*recv_iter).value,io_time);
+			route((*recv_iter).model->typeIsNetwork(),
+			(*recv_iter).model,(*recv_iter).value);
+		}
+	}
+	recvs->clear();
+	recv_pool.destroy_obj(recvs);
+}
+
+template <class X, class T>
+void Simulator<X,T>::getAllChildren(Network<X,T>* model, Set<Devs<X,T>*>& s)
+{
+	Set<Devs<X,T>*> tmp;
+	// Get the component set
+	model->getComponents(tmp);
+	// Add all of the local level elements to s
+	s.insert(tmp.begin(),tmp.end());
+	// Find the components of type network and update s recursively
+	typename Set<Devs<X,T>*>::iterator iter;
+	for (iter = tmp.begin(); iter != tmp.end(); iter++)
+	{
+		if ((*iter)->typeIsNetwork() != NULL)
+		{
+			getAllChildren((*iter)->typeIsNetwork(),s);
+		}
+	}
+}
+
+template <class X, class T>
+Simulator<X,T>::~Simulator()
+{
+	// Clean up the models with stale IO
+	typename Bag<Atomic<X,T>*>::iterator iter;
+	for (iter = activated.begin(); iter != activated.end(); iter++)
+	{
+		clean_up(*iter);
+	}
+}
+
+} // End of namespace
+
+#endif

simulators/adevs-3.3/include/adevs_time.h

@@ -0,0 +1,321 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef __adevs_time_h_
+#define __adevs_time_h_
+#include <cfloat>
+#include <iostream>
+#include <cmath>
+#include <limits>
+
+/// Returns the maximum value for a time type
+template <class T> inline T adevs_inf();
+/// Returns the zero value for a time type
+template <class T> inline T adevs_zero(); 
+/// Returns a value less than adevs_zero()
+template <class T> inline T adevs_sentinel(); 
+/// Returns the interval to the next instant of time
+template <class T> inline T adevs_epsilon(); 
+
+namespace adevs
+{
+
+/**
+ * This time type allows models to evolve on R x Z.
+ */
+template <typename T=double>
+class sd_time
+{
+	public:
+		/// Creates the identify (0,0)
+		sd_time():t(0),k(0){}
+		/// Create a time (t,k)
+		sd_time(T t, int k):t(t),k(k){}
+		/// Copy constructor
+		sd_time(const sd_time& other):t(other.t),k(other.k){}
+		/// Get the real part of time
+		T real() const { return t; }
+		/// Get the logical part of time
+		double integer() const { return k; }
+		/// Assignment operator
+		const sd_time& operator=(const sd_time& other)
+		{
+			t = other.t;
+			k = other.k;
+			return *this;
+		}
+		/// Equivalence
+		bool operator==(const sd_time& t2) const
+		{
+			return (t == t2.t && k == t2.k);
+		}
+		/// Not equal
+		bool operator!=(const sd_time& t2) const
+		{
+			return !(*this == t2);
+		}
+		/// Order by t then by c
+		bool operator<(const sd_time& t2) const
+		{
+			return (t < t2.t || (t == t2.t && k < t2.k));
+		}
+		/// Less than or equal
+		bool operator<=(const sd_time& t2) const
+		{
+			return (*this == t2 || *this < t2);
+		}
+		/// Greater than
+		bool operator>(const sd_time& t2) const
+		{
+			return !(*this <= t2);
+		}
+		/// Greater than or equal
+		bool operator>=(const sd_time& t2) const
+		{
+			return !(*this < t2);
+		}
+		/// Advance this value by a step size t2
+		sd_time operator+(const sd_time& t2) const
+		{
+			sd_time result(*this);
+			result += t2;
+			return result;
+		}
+		/// Advance this value by a step size t2
+		const sd_time& operator+=(const sd_time& t2) 
+		{
+			if (t2.t == 0) k += t2.k;
+			else { t += t2.t; k = t2.k; }
+			return *this;
+		}
+		/// Length of the interval from now to t2
+		sd_time operator-(const sd_time& t2) const
+		{
+			sd_time result(*this);
+			result -= t2;
+			return result;
+		}
+		/// Length of the interval from now to t2
+		const sd_time& operator-=(const sd_time& t2) 
+		{
+			if (t == t2.t) { t = 0; k -= t2.k; }
+			else t -= t2.t; 
+			return *this;
+		}
+		/// Print a time to the output stream
+		friend std::ostream& operator<<(std::ostream& out, const sd_time& t)
+		{
+			out << "(" << t.t << "," << t.k << ")";
+			return out;
+		}
+		/// Read a time from the input stream
+		friend std::istream& operator>>(std::istream& in, sd_time& t)
+		{
+			char junk;
+			in >> junk >> t.t >> junk >> t.k >> junk;
+			return in;
+		}
+	private:
+		T t; int k;
+};
+
+/**
+ * <p>The fcmp() inline function is taken from fcmp.c, which is
+ * Copyright (c) 1998-2000 Theodore C. Belding,
+ * University of Michigan Center for the Study of Complex Systems,
+ * <mailto:Ted.Belding@umich.edu>,
+ * <http://www-personal.umich.edu/~streak/>,
+ * </p>
+ *
+ * <p>This code is part of the fcmp distribution. fcmp is free software;
+ * you can redistribute and modify it under the terms of the GNU Library
+ * General Public License (LGPL), version 2 or later.  This software
+ * comes with absolutely no warranty.</p>
+ */ 
+inline int fcmp(double x1, double x2, double epsilon) 
+{
+    int exponent;
+    double delta;
+    double difference;
+
+    /* Get exponent(max(fabs(x1), fabs(x2))) and store it in exponent. */
+
+    /* If neither x1 nor x2 is 0, */
+    /* this is equivalent to max(exponent(x1), exponent(x2)). */
+
+    /* If either x1 or x2 is 0, its exponent returned by frexp would be 0, */
+    /* which is much larger than the exponents of numbers close to 0 in */
+    /* magnitude. But the exponent of 0 should be less than any number */
+    /* whose magnitude is greater than 0. */
+
+    /* So we only want to set exponent to 0 if both x1 and */
+    /* x2 are 0. Hence, the following works for all x1 and x2. */
+
+    frexp(fabs(x1) > fabs(x2) ? x1 : x2, &exponent);
+
+    /* Do the comparison. */
+
+    /* delta = epsilon * pow(2, exponent) */
+
+    /* Form a neighborhood around x2 of size delta in either direction. */
+    /* If x1 is within this delta neighborhood of x2, x1 == x2. */
+    /* Otherwise x1 > x2 or x1 < x2, depending on which side of */
+    /* the neighborhood x1 is on. */
+
+    delta = ldexp(epsilon, exponent); 
+
+    difference = x1 - x2;
+
+    if (difference > delta)
+        return 1; /* x1 > x2 */
+    else if (difference < -delta)
+        return -1;  /* x1 < x2 */
+    else /* -delta <= difference <= delta */
+        return 0;  /* x1 == x2 */
+}
+
+/**
+ * This is an alternative double that may be used for the simulation clock
+ * (i.e., as the template parameter T for models and simulators). It
+ * uses the fcmp function to check for equality instead of the 
+ * default equality operator. Information on the fcmp function
+ * may be found at http://fcmp.sourceforge.net/
+ */
+class double_fcmp {
+
+private:
+	double d;
+
+public:
+    /**
+	 * The user must instantiate this static variable
+	 * and initialize as required by the fcmp function.
+	 */
+    static double epsilon;
+
+    double_fcmp(double rhs = 0) 
+        : d(rhs) { }
+
+    const double_fcmp& operator=(const double_fcmp& rhs)
+    {
+        d = rhs.d;
+        return *this;
+    }
+    const double_fcmp& operator=(double rhs)
+    {
+        d = rhs;
+        return *this;
+    }
+    operator double()
+    {
+        return d;
+    }
+    bool operator<(double rhs) const
+    {
+        return (fcmp(d, rhs, epsilon) < 0);
+    }
+    bool operator<(const double_fcmp& rhs) const
+    {
+        return (fcmp(d, rhs.d, epsilon) < 0);
+    }
+    bool operator<=(const double_fcmp& rhs) const
+    {
+        return (fcmp(d, rhs.d, epsilon) <= 0);
+    }
+    bool operator>(const double_fcmp& rhs) const
+    {
+        return (fcmp(d, rhs.d, epsilon) > 0);
+    }
+    bool operator>=(const double_fcmp& rhs) const
+    {
+        return (fcmp(d, rhs.d, epsilon) >= 0);
+    }
+    bool operator==(double rhs) const
+    {
+        return (fcmp(d, rhs, epsilon) == 0);
+    }
+    bool operator==(const double_fcmp& rhs) const
+    {
+        return (fcmp(d, rhs.d, epsilon) == 0);
+    }
+};
+
+} // end namespace
+
+template <> inline float adevs_inf() {
+	return std::numeric_limits<float>::max(); }
+template <> inline long double adevs_inf() {
+	return std::numeric_limits<long double>::max(); }
+template <> inline double adevs_inf() {
+	return std::numeric_limits<double>::max(); }
+template <> inline int adevs_inf() {
+	return std::numeric_limits<int>::max(); }
+template <> inline long adevs_inf() {
+	return std::numeric_limits<long>::max(); }
+template <> inline adevs::double_fcmp adevs_inf() {
+	return std::numeric_limits<double>::max(); }
+template <> inline adevs::sd_time<double> adevs_inf() {
+	return adevs::sd_time<double>(std::numeric_limits<double>::max(),std::numeric_limits<int>::max()); }
+template <> inline adevs::sd_time<long> adevs_inf() {
+	return adevs::sd_time<long>(std::numeric_limits<long>::max(),std::numeric_limits<int>::max()); }
+template <> inline adevs::sd_time<int> adevs_inf() {
+	return adevs::sd_time<int>(std::numeric_limits<int>::max(),std::numeric_limits<int>::max()); }
+
+template <> inline float adevs_zero() { return 0.0f; }
+template <> inline long double adevs_zero() { return 0.0L; }
+template <> inline double adevs_zero() { return 0.0; }
+template <> inline int adevs_zero() { return 0; }
+template <> inline long adevs_zero() { return 0; }
+template <> inline adevs::double_fcmp adevs_zero() { return 0.0; }
+template <> inline adevs::sd_time<double> adevs_zero() { return adevs::sd_time<double>(0.0,0); }
+template <> inline adevs::sd_time<long> adevs_zero() { return adevs::sd_time<long>(0,0); }
+template <> inline adevs::sd_time<int> adevs_zero() { return adevs::sd_time<int>(0,0); }
+
+template <> inline float adevs_sentinel() { return -1.0f; }
+template <> inline long double adevs_sentinel() { return -1.0L; }
+template <> inline double adevs_sentinel() { return -1.0; }
+template <> inline int adevs_sentinel() { return -1; }
+template <> inline long adevs_sentinel() { return -1; }
+template <> inline adevs::double_fcmp adevs_sentinel() { return -1.0; }
+template <> inline adevs::sd_time<double> adevs_sentinel() { return adevs::sd_time<double>(-1.0,0); }
+template <> inline adevs::sd_time<long> adevs_sentinel() { return adevs::sd_time<long>(-1,0); }
+template <> inline adevs::sd_time<int> adevs_sentinel() { return adevs::sd_time<int>(-1,0); }
+
+template <> inline float adevs_epsilon() { return 0.0f; }
+template <> inline long double adevs_epsilon() { return 0.0L; }
+template <> inline double adevs_epsilon() { return 0.0; }
+template <> inline int adevs_epsilon() { return 0; }
+template <> inline long adevs_epsilon() { return 0; }
+template <> inline adevs::double_fcmp adevs_epsilon() { return 0.0; }
+template <> inline adevs::sd_time<double> adevs_epsilon() { return adevs::sd_time<double>(0.0,1); }
+template <> inline adevs::sd_time<long> adevs_epsilon() { return adevs::sd_time<long>(0,1); }
+template <> inline adevs::sd_time<int> adevs_epsilon() { return adevs::sd_time<int>(0,1); }
+
+#endif

simulators/adevs-3.3/include/adevs_wrapper.h

@@ -0,0 +1,207 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef __adevs_wrapper_h_
+#define __adevs_wrapper_h_
+#include "adevs_models.h"
+
+namespace adevs
+{
+	/**
+	 * <p>This class wraps a Network or Atomic model with interface type InternalType in an
+	 * Atomic model with interface type ExternalType. Input to the ModelWrapper is passed
+	 * through a user provided input translation method before being handed off
+	 * to the wrapped model for processing. Output from the wrapped model is
+	 * passed through a user provided output translation method before emerging
+	 * as output from the ModelWrapper. If the wrapped model is a Network, the input
+	 * translation method can create inputs for any of its components. Similarly
+	 * the output translation method is provided with every output produced by
+	 * every component in the Network. If the wrapped model is Atomic then there
+	 * is, of course, only one possible destination for incoming events and only
+	 * one source of outgoing events.
+	 * <p>You will need to implement the usual gc_output event for outputs
+	 * produced by the ModelWrapper. You will also need to implement
+	 * gcInputEvents method to clean up events that are created during
+	 * the input translation process.
+	 */
+	template <typename ExternalType, typename InternalType, class T = double> class ModelWrapper:
+		public Atomic<ExternalType,T>,
+		public EventListener<InternalType,T>
+	{
+		public:
+			/**
+			 * Create a wrapper for the specified model. The ModelWrapper takes
+			 * ownership of the supplied model and will delete it when the
+			 * ModelWrapper is deleted.
+			 */
+			ModelWrapper(Devs<InternalType,T>* model);
+			/**
+			 * This method is used to translate incoming input objects into
+			 * input objects that the wrapped model can process. The supplied
+			 * internal_input bag should be filled with Events that contain the targeted
+			 * internal models and the values to supply to them. The external_input
+			 * bag contains the input values supplied to the wrapper's external or
+			 * confluent transition function.
+			 */
+			virtual void translateInput(const Bag<ExternalType>& external_input,
+					Bag<Event<InternalType,T> >& internal_input) = 0;
+			/**
+			 * This method is used to translate outgoing output objects
+			 * into objects that the ModelWrapper can produce. The 
+			 * internal_output bag contains all of the output events that the
+			 * were produced by the wrapped model. The external_output bag
+			 * should be filled with objects of type ExternalType that
+			 * will be produced as output by the ModelWrapper.
+			 */
+			virtual void translateOutput(const Bag<Event<InternalType,T> >& internal_output,
+					Bag<ExternalType>& external_output) = 0;
+			/**
+			 * This is the garbage collection method for internal input events.
+			 * It will be called when the wrapper is done with a set of events
+			 * that you created with the translateInput method. The supplied bag
+			 * is the same one that you filled out in the translateInput method.
+			 */
+			virtual void gc_input(Bag<Event<InternalType,T> >& g) = 0;
+			/// Get the model that is wrapped by this object
+			Devs<InternalType,T>* getWrappedModel() { return model; }
+			/// Atomic internal transition function
+			void delta_int();
+			/// Atomic external transition function
+			void delta_ext(T e, const Bag<ExternalType>& xb);
+			/// Atomic confluent transition function
+			void delta_conf(const Bag<ExternalType>& xb);
+			/// Atomic output function
+			void output_func(Bag<ExternalType>& yb);
+			/// Atomic time advance function
+			T ta();
+			/// EventListener outputEvent method
+			void outputEvent(Event<InternalType,T> y, T t);
+			/// Destructor. This destroys the wrapped model too.
+			~ModelWrapper();
+		private:
+			ModelWrapper(){}
+			ModelWrapper(const ModelWrapper&){}
+			void operator=(const ModelWrapper&){}
+			// Bag of events created by the input translation method 
+			Bag<Event<InternalType,T> > input;
+			// Output from the wrapped model
+			Bag<Event<InternalType,T> > output;
+			// The wrapped model
+			Devs<InternalType,T>* model;
+			// Simulator for driving the wrapped model
+			Simulator<InternalType,T>* sim;
+			// Last event time
+			T tL;
+	};
+
+template <typename ExternalType, typename InternalType, class T> 
+ModelWrapper<ExternalType,InternalType,T>::ModelWrapper(Devs<InternalType,T>* model):
+	Atomic<ExternalType,T>(),
+	EventListener<InternalType,T>(),
+	model(model),
+	tL(adevs_zero<T>())
+{
+	sim = new Simulator<InternalType,T>(model);
+	sim->addEventListener(this);
+}
+
+template <typename ExternalType, typename InternalType, class T> 
+void ModelWrapper<ExternalType,InternalType,T>::delta_int()
+{
+	// Update the internal clock
+	tL = sim->nextEventTime();
+	// Execute the next autonomous event for the wrapped model
+	sim->execNextEvent();
+}
+
+template <typename ExternalType, typename InternalType, class T> 
+void ModelWrapper<ExternalType,InternalType,T>::delta_ext(T e, const Bag<ExternalType>& xb)
+{
+	// Update the internal clock
+	tL += e;
+	// Convert the external inputs to internal inputs
+	translateInput(xb,input);
+	// Apply the input
+	sim->computeNextState(input,tL);
+	// Clean up 
+	gc_input(input);
+	input.clear();
+}
+
+template <typename ExternalType, typename InternalType, class T> 
+void ModelWrapper<ExternalType,InternalType,T>::delta_conf(const Bag<ExternalType>& xb)
+{
+	// Update the internal clock
+	tL = sim->nextEventTime();
+	// Convert the external inputs to internal inputs
+	translateInput(xb,input);
+	// Apply the input
+	sim->computeNextState(input,tL);
+	// Clean up 
+	gc_input(input);
+	input.clear();
+}
+
+template <typename ExternalType, typename InternalType, class T> 
+T ModelWrapper<ExternalType,InternalType,T>::ta()
+{
+	if (sim->nextEventTime() < adevs_inf<T>()) return sim->nextEventTime()-tL;
+	else return adevs_inf<T>();
+}
+
+template <typename ExternalType, typename InternalType, class T> 
+void ModelWrapper<ExternalType,InternalType,T>::output_func(Bag<ExternalType>& yb)
+{
+	// Compute the model's output events; this causes the outputEvent method to be called
+	sim->computeNextOutput();
+	// Translate the output events to external output events
+	translateOutput(output,yb);
+	// Clean up; the contents of the output bag are deleted by the wrapped model's 
+	// gc_output method
+	output.clear();
+}
+
+template <typename ExternalType, typename InternalType, class T> 
+void ModelWrapper<ExternalType,InternalType,T>::outputEvent(Event<InternalType,T> y, T t)
+{
+	// Just save the events for processing by the output_func
+	output.insert(y);
+}
+
+template <typename ExternalType, typename InternalType, class T> 
+ModelWrapper<ExternalType,InternalType,T>::~ModelWrapper()
+{
+	delete sim;
+	delete model;
+}
+
+} // end of namespace
+
+#endif

simulators/adevs-3.3/include/object_pool.h

@@ -0,0 +1,78 @@
+/**
+ * Copyright (c) 2013, James Nutaro
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met: 
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer. 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution. 
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those
+ * of the authors and should not be interpreted as representing official policies, 
+ * either expressed or implied, of the FreeBSD Project.
+ *
+ * Bugs, comments, and questions can be sent to nutaro@gmail.com
+ */
+#ifndef __adevs_object_pool_h_
+#define __adevs_object_pool_h_
+#include "adevs_bag.h"
+
+namespace adevs
+{
+
+/**
+ * A utility class for managing pools of objects that are stored on the heap.
+ * Uses the new and delete operators to create and destroy objects.
+ */
+template <class T> class object_pool
+{
+	public:
+		/// Construct a pool with a specific initial population
+		object_pool(unsigned int pop = 0):
+		pool()
+		{
+			for (unsigned int i = 0; i < pop; i++)
+				pool.insert(new T());
+		}
+		/// Create an object 
+		T* make_obj()
+		{
+			T* obj;
+			if (pool.empty()) obj = new T;
+			else
+			{
+				obj = *((pool.end())--);
+				pool.erase(pool.end()--);
+			}
+			return obj;
+		}
+		/// Return an object to the pool
+		void destroy_obj(T* obj) { pool.insert(obj); }
+		// Delete all objects in the pool
+		~object_pool()
+		{
+			typename Bag<T*>::iterator iter = pool.begin();
+			for (; iter != pool.end(); iter++) delete *iter;
+		}
+	private:
+		Bag<T*> pool;
+};
+
+} // end of namespace
+
+#endif