Skip to content

Instantly share code, notes, and snippets.

@gmalysa
Last active May 3, 2024 08:55
Show Gist options
  • Save gmalysa/94f5b1e61dc3fb060181d3ddc6d10196 to your computer and use it in GitHub Desktop.
Save gmalysa/94f5b1e61dc3fb060181d3ddc6d10196 to your computer and use it in GitHub Desktop.
Calculate Residuals and extra info for Arcosphere Balancing
#include <stdbool.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#define DIM 8
struct vec {
int32_t x[DIM];
};
struct vec a = { .x = {0, -1, 0, 1, 0, 1, 0, -1} };
struct vec b = { .x = {1, -1, 0, 0, 1, -1, 0, 0} };
struct vec c = { .x = {0, 1, -1, -1, 1, 0, 0, 0} };
struct vec d = { .x = {0, 0, 0, -1, -1, 1, 1, 0} };
struct vec e = { .x = {-1, 0, 0, 0, 0, -1, 1, 1} };
struct vec f = { .x = {1, 0, 1, 0, -1, 0, -1, 0} };
struct vec g = { .x = {0, 0, -1, 1, 0, 0, -1, 1} };
struct vec h = { .x = {-1, 1, 1, 0, 0, 0, 0, -1} };
struct vec i = { .x = {-1, -1, 1, -1, 1, 1, -1, 1} };
struct vec u[7] = { 0 };
struct vec zero = {0};
int32_t dots[7] = {0};
static void vec_add(struct vec *d, int32_t c1, struct vec *s1, int32_t c2,
struct vec *s2)
{
for (size_t i = 0; i < DIM; i++)
d->x[i] = c1*s1->x[i] + c2*s2->x[i];
}
static int32_t l1(struct vec *s) {
return abs(s->x[0]) + abs(s->x[1]) + abs(s->x[2]) + abs(s->x[3])
+ abs(s->x[4]) + abs(s->x[5]) + abs(s->x[6]) + abs(s->x[7]);
}
static int32_t dot(struct vec *a, struct vec *b) {
int32_t res = 0;
for (size_t i = 0; i < DIM; ++i)
res += a->x[i]*b->x[i];
return res;
}
static void residual_floor(struct vec *r, struct vec *src) {
vec_add(r, 0, r, 1, src);
for (size_t k = 0; k < 7; ++k) {
int32_t T = dot(src, &u[k]);
vec_add(r, 1, r, -T/dots[k], &u[k]);
}
}
static void residual_round(struct vec *r, struct vec *src) {
vec_add(r, 0, r, 1, src);
for (size_t k = 0; k < 7; ++k) {
int32_t T = dot(src, &u[k]);
if ((T < 0) == (dots[k] < 0)) {
T = T + dots[k]/2;
}
else {
T = T - dots[k]/2;
}
vec_add(r, 1, r, -T/dots[k], &u[k]);
}
}
static void print_vec(struct vec *v) {
printf("[%d, %d, %d, %d, %d, %d, %d, %d]", v->x[0], v->x[1], v->x[2], v->x[3],
v->x[4], v->x[5], v->x[6], v->x[7]);
}
static bool in_residual_region(struct vec *vec) {
for (size_t k = 0; k < 7; ++k) {
int32_t d = dot(&u[k], vec);
if (abs(d) >= abs(dots[k])) {
return false;
}
}
return true;
}
static void (*search_test)(struct vec *s);
/*
* Change this between residual_floor and residual_round to compare the difference
* between the two rounding methods on the max error L1 and max error per arcosphere
*/
static void (*get_residual)(struct vec *r, struct vec *src) = residual_round;
size_t tried = 0;
void iterate_volumes(size_t coord) {
static struct vec s = {0};
if (coord == 7) {
int32_t sum = s.x[0] + s.x[1] + s.x[2] + s.x[3] + s.x[4] + s.x[5] + s.x[6];
s.x[7] = -sum;
tried += 1;
search_test(&s);
}
else {
for (int32_t k = -4; k < 5; ++k) {
s.x[coord] = k;
iterate_volumes(coord+1);
}
}
}
static int32_t max_l1 = 0;
static struct vec max_vec = {0};
static size_t outside = 0;
void max_l1_test(struct vec *s) {
struct vec r;
if (!in_residual_region(s)) {
outside += 1;
return;
}
get_residual(&r, s);
int32_t sum = l1(&r);
if (sum > max_l1) {
max_l1 = sum;
memcpy(&max_vec, &r, sizeof(r));
}
}
static size_t has_max_n = 0;
void equals_max_l1(struct vec *s) {
struct vec r;
if (!in_residual_region(s))
return;
get_residual(&r, s);
int32_t sum = l1(&r);
if (sum == max_l1) {
print_vec(&r);
printf("\n");
has_max_n += 1;
}
}
static int32_t max_coord_val = 0;
static size_t max_coord_coord = 0;
static struct vec max_coord_vec = {0};
void max_coord_test(struct vec *s) {
struct vec r;
if (!in_residual_region(s))
return;
get_residual(&r, s);
if (abs(r.x[max_coord_coord]) > max_coord_val) {
max_coord_val = abs(r.x[max_coord_coord]);
memcpy(&max_coord_vec, &r, sizeof(r));
}
}
int main(void) {
size_t k;
vec_add(&u[0], 1, &zero, 1, &c);
vec_add(&u[1], 1, &zero, 1, &e);
vec_add(&u[2], 1, &c, 1, &e);
vec_add(&u[2], 1, &u[2], 1, &a);
vec_add(&u[2], 1, &u[2], 1, &a);
vec_add(&u[3], 1, &zero, 1, &b);
vec_add(&u[4], 1, &zero, 1, &g);
vec_add(&u[5], 1, &b, 1, &g);
vec_add(&u[5], 1, &u[5], 1, &d);
vec_add(&u[5], 1, &u[5], 1, &d);
vec_add(&u[6], 1, &zero, 1, &i);
for (k = 0; k < 7; ++k) {
dots[k] = dot(&u[k], &u[k]);
}
for (k = 0; k < 7; ++k) {
printf("%zu: ", k);
print_vec(&u[k]);
printf(", len = %d\n", dots[k]);
}
search_test = max_l1_test;
iterate_volumes(0);
printf("Examined %zu vectors\n", tried);
printf("%zu vectors were outside the residual region\n", outside);
printf("-> %zu total points inside the residual region\n", tried - outside);
printf("Max L1: %d at ", max_l1);
print_vec(&max_vec);
printf("\n");
search_test = equals_max_l1;
iterate_volumes(0);
printf("%zu vectors achieve the max L1\n", has_max_n);
search_test = max_coord_test;
for (k = 0; k < DIM; ++k) {
max_coord_val = 0;
max_coord_coord = k;
iterate_volumes(0);
printf("Max imbalance of %d at coord %zu, ex ", max_coord_val, max_coord_coord);
print_vec(&max_coord_vec);
printf("\n");
}
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment