DDA walkthrought example: Chua attractor¶
What follows is an example of the Chua attractor. It is described in the Analog Paradigm Application Note 3 as well as in section 6.15 in Bernd’s new book (Analog Programming II). The attractor is described by a coupled set of three ordinary differential equations,
with \(f(x)\) a simple function decribing the Chua diode (given algebraically) and a number of parameters \(c_{1,2,3}\). What follows is the scaling of these equations. The resulting set of equations is slightly more verbose. It’s implementation is given in the following traditional DDA file:
[1]:
!cat chua.dda
#
# Copyright (c) 2020 anabrid GmbH
# Contact: https://www.anabrid.com/licensing/
#
# This file is part of the examples of the PyAnalog toolkit.
#
# ANABRID_BEGIN_LICENSE:GPL
# Commercial License Usage
# Licensees holding valid commercial anabrid licenses may use this file in
# accordance with the commercial license agreement provided with the
# Software or, alternatively, in accordance with the terms contained in
# a written agreement between you and Anabrid GmbH. For licensing terms
# and conditions see https://www.anabrid.com/licensing. For further
# information use the contact form at https://www.anabrid.com/contact.
#
# GNU General Public License Usage
# Alternatively, this file may be used under the terms of the GNU
# General Public License version 3 as published by the Free Software
# Foundation and appearing in the file LICENSE.GPL3 included in the
# packaging of this file. Please review the following information to
# ensure the GNU General Public License version 3 requirements
# will be met: https://www.gnu.org/licenses/gpl-3.0.html.
# ANABRID_END_LICENSE
#
# Chua attractor, chapter 6.15 from Bernds book ap2.pdf
# Below is the scaled version (equations 6.40-6.51)
x0 = const(0.1)
x1 = mult(-10, neg(sum(x, fx)))
x2 = neg(sum(y, mult(0.5, x1)))
x = neg(sum(mult(3.12, neg(int(x2, dt, 0))), x0))
y1 = neg(sum(z, neg(mult(0.125, y))))
y2 = neg(sum(mult(1.25, x), mult(2, y1)))
y = mult(4, neg(int(y2, dt, 0)))
z = int(mult(3.5, y), dt, 0)
f1 = abs(sum(mult(0.7143,x), 0.2857))
f2 = abs(sum(mult(0.7143,x), -0.2857))
f3 = neg(sum(f1, neg(f2)))
fx = sum(mult(0.714, x), mult(0.3003, f3))
dt = const(0.001)
In the following, we use the PyDDA library to read in this DDA file and demonstrate the internal representation.
[2]:
from dda.dsl import read_traditional_dda
chua_text = open("chua.dda").read()
state = read_traditional_dda(chua_text)
state
[2]:
State({'dt': const(0.001),
'f1': abs(sum(mult(0.7143, x), 0.2857)),
'f2': abs(sum(mult(0.7143, x), -0.2857)),
'f3': neg(sum(f1, neg(f2))),
'fx': sum(mult(0.714, x), mult(0.3003, f3)),
'x': neg(sum(mult(3.12, neg(int(x2, dt, 0))), x0)),
'x0': const(0.1),
'x1': mult(-10, neg(sum(x, fx))),
'x2': neg(sum(y, mult(0.5, x1))),
'y': mult(4, neg(int(y2, dt, 0))),
'y1': neg(sum(z, neg(mult(0.125, y)))),
'y2': neg(sum(mult(1.25, x), mult(2, y1))),
'z': int(mult(3.5, y), dt, 0)})
Syntax trees: Down into the rabbit hole¶
Obviously, the output of the internal data structure state and the DDA file itself does not differ so much. That is by intention, both look quite pythonic. The state itself is basically a dictionary (mapping) from strings (the left hand sides in the DDA file) to the expressions (their right hand sides). Let’s inspect such an expression.
[3]:
state["z"]
[3]:
int(mult(3.5, y), dt, 0)
[4]:
type(state["z"])
[4]:
dda.ast.Symbol
[5]:
print(state["z"].head)
print(state["z"].tail)
int
(mult(3.5, y), dt, 0)
What we are actually looking at is the PyDDA-representation of a mathematical expression tree. We can visualize this tree:
[6]:
state["z"].draw_graph()
[6]:
Given this, we can compute the dependencies of all variables in the sytem:
[7]:
graph = state.draw_dependency_graph(export_dot=False)
graph
[7]:
<networkx.classes.digraph.DiGraph at 0x7f431c5df3d0>
[9]:
# Draw the graph with matlotlib
import networkx as nx
from networkx.drawing.nx_agraph import graphviz_layout
pos = graphviz_layout(graph)
nx.draw(graph, pos=pos)
[10]:
# draw the graph with graphviz (requires pydot,graphviz)
from networkx.drawing.nx_pydot import to_pydot
from graphviz import Source
nx2dot = lambda graph: Source(to_pydot(graph).to_string())
[11]:
nx2dot(graph)
[11]:
Based on this dependency analysis, one can linearize the state, that is, define an ordering how to compute the state numerically:
[12]:
vars = state.variable_ordering()
[13]:
print("The evolved variables are:", vars.evolved)
print("Auxilliary variables are:", vars.aux.all)
The evolved variables are: ['int_1', 'int_2', 'z']
Auxilliary variables are: ['f1', 'f2', 'f3', 'fx', 'mult_10', 'mult_6', 'mult_9', 'sum_1', 'sum_2', 'sum_3', 'sum_4', 'sum_5', 'sum_6', 'sum_7', 'sum_8', 'x', 'x1', 'x2', 'y', 'y1', 'y2']
One sees a number of new variables. They were introduced by naming all intermediate expressions. What are these intermediates? Let’s review again the variable z:
[14]:
state["z"].draw_graph()
[14]:
Here, the left most child is an intermediate expression, since it computes mult(3.5, y). We can give this intermediate result a concrete name and replace the whole subtree with this name:
[15]:
state.name_computing_elements()["z"].draw_graph()
[15]:
Simulating a circuit¶
In the following, we use the C++ code generator to simulate this simple ordinary differential equation:
[29]:
cpp_code = state.export(to="C")
print(cpp_code[:1000])
print("// ... (in total ", cpp_code.count("\n"), " lines of C/C++ code) ...")
// This code was generated by PyDDA.
#include <cmath> /* don't forget -lm for linking */
#include <cfenv> /* for feraisexcept and friends */
#include <limits> /* for signaling NAN */
#include <vector>
#include <string>
#include <sstream>
#include <algorithm>
#include <map>
#include <cstdio>
#include <iostream>
#include <fstream>
bool debug;
constexpr double _nan_ = std::numeric_limits<double>::signaling_NaN();
namespace dda {
/* if you use an old C++ compiler, just remove the newer features */
#define A constexpr double /* constexpr requires C++11 */
#define D template<typename... T> A /* Variadic templates require C++17 */
// Known limitations for div(int, double): If certain arguments appear as integer
// in the code (i.e. 1 instead of 1.0), there is div(int,int) kicking in from
// cstdlib. TODO: Should rename div to Div; following int->Int.
A neg(double a) { return -a; }
A div(double a, double b) { return a/b; }
D Int(T... a) { return -(a + ...); } // int() is res
// ... (in total 481 lines of C/C++ code) ...
We printed the generated C++ code, which is in fact just a string in python, in the cell above. Next come some shortcut functions which call the system C++ compiler and run the binary, all externally on the system shell. The return value is slurped in as CSV data with numpy, so we readily have it for plotting.
[30]:
from dda.cpp_exporter import compile, run
compile(cpp_code)
[31]:
# This shows the stdout of the binary generated by the above C++ code:
print(run(arguments={"max_iterations":10}, fields_to_export=list("xyz"), return_ndarray=False))
Running: ./a.out --max_iterations=10 x y z
x y z
0.100223 0.0005 0
0.100448 0.00100062 -1.75e-06
0.100675 0.00150184 -5.25215e-06
0.100905 0.00200368 -1.05086e-05
0.101136 0.00250611 -1.75215e-05
0.10137 0.00300915 -2.62929e-05
0.101605 0.00351277 -3.68249e-05
0.101843 0.00401699 -4.91196e-05
0.102082 0.0045218 -6.3179e-05
0.102324 0.00502718 -7.90053e-05
[32]:
# We can slurp in the CSV data directly to a numpy recarray:
result = run(arguments={"max_iterations":1000, "modulo_write":10})
Running: ./a.out --max_iterations=1000 --modulo_write=10
[33]:
result["x"]
[33]:
array([0.100223, 0.102568, 0.105126, 0.107902, 0.110904, 0.114137,
0.11761 , 0.121327, 0.125297, 0.129525, 0.134019, 0.138785,
0.143831, 0.149164, 0.154792, 0.16072 , 0.166958, 0.173512,
0.180391, 0.187603, 0.195156, 0.203058, 0.211319, 0.219946,
0.22895 , 0.23834 , 0.248125, 0.258317, 0.268924, 0.279959,
0.291432, 0.303354, 0.315739, 0.328597, 0.341943, 0.355789,
0.370149, 0.385039, 0.400471, 0.415958, 0.430975, 0.445545,
0.459687, 0.473419, 0.486757, 0.499714, 0.512304, 0.524535,
0.536418, 0.54796 , 0.559167, 0.570045, 0.580597, 0.590828,
0.600738, 0.61033 , 0.619603, 0.628559, 0.637196, 0.645513,
0.653509, 0.661182, 0.668528, 0.675547, 0.682235, 0.688589,
0.694606, 0.700283, 0.705617, 0.710605, 0.715244, 0.719531,
0.723464, 0.727041, 0.730258, 0.733115, 0.73561 , 0.737741,
0.739508, 0.74091 , 0.741947, 0.742621, 0.742931, 0.742879,
0.742467, 0.741697, 0.740573, 0.739097, 0.737275, 0.735109,
0.732606, 0.729772, 0.726611, 0.723132, 0.719342, 0.715248,
0.710859, 0.706184, 0.701233, 0.696016])
The above cell output shows a NumPy array, which is a Python-internal representation of the CSV file printed in the cell above. NumPy arrays are suitable for plotting, as we do next.
What is actually called by cpp_exporter.run() is np.genfromtxt(). As our CSV file has a header row with column names, it creates a numpy recarray. Therefore, we can address the column x by writing result["x"].
[34]:
import matplotlib.pyplot as plt
[35]:
plt.plot(result["x"], label="x")
plt.plot(result["y"], label="y")
plt.xlabel("time (iterations)")
plt.ylabel("field values")
plt.legend()
[35]:
<matplotlib.legend.Legend at 0x7f42f5dd2820>
The above plot shows the time evolution of the quantity x and y. The plot is not very meaningful, but at least we see that the values are well within the analog computer bounds \([-1,1]\).
Let’s run the simulation a bit longer and display a phase space plot of x and y:
[78]:
result = run(arguments={"max_iterations":30000, "modulo_write":50})
plt.xlabel("x"); plt.ylabel("y")
plt.plot(result["x"], result["y"], "o-")
Running: ./a.out --max_iterations=30000 --modulo_write=50
[78]:
[<matplotlib.lines.Line2D at 0x7f12dd23d580>]
That’s it for the moment. If you want to see even more advanced plotting, inspect the run-chua.py file in the directory of this notebook file (i.a. in the experiments/ directory).