Merge pull request #41 from siftech/v1.8.0-rc

V1.8.0-rc: WIP: documentation updates for release

.gitignore

@@ -316,3 +316,6 @@ unsat.core
 .gitignore
 core.dimacs
 notebooks/saved-results/out
+**/*.bbl
+**/*.blg
+**/*.out

auxiliary_packages/funman_benchmarks/src/funman_benchmarks/_version.py

@@ -1,4 +1,4 @@
 """Version information."""
 
 # The following line *must* be the last in the module, exactly as formatted:
-__version__ = "1.5.2"
+__version__ = "1.8.0"

auxiliary_packages/funman_demo/src/funman_demo/_version.py

@@ -1,4 +1,4 @@
 """Version information."""
 
 # The following line *must* be the last in the module, exactly as formatted:
-__version__ = "1.7.0"
+__version__ = "1.8.0"

auxiliary_packages/funman_demo/src/funman_demo/example/pde.py

@@ -1,4 +1,5 @@
 import logging
+import os
 from typing import List, Optional
 
 import matplotlib.animation as animation
@@ -15,29 +16,47 @@
 def animate_heat_map(my_df, frames):
     fig = plt.figure()
 
-    data = my_df.loc[0, :]
     vmin = my_df.min().min()
     vmax = my_df.max().max()
     # ax = sns.heatmap(data, vmin=0, vmax=1)
+    # fig, ax = plt.subplots()
+    # sns.heatmap(data, vmin=vmin, vmax=vmax, cmap="crest", ax=ax)
 
-    # def init():
-    #     plt.clf()
-    #     ax = sns.heatmap(data, vmin=0, vmax=1)
+    def init():
+        # plt.clf()
+        fig.clear()
+        data = my_df.loc[0, :]
+        ax = sns.heatmap(data, vmin=vmin, vmax=vmax, cmap="crest")
+        ax.set_xlabel(data.columns[0][0])
+        ax.set_ylabel(data.index.name)
+        ax.set_title(data.columns[0][0])
 
     def animate(i):
-        plt.clf()
+        # plt.clf()
+        fig.clear()
         data = my_df.loc[i, :]
+        # ax.set_data(data)
         ax = sns.heatmap(data, vmin=vmin, vmax=vmax, cmap="crest")
+        ax.set_xlabel(data.columns[0][0])
+        ax.set_ylabel(data.index.name)
+        ax.set_title(f"{data.columns[0][0]}: time = {i}")
 
     anim = animation.FuncAnimation(
-        fig, animate, interval=1000, frames=frames  # init_func=init,
+        fig,
+        animate,
+        interval=1000,
+        frames=frames,
+        init_func=init,
     )
 
     return anim
 
 
 def plot_spatial_timeseries(
-    results: FunmanResults, variables: Optional[List[str]] = None
+    results: FunmanResults,
+    variables: Optional[List[str]] = None,
+    outdir=None,
+    fps=1,
 ):
     logging.getLogger("matplotlib.animation").setLevel(logging.ERROR)
     logging.getLogger("matplotlib.colorbar").setLevel(logging.ERROR)
@@ -62,11 +81,22 @@ def plot_spatial_timeseries(
 
     df.columns = df.columns.str.split("_", expand=True)
     df = df.stack([1])
+    df.index = df.index.set_names(["time"] + [df.columns[0][0]])
 
     anim_h = animate_heat_map(df, steps)
+    if outdir:
+        anim_h.save(
+            os.path.join(outdir, "h.gif"),
+            writer=animation.PillowWriter(fps=fps),
+        )
     hh = HTML(anim_h.to_jshtml())
     dh = df.unstack().diff().fillna(0).stack([1]).rename(columns={"h": "dh"})
     anim_dh = animate_heat_map(dh, steps)
+    if outdir:
+        anim_dh.save(
+            os.path.join(outdir, "dh.gif"),
+            writer=animation.PillowWriter(fps=fps),
+        )
     hdh = HTML(anim_dh.to_jshtml())
 
-    return hh, hdh
+    return hh, hdh, anim_h, anim_dh

auxiliary_packages/funman_dreal/src/funman_dreal/_version.py

@@ -1,4 +1,4 @@
 """Version information."""
 
 # The following line *must* be the last in the module, exactly as formatted:
-__version__ = "1.7.0"
+__version__ = "1.8.0"

auxiliary_packages/funman_dreal/test/test_chime_bilayer_solve.py

@@ -30,7 +30,7 @@
 from pysmt.typing import BOOL, INT, REAL
 
 from funman import Funman
-from funman.model import Model, ModelParameter, QueryLE
+from funman.model import FunmanModel, ModelParameter, QueryLE
 from funman.representation.representation import (
     Point,
     ResultCombinedHandler,

auxiliary_packages/ibex_tests/Makefile

@@ -0,0 +1,11 @@
+.DEFAULT_GOAL := all
+
+num_constraint_test.o: src/num_constraint_test.cc
+	g++ -c -g src/num_constraint_test.cc -I/usr/include/ibex -libex
+
+num_constraint_test: num_constraint_test.o
+	g++ -g num_constraint_test.o -libex -o num_constraint_test
+
+test: num_constraint_test
+
+all: test

auxiliary_packages/ibex_tests/src/num_constraint_test.cc

@@ -0,0 +1,187 @@
+#include "ibex.h"
+#include <iostream>
+#include <chrono>
+
+using namespace ibex;
+using namespace std;
+
+using std::chrono::duration;
+using std::chrono::duration_cast;
+using std::chrono::high_resolution_clock;
+using std::chrono::milliseconds;
+
+NumConstraint func(const Array<const ExprSymbol> &args, const ExprCtr &expr)
+{
+    cout << "size = " << args.size() << endl;
+    switch (args.size())
+    {
+    case 1:
+        return NumConstraint(args[0], expr);
+    case 2:
+        return NumConstraint(args[0], args[1], expr);
+    case 3:
+        return NumConstraint(args[0], args[1], args[2], expr);
+    case 4:
+        return NumConstraint(args[0], args[1], args[2], args[3], expr);
+    case 5:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], expr);
+    case 6:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], expr);
+    case 7:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], expr);
+    case 8:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], expr);
+    case 9:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], expr);
+    case 10:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9], expr);
+    case 11:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9], args[10], expr);
+    case 12:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9], args[10], args[11], expr);
+    case 13:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9], args[10], args[11], args[12], expr);
+    case 14:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9], args[10], args[11], args[12], args[13], expr);
+    case 15:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9], args[10], args[11], args[12], args[13], args[14], expr);
+    case 16:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9], args[10], args[11], args[12], args[13], args[14], args[15], expr);
+    case 17:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9], args[10], args[11], args[12], args[13], args[14], args[15], args[16], expr);
+    case 18:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9], args[10], args[11], args[12], args[13], args[14], args[15], args[16], args[17], expr);
+    case 19:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9], args[10], args[11], args[12], args[13], args[14], args[15], args[16], args[17], args[18], expr);
+    case 20:
+        return NumConstraint(args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9], args[10], args[11], args[12], args[13], args[14], args[15], args[16], args[17], args[18], args[19], expr);
+    default:
+        return NumConstraint(args, expr);
+    }
+}
+
+// void c1()
+// {
+//     int num_vars = 20;
+//     Variable x[num_vars];
+
+//     Array<const ExprSymbol> vars(num_vars);
+
+//     const ExprNode *expr = &(x[0] + 1);
+//     for (int i = 1; i < num_vars; i++)
+//         expr = &(*expr + x[i]);
+//     // expr = &(*expr < 0);
+
+//     for (int i = 0; i < num_vars; i++)
+//     {
+//         vars.set_ref(i, x[i]);
+//         // expr = expr + Function(x[i], x[i]);
+//     }
+//     // Function f(vars, *expr, "f");
+
+//     NumConstraint c1 = func(vars, (*expr < 0)); // x[0] + x[1] + x[2] + x[3] + x[4] = 0);
+//     cout << "Made constraint c1: " << c1 << endl;
+// }
+
+// void c2()
+// {
+//     int num_vars = 20;
+
+//     // Function f(vars, *expr, "f");
+
+//     NumConstraint c2(x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10], x[11], x[12], x[13], x[14], x[15], x[16], x[17], x[18], x[19], *expr < 0); // x[0] + x[1] + x[2] + x[3] + x[4] = 0);
+//     cout << "Made constraint c2: " << c2 << endl;
+// }
+
+duration<double, std::milli> make_constraint(const ExprNode *expr, Array<const ExprSymbol> &vars, bool use_opt)
+{
+    auto t1 = high_resolution_clock::now();
+    // cout << vars.size() << endl;
+    cout << "done expr" << endl;
+
+    NumConstraint c1 = use_opt ? func(vars, (*expr < 0)) : NumConstraint(vars, (*expr < 0));
+
+    cout << "done expr" << endl;
+    auto t2 = high_resolution_clock::now();
+    duration<double, std::milli> ms_double = t2 - t1;
+    return ms_double;
+}
+
+const ExprNode *make_expr(Array<const ExprSymbol> &x)
+{
+
+    int num_vars = x.size();
+    const ExprNode *expr = &(x[0] + 1);
+    for (int i = 1; i < num_vars; i++)
+        expr = &(*expr + x[i]);
+    return expr;
+}
+
+void make_vars(Array<const ExprSymbol> *vars)
+{
+    int num_vars = vars->size();
+    Variable x[num_vars];
+    cout << "num_vars: " << num_vars << endl;
+
+    for (int i = 0; i < num_vars; i++)
+    {
+        cout << "i = " << i << endl;
+        vars->set_ref(i, x[i]);
+    }
+}
+
+void summary(int num_vars, double time_with, double time_without)
+{
+    double speedup = time_without / time_with;
+    cout << num_vars << ": " << time_with << " " << time_without << " (" << speedup << "x)" << endl;
+}
+
+int main()
+{
+    // Variable x;
+    // NumConstraint c(x, x + 1 <= 0);
+    // cout << "Made constraint: " << c << endl;
+    cout << "Starting: " << endl;
+    int num_vars = 20;
+    Variable x[num_vars];
+    Array<const ExprSymbol> vars(num_vars);
+
+    // cout << "Made vars: " << vars.size() << endl;
+    make_vars(&vars);
+    cout << vars[0] << endl;
+    cout << "Made vars: " << vars.size() << endl;
+    const ExprNode *expr = make_expr(vars);
+    cout << "Made expr" << endl;
+    cout << (*expr) << endl;
+    auto time_with = make_constraint(expr, vars, true);
+    auto time_without = make_constraint(expr, vars, false);
+    // summary(vars.size(), time_with.count(), time_without.count());
+    delete expr;
+
+    // auto t1 = high_resolution_clock::now();
+    // c1();
+    // auto t2 = high_resolution_clock::now();
+    // // c2();
+
+    // auto t3 = high_resolution_clock::now();
+
+    // /* Getting number of milliseconds as an integer. */
+    // auto ms_int = duration_cast<milliseconds>(t2 - t1);
+
+    // /* Getting number of milliseconds as a double. */
+    // duration<double, std::milli> ms_double = t2 - t1;
+
+    // std::cout << ms_int.count() << "ms\n";
+    // std::cout << ms_double.count() << "ms\n";
+
+    // /* Getting number of milliseconds as an integer. */
+    // auto ms_int1 = duration_cast<milliseconds>(t3 - t2);
+
+    // /* Getting number of milliseconds as a double. */
+    // duration<double, std::milli> ms_double1 = t3 - t2;
+
+    // // std::cout << ms_int1.count() << "ms\n";
+    // // std::cout << ms_double1.count() << "ms\n";
+
+    return 0;
+}

auxiliary_packages/pde2petri/README.md

@@ -0,0 +1,31 @@
+# PDE to Petri converter
+A package for converting PDEs into Petrinet AMRs
+
+This package contains manual and automated approaches to convert vector space PDEs into Petrinets with custom rate laws.
+
+1) Manual: The script [scripts/generate-pde-amr.py](scripts/generate-pde-amr.py) will write several directories to the current working directory.  Each directory contains alternative formulations of a PDE problem as a Petrinet, including:
+    - `halfar`: ice dome model 
+    - `advection`: advection model for incompressible flows 
+
+    The files in each subdirectory use the naming scheme, as follows:
+
+    `{problem}_{derivative}_{dimensions}_{boundary_slope}_{num_disc}.json`
+
+    where 
+
+    - `problem`: the name of the problem (e.g., `advection`)
+    - `derivative`: the method used to compute spatial derivatives (e.g., `forward`, `backward`, and `centered`)
+    - `dimensions`: the number of spatial dimensions (1, 2, or 3)
+    - `boundary_slope`: the coefficient for boundary conditions, expressed as `u(x, t) = kt`, where `k = boundary_slope`, `t` is the time (relative to starting time at 0), `u` is a state variable, and `x` is a boundary position.
+    - `num_disc`: the number of discrete points in each dimension (e.g., if `dimension = 2` and `num_disc = 5`, then there will be `5^2 = 25` positions, not including boundaries).
+
+    A notebook illustrating the results of FUNMAN analyzing the models is available [here](https://github.com/siftech/funman/blob/pde-amr-examples/notebooks/pde_as_petrinet.ipynb).
+
+2) Automated: The [./test](./test) directory includes tests that illustrate automatically converting a vector space PDE to a Petrinet.  
+
+    The code in [./src](./src) corresponds to the `pde2petri` package, which can be installed with `pip` (e.g., `pip install .` where the current directory is the same directory containing this README).  
+
+    The automated approach involves parsing an expression with the sympy package, discretizing the PDE, and isolating the "next state" variables to determine the Petrinet states and rate laws for the transitions.  The [document](./doc/discretization/main.pdf) describes the general approach taken in more detail, and with examples.
+
+---
+Authored by [email protected] and [email protected]