% Jupyter Notebook code cell colors
\definecolor{nbsphinxin}{HTML}{307FC1}
\definecolor{nbsphinxout}{HTML}{BF5B3D}
\definecolor{nbsphinx-code-bg}{HTML}{F5F5F5}
\definecolor{nbsphinx-code-border}{HTML}{E0E0E0}
\definecolor{nbsphinx-stderr}{HTML}{FFDDDD}
% ANSI colors for output streams and traceback highlighting
\definecolor{ansi-black}{HTML}{3E424D}
\definecolor{ansi-black-intense}{HTML}{282C36}
\definecolor{ansi-red}{HTML}{E75C58}
\definecolor{ansi-red-intense}{HTML}{B22B31}
\definecolor{ansi-green}{HTML}{00A250}
\definecolor{ansi-green-intense}{HTML}{007427}
\definecolor{ansi-yellow}{HTML}{DDB62B}
\definecolor{ansi-yellow-intense}{HTML}{B27D12}
\definecolor{ansi-blue}{HTML}{208FFB}
\definecolor{ansi-blue-intense}{HTML}{0065CA}
\definecolor{ansi-magenta}{HTML}{D160C4}
\definecolor{ansi-magenta-intense}{HTML}{A03196}
\definecolor{ansi-cyan}{HTML}{60C6C8}
\definecolor{ansi-cyan-intense}{HTML}{258F8F}
\definecolor{ansi-white}{HTML}{C5C1B4}
\definecolor{ansi-white-intense}{HTML}{A1A6B2}
\definecolor{ansi-default-inverse-fg}{HTML}{FFFFFF}
\definecolor{ansi-default-inverse-bg}{HTML}{000000}

% Define an environment for non-plain-text code cell outputs (e.g. images)
\makeatletter
\newenvironment{nbsphinxfancyoutput}{%
    % Avoid fatal error with framed.sty if graphics too long to fit on one page
    \let\sphinxincludegraphics\nbsphinxincludegraphics
    \nbsphinx@image@maxheight\textheight
    \advance\nbsphinx@image@maxheight -2\fboxsep   % default \fboxsep 3pt
    \advance\nbsphinx@image@maxheight -2\fboxrule  % default \fboxrule 0.4pt
    \advance\nbsphinx@image@maxheight -\baselineskip
\def\nbsphinxfcolorbox{\spx@fcolorbox{nbsphinx-code-border}{white}}%
\def\FrameCommand{\nbsphinxfcolorbox\nbsphinxfancyaddprompt\@empty}%
\def\FirstFrameCommand{\nbsphinxfcolorbox\nbsphinxfancyaddprompt\sphinxVerbatim@Continues}%
\def\MidFrameCommand{\nbsphinxfcolorbox\sphinxVerbatim@Continued\sphinxVerbatim@Continues}%
\def\LastFrameCommand{\nbsphinxfcolorbox\sphinxVerbatim@Continued\@empty}%
\MakeFramed{\advance\hsize-\width\@totalleftmargin\z@\linewidth\hsize\@setminipage}%
\lineskip=1ex\lineskiplimit=1ex\raggedright%
}{\par\unskip\@minipagefalse\endMakeFramed}
\makeatother
\newbox\nbsphinxpromptbox
\def\nbsphinxfancyaddprompt{\ifvoid\nbsphinxpromptbox\else
    \kern\fboxrule\kern\fboxsep
    \copy\nbsphinxpromptbox
    \kern-\ht\nbsphinxpromptbox\kern-\dp\nbsphinxpromptbox
    \kern-\fboxsep\kern-\fboxrule\nointerlineskip
    \fi}
\newlength\nbsphinxcodecellspacing
\setlength{\nbsphinxcodecellspacing}{0pt}

% Define support macros for attaching opening and closing lines to notebooks
\newsavebox\nbsphinxbox
\newcommand{\nbsphinxstartnotebook}[1]{%
    \par
    % measure needed space
    \setbox\nbsphinxbox\vtop{{#1\par}}
    % reserve some space at bottom of page, else start new page
    \needspace{\dimexpr2.5\baselineskip+\ht\nbsphinxbox+\dp\nbsphinxbox}
    % mimick vertical spacing from \section command
      \addpenalty\@secpenalty
      \@tempskipa 3.5ex \@plus 1ex \@minus .2ex\relax
      \addvspace\@tempskipa
      {\Large\@tempskipa\baselineskip
             \advance\@tempskipa-\prevdepth
             \advance\@tempskipa-\ht\nbsphinxbox
             \ifdim\@tempskipa>\z@
               \vskip \@tempskipa
             \fi}
    \unvbox\nbsphinxbox
    % if notebook starts with a \section, prevent it from adding extra space
    \@nobreaktrue\everypar{\@nobreakfalse\everypar{}}%
    % compensate the parskip which will get inserted by next paragraph
    \nobreak\vskip-\parskip
    % do not break here
    \nobreak
}% end of \nbsphinxstartnotebook

\newcommand{\nbsphinxstopnotebook}[1]{%
    \par
    \setbox\nbsphinxbox\vbox{{#1\par}}
    \nobreak % it updates page totals
    \dimen@\pagegoal
    \advance\dimen@-\pagetotal \advance\dimen@-\pagedepth
    \advance\dimen@-\ht\nbsphinxbox \advance\dimen@-\dp\nbsphinxbox
    \ifdim\dimen@<\z@
      % little space left
      \unvbox\nbsphinxbox
      \kern-.8\baselineskip
      \nobreak\vskip\z@\@plus1fil
      \penalty100
      \vskip\z@\@plus-1fil
      \kern.8\baselineskip
    \else
    \fi
}% end of \nbsphinxstopnotebook

% Ensure height of an included graphics fits in nbsphinxfancyoutput frame
\newdimen\nbsphinx@image@maxheight % set in nbsphinxfancyoutput environment
\newcommand*{\nbsphinxincludegraphics}[2][]{%
    \gdef\spx@includegraphics@options{#1}%
    \setbox\spx@image@box\hbox{\includegraphics[#1,draft]{#2}}%
    \in@false
    \ifdim \wd\spx@image@box>\linewidth
      \g@addto@macro\spx@includegraphics@options{,width=\linewidth}%
      \in@true
    \fi
    % no rotation, no need to worry about depth
    \ifdim \ht\spx@image@box>\nbsphinx@image@maxheight
      \g@addto@macro\spx@includegraphics@options{,height=\nbsphinx@image@maxheight}%
      \in@true
    \fi
    \ifin@
      \g@addto@macro\spx@includegraphics@options{,keepaspectratio}%
    \fi
    \setbox\spx@image@box\box\voidb@x % clear memory
    \expandafter\includegraphics\expandafter[\spx@includegraphics@options]{#2}%
}% end of "\MakeFrame"-safe variant of \sphinxincludegraphics

\renewcommand*\sphinx@verbatim@nolig@list{\do\'\do\`}
\begingroup
\catcode`'=\active
\let\nbsphinx@noligs\@noligs
\g@addto@macro\nbsphinx@noligs{\let'\PYGZsq}
\endgroup
\renewcommand*\sphinxbreaksbeforeactivelist{\do\<\do\"\do\'}
\renewcommand*\sphinxbreaksafteractivelist{\do\.\do\,\do\:\do\;\do\?\do\!\do\/\do\>\do\-}
\fvset{codes*=\sphinxbreaksattexescapedchars\do\^\^\let\@noligs\nbsphinx@noligs}

\setcounter{tocdepth}{2}

% Fix for allowing to show unicode delta character 
% See also option 'textgreek' in latex_elements
% (https://www.sphinx-doc.org/en/master/latex.html#latex-elements-confval)
\DeclareUnicodeCharacter{0394}{\ensuremath{\Delta}}

% another fix for minimal length bug
\setlength{\headheight}{14pt} % https://tex.stackexchange.com/a/198694

% another fix for very long equations
\allowdisplaybreaks % https://tex.stackexchange.com/a/51883

An example on how to use PyDDA from command line: Gradiend descent

Different ways to access the well-known DDA computing elements in Python

I want to simulate electronics or solve an ODE. Do I really need DDA?

#!/usr/bin/env python3
#
# Copyright (c) 2020 anabrid GmbH
# Contact: https://www.anabrid.com/licensing/
#
# This file is part of the DDA module 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.
# For Germany, additional rules exist. Please consult /LICENSE.DE
# for further agreements.
# ANABRID_END_LICENSE
#

"""
The minimalistic pythonic standalone abstract syntax tree (**AST**)
representation in this module is the heart of the PyDDA package.
The code has no external dependencies, especially it does not rely
on a Computer Algebra System or even on SymPy.

The :class:`Symbol` object represents a node in a AST and the edges
to it's children. In order to simplify mass symbol generation,
:meth:`symbols` can be used.

The :class:`State` object represents a list (set) of equations.
It basically maps variables to their expressions. The :class:`State`
represents a (traditional) DDA file. From a python perspective, a
:class:`State` is not much more then a dictionary on stereoids.
"""

# all "batteries included":
import os, sys, pprint, collections, types, warnings, itertools
identity = lambda x:x
flatten = lambda l: [item for sublist in l for item in sublist]
unique = lambda l: list(set(l))

#class VariableTaxonomy(types.SimpleNamespace):
#    """
#    Inherits __init__, __repr__ and __eq__ from the SimpleNamespace.
#    """
#    def __repr__(self):
#        items = ((f"{k}={v!r}" if len(v)<5 else "%s={%s}"%(k,f"{len(v)} items")) for k,v in self.__dict__.items())
#        return f"VariableTaxonomy({', '.join(items)})"
#    def verbose(self):
#        "Calls the 'traditional' full dump of the namespace object"
#        return types.SimpleNamespace.__repr__(self)
#    def __init__(self, **kwargs):
#        super(**kwargs)
#    def find(self,variable):
#        candidates = [ cls for cls, members in self.__dict__.items() if variable in members ]
#        if len(candidates) == 1: return candidates[0]
#        elif len(candidates) > 1:  raise KeyError(f"Variable {variable} is in multiple classes: {candidates}")
#        else: raise KeyError(f"Variable {variable} is in no class.")
#    __contains__ = find
#    def __len__(self): return len(self.__dict__)
#    def where_is(self, variable_list=[]):
#        if not variable_list:
#            variable_list = [ var for var in members for cls,members in self.__dict__.items() ]
#        return { var: self.find(var) for var in variable_list }
#        
    
class Symbol:
    """
    A symbol is similar to a LISP atom which has a Head and a tail,
    where tail is a list. Common notations for such a type are
    
    * ``head[tail]`` in Mathematica,
    * ``(head, tail)`` in Lisp
    * ``head(tail)`` in C-like languages like Python, Perl, Fortran, C
    * Actually ``[head, *tail]`` in Python, but we don't use that.
      
    A symbol also represents a vertex (node) and it's childs in an ordered tree.
    Think of head being the vertex and tail the (edge) list of children.
    We use the Symbol class to represent the abstract syntax tree (AST)
    of the DDA language for describing ODEs and circuitery.
    
    When you call ``str()`` or similar on instances of this class, it will
    print its representation in the C-like notation. This notation is
    identical to the "classical" DDA language.
    
    There are two types of Symbols: **Variables** have no tail, they just
    consist of a head:
    
    >>> x = Symbol("x")
    >>> print(x)
    x
    >>> x.head
    'x'
    >>> x.tail
    ()
    
    In contrast, **Terms** have a tail:
    
    >>> f = Symbol("f", Symbol("x"), Symbol("y"))
    >>> print(f)
    f(x, y)
    >>> f.head
    'f'
    >>> f.tail
    (x, y)
    
    Variables can be used to create complex expressions for which
    they then serve for as a head:
    
    >>> f,x,y,z = Symbol("f"), Symbol("x"), Symbol("y"), Symbol("z")
    >>> f(x,y)
    f(x, y)
    >>> # example for kind of nonsensical terms:
    >>> x(x,f,x)
    x(x, f, x)
    
    Calling a symbol will always replace it's tail:
    
    >>> f(x)(y)
    f(y)
    
    Symbols are equal to each other if their head and tail equals:
    
    >>> a1, a2 = Symbol("a"), Symbol("a")
    >>> a1 == a2
    True
    >>> f(x) == f(x)
    True
    >>> f(x) == f(x,x)
    False
    
    Symbols can be used as dictionary keys, since they hash trivially
    due to their unique canonical (pythonic) string interpretation.
    
    .. note::
    
       In order to avoid confusion between Python Strings and Symbols,
       you should
       
       - *always* use strings as Symbol heads but
       - *never* use strings in Symbol tails. Instead, use there Symbols
         only.
       
       Think of Symbol implementing the following type (hint):
       ``Tuple[str, List[Symbol]]``.
       
       The DDA code helps you to follow this guide. For instance, the
       representation of ``f1`` shows that it is a symbol with two
       string arguments, while ``f2`` has symbol arguments:

       >>> f1 = Symbol("f", "x", "y")
       >>> f2 = Symbol("f", Symbol("x"), Symbol("y"))
       >>> f1
       f('x', 'y')
       >>> f2
       f(x, y)
       
       And DDA prevents you from shooting in your foot:
       
       >>> f, x, y = symbols("f,x,y")
       >>> f3 = Symbol(f,x,y)
       Traceback (most recent call last):
       ...
       TypeError: Trying to initialize Symbol f(x, y) but head f is a Symbol, not a String.
       
       In previous versions of DDA, the thin line between strings and
       symbols hasn't been made so clear and tracing errors was harder.
       
       Summing up, it is a good convention to *only* have Symbols and
       floats/integers being part of the Symbol tails.
    
    """
    def __init__(self, head, *tail):
        self.head = head
        self.tail = tail
        if isinstance(self.head, Symbol):
            raise TypeError(f"Trying to initialize Symbol {self} but head {head} is a Symbol, not a String.")
    def __call__(self, *tail):
        return Symbol(self.head, *tail)
    def __str__(self):
        # could also say tailstr = str(self.tail) if len(self.tail) else ""
        tailstr = "(" + ", ".join(map(repr, self.tail)) + ")" if self.tail else ""
        return str(self.head) + tailstr
    def __eq__(self, other):
        "Allow to check Symbol('foo') == Symbol('foo')"
        if isinstance(other, str): raise TypeError(f"Convert the string '{other}' to Symbol before comparing with Symbol {self}.")
        if not isinstance(other, Symbol): raise TypeError(f"{other} is not a Symbol. Cannot compare with {self}.")
        return self.head == other.head and \
               len(self.tail) == len(other.tail) and \
               all([ a == b for a,b in zip(self.tail,other.tail) ])
    def __hash__(self):
        "Allows Symbol to be keys in dictionaries. Also used for persistence"
        return hash(str(self)) # not yet tested, use with care
    
    # Don't repeat yourself:
    __repr__ = __str__
    __unicode__ = __str__
    
    # Here is a bunch of handy methods to deal with Symbols:
    
    def is_variable(self):
        "A variable is a symbol without a tail."
        return len(self.tail) == 0

    def is_term(self):
        "A term is a symbol with a tail."
        return not self.is_variable()

    def variables(self):
        """Compute the direct dependencies of this symbol, i.e. other variables
        directly occuring in the tail."""
        return [ el for el in self.tail if el.is_variable() ]

    def all_variables(self):
        "Like :meth:variables, but also find variables in all children."
        return unique(flatten([ el.all_variables() if el.tail else [el] for el in self.tail if is_symbol(el)]))

    def all_terms(self):
        "Like :meth:all_variables, but for terms: Returns a list of all terms in all children of this node."
        return unique(flatten([ el.all_terms() for el in self.tail if is_symbol(el) and el.tail() ]))

    def map_heads(self, mapping):
        """
        Call a mapping function on all heads in all (nested) subexpressions.
        The mapping is effectively carried out on the head (ie. maps strings)
        Returns a new mapped Symbol. This routine is suitable for renaming
        variables and terms within an AST. Example usage:
        
        >>> Symbol("x", Symbol("y"), 2).map_heads(lambda head: head+"foo")
        xfoo(yfoo, 2)

        The mapping is unaware of the AST context, so you have to distinguish
        between variables and terms yourself if you need to. See also
        :meth:`map_variables` for context-aware head mapping. Compare these
        examples to the ones given for :meth:`map_variables`:

        >>> x, map, y = Symbol("x"), lambda _: "y", Symbol("y")
        >>> x.map_heads(map) == x.map_variables(map)  # == y
        True
        >>> x(x,x).map_heads(map) == y(y,y)
        True
        >>> x(x, x(x)).map_heads(map) == y(y, y(y))
        True
        """
        return Symbol(mapping(self.head), *[(el.map_heads(mapping) if is_symbol(el) else el) for el in self.tail])

    def map_variables(self, mapping, returns_symbol=False):
        """
        Calls a mapping function on all variables within the (nested) subexpressions.
        The mapping is effectively carried out on the head (ie. maps strings). This
        is a mixture between :meth:`map_heads` and :meth:`map_tails`.
        
        variables but not terms within the AST. Examples:
        
        >>> x.map_variables(map) == y
        True
        >>> x(x,x).map_variables(map) == x(y,y)
        True
        >>> x(x, x(x)).map_variables(map) == x(y, x(y))
        True
        
        This function ignores non-symbols as they cannot be variables. This is the same
        as :meth:`map_tails` does and is handy when you have numbers within your expressions:
        
        >>> x = Symbol("x")
        >>> expr = x(123, x(9.1), x, x(x, 0.1, x))
        >>> res1 = expr.map_variables(lambda xx: "y")
        >>> res2 = expr.map_variables(lambda xx: Symbol("y"), returns_symbol=True)
        >>> res1 == res2
        True
        >>> res1
        x(123, x(9.1), y, x(y, 0.1, y))
        
        If you want to use ``map_variables`` to change a variable to a term, and/or if your
        mapping function does not return strings but Symbols, use ``returns_symbol=True``:
        
        >>> Symbol("x").map_variables(lambda x: Symbol("y", 123), returns_symbol=True)
        y(123)
        >>> Symbol("x").map_variables(lambda x: Symbol("y", 123)) # this won't work
        Traceback (most recent call last):
        ...
        TypeError: Trying to initialize Symbol y(123) but head y(123) is a Symbol, not a String.
        """
        return (identity if returns_symbol else Symbol)(mapping(self.head)) if self.is_variable() else \
            Symbol(self.head, *[ (el.map_variables(mapping, returns_symbol=returns_symbol) if is_symbol(el) else el) for el in self.tail ])

    def map_tails(self, mapping, map_root=False):
        """
        Calls a mapping function on all tails in all (nested) subexpressions.
        The mapping is carried out on the tail symbols (ie. maps Symbols).
        Returns a new mapped Symbol. The routine is suitable for AST walking,
        adding/removing stuff in the tails while preserving the root symbol.
        This could also be called ``map_symbols``, c.f. :meth:`map_terms`.
        
        Example for recursively wrapping all function calls:
        
        >>> x,y,z = symbols("x,y,z")
        >>> x(y, z(x), x(y)).map_tails(lambda smb: Symbol("foo")(smb))
        x(foo(y), foo(z(foo(x))), foo(x(foo(y))))
        >>> x(y, z(x), x(y)).map_tails(lambda smb: Symbol("foo")(smb), map_root=True)
        foo(x(foo(y), foo(z(foo(x))), foo(x(foo(y)))))
        
        Example for recursively removing certain unary functions ``z(x)`` for any ``x``:

        >>> remover = lambda head: lambda x: x.tail[0] if isinstance(x,Symbol) and x.head==head else x
        >>> x(y, z(x), x(z(y),x)).map_tails(remover("z"))
        x(y, x, x(y, x))
        
        The argument ``map_root`` decides whether the map is run on the
        root node or not. It will be ``map_root=False`` in any recursive use.
        In former instances of this code, it was always ``map_root=False``.
        Example:
        
        >>> (a, b), flip = symbols("a,b"), lambda smb: b if smb.head==a.head else a
        >>> a(b,a,b).map_tails(flip, map_root=True)
        b
        >>> a(b,a,b).map_tails(flip, map_root=False)
        a(a, b, a)
        
        Note how the ``flip`` function cuts every tail and returns variables only.
        Here is a variant which perserves any tail:
        
        >>> (a, b) = symbols("a,b")
        >>> flipper = lambda smb: (b if smb.head==a.head else a)(*smb.tail)
        >>> a(b,a,b).map_tails(flipper, map_root=True)
        b(a, b, a)
        
        Here is another example which highlights how ``map_tails`` can convert terms to
        variables:
        
        >>> x, map, y = Symbol("x"), lambda _: Symbol("y"), Symbol("y")
        >>> x(x,x).map_tails(map, map_root=False)
        x(y, y)
        >>> x(x,x(x,x)).map_tails(map, map_root=False)
        x(y, y)
        
        For real-life examples, study for instance the source code of :mod:`cpp_exporter` or
        grep any DDA code for ``map_tails``.
        
        See also :meth:`map_heads` and :meth:`map_variables` for variants.
        """
        r = Symbol(self.head, *[(mapping(el.map_tails(mapping, map_root=False)) if is_symbol(el) else el) for el in self.tail])
        return mapping(r) if map_root else r

    def map_terms(self, mapping, returns_symbol=False):
        """
        Calls a mapping function on all terms within the (nested) subexpressions.
        The mapping is effectively carried out on the term head (ie. maps strings).
        See :meth:`map_variables` for the similar-minded antoganist as well as
        :meth:`map_heads` and :meth:`map_tails` for more *low level* minded variants.
        
        Returns a new mapped Symbol. This routine is suitable for renaming terms
        but not variables within the AST. Examples:
        
        >>> x.map_terms(map) == x
        True
        >>> x(x,x).map_terms(map) == y(x,x)
        True
        >>> x(x, x(x)).map_terms(map) == y(x, y(x))
        True
        
        This function ignores non-symbols as they cannot be variables, similar to
        :meth:`map_variables`.
        
        It is basically ``map_terms(map) = map_tails(lambda smb: Symbol(map(smb.head))(smbl.tail) if symb.is_variable() else smb)``.
        
        The argument ``returns_symbol`` allows to discriminate between mapping
        functions which return strings (for symbol heads) or symbols. The later
        allows for manipulating expressions.
        """
        if self.is_variable(): return self
        r = mapping(self.head)
        if not returns_symbol: r = Symbol(r)
        return r(*[ (el.map_terms(mapping) if is_symbol(el) else el) for el in self.tail ])
    
    def draw_graph(self, graph=None):
        """
        Uses graphviz to draw the AST down from this symbol.
        
        See also :method:`State.draw_dependency_graph` for similar draph
        drawing code and notes on python library dependencies.
        
        .. note:: This method constructs the graph by drawing edges between similar
           named symbols. This will *not* represent the abstract syntax tree if
           a single symbol head, regardless of whether variable or term, appears
           twice.
           
           If you want to draw the actual AST with this function, you have to make
           each symbol (head) unique by giving them distinct names.
        
        Simple usage example:
        
        >>> expression = x(1,y,2,z(3,4))
        >>> graph = expression.draw_graph()
        >>> print(graph)  # doctest: +NORMALIZE_WHITESPACE
        digraph "DDA-Symbol" {
           node [shape=doublecircle]
           x
           node [shape=circle]
           x -> 1
           y
           x -> y
           x -> 2
           z
           z -> 3
           z -> 4
           x -> z
        }
        >>> graph.view() # Call this to draw the graph  # doctest: +SKIP
        """
        if not graph:
            from graphviz import Digraph
            graph = Digraph("DDA-Symbol", filename="symbol.gv")
            #graph.attr(size="8,5")  # rankdir="LR"
        
        graph.attr("node", shape="doublecircle")
        graph.node(self.head)
        
        graph.attr("node", shape="circle")
        for k in self.tail:
            if is_symbol(k):
                graph = k.draw_graph(graph)
                graph.edge(self.head, k.head)
            else:
                graph.edge(self.head, str(k))

        return graph

def is_symbol(smbl):
    "Convenience function"
    return isinstance(smbl, Symbol)

#def convenience(method, cls=Symbol):
#    def fun(inst):
#        if not isinstance(inst, cls): raise ValueError(f"'{str(inst)}' is of type {type(inst)} but {str(cls)} required")
#        return getattr(inst, method)
#    return fun
#map_heads = convenience("map_heads")

def symbols(*query):
    """
    Quickly make symbol objects. Usage similar to sympy's symbol function:
    
    >>> a, b = symbols("a", "b")
    >>> x, y, z = symbols("x, y, z")
    """
    symbols = [ Symbol(p.strip()) for q in query for p in q.split(",") ]
    return symbols if len(symbols)>1 else symbols[0]

def topological_sort(dependency_pairs):
    """Sort a graph (given as edge list) subject to dependency constraints.
    The result are two lists: One for the sorted nodes, one for the
    unsortable (cyclically dependent) nodes.
    
    Implementation shamelessly stolen from https://stackoverflow.com/a/42359401
    """
    num_heads = collections.defaultdict(lambda: 0)   # num arrows pointing in
    tails = collections.defaultdict(list)      # list of arrows going out
    heads = []                     # unique list of heads in order first seen
    for h, t in dependency_pairs:
        num_heads[t] += 1
        if h in tails:
            tails[h].append(t)
        else:
            tails[h] = [t]
            heads.append(h)

    ordered = [h for h in heads if h not in num_heads]
    for h in ordered:
        for t in tails[h]:
            num_heads[t] -= 1
            if not num_heads[t]:
                ordered.append(t)
    cyclic = [n for n, heads in num_heads.items() if heads]
    Results = collections.namedtuple('topological_sort', ['sorted', 'cyclic'])
    ordered, cyclic = ordered[::-1], cyclic[::-1] # I find reverse logic more logical
    return Results(ordered, cyclic)

class State(collections.UserDict):
    """
    A state is a dictionary which is by convention a mapping from variable
    names (as *strings*) to their symbolic meaning, i.e. a ``Symbol()``.
    We refer to the keys in the dictionary as the *Left Hand Side* (LHS) and
    the values in the dictionary as the *Right Hand Side* (RHS), in analogy
    to an Equation.

    .. note::

       Since ``Symbol()`` spawns an AST, a state is a list of variable definitions.    
       A DDA file is a collection of equations. Therefore, a state holds the
       content of a DDA file.
   
    This class collects a number of basic helper routines for dealing with
    states.
    
    In order to simplify writing DDA files in Python, this class extends
    the dictionary idiom with the following optional features, which
    are turned on by default (but can be disabled by constructor
    arguments ``type_peacemaking`` and ``default_symbol``).
    
    - **Type peacemaking**: Query a ``Symbol()``, get translated to
      ``str()``:
      
      >>> State({ "foo": Symbol("bar")})[Symbol("foo")] == Symbol("bar")
      True
      
    - **Default Symbol**: Automatically add an entry when unknown:
      
      >>> State()["foo"] == Symbol("foo")
      True
    
    .. note::
    
       By intention, the keys of the State are always strings, never Symbols.
       This also should make sure you don't use complex ASTs for keys,
       such as ``Symbol("foo", "bar")``.
    
    As ``State``` extends ``collections.UserDict``, you can access the underlying
    dictionary:
    
    >>> x,y = symbols("x,y")
    >>> add, integrate = symbols("add", "integrate")
    >>> eqs = { x: add(y,y), y: integrate(x) }
    >>> state = State(eqs); print(state)
    State({'x': add(y, y), 'y': integrate(x)})
    >>> state.data
    {'x': add(y, y), 'y': integrate(x)}
    
    .. warning::

       Don't be fooled by refering to the state while constructing the state. This
       will end up in overly complex expressions. By rule of thumb, only use
       :class:`Symbols` at the state definition (or in particular on the right hand side).
       For instance, you do want to construct a state like
       
       >>> state = State()
       >>> state["x"] = Symbol("add", Symbol("y"), Symbol("y"))
       >>> state["y"] = Symbol("int", Symbol("x"))
       >>> state
       State({'x': add(y, y), 'y': int(x)})
       
       In contrast, this is most likely not what you want:
       
       >>> state["y"] = Symbol("int", state["x"])
       >>> state
       State({'x': add(y, y), 'y': int(add(y, y))})
       
       This time, you did not exploit the definition of ``state["x"]`` by referencing
       on DDA level but instead inserted the expression by referencing on Python
       level. This is like *compile-time evaluation* versus *runtime evaluation*,
       when *compile-time* is at python and *runtime* is when evaluating the DDA
       expressions in some time evolution code.
       
       Summing up, the mistake above is to reference to ``state`` while constructing7
       the ``state``. You should not do that. You go best by defining the ``Symbol``
       instances before and then only using them all over the place:
       
       >>> x, y, add, int = symbols("x, y, add, int")
       >>> state[x] = add(y,y)
       >>> state[y] = int(x)
       >>> state

    .. note::

       Why the name? The class name ``State`` seems arbitrary and quirky, ``System`` may
       be a better choice (given that the class instances hold an equation system).
       However, one could also argue that the class instances hold the definition for
       a system in a particular state. especially, ``State.keys()`` are the state
       variables which undergo a definition by their corresponding ``State.values()``.
       Most CAS do not have a special class for collections of equations. Instead,
       they typically have some *equation* type and equation systems are sets or lists
       of equations. In PyDDA, we don't have an equation type because the DDA domain
       specific language (see :mod:`dsl`) doesn't provide advanced treatments of
       equations but is basically only a lengthy definition of a set of equations,
       which you could understand as a mapping/dictionary data type defining the state
       of the system. That's why `State` is actually an enriched `dict`.

    """    
    def __init__(self, initialdata=dict(), type_peacemaking=True, default_symbol=True):
        self.type_peacemaking = type_peacemaking
        self.default_symbol = default_symbol
        super().__init__(initialdata)
        
    @classmethod
    def from_string(cls, *string_or_list_of_strings):
        "Shorthand for :meth:`.dsl.read_traditional_dda`. Returns new instance."
        from .dsl import read_traditional_dda
        return read_traditional_dda("\n".join(string_or_list_of_strings))
    
    def to_string(self):
        "Shorthand for :meth:`.dsl.to_traditional_dda`. Returns a string (representation)."
        from .dsl import to_traditional_dda
        return to_traditional_dda(self, cleanup=False, prefix=None, suffix=None)

    def __getitem__(self, name):
        if self.type_peacemaking and isinstance(name, Symbol):
            name = name.head
        if self.default_symbol and not name in self.data:
            self.data[name] = Symbol(name)
        return super().__getitem__(name) # or just self.data[name] 
    
    def __setitem__(self, name, value):
        return super().__setitem__(name, value)
            
    def update(self, other_dict):
        if isinstance(other_dict, State):
            other_dict = other_dict.data
        return super().update(other_dict) # always returns None!
    
    def __repr__(self):
        return "{}({})".format(type(self).__name__, pprint.pformat(self.data))
        
    # Actual useful methods come here:
    
    def equation_adder(self):
        """
        Syntactic sugar for adding new equations to the system. Usage:

        >>> state = State()
        >>> x,y,z,add,int = symbols("x,y,z,add,int")
        >>> eq = state.equation_adder()
        >>> eq(y=int(x))
        >>> eq(x=add(y,z), z=int(x,0,0.1))
        >>> state
        State({'x': add(y, z), 'y': int(x), 'z': int(x, 0, 0.1)})
        
        Known limitations: This doesn't work any better then the ``BreveState``
        below because keywords must not be variables, they will always resolve
        to their string representation.
        
        >>> foo = Symbol("bar")
        >>> s1, s2 = State(), State()
        >>> eq1 = s1.equation_adder()
        >>> eq1(foo=42)
        >>> s2[foo] = 42
        >>> s1
        State({'foo': 42})
        >>> s2
        State({'bar': 42})
        """
        def adder(**dct):
            for k,v in dct.items():
                self[k] = v
        return adder
        
    def map_tails(self, mapper, map_root=True):
        "Apply :meth:`Symbol.map_tails` on all right hand sides."
        apply_mapper = lambda el: el.map_tails(mapper, map_root=map_root) if el.is_term() else mapper(el)
        return State({var: apply_mapper(self[var]) for var in self })
    
    def map_heads(self, mapper):
        """
        This function is suitable for renaming variables.
        mapper is always executed on the string variable names (Symbol heads)
        """
        return State({ mapper(var): self[var].map_heads(mapper) for var in self })
    
    def symbols(self, *query):
        "Same as :meth:`symbols()` above, but register at self (state)"
        return [ self[x.head] for x in symbols(*query) ]
        
    def constant_validity(self):
        """
        Check validity of numeric constants in the state.
        Depending on context, values -1 < t < +1 are illegal.
        
        (Not yet implemented!)
        """
        pass # TBD, probably not here
    
    def dependency_graph(self):
        """
        Returns the edge list of the variable dependency graph of this state.
        We can call :meth:`topological_sort()` on the result of this method.
        
        A weird example including some corner cases:

        >>> s1 = State.from_string("foo = const(0.7)", "bar=mult(bar,baz)", "baz=f(bar)")
        >>> s1.dependency_graph()
        [('bar', 'bar'), ('bar', 'baz'), ('baz', 'bar')]
        
        Another even more weird example which exploits raw value assignment,
        something which is not following the ``foo=call(bar)`` requirement for DDA
        files:
        
        >>> a,b,c,d,f = symbols("a,b,c,d,f")
        >>> s2 = State({ a: f(0.7), b: c(b,b), c: 42, d: c })
        >>> s2.dependency_graph()
        [('b', 'b'), ('d', 'c')]
        
        Note that this function always returns list of tuples of strings. No more
        symbols. See also :meth:`draw_dependency_graph` for a quick way of exporting
        or plotting this graph.
        """
        def dependent(rhs):
            if not isinstance(rhs, Symbol): return []
            if rhs.is_term(): return map(str, rhs.all_variables())
            if rhs.is_variable(): return [str(rhs)]
        # Comptue adjacency list of dependencies. All is strings, no more symbols
        adjacency_list = { k: list(dependent(self[k])) for k in self }
        edge_list = [ (k,e) for k,dep in adjacency_list.items() for e in dep ]
        return sorted(edge_list)
        
    def draw_dependency_graph(self, export_dot=True, dot_filename="test.dot"):
        """
        If you have ``networkx`` and ``pyGraphViz`` installed, you can use this method
        to draw the *variable dependency graph* (see method :meth:`dependency_graph()`)
        with ``Dot``/``GraphViz``. This method will return the ``nx.DiGraph()`` instance.
        If ``export_dot`` is set, it will also write a dotfile, call ``dot`` to
        render it to a bitmap and open that bitmap.
        
        .. note::
        
           Your distribution package ``python-graphviz`` is probably not ``pygraphviz``.
           You are on the safe side if you run: ``pip install pygraphviz``
        """
        import networkx as nx
        from networkx.drawing.nx_agraph import graphviz_layout
        
        G = nx.DiGraph()
        G.add_edges_from(self.dependency_graph())
        if export_dot:
            nx.nx_agraph.write_dot(G, dot_filename)
            print(f"Exported to {dot_filename}. Running dot...")
            os.system(f"dot -Tpng {dot_filename} > {dot_filename}.png && open {dot_filename}.png")
        return G

    def name_computing_elements(self, strict=False):
        """
        Name all computing elements / intermediate expressions. Returns
        a new State which is *linearized* in a way that the numbering proposes a computing
        order.
        
        Linearization is an idempotent operation, i.e. for any
        ``lin = state.name_computing_elements()`` it is ``lin == lin.name_computing_elements()``.
        Mathematically, it is a projection of the state on its linearized one.
        
        Linearization means to define a evaluation order and to give unique names
        to all terms occuring. Note that *all* depends on the strictness (`strict=True` vs.
        the default `strict=False`):
        
        >>> x, y, sum, mult = symbols("x, y, sum, mult")
        >>> ns = State({ x: sum(x,y, sum(y, mult(y,x))), y: mult(x) })
        >>> print(ns.name_computing_elements().to_string())              # doctest: +NORMALIZE_WHITESPACE
        mult_1 = mult(y, x)
        sum_1 = sum(y, mult_1)
        x = sum(x, y, sum_1)
        y = mult(x)
        >>> print(ns.name_computing_elements(strict=True).to_string())   # doctest: +NORMALIZE_WHITESPACE
        mult_2 = mult(x)
        sum_2 = sum(x, y, sum_1)
        x = sum_2
        y = mult_2
        
        Here, *strict* means that really *every* term is labeled, even if this yields in "dumb"
        assignments such as ``x = sum_2``. You want a strict naming when enumerating computing elements,
        while a non-strict naming is preferable for brief evaluation. Also note that
        
        >>> x,const = symbols("x,const")
        >>> State({ x: const(42) }).name_computing_elements(strict=False)
        State({'x': const(42)})
        >>> State({ x: const(42) }).name_computing_elements(strict=True)
        State({'const_1': const(42), 'x': const_1})
        >>> State({ x: const(42) }).name_computing_elements(strict=True).name_computing_elements(strict=True) # doctest:+ELLIPSIS +NORMALIZE_WHITESPACE +SKIP
        /.../ast.py:819: UserWarning: State.named_computing_elements(): While counting const, I notice that  const_1 is already part of the state. Maybe you want to run name_computing_elements(strict=False) for idempotence.
          warnings.warn(...)
        State({'const_1': const_1_, 'const_1_': const(42), 'x': const_1})
        
        On this mini example, one especially sees that idempotence is only given when ``strict=False``. 
        It is ``state.name_computing_elements(True) == state.name_computing_elements(s[0]).name_computing_elements(s[1])....name_computing_elements(s[n])``
        when ``s`` is a boolean array of ``len(s)==n`` and ``sum(s) == 1``, i.e. only one
        occurence of ``strict=True`` and all other ``False``.

        The linearized state only has entries of a *normal form* ``state[f_i] = f(v1,v2,...)``
        for a function (term) ``f`` and some variables ``v_j``. Furthermore note how even
        ``x = sum_2`` in the above example indirects the former assignment of ``x = sum(x,y...)``.
        Again, for any value in the linearized state, the tail only contains variables,
        no terms. This is handy for many things, such as circuit drawing, imperative
        evaluation (in combination with :meth:`variable_ordering`, cf. :mod:`cpp_exporter`)
        and determination of integrands/actual variables. For instance
        
        >>> s = State.from_string("foo = const(0.7)", "baz=mult(bar,bar)", "bar = neg(int(neg(baz), foo, 0.3))")
        >>> s
        State({'bar': neg(int(neg(baz), foo, 0.3)), 'baz': mult(bar, bar), 'foo': const(0.7)})
        >>> print(s.name_computing_elements(strict=True).to_string())    # doctest: +NORMALIZE_WHITESPACE
        bar = neg_2
        baz = mult_1
        const_1 = const(0.7)
        foo = const_1
        int_1 = int(neg_1, foo, 0.3)
        mult_1 = mult(bar, bar)
        neg_1 = neg(baz)
        neg_2 = neg(int_1)
        >>> print(s.name_computing_elements(strict=False).to_string())    # doctest: +NORMALIZE_WHITESPACE
        bar = neg(int_1)
        baz = mult(bar, bar)
        foo = const(0.7)

        Here one sees immediately that ``int_1`` is the actual integral solution while
        ``bar`` is only a derived quantity. Calls like ``const(float)`` remain unchanged
        since they are already in the normal form ``f(v1,v2,...)``.
        
        Here is a more complex example:
        
        >>> from dda.computing_elements import neg,int,mult
        >>> dda_state = State({"x": neg(int(neg(int(neg(mult(1, Symbol("x")), 0.005, 1)), 0.005, 0))) })
        >>> dda_state.name_computing_elements().variable_ordering().where_is # doctest: +NORMALIZE_WHITESPACE
        {'x': 'vars.aux.sorted',
        'mult_1': 'vars.aux.sorted',
        'neg_2': 'vars.aux.sorted',
        'neg_1': 'vars.aux.sorted',
        'int_1': 'vars.evolved',
        'int_2': 'vars.evolved'}

        """
        symbol_counter = collections.defaultdict(lambda:0)
        intermediates = {}
        #NamedSymbol = collections.namedtuple('NamedSymbol', ['symbol', 'number'])
        def register_computing_element(el):
            if len(el.tail) and not el in intermediates:
                symbol_counter[el.head] += 1
                named_symbol = "%s_%d" % (el.head, symbol_counter[el.head])
                while named_symbol in self:
                    # Note that this naming clash could also happen by pure coincidence.
                    warnings.warn(f"State.named_computing_elements(): While counting {el.head}, I notice that  {named_symbol} is already part of the state. Maybe you want to run name_computing_elements(strict=False) for idempotence.")
                    named_symbol += "_"
                intermediates[named_symbol] = el
                return Symbol(named_symbol)
            return el
        
        linearized_state = self.map_tails(register_computing_element, map_root=strict)
        linearized_state.update(intermediates)
        return linearized_state
    
    def variable_ordering(self):
        """
        Will perform an analysis of all variables occuring in this state (especially in the RHS).
        This is based on the linarized variant of this state (see :meth:`name_computing_elements()`).
        
        The return value is an object (actually a types.SimpleNamespace instance) which contains
        lists of variable names (as strings). The properties (categories) are primarly
        
          * explicit constants: Any entry ``state["foo"] = const(1.234)``
          * State variables/evolved variables: Any outcome of a time integration, i.e. ``int(...)``,
            i.e. ``Symbol("int")``. This can be as simple as ``state["foo"] = int(Symbol("foo"),...)``.
            Complex terms such as ``state["foo"] = mult(int(foo), int(bar))`` will result in
            intermediate variables called like ``int_0``, ``int_1``(see :meth:`name_computing_elements`
            for the code which invents these names), which are the actual evolved variables.
          * Auxilliary variables: Any other variables which are required to compute evolved
            variables.
            
        By intention, we **sort only the aux. variables**. One should check that they DO NOT have any
        cyclic dependency, because feedback loops are only useful on integrators at this level
        of circuit modeling.

        We differntiate the auxilliaries further into:
        
          * ``sorted_aux_vars``: Auxilliaries required to compute the state variable changes
          * ``cyclic_aux_vars``: Auxilliaries which have a cyclic dependency on each other
            (this should not happen as it won't lead to a stable circuit)
          * ``unneeded_auxers``: Auxilliaries which are not required to compute the state.
            These are probably used in postprocessing. If they depend on the state variables,
            further work is neccessary.
        
        Given an ODE problem ``dq/dt = f(q)``, an imperative code for evolving the state
        ``q`` in time  should compute all auxillairy variables in the respective order
        before computing the actual ``dq/dt``. The dependency is basically, in pseudo code:
        
        >>> aux  = function_of(aux, state)                    # doctest: +SKIP
        >>> dqdt = function_of(aux, state)                    # doctest: +SKIP
        
        and in the numerical integration schema step
        
        >>> state = function_of(dqdt)                         # doctest: +SKIP
        
        This method returns a namespace object, which is basically a fancy dictionary. It
        is used over a simple dictionary just for shorter syntax.
        
        The following examples demonstrate a deeply nested corner case, i.e. a 
        compute graph consisting of a single "long" Euler cycle. By breaking up this
        cycle at the integrations, ``variable_ordering()`` can linearize these cycles
        correctly. This works both for non-strict and strict element naming.
        
        >>> # variable ordering is made based on non-strict naming:
        >>> dda_state.name_computing_elements(strict=False) # doctest: +NORMALIZE_WHITESPACE
        State({'int_1': int(neg_1),
               'int_2': int(neg_2),
               'mult_1': mult(1, x),
               'neg_1': neg(mult_1, 0.005, 1),
               'neg_2': neg(int_1, 0.005, 0),
               'x': neg(int_2)})
        >>> dda_state.name_computing_elements(strict=True) # doctest: +NORMALIZE_WHITESPACE
        'int_2': int(neg_2),
        'mult_1': mult(1, x),
        'neg_1': neg(mult_1, 0.005, 1),
        'neg_2': neg(int_1, 0.005, 0),
        'neg_3': neg(int_2),
        'x': neg_3})
        >>> # This is how the full output looks like
        >>> dda_state.variable_ordering()  # doctest: +SKIP
        namespace(aux=namespace(all=['mult_1', 'neg_1', 'neg_2', 'x'],
                                sorted=['x', 'mult_1', 'neg_2', 'neg_1'],
                                cyclic=[],
                                unneeded=set()),
                  evolved=['int_1', 'int_2'],
                  explicit_constants=[],
                  all=['int_1', 'int_2', 'mult_1', 'neg_1', 'neg_2', 'x'],
                  ordering=OrderedDict([('vars.explicit_constants', []),
                                        ('vars.aux.sorted',
                                         ['x', 'mult_1', 'neg_2', 'neg_1']),
                                        ('vars.aux.cyclic', []),
                                        ('vars.evolved', ['int_1', 'int_2']),
                                        ('vars.aux.unneeded', set())]),
                  where_is={'int_1': 'vars.evolved',
                            'int_2': 'vars.evolved',
                            'mult_1': 'vars.aux.sorted',
                            'neg_1': 'vars.aux.sorted',
                            'neg_2': 'vars.aux.sorted',
                            'x': 'vars.aux.sorted'})
        >>> # Compare the strict and nonstrict orderings:
        >>> for k, v in dda_state.variable_ordering().ordering.items(): print(f"{k:25s}: {v}")
        vars.explicit_constants  : []
        vars.aux.sorted          : ['x', 'mult_1', 'neg_2', 'neg_1']
        vars.aux.cyclic          : []
        vars.evolved             : ['int_1', 'int_2']
        vars.aux.unneeded        : set()
        >>> for k, v in dda_state.name_computing_elements(strict=True).variable_ordering().ordering.items(): print(f"{k:25s}: {v}")
        vars.aux.sorted          : ['neg_3', 'x', 'mult_1', 'neg_2', 'neg_1']
        
        """
        
        lin = self.name_computing_elements(strict=False)
        vars = types.SimpleNamespace() # VariableTaxonomy()
        vars.aux = types.SimpleNamespace() # VariableTaxonomy()
        
        # Thanks to named computing elements, can find all int(...) expressions
        # without searching, since they must come first.
        vars.evolved  = sorted(filter(lambda var: lin[var].head in ["int", "diff"],   lin))
        vars.explicit_constants = sorted(filter(lambda var: lin[var].head == "const", lin))
        
        # Then compute ALL aux variables BEFORE computing dqdt.
        # The order of dqdt should not be touched, as there CAN NOT be any
        # dependency, since dqdt.foo = int(lin).

        vars.all = sorted(set.union(*[set(map(str, lin[k].all_variables())) for k in lin], set(lin)))
        vars.aux.all = [ v for v in vars.all if not v in vars.evolved and not v in vars.explicit_constants  ]

        # Linearize aux expressions by dependency sorting.
        dep_edge_list = lin.dependency_graph()
        # Edge direction: (a,b) = [dependent variable a]--[depends on]-->[dependency b]
        aux_dep_edges = [ (a,b) for a,b in dep_edge_list if b in vars.aux.all ]  # a in vars.aux.all and
        # Note:
        #  In order to detect auxers dependent on evolved vars, we only filter for b being aux and let a free.
        #  We cannot feed the dep_edge_list directly into topological_sort since it would yield cyclic
        #  dependency for any meaningful compute circuit, which always has feedback loops.
        #  These loops are breaked by the evolved variables.
        dep_sorted, dep_cyclic = topological_sort(aux_dep_edges)
        vars.aux.sorted = [ v for v in dep_sorted if v in vars.aux.all ]
        vars.aux.cyclic = [ v for v in dep_cyclic if v in vars.aux.all ]

        vars.aux.unneeded = set(vars.aux.all) - (set(vars.aux.sorted) | set(vars.aux.cyclic) | set(vars.explicit_constants))
        
        # The following is an evil eval trick to avoid repetition. It just dumps the content of the mentioned vars.[...] arrays.
        ordering = ["vars.explicit_constants", "vars.aux.sorted", "vars.aux.cyclic", "vars.evolved", "vars.aux.unneeded"]
        vars.ordering = collections.OrderedDict([ (name, eval(name, {"vars":vars})) for name in ordering ])
        
        # allow inverse lookup for debugging
        vars.where_is = { i: k  for k,l in vars.ordering.items() for i in l }
        
        # In the end, make sure we haven't lost something!
        residium = set(lin.keys()) - set(vars.where_is.keys())
        if len(residium) != 0: 
                warnings.warn("State.variable_ordering(): Lost these variables in ordering analysis: " + str(residium))
                # could also raise an issue, because it is kind of serious.
        
        return vars

    def remove_duplicates(self):
        """
        Assuming a linearized state, this function simplifies the system by removing/resolving
        duplicate entries. No further renaming takes place: Always the first encounter of a
        term determines the name for all equivalent terms.
        
        Returns a new state.
        """
        # 1. equivalence classes (master_key: [equivalents]} and reverse lookup values {term: master_key}
        classes, values = {}, {}
        for i,iv in self.items():
            if iv in values:
                classes[ values[iv] ].append(i)
            else: # new equivalence class
                classes[ i ] = []
                values[ iv ] = i
        # 2. Create the new state based on the equivalence classes
        #    Start with a straight copy of the master classes
        ret = State({ master: self[master] for master in classes.keys() })
        # 3. Recurvsively replace equivalent class members in any term in ret
        for master, equivalents in classes.items():
            slave2master = lambda varname: master if varname in equivalents else varname
            for testkey, testval in ret.items():
                ret[testkey] = testval.map_variables(slave2master, returns_symbol=False)
        return ret
    
    def term_statistics(self):
        """
        Assuming a linearized state, tells how much each term occurs within this state.
        Returns instance of ``collections.Counter``. For simplicity, the term heads
        are returned as strings.
        
        Typical usage is like ``state.name_computing_elements().remove_duplicates().term_statistics()``
        """
        return collections.Counter([ term.head for term in self.values() ])
        
    def export(self, to, **passed_args):
        "Syntactic sugar for dda.export(), for convenience"
        from . import export
        return export(self, to, **passed_args)

class BreveState(State):
    """
    This subclass of a state adds *syntactic sugar* by allowing attribute/member access
    notation. Instead of ``state["foo"]`` you can write ``state.foo`` on instances of
    this class. Example usage:
    
    >>> x,y,z = symbols("x,y,z")
    >>> state = BreveState()
    >>> state.x = x(y,z)
    >>> state.y = y(x,z)
    >>> state.z = z(x,y)
    >>> print(state)
    BreveState({'x': x(y, z), 'y': y(x, z), 'z': z(x, y)})
    
       Known limitations:
    
       * Breaks Python class introspection (for instance tab completion in *iPython*)
       * Of course users cannot add any non-data related attribute (or method)

    See also :meth:`State.equation_adder()` for similar sugar which might have
    unexpected effects:
    
    >>> s, b = State(), BreveState()
    >>> foo = Symbol("bar") # in this context, similar to a string foo = "bar"
    >>> s[foo] = 42         # foo resolves to string represntation "bar"
    >>> b.foo = 42          # equals b["foo"], thus has nothing to do with variable foo
    >>> s
    State({'bar': 42})
    >>> b
    BreveState({'foo': 42})
    
    If you find this class useful, you also might like ``types.SimpleNamespace`` or
    ``collections.namedtuple``. Both are basically immutable, while this object type
    is mutable by intention.
    """    
    wellknown = ("data", "type_peacemaking", "default_symbol")
    isinternal = lambda name: (not is_symbol(name)) and (name in BreveState.wellknown or name.startswith("_"))

    def __getattr__(self, name):
        return self.__dict__[name] if BreveState.isinternal(name) else self.__getitem__(name)
    def __setattr__(self, name, value):
        if BreveState.isinternal(name):
            self.__dict__[name] = value
        else:
            self.__setitem__(name, value) # or self.data[name] = value

#
#
#
# 
#

"""
This module provides interplay with the SymPy package. SymPy is a lightweight
pure-python computer algebra system which is bundled with SciPy.
An adapter to/from SymPy allows to use powerful Computer Algebra basic
functions such as expression simplification.

We use this currently to provide a lean latex representation for the
cumbersome DDA expressions.
"""

from . import Symbol, dda
import builtins

def from_sympy(sympy_equation_list):
    """
    Import a state from a set of equations from SymPy.
    
    This function expects a python list of sympy equations where there
    is a single sympy symbol on one hand and an expression on the other
    hand (see examples below).
    
    The mapping basically follows the `SymPy key invariant
    <https://docs.sympy.org/latest/tutorial/manipulation.html#args>`_:
    "Every well-formed SymPy expression must either have empty ``args``
    or satisfy ``expr == expr.func(*expr.args)``".

    Therefore we basically map a sympy expression ``(expr.func,
    expr.args)`` to the DDA ``(head, tail)`` notation. While the heads
    are easy to map (for instance, ``sympy.Mul`` equals ``dda.mult``),
    special attention must be given to the tails, for instance SymPys
    ``Mul(a,b,c)`` translates to DDAs ``mult(mult(a,b),c)`` (in DDA
    we always assume commutative real-valued variables). Also in DDA
    there is ``neg(x)`` or ``div(x,y)`` which is represented in SymPy as
    ``Mul(Integer(-1), x)`` and ``Mul(Symbol('x'), Pow(Symbol('y'), Integer(-1)))``,
    respectively.
    """
    import sympy
    
    raise ValueError("Not yet implemented!")

    # This is terribly nontrivial.
    
    # final mappings without simplifications
    sympy2dda = {
        sympy.Mul:   dda.mult,
        sympy.sqrt:  dda.sqrt,
        sympy.Abs:   dda.abs,
        sympy.Add:   lambda *x: dda.neg(dda.sum(*x)),
        sympy.exp:   dda.exp,
        sympy.floor: dda.floor,
        sympy.integral: lambda function, *symbols: dda.int(function)
    }
    
    dda_Symbol = Symbol # just to be verbose
    atom = sympy.Wild("x", properties=[lambda k: k.is_Symbol])
    
    # Do a depth-first-traversal mapping sympy to dda expressions
    def expr2dda(expr):
        if len(expr.args) == 0: # reached leaf
            if expr.is_number:
                return float(expr)
            else:
                raise ValueError(f"Found {expr} but don't know how to handle.\n(srepr: {srepr(expr)})")
        else: # map compound expression
            # First, try shorthands:
            
            # expressions 1/x for atom x:
            #( sympy.Mul(sympy.Integer(-1), atom), lambda match: dda.inv(
                
            # border cases which are really mad to catch:
            # srepr(-3*b) -> Mul(Integer(-1), Integer(3), Symbol('b'))
            
            # Last, try anything else
            if type(expr) in sympy2dda:
                return sympy2dda[type(expr)](*expr.args)
            else:
                raise ValueError(f"Found compound expression of type {type(expr)}, don't know how to handle. It is: {expr}\n(srepr: {srepr(expr)})")
            
        for arg in expr.args:
            expr2dda(arg)
    
    for eq in sympy_equation_list:
        eq = eq.canonical # ensure symbol on the left
        lhs, rhs = eq.args
        if not isinstance(lhs, sympy.Symbol):
            raise ValueError(f"Missing single symbol on LHS of Sympy equation {eq}")
        
        dda_lhs = dda_Symbol(lhs.name)
        
def to_sympy(state, symbol_mapper=identity, round_n=15):
    """
    Export a state to a set of equations for SymPy.
    Returns a list of ``sympy.Eq`` objects.
    Of course it requires Sympy installed/available.
    
    .. note::

       The heart of this function is a mapping from :class:`ast.Symbol`
       terms (functions) to Sympy functions, for instance by mapping
       ``Symbol("int")(...)`` to ``-sympy.Integral(sympy.Add(...), t)``.

       Thanks to the ease of the computing elements, this mapping does
       not require pattern matching but can be performed on a basic
       level. However, not all terms are yet supported.

    The argument `symbol_mapper` allows to apply another mapping
    on the DDA Symbol heads. By default, it is the identity
    function.
    
    With Sympy, you can do all funny things, such as:
    
    >>> from dda import *
    >>> x,int,neg=symbols("x,int,neg")
    >>> state = State({'x': int(neg(x), 0.2, 1)})
    >>> to_sympy(state)
    [Eq(x, -Integral(1.2 - x, t))]
    """
    
    sympy_Symbol = lambda label: sympy.Symbol(symbol_mapper(label))
    
    # round large floats. Doesn't work yet.
    rhs_rounder = lambda rhs: sympy.N(rhs, n=round_n)
    
    int = builtins.int # just to go sure
    t = sympy_Symbol("t")
    
    def todo(*x):
        raise ValueError("Requried DDA 2 Sympy function yet implemented")
    
    dda2sympy = {
        "const": lambda x: x,
        "neg":   lambda x: -x,
        "div":   lambda x,y: x/y,
        "int":   lambda *x: - sympy.Integral(sympy.Add(*x), t),
        "sum":   lambda *x: - sympy.Add(*x),
        "mult":  sympy.Mul,
        "sqrt":  sympy.sqrt,
        "abs":   sympy.Abs,
        "exp":   sympy.exp,
        "floor": sympy.floor,
    }
    
    def symbol2sympy(smbl):
        if isinstance(smbl, float) or isinstance(smbl, int):
            return smbl # let Sympy handle the numbers
        if not isinstance(smbl, Symbol): raise TypeError(f"Expecting symbol, got {smbl}")
        if smbl.is_variable():
            return sympy_Symbol(smbl.head) # just a variable
        else: # symbl.is_term()
            if smbl.head in dda2sympy:
                return dda2sympy[smbl.head](*map(symbol2sympy, smbl.tail))
            else:
                raise ValueError(f"DDA Symbol {smbl.head} in expression {smbl} not (yet) implemented.")
    
    equation_list = [ sympy.Eq(sympy_Symbol(lhs), rhs_rounder(symbol2sympy(state[lhs])))  for lhs in sorted(state) ]
    return equation_list

def to_latex(state, chunk_n=None):
    """
    Export to latex, using sympy.
    
    This mostly differs from ``sympy.latex`` for large equation
    systems where we use the align latex environment instead of
    a single equation. For the above example:
    
    >>> import sympy, dda
    >>> x,int,neg=dda.symbols("x,int,neg")
    >>> state = dda.State({'x': int(neg(x), 0.2, 1)})
    >>> print(sympy.latex(to_sympy(state)))  
    \\left[ x = - \\int \\left(1.2 - x\\right)\\, dt\\right]
    >>> print(to_latex(state))               
    \\begin{align}
    x &= - \\int \\left(1.2 - x\\right)\\, dt
    \\end{align}
    """
    
    equation_list = to_sympy(state)

    # Pretty-print set of equations. That's nicer then just
    # sympy.latex(self.equation_list)
    s = [ f"{sympy.latex(eq.lhs)} &= {sympy.latex(eq.rhs)}" for eq in equation_list ]
    equation_list_to_align = lambda t: r"\begin{align}"+"\n" + (r" \\"+"\n").join(t) + "\n" + r"\end{align}"
    
    def chunks(lst, n):
        """Yield successive n-sized chunks from lst."""
        for i in range(0, len(lst), n):
            yield lst[i:i + n]
    
    latex = "\n".join(map(equation_list_to_align, chunks(s,chunk_n) if chunk_n else [s]))
    return latex
 
#
#
#
# 
#

"""
C++ code generation is a major feature of the PyDDA code. The generated
code is built from a string template and has the following features:

* Standalone code: No further dependencies (beyond standard libc, lm and STL).
* Lightweight object oriented: Uses classes (structures) to hold the
  different variables (basically AoS instead of SoA approach). Little
  C++ templating.
* Organization in a few functions which allows to edit the generated 
  C++ code manually afterwards without going mad.
* CSV or binary output, or no output at all. Output is always made
  to stdout. Information messages are always sent to stderr.
* Debugging facilities built right into the code for setting NaNs and
  abortion in case of floating point exceptions.
* Runtime arguments via the commandline (argv): Parsing and passing.

.. note::

   **C++17** is required for building the C++ code. This is because
   we use variadic templates.

For the C++ runtime arguments, we support so far:

- *Simulation steering*: Selection of number of integration iterations
  and frequency of dumping the solution.
- *Query based plotting*: Selection which variables shall be outputted
  at runtime.
- Further *Flags* and *Numeric arguments* as well as a useful ``--help``
  message.
- *Initial data* and *time step sizes* can also be chosen at run time.
- *Introspection capabilities*, for instance one can ask the binary
  about the evolution quantities built in.

Basically the equation structure is the only thing left hardcoded at
C++ code generation time.

"""

from . import clean
from .ast import State, Symbol, topological_sort
from .computing_elements import cpp_impl

# could also say is_number = lambda obj: str(obj).is_digit()
is_number = lambda obj: isinstance(obj, float) or isinstance(obj, int)

import textwrap, itertools, math, functools, os

cpp_template = """\
// 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_name)s = std::numeric_limits<double>::signaling_NaN();

namespace dda {
%(cpp_impl)s
} // end namespace dda

using namespace dda;

// Time-evolved variables, actual "state" (in general dq/dt!=0)
struct %(state_type)s {
    %(state_var_definition)s
    %(dqdt_operators)s
    %(diff_or_integrate)s
    const double* byName(std::string %(lookup_name)s) const {
        %(state_vars_by_name)s
        return nullptr;
    }
};

// Auxiliary variables, derived from %(state_type)s, not evolved in time (dqdt=0)
struct %(aux_type)s  {
    %(aux_var_definition)s
    void set_to_nan() {
        %(aux_var_set_to_nan)s
    }
        %(aux_vars_by_name)s
    }
};

// Explicit constants
// Originally, in this code, they just have been constexprs.
// For some applications, it is useful to have them being runtime arguments. This has to
// be chosen at code generation time, they are then compile time constants.
#define %(const_type)s_is_really_const %(constexpr_consts)d
static %(const_type_qualifier)s struct %(const_type)s {
    %(explicit_constant_definitions)s
    %(const_type_qualifier)s double* byName(std::string %(lookup_name)s) %(const_type_qualifier)s {
        %(explicit_constants_by_name)s
    }
} %(const_name)s;

%(explicit_constant_initializations)s

// All variables, time evolveds, constants, as well-ordered strings:
%(all_variables_as_string)s

/// Compute the equations as given in the dda file
/// It is actually %(dqdt_name)s = f(%(state_name)s, %(aux_name)s), or at least
/// (%(dqdt_name)s,%(aux_name)s) = f(%(state_name)s).
/// %(aux_name)s is only returned for debugging and controlling purposes of intermediate results

void f(%(state_type)s const &%(state_name)s, %(state_type)s &%(dqdt_name)s, %(aux_type)s &%(aux_name)s) {
    if(debug) %(aux_name)s.set_to_nan(); // only for debugging: ensure no use of uninitialized variables

    %(equations)s
}

/// Helper routine: Compute aux variables from a state. Helpful for consistent printing.
void compute_aux(%(state_type)s const& %(state_name)s, %(aux_type)s &%(aux_name)s) {
    %(state_type)s dummy_%(dqdt_name)s;
    f(%(state_name)s, dummy_%(dqdt_name)s, %(aux_name)s);
}

%(state_type)s
    initial_data{ %(initial_data)s },
    dt{ %(timestep_data)s };
    
// not used:
// bool
//    shall_differentiate{ %(differentiate_data)s };
    
%(state_type)s integrate(%(state_type)s& %(state_name)s, %(aux_type)s& %(aux_name)s, int rk_order) {
    %(state_type)s k1, k2, k3, k4;
    %(state_type)s updated;

    switch(rk_order) {
        case 1:
            // Explicit Euler scheme
            f(%(state_name)s, k1, %(aux_name)s);
            updated = %(state_name)s + k1*dt;
            break;
        case 2:
            // RK2 scheme
            f(%(state_name)s + k1*dt, k2, %(aux_name)s);
            updated = %(state_name)s + (k1+k2)*dt*0.5;
            break;
        case 3:
            // Kutta's third order scheme 
            f(%(state_name)s + dt*k1*0.5, k2, %(aux_name)s);
            f(%(state_name)s + dt*k1*(-1.0) + dt*k2*2.0, k3, %(aux_name)s);
            updated = %(state_name)s + (k1 + k2*4.0 + k3*1.0)*dt*(1./6.);
            break;
        case 4:
            // Classical RK4 scheme 
            f(%(state_name)s + dt*k2*0.5, k3, %(aux_name)s);
            f(%(state_name)s + dt*k3*1.0, k4, %(aux_name)s);
            updated = %(state_name)s + (k1 + k2*2.0 + k3*2.0 + k4)*dt*(1./6.);
            break;
        default:
            exit(-42);
    }
    
    return updated;
}

struct csv_writer {
    std::vector<std::string> query_variables;
    bool write_initial_conditions, always_compute_aux_before_printing, binary_output, skip_header;

    const char* sep(size_t i) const { return (i!=query_variables.size() ? "\\t" : "\\n"); }
    
    void write_header() const {
        if(binary_output || skip_header) return;
        for(size_t i=0; i<query_variables.size(); i++)
            std::cout << query_variables[i] << sep(i+1);
    }
    
    void write_line(const %(state_type)s& %(state_name)s, const %(aux_type)s& %(aux_name)s, %(const_type_qualifier)s %(const_type)s& %(const_name)s) const {
        %(aux_type)s recomputed__%(aux_name)s;
        const %(aux_type)s *actual__%(aux_name)s = &%(aux_name)s;
        if(always_compute_aux_before_printing) {
            compute_aux(%(state_name)s, recomputed__%(aux_name)s);
            actual__%(aux_name)s = &recomputed__%(aux_name)s;
        }
        for(size_t i=0; i<query_variables.size(); i++) {
            const double* lookup=nullptr; std::string var = query_variables[i];
            if(!lookup) lookup = %(state_name)s.byName(var);
            if(!lookup) lookup = actual__%(aux_name)s->byName(var);
            if(!lookup) lookup = %(const_name)s.byName(var);
            if(!lookup) { std::cerr << "Lookup failed" << std::endl; exit(-123); }
            if(binary_output) {
                std::cout.write(reinterpret_cast<const char*>( lookup ), sizeof(double));
            } else {
                std::cout << *lookup << sep(i+1);
            }
        }
    }
} writer;

%(state_type)s simulate_dda(%(state_type)s initial, int max_iterations, int modulo_writer, int rk_order, int modulo_progress) {
    %(state_type)s %(state_name)s = initial, old = initial;
    %(aux_type)s %(aux_name)s;
        
    if(debug) {
#ifdef _GNU_SOURCE
        feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); // glibc-specific
#else
        // Debugging for non-glibc (ie. non-Linux), see https://stackoverflow.com/a/12762999
        // If you want do debugging on some other platform, please implement it on your own...
        std::cerr << "Warning: Debugging floating point not yet supported" << std::endl;
#endif
    }

    writer.write_header();

    if(writer.write_initial_conditions) {
        %(aux_name)s.set_to_nan(); // definite state for auxers at printing: not yet computed.
        writer.write_line(%(state_name)s, %(aux_name)s, %(const_name)s);
    }
    
    for(int iter = 0; iter < max_iterations; iter++) {
        %(state_type)s integrated = integrate(%(state_name)s, %(aux_name)s, rk_order);
        
        %(aux_type)s dummy; %(state_type)s dqdt;
        f(%(state_name)s, dqdt, dummy); // etxra round for dqdt only neccessary for first timestep
        %(state_type)s differences = (dqdt + old*(-1.0)) / dt;
        
        old = dqdt;
        %(state_name)s = %(state_type)s::diff_or_integrate(integrated, differences);
            
        // TODO: Currently, differentiation is always first order forward
        //       in time. Could also do higher order.

        if(iter %% modulo_writer == 0)
            writer.write_line(%(state_name)s, %(aux_name)s, %(const_name)s);
            
        if(modulo_progress >= 0 && iter %% modulo_progress == 0) {
            fprintf(stderr, "%%.0d%%%%... ", 100*iter / max_iterations );
            fflush(stderr);
        }
    }
    
    return %(state_name)s;
}

using namespace std;

vector<string> split(const string& input, const char seperator='=') {
    stringstream s(input);
    string segment;
    vector<string> seglist;
    while(getline(s, segment, seperator))
        seglist.push_back(segment);
    return seglist;
}

template<typename T>  T extract(const string& input) {
    T target;
    istringstream ss(input); ss >> target;
    if(ss.fail()) {
        cerr << "Failure: Not a number: " << input << endl;
        exit(-4);
    }
    return target;
}

bool unprefix(string& str, const string prefix) {
    bool is_prefixed = !str.compare(0, prefix.size(), prefix);
    if(is_prefixed) str = str.substr(prefix.size());
    return is_prefixed;
}

bool contains(const std::vector<std::string>& haystack, const std::string& needle) {
    return std::find(haystack.begin(), haystack.end(), needle) != haystack.end();
}

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main(int argc, char** argv) {
    // cstd rand 
    srand((unsigned int)time(NULL));

    // runtime-definable arguments and their values:
    
    map<string, int> numbers;
    numbers["modulo_write"] = 1;
    numbers["max_iterations"] = 100;
    numbers["rk_order"] = 1;
    numbers["number_precision"] = 5;
    numbers["modulo_progress"] = -1;
    
    map<string, bool> flags;
    flags["debug"] = false;
    flags["list_all_variables"] = false;
    flags["write_initial_conditions"] = false;
    flags["always_compute_aux_before_printing"] = true;
    flags["binary_output"] = false;
    flags["skip_header"] = false;

    // Our primitive argument processing:
    vector<string> args(argv + 1, argv + argc);
    for(auto arg : args) {
        if(arg == "--help") {
            string ind = "  ";
            cerr << "Usage: " << argv[0] << " [arguments] <variables_to_print>" << endl;
            cerr << "This is an ODE integrator generated by PyDDA." << endl << endl;
            cerr << "* Boolean arguments (Usage --foo or --foo=1 or --foo=0)" << endl;
            for(auto const& [key, val] : flags) cerr << ind << key << " (default value: " << val << ")" << endl;
            cerr << "* Numeric arguments: (Usage --foo=123)" << endl;
            for(auto const& [key, val] : numbers) cerr << ind << key << " (default value: " << val << ")" << endl;
            cerr << "* Overwrite initial conditions (initial data): (Usage --initial:foo=1.23)" << endl;
            for(auto const& key : all_%(state_type)s) cerr << ind << key << " (default value: " << *initial_data.byName(key) << ")" << endl;
            cerr << "* Overwrite (per-variable) time step size: (Usage --dt:foo=0.001 or even --dt:foo=1e-6)" << endl;
            for(auto const& key : all_%(state_type)s) cerr << ind << key << " (default value: " << *dt.byName(key) << ")" << endl;
            #if %(const_type)s_is_really_const == 0
            cerr << "* Overwrite (runtime) constants: (Usage --const:foo=1.23)" << endl;
            for(auto const& key : all_%(const_type)s) cerr << ind << key << " (default value: " << *%(const_name)s.byName(key) << ")" << endl;
            #else
            cerr << "* Overwrite constants: Not possible since compiled as compile time constants." << endl;
            #endif
            cerr << "* Query fields: (if none given, all are dumped)";
            for(size_t i=0; i<all_variables.size(); i++) { cerr << endl << ind;
                for(size_t j=0;j<5 && i<all_variables.size();j++) cerr << all_variables[i] << (i+1!=all_variables.size() ? ", " : ""); }
            cerr << endl << endl << "Exemplaric usage:" << endl;
            cerr << ind << argv[0] << " --foo=1 --bar=0 --baz=7 --ic:var1=0.5 --dt:var2=0.01 --const:something=42 var1 var2 var3" << endl;
            cerr << endl << "For more options and help, see the PyDDA code documentation at https://pyanalog.readthedocs.io/" << endl;
            exit(-1);
        }
        
        // TODO: Compress the following repetitive lines with some function.
        if(unprefix(arg, "--dt:")) {
            vector<string> splitted = split(arg);
            if(splitted.size() != 2) { cerr << "ERR: Usage: --dt:key=value"; exit(-2); }
            string key=splitted[0]; double value=extract<double>(splitted[1]);
            if(!dt.byName(key)) { cerr << "ERR: --dt:"<<key<<" is not a valid argument. Check with --help"; exit(-2); }
            *(const_cast<double*>(dt.byName(key))) = value;
        } else if(unprefix(arg, "--initial:")) {
            if(splitted.size() != 2) { cerr << "ERR: Usage: --initial:key=value"; exit(-2); }
            if(!initial_data.byName(key)) { cerr << "ERR: --initial:"<<key<<" is not a valid argument. Check with --help"; exit(-2); }
            *(const_cast<double*>(initial_data.byName(key))) = value;
        } else if(unprefix(arg, "--const:")) {
            if(splitted.size() != 2) { cerr << "ERR: Usage: --const:key=value"; exit(-2); }
            if(!%(const_name)s.byName(key)) { cerr << "ERR: --const:"<<key<<" is not a valid argument. Check with --help"; exit(-2); }
                *(%(const_name)s.byName(key)) = value;
            #else
                cerr << "ERR: Cannot overwrite compile time constant." << endl;
                cerr << "ERR: Regenerate code with constexpr_consts=True in order to use this feature." << endl;
                exit(-3);
            #endif
        } else if(unprefix(arg, "--")) {
            if(splitted.size() == 1) splitted.push_back("1"); // for the booleans
            if(splitted.size() != 2) {
                cerr << "ERR: Illegal key-value argument " << arg << endl;
                cerr << "ERR: Try --help" << endl;
                exit(-3);
            }
            string key=splitted[0], value=splitted[1];
            if(numbers.count(key))    numbers[key] = extract<int>(value);
            else if(flags.count(key)) flags[key] = extract<bool>(value);
            else {
                cerr << "ERR: Illegal argument: " << arg << endl;
                exit(-3);
            }
        } else {
            if(contains(all_variables, arg)) {
                writer.query_variables.push_back(arg);
            } else {
                cerr << "ERR: It is not a member of the variables. Try --help." << endl;
                exit(-2);
            }
        }
    } // for args
    
    if(flags["list_all_variables"]) {
        for(auto var : all_variables) puts(var.c_str());
        exit(0);
    }
    
    // cout.set_precision(numbers["number_precision"]); // TODO, move above correctly, before every number
    
    if(writer.query_variables.empty())
        writer.query_variables = all_variables;

    debug = flags["debug"];
    writer.write_initial_conditions = flags["write_initial_conditions"];
    writer.always_compute_aux_before_printing = flags["always_compute_aux_before_printing"];
    writer.binary_output = flags["binary_output"];
    writer.skip_header = flags["skip_header"];
    
    simulate_dda(initial_data, numbers["max_iterations"], numbers["modulo_write"], numbers["rk_order"], numbers["modulo_progress"]);
}

"""

def to_cpp(state, number_precision=math.inf, constexpr_consts=False):
    """
    Given a state, returns standalone C++ code as string.
    
    This code can be written to a file, compiled with a recent 
    C++ compiler and then solves the differential equation system when executed.
    
    The algorithm is basically:
    
    1. linearize the state (this can raise)
    2. determine all the C++ template fields
    3. Return the filled out template
    
    We plan to add logging for non-fatal information about the C++ code quality
    (see TODOs in the code).
    
    The argument ``number_precision`` currently has no effect.
    """

    indent = " "*5 # or tab, whatever you prefer - should be fit to cpp_template

    # Despite all user-chosen variable names are scoped within structs/classes, name
    # clashes are possible in some contexts. Therefore, the following names should be
    # chosen carefully.
    state_type, aux_type, const_type = "state_variables", "auxillaries", "constants"
    state_name, dqdt_name, aux_name, other_name, const_name = "_state", "_dqdt", "_aux", "_other", "_constants"
    lookup_name, nan_name = "_name_", "_nan_"

    state = clean(state, target="C").name_computing_elements()
    vars = state.variable_ordering()
    
    # NB: One could require that there no implicit constants, i.e. numbers appearing randomly
    #     on the RHS of expressions. One could even go further and allow numbers only in
    #     potentiometer settings
    
    if not vars.evolved:
        pass # TODO: Warn that there are (interesting) variables at all.
    if vars.aux.cyclic:
        pass # TODO: Warn that there are unwanted feedback loops

    # prepare for prefixing all RHS variables
    pystruct2cprefix = { state_name: vars.evolved, aux_name: vars.aux.all, const_name: vars.explicit_constants  }
    struct_lookup = { var: cstruct for (cstruct, pystruct) in pystruct2cprefix.items()  for var in pystruct  }
    csymbol = lambda var: f"{struct_lookup[var]}.{var}"
    prefix_rhs = lambda el: Symbol(csymbol(el.head)) if el.is_variable() else el
    
    def number2cstr(number):
        if not is_number(number):
            raise ValueError(f"Was expecting a number, but got {type(number)} with value {number}.")
        formatstring = f"%{number_precision}f" if number_precision < math.inf else "%f"
        return formatstring % number
    
    # Remove any const() which remained. Would be nicer to assert not having consts() at all.
    remove_const = lambda x: number2cstr(x.tail[0]) if isinstance(x,Symbol) and x.head=="const" else x
    
    # Translate const(foo) by stripping foo or some bar by looking up if it is an explicit constant.
    # Dynamical variables are not allowed here.
    def lookup_const(var):
        if isinstance(var,Symbol):
            if var.head == "const":
                var = var.tail[0] # continue
            elif var.is_variable():
                if not var.head in vars.explicit_constants:
                    raise ValueError(f"Only constants allowed in this context. {var} however refers to {csymbol(var.head)}.")
                return csymbol(var.head)
            else: # remaining case: var.is_term()
                raise ValueError(f"Was expecting const(foo) or so, but got term {var}.")
        if not is_number(var): raise ValueError(f"Got a weird {type(var)} in a constant context: {var}")
        return number2cstr(var)

    # HOTFIX rename reserved keywords in the C language
    # FIXME: Should have a better solution valid for *every* possible collision. 
    # C++ namespaces are apparently not enough.
    c_names = { "div": "Div", "abs": "Abs" }
    c_substitute = lambda head: c_names.get(head, head)
    state = State({ var: state[var].map_heads(c_substitute) for var in state })

    # Extract int(..., timestep, initial_data) and rewrite reserved C keyword
    # Also extract diff(..., timestep) and rewrite symbol
    timesteps = {}
    initial_data = {}
    differentiate = {}
    def map_and_treat_integrals(var):
        if not var in vars.evolved: return state[var]
        assert state[var].head in ["int", "diff"], f"Expected int(...) or diff(...) but got {state}"
        head, tail = state[var].head, state[var].tail
        if not len(tail) >= 3: raise ValueError(f"{head}(...) requires at least {head}(value, dt, initial_value)")
        timesteps[var] = lookup_const(tail[-2])
        initial_data[var] = lookup_const(tail[-1])
        differentiate[var] = head == "diff"
        return Symbol(head.capitalize(), *tail[0:len(tail)-2])
    state = State({ var: map_and_treat_integrals(var) for var in state })

    # do the renaming *after* variable dependency analysis
    state = state.map_tails(remove_const, map_root=False)
    state = state.map_tails(prefix_rhs,   map_root=False)
        
    # C-format lists of statements or so. Do indentation.
    J = lambda whatever: ", ".join(whatever)
    C = lambda whatever: textwrap.indent(whatever if isinstance(whatever, str) else  "\n".join(whatever), indent)
    CC = lambda whatevr: C(C(whatevr)) # two indentations ;-)
    Cw = lambda whatever: C(textwrap.wrap(whatever, width=50))
    enc = lambda ws, inner: ws[0] + inner + (ws[1] if len(ws)==2 else ws[0]) # enc('"','foo') or enc('{}','foo')
    varlist = lambda ctype, lst: Cw(f"{ctype} {', '.join(lst)};") if lst else ""

    state_var_definition = varlist("double", vars.evolved)
    aux_var_definition = varlist("double", vars.aux.all)
    
    # For runtime (i.e. CLI invocation) introspection capabilities
    lookup = lambda struct, prefix="": CC(f'if({lookup_name} == "{var}") return &{prefix}{var};' for var in struct)
    state_vars_by_name = lookup(vars.evolved)
    aux_vars_by_name = lookup(vars.aux.all)
    explicit_constants_by_name = lookup(vars.explicit_constants, prefix=const_type+"::")
    
    # For even more runtime (i.e. CLI invocation) introspection capabilities
    variable_namings = {
        "variables": vars.all,
        state_type:  vars.evolved,
        aux_type:    vars.aux.all,
        const_type:  vars.explicit_constants,
    }
    all_variables_as_string = varlist("static const std::vector<std::string>",
        ["all_"+k+" = "+enc("{}", Cw(', '.join([ enc('"',vi) for vi in v ])))
            for k,v in variable_namings.items() ]).lstrip()
    
    # the above replaces repepetitive stuff such as:
    #all_variables_as_string = C(f'"{var}",' for var in vars.all)
    #all_state_vars_by_name = C(f'"{var}",' for var in vars.evolved)
    #all_constants_by_name = C(f'"{var}",' for var in vars.explicit_constants)

    # For debugging:
    aux_var_set_to_nan = C(f"{var} = {nan_name};" for var in vars.aux.all)
    
    # Workaround for avoiding logic in template...
    either = lambda lst, otherwise: lst if lst else [otherwise]
    nonemptyvarlist = lambda lst: either(lst, otherwise="_none=-1")
    
    const_assignments = [f"{var}={state[var].tail[0]}" for var in vars.explicit_constants] # extract x in const(x)
    if constexpr_consts:
        const_type_qualifier = "const"
        # generates a member in const struct foo { constexpr static double bla = 1, blo = 2, ...; }
        explicit_constant_definitions = varlist("constexpr static double", nonemptyvarlist(const_assignments))
        explicit_constant_initializations = ""
    else:
        const_type_qualifier = ""
        # generates a member in struct foo { static double bla, blo, ...; }
        explicit_constant_definitions = varlist("static double", either(vars.explicit_constants, otherwise="_none"))
        # and subsequently foo::bla=1, foo::blo=2, ...
        explicit_constant_initializations = C(f"double {const_type}::{ass};" for ass in const_assignments)
    
    initial_data = J(str(initial_data[var]) for var in vars.evolved)
    timestep_data = J(str(timesteps[var]) for var in vars.evolved)
    
    diff_or_integrate = C([ \
        f"static {state_type} diff_or_integrate(const {state_type} &integrals, const {state_type} &difference) "+'{', \
        C(f"{state_type} {state_name};"), \
        C(f"{state_name}.{var} = {'difference' if differentiate[var] else 'integrals'}.{var};" for var in vars.evolved), \
        C(f"return {state_name};"), \
        "}"])
    
    # leftover/bookkeeping/better doing    
    differentiate_data = J(str(differentiate[var]) for var in vars.evolved)

    #state_assignments = lambda lst: C(f"{v} = {state[v]};" for v in lst)) if lst else C("/* none */")
    state_assignments = lambda lhs_struct,lst: [f"{lhs_struct}.{var} = {state[var]};" for var in lst] if lst else ["/* none */"]

    # The dependency is basically
    #   aux = function_of(aux, state)
    #   dqdt = function_of(aux, state)
    # and in the integration schema step
    #   state = function_of(dqdt)
    
    equations = []
    equations.append("// 1. Topologically sorted aux variables")
    equations += state_assignments(aux_name, vars.aux.sorted)
    equations.append("// 2. Cyclic aux variables")
    equations += state_assignments(aux_name, vars.aux.cyclic)
    # FIXME:
    # There is a special case like the following where unneeded auxers are wrongly
    # detected as unneeded, depsite urgently needed:
    """
    # CONTEXT: dipole_water.py
     // 3. State variable changes (dqdt), finally
     _dqdt.mdm0_0 = Int(_aux.Fm0_01, _aux.Fn0_00, _aux.Fn0_01);
     _dqdt.mdm0_1 = Int(_aux.Fm0_01, _aux.Fn0_10, _aux.Fn0_11);
     _dqdt.mdm1_0 = Int(_aux.Fm1_01, _aux.Fn1_00, _aux.Fn1_01);
     _dqdt.mdm1_1 = Int(_aux.Fm1_01, _aux.Fn1_10, _aux.Fn1_11);
     _dqdt.mdr0_0 = Int(_aux.Fr0_01, _aux.Fn0_00, _aux.Fn0_01);
     _dqdt.mdr0_1 = Int(_aux.Fr0_01, _aux.Fn0_10, _aux.Fn0_11);
     _dqdt.mdr1_0 = Int(_aux.Fr1_01, _aux.Fn1_00, _aux.Fn1_01);
     _dqdt.mdr1_1 = Int(_aux.Fr1_01, _aux.Fn1_10, _aux.Fn1_11);
     _dqdt.mu0_0 = Int(_aux.mult_19);
     _dqdt.mu0_1 = Int(_aux.mult_37);
     _dqdt.mu1_0 = Int(_aux.mult_45);
     _dqdt.mu1_1 = Int(_aux.mult_51);
     _dqdt.r0_0 = Int(_aux.mult_20);
     _dqdt.r0_1 = Int(_aux.mult_38);
     _dqdt.r1_0 = Int(_aux.mult_46);
     _dqdt.r1_1 = Int(_aux.mult_52);
     // 4. Unneeded auxilliary variables (maybe postprocessing, etc.)
     _aux.mult_19 = mult(0.5, _state.mdm0_0);
     _aux.mult_46 = mult(0.0625, _state.mdr1_0);
     _aux.mult_52 = mult(0.0625, _state.mdr1_1);
     _aux.mult_51 = mult(0.5, _state.mdm1_1);
     _aux.mult_37 = mult(0.5, _state.mdm0_1);
     _aux.mult_45 = mult(0.5, _state.mdm1_0);
     _aux.mult_20 = mult(0.0625, _state.mdr0_0);
     _aux.mult_38 = mult(0.0625, _state.mdr0_1);
     """
    equations.append("// 3. Unneeded auxilliary variables (maybe postprocessing, etc.)")
    equations += state_assignments(aux_name, vars.aux.unneeded)
    equations.append("// 4. State variable changes (dqdt), finally")
    equations += state_assignments(dqdt_name, vars.evolved)
    equations = C(equations)

    make_operator = lambda operator_symbol, other_type, a=True: [ \
        f"{state_type} operator{operator_symbol}(const {other_type} &{other_name}) const "+'{', \
        C(f"{state_name}.{var} = {var} {operator_symbol} {other_name}{'.'+var if a else ''};" for var in vars.evolved), \
        "}"]

    dqdt_operators = C(C(make_operator(s,o,a)) for s,(o,a) in
        itertools.product("*+/", zip((state_type, "double"), (True,False,False))))

    output = cpp_template % {**locals(), **globals()}
    return output

# What follows are a few helper functions to make the usage nicer

default_c_filename = "generated.cc"
default_compiler_output = "./a.out"

def compile(code, c_filename=default_c_filename, compiler="g++", compiler_output=default_compiler_output, options="--std=c++17 -Wall"):
    """
    Small helper function to compile C++ code from python.
    
    Write string *code* to *c_filename* and run the *compiler* on that, afterwards.
    Will raise an error if compilation fails.
    """
    from os import system
    with open(c_filename, "w", encoding="utf-8") as fh:
        print(code, file=fh)
    if system(f"{compiler} -o{compiler_output} {options} {c_filename}"):
        raise ValueError("Could not compile C source!")

def runproc(command, decode=False):
    "Helper to run external command and slurp its output to a binary array"
    import subprocess # internals
    # Popen: Python 3.6 does not yet have the text keyword. Newer Pythons have.
    # We do a decode("utf-8") later for the same purpose.
    proc = subprocess.Popen(command, stdout=subprocess.PIPE)#, stderr=subprocess.PIPE) # text=True)        
    # pass stderr to user for progress bar.
    (stdout,stderr) = proc.communicate()
    #stderr += proc.stderr.read()
    if decode:
        stdout = stdout.decode("utf-8")
        #stderr = stderr.decode("utf-8")
    return proc, stdout#, stderr

def list_all_variables(command=default_compiler_output):
    "in order to know which fields have been read, slurp all variables"
    fproc, fields = runproc([command, "--list_all_variables"], decode=True)
    if fproc.returncode:
        #warnings.warn("Could not extract evolution variables automagically from command: "+stderr+", Stdout: "+stdout)
        #warnings.warn("Could not extract evolution variables automagically from command (see stderr above). Stdout: "+stdout)
        raise ValueError("Could not extract evolution variables automagically from command (see stderr above). Stdout: "+stdout)
    else:
        fields_to_export = fields.strip().split("\n")
    return fields_to_export

def run(command=default_compiler_output, binary=False, arguments={}, fields_to_export=[]):
    """
    Small helper function to execute a code generated by this module.
    
    Runs ``command`` on the command line, with given dict ``arguments`` in ``--foo=bar``
    fashion and ``fields_to_export`` just as a sequential argument list. If no
    ``fields_to_export`` is given, ``command --list_all_variables`` will be run to query
    all default fields.
    
    Pipes stdout to a string, which is returned. Stderr will just be passed.
    The function will return once the binary finished or raise in case of error.
    
    If you set ``binary=True``, raw data instead of CSV files will be passed between the spawned
    command and this python program. This decreases the runtime significantly if you write a lot
    of data (since CSV generating and parsing overload is gone). 
    
    Example usage:
    
    >>> state = State()
    >>> state["x"] = Symbol("int", Symbol("neg", state["x"]), 0.2, 1)
    >>> state
    State({'x': int(neg(x), 0.2, 1)})
    >>> cpp_code = to_cpp(state)
    >>> print(cpp_code)  # doctest: +SKIP
    // 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>
    ....
    >>> compile(cpp_code, compiler_output="foo.exe")
    >>> res = run("./foo.exe", arguments={'max_iterations':10}, fields_to_export=['x']) # doctest: +SKIP
    Running: ./foo.exe --max_iterations=10 x
    TODO: Doctesting this doesn't work good due to stderr (cf https://stackoverflow.com/a/61533524)
    >>> print(res) # doctest: +NORMALIZE_WHITESPACE +SKIP
    x
    1.2
    1.44
    1.728
    2.0736
    2.48832
    2.98598
    3.58318
    4.29982
    5.15978
    6.19174

    """
    if binary: arguments["binary_output"] = 1
    command_with_args = [command] + [ f"--{k}={v}" for k,v in arguments.items() ] + fields_to_export
    print(f"Running: {' '.join(command_with_args)}")
    proc, stdout = runproc(command_with_args, decode=not binary)

    if proc.returncode:
        print(stdout)
        #print("STDERR:", stderr)
        raise ValueError(f"Could not execute '{command_with_args}'. Please run on the command line for inspection. Probably use gdb.")
    else:
        return stdout

def numpy_read(stdout, binary=False, return_ndarray=True, return_recarray=False, fields_to_export=[]):
    """
    Postprocessing to fill the gap between the C++ output and a suitable numpy array.
    In order to so, this function has to know whether your output was binary or text.
    Furthermore, you need to tell him how many fields you had. You can use :meth:`list_all_variables`
    
    Old Text:
    
    This option only makes real sense if
    you set (the default) ``return_ndarray=True``. Note that if you don't pass the  ``fields_to_export``
    option but set ``binary=True``, in the moment the returned array is one-dimensional (a warning
    will be printed). If you like even more structured data be returned, turn on ``return_recarray=True``.
    It will return a numpy.recarray, the same data type which you get when you read CSV data
    with column headers. ``return_recarray=True`` implies ``return_ndarray=True``.
    """
    if return_recarray:
        return_ndarray = True
    if return_ndarray:
        import warnings # bultin
        import numpy as np # external
        
        if binary:
            # note that `binary` will give you a long array of numbers instead of a two-dimensional array.
            # Any structural information is lacking and has to be introduced again.
            
            if not fields_to_export:
                raise ValueError("No fields given. Please pass fields_to_export=list_all_variables().")
            
            data = np.frombuffer(stdout, dtype="<f8")
            if fields_to_export:
                N = len(data) / len(fields_to_export)
                assert N.is_integer(), "1d number sequence from C run is not divisible by number of fields."
                data = data.reshape(int(N), len(fields_to_export))
                if return_recarray:
                    dt = [ (f, data.dtype) for f in fields_to_export ]
                    data = data.view(dt)  # this will be of shape (N,1)
                    data = data[:,0]      # strip this unneccessary extra axis!
            else:
                warnings.warn("Returning one-dimensional numpy array instead of a structured two-dimensional array. Pass fields_to_export=['foo','bar','baz'] in order to recieve two-dimensional array.")
                pass
        else:
            from io import StringIO
            if return_recarray:
                data = np.genfromtxt(StringIO(stdout), names=True)
            else:
                data = np.genfromtxt(StringIO(stdout), skip_header=1)
        return data
    else:

class Solver:
    """
    Syntactic sugar for a more concise OOP feeling. Instead of calling ``export(to="C")``,
    ``compile()`` and ``run()`` you can just write ``Solver(state, runtime_arguments)``.
    This object will even clean up after running.
    """
    def __init__(self, dda_state_or_code, *runtime_fields_to_export, constexpr_consts=True, **runtime_arguments):
        self.c_code = to_cpp(dda_state_or_code, constexpr_consts=constexpr_consts) if isinstance(dda_state_or_code, State) else dda_state_or_code
        
        self.code_name = "cpp_generated.cc"
        self.output_name = "cpp_generated.exe"
        self.default_runtime_fields_to_export = runtime_fields_to_export
        self.default_runtime_arguments = runtime_arguments
        
        compile(self.c_code, c_filename=self.code_name, compiler_output=self.output_name)
        
    def run(self, *runtime_fields_to_export, binary=False, cleanup=True, **runtime_arguments):
        "Chaining and Syntactic sugar for delayed argument setting/overwriting"
        self.runtime_arguments = { **runtime_arguments, **self.default_runtime_arguments }
        self.runtime_fields_to_export = list(runtime_fields_to_export if runtime_fields_to_export else self.default_runtime_fields_to_export)
        self.binary = binary
    
        self.stdout = run(command="./" + self.output_name, binary=self.binary,
                  arguments=self.runtime_arguments,
                  fields_to_export=self.runtime_fields_to_export)
        
        if self.binary and not self.runtime_fields_to_export:
            # Could also extract field values from dda state...
            self.fields = list_all_variables("./" + self.output_name)
        else:
            self.fields = self.runtime_fields_to_export
        
        if cleanup:
            os.unlink(self.code_name)
            os.unlink(self.output_name)
        return self # chainable!

    def as_ndarray(self):
        "Return run results as a ``np.ndarray`` (i.e. like a table without headers, typically 2D data)"
        if not hasattr(self, "stdout"): raise LogicError("First call run() on the simulator before data readout.")
        return numpy_read(self.stdout, binary=self.binary, return_ndarray=True, return_recarray=False, fields_to_export=self.fields)
        
    def as_recarray(self):
        "Return run results as a ``np.recarray`` (i.e. like CSV table with named headers)"
        return numpy_read(self.stdout, binary=self.binary, return_ndarray=True, return_recarray=True,  fields_to_export=self.fields)

#
#
#
# 
#

"""
The DDA domain specific language, also refered to as "traditional dda",
is the C-like language invented by Bernd for his Perl'ish dda code.

It basically reads as the following snippet:

::

    # Single-line comments are written like this

    dt = const(0.5)  # constants are defined like this
    y0 = const(1)

    z = mult(y, y)
    y = int(y, dt, y0)

As you see, variables do not even have to be introduced and can be
used in any order. There is only one data type, the *analog line signal*
which is basically a real number within a fixed interval.

.. note::

   Interestingly, traditional DDA files are a python subset and thus can
   be easily parsed and generated by python syntax. That's why the code
   of this PyDDA module is so short. That's also a primary reason why
   DDA was rewritten in Python.
   
Demonstration usage of this module:

>>> traditional_dda_document='''
... # Single-line comments are written like this
... 
... dt = const(0.5)  # constants are defined like this
... y0 = const(1)
... 
... z = mult(y, y)
... y = int(y, dt, y0)
... '''
>>> state = read_traditional_dda(traditional_dda_document)
>>> state
State({'dt': const(0.5), 'y': int(y, dt, y0), 'y0': const(1), 'z': mult(y, y)})
>>> print(to_traditional_dda(state))   # doctest: +NORMALIZE_WHITESPACE
# Canonical DDA file generated by PyDDA
<BLANKLINE>
dt = const(0.5)
y = int(y, dt, y0)
y0 = const(1)
z = mult(y, y)

"""

from . import clean, ast as dda # in order to explicitely write dda.Symbol, dda.State
import collections, sys, argparse, builtins, os, inspect # python included

def to_traditional_dda(state, cleanup=True, prefix="# Canonical DDA file generated by PyDDA\n", suffix="\n"):
    """
    Export state to canonical dda file format (i.e. without all the python).
    
    Returns the generated DDA file as string
    """
    
    if cleanup:
        # Not sure about scattered consts, maybe just remove them for the time being.
        remove_const = lambda x: x.tail[0] if isinstance(x,dda.Symbol) and x.head=="const" else x
        state = clean(state, target="dda").map_tails(remove_const, map_root=False)
        # TODO: Treat constants better. They have deserved it!

    dda_lines = []
    if prefix: dda_lines.append(prefix)
    for k in sorted(state):
        dda_lines.append(f"{k} = {state[k]}")
    if suffix: dda_lines.append(suffix)

    output = "\n".join(dda_lines)

def read_traditional_dda(content, return_ordered_dict=False):
    """
    Read some traditional dda file. We use the Python parser (``ast`` builtin)
    for this job. This is possible because the DDA syntax is a python subset and
    the parser doesn't care about semantics, only syntax.
    
    Thanks to the ``ast`` builtin package, we can just transform the python AST to
    the Symbolic/State class data structures used in this module.
    
    .. note::
    
       If some of the assertions fail, you can debug your DDA file by inspecting
       the output of ast.parse(content) on iPython. You can also run the Python
       debugger (pdb) on this function, for instance in iPython:
    
       >>> %pdb                                                  # doctest: +SKIP
       >>> read_traditional_dda(file("foo.dda").read())          # doctest: +SKIP
    
    Returns a state instance or OrderedDict, on preference.
    """
    import ast # python builtin
    tree = ast.parse(content)
    
    assert type(tree) == ast.Module, "I was expecting a whole file as content"
    assert type(tree.body) == list, "DDA file malformed, I was expecting a list of statements"
    assert all(type(f) == ast.Assign for f in tree.body), "DDA file malformed, I was expecting a list of assignments only"
    
    def expr2str(ast_obj):
        "Get source back; https://stackoverflow.com/questions/32146363/python-ast-abstract-syntax-trees-get-back-source-string-of-subnode"
        if hasattr(ast, "get_source_segment"): # Python >=3.8
            return ast.get_source_segment(content, ast_obj)
        else:
            return ast.dump(ast_obj) + " (Hint: Use Python >=3.8 to get source segment)"

    def arg2symbol(argument):
        "Transform some DDA function argument to the Symbol hierarchy"
        if isinstance(argument, ast.Constant):
            return argument.value
        elif isinstance(argument, ast.Name):
            return dda.Symbol(argument.id)
        elif isinstance(argument, ast.Call):
            return call2symbol(argument)
        elif isinstance(argument, ast.UnaryOp):
            # something like -1 is represented as ~ ast.USub(-1)
            return ast.literal_eval(argument)
        elif isinstance(argument, ast.Num):
            return argument.n
        else:
            raise TypeError(f"Don't understand argument type '{expr2str(argument)}'")
    
    def call2symbol(statement):
        "Transform some Right Hand Side nested function call to Symbol hierarchy"
        assert type(statement) == ast.Call, f"Was expecting a simple f(x) call but got '{expr2str(statement)}'"
        assert len(statement.keywords) == 0, f"Did not expect pythonic keyword arguments f(x=bar) in '{expr2str(statement)}'"
        assert type(statement.func) == ast.Name, f"Dunno, what is {statement.func}?"
        head = statement.func.id
        tail = map(arg2symbol, statement.args)
        return dda.Symbol(head, *tail)
    
    def ast_assignment_to_tuple(assign):
        assert len(assign.targets)==1, f"Was expecting only a single assignment, but got '{expr2str(assign)}'"
        assert type(assign.value) == ast.Call, f"DDA file malformed, expecting foo=call(bar), but got '{expr2str(assign)}'."
        variable_name = assign.targets[0].id
        rhs = call2symbol(assign.value)
        return (variable_name, rhs)
    
    result = map(ast_assignment_to_tuple, tree.body)
    mapping = collections.OrderedDict(result)
    return mapping if return_ordered_dict else dda.State(mapping)

def read_traditional_dda_file(filename, **kwargs):
    """
    Syntactic sugar for :meth:`read_traditional_dda`, so users can directly pass a filename
    if the have their DDA code in a file.
    """
    with open(filename, "r") as fh:
        content = fh.read()
    return read_traditional_dda(content, **kwargs)

def cli_exporter():
    """
    A Command Line Interface (CLI) for PyDDA.

    This CLI API does mainly what the old dda2c.pl script did, i.e.
    translating a (traditional) DDA file to C code. There are fewer
    options, because --iterations, --modulus and --variables are
    now runtime options for the generated C program.
    
    However, we can generate much more then C. Output is always text.
    
    Invocation is either ``python -m dda --help`` or ``python -m dda.dsl --help``
    anywhere from the system. ``setup.py`` probably also installed a
    ``pydda`` binary somewhere calling the same. You can also just
    call ``./dsl.py --help``.
    """
    
    parser = argparse.ArgumentParser(description="PyDDA, the AST-based DDA compiler", epilog=inspect.getdoc(cli_exporter))

    parser.add_argument("circuit_file", nargs='?', type=argparse.FileType('r'), default=sys.stdin, help="DDA setup (traditional file). Default is stdin.")
    parser.add_argument("-o", "--output", nargs='?', type=argparse.FileType('w'), default=sys.stdout, help="Where to write exported code to. Default is stdout.")
    
    # case insensitive choices. Will be parsed by export() anyway
    parser.add_argument("format", choices=["c", "dda", "dot", "latex"], type=str.lower, help="File formats which can be generated")
    
    # The following is deprecated because it is now runtime information for the C++ code.
    #C = parser.add_argument_group(title="Arguments for C++ code generation (Only apply if --export=C)")
    #C.add_argument("-N", "--iterations", type=int, help="Number of integration steps to be performed")
    #C.add_argument("-m", "--modulus", type=int, help="Output a value every <modulus> iteration steps")
    #C.add_argument("-v", "--variables", nargs="*", help="List of variables to be plotted (comma seperated)")
    
    arg = parser.parse_args()
    
    dda_text = arg.circuit_file.read()
    state = read_traditional_dda(dda_text)
    exported_code = state.export(to=arg.format)
    print(exported_code, file=arg.output)
    
if __name__ == "__main__":
    cli_exporter()

#
#
#
# 
#

"""
PyDDA is a small library to write and generate DDA code in Python.
DDA stands for *digital differential analyzer*. In this context, it is
a code for solving ordinary differential equations given in a domain
specific language description (i.e. an electrical circuit).

For further details, please see the ``doc/`` directory (Sphinx
documentation).
""" 

import sys

assert sys.version_info >= (3, 6), \
   "PyDDA uses f-strings all the ways, which require Python 3.6."

def export(state, to, **kw):
    """
    Convenience function to export (transform) a state to some other
    programming language.

    Possible formats (allowed values for ``to``) supported so far are:
    
    * C/C++ (via :py:mod:`dda.cpp_exporter`)
    * DDA   (via :py:mod:`dda.dsl`)
    * SymPy (via :py:mod:`dda.sympy`)
    * Latex (via :py:mod:`dda.sympy`)

    This function shall be nice, so it accepts many spelling/notation
    of these language names.

    The return value are typically strings or tuples, dicts. There should
    be no side effects.
    """
    import re
    from .cpp_exporter import to_cpp
    from .cpp_exporter import Solver as CppSolver
    from .dsl import to_traditional_dda
    from .sympy import to_sympy, to_latex
    from .scipy import to_scipy
    exporters = {
        "CppSolver": CppSolver,
        r"c(\+\+|pp)?": to_cpp,
        "dda": to_traditional_dda,
        "sympy": to_sympy,
        "latex": to_latex,
        "scipy": to_scipy,
    }
    for k,v in exporters.items():
        if re.match(k, to, re.IGNORECASE):
            return v(state, **kw)
    raise ValueError(f"Export format {to} not known. Valid (regexps) are {list(exporters.keys())}.")

def clean(thing, target="C"):
    """
    Cleans an identifier for being compatible with the `target` language. This can be something like
    `C`, `python` or `dda` (cf. languages supported by :py:meth:`dda.export`) or also `latex`.
    
    It will basically try to transliterate all Unicode to ASCII and then try to ensure that the
    identifier is a valid C variable name (i.e. don't start with numbers, etc.).
    
    This function is nice, if you pass a :py:class:`dda.State` or :py:class:`dda.Symbol`, it will
    map the whole State/Symbol. Otherwise, it expects a string.
    
    Examples:
    
    >>> clean("\\frac{x}{y}")  # backslashes are just removed   # doctest: +SKIP
    'fracxy'
    >>> clean(r'a^{-1}')
    'a__1'
    >>> clean('a^b_c^{ef}')
    'a_b_c_ef'
    >>> clean(u'µ²')   # only if python package "unidecode" is installed # doctest: +SKIP
    'u2'
    >>> clean('77%alc')   # well, you can use numbers at the beginning of strings
    '_77alc'
    """
    import re, functools
    wants = lambda expr: re.match(expr, target, re.IGNORECASE)

    if isinstance(thing, State) or isinstance(thing, Symbol):
        return thing.map_heads(functools.partial(clean, target=target))
    elif not isinstance(thing, str):
        raise TypeError(f"Can only operate on strings (or States, Symbols, for convenience), but got '{thing}', which is of type {type(thing)}.")
    else:
        identifier = thing
        original = identifier # keep a copy

    # First, compose a decoder to remove any non-ascii stuff:
    try:
        from unidecode import unidecode
        identifier = unidecode(identifier)
    except ImportError as e:
        identifier = identifier.encode("ascii", errors="ignore").decode()

    if wants(r"c(\+\+|pp)?|python") or wants("dda"):
        if re.match(r"^\d", identifier):
            identifier = "_" + identifier

        identifier = identifier.replace("{","").replace("}","").replace("^","_").replace("\\","")

        # a^{-1} gets a_m1 and -a gets _a but a-b-c gets a_b_c
        identifer = re.sub(r"(\W)-", r"\1m", identifier)
        identifer = re.sub(r"^-", r"m", identifier)
        identifier = identifier.replace("-","_")

        # Remove anything left
        allowed = "_0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXY"
        identifier = "".join(c for c in identifier if c in allowed)
    elif wants("tex"):
        # Don't do anything, expect the user to be able to write valid latex identifiers.
        pass

    if not identifier: # empty string?
        raise ValueError(f"Identifier '{original}' is so weird that it was chopped to an empty string by the clean() function. Please choose a more subtle identifier.")

    return identifier

# populate namespace with some useful objects which should
# live in the package namespace.

from .ast import Symbol, State, BreveState, symbols, is_symbol
from .computing_elements import dda, dda_functions, dda_symbols
from .sympy import to_sympy
from .dsl import read_traditional_dda, to_traditional_dda, read_traditional_dda_file

#
#
#
# 
#
 
"""
The ``dda.scipy`` module allows in-python evaluation of DDA systems
as well as their solution with ODE Integrators in Python,
such as `scipy.integrate <https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html>`_.
For the usage and examples, see the main class :class:`to_scipy`.

.. note::

   In order to run this code, you need, obviously, `SciPy <https://www.scipy.org/>`_ next 
   to `NumPy <https://numpy.org/>`_.
   
.. warning::

   This module exposes a solver of a DDA system which is quite different to the
   :mod:`cpp_exporter`. In particular,
   
   * The solver is required to be told the solution time span or final time in
     time units, not iteration indices.
   * The solver only spills out the evolved (integration) quantities and not any
     derived quantities. You can recompute them at any timestep, but there is currently
     no code helping you in order to archieve this result. This can result in confusion
     when you cannot query the fields you asked for (in particular in the CLI frontend).
   * The SciPy time integrator tries to find an optimal (and minimal) time step, yielding
     in a quite "rough" solution. You can turn on *dense output* in order to tell the
     SciPy solver to integrate between these time steps, yielding a more smooth output
     with more datapoints.   
"""

from . import State, Symbol, symbols, clean
import math, functools, itertools, operator, copy, argparse, inspect, sys # builtins
import numpy as np # externals
from builtins import sum # just to be explicit, sum is not overwritten
from collections import namedtuple, OrderedDict

first = lambda iterable: next(x for x in iterable if x)  # first element of iterator/list

dda2python = {  #: Mapping directory of Symbol heads to their pure-python implementation.
    "const": lambda x: x,
    "neg":   lambda x: -x,
    "div":   lambda x,y: x/y,
    "Int":   lambda *x: - sum(x), # mangled integral
    "sum":   lambda *x: - sum(x),
    "mult":  lambda x,y: x*y,
    "sqrt":  lambda x: math.sqrt(x),
    "abs":   lambda x: math.abs(x),
    "exp":   lambda x: math.exp(x),
    "floor": lambda x: math.floor(x),
}

def evaluate_values(smbl, values):
    "Evaluate a symbol within the context of an already evaluated values dictionary given."
    if isinstance(smbl, float) or isinstance(smbl, int):
        return smbl # usable node
    if not isinstance(smbl, Symbol): raise TypeError(f"Expecting symbol, got {smbl}")
    if smbl.is_variable():
        return values[smbl.head] # can raise KeyError if variable not found.
    else: # symbl.is_term()
        if smbl.head in dda2python:
            return dda2python[smbl.head](*(evaluate_values(t, values) for t in smbl.tail))
        else:
            raise ValueError(f"DDA Symbol {smbl.head} in expression {smbl} not (yet) implemented.")

class to_scipy:
    """
    The SciPy exporter. When initializing this class with your DDA system, it will setup
    a function ``f(y)`` which can be evaluated as any right hand side in the ordinary
    differential equation (ODE) system ``dy/dt = f(y)``. Here, ``y`` are the evolution
    quantities, i.e. a vector which is composed automatically from the linearized DDA system
    (see :meth:`ast.State.name_computing_elements` and :meth:`ast.State.variable_ordering`).
    Furthermore, this class prepares the initial values ``y0`` for the integration.
    
    You can evaluate these quantities in any python context, i.e. with any scientific python
    ODE solver library. For the time being, this class provides a convenience method
    :meth:`solve` which calls `scipy.integrate.solve_ivp <https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html>`_.
    There is no other scipy dependence in this code.
    
    Usage example:
    
    >>> from dda.computing_elements import neg,int,mult
    >>> from dda import State, symbols
    >>> dda_state = State({x: int(y, 1, 1), y: mult(2,x), z: neg(int(z, 1, -1)) })
    >>> py_state = to_scipy(dda_state)
    >>> py_state.state # state has been linearized before processing
    State({'int_1': Int(z), 'x': Int(y), 'y': mult(2, x), 'z': neg(int_1)})
    >>> py_state.vars.evolved # evolved variables are therefore not [x,z] but [int_1,x]
    ['int_1', 'x']
    >>> py_state.y0  # initial values
    array([-1,  1])
    >>> py_state.dt # same timestep for all integrals (but see note below)
    1
    >>> py_state.rhs(py_state.y0)  # evaluation of f(y) on y0
    array([-1, -2])
    >>> y1 = py_state.rhs(py_state.y0) * py_state.dt  # a single Euler integration timestep, for instance
    >>> y1
    >>> sol = py_state.solve(10)  # ODE integration with SciPy
    >>> sol.t  # integration went from 0->10 with 17 timesteps     # doctest:+ELLIPSIS
    array([ 0.        , ..., 10.        ])
    >>> sol.y[:,-1] # the first solution is ~exp(+t)->inf, the second exp(-t)->0
    array([-2.20269685e+04,  1.94334984e-08])
    >>> from pylab import plot, legend, show  # plotting example   # doctest:+SKIP
    >>> for i,fieldname in enumerate(py_state.vars.evolved):       # doctest:+SKIP
            plot(sol.t, sol.y[i], label=fieldname)
    >>> legend(); show()                                           # doctest:+SKIP
    
    .. warning:: Due to the way how widespread ODE integrators work, the per-integral
       step size ``dt`` is required to be the same for every integration which is part
       of the DDA system. That is, the following generates an error:
       
       >>> from dda import dda, symbols
       >>> a,b=symbols("a,b")
       >>> state = State({ a: dda.int(a, 0.2, 0), b: dda.int(b, 0.5, 0) })
       >>> to_scipy(state)
       ...
       ValueError: Scipy requires all timesteps to be the same, however dt_(['a', 'b']) = [0.2 0.5]
       
       Most high-level integrators available in scientific Python toolkits (such as scipy)
       assume the overall system to have a single timestep size ``Δt`` (which is also quite
       natural from a mathematical perspective). The signature ``dda.int(f,dt,ic)`` is thus
       quirky from a mathematical or numerical viewpoint. It is written in such a way because
       analog computing integrators have a tunable time scale ``k_0 ~= 1/dt`` which however
       can also be consumed in the integrand itself: ``dda.int(f,1/k_0,ic) == dda.int(f/k_0,1,ic)``.
       
       Furthermore, most high-level integrators do adaptive timestepping anyway. The fine-tuning
       of timestep sizes is something which is only paid respect to in the :mod:`cpp_exporter`
       module.
    """
    
    def __init__(self, state):
        "state needs to be an instance of :class:`ast.State`"
        self.state = clean(state, target="python").name_computing_elements()
        self.vars = self.state.variable_ordering()
        
        if not self.vars.evolved:
            raise ValueError("Nothing to evolve. Lacking some dda.int(f,dt,ic) integrators")
        
        # by intention, dt and initial condition are ignored

        # Extract int(..., timestep, initial_data) and mark integrator as visited (mangled)
        timesteps = {}
        initial_data = {}
        def map_and_treat_integrals(var):
            if not var in self.vars.evolved: return self.state[var]
            tail = self.state[var].tail
            if not len(tail) >= 3: raise ValueError("int(...) requires at least int(value, dt, ic)")
            timesteps[var] = self.evaluate_const(tail[-2])
            initial_data[var] = self.evaluate_const(tail[-1])
            return Symbol("Int", *tail[0:len(tail)-2])
        self.state = State({ var: map_and_treat_integrals(var) for var in self.state })
        
        # Scipy needs numpy arrays for further processing
        timesteps = np.array([timesteps[k] for k in self.vars.evolved])
        self.y0 = np.array([initial_data[k] for k in self.vars.evolved])

        # The Scipy integrate methods cannot treat different timestep sizes.
        self.dt = timesteps[0]
        if not np.all(timesteps == self.dt):
            raise ValueError(f"Scipy requires all timesteps to be the same, however dt_({self.vars.evolved}) = {timesteps}")
        
        # set beginning values for each call of f(y). This dictionary can be
        # mutated over repeated calls, but if you want to go sure, make a 
        # deep copy at every call.
        self.evaluation_default_values = { k: np.nan for k in self.state.keys() }
        for var in self.vars.explicit_constants:
            self.evaluation_default_values[var] = self.state[var].tail[0]
            
        self.debug = False
    
    def evaluate_const(self, var):
        """
        Translate ``const(foo)`` by stripping foo or some bar by looking up if it is an explicit constant.
        Dynamical variables are not allowed here. This is somewhat similar but different to
        ``cpp_exporter.lookup_const(var)``.
        """
                if not var.head in self.vars.explicit_constants:
                    raise ValueError(f"Only constants allowed in this context. {var} however refers to {var.head}.")
                return self.evaluate_const(self.state[var.head])
        return var
    
    def evaluate_state(self, evolution_vector, copy=False):
        """
        Recomputes the full state from the evolution state vector.
        Returns a dictionary with same keys as ``self.state`` and scalars (floats) as
        values.
        
        This will especially compute the aux variables, while for the evaluation
        variables the RHS of ``dy / dt = f(y)`` is computed.
        
        .. note:: As a user, you most likely want to call :meth:`reconstruct_state` or
           :meth:`rhs` instead of this function.
        
        For optimization purpose, numerical state evaluation is always carried
        out on the ``evaluation_default_values`` member (which also hold the initial
        values for the first ``rhs`` evaluation). If you set ``copy=True``,
        a shallow copy (which is equal to a deep copy for a dict holding floats)
        is returned. In external calls, you should probably always set ``copy=True``.
        
        .. note:: The implementation of this function currently evaluates the (prepared) DDA 
            sytem by recursive calls with the help of a variable assignment directory. This is
            basically a *run-time compilation* (JIT/VM) in pure python. Needless to say, this
            won't give a great performance!
            
            There are plenty of low-hanging fruits to provide optimized versions of this code:
            One could call the efficient (but still scalar) C++ implementation which
            :mod:`cpp_exporter` provides by `methods provided by Cython
            <https://cython.readthedocs.io/en/latest/src/userguide/external_C_code.html>`_.
            One could also map the DDA abstract syntax tree (AST) to the python one and
            use some unparser code to evaluate the DDA code as pure python
            (see for instance
            `Python: Modify AST and write back python code 
            <https://stackoverflow.com/questions/768634/parse-a-py-file-read-the-ast-modify-it-then-write-back-the-modified-source-c>`_).
            
            For the time being, this code remains as a pure demonstration code.
            Thanks to using the linearized state, there should be no troubles with call stack
            overflows, however cyclic dependencies may not be properly resolved
            and can result in infinite recursions (stack overflow).
        """
        values = self.evaluation_default_values
        # A shallow copy should be fine, since self.evaluation_default_values 
        # only holds immutable datat types (floats) and never mutable ones
        # (such as lists). We allow for *not* doing this copy only for performance
        # reasons, as calls to self.rhs() basically overwrite the values dict
        # completely and thus don't need to reallocate all the dict on every call.
        if copy or self.debug: values = values.copy() #copy.deepcopy(values)
        for i,var in enumerate(self.vars.evolved):
            values[var] = evolution_vector[i]  # scatter
        # TODO: md3.py doesnt work when vars.aux.unneeded is skipped, because
        #       Int(...) depends on such an unneeded guy. This is wrong.
        for var in itertools.chain(self.vars.aux.sorted, self.vars.aux.cyclic, self.vars.aux.unneeded):
            values[var] = evaluate_values(self.state[var], values)
        nanmask = np.isnan(list(values.values()))
        if any(nanmask):
            nankeys = np.array(list(values.keys()))[nanmask]
            print(f"to_scipy(state).evaluate_state({evolution_vector}): NaN detected in evaluation of {nankeys}")
            print("Here is a PDB shell for debugging:")
            from pdb import set_trace as bp
            bp()
        return values
    
    def reconstruct_state(self, evolution_vector, copy=True):
        """
        Given the evolution vector sizes, this computes the full state. That is, this function
        differs from :meth:`evaluate` at all evaluation quantities where the values of the
        evolution vector itself are put in place.
        """
        values = self.evaluate_state(evolution_vector, copy=copy)
        for i,var in self.vars.evolved:
            values[var] = evolution_vector[i]

    def rhs(self, evolution_vector):
        """
        The ODE *right hand side* in ``dy / dt = f(y)``. ``y`` is a numpy vector,
        and ``f(y)`` returns a similarly sized (numpy) vector which we call ``rhs`` here:
        
        >>> ode = to_scipy(State({ Symbol("x"): Symbol("int")(Symbol("x"),0.1,1) }))
        >>> y1 = ode.y0 + ode.rhs(ode.y0) * ode.dt   # perform some euler integration step
        >>> y1
        array([0.9])
        
        Usually, you want to pass this function to some scipy integrator. See
        also :meth:`ft`.
        """
        values = self.evaluate_state(evolution_vector)
        dqdt = [ evaluate_values(self.state[var], values) for var in self.vars.evolved ] # gather
        if np.any(np.isnan(np.array(dqdt))):
            bp()
        return np.array(dqdt)

    def rhst(self, t, evolution_vector):
        """Syntactic sugar for scipy integrators who want a signature ``rhst(t,y)``. Will
        just call ``rhs(y)`` instead."""
        return self.rhs(evolution_vector)

    def solve(self, tfinal, **kwargs):
        """
        Basically passes all arguments to ``scipy.integrate.solve_ivp``. See
        documentation for :class:`to_scipy` for usage example.
        
        Currently, it is hardcodedly ``tspan=[0,tfinal]``. All other (keyword) arguments
        are passed to ``solve_ivp``.
        """
        from scipy.integrate import solve_ivp
        return solve_ivp(self.rhst, [0, tfinal], self.y0, **kwargs)

def cli_scipy():
    """
    A Command Line Interface (CLI) for :mod:`dda.scipy`.

    This CLI API basically solves a DDA file (see :mod:`dda.dsl` for the syntax).
    This is a different approach then using the :mod:`dda.cpp_exporter` Instead
    of code generation (and the need for a C++ compiler), this evaluates the
    DDA file within python. The disadvantage is that this is damned slow, the
    advantage is that the time integrator is much better then the selfmade one
    in the :mod:`dda.cpp_exporter` module. And there is no need for a C++
    compiler at all, all is (more or less) pure python.
    
    Invocation is like  ``python -m dda.scipy --help`` anywhere from the system.
    Run this to explore the usage of this command line interface.
    
    The output will be CSV (file or stdout), in terms of one line per integration
    step (called *dense solution* in scipy ODESolver language).
    """
    from .dsl import read_traditional_dda
    
    parser = argparse.ArgumentParser(description="PyDDA's scipy interface simulation runner", epilog=inspect.getdoc(cli_scipy))

    parser.add_argument("-o", "--output", nargs='?', type=argparse.FileType('w'), default=sys.stdout, help="Where to write output CSV to. Default is stdout.")
    parser.add_argument("-q", "--query-fields", nargs='*', help="List of fields to plot. Just pass whitespace seperated (i.e. -q a b c). Also add 't' if you want to have the solution time (recommended).")
    
    g = parser.add_argument_group(title="Arguments passed to scipy.integrate.solve_ivp")
    g.add_argument("-t", "--tfinal", required=True, type=float, help="Time (in simulation units) to run up to. Do not confuse this with some iteration counter.")
    g.add_argument("-m", "--method", nargs="?", default=False, help="Integration method to use")
    g.add_argument("-d", "--dense-output", action='store_true', help="Dense Output (default is not dense)")
    # add whatever you want to expose and add them to the passing list below:
    scipy_args = [ "tfinal", "method", "dense_output" ]
    
    # FIXME: Dense output does not yet seem to work.
    # FIXME: Leaving --method out does not yet work!
    
    dda_state = read_traditional_dda(dda_text)
    scipy_state = to_scipy(dda_state)
    
    scipy_args = { k: vars(arg)[k] for k in scipy_args if k }
    sol = scipy_state.solve(**scipy_args)
    
    named_sol = OrderedDict()
    named_sol["t"] = sol.t
    for i,fieldname in enumerate(sorted(scipy_state.vars.evolved)):
        named_sol[fieldname] = sol.y[i]
        
    if arg.query_fields:
        try:
            writeout = OrderedDict((q, named_sol[q]) for q in arg.query_fields)
        except KeyError as e:
            raise KeyError(f"Make sure that your query fields {arg.query_fields} are within the available keys {list(named_sol)}.") from e
    else:
        writeout = named_sol

    sep = "\t"
    # don't using np.savetxt(header=...) because it prepends a comment sign "# "
    print(sep.join(writeout.keys()), file=arg.output)
    np.savetxt(arg.output, np.array(list(writeout.values())).T, delimiter=sep)

    cli_scipy()

#
#
# This file is part of the FPAA module of the PyAnalog toolkit.
#
# 
#

# Python-included
import sys, os, argparse, glob, logging
from math import ceil
from copy import deepcopy
from collections import OrderedDict, namedtuple, defaultdict
from pprint import pprint, pformat
from collections.abc import Iterable
from pathlib import Path
from functools import reduce
from numbers import Number

# Helper routines for simply nested dictionaries:
flatten_dict = lambda dct: dict(reduce(lambda a,b: a+b, [list(dict(k).items()) for k in dct.values()]))
filter_dict = lambda dct: { k:v for k,v in dct.items() if v }
# Map [{'I1':'a'},{'I2':'b'},...] -> [('I1','a'),('I2','b'),...]
Target = namedtuple("Target", ['part','pin'])
pin2tuple = lambda dct: [ Target(t,p) for (t,p) in dct.items() ][0]
pins2tuples = lambda pl: [ Target(t,p) for tp in pl for (t,p) in tp.items() ]
# no need for bitarray
bool2bin = lambda boolean: '1' if boolean else '0'
int2bin = lambda number: bin(number)[2:] # cutting away the 0b from 0b10101
boolList2BinString = lambda lst: ''.join(map(bool2bin, lst))
bitstring2bin = lambda s: int('0b'+s, base=2)

def chunks(lst, n):
    for i in range(0, len(lst),n):
        yield lst[i:i+n]

module = __file__
log = logging.getLogger(module)

info = log.info
debug = log.debug

machines_from_list = True # choose Machines from current directory instead of providing YAML file

def yaml_load(fname):
    # external dependencies, install with "pip install pyyaml"
    # If you don't have pip, install pip with "easy_install pip"
    import yaml # PyYAML
    
    with open(fname, "r") as cfh:
        return yaml.load(cfh, Loader=yaml.SafeLoader) # may rise ScannerError

architectures_basedir = os.path.dirname(os.path.realpath(__file__))
available_architectures = {Path(fn).stem: fn for fn in glob.glob(architectures_basedir+"/*.yml")}

def load_from_yaml(circuit, arch):
    """
    Expects arch and circuit to be strings.
    """
    if not Path(arch).is_file():
        if arch in available_architectures:
            arch = available_architectures[arch]
        else:
            raise ValueError(f"{arch} is neither a file nor a registered architecture.")
    arch = yaml_load(arch)
    circuit = yaml_load(circuit)
    return (circuit, arch)

def synthesize(circuit, arch):
    """
    Translate a circuit to a netlist for a given target architecture.
    
    This routine is the heart of the FPAA compiler. It mainly
    
    * Allocates available hardware to the requested ones by the user circuits
      and allows for book keeping between the user named and architecture named computing
      components.
    * Ensures consistency of the resulting circuit (there are no dangling wires,
      no over-allocation, etc.)
      
    There is a lot of *info* and *debug* output available if turned on via Python
    logging.
    
    Expects ``arch`` and ``circuit`` to be nested data structures (dicts and lists holding
    strings and numbers), similar to their YAML representation. The documentation does
    not yet cover an in-depth description of these data structures, but there are tons
    of example YAML files which are straightforward to understand.
    
    Returns ``wired_circuit``, a list of computing components (which itself are again
    "PODs", i.e. dicts with nested data structures).
    """

    info("Welcome to the HyConAVR.ino Program Compiler.")
    info("Input program: "  + circuit['title'])
    info("Target machine: " + arch['title'])

    debug("arch:")
    debug(pformat(arch))
    debug("-----------------------------------------------------")
    debug("circuit:")
    debug(pformat(circuit))

    assigned_parts_by_entity = filter_dict({
        entity: OrderedDict({ part: None for part,parch in arch['configurable_parts'].items()
            if parch['type']==entity and not 'cannot_be_allocated' in parch })
        for entity in arch['entities'].keys() })

    # Determine mapping of part names (circuit -> arch)
    for part, spec in circuit['program'].items():
        t = spec['type']
        if not t in arch['entities']:
            raise ValueError(f"Invalid type {t} for Part {part} in Netlist {circuit['title']}. Available types for given architecture {arch['title']} are: {', '.join(arch['entities'].keys())}")
        assigned_parts = assigned_parts_by_entity[t]
        none_allocated_parts = [k for k,v in assigned_parts.items() if not v ]
        #import ipdb; ipdb.set_trace() 
        if len(none_allocated_parts) == 0:
            raise ValueError(f"Have used up all {len(assigned_parts)} parts of type {t} in architecture {arch['title']}! Cannot allocate another one.")
        target = none_allocated_parts[0]
        info(f"Allocating Type {t}: Mapping circuit part {part} onto architecture part {target}")
        assigned_parts[target] = part

    # Mapping from architectured parts to user-named parts (having also None's for unallocated parts)
    assigned_parts = flatten_dict(assigned_parts_by_entity)
    # Identity mapping for architectured parts
    default_mapping = { k:k for k in arch['configurable_parts'].keys() }
    # Extend the architectured parts with basically the nonallocable ones
    arch2user = {**default_mapping, **assigned_parts}
    user2arch = { v:k for k,v in arch2user.items() }
    # Note that arch2user is not always invertible. Several arch parts could
    # not be used and resolve to None.

    # Setup the wired circuit. In this dictionary, the parts are named as in
    # the architecture and *not* as from the user view. Use arch2user to translate
    # the user view, i.e. access like wired_circuit[arch2user[userpartname]]...
    wired_circuit = deepcopy(arch['configurable_parts'])
    for part in filter_dict(assigned_parts).keys():
        wired_circuit[part].update(circuit['program'][assigned_parts[part]])

    def resolve_user_pin(item): # closure over arch
        """
        Excepts a user-named part and always returns an architecture part name.
        Expands something like 'M2' to {'M2':'o'}, i.e. part -> {part:signal},
        where the default signal is the first output signal.
        (Note: This is a closoure over wired_circuit in a pseudo-OOP fashion.)
        """
        if isinstance(item, str):
            item = user2arch[item]
            return {item: arch['entities'][wired_circuit[item]['type']]['output'][0]['name'] }
        elif (isinstance(item, dict) and len(item)==1):
            return { user2arch[k]: v for k,v in item.items() } # could also convert to Target() at this place
        elif isinstance(item, Number):
            return item
        else: raise ValueError(f"Malformed target notation: {item}")

    # Sweep over all configurable parts
    for pname, part in wired_circuit.items():
        # Prepare for giving output information in next sweep
        part['output'] = defaultdict(list)

        if not 'input' in part:
            part['input'] = []

        reference = arch['entities'][part['type']]

        # ensure the inputs are dictionaries
        if isinstance(part['input'], list):
            new_input = {}
            for i, inputdesc in enumerate(reference['input']):
                if i < len(part['input']):
                    new_input[ inputdesc['name'] ] = part['input'][i]
            part['input'] = new_input

        # fill up defaults
        for name, value in reference.get('default_inputs', {}).items():
            if name not in part['input']:
                part['input'][name] = value

        # Resolve variables (inteded for numerical values):
        if 'coefficients' in circuit:
            for name, target in part['input'].items():
                if isinstance(target, str) and target in circuit['coefficients'].keys():
                    part['input'][name] = circuit['coefficients'][target]
                    info(f"Resolving variable {target}={part['input'][name]} at architecture part {pname}/{name}")

        # Name implicit (first) output lines
        for name, target in part['input'].items():
            try:
                part['input'][name] = resolve_user_pin(target)
            except KeyError: # thrown by user2arch
                #debug(f"Will pass malformed user target '{target}' at users {pname} to next sweep")
                userdesc = f"Architecture part {pname} (User part {arch2user[pname]})"
                raise ValueError(f"{userdesc}, input {name}: Cannot understand {target}, certainly because it is nonexistent")
            #if isinstance(target, str):
            #    part['input'][name] = {target: arch['entities'][wired_circuit[target]['type']]['output'][0]['name'] }

    # final sweep over all configurable parts:
        userdesc = f"Architecture part {pname} (User part {arch2user[pname]})"
        # Check wire types
            adesc = { dct['name']: dct for dct in arch['entities'][part['type']]['input'] }
            if not name in adesc:
                raise ValueError(f"{userdesc} constructs input line {name} which doesn't exist for type {part['type']}")
            adesc = adesc[name]
            if adesc['type'] == 'numeric' and not isinstance(target, Number):
                raise ValueError(f"{userdesc} requires a number, but {target} given. (Hint: Maybe you used an undefined variable)")
            if isinstance(target, dict):
                if len(target) > 1:
                    raise ValueError(f"{userdesc} contains too many information. {target} given")
                (tpart,tline), = target.items()
                tpart_reference = {dct['name']:dct for dct in arch['entities'][wired_circuit[tpart]['type']]['output'] }
                if not tline in tpart_reference:
                    usertarget = { arch2user[k]:v for k,v in target.items() }
                    raise ValueError(f"{userdesc} wires to nonexisting target in input {target} (User provided {usertarget})")
                tpart_reference = tpart_reference[tline]
                if tpart_reference['type'] != adesc['type']:
                    raise ValueError(f"I{userdesc}: Incompatible target line {name}. Required type: {tpart_reference['type']}, but lined to {adesc}")

                # Give output information, because we can.
                wired_circuit[tpart]['output'][tline].append({pname:name})

        # Check if everything is given
        missing_keys = set([i['name'] for i in arch['entities'][part['type']]['input']]) \
                    - set(part['input'].keys())
        if missing_keys:
            raise ValueError(f"{userdesc}: Too few input lines given: Missing keys {missing_keys}")

    # really final sweep: Ensure nonused parts have no output
        if len(part['input']) == 0 and len(part['output']) > 0 and not 'cannot_be_allocated' in part:
            raise ValueError(f"{userdesc} has no input but is wired to {part['output']}. The universe will collapse into a black hole!")

    # As a service, for the time being, provide the readout/measurement positions as info
    if "observables" in circuit:
        for var, utarget in circuit["observables"].items():
            (tpart,tline), = resolve_user_pin(utarget).items()
            try: address = arch["configurable_parts"][tpart]['address'][tline]
            except KeyError:
                # should check before if tpart actually has output named tline.
                raise ValueError(f"At requesting output for {var}: Archtecture missing address for output line {tline} in part {tpart} in section configurable_parts")
            info(f"Observable {var} (user defined part {utarget}) can be measured at architecture part output {tpart}:{tline}, adress 0x{address:04x}")

    # yay
    #info("wired_circuit: ")
    #info(pformat(wired_circuit))
    #sys.exit(-1)
    return wired_circuit

def normalize_potentiometer(value, resolution_bits=10):
    "Map a real value [0..1] to Potentiometer value [0..1023]"
    maxval = 2**resolution_bits - 1 # 2**10-1 = 1023
    value = float(value)
    if value < 0 or value > 1:
        raise ValueError(f"Digital potentiometer value {value} out of bounds")
    return int(round(value * maxval))

last_seen_xbars = [] # an ugly global, filled by compile_instructions; for later plotting

def compile_instructions(wired_circuit, arch):
    """
    Compile a netlist (``wired_circuit``) to configuration *instructions* for setting
    up the analog computer hybrid controller.
    
    This routine basically loops over the *hardwired parts* of the given architecture,
    i.e. built-in
    
    * potentiometers (DPT-24 and HC)
    * cross bar arrays (XBAR)
    
    and configures them according to the given wired_circuit. This means that relevant
    allocated potentiometers will be set and the XBAR configuration bitmask will be
    computed from the hardware description provided by the circuit and the architecture.
    
    Currently returns a list of instructions (tuples) which could be directly be written
    out to serial or passed to PyHyCon.
    
    .. note::
       This method will change in near time and talk directly to a PyHyCon instance.
    """
    instructions = []

    # should actually call write(tpl) or directly PyHyCon
    instruct = lambda *tpl: instructions.append(tpl)

    def resolve_machine_pin(item): # closure over arch
        """
        Expects an architecture part name. Never comes in touch with user names.
        """
        elif (isinstance(item, dict) and len(item)==1) or isinstance(item, Number):

    # Go over hardwired parts
    for hwname, hw in arch['wired_parts'].items():
        if hw['type'] == "DPT24":
            # DPT24 Potentiometers
            assert len(hw['enumeration']) <= 24, "DPT24 has only 24 digital potentiometers"
            for port, t in enumerate(pins2tuples(map(resolve_machine_pin,hw['enumeration']))):
                numeric_value = wired_circuit[t.part]['input'][t.pin]
                normalized_value = normalize_potentiometer(numeric_value)
                info(f"DPT24@{hw['address']:x}: Storing value {'%4d'%normalized_value} at DPT port {port:2} (corresponding to {t.part}:{t.pin})")
                #print(f'$ac->set_pt("DPT24-{port}", {numeric_value});')
                instruct("P", hw['address'], "%02X"%port, "%04d"%normalized_value)
        elif hw['type'] == 'HC':
            # Hybrid controller: DPTs (same code as DPT24)
            assert len(hw['dpt_enumeration']) <= 8, "HC has only eight digital potentiometers"
            for port, t in enumerate(pins2tuples(map(resolve_machine_pin,hw['dpt_enumeration']))):
                info(f"HC@{hw['address']:x}: Storing value {normalized_value:4} at DPT port {port:2} (corresponding to {t.part}:{t.pin})")
                #print(f'$ac->set_pt("HCDPT-{port}", {numeric_value});')
            # Hybrid controller: Digital output
            assert len(hw['digital_output']) <= 8, "HC has only eight digital outputs"
            for port, t in enumerate(pins2tuples(map(resolve_machine_pin,hw['digital_output']))):
                value = wired_circuit[t.part]['input'][t.pin]
                info(f"HC@{hw['address']:x}: Storing {value} at digital output port {port} (corresponding to {t.part}:{t.pin})")
                instruct("D" if value else "d", hw['address'], "%1d"%port)
                #print(f'$ac->digital_output({port}, %d);' % (1 if value else 0))
        elif hw['type'] == 'XBAR':
            # XBAR matrix
            N,M = len(hw['output_rows']), len(hw['input_columns'])
            assert N==16 and M==16, "XBAR only implemented for 16x16"
            info(f"XBAR@{hw['address']:x}: Computing XBAR of size NxM={N}x{M}")
            cols = pins2tuples(map(resolve_machine_pin, hw['input_columns']))
            rows = pins2tuples(map(resolve_machine_pin, hw['output_rows']))
            #outputs = { pname: dict(part['output']) for pname,part in wired_circuit.items() }
            inputs = { pname: { line: pin2tuple(target) for line,target in part['input'].items()
                    if isinstance(target,dict) } #and not "None" in target } # filter out empty inputs
                    for pname,part in wired_circuit.items() }

            # The AD8113 enforces that there is only one connection per (output) row.
            # In other words: In the XBAR, an output line can be connected only to one input
            #   line, but an input line in the XBAR can connect up to 16 outputs.
            # This is realized by having an output row being encoded in only 4 bits instead of 16.

            boolean_matrix = [[ ol in inputs[op] and Target(ip,il) == inputs[op][ol] and op!="None" and ip!="None" for (ip,il) in cols] for (op,ol) in rows]
            row_bitstrings = list(map(boolList2BinString, boolean_matrix))
            row_numbers = [ row.index(True) if sum(row) else 0 for row in boolean_matrix ]
            row_active = [sum(row)==1 for row in boolean_matrix]
            row_bitstring = [ f"{active:b}{num:04b}" for num,active in zip(row_numbers, row_active) ]
                    
            for i,(bitvec,num,active,bitvec2,(op,ol)) in enumerate(zip(row_bitstrings,row_numbers,row_active,row_bitstring,rows)):
                info(f"XBAR@{hw['address']:x}: Writing bitmatrix[row {i:2}]:",
                    f"{bitvec}={num:2d}=0x{num:1x} -> {op}:{ol}       [sending {bitvec2}]" if active else
                    f"{bitvec} [output not enabled] [sending {bitvec2}]")

            if not all([sum(row) in (0,1) for row in boolean_matrix ]):
                raise ValueError("XBAR matrix is unsuitable. See info output for it's values. Only a maximum of one `True` bit per row allowed.")
            
            last_seen_xbars.append( (cols,rows,boolean_matrix) ) # for later plotting...

            # Caveat 1: Chip expects rows in order row15...row0
            # Caveat 2: Don't try to convert a bitstring of length 80 to a single int, it will overflow.
            bitstring = "".join(row_bitstring[::-1])
            assert len(bitstring)==80, "XBAR bitstring has wrong length"
            bitstring_hex = "".join(["%02X"%int(x, base=2) for x in chunks( "".join(row_bitstring[::-1]), 8)])
            info(f"Bitstream to send ({len(bitstring)} characters): {bitstring}")
            info(f"Hextream  to send ({len(bitstring_hex)} characters): {bitstring_hex}")
            instruct("X", hw['address'], bitstring_hex)
            #print(f"$ac->set_xbar('XBAR16', '{bitstring_hex}');")
        else:
            raise ValueError(f"Wired part {hwname}: Don't know what to do with type {hw['type']}.")
        
    return instructions

def plot_xbar(target_file, circuit_title, xbar_config=None, interactive_plotting=False):
    """
    Draw an the allocation of a crossbar switch array (xbar) matrix.
    
    xbar_config is a tuple with (cols,rows,boolean_matrix), and by default the last one
    from a global registry (last_seen_xbars) is taken, which is what you want.
    """
    info(f"Drawing the XBAR to {target_file}...")
    (cols,rows,boolean_matrix) = xbar_config if xbar_config else last_seen_xbars[-1]
    M,N = 16,16

    import numpy as np
    import matplotlib as mpl
    import matplotlib.pyplot as plt
    if interactive_plotting: plt.ion()

    mpl.rcParams['font.family'] = ['monospace'] # default is sans-serif
    fig = plt.figure(figsize=[7.,7.5])
    ax = fig.gca()
    plt.xticks(np.arange(M), [ f'{i} ({part}:{pin})' for i,(part,pin) in enumerate(cols) ], rotation=30, ha="left")
    plt.yticks(np.arange(N), [ f'{i} ({part}:{pin})' for i,(part,pin) in enumerate(rows) ])
    
    ax.xaxis.tick_top()
    ax.yaxis.tick_right()
    ax.set_axisbelow(True) # grid in background
    ax.set_aspect(1)
    plt.grid()
    
    Y,X = np.where(np.array(boolean_matrix))
    plt.scatter(X,Y)
    
    plt.xlim(-1,M)
    plt.ylim(-1,N)
    ax.invert_yaxis()

    plt.title(f"XBAR for {circuit_title}", y=1.18, fontweight="bold")
    plt.tight_layout()
    plt.subplots_adjust(top=0.82)

    plt.savefig(target_file)
    
def cli():
    """
    This module is callable via ``python -m fpaa`` or ``./fpaa.py``. It exposes the main
    functions on the command line which is especially helpful for debugging or
    interactively programming an analog computer from the command line.

    Call ``--help`` for all possible command line options.
    """
    parser = argparse.ArgumentParser(description="A circuit synthesizer for the HyConAVR.", epilog=__doc__)#, formatter_class=argparse.RawTextHelpFormatter
    parser.add_argument("-v", "--verbose", dest="verbose_count",
                        action="count", default=0,
                        help="increases log verbosity for each occurence.")
    parser.add_argument("-o", "--output", default="-", metavar="OUTPUT.txt", help="Put output string into file (default is '-' and means stdout)")
    parser.add_argument("-p", "--plot", metavar="OUTPUT.pdf", help="Plot crossbar switch")
    
    arch_group = parser.add_mutually_exclusive_group(required=True)
    arch_group.add_argument("-A", "--registered-arch", choices=available_architectures.keys(), default="AP-M1-Mini", help=f"Target machine architecture description: Any YAML file in directory {architectures_basedir} is available as machine") 
    arch_group.add_argument("-a", "--arch", metavar="path/to/my/MACHINE.yml", help="Target machine architecture description (any valid filename)")
    
    parser.add_argument("circuit", metavar="CIRCUIT.yml", help="The YAML file holding the circuit description")
    args = parser.parse_args()
    
    logging.basicConfig(stream=sys.stderr, level=logging.DEBUG,
                        format='%(name)s (%(levelname)s): %(message)s')
    # Sets log level to WARN going more verbose for each new -v.
    log.setLevel(max(3 - args.verbose_count, 0) * 10)
    
    chipout = open(args.output, "w") if args.output != "-" else sys.stdout
    
    circuit, arch = load_from_yaml(args.circuit, args.arch or args.registered_arch)
    wired_circuit = synthesize(circuit, arch)
    instructions = compile_instructions(wired_circuit, arch)
    
    # That was used for raw writing once. Could just write map(write, instructions).
    # However, we want to use the HyCon interface anyway in the future, so this code
    # will be superseded.
    def write(command_letter, address, *data):
        if not isinstance(address,int):
            raise ValueError("Need address as integer (may specify as 0x123)")
        command = f"{command_letter}{address:04X}" + "".join(data)
        debug(f"Writing out: {command}")
        print(command, file=chipout)
    
    print(instructions)

    if args.plot:
        plot_xbar(args.plot, circuit_title=circuit['title'])#, interactive_plotting=args.debug)
    
    try:
        cli()
    except Exception as e:
        log.error(e)  # will show stack trace only
        sys.exit(-1)  # when at least one -v is given.
    
"""
Pythonic API for Siglent Digital Oscilloscope in order to perform data aquisition.
This is a python standalone code implementing the SCPI protocol (without further
Python libraries) over sockets. The connection is assumed to be established over
Ethernet (not USB, despite this should also work, in principle, as the siglent
registers as UART device over USB). The code can work in blocking and non-blocking
fashion (but is not fully async).

See the programming guide
https://siglentna.com/wp-content/uploads/dlm_uploads/2021/01/SDS1000-SeriesSDS2000XSDS2000X-E_ProgrammingGuide_PG01-E02D.pdf
for the reference. This code was tested against an Siglent SDS1104X-E with four 10bit channels.

Here is a picture of such an oscilloscope. On the back, it has an RJ45 ethernet
connector and a USB client socket. Whether over USB-UART or TCP, it allows to be
programmed with the
`Standards Commands for Programmable Instruments (SCPI) <https://en.wikipedia.org/wiki/Standard_Commands_for_Programmable_Instruments>`_
standard, which is basically a simple command-reply text protocol similar to the
one we implement with the :class:`hycon.HyCon`.

.. image:: https://www.welectron.com/media/image/product/9607/lg/siglent-sds1104x-e-oszilloskop~3.jpg
   :width: 80%
   :alt: Photo of Siglent

Over ethernet, one can also access the oscilloscope screen via
`VNC <https://en.wikipedia.org/wiki/Virtual_Network_Computing>`_. 
The oscilloscope provides a web interface which includes an in-browser
VNC client and models the interface controls. Both is shown side on side in
the UI. The interface knobs are translated to SCPI commands via Javascript.
The status indicators are polled. The overall UI is quite heavy, rather slow
and generates, in particular, high load on the oscilloscope.

.. image:: https://www.eevblog.com/forum/testgear/siglent-sds1104x-e-in-depth-review/?action=dlattach;attach=438190;image
   :width: 100%
   :alt: Screenshot of Web Interface

    Make sure you **never ever run the browser interface at the same time as 
    a telnet session or this script**. This will break the oscilloscope randomly
    and requires you to cold start it. This is a reproducable bug in the Siglent
    firmware.
    
    It is, however, fine to run a VNC connection to the oscilloscope while
    running this script.
    
This script is also executable as standalone and offers an API to contact the siglent
and do a single readout, assuming the triggering took place independently.
"""

import socket, time, re, numpy as np, sys, argparse

# TODO: Improve logging
log = lambda *a,**b: print(*a,**b,file=sys.stderr)

def downsample(data, factor):
    """
    Downsample all numpy 1d arrays, since ``query("MEMORY_SIZE 70K")`` has no effect...
    TODO: MEMSIZE only works in STOP mode.
    """
    return { k: v[::factor] for k,v in data.items() }

class siglent:
    """
    This is a classy API to the oscilloscope. A usage example intermixed with
    HyCon controlling may look like this:
    
    .. code-block:: python
    
       import hycon.aquisition.siglent_scpi as sig
       
       hc = AutoConfHyCon("something.yml")
       print(hc.get_status())
       
       siglent = sig.siglent("192.168.32.68")
       print("Siglent Status: ", siglent.status())
       
       hc.ic()
       print("ICs freely running, preparing Siglent trigger...")
       #input("Press enter to continue")
       prepped = siglent.trigger_single()
       siglent.log("Siglent Status: ", siglent.status())
       if not prepped:
           raise ValueError("Could not properly prepare the Siglent.")
       time.sleep(1)
       hc.single_run_sync()
       time.sleep(2) # 2 seconds "OP time".
       print("Finished with single run sync. Waiting another second for siglent.")
       hc.ic()
       #input("HC stopped. Press enter to readout siglent.")
       data = siglent.read_all()
       sig.write_npz(data, f"output.npz")
    
       While this code basically works, the refactoring at PyAnalog-introduction
       time was not tested. That is, this code is not yet fully tested.

    """
    
    def send(self, md, log_end="\n", blocking=False, sleep_sec=0.5):
        "Send SCPI command to siglent."
        if blocking:
            # to ensure the oscilloscope is finished with the command...
            if not query("*OPC?", sleep_sec=0) == "1":
                raise ValueError("OPC before blocking send did not return 1")
        binary = (cmd + "\n").encode("ascii")
        log("SCPI> ", cmd, end=log_end, flush=True)
        sent_len = self.s.send(binary)
        time.sleep(sleep_sec) # 1sec suggested by the programming examples...
        if len(binary) != sent_len:
            log(f"Tried to send {len(binary)} bytes but could send only {sent_len}.")
                raise ValueError("OPC after blocking send did not return 1")

    def query(cmd, sleep_sec=0.5):
        """
        Query the siglent for something.
        This is basically a :meth:`send`/:meth:`recv` cycle, according to the
        SCPI standard.
        """
        send(cmd, log_end="", sleep_sec=sleep_sec)
        return s.recv(self.bufsize).decode("ascii", "ignore").strip() # ignore binary non-decodable

    def query_num(cmd, pat):
        "A shorthand which returns the number in some readout which looks like ``a = 123``."
        r = query(cmd, sleep_sec=0) # we never experienced race conditions here, so don't sleep
        readout = float(re.match(pat, r).group(1))
        log(" = ",readout)
        return readout
    
    def __init__(self, host, port=5025):
        """
        This class connects to SCPI Raw TCP, so use port 5025 instead of Telnet 5024.
        Insert the IPv4 address (or resolvable domain, or IPv6) of the SIGLENT oscilloscope
        as string, for instance as in ``osci = siglent("192.168.32.68")``. It is
        passed to ``socket.connect()``.
        """
        
        s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)

        # Try out if TCP_NODELAY improves something.
        s.setsockopt(socket.IPPROTO_TCP, socket.TCP_NODELAY, 1)

        address = (host, port)
        s.connect(address)
        log("SCPI Raw TCP Py3 Socket reader connected to ", target)

        self.bufsize = 4 * 1024 # bytes
        self.channels = list(range(1,5))

        # first check whether we can work
        log("= ", query("*IDN?"))

        run_full_setup = True
        if run_full_setup:
            send("*RST", blocking=True)

            send("stop", blocking=True) # attention, STOP seems to reset the oscilloscope :-(

            # Correctly setup for aquiring 4 channels and so on
            for c in self.channels:
                send(f"C{c}:TRACE ON", blocking=True)
                send(f"C{c}:VDIV 1V", blocking=True)

            send("tdiv 100ms", blocking=True)

            # set memory size *after* setting the number of channels.
            # choices: 7k, 70k, 700k, 7M
            send("MEMORY_SIZE 70K", blocking=True) # only works in STOP mode.
            log(" = ", query("MEMORY_SIZE?"))
            # ALSO: query("TRIG_DELAY -600ms") is reasonable with 100ms time division

        # TODO: all this could be encapsuled within a class.

        self.vdiv = { c: query_num("c%d:vdiv?"%c, "C%d:VDIV ([\d\.E+-]*)V"%c) for c in self.channels }
        self.ofst = { c: query_num("c%d:ofst?"%c, "C%d:OFST ([\d\.E+-]*)V"%c) for c in self.channels }
        self.tdiv = query_num("tdiv?", "TDIV ([\d\.E+-]*)S")
        self.sara = query_num("sara?", "SARA ([\d\.E+-]*)Sa/s")

    def status(self):
        """
        Reads out the status bit (``INR?``) and translates it according to the manual,
        thus returning a string.
        
        .. warning::
        
           Attention, Querying INR? is *NOT* a read-only operation but changes
           the oscilloscope status.
        """
        value = self.query("INR?")
        try:
            value = int(value.split()[-1]) # "INR 0" -> int("0")
        except ValueError:
            raise ValueError(f"Expected INR status msg but got {value=}")
        # hint: INR only returns bit 0 and bit 13...
        bits = [
            "15 Not used (always 0)",
            "14 Not used (always 0)",
            "Trigger is ready", # bit 13
            "Pass/Fail test detected desired outcome",
            "Waveform processing has terminated in Trace D",
            "Waveform processing has terminated in Trace C",
            "Waveform processing has terminated in Trace B",
            "Waveform processing has terminated in Trace A",
            "A memory card, floppy or hard disk exchange has been detected",
            "Memory card, floppy or hard disk has become full in “AutoStore Fill” mode",
            "5 Not use(always 0)",
            "A segment of a sequence waveform has been acquired",
            "A time-out has occurred in a data block transfer",
            "A return to the local state is detected",
            "A screen dump has terminated",
            "A new signal has been acquired" # bit 0
        ]
        meanings = [ desc for i,desc in enumerate(bits[::-1]) if value & (2**i) ]
        if not meanings:
            meanings = ["Zero"]
        return meanings

    def trigger_single(self):
        """
        Brings the Oscilloscope in single trig mode and checks the success of this
        operation (returns ``True`` in case of passed check).
        """
        #send("STOP", blocking=True) # somehow this resets the oscilloscope...
        print("I do *not* stop or so. Just asking for trigger mode single...")
        self.send("TRIG_MODE SINGLE", blocking=True)
        q = self.query("TRIG_MODE?")
        log(" = ",q)
        
        # Again, wait wait wait, why not wait a bit more? The oscilloscope otherwise
        # doesn't catch up...
        time.sleep(1)
        if not self.query("*OPC?") == "1":
            raise ValueError("OPC before blocking send did not return 1")
        
        return q == "TRMD SINGLE"
        """
        inr = query("INR?")
        if inr == "INR 8193":
            inr = query("INR?")
        if inr != "INR 0":
            log(f"Cannot ARM, since {inr=} but expected 0")
            return False
        else:
            log("=", query("ARM")) # ARM is send only, not query
            return True
        """

    def read_wf(self, channel, try_multiples=2):
        """
        Read the actual waveform. This can be called after an aquisition took place
        and the oscilloscope trigger fired. In many times, this will fail and return
        no data. This results in a thrown ``ValueError``. Otherwise, data recieving
        is looped (with lot's of logging since it takes some seconds). When some 
        consistency is reached (all announced data is recieved), will return the 
        *raw recieved channel data* as 1d list of numbers. That is, these are 10bit
        numbers.
        
        This function will by default try several times to ask for data, in case the
        siglent did not properly answer.
        """
        assert channel in self.channels
        
        # Wait until for instance some slow data aquisition is finished
        
        self.send("c%d:wf? dat2" % channel, log_end="")
        data = s.recv(bufsize)
        log(" = ", data[:25], "...")
        prefix = "Channel %d | " % channel
        log(prefix, f"Read (first) {len(data)} bytes.")
        
        if(not data):
            raise ValueError("Did not recieve any data.")

        # parse data format, begins like "C1:WF DAT2,#9007000000\x1e...."
        ascii_reply = data.decode("ascii", "ignore") # ignore binary non-decodable
        try:
            num_digits_idx = ascii_reply.find("#")
            num_digits = int(ascii_reply[num_digits_idx+1])
            num_bytes = int(ascii_reply[num_digits_idx+2:num_digits_idx+2+num_digits])
        except ValueError as e:
            print(f"Channel {channel}: Error at parsing data, Full reply is {ascii_reply=}")
            if try_multiples:
                print(f"Trying again for {try_multiples-1} attempts...")
                return self.read_wf(channel, try_multiples - 1)
            else:
                raise e
            
        start_payload_idx = num_digits_idx+2+num_digits

        # byte array that will grow
        all_data = data[start_payload_idx:]
        
        if not num_bytes:
            print(f"Got a zero answer ({num_bytes=}). Input was {ascii_reply=}")
            else:
                raise ValueError("Could not read any data (got 0 length).")

        log(prefix, "Recieved %9d / %d bytes (%.2f%%)" % (len(all_data), num_bytes, len(all_data)/num_bytes*100))

        # now loop until all data are recieved:
        while len(all_data) < num_bytes:
            data = s.recv(self.bufsize)
            all_data += data
            log(prefix, "Recieved %9d / %d bytes (%.2f%%)" % (len(all_data), num_bytes, len(all_data)/num_bytes*100))

        return all_data

    def read_all(self, channels=None):
        """
        Read the actual waveforms by querying the passed channel list. A channel is
        an integer between ``0`` and ``4``, see ``siglent.channels`` for the valid
        channel list. If not provided, all channels are read.
        
        Will do the conversion to voltages, that is, returns proper floats representing
        a voltage in unit V.  Returns a dictionary mapping the channel name to the 1d data.
        """
        if not channels: channels = self.channels
        channel_data = { d: self.read_wf(d) for d in self.channels }

        log("Finished reading channels:")
        for k,v in channel_data.items(): log(f"{k}: {len(v)} entries ({type(v)})")
        data_len = len(channel_data[next(iter(channel_data))])
        assert all([data_len == len(v) for k,v in channel_data.items() ]), "Non-uniform data recieved."

        log("Conversion to voltages...")
        def adc2volt(channel):
            unsigned = channel_data[channel]
            #signed = np.array([ (256-d if d > 127 else d) for d in unsigned])
            
            signed = np.frombuffer(unsigned, dtype=np.uint8).astype(np.int16)
            signed[signed > 127] = signed[signed > 127] - 256

            weird_factor = 25 # comes from documentation!
            volt = signed / weird_factor * self.vdiv[channel] - self.ofst[channel]
            
            # unfortunately, computing with floats blows up the data size in bytes
            # a lot. We basically still only have 8 bit Resolution. Floating point
            # single precision is 32 bit. Numpy supports half precision, thought
            # (16bit). Some modern numpy's also support quad precision (8bit).
            volt = volt.astype(np.float16)
            return volt

        # voltages in volt
        data = { f"channel{k}_volt": adc2volt(k) for k,v in channel_data.items() }

        # time axis in seconds
        weird_factor = 14/2
        time = -self.tdiv * weird_factor + np.arange(data_len) / self.sara
        
        data["time_sec"] = time
        
def write_npz(data, fname=None):
    """
    Write out NPZ files. See https://numpy.org/doc/stable/reference/generated/numpy.savez.html
    for this numpy specific binary file format.
    Expects data to be a python dictionary (ideally some ``np.array``).
    """
    if not fname and not sys.stdout.isatty():
        log("Writing compressed .npz file to stdout...")
        target = sys.stdout.buffer # sys.stdout does not accept binary, .buffer does.
    elif not fname:
        log("Stdout is a Terminal, writing to output.npz.")
        log("Hint: Call like ./scpi-socket.py > target.npz for other target.")
        target = open("output.npz", "wb")
    else:
        log(f"Writing npz file to {fname}...")
        target = open(fname, "wb")
        
    np.savez_compressed(target, **data) # writes binary

    parser = argparse.ArgumentParser(description="Make a read swipe on a Siglent oscilloscope and write to stdout.")
    parser.add_argument("--host", help="Host/IP address of siglent")
    siglent = siglent(args.host)
    write_npz(siglent.read_all())
 
#
#
# This file is part of the HyCon module of the PyAnalog toolkit.
#
# 
#

"""
PyHyCon -- a Python Hybrid Controller interface.

Note that the ``IO::HyCon`` Perl 5 module is the reference implementation
that is maintained by the HyConAVR firmware author (Bernd). You can find
the ``IO::HyCon`` Perl module as well as the mentioned Arduino firmware
at https://github.com/anabrid/Model-1/tree/main/software/HC. The Perl
code is also part of the CPAN at https://metacpan.org/pod/IO::HyCon.

While this implementation tries to be API-compatible with the reference
implementation, it tries to be minimal/low-level and won't implement any
client-side luxury (such as address mapping). It is the task of the
user to implement something high-level ontop of this.

Especially, the following tasks are implemented by different modules which
can but do not needed to be used:

* Connection managament: HyCon just assumes a file handle, but different
  connection types are proposed in the :mod:`connections` module.
* Autosetup: PyHyCon is plain python and has no dependency, for instance 
  on YAML. There is :mod:`autosetup` which implements the "autosetup"
  functionality of Perl-HyCon.
* High level functionality is implemented on top of HyCon and not within.
  See for instance :mod:`fpaa` for an abstraction which can generate
  HyCon instructions and is aware of the circuit design at the same time.
  
The hycon module also includes an interpreter for the HyCon serial stream
"protocol". See :mod:`replay` for further details.

Logging and Debugging
^^^^^^^^^^^^^^^^^^^^^

There are several ways to inspect what HyCon is doing. One of the simplest
is to activate logging on ``INFO`` level:

>>> import logging
>>> logging.basicConfig(level=logging.INFO)
>>> # proceed here as usual, i.e.: hc = HyCon(...)

"""

# All these modules are Python internals
import sys, re, logging, time
from collections import namedtuple

#logging.basicConfig(level=logging.INFO) # only for testing
log = logging.getLogger('HyCon') # or __name__

def ensure(var, **q):
    """
    This is our assert function which is used widely over the code for dynamic
    parameter checking. ``q`` stands for *query*. The function will return
    silently if the query is fullfilled and raise a ``ValueError`` otherwise.
    Examples for success:
    
    >>> ensure(42, eq=42)
    >>> ensure("foo bar", re="fo+.*")
    >>> ensure(17, inrange=(0,20))
    >>> ensure("x", within="xyz")
    >>> ensure("bla", length=3)
    >>> ensure("blub", isa=str, length=4, re="b.*")
    
    And in case of failure:
    >>> ensure(3, within=[1,2,9])
    Traceback (most recent call last):
    ...
    ValueError: Got var=3, but it is none of [1, 2, 9].
    """
    # Should probably use inspect.stack() or traceback.extract_stack() to get original varname
    basemsg=f"Got {var=}" if not 'basemsg' in q else q['basemsg']
    if 'eq' in q and not var == q['eq']: raise ValueError(f"{basemsg}, but should be '{q['eq']}'")
    if 're' in q and not re.match(q['re'], var): raise ValueError(f"{basemsg}, but that doesn't match regexp '{q['re']}'")
    if 'inrange' in q and not (var >= q['inrange'][0] and var <= q['inrange'][1]): raise ValueError(f"{basemsg}, but it is not in range{q['inrange']}.")
    if 'within' in q and not var in q['within']: raise ValueError(f"{basemsg}, but it is none of {q['within']}.")
    if 'length' in q and not len(var)==q['length']: raise ValueError(f"{basemsg}, of {len(var)=} but expected to be len(var)={q['length']}.")
    if 'isa' in q and not isinstance(var, q['isa']): raise ValueError(f"{basemsg}, which is of {type(var)=} but expected  type {q['isa']}.")

class expect:
    """
    ``ensure`` delayed and on stereoids: Can be initialized with a ``query`` (but with
    further options) and then called with a ``HyConRequest``. Will check the *response*
    and also allows *return value mapping*, for instance with regexpes or by splitting.
    Example:
    
    >>> R = HyConRequest("dummy")
    >>> R.response = "1,2,3"
    >>> E = expect(split=",", type=int)
    >>> print(list(E(R)))
    [1, 2, 3]
    """
    def __init__(self, **q): self.q=q
    def __call__(self, r): # r: HyConRequest
        q = deepcopy(self.q);
        q['basemsg'] = f"Unexpected response: Command {r.command} yielded '{r.response}'"
        ensure(r.response, **q)
        mapper = q['type'] if 'type' in q else lambda x:x # id
        try:
            if 'ret' in q: return mapper(re.match(q['re'], r.response).groupdict()[ q['ret'] ])
            if 'split' in q: return map(mapper, re.split(q['split'], r.response))
            if 're' in q and not 'type' in q: return re.match(q['re'], r.response)
            return mapper(r.response)
        except ValueError:  raise ValueError(f"{basemsg} but cannot be casted/mapped to {mapper}")
    def __str__(self): return "expect(%s)"%str(self.q)[1:-1].replace("'",'')

def wont_implement(reason):
    "Will not implement: Returns a function which raises ``NotImplementedError(reason)`` when called."
    def not_implemented(*v,**kw): raise NotImplementedError(reason)
    not_implemented.__doc__ = f"*Not implemented* because {reason}"
    return not_implemented

class HyConRequest:
    """
    A HyConRequest models a single *request* and *response* cycle. It stores the ASCII ``command``
    emitted by the HyCon and can save a *expected response* future/promise (see ``expect`` class).
    A HyConRequest can only be made once. If you want to do it several times, you have to (deep)
    copy the instance.
    """
    def __init__(self, command, expected_response=None):
        self.executed = False
        self.command = command
        self.expected_response = expect(re=expected_response) if isinstance(expected_response, str) else expected_response
        
        return f"HyConRequest({self.command}, {self.expected_response}, {self.executed=}, "+\
            f"response={self.response if hasattr(self, 'response') else 'n.a.'}, "+\
            f"reply={self.reply if hasattr(self, 'reply') else 'n.a.'})"
    
    def write(self, hycon):
        "Run this request. Can only be executed once. Can be chained."
        if self.executed:
            raise ValueError("Shall not execute same command twice.")
        self.executed = True
        log.info(f"Sending [{self.command}]")
        hycon.fh.write(self.command)
        return self # chainable
    
    def read(self, hycon, expected_response=None, read_again=False):
        "Read response from HyCon. If *read_again* is given, will read several times. Can be chained."
        if not expected_response: expected_response = self.expected_response
        if hycon.unidirectional:
            log.debug("Unidirectional channel, skipping reading from HyCon...")
            return self # chainable
        if not expected_response:
            log.debug("No response expected, skipping reading from HyCon...")
        log.debug(f"Waiting for response {expected_response} ... ")
        if read_again or not hasattr(self, "response"):
            self.response = hycon.fh.readline().strip() # Note: The HyConAVR always answers with a full line.
        if not self.response:
            raise ValueError(f"No Response from Hybrid controller! Command was '{self.command}'")
        self.reply = expected_response(self) # "Reply" is the highlevel, mapped answer
        return self # always chainable

class HyCon:
    """
    Low-Level Hybrid Controller OOP interface, similar to the Perl Hybrid controller.
    
    This is a minimalistic implementation which tries to implement all neccessary checking
    of input/output request/reply structure correctness, but won't do any *high level*
    support for applications. Users are assumed to write such code on themselves. The PyFPAA
    library is an example for a high level "frontend" against HyCon, which includes a
    circuit understanding, etc.
    """
    
    DIGITAL_OUTPUT_PORTS = 8
    DIGITAL_INPUT_PORTS = 8
    DPT_RESOLUTION = 10
    DPT_MAX_INT_VALUE = (2 ** DPT_RESOLUTION - 1)  # 0 <= value <= 1023
    XBAR_CONFIG_BYTES = 10
    MAX_RO_GROUP_SIZE = 500 # hardcoded in HyConAVR.ino
    
    def __init__(self, fh, unidirectional=False):
        """
        Expects fh to be an IOHandler. This could be an open file,
        a (unix/inet) domain socket, a special device (serial port),
        some serial port library, etc.
        
        .. note:: We expect the `fh` to have an API compatible to
            `ABC.IOBase <https://docs.python.org/3/library/io.html>`_. The two methods
            `fh.readline()` and `fh.write()` are required. You find some example "polyfills"
            in the :mod:`connections` module.
        
        Make sure you disable output buffering (on print/write to fh), since commands do
        not end with newlines. Responses from the uC always end with newlines.
        
        If unidirectional is set, PyHyCon won't try to read from fh. Any
        read checks and mappings will be skipped and HyConRequest.read()
        will always return None.
        
        Important: We don't send anything at construction time, unless the
        reference implementation. Call reset() on yourself if you feel to do
        so.
        """
        self.fh = fh
        self.unidirectional = unidirectional
        if unidirectional: log.info(f"Won't make any attempt to read from the HyCon at {fh}")
    
    def query(self, *args, **kwargs):
        "Create a request, run it and check the reply"
        return HyConRequest(*args, **kwargs).write(self).read(self)
    
    def command(*args, help=None, **kwargs):
        "Return a method which, when called, creates a request, runs it and checks the reply"
        method = lambda self: self.query(*args, **kwargs)
        method.__doc__ = help
        return method
    
    ic               = command('i', '^IC',            help="Switch AC to IC-mode")
    op               = command('o', '^OP',            help="Switch AC to OP-mode")
    halt             = command('h', '^HALT',          help="Switch AC to HALT-mode")
    disable_ovl_halt = command('a', '^OVLH=DISABLED', help="Disable HALT-on-overflow")
    enable_ovl_halt  = command('A', '^OVLH=ENABLED',  help="Enable HALT-on-overflow")
    disable_ext_halt = command('b', '^EXTH=DISABLED', help="Disable external HALT")
    enable_ext_halt  = command('B', '^EXTH=ENABLED',  help="Enable external HALT")
    repetitive_run   = command('e', '^REP-MODE',      help="Switch to RepOp")
    single_run       = command('E', '^SINGLE-RUN',    help="One IC-OP-HALT-cycle")
    pot_set          = command('S', '^PS',            help="Activate POTSET-mode")

    def repeated_reset(self):
        from .connections import repeated_reset
        repeated_reset(self.fh)
        
    def single_run_sync(self):
        "Synchronous run (finishes after a single run finished). Return value is true if terminated by ext. halt condition"
        q = self.query('F', '^SINGLE-RUN')
        timeout = 1.1 * (self.ictime + self.optime) # raises error if set_ic_time/set_op_time not called before
        log.info(f"Waiting for single run sync to finish for {timeout} ms...")
        time.sleep(timeout / 1000)  # converting timeout (microseconds) to seconds
        res = q.read(self, expect(re='^EOSR(HLT)?'), read_again=True)
        was_terminated_by_ext_halt_condition = res.reply.groups()[0]=="HLT" # EOSRHLT
        return was_terminated_by_ext_halt_condition
    
    def set_ic_time(self, ictime):
        "Sets IC (initial condition) time in MILLISECONDS."
        ensure(ictime, within=range(0,999999)); self.ictime = ictime
        return self.query('C%06d' % ictime, expect(eq=f"T_IC={ictime}"))
    
    def set_op_time(self, optime):
        "Sets OP (operation mode) time in MILLISECONDS"
        ensure(optime, inrange=(0,999999)); self.optime = optime
        return self.query('c%06d' % optime, expect(eq=f"T_OP={optime}"))
    
    def get_data(self):
        "Supposed to be called when a read out group is defined and the machine is in (synchronous) OP mode. "
        q = self.query('l', "^No data!|.*$")
        if q.response == "No data!": return None
        data = []
        while True:
            resp = q.read(self, expect(re=r"^([-\d\.\s]*|EOD)*$"), read_again=True).response.split()
            if "EOD" in resp: break
            data.append(list(map(float, resp)))
    
    def read_element_by_address(self, address):
        "Read any machine element voltage. Expecting 16-bit element address as integer."
        ensure(address, isa=int)
        response_match = self.query("g%04X" % address, r"(?P<value>.+)\s+(?P<id>.+)").reply
        return response_match.groupdict() # return the dictionary value-> ..., id-> ..., caveat, should be all numeric!?
    
    def set_ro_group(self, addresses):
        "Defines a read out group, expects addresses to be an integer list of 16-bit element addresses."
        ensure(len(addresses), inrange=(1,5*self.MAX_RO_GROUP_SIZE))
        for a in addresses: ensure(a, isa=int)
        return self.query("G" + ";".join([f"{a:04X}" for a in addresses]) + ".")
    
    read_ro_group = command('f', expect(split=";", type=float), help="Query for currently set read out group")
    read_digital = command("R", expect(re="^"+r"\d\s"*(DIGITAL_INPUT_PORTS-1), split='', type=bool), help="Read digital inputs")
    
    def digital_output(self, port, state):
        "Set digital output pins of the Hybrid Controller"
        ensure(port, inrange=(0, self.DIGITAL_OUTPUT_PORTS)); ensure(state, isa=bool)
        return self.query(f"{'D' if state else 'd'}{port:1d}")

    def set_xbar(self, address, config):
        """Exactly {self.XBAR_CONFIG_BYTES*2} HEX-nibbles are required to config data."""
        ensure(address, isa=int); ensure(config, length=self.XBAR_CONFIG_BYTES*2)
        return self.query(f"X{address:04X}{config}", expect(eq="XBAR READY"))
    
    read_mpts = wont_implement("because it is just a high-level function which calls pot_set and iterates a list of potentiometer address/names.")
    
    def set_pt(self, address, number, value):
        "Set a digital potentiometer by address/number."
        ensure(value, inrange=(0,1))
        value = int(value * self.DPT_MAX_INT_VALUE) # 0000 <= value <= 1023
        return self.query(f"P{address:04X}{number:02X}{value:04d}", expect(eq=f"P{address:X}.{number:X}={value:d}"))
    
    def read_dpts(self):
        """
        Asks the Hybridcontroller for reading out *all* DPTs in the machine (also DPT24 modules).
        Returns mapping of PT module to list of values in that module.
        """
        response = self.query("q", '.*,.*').reply.string
        # strings looks like:
        #{'200': '204,44,244,102,89,244,89,639,639,1023,1023,1023,639,639,1023,1023,1023,300,300,79,79,102,102,300',
        # '300': '0,0,0,0,0,0,0,306'} 
        # numbers is the same parsed as hex/int.
        strings = dict([ items.split(":") for items in response.strip().split(";") ])
        numbers = { int(k, base=16): list(map(int, v.split(","))) for k,v in strings.items() }
        dpt2float = lambda val: val / 2**self.DPT_RESOLUTION
        floats = { k: list(map(dpt2float, v)) for k,v in numbers.items() }
        return floats
    
    def get_status(self):
        "Queries the HybridController about it's current status. Will return a dictionary."
        response = self.query('s', '.*=.*,.*').reply.string # expecting something with commas :-)
        state = dict(items.split("=") for items in response.split(","))
        state['RO-GROUP'] = state['RO-GROUP'].split(";")
        state['DPTADDR'] = state['DPTADDR'].split(";") # don't resolve mapping
        return state

    get_op_time = command('t', expect(re=r"t_OP=(?P<time>-?\d*)", ret='time', type=float), help="Asks about current OP time")
    reset = command('x', expect(eq='RESET'), help="Resets the HybridController (has no effect on python instance itself)")

#
#
#
# 
#

"""
Connection or "backends" for the PyHyCon.

The :class:`HyCon.HyCon` class requires a file handle to be passed. Usually, file
APIs are cursed in many languages (also python), but you can get your way out
with the following examples and also classes in this module.

Tested or "proven" connection interfaces are:

* :class:`tcpsocket`: A small adapter for the :meth:`socket` python builtin.
* :class:`human`: A small dummy adapter which prints to the interactive user terminal session
  and expects commands from there (the naming is ironically pointing to the human acting as
  actual Hybrid controller hardware endpoint).

Somewhat experimental but known to work is especially for unidirectional access:

* :class:`StringIO.StringIO`: Circumventing file access by reading from/to strings.
* :class:`sys.stdout` for just dumping HyCon-generated instructions

.. note:: All functions in this module do some progress reporting if you enable python logging
   facilities. Do so with
   
   >>> import logging
   >>> logging.basicConfig(level=logging.INFO)

Usage Examples
--------------

The following examples are suitable to be run in an interactive python REPL to explore
the APIs.

Using PyHyCon with a microcontroller "simulator"
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

>>> from hycon import HyCon
>>> ac = HyCon(human())                                                                                            
>>> ac.set_ic_time(1234)                                  # doctest: +SKIP                                                                                  
<< Sending [C001234] to uC
[type reply of uC]>> T_IC=1234
HyConRequest(C001234, expect(eq: T_IC=1234), self.executed=True, response=T_IC=1234, reply=T_IC=1234)

Using PyHyCon only for writing firmware commands

>>> import hycon, sys
>>> ac = hycon.HyCon(sys.stdout, unidirectional=True)
>>> ac.set_ic_time(234)
C000234HyConRequest(C000234, expect(eq: T_IC=234), self.executed=True, response=n.a., reply=n.a.)

Such a unidirectional approach can be interesting when generating bitstreams, for
larger integration tests, etc.

Using PyHyCon over TCP/IP
^^^^^^^^^^^^^^^^^^^^^^^^^

>>> sock = tcpsocket("localhost", 12345)                 # doctest: +SKIP
>>> ac = HyCon(sock)                                     # doctest: +SKIP
>>> ac.reset()                                           # doctest: +SKIP
>>> ac.digital_output(3, True)                           # doctest: +SKIP
>>> ac.set_op_time(123)                                  # doctest: +SKIP
>>> ac.set_xbar(0x0040, "0000000210840000781B")          # doctest: +SKIP

This setup is particularly interesting when connecting network-transparently to
actual hardware. The target TCP server is expected to route the contents to a
serial port/USB UART without introducing buffering. Examples for this kind of
stub servers are given at networking-hc_.

Using PyHyCon over Serial

>>> fh = serial("/dev/ttyUSB0", 115200)                     # doctest: +SKIP
>>> ac = HyCon(fh)                                          # doctest: +SKIP
>>> ac.digital_output(3, True)                              # doctest: +SKIP
>>> # etc.

You are encouraged to use the :class:`serial` class, which uses
`PySerial <https://pythonhosted.org/pyserial/>`_ under the hood and does the
clearing/resetting of the stream for you (something which is more cumbersome over serial
then over TCP).

If you really want, you can also use PySerial directly:

>>> import Serial from serial                               # doctest: +SKIP
>>> fh = Serial("/dev/ttyUSB0", 115200)                     # doctest: +SKIP

Note that this approach suffers from binary/string conversions, but you could
probably wrap ``open(fh)`` in some text mode.

If you (also) do not like PySerial, you can connect to a char device on a unixoid
operating system with vanilla python (this example is kind-of-untested):

>>> import os
>>> fd = os.open("/dev/ttyUSB0", os.O_RDWR | os.O_NOCTTY)      # doctest: +SKIP
>>> fh = os.fdopen(self.fd, "wb+", buffering=0)                # doctest: +SKIP
>>> ac = HyCon(fh)                                             # doctest: +SKIP

In this case, you certainly want to set the connection parameters (baud rate, etc.) by
``ioctl``, for instance in before on your linux terminal using a command
like ``stty -F /dev/ttyUSB0 115200``, or with ``stty ospeed 115200``
and ``stty ispeed 115200`` on Mac. Furthermore, when using this approach, consider
writing a small wrapper which runs ``fh.flush()`` after writing.
"""

import logging, time
log = logging.getLogger('connections')

def repeated_reset(fh):
    """
    This routine tries to *clear* output buffers of the hycon UART by sending repeated
    reset instructions and waiting until the reply "RESET" appears on the line. Doing this,
    it implements the HyCon protocol, but the :mod:`HyCon` code does not deal with
    connection issues, which is why this function is aprt of :mod:`connections`.
    
    Calling this function on beginning the setup is recommended for direct serial
    connections.
    
    You can also call to this method with the :fun:`HyCon.HyCon.repeated_reset()` shorthand.
    
    This function returns ``True`` when the connection suceeded, else ``False``.
    """
    max_reset_attempts = 10
    for i in range(max_reset_attempts):
        log.info(f"Attempt {i}/{max_reset_attempts} to reset the controller/connection...")
        fh.write("x")               # HyCon protocol reset instruction
        if max_reset_attempts > 0:  # save time: If connection is stable, do not sleep...
            time.sleep(1)           # ...otherwise sleep one second to allow reset
        res = fh.readline()
        if res == "RESET\n":
    log.warn("Could not properly reset the controller!")
    return False

class human:
    "Dummy IOWrapper for testing HyCon.py without the actual hardware"
    def write(self, sth):      print(f"<< Sending [{sth}] to uC")
    def readline(self): return input("[type reply of uC]>> ")

class tcpsocket:
    "Wrapper for communicating with HyCon over TCP/IP. See also HyCon-over-TCP.README for further instructions"
    def __init__(self, host, port):
        from socket import socket # builtin
        self.s = socket()
        self.s.connect((host,port))
        self.fh = self.s.makefile(mode="rw", encoding="utf-8")
        log.info(f"Connected to TCP {host}:{port}")
        
        repeated_reset(self)
    def write(self, sth):
        "Expects sth to be a string"
        self.s.sendall(sth.encode("ascii"))
    def readline(self):
        # instead also: self.s.recv(123)
        return self.fh.readline()
    
class serial:
    "Small wrapper for making the use of PySerial more handy (no need for extra import)"
    def __init__(self, port, baudrate, **passed_options):
        try:
            from serial import Serial  # requires pyserial
        except ImportError:
            raise ImportError("Please install PySerial in order to use it.")
        # in the following, some arguments are added for debugging...
        self.s = Serial(port, baudrate, 
            timeout=0.2, # in seconds            
            **passed_options)
        log.info(f"Connected to serial port {port} with {baudrate} baud")

        repeated_reset(self) # this is crucial for direct serial connections
        
        self.s.write(sth.encode("ascii"))
        return self.s.readline().decode("ascii")
            
"""
The autosetup module of the hycon package is the python implementation of the
similar named feature of the Perl HyCon library.

It is used to setup a hybrid controller based from a YAML file which includes
a mapping from names to computing element and potentiometer addresses and
a problem description containing information about timing, potentiometer
values (coefficients) and a readout group of interest. It can also describe
the configuration of an XBAR module.

The idea of this YAML file is to describe the analog circuit as complete as
possible, to keep the steering hycon code in perl (or python, respectively)
short. Furthermore, it brings some kind of highlevel description of the
circuit, since many parts of the circuit are given names.

This idea is some intermediate idea to the pyFPAA code which I wrote. It can
be seen as an alternative high-level frontend to pyHyCon. Remember, the idiom
of pyHyCon is to provide only a lowest level API for interaction with the
hardware hybrid controller.

About the history of this code: Bernd started to write his auto-setup code
at 25-DEC-2019. I started to write my pyFPAA code at the same time. Roughly
one year later, where most of the time was spent at other stuff, I now port
parts of Bernds auto-setup code to python in order to be able to use the
same YAML files.
"""

from .HyCon import HyCon

class DotDict(dict):
    """
    Small syntactic sugar: Dot notation to access to dictionary attributes, which is especially
    handy for deeply nested dicts. There are plenty of similar library for python around, but
    this implementation is only five lines (yes, five). The following usage example is longer
    then the implementation:
    
    >>> a = { "b": 42, "non-identifier": 3, "foo": { "bar": { 3: 123 } }}
    >>> a = DotDict(a)
    >>> a.b
    42
    >>> a.foo.bar
    {3: 123}
    >>> a["non-identifier"]   # traditional __getitem__ access is still possible
    3
    >>> a.foo.bar[3]          # especially hand for non-pure-ascii identifiers
    123
    >>> DotDict(DotDict(a)).b # DotDict can be applied repeatedly without loss of functionality
    42
    >>> c = DotDict()
    >>> c.foo = "b"
    >>> c                     # also works for setting, not only reading
    {'foo': 'b'}
    >>> c.bar = DotDict()
    >>> c.bar.baz = "bla"     # Limitation for nested setting: create nested DotDicts first.
    """
    def __getattr__(*args):
        val = dict.get(*args)
        return DotDict(val) if type(val) is dict else val
    __setattr__ = dict.__setitem__
    __delattr__ = dict.__delitem__ 

# This is tolen from ../fpaa/fpaa.py  . TODO: just import or so.
    
class PotentiometerAddress(namedtuple("PotentiometerAddress", ["address", "number"])):
    """
    Stores a potentiometer address, which is a tuple of a (typically hex-given) bus address
    of the hardware element and an element-internal number. Example:
    
    >>> a = PotentiometerAddress(0x200, 0x20)
    >>> b = PotentiometerAddress.fromText("0x200/20")  # FIXME: Is number really base 16?
    >>> a == b
    True
    >>> a.address   # Don't forget that python standard numeric output is in decimal
    512
    >>> a.toText()
    '0x200/20'
    """
    
    def fromText(cls, text):
        "Parses something like 0x200/2 to (0x200, 2). Will also accept 0200/2 as hex."
        if not cls.isPotentiometerAddress(text):
            raise ValueError("'%s' doesn't look like a valid potentiometer address. Should be like 0x200/2 or 0200/2", text)
        address,number = text.split("/")
        return cls(int(address,16), int(number,16))
    def isPotentiometerAddress(cls, text):
        # should rather have a regexp or so...
        return isinstance(text,str) and text.count("/") == 1 # or not "0x" in text:
    def toText(self):
        return "0x%x/%x" % (self.address,self.number)

def autosetup(hycon, conf, reset=True):
    """
    hycon is expected to be an instance of HyCon.
    conf is expected to be a dictionary.
    
    If you want to load from a YAML file, use the yaml_load function.
    
    TODO: XBAR support not yet implemented.
    """
    
    if not 'problem' in conf: raise ValueError("No problem section defined!")
    problem = DotDict(conf["problem"]) # syntactic sugar
    elements = DotDict(conf["elements"]) # syntactic sugar
    
    if reset: hycon.reset()
    
    if "times" in problem:
        # Note: All these times are in MILLISECONDS
        if "ic" in problem.times: 
            hycon.set_ic_time(problem.times.ic)
        if "op" in problem.times:
            hycon.set_op_time(problem.times.op)

    # Initial Conditions
    # TODO: Skipping here, because mainly used in XBAR-relevant code...
    #for element in problem.get("IC",[]): # aka is dict
    #    value = problem.IC[element]
    #    sign = (value < 0)
    #    number = ...

    try:
        # Set potentiometer values (coefficients):
        for name, value in problem.get("coefficients", {}).items():
            dp = PotentiometerAddress.fromText(elements[name])
            hycon.set_pt(dp.address, dp.number, value)

        # Define read out group if specified:
        if "ro-group" in problem:
            # Remark: This will fail with weird "None" errors for elements in ro-group which are not defined.
            addresses = list(map(elements.get, problem["ro-group"]))
            hycon.set_ro_group(addresses)
    except KeyError as e:
        raise KeyError("Unknown coefficient or computing element: '%s'. It is not part of the element map." % e)
        
    # Derive the required XBAR setup:
    #if (defined($self->{problem}) and defined($xbar_address)) {
    # TODO: Skipping here...

def autoconnect(conf):
    """
    Opens a file handle to the target position found in the YAML file. Follows the same
    rules as the perl routine, i.e. looks for *serial* or *tcp* key and connects
    according to the parameters.
    
    Example serial port configuration:
    
    ::
    
        serial:
          port: /dev/cu.usbserial-DN050L21
          bits: 8
          baud: 250000
          parity: none
          stopbits: 1
          poll_interval: 10
          poll_attempts: 20000
          
    Note that we only take into account port and baud rate, since everything else looks standard
    and the pySerial port cannot deal with an integer stopbit ``1`` but expects something like ``serial.STOPBITS_ONE``.
    As a note to the future, https://tools.ietf.org/html/rfc2217.html is supported by pySerial and
    should be adopted by the YAML definition.,
    
    Example TCP port configuration:
    
    ::
    
        tcp:
            addr: 192.168.31.190
            port: 12345
            connection_timeout: 2
            timeout: 0.1
            poll_interval: 10
            poll_attempts: 2000
            quick_start: False
            
    Again, we only take into account the IP address and the TCP port, everything else is ignored for the
    time being.
    """
    from .connections import tcpsocket, serial
    
    c = DotDict(conf)
    if 'tcp' in c:
        return tcpsocket(c.tcp.addr, c.tcp.port)
    elif 'serial' in c:
        return serial(c.serial.port, c.serial.baud)
    else:
        raise ValueError("Could neither find `tcp` nor `serial` section in configuration. Available keys: {list(c)}")
    
class AutoConfHyCon(HyCon):
    """
    Syntactic sugar to provide a "setup" method similar to the perl
    HyCon API.
    
    TODO: Should also provide other methods for high level value read
    and set access.
    
    ``conf`` can either be a string holding the YAML filename  or a dictionary (holding the
    configuration content, i.e. parsed YAML file).
    """
    def __init__(self, conf):
        # Will store conf as DotDict for easier later access
        self.conf = DotDict(yaml_load(conf) if isinstance(conf, str) else conf)
        self.fh = autoconnect(self.conf)
        self.unidirectional = False # could be steered by conf, too
        self.autosetup()
        
    def autosetup(self, conf=None, reset=False):
        if not conf: conf = self.conf
        autosetup(self, conf, reset)
        
    def get_data_by_name(self):
        "Get readout group data handily labeled by name"
        if not "ro-group" in self.conf.problem: raise ValueError("No Read-out group defined")
        data = self.get_data() # shape (sample_points, readout_group_length)
        if not data: return None
        return OrderedDict( (name, [line[i] for line in data]) for i,name in enumerate(self.conf.problem["ro-group"]) )
    
    def set_pt_by_name(self, name, value):
        "Set a digital potentiometer by name"
        dp = PotentiometerAddress.fromText(self.conf.elements[name])
        return self.set_pt(dp.address, dp.number, value)

    def read_element_by_name(self, name):
        "Reads element by name"
        if not name in self.conf.elements:
            raise ValueError(f"{name} not defined in yaml element dict")
        address = self.conf.elements[name]
        # address is already an int, correctly parsed...
        return self.read_element_by_address(address)

    def read_dpts_by_name(self):
        """
        Returns single map of DPT name to value (as float).
        """
        floats = self.read_dpts()
        interesting_potentiometers = \
          { PotentiometerAddress.fromText(addr): name for name,addr in self.conf.elements.items() \
            if PotentiometerAddress.isPotentiometerAddress(addr) }
        named_values = {}
        for p, name in interesting_potentiometers.items():
            if p.address in floats and p.number < len(floats[p.address]):
                named_values[name] = floats[p.address][p.number]
        return named_values
    
    def read_ro_group_by_name(self):
        "Returns an OrderedDict of read-out group elements, with names"
        req = self.read_ro_group()
        return OrderedDict(zip(self.conf.problem["ro-group"], map(float,req.response.split(";"))))

#
#
#
# 
#

"""
A HyCon command stream interpreter.

Will spill out LISP-like commands which can be fed into the hycon
again. This allows for replaying, which is helpful for a number of 
special scenarios such as:

* Offline-validating a HyCon instruction stream
* Man-in-the-middle inspecting an HyCon instruction stream
* Validating the correctness of high-level HyCon instructions
  (such as emitted by the PyFPAA or autosetup codes)

The code basically implements a character-by-character tokenizer/parser.
It is built based on a simple mapping datastructure which assigns each
one-letter command the respective PyHyCon method name. Furthermore,
method arguments can be read and converted.

It would be nice to join HyCon.py and replay.py to a single file which
translates between the serial protocol and the OOP API calls. The transformation
is quite trivial, but now we have a lot of code doing nothing of 
bigger interest.

The ordering follows the AVR Ino code.
"""

import sys, io, functools

def delayed(static_method):
    """
    This is a decorator for a static method in a class. It will
    make the function body "delayed", i.e. when calling the function,
    a future/promise/delay/deferred element is returned.
    Some parameter bounding (closure) happens: The decorated arguments
    are evaluated early while the later execution expects only a single
    reader argument. Example to follow the logic:
    
    >>> f = lambda a,b,c,d: print(a,b,c,d)
    >>> g = delayed(f)
    >>> h = g(1,2,3)
    >>> h(4)
    4 1 2 3
    """
    @functools.wraps(static_method)
    def decorated(*w,**kw):
        def deferred(reader):
            return static_method(reader, *w, **kw)
        #deferred.__name__ = f"{static_method.__name__}({str(w)[1:-1]}{',' if kw else ''}{str(kw)[1:-1]})"
        return deferred
    return decorated

class consume:
    """
    Consume is an ugly namespace and not a class, actually.
    The basic idea of these (static!) functions is to be
    called delayedly with a function as it's argument which
    acts like the io.IOBase.read() function, i.e. advances
    an internal cursor (side effect) and returns *n* characters
    from the stream. Crude Example:
    
    >>> tokenizer = [consume.exact("test"), consume.decimals(3), consume.exact("foo"), consume.hex(2)]
    >>> test = io.StringIO("test123fooAA")
    >>> reader = test.read
    >>> [ token(reader) for token in tokenizer ]
    ['test', 123, 'foo', 170]
    """
    @delayed
    def exact(reader, text):
        proof = reader(len(text))
        if proof != text: raise ValueError(f"Expected '{text}' but got '{proof}'.")
        return text
    @delayed
    def number(reader, digits, base, multiply=1):
        """
        Reads a number with #digit digits in some base. Can perform a multiplication afterwards.
        
        >>> consume.number(8,16)(io.StringIO("deadbeef").read)
        3735928559
        >>> consume.number(2,10,multiply=2)(io.StringIO("42").read)
        84
        """
        text = reader(digits)
        if not text: raise ValueError(f"Encountered end of input but expected {digits} digits in base {base}")
        try: return int(text,base) * multiply
        except ValueError as e: raise ValueError(f"Could not cast '{text}' as a base-{base} number.")
    @delayed
    def decimals(reader, number_of_digits):
        return consume.number(number_of_digits, 10)(reader)
    @delayed
    def hex(reader, number_of_digits):
        return consume.number(number_of_digits, 16)(reader)
    @delayed
    def list(reader, split, end, digits, base):
        """
        Reads a list of numbers. Limitations:
        * Always expects end token to come
        * All numbers must have same number of digits (and same base)
        * Cannot handle empty list or and won't accept end-of-file before end token.
        
        Examples:
        
        >>> consume.list(split=",",digits=1,base=10,end=".")(io.StringIO("1,5,2,3,9.").read)
        [1, 5, 2, 3, 9]
        >>> consume.list(split=":",digits=2,base=16,end=";")(io.StringIO("5a:88:ff:ff;").read)
        [90, 136, 255, 255]
        """
        assert len(split) == len(end)
        retlst = []
            retlst.append( consume.number(digits, base)(reader) )
            #print(f"List: {retlst}")
            sep = reader(len(split))
            #print(f"Sep: '{sep}'")
            if sep == end: break
            elif sep != split: raise ValueError(f"Unexpected character(s) or end of stream: '{sep}'")
        return retlst

not_implemented = lambda comment: ("NOT_IMPLEMENTED", comment)

mapping = {
    "a": "disable_ovl_halt",
    "A": "enable_ovl_halt",
    "b": "disable_ext_halt",
    "B": "enable_ext_halt",
    "c": ("set_op_time", consume.decimals(6)),
    "C": ("set_ic_time", consume.decimals(6)),
    "d": ("digital_output", consume.decimals(1), False),
    "D": ("digital_output", consume.decimals(1), True),
    "e": "repetitive_run",
    "E": "single_run",
    "F": "single_run_sync",
    "f": "read_ro_group",
    "G": ("set_ro_group", consume.list(split=";", digits=4, base=16, end=".")),
    "g": ("read_element_by_address", consume.hex(4)),
    "h": "halt",
    "i": "ic",
    "L": not_implemented("Locate a computing element"),
    "l": "get_data",
    "o": "op",
    "P": ("set_pt", consume.hex(4), consume.hex(2), consume.number(digits=4, base=10, multiply=1./HyCon.DPT_MAX_INT_VALUE)),
    "q": "read_dpts", # is not_implemented()
    "R": "read_digital",
    "s": "get_status",
    "S": "pot_set",
    "t": "get_op_time",
    "x": "reset",
    "X": ("set_xbar", consume.hex(4), consume.hex(HyCon.XBAR_CONFIG_BYTES*2)),
    "?": not_implemented("Prints help")
}

class HyConRequestReader:
    """
    Converts HyCon "configuration strings" to high level API calls.
    This can be seen as the inverse operation to calling the HyCon.
    
    Instances of this class act as iterator. It will consume the
    incoming stream character by character (or the whole string, if
    given as a string).
    
    Example:
    
    >>> instructions = 'C000100c015000P0200000204P0300030000G0362;0363;0220;0221;0222;0223.'
    >>> commands = list(HyConRequestReader(instructions))
    >>> for c in commands: print(c)
    ('set_ic_time', 100)
    ('set_op_time', 15000)
    ('set_pt', 512, 0, 0.19941348973607037)
    ('set_pt', 768, 3, 0.0)
    ('set_ro_group', [866, 867, 544, 545, 546, 547])
    >>> replayed = io.StringIO()
    >>> hc = HyCon(replayed, unidirectional=True)
    >>> replay(hc, commands)
    >>> replayed.getvalue() == instructions
    True
    """
    def __init__(self, stream_or_string, mapping=mapping):
        self.subject = io.StringIO(stream_or_string) if isinstance(stream_or_string, str) else stream_or_string
        self.shadow = io.StringIO()
        self.mapping = mapping
    
    def read(self, size=-1):
        ret = self.subject.read(size)
        self.shadow.write(ret)
        return ret
    
    def debugline(self):
        read_in = self.shadow.getvalue()
        return f"Encountered after having read {len(read_in)} characters '{read_in}'"
        
    def read_arguments(self, rhs):
        if isinstance(rhs, tuple):
            return tuple(map(self.read_arguments, rhs))
        if callable(rhs):
            res = rhs(self.read)
            return res
        return rhs

    def read_command(self):
        command = self.read(1)
        if not command:
            raise StopIteration # end of input
        if command not in self.mapping:
            raise ValueError(f"Command '{command}' not part of the valid HyCon command characters. {self.debugline()}")
        return self.read_arguments(self.mapping[command])

    # act as an iterator
    __iter__ = identity
    __next__ = read_command
    
def replay(hycon, commands):
    """
    Given commands a list of tuples, this will mostly act like operator.methodcaller
    on them. If no arguments are given, the tuple can be omitted.
    
    Basic example:
    
    >>> class WannaBeHyCon:
    ...     def toot(self,x): print("too(%s)ooot" % x)
    ...     def bar(self): print("this is bar")
    ...     def buz(self,a,b,c): print(f"a*b = {a*b} but what is {c}")
    >>> replay(WannaBeHyCon(), [ "bar", ("toot", "fuz"), ("buz", 1,2,3) ])
    this is bar
    too(fuz)ooot
    a*b = 2 but what is 3

    The command format is produced by the HyConRequestReader and thus can be fed
    into this replay function:
    
    >>> replay(HyCon(sys.stdout, unidirectional=True), HyConRequestReader("xiohaARt"))
    xiohaARt
    
    This works for almost any useful instruction stream.
    """
    for command in commands:
        if not isinstance(command, tuple): command = (command,)
        method,*args = command        # LISP-like unpacking (command, *args)
        getattr(hycon, method)(*args) # calling