// 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