From 6e1eaaf13a13584294c6a090ef40500cd786f56f Mon Sep 17 00:00:00 2001 From: Jesper Jensen Date: Sun, 23 Mar 2025 11:26:53 +0100 Subject: Initial commit --- src/leven.c | 331 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/leven.h | 57 +++++++++++ src/log.h | 13 +++ src/mat.h | 43 ++++++++ src/parse.c | 172 +++++++++++++++++++++++++++++++ src/parse.h | 7 ++ 6 files changed, 623 insertions(+) create mode 100644 src/leven.c create mode 100644 src/leven.h create mode 100644 src/log.h create mode 100644 src/mat.h create mode 100644 src/parse.c create mode 100644 src/parse.h (limited to 'src') diff --git a/src/leven.c b/src/leven.c new file mode 100644 index 0000000..9ed2ab3 --- /dev/null +++ b/src/leven.c @@ -0,0 +1,331 @@ +#include +#include +#include +#include +#include + +#include "leven.h" +#include "log.h" + +#if 1 +#define debug2D(array, len) debug2D_impl(#array, array, len) +static void debug2D_impl(char* name, mat_uint32_t mat, uint32_t len) { + logb("2D Array %s [%ldx%d]", name, mat.stride, len); + for(uint32_t y = 0; y < len; y++) { + lognl(); + for(uint32_t x = 0; x < mat.stride; x++) { + logc("%03d ", *imat_uint32_t(mat, x, y)); + } + } + loge(); +} +#endif + +static inline int min(uint32_t a, uint32_t b) { + return (a < b) ? a : b; +} + +void string_edit_distance(const nid *a, size_t la, const nid *b, size_t lb, const mat_uint32_t weight, mat_uint32_t cost) { + *imat_uint32_t(cost, 0, 0) = 0; + + for (int i = 1; i <= la; i++) { + nid ac = a[i-1]; + assert(ac != 0); + *imat_uint32_t(cost, i, 0) = *imat_uint32_t(cost, i-1, 0) + *imat_uint32_t(weight, ac, 0); + } + + for (int i = 1; i <= lb; i++) { + nid bc = b[i-1]; + assert(bc != 0); + *imat_uint32_t(cost, 0, i) = *imat_uint32_t(cost, 0, i-1) + *imat_uint32_t(weight, 0, bc); + } + + for(int j = 1; j <= lb; j++) { + for(int i = 1; i <= la; i++) { + nid ac = a[i-1]; + assert(ac != 0); + nid bc = b[j-1]; + assert(bc != 0); + + // Remove + *imat_uint32_t(cost, i, j) = *imat_uint32_t(cost, i-1, j ) + *imat_uint32_t(weight, ac, 0); + // Add + *imat_uint32_t(cost, i, j) = min(*imat_uint32_t(cost, i , j-1) + *imat_uint32_t(weight, 0, bc), *imat_uint32_t(cost, i, j)); + // Replace + *imat_uint32_t(cost, i, j) = min(*imat_uint32_t(cost, i-1, j-1) + *imat_uint32_t(weight, ac, bc), *imat_uint32_t(cost, i, j)); + } + } + debug2D(cost, lb+1); +} + +void string_edit_alignment(const nid *a, size_t la, const nid *b, size_t lb, const mat_uint32_t weight, const mat_uint32_t cost, uint32_t *alignment) { + int i = la; + int j = lb; + while(j > 0) { + // The order of these operations dictate our preference if we have multiple equivalent answers + if(i == 0) { + // If we run out of characters the rest must be added + alignment[j-1] = -1; + j--; + } else if(*imat_uint32_t(cost, i, j) == *imat_uint32_t(cost, i-1, j-1) + *imat_uint32_t(weight, a[i-1], b[j-1])) { + alignment[j-1] = i-1; + i--; + j--; + } else if(*imat_uint32_t(cost, i, j) == *imat_uint32_t(cost, i-1, j) + *imat_uint32_t(weight, a[i-1], 0)) { + i--; + } else if(*imat_uint32_t(cost, i, j) == *imat_uint32_t(cost, i , j-1) + *imat_uint32_t(weight, 0, b[j-1])) { + alignment[j-1] = -1; + j--; + } else { + abort(); + } + } +} + +void constrained_tree_distance( + const struct Tree a, + const struct Tree b, + const mat_uint32_t cost, // The cost matrix to map a node from a (x-axis) to a node from b (y-axis) + mat_uint32_t cost_n, // The resulting computed cost matrixes node and forest. Size a.len x b.len + mat_uint32_t cost_f, + mat_uint32_t cost_s // Scratch space to calculate the edit distance between subtrees. Size a.adj.stride x b.adj.stride. +) { + *imat_uint32_t(cost_n, 0, 0) = 0; + *imat_uint32_t(cost_f, 0, 0) = 0; + + // Outer edges is creating/deleting the node + for(size_t i = a.len; i > 0; i--) { + *imat_uint32_t(cost_f, i, 0) = 0; + for(size_t j = 0; j < a.adj.stride && *imat_nid(a.adj, j, i-1) != 0; j++) { + *imat_uint32_t(cost_f, i, 0) += *imat_uint32_t(cost_n, *imat_nid(a.adj, j, i-1), 0); + } + *imat_uint32_t(cost_n, i, 0) = *imat_uint32_t(cost_f, i, 0) + *imat_uint32_t(cost, i, 0); + } + + for(size_t i = b.len; i > 0; i--) { + *imat_uint32_t(cost_f, 0, i) = 0; + for(size_t j = 0; j < b.adj.stride && *imat_nid(b.adj, j, i-1) != 0; j++) { + *imat_uint32_t(cost_f, 0, i) += *imat_uint32_t(cost_n, 0, *imat_nid(b.adj, j, i-1)); + } + *imat_uint32_t(cost_n, 0, i) = *imat_uint32_t(cost_f, 0, i) + *imat_uint32_t(cost, 0, i); + } + + for(size_t j = b.len; j > 0; j--) { + size_t b_adj_len = 0; + while(b_adj_len < b.adj.stride && *imat_nid(b.adj, b_adj_len, j-1) != 0) b_adj_len++; + + for(size_t i = a.len; i > 0; i--) { + size_t a_adj_len = 0; + while(a_adj_len < a.adj.stride && *imat_nid(a.adj, a_adj_len, i-1) != 0) a_adj_len++; + + string_edit_distance(imat_nid(a.adj, 0, i-1), a_adj_len, imat_nid(b.adj, 0, j-1), b_adj_len, cost_n, cost_s); + uint32_t min_cost = *imat_uint32_t(cost_s, a_adj_len, b_adj_len); + + if(a_adj_len > 0) { + uint32_t temp_min = UINT32_MAX; + for(uint32_t k = 0; k < a_adj_len; k++) { + uint32_t cost = *imat_uint32_t(cost_f, *imat_nid(a.adj, k, i-1), j) - *imat_uint32_t(cost_f, *imat_nid(a.adj, k, i-1), 0); + if(temp_min > cost) + temp_min = cost; + } + min_cost = min(min_cost, *imat_uint32_t(cost_f, i, 0) + temp_min); + } + + if(b_adj_len > 0) { + uint32_t temp_min = UINT32_MAX; + for(uint32_t k = 0; k < b_adj_len; k++) { + uint32_t cost = *imat_uint32_t(cost_f, i, *imat_nid(b.adj, k, j-1)) - *imat_uint32_t(cost_f, 0, *imat_nid(b.adj, k, j-1)); + if(temp_min > cost) + temp_min = cost; + } + min_cost = min(min_cost, *imat_uint32_t(cost_f, 0, j) + temp_min); + } + + *imat_uint32_t(cost_f, i, j) = min_cost; + + min_cost = *imat_uint32_t(cost_f, i, j) + *imat_uint32_t(cost, i, j); + + if(a_adj_len > 0) { + uint32_t temp_min = UINT32_MAX; + for(uint32_t k = 0; k < a_adj_len; k++) { + uint32_t cost = *imat_uint32_t(cost_n, *imat_nid(a.adj, k, i-1), j) - *imat_uint32_t(cost_n, *imat_nid(a.adj, k, i-1), 0); + if(temp_min > cost) + temp_min = cost; + } + min_cost = min(min_cost, *imat_uint32_t(cost_n, i, 0) + temp_min); + } + + if(b_adj_len > 0) { + uint32_t temp_min = UINT32_MAX; + for(uint32_t k = 0; k < b_adj_len; k++) { + uint32_t cost = *imat_uint32_t(cost_n, i, *imat_nid(b.adj, k, j-1)) - *imat_uint32_t(cost_n, 0, *imat_nid(a.adj, k, j-1)); + if(temp_min > cost) + temp_min = cost; + } + min_cost = min(min_cost, *imat_uint32_t(cost_n, 0, j) + temp_min); + } + + *imat_uint32_t(cost_n, i, j) = min_cost; + } + } +} + +void constrained_tree_alignment ( + const struct Tree a, + const struct Tree b, + const mat_uint32_t cost, + const mat_uint32_t cost_n, // The resulting computed cost matrixes node and forest. Size a.len x b.len + const mat_uint32_t cost_f, + mat_uint32_t cost_s, + uint32_t *adj_alignment, + mat_uint32_t alignment +) { + + mat_uint32_t to_compute = { + .data = malloc((2 * (a.len * b.len)) * sizeof(uint32_t)), + .stride = 2 + }; + *imat_uint32_t(to_compute, 0, 0) = 1; + *imat_uint32_t(to_compute, 1, 0) = 1; + size_t to_compute_head = 1; + + while(to_compute_head > 0) { + to_compute_head--; + uint32_t i = *imat_uint32_t(to_compute, 0, to_compute_head); + uint32_t j = *imat_uint32_t(to_compute, 1, to_compute_head); + + if(i == -1) { + size_t remain = 1; + size_t cursor = j-1; + while(remain > 0) { + log("ADD %ld", cursor); + + size_t adj_len = 0; + while(adj_len < b.adj.stride && *imat_nid(b.adj, adj_len, cursor) != 0) adj_len++; + remain += adj_len; + + cursor++; + remain--; + } + continue; + } else if(j == -1) { + size_t remain = 1; + size_t cursor = i-1; + while(remain > 0) { + log("REMOVE %ld", cursor); + + size_t adj_len = 0; + while(adj_len < a.adj.stride && *imat_nid(a.adj, adj_len, cursor) != 0) adj_len++; + remain += adj_len; + + cursor++; + remain--; + } + continue; + } + + size_t a_adj_len = 0; + while(a_adj_len < a.adj.stride && *imat_nid(a.adj, a_adj_len, i-1) != 0) a_adj_len++; + + size_t b_adj_len = 0; + while(b_adj_len < b.adj.stride && *imat_nid(b.adj, b_adj_len, j-1) != 0) b_adj_len++; + + uint32_t min_cost_a = UINT32_MAX; + ssize_t min_a = -1; + if(a_adj_len > 0) { + for(uint32_t k = 0; k < a_adj_len; k++) { + uint32_t cost = *imat_uint32_t(cost_n, *imat_nid(a.adj, k, i-1), j) - *imat_uint32_t(cost_n, *imat_nid(a.adj, k, i-1), 0); + if(min_cost_a > cost) { + min_cost_a = cost; + min_a = k; + } + } + } + + uint32_t min_cost_b = UINT32_MAX; + ssize_t min_b = -1; + if(b_adj_len > 0) { + for(uint32_t k = 0; k < b_adj_len; k++) { + uint32_t cost = *imat_uint32_t(cost_n, i, *imat_nid(b.adj, k, j-1)) - *imat_uint32_t(cost_n, 0, *imat_nid(b.adj, k, j-1)); + if(min_cost_b > cost) { + min_cost_b = cost; + min_b = k; + } + } + } + + if(*imat_uint32_t(cost_n, i, j) == *imat_uint32_t(cost_f, i, j) + *imat_uint32_t(cost, i, j)) { + log("MATCH %d %d", i, j); + // Compare the forests underneath this node + string_edit_distance(imat_nid(a.adj, 0, i-1), a_adj_len, imat_nid(b.adj, 0, j-1), b_adj_len, cost_n, cost_s); + assert(*imat_uint32_t(cost_s, a_adj_len, b_adj_len) == *imat_uint32_t(cost_f, i, j)); + string_edit_alignment(imat_nid(a.adj, 0, i-1), a_adj_len, imat_nid(b.adj, 0, j-1), b_adj_len, cost_n, cost_s, adj_alignment); + + ssize_t a_cursor = a_adj_len-1; + ssize_t b_cursor = b_adj_len-1; + + while(b_cursor >= 0) { + if(adj_alignment[b_cursor] == -1) { + uint32_t *slot = imat_uint32_t(to_compute, 0, to_compute_head); + to_compute_head++; + slot[0] = -1; + slot[1] = *imat_nid(b.adj, b_cursor, j-1); + + b_cursor--; + continue; + } + + while(a_cursor > adj_alignment[b_cursor]) { + uint32_t *slot = imat_uint32_t(to_compute, 0, to_compute_head); + to_compute_head++; + slot[0] = *imat_nid(a.adj, a_cursor, i-1); + slot[1] = -1; + a_cursor--; + + continue; + } + + uint32_t *slot = imat_uint32_t(to_compute, 0, to_compute_head); + to_compute_head++; + slot[0] = *imat_nid(a.adj, a_cursor, i-1); + slot[1] = *imat_nid(b.adj, b_cursor, j-1); + + b_cursor--; + a_cursor--; + } + } else if(a_adj_len > 0 && *imat_uint32_t(cost_n, i, j) == *imat_uint32_t(cost_n, i, 0) + min_cost_a) { + // Remove this node and replace it with one of its children + log("REMOVE %ld %ld", i, j); + for(size_t x = 0; x < a_adj_len; x++) { + uint32_t *slot = imat_uint32_t(to_compute, 0, to_compute_head); + to_compute_head++; + if(x == min_a) { + slot[0] = *imat_nid(a.adj, x, i-1); + slot[1] = -1; + } else { + slot[0] = *imat_nid(a.adj, x, i-1); + slot[1] = i; + } + } + } else if(b_adj_len > 0 && *imat_uint32_t(cost_n, i, j) == *imat_uint32_t(cost_n, 0, j) + min_cost_b) { + // Inject a node here, moving the current node (from a) into the child forest. + log("ADD %ld %ld", i, j); + for(size_t x = 0; x < b_adj_len; x++) { + uint32_t *slot = imat_uint32_t(to_compute, 0, to_compute_head); + to_compute_head++; + if(x == min_b) { + slot[0] = i; + slot[1] = *imat_nid(b.adj, x, j-1); + } else { + slot[0] = -1; + slot[1] = *imat_nid(b.adj, x, j-1); + } + } + } else { + abort(); + } + } + + log("END"); + free(to_compute.data); +} diff --git a/src/leven.h b/src/leven.h new file mode 100644 index 0000000..318eaab --- /dev/null +++ b/src/leven.h @@ -0,0 +1,57 @@ +#pragma once + +#include +#include + +#define MAT_TYPE uint32_t +#include "mat.h" + +typedef struct { + uint32_t *data; + size_t stride[2]; +} mat3_uint32_t; + +static inline uint32_t* imat3_uint32_t(mat3_uint32_t mat, size_t x, size_t y, size_t z) { + assert(x >= 0); + assert(y >= 0 && y < mat.stride[0]); + assert(z >= 0 && z < mat.stride[1]); + return &mat.data[x + y*mat.stride[0] + z*(mat.stride[0]*mat.stride[1])]; +} +#define DECL_MAT3(NAME, T, X, Y, Z) \ + mat3_uint32_t NAME = { \ + .data = (T[X * Y * Z]){}, \ + .stride = {X, Y}, \ + } \ + +typedef uint16_t nid; + +#define MAT_TYPE nid +#include "mat.h" + +void string_edit_distance(const nid *a, size_t la, const nid *b, size_t lb, const mat_uint32_t weight, mat_uint32_t cost); +void string_edit_alignment(const nid *a, size_t la, const nid *b, size_t lb, const mat_uint32_t weight, const mat_uint32_t cost, uint32_t *alignment); + +struct Tree { + mat_nid adj; + size_t len; +}; + +void constrained_tree_distance( + struct Tree a, + struct Tree b, + mat_uint32_t cost, + mat_uint32_t cost_n, + mat_uint32_t cost_f, + mat_uint32_t cost_s +); + +void constrained_tree_alignment( + const struct Tree a, + const struct Tree b, + const mat_uint32_t cost, + const mat_uint32_t cost_n, + const mat_uint32_t cost_f, + mat_uint32_t cost_s, + uint32_t *adj_alignment, + mat_uint32_t alignment +); diff --git a/src/log.h b/src/log.h new file mode 100644 index 0000000..6954106 --- /dev/null +++ b/src/log.h @@ -0,0 +1,13 @@ +#pragma once + +#include + +#define logb(FORMAT, ...) fprintf(stderr, "%s:%d: " FORMAT, __FILE__, __LINE__ __VA_OPT__(,) __VA_ARGS__) +#define logc(FORMAT, ...) fprintf(stderr, FORMAT __VA_OPT__(,) __VA_ARGS__) +#define lognl() fprintf(stderr, "\n") +#define loge() lognl() + +#define log(FORMAT, ...) do { \ + logb(FORMAT __VA_OPT__(,) __VA_ARGS__); \ + loge(); \ +} while(0) diff --git a/src/mat.h b/src/mat.h new file mode 100644 index 0000000..2f116ea --- /dev/null +++ b/src/mat.h @@ -0,0 +1,43 @@ +#ifndef MAT_TYPE +#error MAT_TYPE must be defined +#endif + +#include +#include + +#define PP_CONCAT1(x,y) x ## y +#define PP_CONCAT(x,y) PP_CONCAT1(x,y) + +#ifndef MAT_H_SHARED +#define MAT_H_SHARED + +#define DECL_MAT(NAME, TYPE, X, Y) \ + PP_CONCAT(mat_, TYPE) NAME = { \ + .data = (TYPE[X * Y]){}, \ + .stride = X, \ + } + +#define DECL_MAT_DATA(NAME, TYPE, X, Y, ...) \ + PP_CONCAT(mat_, TYPE) NAME = { \ + .data = (TYPE[X * Y]){__VA_ARGS__}, \ + .stride = X, \ + } + +#endif + +#define SELF_T PP_CONCAT(mat_, MAT_TYPE) +#define F(NAME) PP_CONCAT(NAME##_, MAT_TYPE) + +typedef struct { + MAT_TYPE *data; + size_t stride; +} SELF_T; + +static inline MAT_TYPE *F(imat)(SELF_T mat, size_t x, size_t y) { + assert(x >= 0 && x < mat.stride); + assert(y >= 0); + return &mat.data[x + y*mat.stride]; +} + +#undef SELF_T +#undef MAT_TYPE diff --git a/src/parse.c b/src/parse.c new file mode 100644 index 0000000..6ab2007 --- /dev/null +++ b/src/parse.c @@ -0,0 +1,172 @@ +#include "parse.h" + +#include +#include +#include +#include +#include + +#include "log.h" + +enum ParsePhase { + PHASE_COUNT, + PHASE_BUILD, +}; + +struct ParseCtx { + char * const str; + char *cursor; + struct Tree *tree; + size_t tree_cursor; + size_t *chunks; + + enum ParsePhase phase; +}; + +static uint64_t max(uint64_t a, uint64_t b) { + return a > b ? a : b; +} + +static bool alnum(struct ParseCtx *ctx) { + return isalnum(ctx->cursor[0]); +} + +static int read_STag(struct ParseCtx *ctx, bool *self_close, size_t nodeId) { + // Record start position of the tag + if (ctx->phase == PHASE_BUILD) { + ctx->chunks[nodeId] = ctx->cursor - ctx->str; + } + + if(*ctx->cursor != '<') return 1; + ctx->cursor++; + + while(isalnum(*ctx->cursor) || *ctx->cursor == ' ' || *ctx->cursor == '"' || *ctx->cursor == '=') ctx->cursor++; + + // Check for self-closing tag + if(*ctx->cursor == '/') { + *self_close = true; + ctx->cursor++; + } + + if(*ctx->cursor != '>') return 1; + ctx->cursor++; + + return 0; +} + +static int read_Element(struct ParseCtx *ctx, size_t nodeId); + +static int read_Content(struct ParseCtx *ctx, size_t nodeId) { + int err; + uint64_t children = 0; + while(true) { + if(ctx->cursor[0] == '<' && ctx->cursor[1] != '/') { + size_t childId = ctx->tree_cursor++; + + err = read_Element(ctx, childId); + if(err != 0) return err; + + if(ctx->phase == PHASE_BUILD) + *imat_nid(ctx->tree->adj, children, nodeId) = childId+1; + + children++; + } else if(alnum(ctx)) { + ctx->cursor++; + } else { + break; + } + } + + if(ctx->phase == PHASE_COUNT) + ctx->tree->adj.stride = max(ctx->tree->adj.stride, children); + + if(ctx->phase == PHASE_BUILD) { + for(size_t i = children; i < ctx->tree->adj.stride; i++) { + *imat_nid(ctx->tree->adj, i, nodeId) = 0; + } + } + + return 0; +} + +static int read_ETag(struct ParseCtx *ctx) { + if(*ctx->cursor != '<') return 1; + ctx->cursor++; + + if(*ctx->cursor != '/') return 1; + ctx->cursor++; + + while(isalnum(*ctx->cursor)) ctx->cursor++; + + if(*ctx->cursor != '>') return 1; + ctx->cursor++; + + return 0; +} + +static int read_Element(struct ParseCtx *ctx, size_t nodeId) { + int err; + bool self_close = false; + + err = read_STag(ctx, &self_close, nodeId); + if(err != 0) return err; + + if(!self_close) { + err = read_Content(ctx, nodeId); + if(err != 0) return err; + err = read_ETag(ctx); + if(err != 0) return err; + } else if(ctx->phase == PHASE_BUILD) { + for(size_t i = 0; i < ctx->tree->adj.stride; i++) + *imat_nid(ctx->tree->adj, i, nodeId) = 0; + } + + if(ctx->phase == PHASE_COUNT) ctx->tree->len++; + return 0; +} + +// This follows a 2 phase approach. First we size out the tree to figure out how +// much data we are going to store. Then we allocate the space based on that +// count pass before then doing a second pass in the BUILD phase which fill in +// those data structures. +int parse_string(char *str, struct Tree *tree, size_t **chunks) { + assert(chunks != NULL); // Chunks parameter is required + + tree->adj.stride = 1; + tree->len = 0; + struct ParseCtx ctx = { + .str = str, + .cursor = str, + .tree = tree, + .tree_cursor = 1, + .phase = PHASE_COUNT, + }; + int err; + + err = read_Element(&ctx, 0); + if(err != 0) { + log("Parse error at %ld", ctx.cursor - ctx.str); + return err; + } + + tree->adj.data = malloc(tree->len * tree->adj.stride * sizeof(*tree->adj.data)); + + // Allocate memory for chunks + *chunks = malloc(tree->len * sizeof(size_t)); + if (*chunks == NULL) { + free(tree->adj.data); + return 1; + } + + // Reset the context for the next phase + ctx.phase = PHASE_BUILD; + ctx.cursor = str; + ctx.tree_cursor = 1; + ctx.chunks = *chunks; + + err = read_Element(&ctx, 0); + // The second phase we assume the tree is valid since it passed the first phase + assert(err == 0); + + return 0; +} diff --git a/src/parse.h b/src/parse.h new file mode 100644 index 0000000..516b171 --- /dev/null +++ b/src/parse.h @@ -0,0 +1,7 @@ +#include "leven.h" + +// Chunk offsets for each node in the tree +// chunks[nodeId] = start position of the tag in the original string +// The end position is implicitly the start of the next tag or the end of the string + +int parse_string(char *str, struct Tree *tree, size_t **chunks); -- cgit v1.2.3