Cactus Code Thorn Forms
Author(s)    : Erik Schnetter <schnetter@gmail.com>
Maintainer(s): Erik Schnetter <schnetter@gmail.com>
Licence      : LGPL
--------------------------------------------------------------------------
1. Purpose
Discrete Differential Forms

# Configuration definitions for thorn Forms
REQUIRES googletest

# Interface definition for thorn Forms
IMPLEMENTS: Forms

ActiveThorns = "Forms"
Forms::test_forms = yes

# Parameter definitions for thorn Forms
BOOLEAN test_forms "Test Forms at startup"
{
} "no"

# Schedule definitions for thorn Forms
if (test_forms) {
  SCHEDULE TestForms_Test AT basegrid
  {
    LANG: C
    OPTIONS: meta
  } "Test Forms"
}

#include "forms.hxx"
namespace Forms {
void q() {
  form<double, 3, 0> a0({0});
  form<double, 3, 1> a1({1, 2, 3});
  form<double, 3, 2> a2({4, 5, 6});
  form<double, 3, 3> a3({7});
  +a0;
  +a1;
  +a2;
  +a3;
  -a0;
  -a1;
  -a2;
  -a3;
  a0 + a0;
  a1 + a1;
  a2 + a2;
  a3 + a3;
  a0 - a0;
  a1 - a1;
  a2 - a2;
  a3 - a3;
  2 * a0;
  2 * a1;
  2 * a2;
  2 * a3;
  a0 * 2;
  a1 * 2;
  a2 * 2;
  a3 * 2;
  a0 / 2;
  a1 / 2;
  a2 / 2;
  a3 / 2;
  wedge(a0, a0);
  wedge(a0, a1);
  wedge(a0, a2);
  wedge(a0, a3);
  wedge(a1, a0);
  wedge(a1, a1);
  wedge(a1, a2);
  wedge(a2, a0);
  wedge(a2, a1);
  wedge(a3, a0);
}
} // namespace Forms

#include <cxxabi.h>
#include <array>
#include <cassert>
#include <cstddef>
#include <initializer_list>
#include <memory>
#include <ostream>
#include <string>
#include <type_traits>
namespace Forms {
using namespace std;
constexpr int bitsign(bool s) { return s ? -1 : +1; }
constexpr int bitsign(int s) { return bitsign(bool(s & 1)); }
constexpr size_t count_bits(size_t x) {
  size_t n = 0;
  while (x > 0) {
    n += x & 1;
    x >>= 1;
  }
  return n;
}
// Taken from
// <https://stackoverflow.com/questions/9640256/define-multiple-methods-with-parameters-from-variadic-templates>
inline string demangle(const char *mangled) {
  int status;
  unique_ptr<char[], void (*)(void *)> result(
      abi::__cxa_demangle(mangled, 0, 0, &status), std::free);
  return result.get() ? string(result.get()) : "<unknown type>";
}
template <typename T, size_t D, size_t R> struct form;
template <typename T, size_t D, size_t R> struct form {
  static_assert(D >= 0, "");
  static_assert(R >= 0 && R <= D, "");
  constexpr static size_t ncomponents() {
    size_t ncomps = 0;
    for (size_t idx = 0; idx < (1 << D); ++idx)
      ncomps += count_bits(idx) == R;
    return ncomps;
  }
  constexpr static size_t idx2comp(const size_t idx0) {
    assert(count_bits(idx0) == R);
    size_t comp = 0;
    for (size_t idx = 0; idx < (1 << D); ++idx) {
      if (idx == idx0)
        break;
      comp += count_bits(idx) == R;
    }
    return comp;
  }
  constexpr static size_t comp2idx(const size_t comp0) {
    size_t comp = 0;
    size_t idx;
    for (idx = 0; idx < (1 << D); ++idx) {
      if (count_bits(idx) == R) {
        if (comp == comp0)
          break;
        ++comp;
      }
    }
    assert(count_bits(idx) == R);
    return idx;
  }
  static constexpr size_t N = ncomponents();
  array<T, N> comps;
  form() : comps() {}
  form(const array<T, N> &comps) : comps(comps) {}
  static constexpr form pure(const T &rc) {
    form r;
    for (size_t n = 0; n < N; ++n)
      r.comps[n] = rc;
    return r;
  }
  template <typename T1, size_t D1, size_t R1>
  friend constexpr form<T1, D1, R1> make_pure(const T1 &rc);
  template <typename T1, size_t D1, size_t R1, typename F, typename U>
  friend constexpr form<U, D1, R1> fmap(const F &f, const form<T1, D1, R1> &x);
  template <typename T1, size_t D1, size_t R1, typename F, typename U>
  friend constexpr form<U, D1, R1> fmap(const F &f, const form<T1, D1, R1> &x,
                                        const form<T1, D1, R1> &y);
  template <typename T1, size_t D1, size_t R1, typename F>
  friend constexpr T1 foldl(const F &f, const T1 &z, const form<T1, D1, R1> &x);
  template <typename T1, size_t D1, size_t R1, typename F, typename G,
            typename U>
  friend constexpr U foldmap(const F &f, const G &g, const U &z,
                             const form<T1, D1, R1> &x);
  template <typename T1, size_t D1, size_t R1, typename F, typename G,
            typename U>
  friend constexpr U foldmap(const F &f, const G &g, const U &z,
                             const form<T1, D1, R1> &x,
                             const form<T1, D1, R1> &y);
  friend constexpr form operator+(const form &x) {
    return fmap([](const T &xc) { return +xc; }, x);
  }
  friend constexpr form operator-(const form &x) {
    return fmap([](const T &xc) { return -xc; }, x);
  }
  friend constexpr form operator+(const form &x, const form &y) {
    return fmap([](const T &xc, const T &yc) { return xc + yc; }, x, y);
  }
  friend constexpr form operator-(const form &x, const form &y) {
    return fmap([](const T &xc, const T &yc) { return xc - yc; }, x, y);
  }
  friend constexpr form operator*(const T &a, const form &x) {
    return fmap([&](const T &xc) { return a * xc; }, x);
  }
  friend constexpr form operator*(const form &x, const T &a) {
    return fmap([&](const T &xc) { return xc * a; }, x);
  }
  friend constexpr form operator/(const form &x, const T &a) {
    return fmap([&](const T &xc) { return xc / a; }, x);
  }
  friend constexpr bool operator==(const form &x, const form &y) {
    return fmap([](const T &xc, const T &yc) { return xc == yc; }, x, y).all();
  }
  friend constexpr bool operator!=(const form &x, const form &y) {
    return !(x == y);
  }
  constexpr bool all() const {
    return foldl([](const bool x, const bool y) { return x && y; }, true,
                 *this);
  }
  friend ostream &operator<<(ostream &os, const form &x) {
    os << "form<" << demangle(typeid(T).name()) << "," << D << "," << R
       << ">{comps:[";
    for (size_t n = 0; n < N; ++n) {
      if (n > 0)
        os << ",";
      const size_t idx = comp2idx(n);
      os << "(" << n << "/" << idx << ")";
      for (size_t d = 0; d < D; ++d)
        os << int((idx & (1 << d)) != 0);
      os << ":" << x.comps[n];
    }
    os << "]}";
    return os;
  }
};
template <typename T, size_t D, size_t R>
constexpr form<T, D, R> make_pure(const T &rc) {
  constexpr size_t N = form<T, D, R>::N;
  form<T, D, R> r;
  for (size_t n = 0; n < N; ++n)
    r.comps[n] = rc;
  return r;
}
template <typename T, size_t D, size_t R, typename F,
          typename U = result_of_t<F(T)> >
constexpr form<U, D, R> fmap(const F &f, const form<T, D, R> &x) {
  constexpr size_t N = form<T, D, R>::N;
  form<U, D, R> r;
  for (size_t n = 0; n < N; ++n)
    r.comps[n] = f(x.comps[n]);
  return r;
}
template <typename T, size_t D, size_t R, typename F,
          typename U = result_of_t<F(T, T)> >
constexpr form<U, D, R> fmap(const F &f, const form<T, D, R> &x,
                             const form<T, D, R> &y) {
  constexpr size_t N = form<T, D, R>::N;
  form<U, D, R> r;
  for (size_t n = 0; n < N; ++n)
    r.comps[n] = f(x.comps[n], y.comps[n]);
  return r;
}
template <typename T, size_t D, size_t R, typename F>
constexpr T foldl(const F &f, const T &z, const form<T, D, R> &x) {
  constexpr size_t N = form<T, D, R>::N;
  T r(z);
  for (size_t n = 0; n < N; ++n)
    r = f(r, x.comps[n]);
  return r;
}
template <typename T, size_t D, size_t R, typename F, typename G,
          typename U = result_of_t<F(T)> >
constexpr U foldmap(const F &f, const G &g, const U &z,
                    const form<T, D, R> &x) {
  constexpr size_t N = form<T, D, R>::N;
  U r(z);
  for (size_t n = 0; n < N; ++n)
    r = g(r, f(x.comps[n]));
  return r;
}
template <typename T, size_t D, size_t R, typename F, typename G,
          typename U = result_of_t<F(T, T)> >
constexpr U foldmap(const F &f, const G &g, const U &z, const form<T, D, R> &x,
                    const form<T, D, R> &y) {
  constexpr size_t N = form<T, D, R>::N;
  U r(z);
  for (size_t n = 0; n < N; ++n)
    r = g(r, f(x.comps[n], y.comps[n]));
  return r;
}
template <typename T, size_t D, size_t Rx, size_t Ry, size_t R = Rx + Ry>
constexpr form<T, D, R> wedge(const form<T, D, Rx> &x,
                              const form<T, D, Ry> &y) {
  constexpr size_t Nx = form<T, D, Rx>::N;
  constexpr size_t Ny = form<T, D, Ry>::N;
  form<T, D, R> r;
  for (size_t ny = 0; ny < Ny; ++ny) {
    for (size_t nx = 0; nx < Nx; ++nx) {
      const auto idxx = form<T, D, Rx>::comp2idx(nx);
      const auto idxy = form<T, D, Ry>::comp2idx(ny);
      const bool skip = (idxx & idxy) != 0;
      if (!skip) {
        const size_t idx = idxx | idxy;
        const size_t n = form<T, D, R>::idx2comp(idx);
        bool parity = false;
        bool count = false;
        for (size_t d = 0; d < D; ++d) {
          parity ^= (idxx & (1 << d)) != 0 && count;
          count ^= (idxy & (1 << d)) != 0;
        }
        r.comps[n] += bitsign(parity) * x.comps[nx] * y.comps[ny];
      }
    }
  }
  return r;
}
} // namespace Forms

# Main make.code.defn file for thorn Forms
# Source files in this directory
SRCS = forms.cxx test.cxx
# Subdirectories containing source files
SUBDIRS =

#include "forms.hxx"
#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>
#include <gtest/gtest.h>
#include <cstdint>
#include <iostream>
#include <random>
#include <vector>
namespace TestForms {
using namespace Forms;
using namespace std;
template <typename T> class arbitrary {
  vector<T> vals;
  size_t pos;
  random_device rd;
  default_random_engine gen;
  uniform_int_distribution<> dist;
public:
  arbitrary() : pos(0), gen(rd()), dist(-1000 * 1000, 1000 * 1000) {
    for (int i = 0; i <= 10; ++i)
      vals.push_back(i);
    {
      int n = vals.size();
      for (int i = 0; i < n; ++i)
        vals.push_back(-vals.at(i));
    }
    {
      int n = vals.size();
      for (int i = 0; i < n; ++i)
        vals.push_back(1000 * vals.at(i));
    }
  }
  T next() {
    if (pos < vals.size())
      return vals.at(pos++);
    return dist(gen);
  }
};
template <typename T, int D, int R> void TestVector() {
  arbitrary<T> arb;
  const int n = 100;
  for (int i = 0; i < n; ++i) {
    constexpr int N = form<T, D, R>::N;
    form<T, D, R> x, y, z;
    for (int c = 0; c < N; ++c) {
      x.comps[c] = arb.next();
      y.comps[c] = arb.next();
      z.comps[c] = arb.next();
    }
    const auto zero = form<T, D, R>();
    const T a = arb.next();
    const T b = arb.next();
    const T szero = 0;
    const T sone = 1;
    // Addition:
    // Associativity
    EXPECT_EQ((x + y) + z, x + (y + z));
    // Left zero
    EXPECT_EQ(zero + x, x);
    // Right zero
    EXPECT_EQ(x + zero, x);
    // Commutativity
    EXPECT_EQ(x + y, y + x);
    // Left inverse
    EXPECT_EQ((-x) + x, zero);
    // Right inverse
    EXPECT_EQ(x + (-x), zero);
    // Subtraction
    EXPECT_EQ(x - y, x + (-y));
    // Scaling:
    // Associativity
    EXPECT_EQ((a * b) * x, a * (b * x));
    // Left Unit
    EXPECT_EQ(sone * x, x);
    // Right Unit
    EXPECT_EQ(x * sone, x);
    // Left Zero
    EXPECT_EQ(szero * x, zero);
    // Right Zero
    EXPECT_EQ(x * szero, zero);
    // Commutativity
    EXPECT_EQ(a * x, x * a);
    if (false) {
      if (a != szero) {
        // Left inverse
        EXPECT_EQ((sone / a) * (a * x), x);
        // Right inverse
        EXPECT_EQ((x * a) * (sone / a), x);
        // Division
        EXPECT_EQ(x / a, x * (sone / a));
      }
    }
    // Vector distributivity
    EXPECT_EQ(a * (x + y), a * x + a * y);
    // Scalar distributivity
    EXPECT_EQ((a + b) * x, a * x + b * x);
  }
}
template <typename T, int D, int Rx, int Ry> void TestWedge() {
  constexpr int R = Rx + Ry;
  arbitrary<T> arb;
  const int n = 100;
  for (int i = 0; i < n; ++i) {
    constexpr int Nx = form<T, D, Rx>::N;
    constexpr int Ny = form<T, D, Ry>::N;
    form<T, D, Rx> x;
    for (int c = 0; c < Nx; ++c)
      x.comps[c] = arb.next();
    form<T, D, Ry> y;
    for (int c = 0; c < Ny; ++c)
      y.comps[c] = arb.next();
    const auto xzero = form<T, D, Rx>();
    const auto yzero = form<T, D, Ry>();
    const auto rzero = form<T, D, R>();
    const T a = arb.next();
    const T b = arb.next();
    const T szero = 0;
    const T sone = 1;
    // Scaling:
    // Associativity
    EXPECT_EQ(a * wedge(x, y), wedge(a * x, y));
    // Left Zero
    EXPECT_EQ(wedge(xzero, y), rzero);
    // Right Zero
    EXPECT_EQ(wedge(x, yzero), rzero);
    // Anti-commutativity
    constexpr bool s = (Rx & Ry) & 1;
    EXPECT_EQ(wedge(x, y), bitsign(s) * wedge(y, x));
  }
}
extern "C" void TestForms_Test(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  CCTK_VINFO("Testing...");
  typedef int64_t val_t;
  TestVector<val_t, 0, 0>();
  TestVector<val_t, 1, 0>();
  TestVector<val_t, 1, 1>();
  TestVector<val_t, 2, 0>();
  TestVector<val_t, 2, 1>();
  TestVector<val_t, 2, 2>();
  TestVector<val_t, 3, 0>();
  TestVector<val_t, 3, 1>();
  TestVector<val_t, 3, 2>();
  TestVector<val_t, 3, 3>();
  TestWedge<val_t, 0, 0, 0>();
  TestWedge<val_t, 1, 0, 0>();
  TestWedge<val_t, 1, 0, 1>();
  TestWedge<val_t, 1, 1, 0>();
  TestWedge<val_t, 2, 0, 0>();
  TestWedge<val_t, 2, 0, 1>();
  TestWedge<val_t, 2, 0, 2>();
  TestWedge<val_t, 2, 1, 0>();
  TestWedge<val_t, 2, 1, 1>();
  TestWedge<val_t, 2, 2, 0>();
  TestWedge<val_t, 3, 0, 0>();
  TestWedge<val_t, 3, 0, 1>();
  TestWedge<val_t, 3, 0, 2>();
  TestWedge<val_t, 3, 0, 3>();
  TestWedge<val_t, 3, 1, 0>();
  TestWedge<val_t, 3, 1, 1>();
  TestWedge<val_t, 3, 1, 2>();
  TestWedge<val_t, 3, 2, 0>();
  TestWedge<val_t, 3, 2, 1>();
  TestWedge<val_t, 3, 3, 0>();
  CCTK_VINFO("Done.");
}
} // namespace TestForms

CarpetX/Forms