Initial commit

1 files changed, 374 insertions, 0 deletions

diff --git a/c/importgfa.c b/c/importgfa.c
new file mode 100644
index 0000000..43cc948
--- /dev/null
+++ b/c/importgfa.c
@@ -0,0 +1,374 @@
+/// domagi --- DuckDB-powered pangenome Swiss Army knife
+/// Copyright © 2026 Arun Isaac <arunisaac@systemreboot.net>
+///
+/// This file is part of domagi.
+///
+/// domagi is free software: you can redistribute it and/or modify it under the
+/// terms of the GNU General Public License as published by the Free Software
+/// Foundation, either version 3 of the License, or (at your option) any later
+/// version.
+///
+/// domagi is distributed in the hope that it will be useful, but WITHOUT ANY
+/// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+/// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+/// details.
+///
+/// You should have received a copy of the GNU General Public License along with
+/// domagi. If not, see <https://www.gnu.org/licenses/>.
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <duckdb.h>
+
+#define NAME str_int_map
+#define KEY_TY const char*
+#define VAL_TY int
+#include <verstable.h>
+
+void die(const char *message)
+{
+  fprintf(stderr, "%s\n", message);
+  exit(EXIT_FAILURE);
+}
+
+static int ceildiv (int x, int y)
+{
+  return (x / y) + (x % y ? 1 : 0);
+}
+
+static int hashtable_get (str_int_map *table, const char *key)
+{
+  // Look up key in hash table and return the corresponding key. The program is
+  // aborted if the key is not found. This is ok for our limited purposes since
+  // all keys we look up are guaranteed to be present.
+  str_int_map_itr itr = vt_get(table, key);
+  assert(!vt_is_end(itr));
+  return itr.data->val;
+}
+
+static uint8_t orientation2int (char c)
+{
+  return c == '+' ? 0 : 1;
+}
+
+static void pass1_handle_s_line (char *line, str_int_map *segment_id_table, duckdb_appender *appender)
+{
+  static int segment_id = 0;
+
+  // Split S line.
+  strsep(&line, "\t");
+  char *segment_name = strsep(&line, "\t");
+  char *sequence = strsep(&line, "\t\n");
+
+  // Map segment name to ID in table.
+  vt_insert(segment_id_table, strdup(segment_name), segment_id);
+  // Append to database.
+  duckdb_append_int32(*appender, segment_id);
+  duckdb_append_varchar(*appender, segment_name);
+  duckdb_append_varchar(*appender, sequence);
+  duckdb_appender_end_row(*appender);
+  // Increment segment ID for next segment.
+  segment_id++;
+}
+
+static void pass1_handle_p_line (char *line, duckdb_appender *appender)
+{
+  static int path_id = 0;
+
+  // Split P line.
+  strsep(&line, "\t");
+  char *path_name = strsep(&line, "\t");
+
+  // Append to database. There are usually relatively few P lines. So, we're
+  // content with a simple row-wise appender.
+  duckdb_append_int32(*appender, path_id);
+  duckdb_append_varchar(*appender, path_name);
+  duckdb_appender_end_row(*appender);
+  // Increment path ID for next path.
+  path_id++;
+}
+
+#define LINK_COLUMN_COUNT 4
+static void process_l_lines (char **lines, size_t line_count, str_int_map *segment_id_table, duckdb_database *db)
+{
+  idx_t vector_size = duckdb_vector_size();
+  duckdb_logical_type int_type = duckdb_create_logical_type(DUCKDB_TYPE_INTEGER);
+  duckdb_logical_type orientation_type = duckdb_create_enum_type((const char *[]){"+", "-"}, 2);
+  duckdb_logical_type types[LINK_COLUMN_COUNT]
+    = {int_type, orientation_type, int_type, orientation_type};
+#pragma omp parallel
+  {
+    duckdb_connection con;
+    if (duckdb_connect(*db, &con) == DuckDBError)
+      die("Unable to connect to DuckDB database");
+    duckdb_appender appender;
+    if (duckdb_appender_create(con, NULL, "link", &appender) == DuckDBError)
+      die("Unable to create appender for link table");
+    duckdb_data_chunk chunk = duckdb_create_data_chunk(types, LINK_COLUMN_COUNT);
+    duckdb_data_chunk_set_size(chunk, vector_size);
+    int32_t *from_segment_data = duckdb_vector_get_data(duckdb_data_chunk_get_vector(chunk, 0));
+    uint8_t *from_orientation_data = duckdb_vector_get_data(duckdb_data_chunk_get_vector(chunk, 1));
+    int32_t *to_segment_data = duckdb_vector_get_data(duckdb_data_chunk_get_vector(chunk, 2));
+    uint8_t *to_orientation_data = duckdb_vector_get_data(duckdb_data_chunk_get_vector(chunk, 3));
+#pragma omp for
+    for (size_t slice_index=0; slice_index<ceildiv(line_count, vector_size); slice_index++) {
+      size_t chunk_index, line_index;
+      for (chunk_index=0, line_index=slice_index*vector_size;
+           (chunk_index<vector_size) && (line_index<line_count);
+           chunk_index++, line_index++) {
+        char *original_line = lines[line_index];
+        // Split L line and write results into data chunk.
+        strsep(&lines[line_index], "\t");
+        from_segment_data[chunk_index]
+          = hashtable_get(segment_id_table, strsep(&lines[line_index], "\t"));
+        from_orientation_data[chunk_index]
+          = orientation2int(*strsep(&lines[line_index], "\t"));
+        to_segment_data[chunk_index]
+          = hashtable_get(segment_id_table, strsep(&lines[line_index], "\t"));
+        to_orientation_data[chunk_index]
+          = orientation2int(*strsep(&lines[line_index], "\t\n"));
+        free(original_line);
+      }
+      // Write chunk to database.
+      duckdb_data_chunk_set_size(chunk, chunk_index);
+      if (duckdb_append_data_chunk(appender, chunk) == DuckDBError)
+        die("Unable to append chunk to link table");
+    }
+    duckdb_destroy_data_chunk(&chunk);
+    duckdb_appender_destroy(&appender);
+    duckdb_disconnect(&con);
+  }
+  duckdb_destroy_logical_type(&int_type);
+  duckdb_destroy_logical_type(&orientation_type);
+}
+#undef LINK_COLUMN_COUNT
+
+#define PATH_SEGMENT_COLUMN_COUNT 4
+#define SLICE_SIZE 65536
+static void process_p_line (char *line, size_t line_length, str_int_map *segment_id_table, int path_id, duckdb_database *db)
+{
+  // Split P line.
+  strsep(&line, "\t");
+  char *path_name = strsep(&line, "\t");
+  char *segment_names = strsep(&line, "\t\n");
+  // segment_names is likely a very long string, and we want to avoid calling
+  // strlen on it. Hence, we compute its length in this roundabout way from the
+  // line length.
+  size_t segment_names_length = line_length
+    // This is the stuff before segment names.
+    - strlen("P\t") - strlen(path_name) - strlen("\t")
+    // And, this is the stuff after it.
+    - (line ? strlen(line) : 0);
+
+  // Split segment_names string into slices.
+  // Allocate the slices array to an approximate upper bound; we need not bother
+  // with ceil.
+  char **slices = malloc((1+segment_names_length/SLICE_SIZE)*sizeof(char*));
+  size_t slice_count;
+  {
+    slices[0] = segment_names;
+    int i = 1;
+    while (slices[i-1] + SLICE_SIZE < segment_names + segment_names_length) {
+      // Start the ith slice SLICE_SIZE bytes away from the beginning of the
+      // (i-1)th slice.
+      slices[i] = slices[i-1] + SLICE_SIZE;
+      // But, extend the (i-1)th slice a bit more so that its end lines up with
+      // a comma.
+      char *next_comma;
+      if ((next_comma = strchr(slices[i], ','))) {
+        *next_comma = '\0';
+        slices[i] = next_comma + 1;
+      }
+      i++;
+    }
+    slice_count = i;
+  }
+
+  // Count commas in each slice.
+  int *commas = calloc(slice_count, sizeof(int));
+#pragma omp parallel for
+  for (size_t i=0; i<slice_count; i++) {
+    char *slice = slices[i];
+    char *comma;
+    while ((comma=strchr(slice, ','))) {
+      commas[i]++;
+      slice = comma + 1;
+    }
+  }
+  // Prefix-sum to get the number of commas in segment_names up to that slice.
+  for (size_t i=1; i<slice_count; i++)
+    commas[i] += commas[i-1];
+
+  // Split path segments and append them to the database.
+  {
+    idx_t vector_size = duckdb_vector_size();
+    duckdb_logical_type int_type = duckdb_create_logical_type(DUCKDB_TYPE_INTEGER);
+    duckdb_logical_type orientation_type = duckdb_create_enum_type((const char *[]){"+", "-"}, 2);
+    duckdb_logical_type types[PATH_SEGMENT_COLUMN_COUNT]
+      = {int_type, int_type, int_type, orientation_type};
+#pragma omp parallel
+    {
+      duckdb_connection con;
+      if (duckdb_connect(*db, &con) == DuckDBError)
+        die("Unable to connect to DuckDB database");
+      duckdb_appender appender;
+      if (duckdb_appender_create(con, "internal", "path_segment", &appender) == DuckDBError)
+        die("Unable to create appender for path_segment table");
+      duckdb_data_chunk chunk = duckdb_create_data_chunk(types, PATH_SEGMENT_COLUMN_COUNT);
+      int32_t *position_data = duckdb_vector_get_data(duckdb_data_chunk_get_vector(chunk, 0));
+      int32_t *path_id_data = duckdb_vector_get_data(duckdb_data_chunk_get_vector(chunk, 1));
+      int32_t *segment_id_data = duckdb_vector_get_data(duckdb_data_chunk_get_vector(chunk, 2));
+      uint8_t *segment_orientation_data = duckdb_vector_get_data(duckdb_data_chunk_get_vector(chunk, 3));
+#pragma omp for
+      for (size_t i=0; i<slice_count; i++) {
+        char orientation;
+        int position = i ? commas[i-1] : 0;
+        char *comma;
+        do {
+          size_t chunk_index = 0;
+          do {
+            char *orientation_ptr;
+            if ((comma = strchr(slices[i], ','))) {
+              // Comma is found; pick up the orientation character and blot out
+              // the comma.
+              orientation_ptr = comma - 1;
+              *comma = '\0';
+            } else
+              // Comma is not found; this must be the last segment—the character
+              // before the terminating null byte is the orientation character.
+              orientation_ptr = strchr(slices[i], '\0') - 1;
+            // Copy the orientation character and blot it out from the path
+            // segment name.
+            orientation = *orientation_ptr;
+            *orientation_ptr = '\0';
+            // Look up segment ID in hash table.
+            int path_segment_id = hashtable_get(segment_id_table, slices[i]);
+            // Append to data chunk.
+            position_data[chunk_index] = position;
+            path_id_data[chunk_index] = path_id;
+            segment_id_data[chunk_index] = path_segment_id;
+            segment_orientation_data[chunk_index] = orientation2int(orientation);
+            // Prepare for the next iteration.
+            position++;
+            chunk_index++;
+            if (comma)
+              slices[i] = comma + 1;
+          } while (chunk_index<vector_size && comma);
+          duckdb_data_chunk_set_size(chunk, chunk_index);
+          if (duckdb_append_data_chunk(appender, chunk) == DuckDBError)
+            die("Unable to append chunk to path_segment table");
+        } while (comma);
+      }
+      duckdb_destroy_data_chunk(&chunk);
+      duckdb_appender_destroy(&appender);
+      duckdb_disconnect(&con);
+    }
+    duckdb_destroy_logical_type(&int_type);
+    duckdb_destroy_logical_type(&orientation_type);
+  }
+
+  free(commas);
+  free(slices);
+}
+#undef SLICE_SIZE
+#undef PATH_SEGMENT_COLUMN_COUNT
+
+int import_gfa (const char *gfa_path, const char *duckdb_path)
+{
+  str_int_map segment_id_table;
+  vt_init(&segment_id_table);
+
+  FILE *gfa_fp = fopen(gfa_path, "r");
+  if (!gfa_fp) {
+    // TODO: Replace with die.
+    fprintf(stderr, "Unable to open GFA: %s\n", gfa_path);
+    exit(EXIT_FAILURE);
+  }
+
+  duckdb_database db;
+  duckdb_connection con;
+  if (duckdb_open(duckdb_path, &db) == DuckDBError) {
+    // TODO: Replace with die.
+    fprintf(stderr, "Unable to open DuckDB database: %s\n", duckdb_path);
+    exit(EXIT_FAILURE);
+  }
+  if (duckdb_connect(db, &con) == DuckDBError)
+    die("Unable to connect to DuckDB database");
+
+  char *line = NULL;
+  size_t n;
+  ssize_t line_length;
+  size_t l_line_count = 0;
+
+  // Pass 1
+  {
+    duckdb_appender segment_appender, path_appender;
+    if (duckdb_appender_create(con, NULL, "segment", &segment_appender) == DuckDBError)
+      die("Unable to create appender for segment table");
+    if (duckdb_appender_create(con, NULL, "path", &path_appender) == DuckDBError)
+      die("Unable to create appender for path table");
+    while ((getline(&line, &n, gfa_fp)) > 0) {
+      switch (line[0]) {
+      case 'S':
+        pass1_handle_s_line(line, &segment_id_table, &segment_appender);
+        break;
+      case 'P':
+        pass1_handle_p_line(line, &path_appender);
+        break;
+      case 'L':
+        l_line_count++;
+        break;
+      default:
+        break;
+      }
+    }
+    if (duckdb_appender_destroy(&segment_appender) == DuckDBError)
+      die("Unable to destroy segment appender");
+    if (duckdb_appender_destroy(&path_appender) == DuckDBError)
+      die("Unable to destroy path appender");
+  }
+  rewind(gfa_fp);
+  // Pass 2
+  {
+    char **l_lines = malloc(l_line_count*sizeof(char*));
+    size_t l_line_index = 0, path_id = 0;
+    while ((line_length = getline(&line, &n, gfa_fp)) > 0) {
+      switch (line[0]) {
+      case 'L':
+        l_lines[l_line_index] = strdup(line);
+        l_line_index++;
+        break;
+      case 'P':
+        process_p_line(line, line_length, &segment_id_table, path_id, &db);
+        path_id++;
+        break;
+      default:
+        break;
+      }
+    }
+    process_l_lines(l_lines, l_line_count, &segment_id_table, &db);
+    free(l_lines);
+  }
+
+  free(line);
+  fclose(gfa_fp);
+  duckdb_disconnect(&con);
+  duckdb_close(&db);
+  for (str_int_map_itr itr=vt_first(&segment_id_table); !vt_is_end(itr); itr=vt_next(itr))
+    free((char*)itr.data->key);
+  vt_cleanup(&segment_id_table);
+
+  return 0;
+}
+
+int main (int argc, char *argv[])
+{
+  if (argc != 3) {
+    fprintf(stderr, "Usage: %s <gfa_path> <duckdb_path>\n", argv[0]);
+    exit(EXIT_FAILURE);
+  }
+  return import_gfa(argv[1], argv[2]);
+}