add forward autodiff

lib/gpt/ad/__init__.py

@@ -17,5 +17,4 @@
 #    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 #
 import gpt.ad.reverse
-
-# import gpt.ad.forward
+import gpt.ad.forward

lib/gpt/ad/forward/__init__.py

@@ -0,0 +1,22 @@
+#
+#    GPT - Grid Python Toolkit
+#    Copyright (C) 2023  Christoph Lehner ([email protected], https://github.com/lehner/gpt)
+#
+#    This program is free software; you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation; either version 2 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License along
+#    with this program; if not, write to the Free Software Foundation, Inc.,
+#    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+from gpt.ad.forward.infinitesimal import infinitesimal
+from gpt.ad.forward.landau import landau
+from gpt.ad.forward.series import series
+from gpt.ad.forward.transform import norm2, inner_product, cshift

lib/gpt/ad/forward/infinitesimal.py

@@ -0,0 +1,73 @@
+#
+#    GPT - Grid Python Toolkit
+#    Copyright (C) 2023  Christoph Lehner ([email protected], https://github.com/lehner/gpt)
+#
+#    This program is free software; you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation; either version 2 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License along
+#    with this program; if not, write to the Free Software Foundation, Inc.,
+#    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+
+
+class infinitesimal:
+    def __init__(self, value):
+        if isinstance(value, str):
+            value = {value: 1}
+        self.value = value
+
+    def __pow__(self, n):
+        value = {}
+        for v in self.value:
+            value[v] = self.value[v] * n
+        return infinitesimal(value)
+
+    def __mul__(self, other):
+        value = {}
+        for v1 in self.value:
+            value[v1] = self.value[v1]
+        for v2 in other.value:
+            if v2 in value:
+                value[v2] += other.value[v2]
+            else:
+                value[v2] = other.value[v2]
+        return infinitesimal(value)
+
+    def __str__(self):
+        r = ""
+        for v in sorted(self.value):
+            if r != "":
+                r = r + "*"
+            if self.value[v] == 1:
+                r = r + v
+            else:
+                r = r + f"{v}**{self.value[v]}"
+        return r
+
+    def __hash__(self):
+        return hash(self.__str__())
+
+    def __eq__(self, other):
+        return self.__str__() == other.__str__()
+
+    def __cmp__(self, other):
+        return self.__str__().__cmp__(other.__str__())
+
+    def symbols(self):
+        return tuple(sorted(list(self.value.keys())))
+
+    def behaves_as(self, other):
+        for s in other.value:
+            n0 = self.value[s] if s in self.value else 0
+            n1 = other.value[s]
+            if n0 < n1:
+                return False
+        return True

lib/gpt/ad/forward/landau.py

@@ -0,0 +1,55 @@
+#
+#    GPT - Grid Python Toolkit
+#    Copyright (C) 2023  Christoph Lehner ([email protected], https://github.com/lehner/gpt)
+#
+#    This program is free software; you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation; either version 2 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License along
+#    with this program; if not, write to the Free Software Foundation, Inc.,
+#    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+
+
+class landau:
+    def __init__(self, *infinitesimals):
+        self.infinitesimals = infinitesimals
+
+    def accept(self, i):
+        for j in self.infinitesimals:
+            if i.behaves_as(j):
+                return False
+        return True
+
+    def __add__(self, other):
+        if self is other:
+            return self
+        infinitesimals = []
+        for i in self.infinitesimals + other.infinitesimals:
+            keep = True
+            for n, j in enumerate(infinitesimals):
+                if i.behaves_as(j):
+                    keep = False
+                elif j.behaves_as(i):
+                    infinitesimals[n] = i
+            if keep:
+                infinitesimals.append(i)
+        infinitesimals = list(set(infinitesimals))
+        return landau(*infinitesimals)
+
+    def __str__(self):
+        a = []
+        for i in self.infinitesimals:
+            a.append(str(i))
+        r = ",".join(sorted(a))
+        return f"O({r})"
+
+    def __eq__(self, other):
+        return str(self) == str(other)

lib/gpt/ad/forward/series.py

@@ -0,0 +1,168 @@
+#
+#    GPT - Grid Python Toolkit
+#    Copyright (C) 2023  Christoph Lehner ([email protected], https://github.com/lehner/gpt)
+#
+#    This program is free software; you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation; either version 2 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License along
+#    with this program; if not, write to the Free Software Foundation, Inc.,
+#    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+import gpt as g
+from gpt.ad.forward import infinitesimal
+
+
+def promote(other, landau_O):
+    if isinstance(other, infinitesimal):
+        other = series({other: 1}, landau_O)
+    elif g.util.is_num(other):
+        other = series({infinitesimal({}): other}, landau_O)
+    return other
+
+
+class series:
+    def __init__(self, terms, landau_O):
+        self.landau_O = landau_O
+        if not isinstance(terms, dict):
+            i0 = infinitesimal({})
+            terms = {i0: terms}
+        self.terms = terms
+
+    def __str__(self):
+        r = ""
+        for t in self.terms:
+            if r != "":
+                r = r + " + "
+            r = r + "(" + str(self.terms[t]) + ")"
+            si = str(t)
+            if si != "":
+                r = r + "*" + si
+        return r
+
+    def distribute2(self, other, functional):
+        other = promote(other, self.landau_O)
+        # first merge landau_Os
+        landau_O = self.landau_O + other.landau_O
+        # then merge terms
+        terms = {}
+        for t1 in self.terms:
+            for t2 in other.terms:
+                i = t1 * t2
+                if not landau_O.accept(i):
+                    continue
+                if i not in terms:
+                    terms[i] = g(functional(self.terms[t1], other.terms[t2]))
+                else:
+                    terms[i] += functional(self.terms[t1], other.terms[t2])
+        return series(terms, landau_O)
+
+    def distribute1(self, functional):
+        # then merge terms
+        terms = {}
+        for t1 in self.terms:
+            if t1 not in terms:
+                terms[t1] = g(functional(self.terms[t1]))
+            else:
+                terms[t1] += functional(self.terms[t1])
+        return series(terms, self.landau_O)
+
+    def function(self, functional):
+        root = self[1]
+        # get nilpotent power
+        nilpotent = self - root
+        maxn = 0
+        i0 = infinitesimal({})
+        for t in nilpotent.terms:
+            if t == i0:
+                continue
+            n = 1
+            tn = t
+            while self.landau_O.accept(tn):
+                tn = tn * t
+                n += 1
+            maxn = max([maxn, n])
+        res = series({i0: functional(root, 0)}, self.landau_O)
+        delta = nilpotent
+        nfac = 1.0
+        for i in range(1, maxn):
+            nfac *= i
+            res += delta * functional(root, i) / nfac
+            if i != maxn - 1:
+                delta = delta * nilpotent
+        return res
+
+    def __iadd__(self, other):
+        res = self + other
+        self.landau_O = res.landau_O
+        self.terms = res.terms
+        return self
+
+    def __mul__(self, other):
+        return self.distribute2(other, lambda a, b: a * b)
+
+    def __rmul__(self, other):
+        if g.util.is_num(other):
+            return self.__mul__(other)
+        raise Exception("Not implemented")
+
+    def __add__(self, other):
+        other = promote(other, self.landau_O)
+        # first merge landau_Os
+        landau_O = self.landau_O + other.landau_O
+        # then merge terms
+        terms = {}
+        for t1 in self.terms:
+            if not landau_O.accept(t1):
+                continue
+            terms[t1] = self.terms[t1]
+        for t2 in other.terms:
+            if not landau_O.accept(t2):
+                continue
+            if t2 not in terms:
+                terms[t2] = other.terms[t2]
+            else:
+                terms[t2] = g(terms[t2] + other.terms[t2])
+        return series(terms, landau_O)
+
+    def __sub__(self, other):
+        other = promote(other, self.landau_O)
+        # first merge landau_Os
+        landau_O = self.landau_O + other.landau_O
+        # then merge terms
+        terms = {}
+        for t1 in self.terms:
+            if not landau_O.accept(t1):
+                continue
+            terms[t1] = self.terms[t1]
+        for t2 in other.terms:
+            if not landau_O.accept(t2):
+                continue
+            if t2 not in terms:
+                terms[t2] = other.terms[t2]
+            else:
+                terms[t2] = g(terms[t2] - other.terms[t2])
+        return series(terms, landau_O)
+
+    def __truediv__(self, other):
+        return (1.0 / other) * self
+
+    def __radd__(self, other):
+        return self.__add__(other, self)
+
+    def __getitem__(self, tag):
+        if tag == 1:
+            tag = infinitesimal({})
+        return self.terms[tag]
+
+    def __setitem__(self, tag, value):
+        if tag == 1:
+            tag = infinitesimal({})
+        self.terms[tag] = value

lib/gpt/ad/forward/transform.py

@@ -0,0 +1,31 @@
+#
+#    GPT - Grid Python Toolkit
+#    Copyright (C) 2023  Christoph Lehner ([email protected], https://github.com/lehner/gpt)
+#
+#    This program is free software; you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation; either version 2 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License along
+#    with this program; if not, write to the Free Software Foundation, Inc.,
+#    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+import gpt as g
+
+
+def inner_product(sx, sy):
+    return sx.distribute2(sy, lambda a, b: g.inner_product(a, b))
+
+
+def norm2(sx):
+    return inner_product(sx, sx)
+
+
+def cshift(sx, mu, disp):
+    return sx.distribute1(lambda a: g.cshift(a, mu, disp))