CG-P-CG conversion working correctly

Tested the correctness by converting the rows of data/2222_inq.txt from
CG to P and back and verifying the start and end result matched.
This commit is contained in:
Miguel M 2023-04-30 20:18:39 +01:00
parent 8c1bf8a57f
commit 216836f7f9
2 changed files with 92 additions and 201 deletions

138
src/cg.c
View File

@ -2,6 +2,33 @@
#include <stdio.h>
static inline size_t CgJointIndex(size_t a_out, size_t b_out, size_t a_in __attribute__((unused)),
size_t b_in, size_t a_out_i, size_t b_out_i,
size_t a_in_i, size_t b_in_i) {
return b_out_i +
(b_out - 1) * (b_in_i + b_in * (a_out_i + (a_out - 1) * a_in_i));
}
static inline size_t CgAMarginalIndex(size_t a_out, size_t b_out, size_t a_in,
size_t b_in, size_t a_out_i,
size_t a_in_i) {
return (a_out - 1) * (b_out - 1) * a_in * b_in + a_out_i +
(a_out - 1) * a_in_i;
}
static inline size_t CgBMarginalIndex(size_t a_out, size_t b_out, size_t a_in,
size_t b_in, size_t b_out_i,
size_t b_in_i) {
return (a_out - 1) * (b_out - 1) * a_in * b_in + (a_out - 1) * a_in +
b_out_i + (b_out - 1) * b_in_i;
}
static inline size_t PIndex(size_t a_out, size_t b_out, size_t a_in __attribute__((unused)),
size_t b_in, size_t a_out_i, size_t b_out_i,
size_t a_in_i, size_t b_in_i) {
return b_out_i + b_out * (b_in_i + b_in * (a_out_i + a_out * a_in_i));
}
void FromCgToP(size_t a_out, size_t b_out, size_t a_in, size_t b_in,
const data_t *restrict cg, data_t *restrict p) {
// size_t p_row_size = a_out * a_in * b_out * b_in;
@ -15,25 +42,23 @@ void FromCgToP(size_t a_out, size_t b_out, size_t a_in, size_t b_in,
for (size_t b_in_i = 0; b_in_i < b_in; b_in_i++) {
for (size_t a_out_i = 0; a_out_i < a_out - 1; a_out_i++) {
for (size_t b_out_i = 0; b_out_i < b_out - 1; b_out_i++) {
size_t cg_index =
b_out_i +
(b_out - 1) * (b_in_i + b_in * (a_out_i + (a_out - 1) * a_in_i));
size_t p_index =
b_out_i + b_out * (b_in_i + b_in * (a_out_i + a_out * a_in_i));
size_t cg_index = CgJointIndex(a_out, b_out, a_in, b_in, a_out_i,
b_out_i, a_in_i, b_in_i);
size_t p_index = PIndex(a_out, b_out, a_in, b_in, a_out_i, b_out_i,
a_in_i, b_in_i);
p[p_index] = cg[cg_index];
}
{
size_t p_index = (b_out - 1) +
b_out * (b_in_i + b_in * (a_out_i + a_out * a_in_i));
size_t p_index = PIndex(a_out, b_out, a_in, b_in, a_out_i, b_out - 1,
a_in_i, b_in_i);
p[p_index] = 0;
}
}
for (size_t b_out_i = 0; b_out_i < b_out; b_out_i++) {
size_t p_index =
(b_out - 1) +
b_out * (b_in_i + b_in * ((a_out - 1) + a_out * a_in_i));
PIndex(a_out, b_out, a_in, b_in, a_in - 1, b_out_i, a_in_i, b_in_i);
p[p_index] = 0;
}
}
@ -45,7 +70,7 @@ void FromCgToP(size_t a_out, size_t b_out, size_t a_in, size_t b_in,
for (size_t a_out_i = 0; a_out_i < a_out; a_out_i++) {
for (size_t b_out_i = 0; b_out_i < b_out; b_out_i++) {
size_t p_index = b_out_i + b_out * (0 + b_in * (a_out_i + a_out * 0));
size_t p_index = PIndex(a_out, b_out, a_in, b_in, a_out_i, b_out_i, 0, 0);
// Negative sign comes from moving L to the LHS of the inequality
p[p_index] -= l;
}
@ -58,19 +83,16 @@ void FromCgToP(size_t a_out, size_t b_out, size_t a_in, size_t b_in,
// A's marginals:
for (size_t a_in_i = 0; a_in_i < a_in; a_in_i++) {
for (size_t a_out_i = 0; a_out_i < a_out - 1; a_out_i++) {
size_t cg_index = (a_out - 1) * (b_out - 1) * a_in * b_in + a_out_i +
(a_out - 1) * a_in_i;
size_t cg_index =
CgAMarginalIndex(a_out, b_out, a_in, b_in, a_out_i, a_in_i);
data_t a_marginal = cg[cg_index];
if (a_marginal == 0) {
continue;
}
for (size_t b_out_i = 0; b_out_i < b_out; b_out_i++) {
// printf("A(%zu|%zu) contributing to (%zu,%zu|%zu,0) with %" PRImDATA
//"\n",
// a_out_i, a_in_i, a_out_i, b_out_i, a_in_i, a_marginal);
size_t p_index =
b_out_i + b_out * (0 + b_in * (a_out_i + a_out * a_in_i));
PIndex(a_out, b_out, a_in, b_in, a_out_i, b_out_i, a_in_i, 0);
p[p_index] += a_marginal;
// PrintP(a_out, b_out, a_in, b_in, p);
}
@ -88,11 +110,8 @@ void FromCgToP(size_t a_out, size_t b_out, size_t a_in, size_t b_in,
}
for (size_t a_out_i = 0; a_out_i < a_out; a_out_i++) {
// printf("B(%zu|%zu) contributing to (%zu,%zu|0,%zu) with %" PRImDATA
//"\n",
// b_out_i, b_in_i, a_out_i, b_out_i, b_in_i, b_marginal);
size_t p_index =
b_out_i + b_out * (b_in_i + b_in * (a_out_i + a_out * 0));
PIndex(a_out, b_out, a_in, b_in, a_out_i, b_out_i, 0, b_in_i);
p[p_index] += b_marginal;
// PrintP(a_out, b_out, a_in, b_in, p);
}
@ -112,20 +131,18 @@ void FromPToCg(size_t a_out, size_t b_out, size_t a_in, size_t b_in,
for (size_t b_out_i = 0; b_out_i < b_out - 1; b_out_i++) {
for (size_t b_in_i = 0; b_in_i < b_in; b_in_i++) {
{ // Common factors
size_t p_index =
b_out_i + b_out * (b_in_i + b_in * (a_out_i + a_out * a_in_i));
size_t cg_index =
b_out_i +
(b_out - 1) *
(b_in_i + b_in * (a_out_i + (a_out - 1) * a_in_i));
size_t p_index = PIndex(a_out, b_out, a_in, b_in, a_out_i, b_out_i,
a_in_i, b_in_i);
size_t cg_index = CgJointIndex(a_out, b_out, a_in, b_in, a_out_i,
b_out_i, a_in_i, b_in_i);
cg[cg_index] = p[p_index];
}
}
}
{ // Zero A's marginals
size_t cg_marginal_index = (a_out - 1) * (b_out - 1) * a_in * b_in +
(a_out - 1) * a_in_i + a_out_i;
size_t cg_marginal_index =
CgAMarginalIndex(a_out, b_out, a_in, b_in, a_out_i, a_in_i);
cg[cg_marginal_index] = 0;
}
}
@ -133,9 +150,8 @@ void FromPToCg(size_t a_out, size_t b_out, size_t a_in, size_t b_in,
// Zero B's marginals
for (size_t b_out_i = 0; b_out_i < b_out - 1; b_out_i++) {
for (size_t b_in_i = 0; b_in_i < b_in; b_in_i++) {
size_t cg_marginal_index = (a_out - 1) * (b_out - 1) * a_in * b_in +
(a_out - 1) * a_in + (b_out - 1) * b_in_i +
b_out_i;
size_t cg_marginal_index =
CgBMarginalIndex(a_out, b_out, a_in, b_in, b_out_i, b_in_i);
cg[cg_marginal_index] = 0;
}
}
@ -147,9 +163,8 @@ void FromPToCg(size_t a_out, size_t b_out, size_t a_in, size_t b_in,
for (size_t a_in_i = 0; a_in_i < a_in; a_in_i++) {
for (size_t b_in_i = 0; b_in_i < b_in; b_in_i++) {
size_t oob_p_index =
(b_out - 1) +
b_out * (b_in_i + b_in * ((a_out - 1) + a_out * a_in_i));
size_t oob_p_index = PIndex(a_out, b_out, a_in, b_in, a_out - 1,
b_out - 1, a_in_i, b_in_i);
// L contribution
// The minus side comes from considering L to be on the RHS of the
@ -159,25 +174,22 @@ void FromPToCg(size_t a_out, size_t b_out, size_t a_in, size_t b_in,
for (size_t a_out_i = 0; a_out_i < a_out - 1; a_out_i++) {
// Contributions to terms with no output labels OOBs
for (size_t b_out_i = 0; b_out_i < b_out - 1; b_out_i++) {
size_t cg_index =
b_out_i +
(b_out - 1) *
(b_in_i + b_in * (a_out_i + (a_out - 1) * a_in_i));
size_t cg_index = CgJointIndex(a_out, b_out, a_in, b_in, a_out_i,
b_out_i, a_in_i, b_in_i);
cg[cg_index] += p[oob_p_index];
}
{ // Contribution to A's marginal
size_t cg_marginal_index = (a_out - 1) * (b_out - 1) * a_in * b_in +
(a_out - 1) * a_in_i + a_out_i;
size_t cg_marginal_index =
CgAMarginalIndex(a_out, b_out, a_in, b_in, a_out_i, a_in_i);
cg[cg_marginal_index] -= p[oob_p_index];
}
}
// Contributions to B's marginal
for (size_t b_out_i = 0; b_out_i < b_out - 1; b_out_i++) {
size_t cg_marginal_index = (a_out - 1) * (b_out - 1) * a_in * b_in +
(a_out - 1) * a_in + (b_out - 1) * b_in_i +
b_out_i;
size_t cg_marginal_index =
CgBMarginalIndex(a_out, b_out, a_in, b_in, b_out_i, b_in_i);
cg[cg_marginal_index] -= p[oob_p_index];
}
}
@ -188,20 +200,18 @@ void FromPToCg(size_t a_out, size_t b_out, size_t a_in, size_t b_in,
for (size_t a_in_i = 0; a_in_i < a_in; a_in_i++) {
for (size_t a_out_i = 0; a_out_i < a_out - 1; a_out_i++) {
size_t cg_marginal_index = (a_out - 1) * (b_out - 1) * a_in * b_in +
(a_out - 1) * a_in_i + a_out_i;
(a_out - 1) * a_in_i + a_out_i;
for (size_t b_in_i = 0; b_in_i < b_in; b_in_i++) {
size_t oob_p_index =
(b_out - 1) +
b_out * (b_in_i + b_in * (a_out_i + a_out * a_in_i));
size_t oob_p_index = PIndex(a_out, b_out, a_in, b_in, a_out_i,
b_out - 1, a_in_i, b_in_i);
// A's marginals contributions
cg[cg_marginal_index] += p[oob_p_index];
// Contributions to terms with no labels OOBs
for (size_t b_out_i = 0; b_out_i < b_out - 1; b_out_i++) {
size_t cg_index =
b_out_i +
(b_out - 1) * (b_in_i + b_in * (a_out_i + (a_out - 1) * a_in_i));
size_t cg_index = CgJointIndex(a_out, b_out, a_in, b_in, a_out_i,
b_out_i, a_in_i, b_in_i);
cg[cg_index] -= p[oob_p_index];
}
}
@ -211,27 +221,19 @@ void FromPToCg(size_t a_out, size_t b_out, size_t a_in, size_t b_in,
// 4. Contributions from the A output label OOBs.
for (size_t b_in_i = 0; b_in_i < b_in; b_in_i++) {
for (size_t b_out_i = 0; b_out_i < b_out - 1; b_out_i++) {
{ // B's marginals contributions
size_t cg_marginal_index = (a_out - 1) * (b_out - 1) * a_in * b_in +
(a_out - 1) * a_in + (b_out - 1) * b_in_i +
b_out_i;
for (size_t a_in_i = 0; a_in_i < a_in; a_in_i++) {
size_t oob_p_index =
(b_out - 1) +
b_out * (b_in_i + b_in * ((a_out - 1) + a_out * a_in_i));
cg[cg_marginal_index] += p[oob_p_index];
}
}
size_t cg_marginal_index =
CgBMarginalIndex(a_out, b_out, a_in, b_in, b_out_i, b_in_i);
// Other contributions
for (size_t a_in_i = 0; a_in_i < a_in; a_in_i++) {
size_t oob_p_index =
b_out_i + b_out * (b_in_i + b_in * ((a_out - 1) + a_out * a_in_i));
// B's marginals contributions
size_t oob_p_index = PIndex(a_out, b_out, a_in, b_in, a_out - 1,
b_out_i, a_in_i, b_in_i);
cg[cg_marginal_index] += p[oob_p_index];
// Other contributions
for (size_t a_out_i = 0; a_out_i < a_out - 1; a_out_i++) {
size_t cg_index =
b_out_i +
(b_out - 1) * (b_in_i + b_in * (a_out_i + (a_out - 1) * a_in_i));
size_t cg_index = CgJointIndex(a_out, b_out, a_in, b_in, a_out_i,
b_out_i, a_in_i, b_in_i);
cg[cg_index] -= p[oob_p_index];
}
}

View File

@ -3,21 +3,21 @@
#include "../includes/matrix.h"
#include "../includes/permutation.h"
#include "../includes/cg.h"
static void ReadArguments(char *from, size_t *restrict into, const char *varname)
{
if (!sscanf(from, "%zu", into))
{
static void ReadArguments(char *from, size_t *restrict into,
const char *varname) {
if (!sscanf(from, "%zu", into)) {
fprintf(stderr, "Failed to read %s.", varname);
exit(EXIT_FAILURE);
}
}
int main(int argc, char *argv[])
{
if (argc < 5)
{
fprintf(stderr, "Usage:\n ./main.exe <A outputs> <B outputs> <A inputs> <B inputs>");
int main(int argc, char *argv[]) {
if (argc < 5) {
fprintf(
stderr,
"Usage:\n ./main.exe <A outputs> <B outputs> <A inputs> <B inputs>");
exit(EXIT_FAILURE);
}
@ -27,133 +27,22 @@ int main(int argc, char *argv[])
ReadArguments(argv[3], &a_in, "A inputs");
ReadArguments(argv[4], &b_in, "B inputs");
size_t row_len = (a_out - 1) * (b_out - 1) * a_in * b_in + (a_out - 1) * a_in + (b_out - 1) * b_in + 1;
size_t row_len = (a_out - 1) * (b_out - 1) * a_in * b_in +
(a_out - 1) * a_in + (b_out - 1) * b_in + 1;
matrix_t matrix = ParseMatrix(row_len);
DebugPrintMatrix(&matrix);
//DebugPrintMatrix(&matrix);
size_t row_count = matrix.len / matrix.row_len;
permutation_generator_t perm_a_out = PermutationNewGenerator(a_out);
permutation_generator_t perm_b_out = PermutationNewGenerator(b_out);
permutation_generator_t perm_a_in = PermutationNewGenerator(a_in);
permutation_generator_t perm_b_in = PermutationNewGenerator(b_in);
_Bool *seen = calloc(row_count, sizeof(_Bool));
for (size_t base_row_i= 0; base_row_i < row_count - 1; base_row_i++)
{
if (seen[base_row_i])
{
continue;
}
data_t *base_row = matrix.head + base_row_i * matrix.row_len;
data_t base_l = base_row[matrix.row_len - 1];
for (size_t cmp_row_i = base_row_i + 1; cmp_row_i < row_count; cmp_row_i++)
{
data_t *cmp_row = matrix.head + cmp_row_i * matrix.row_len;
data_t cmp_l = cmp_row[matrix.row_len - 1];
if (base_l != cmp_l)
{
// If the L values are different, the rows cannot be different under permutation.
continue; // To next cmp_row_i
}
PermutationReset(&perm_a_in);
while (!perm_a_in.exhausted)
{
PermutationReset(&perm_b_in);
while (!perm_b_in.exhausted)
{
PermutationReset(&perm_a_out);
while (!perm_a_out.exhausted)
{
PermutationReset(&perm_b_out);
while (!perm_b_out.exhausted)
{
// Compare the rows while applying this permutation.
// If every entry is the same, the rows are equivalent.
_Bool any_different = 0;
for (size_t a_in_i = 0; a_in_i < a_in; a_in_i++)
{
for (size_t a_out_i = 0; a_out_i < a_out; a_out_i++)
{
for (size_t b_in_i = 0; b_in_i < b_in; b_in_i++)
{
for (size_t b_out_i = 0; b_out_i < b_out; b_out_i++)
{
size_t base_idx = b_out_i + b_out * (b_in_i + b_in * (a_out_i + a_out * a_in_i));
size_t cmp_idx = (perm_b_out.permutation[b_out_i] +
b_out * (perm_b_in.permutation[b_in_i] +
b_in * (perm_a_out.permutation[a_out_i] +
a_out * perm_a_in.permutation[a_in_i])));
if (base_row[base_idx] != cmp_row[cmp_idx])
{
any_different = 1;
goto stop_comparison;
}
}
}
}
}
stop_comparison:
if (!any_different)
{
// The two rows are equivalent under this particular permutation.
seen[cmp_row_i] = 1;
goto skip_permutations;
}
PermutationNext(&perm_b_out);
}
PermutationNext(&perm_a_out);
}
PermutationNext(&perm_b_in);
}
PermutationNext(&perm_a_in);
}
skip_permutations:;
} // End of loop over comparison rows
} // End of loop over base rows
// Print the non-equivalent rows to stdout.
for (size_t i = 0; i < row_count; i++)
{
if (!seen[i])
{
#if RELEASE
PrintMatrixRow(&matrix, i);
#else
printf("%zu:\n ", i);
data_t *row = matrix.head + i * matrix.row_len;
for (size_t a_in_i = 0; a_in_i < a_in; a_in_i++)
{
for (size_t a_out_i = 0; a_out_i < a_out; a_out_i++)
{
for (size_t b_in_i = 0; b_in_i < b_in; b_in_i++)
{
for (size_t b_out_i = 0; b_out_i < b_out; b_out_i++)
{
size_t idx = b_out_i + b_out * (b_in_i + b_in * (a_out_i + a_out * a_in_i));
printf("(%zu,%zu|%zu,%zu): %" PRImDATA ", ", a_out_i, b_out_i, a_in_i, b_in_i, row[idx]);
}
}
}
}
printf("L: %" PRImDATA, row[matrix.row_len - 1]);
printf("\n");
#endif
}
data_t *p = malloc(a_out * b_out * a_in * b_in * sizeof(data_t));
for (size_t i = 0; i < row_count; i++) {
data_t *row = matrix.head + i * row_len;
printf("%zu:--------------------------------\n", i);
PrintCg(a_out, b_out, a_in, b_in, row);
FromCgToP(a_out, b_out, a_in, b_in, row, p);
PrintP(a_out, b_out, a_in, b_in, p);
FromPToCg(a_out, b_out, a_in, b_in, row, p);
PrintCg(a_out, b_out, a_in, b_in, row);
}
PermutationFree(&perm_a_in);
PermutationFree(&perm_b_in);
PermutationFree(&perm_a_out);
PermutationFree(&perm_b_out);
free(seen);
return EXIT_SUCCESS;
}