refactor

1 files changed, 730 insertions, 0 deletions

diff --git a/include/ufo.h b/include/ufo.h
new file mode 100755
index 0000000..a92cff1
--- /dev/null
+++ b/include/ufo.h
@@ -0,0 +1,730 @@
+// https://www.sciencedirect.com/science/article/abs/pii/S1568494620303574
+
+// ufo.h - Universal Function Originator
+//
+// Thread safety: Not thread-safe. External synchronization required.
+// Memory: Single-allocation design — struct and all arrays stored in one
+//   contiguous buffer for easy persistence. Use ufo_copy for deep clone;
+//   ufo_free to release. Serialization via ufo_serialize / ufo_deserialize.
+
+// f_j(X) = f_j(a_0j (⊙_1j) a_1j*x_1^b_1j (⊙_2j) a_2j*x_2^b_2j (⊙_3j) ... (⊙_nj)
+// a_nj*x_2^b_nj) g_j(X) = w_j * [f_j(X)]^(c_j) y = w_0 + g_1(X) ⊚_1 g_1(X) ⊚_2
+// ... ⊚_v g_v(X)
+
+#ifndef UFO_H
+#define UFO_H
+
+#include "type.h"
+
+struct ufo_t;
+
+typedef struct ufo_t {
+  // Input count
+  u32 n;
+  // The number of internal blocks
+  u32 v;
+  // Internal weights [v][n+1]
+  f64 *a;
+  // Internal exponents [v][n]
+  f64 *b;
+  // External weights [v+1]
+  f64 *w;
+  // External exponents [v]
+  f64 *c;
+  // Unary functions [v]
+  f64 (**f)(f64);
+  // External binary functions [v]
+  f64 (**fe)(f64, f64);
+  // Internal binary functions [v][n]
+  f64 (**fi)(f64, f64);
+} ufo_t;
+
+#define ufo_ajk(ufo, j, k) (*(ufo->a + (j) * (ufo->n + 1) + (k)))
+#define ufo_bjk(ufo, j, k) (*(ufo->b + (j) * (ufo->n) + (k)))
+#define ufo_fijk(ufo, j, k) (*(ufo->fi + (j) * (ufo->n) + (k)))
+
+ufo_t *ufo_init(u32 n, u32 v);
+void ufo_free(ufo_t *ufo);
+f64 ufo_evaluate(ufo_t *ufo, f64 *x);
+void ufo_print(ufo_t *ufo);
+ufo_t *ufo_copy(const ufo_t *src);
+void ufo_randomize(ufo_t *ufo, f64 min_w, f64 max_w);
+void ufo_mutate(ufo_t *ufo, f64 rate, f64 min_w, f64 max_w);
+void ufo_crossover(ufo_t *parent1, ufo_t *parent2, u32 split);
+u32 ufo_complexity(const ufo_t *ufo);
+void ufo_set_unary(ufo_t *ufo, u32 j, f64 (*fn)(f64));
+void ufo_set_external_binary(ufo_t *ufo, u32 j, f64 (*fn)(f64, f64));
+void ufo_set_internal_binary(ufo_t *ufo, u32 j, u32 k, f64 (*fn)(f64, f64));
+void ufo_evaluate_jacobian(ufo_t *ufo, f64 *x, f64 *jacobian);
+void ufo_evaluate_batch(ufo_t *ufo, f64 *x, u32 count, f64 *out);
+void ufo_reset(ufo_t *ufo);
+void ufo_simplify(ufo_t *ufo);
+f64 ufo_similar(const ufo_t *a, const ufo_t *b);
+void ufo_serialize(const ufo_t *ufo, u8 *buffer, u64 *size);
+int ufo_deserialize(ufo_t **out, const u8 *buffer, u64 size);
+u32 ufo_serialize_size(u32 n, u32 v);
+void ufo_gradient_weights(ufo_t *ufo, f64 *x, f64 *weights_grad);
+void ufo_to_formula_string(const ufo_t *ufo, char *buf, u32 buf_size);
+
+#endif // UFO_H
+
+#ifdef UFO_IMPLEMENTATION
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <tgmath.h>
+
+static f64 ufo_identity(f64);
+static f64 ufo_increment(f64);
+static f64 ufo_decrement(f64);
+static f64 ufo_reciprocal(f64);
+static f64 ufo_root(f64);
+static f64 ufo_absolute_value(f64);
+static f64 ufo_floor(f64);
+static f64 ufo_ceiling(f64);
+static f64 ufo_round(f64);
+static f64 ufo_exponent(f64);
+static f64 ufo_log_natural(f64);
+static f64 ufo_log_decimal(f64);
+static f64 ufo_log_binary(f64);
+static f64 ufo_sine(f64);
+static f64 ufo_cosine(f64);
+static f64 ufo_tangent(f64);
+
+static f64 ufo_addition(f64, f64);
+static f64 ufo_subtraction(f64, f64);
+static f64 ufo_division(f64, f64);
+static f64 ufo_multiplication(f64, f64);
+static f64 ufo_power(f64, f64);
+
+ufo_t *ufo_init(u32 n, u32 v) {
+  if (n == 0 || v == 0)
+    return NULL;
+  if (n > UINT32_MAX / 5)
+    return NULL;
+  {
+    u32 factor_a = 3 * n + 5;
+    if (factor_a < 3 * n)
+      return NULL;
+    if (v > UINT32_MAX / factor_a)
+      return NULL;
+  }
+  u32 N = 3 * v * n + 5 * v;
+  {
+    u32 factor_d = 2 * n + 3;
+    if (factor_d < 2 * n)
+      return NULL;
+    if (v > UINT32_MAX / factor_d)
+      return NULL;
+  }
+  u32 D = 2 * v * n + 3 * v;
+
+  size_t bytes_a = sizeof(f64) * (size_t)N;
+  size_t bytes_d = sizeof(f64 (*)(f64)) * (size_t)D;
+  if (N > 0 && bytes_a / N != sizeof(f64))
+    return NULL;
+  if (D > 0 && bytes_d / D != sizeof(f64 (*)(f64)))
+    return NULL;
+  if (bytes_a > SIZE_MAX - bytes_d)
+    return NULL;
+  size_t bytes = bytes_a + bytes_d;
+
+  void *buffer = malloc(bytes);
+  if (!buffer)
+    return NULL;
+
+  ufo_t *ufo = (ufo_t *)buffer;
+  ufo->n = n;
+  ufo->v = v;
+
+  // Assign everything
+  ufo->a = (f64 *)(ufo + 1);
+  ufo->b = ufo->a + v * (n + 1);
+  ufo->w = ufo->b + v * n;
+  ufo->c = ufo->w + v + 1;
+  ufo->f = (f64(**)(f64))(ufo->c + v);
+  ufo->fe = (f64(**)(f64, f64))(ufo->f + v);
+  ufo->fi = (f64(**)(f64, f64))(ufo->fe + v * n);
+
+  // Initialize all weights to one
+  f64 *tmp_fp = ufo->a;
+  while ((void *)tmp_fp != (void *)ufo->f)
+    *tmp_fp++ = 1.0;
+
+  // Initialize all unary functions to identity
+  f64 (**tmp_unary)(f64) = ufo->f;
+  while ((void *)tmp_unary < (void *)ufo->fe)
+    *tmp_unary++ = ufo_identity;
+
+  // Initialize all binary functions to addition
+  f64 (**tmp_binary)(f64, f64) = ufo->fe;
+  while ((u8 *)tmp_binary != (u8 *)buffer + bytes)
+    *tmp_binary++ = ufo_addition;
+
+  // Set bias to zero
+  ufo->w[v] = 0;
+
+  for (u32 j = 0; j < v; j++)
+    ufo_ajk(ufo, j, n) = 0;
+
+  return ufo;
+}
+
+void ufo_free(ufo_t *ufo) {
+  free(ufo);
+}
+
+static f64 ufo_ujk(ufo_t *ufo, f64 *x, u32 j, u32 k) {
+  (void)k;
+  if (!ufo || !x)
+    return 0;
+  f64 acc = ufo_ajk(ufo, j, ufo->n);
+  for (u32 i = ufo->n; i > 0; i--)
+    acc = ufo_fijk(ufo, j, i - 1)(acc, x[i - 1]);
+  return acc;
+}
+
+static f64 ufo_gj(ufo_t *ufo, f64 *x, u32 j) {
+  if (!ufo || !x)
+    return 0;
+  f64 base = ufo->f[j](ufo_ujk(ufo, x, j, 0));
+  if (base < 0 && fmod(ufo->c[j], 1.0) != 0.0)
+    return 0;
+  if (base == 0.0 && ufo->c[j] <= 0.0)
+    return 0.0;
+  return ufo->w[j] * pow(base, ufo->c[j]);
+}
+
+f64 ufo_evaluate(ufo_t *ufo, f64 *x) {
+  if (!ufo || !x)
+    return 0;
+  f64 y = 0;
+
+  // Calculate g_j(X) for every j (excluding constant term)
+  u32 j;
+  for (j = 0; j < ufo->v; j++)
+    y = ufo->fe[j](y, ufo_gj(ufo, x, j));
+
+  // Constant term
+  y = ufo->fe[0](y, ufo->w[ufo->v]);
+
+  return y;
+}
+
+void ufo_print(ufo_t *ufo) {
+  if (!ufo)
+    return;
+  printf("UFO: %p\n"
+         "\n"
+         "N: %u\n"
+         "V: %u\n"
+         "A: %p\n"
+         "B: %p\n"
+         "W: %p\n"
+         "C: %p\n"
+         "F: %p\n"
+         "FE: %p\n"
+         "FI: %p\n",
+         (void *)ufo, (unsigned)ufo->n, (unsigned)ufo->v,
+         (void *)ufo->a, (void *)ufo->b, (void *)ufo->w, (void *)ufo->c,
+         (void *)ufo->f, (void *)ufo->fe, (void *)ufo->fi);
+}
+
+static f64 ufo_identity(f64 x) { return x; }
+static f64 ufo_increment(f64 x) { return x + 1; }
+static f64 ufo_decrement(f64 x) { return x - 1; }
+static f64 ufo_reciprocal(f64 x) { return (x == 0) ? 0 : 1 / x; }
+static f64 ufo_root(f64 x) { return (x < 0) ? 0 : sqrt(x); }
+static f64 ufo_absolute_value(f64 x) { return fabs(x); }
+static f64 ufo_floor(f64 x) { return floor(x); }
+static f64 ufo_ceiling(f64 x) { return ceil(x); }
+static f64 ufo_round(f64 x) { return round(x); }
+static f64 ufo_exponent(f64 x) { return exp(x); }
+static f64 ufo_log_natural(f64 x) { return (x <= 0) ? 0 : log(x); }
+static f64 ufo_log_decimal(f64 x) { return (x <= 0) ? 0 : log10(x); }
+static f64 ufo_log_binary(f64 x) { return (x <= 0) ? 0 : log2(x); }
+static f64 ufo_sine(f64 x) { return sin(x); }
+static f64 ufo_cosine(f64 x) { return cos(x); }
+static f64 ufo_tangent(f64 x) { return tan(x); }
+
+__attribute__((unused)) static f64 (*UNARY[])(f64) = {
+    ufo_identity,   ufo_increment,      ufo_decrement,   ufo_reciprocal,
+    ufo_root,       ufo_absolute_value, ufo_floor,       ufo_ceiling,
+    ufo_round,      ufo_exponent,       ufo_log_natural, ufo_log_decimal,
+    ufo_log_binary, ufo_sine,           ufo_cosine,
+    ufo_tangent
+};
+
+f64 ufo_addition(f64 x0, f64 x1) { return (x0 + x1); }
+f64 ufo_subtraction(f64 x0, f64 x1) { return (x0 - x1); }
+f64 ufo_division(f64 x0, f64 x1) { return (x1 == 0) ? 0 : (x0 / x1); }
+f64 ufo_multiplication(f64 x0, f64 x1) { return (x0 * x1); }
+f64 ufo_power(f64 x0, f64 x1) { return (x0 >= 0 || fmod(x1, 1.0) == 0.0) ? pow(x0, x1) : 0; }
+
+__attribute__((unused)) static f64 (*BINARY[])(f64, f64) = {ufo_addition, ufo_subtraction, ufo_division,
+                                   ufo_multiplication, ufo_power};
+
+ufo_t *ufo_copy(const ufo_t *src) {
+  if (!src)
+    return NULL;
+  ufo_t *dst = ufo_init(src->n, src->v);
+  if (!dst)
+    return NULL;
+  memcpy(dst->a, src->a, sizeof(f64) * (src->v * (src->n + 1)));
+  memcpy(dst->b, src->b, sizeof(f64) * (src->v * src->n));
+  memcpy(dst->w, src->w, sizeof(f64) * (src->v + 1));
+  memcpy(dst->c, src->c, sizeof(f64) * src->v);
+  memcpy(dst->f, src->f, sizeof(f64 (*)(f64)) * src->v);
+  memcpy(dst->fe, src->fe, sizeof(f64 (*)(f64, f64)) * (src->v + src->v * src->n));
+  return dst;
+}
+
+void ufo_randomize(ufo_t *ufo, f64 min_w, f64 max_w) {
+  if (!ufo)
+    return;
+  f64 range = max_w - min_w;
+  for (u32 j = 0; j < ufo->v; j++) {
+    for (u32 k = 0; k < ufo->n + 1; k++)
+      ufo_ajk(ufo, j, k) = min_w + (f64)rand() / RAND_MAX * range;
+    for (u32 k = 0; k < ufo->n; k++)
+      ufo_bjk(ufo, j, k) = min_w + (f64)rand() / RAND_MAX * range;
+    ufo->w[j] = min_w + (f64)rand() / RAND_MAX * range;
+    ufo->c[j] = min_w + (f64)rand() / RAND_MAX * range;
+    ufo->f[j] = UNARY[(rand() % 16)];
+    u32 bi_idx = (rand() % 5);
+    ufo->fe[j] = BINARY[bi_idx];
+    for (u32 k = 0; k < ufo->n; k++)
+      ufo_fijk(ufo, j, k) = BINARY[(rand() % 5)];
+  }
+  ufo->w[ufo->v] = min_w + (f64)rand() / RAND_MAX * range;
+}
+
+void ufo_mutate(ufo_t *ufo, f64 rate, f64 min_w, f64 max_w) {
+  if (!ufo || rate <= 0.0)
+    return;
+  f64 range = max_w - min_w;
+  for (u32 j = 0; j < ufo->v; j++) {
+    for (u32 k = 0; k < ufo->n + 1; k++)
+      if ((f64)rand() / RAND_MAX < rate)
+        ufo_ajk(ufo, j, k) = min_w + (f64)rand() / RAND_MAX * range;
+    for (u32 k = 0; k < ufo->n; k++)
+      if ((f64)rand() / RAND_MAX < rate)
+        ufo_bjk(ufo, j, k) = min_w + (f64)rand() / RAND_MAX * range;
+    if ((f64)rand() / RAND_MAX < rate)
+      ufo->w[j] = min_w + (f64)rand() / RAND_MAX * range;
+    if ((f64)rand() / RAND_MAX < rate)
+      ufo->c[j] = min_w + (f64)rand() / RAND_MAX * range;
+    if ((f64)rand() / RAND_MAX < rate)
+      ufo->f[j] = UNARY[(rand() % 16)];
+    if ((f64)rand() / RAND_MAX < rate)
+      ufo->fe[j] = BINARY[(rand() % 5)];
+    for (u32 k = 0; k < ufo->n; k++)
+      if ((f64)rand() / RAND_MAX < rate)
+        ufo_fijk(ufo, j, k) = BINARY[(rand() % 5)];
+  }
+  if ((f64)rand() / RAND_MAX < rate)
+    ufo->w[ufo->v] = min_w + (f64)rand() / RAND_MAX * range;
+}
+
+void ufo_crossover(ufo_t *parent1, ufo_t *parent2, u32 split) {
+  if (!parent1 || !parent2 || parent1->n != parent2->n || parent1->v != parent2->v)
+    return;
+  if (split > 0 && split < parent1->v) {
+    u32 n = parent1->n;
+    for (u32 j = split; j < parent1->v; j++) {
+      u32 offset = j * (n + 1);
+      f64 tmp_a;
+      for (u32 k = 0; k < n + 1; k++) {
+        tmp_a = parent1->a[offset + k];
+        parent1->a[offset + k] = parent2->a[offset + k];
+        parent2->a[offset + k] = tmp_a;
+      }
+      {
+        u32 bob = j * n;
+        f64 tmp_b;
+        for (u32 k = 0; k < n; k++) {
+          tmp_b = parent1->b[bob + k];
+          parent1->b[bob + k] = parent2->b[bob + k];
+          parent2->b[bob + k] = tmp_b;
+        }
+      }
+      {
+        f64 tmp = parent1->w[j];
+        parent1->w[j] = parent2->w[j];
+        parent2->w[j] = tmp;
+        tmp = parent1->c[j];
+        parent1->c[j] = parent2->c[j];
+        parent2->c[j] = tmp;
+        f64 (*tmp_fn)(f64) = parent1->f[j];
+        parent1->f[j] = parent2->f[j];
+        parent2->f[j] = tmp_fn;
+        f64 (*tmp_fn2)(f64, f64) = parent1->fe[j];
+        parent1->fe[j] = parent2->fe[j];
+        parent2->fe[j] = tmp_fn2;
+        for (u32 k = 0; k < n; k++) {
+          u32 fi_off = j * n + k;
+          f64 (*tmp_fi)(f64, f64) = parent1->fi[fi_off];
+          parent1->fi[fi_off] = parent2->fi[fi_off];
+          parent2->fi[fi_off] = tmp_fi;
+        }
+      }
+    }
+  }
+}
+
+u32 ufo_complexity(const ufo_t *ufo) {
+  if (!ufo)
+    return 0;
+  u32 score = 0;
+  for (u32 j = 0; j < ufo->v; j++) {
+    score++;
+    f64 c_val = ufo->c[j];
+    if (c_val != 1.0 && c_val != 0.0)
+      score++;
+    f64 test = ufo->f[j](2.0);
+    if (test != 2.0)
+      score++;
+    for (u32 k = 0; k < ufo->n; k++) {
+      f64 b_val = ufo_bjk(ufo, j, k);
+      if (b_val != 1.0 && b_val != 0.0)
+        score++;
+    }
+  }
+  return score;
+}
+
+void ufo_set_unary(ufo_t *ufo, u32 j, f64 (*fn)(f64)) {
+  if (!ufo || j >= ufo->v)
+    return;
+  ufo->f[j] = fn;
+}
+
+void ufo_set_external_binary(ufo_t *ufo, u32 j, f64 (*fn)(f64, f64)) {
+  if (!ufo || j >= ufo->v)
+    return;
+  ufo->fe[j] = fn;
+}
+
+void ufo_set_internal_binary(ufo_t *ufo, u32 j, u32 k, f64 (*fn)(f64, f64)) {
+  if (!ufo || j >= ufo->v || k >= ufo->n)
+    return;
+  ufo_fijk(ufo, j, k) = fn;
+}
+
+void ufo_evaluate_jacobian(ufo_t *ufo, f64 *x, f64 *jacobian) {
+  if (!ufo || !x || !jacobian)
+    return;
+  f64 y0 = ufo_evaluate(ufo, x);
+  for (u32 i = 0; i < ufo->n; i++)
+    jacobian[i] = 0;
+  f64 eps = 1e-6;
+  f64 *xp = malloc(sizeof(f64) * ufo->n);
+  if (!xp)
+    return;
+  for (u32 i = 0; i < ufo->n; i++) {
+    memcpy(xp, x, sizeof(f64) * ufo->n);
+    xp[i] += eps;
+    f64 yp = ufo_evaluate(ufo, xp);
+    jacobian[i] = (yp - y0) / eps;
+  }
+  free(xp);
+}
+
+void ufo_evaluate_batch(ufo_t *ufo, f64 *x, u32 count, f64 *out) {
+  if (!ufo || !x || !out || count == 0)
+    return;
+  for (u32 i = 0; i < count; i++)
+    out[i] = ufo_evaluate(ufo, x + (i64)i * (i64)ufo->n);
+}
+
+void ufo_reset(ufo_t *ufo) {
+  if (!ufo)
+    return;
+  f64 *tmp_fp = ufo->a;
+  while ((void *)tmp_fp != (void *)ufo->f)
+    *tmp_fp++ = 1.0;
+  f64 (**tmp_unary)(f64) = ufo->f;
+  while ((void *)tmp_unary < (void *)ufo->fe)
+    *tmp_unary++ = ufo_identity;
+  f64 (**tmp_binary)(f64, f64) = ufo->fe;
+  u32 total_bin = ufo->v + ufo->v * ufo->n;
+  for (u32 i = 0; i < total_bin; i++)
+    *tmp_binary++ = ufo_addition;
+  ufo->w[ufo->v] = 0;
+  for (u32 j = 0; j < ufo->v; j++)
+    ufo_ajk(ufo, j, ufo->n) = 0;
+}
+
+static int ufo_fn_index_unary(f64 (*fn)(f64)) {
+  for (int i = 0; i < 16; i++)
+    if (UNARY[i] == fn)
+      return i;
+  return -1;
+}
+
+static int ufo_fn_index_binary(f64 (*fn)(f64, f64)) {
+  for (int i = 0; i < 5; i++)
+    if (BINARY[i] == fn)
+      return i;
+  return -1;
+}
+
+void ufo_simplify(ufo_t *ufo) {
+  if (!ufo || ufo->v < 2)
+    return;
+  u32 write = 0;
+  for (u32 read = 0; read < ufo->v; read++) {
+    if (ufo->w[read] == 0.0 && ufo->c[read] == 0.0)
+      continue;
+    if (read != write) {
+      u32 n = ufo->n;
+      memcpy(ufo->a + write * (n + 1), ufo->a + read * (n + 1), sizeof(f64) * (n + 1));
+      memcpy(ufo->b + write * n, ufo->b + read * n, sizeof(f64) * n);
+      ufo->w[write] = ufo->w[read];
+      ufo->c[write] = ufo->c[read];
+      ufo->f[write] = ufo->f[read];
+      ufo->fe[write] = ufo->fe[read];
+      {
+        size_t sz = sizeof(void *) * (size_t)n;
+        memcpy((void *)(ufo->fi + write * n), (void *)(ufo->fi + read * n), sz);
+      }
+    }
+    write++;
+  }
+  for (u32 j = write; j < ufo->v; j++) {
+    ufo->w[j] = 0;
+    ufo->c[j] = 0;
+    ufo->f[j] = ufo_identity;
+    ufo->fe[j] = ufo_addition;
+    for (u32 k = 0; k < ufo->n + 1; k++)
+      ufo_ajk(ufo, j, k) = 0;
+    for (u32 k = 0; k < ufo->n; k++) {
+      ufo_bjk(ufo, j, k) = 0;
+      ufo_fijk(ufo, j, k) = ufo_addition;
+    }
+  }
+}
+
+f64 ufo_similar(const ufo_t *a, const ufo_t *b) {
+  if (!a || !b || a->n != b->n || a->v != b->v)
+    return 0;
+  u32 n = a->n;
+  u32 v = a->v;
+  f64 total = 0;
+  u32 count = 0;
+  for (u32 j = 0; j < v; j++) {
+    for (u32 k = 0; k < n + 1; k++) {
+      total += (a->a[j * (n + 1) + k] == b->a[j * (n + 1) + k]) ? 1 : 0;
+      count++;
+    }
+    for (u32 k = 0; k < n; k++) {
+      total += (a->b[j * n + k] == b->b[j * n + k]) ? 1 : 0;
+      count++;
+    }
+    total += (a->w[j] == b->w[j]) ? 1 : 0;
+    count++;
+    total += (a->c[j] == b->c[j]) ? 1 : 0;
+    count++;
+    total += (a->f[j] == b->f[j]) ? 1 : 0;
+    count++;
+    total += (a->fe[j] == b->fe[j]) ? 1 : 0;
+    count++;
+    for (u32 k = 0; k < n; k++) {
+      total += (a->fi[j * n + k] == b->fi[j * n + k]) ? 1 : 0;
+      count++;
+    }
+  }
+  total += (a->w[v] == b->w[v]) ? 1 : 0;
+  count++;
+  return count ? total / count : 1;
+}
+
+u32 ufo_serialize_size(u32 n, u32 v) {
+  return sizeof(u32) * 2
+         + sizeof(f64) * v * (n + 1)
+         + sizeof(f64) * v * n
+         + sizeof(f64) * (v + 1)
+         + sizeof(f64) * v
+         + sizeof(u8) * v
+         + sizeof(u8) * v
+         + sizeof(u8) * v * n;
+}
+
+void ufo_serialize(const ufo_t *ufo, u8 *buffer, u64 *size) {
+  if (!ufo || !buffer || !size)
+    return;
+  u32 n = ufo->n;
+  u32 v = ufo->v;
+  *size = ufo_serialize_size(n, v);
+  u8 *p = buffer;
+  memcpy(p, &n, sizeof(u32));
+  p += sizeof(u32);
+  memcpy(p, &v, sizeof(u32));
+  p += sizeof(u32);
+  memcpy(p, ufo->a, sizeof(f64) * v * (n + 1));
+  p += sizeof(f64) * v * (n + 1);
+  memcpy(p, ufo->b, sizeof(f64) * v * n);
+  p += sizeof(f64) * v * n;
+  memcpy(p, ufo->w, sizeof(f64) * (v + 1));
+  p += sizeof(f64) * (v + 1);
+  memcpy(p, ufo->c, sizeof(f64) * v);
+  p += sizeof(f64) * v;
+  for (u32 j = 0; j < v; j++) {
+    p[0] = (u8)ufo_fn_index_unary(ufo->f[j]);
+    p++;
+  }
+  for (u32 j = 0; j < v; j++) {
+    p[0] = (u8)ufo_fn_index_binary(ufo->fe[j]);
+    p++;
+  }
+  for (u32 j = 0; j < v; j++)
+    for (u32 k = 0; k < n; k++) {
+      p[0] = (u8)ufo_fn_index_binary(ufo_fijk(ufo, j, k));
+      p++;
+    }
+}
+
+int ufo_deserialize(ufo_t **out, const u8 *buffer, u64 size) {
+  if (!out || !buffer)
+    return 0;
+  u64 min = sizeof(u32) * 2 + 1;
+  if (size < min)
+    return 0;
+  u32 n, v;
+  memcpy(&n, buffer, sizeof(u32));
+  memcpy(&v, buffer + sizeof(u32), sizeof(u32));
+  u64 needed = ufo_serialize_size(n, v);
+  if (size < needed)
+    return 0;
+  ufo_t *ufo = ufo_init(n, v);
+  if (!ufo)
+    return 0;
+  *out = ufo;
+  const u8 *p = buffer + sizeof(u32) * 2;
+  memcpy(ufo->a, p, sizeof(f64) * v * (n + 1));
+  p += sizeof(f64) * v * (n + 1);
+  memcpy(ufo->b, p, sizeof(f64) * v * n);
+  p += sizeof(f64) * v * n;
+  memcpy(ufo->w, p, sizeof(f64) * (v + 1));
+  p += sizeof(f64) * (v + 1);
+  memcpy(ufo->c, p, sizeof(f64) * v);
+  p += sizeof(f64) * v;
+  for (u32 j = 0; j < v; j++) {
+    int idx = (int)p[0];
+    p++;
+    ufo->f[j] = (idx >= 0 && idx < 16) ? UNARY[(size_t)idx] : ufo_identity;
+  }
+  for (u32 j = 0; j < v; j++) {
+    int idx = (int)p[0];
+    p++;
+    ufo->fe[j] = (idx >= 0 && idx < 5) ? BINARY[(size_t)idx] : ufo_addition;
+  }
+  for (u32 j = 0; j < v; j++)
+    for (u32 k = 0; k < n; k++) {
+      int idx = (int)p[0];
+      p++;
+      ufo_fijk(ufo, j, k) = (idx >= 0 && idx < 5) ? BINARY[(size_t)idx] : ufo_addition;
+    }
+  return 1;
+}
+
+void ufo_gradient_weights(ufo_t *ufo, f64 *x, f64 *weights_grad) {
+  if (!ufo || !x || !weights_grad)
+    return;
+  f64 y0 = ufo_evaluate(ufo, x);
+  f64 eps = 1e-6;
+  u32 n = ufo->n;
+  u32 v = ufo->v;
+  for (u32 j = 0; j < v; j++)
+    for (u32 k = 0; k < n + 1; k++) {
+      f64 orig = ufo_ajk(ufo, j, k);
+      ufo_ajk(ufo, j, k) = orig + eps;
+      weights_grad[j * (n + 1) + k] = (ufo_evaluate(ufo, x) - y0) / eps;
+      ufo_ajk(ufo, j, k) = orig;
+    }
+  for (u32 j = 0; j < v; j++) {
+    f64 orig = ufo->w[j];
+    ufo->w[j] = orig + eps;
+    weights_grad[v * (n + 1) + j] = (ufo_evaluate(ufo, x) - y0) / eps;
+    ufo->w[j] = orig;
+  }
+  {
+    u32 off = v * (n + 1) + v;
+    f64 orig = ufo->w[v];
+    ufo->w[v] = orig + eps;
+    weights_grad[off] = (ufo_evaluate(ufo, x) - y0) / eps;
+    ufo->w[v] = orig;
+  }
+  for (u32 j = 0; j < v; j++) {
+    u32 off = v * (n + 1) + v + 1 + j;
+    f64 orig = ufo->c[j];
+    ufo->c[j] = orig + eps;
+    weights_grad[off] = (ufo_evaluate(ufo, x) - y0) / eps;
+    ufo->c[j] = orig;
+  }
+}
+
+static void ufo_unary_name(char *buf, u32 sz, f64 (*fn)(f64)) {
+  int idx = ufo_fn_index_unary(fn);
+  if (idx < 0) {
+    snprintf(buf, sz, "<unknown>");
+    return;
+  }
+  static const char *names[] = {"id", "inc", "dec", "recip", "root", "abs",
+                                 "floor", "ceil", "round", "exp", "log", "log10",
+                                 "log2", "sin", "cos", "tan"};
+  snprintf(buf, sz, "%s", names[idx]);
+}
+
+static void ufo_binary_name(char *buf, u32 sz, f64 (*fn)(f64, f64)) {
+  int idx = ufo_fn_index_binary(fn);
+  if (idx < 0) {
+    snprintf(buf, sz, "?");
+    return;
+  }
+  static const char *names[] = {"+", "-", "/", "*", "^"};
+  snprintf(buf, sz, "%s", names[idx]);
+}
+
+void ufo_to_formula_string(const ufo_t *ufo, char *buf, u32 buf_size) {
+  if (!ufo || !buf || buf_size < 2) {
+    if (buf && buf_size >= 1) buf[0] = '\0';
+    return;
+  }
+  u32 n = ufo->n;
+  u32 v = ufo->v;
+  u32 pos = 0;
+  buf[pos++] = '(';
+  for (u32 j = 0; j < v; j++) {
+    if (j > 0) {
+      char bn[8];
+      ufo_binary_name(bn, sizeof(bn), ufo->fe[j]);
+      pos += (u32)snprintf(buf + pos, buf_size - pos, " %s ", bn);
+      if (pos >= buf_size - 1) break;
+    }
+    char un[16];
+    ufo_unary_name(un, sizeof(un), ufo->f[j]);
+    pos += (u32)snprintf(buf + pos, buf_size - pos, "%s(", un);
+    if (pos >= buf_size - 1) break;
+    pos += (u32)snprintf(buf + pos, buf_size - pos, "%g", ufo_ajk(ufo, j, n));
+    if (pos >= buf_size - 1) break;
+    for (u32 k = 0; k < n; k++) {
+      char bn[8];
+      ufo_binary_name(bn, sizeof(bn), ufo_fijk(ufo, j, k));
+      pos += (u32)snprintf(buf + pos, buf_size - pos, " %s ", bn);
+      if (pos >= buf_size - 1) break;
+      pos += (u32)snprintf(buf + pos, buf_size - pos, "(%g*x_%u)^%g",
+                            ufo_ajk(ufo, j, k), (unsigned)k, ufo_bjk(ufo, j, k));
+      if (pos >= buf_size - 1) break;
+    }
+    pos += (u32)snprintf(buf + pos, buf_size - pos, "))^%g", ufo->c[j]);
+    if (pos >= buf_size - 1) break;
+  }
+  if (pos < buf_size - 1)
+    pos += (u32)snprintf(buf + pos, buf_size - pos, ") + %g", ufo->w[v]);
+  buf[pos < buf_size ? pos : buf_size - 1] = '\0';
+}
+
+#endif // UFO_IMPLEMENTATION