From 4b2ea2ca00eea7feea036b7642c0c1443b8f77a1 Mon Sep 17 00:00:00 2001 From: Nicholas Noll Date: Sun, 3 May 2020 20:20:22 -0700 Subject: removed buggy qsort header and implemented myself --- .gitignore | 5 + Makefile | 6 +- bin/gentags | 3 +- compile_commands.json | 478 ++++++++++++++++++++++- include/libbio.h | 16 +- include/libn.h | 8 + include/libn/macro/map.h | 929 ++++++++++----------------------------------- include/libn/macro/qsort.h | 261 ++++--------- rules.mk | 1 + sys/libbio/align.c | 70 ++-- sys/libbio/phylo.c | 255 +++++++++++-- sys/libbio/rules.mk | 2 +- sys/libbio/simulate.c | 117 ++++++ sys/libbio/test.c | 18 +- sys/libmath/blas.c | 2 - sys/libn/rules.mk | 4 +- sys/libn/test.c | 53 ++- 17 files changed, 1253 insertions(+), 975 deletions(-) create mode 100644 sys/libbio/simulate.c diff --git a/.gitignore b/.gitignore index fc20b4d..941439d 100644 --- a/.gitignore +++ b/.gitignore @@ -10,5 +10,10 @@ sys/nixos bin/fasttree bin/mafft +GRTAGS +GTAGS +GPATH + .cscope +cscope.out tags diff --git a/Makefile b/Makefile index 2d6485c..e09510d 100644 --- a/Makefile +++ b/Makefile @@ -11,13 +11,13 @@ LIB_DIR := lib OBJ_DIR := build # C runtime library -CINIT := $(LIB_DIR)/crt/crt1.o $(LIB_DIR)/crt/x86_64/crti.o `gcc --print-file-name=crtbeginT.o` +CINIT := $(LIB_DIR)/crt/crt1.o $(LIB_DIR)/crt/x86_64/crti.o `gcc --print-file-name=crtbeginS.o` CFINI := `gcc --print-file-name=crtendS.o` $(LIB_DIR)/crt/x86_64/crtn.o # Flags, Libraries and Includes CFLAGS := -g -O3 -march=native \ - -ffast-math -fno-strict-aliasing -fwrapv -fms-extensions \ - -Wno-microsoft-anon-tag -Wno-incompatible-function-pointer-types + -ffast-math -fno-strict-aliasing -fwrapv -fms-extensions \ + -Wno-microsoft-anon-tag -Wno-incompatible-function-pointer-types STATIC := -static -nodefaultlibs -nostartfiles AFLAGS := -f elf64 INCS := -isystem $(INC_DIR)/vendor/libc -I $(INC_DIR) diff --git a/bin/gentags b/bin/gentags index f0494f0..f398349 100755 --- a/bin/gentags +++ b/bin/gentags @@ -10,10 +10,11 @@ echo "finding files ..." ROOT=/home/nolln/root find $ROOT \ -path "$ROOT/sys/*.[chs]" -prune -o \ + -path "$ROOT/vendor/musl/src/*.[chs]" -prune -o \ -path "$ROOT/include/*.h" > "$CSCOPE_DIR/files" echo "adding files to cscope db: $ROOT/cscope.db ..." -cscope -b -i "$CSCOPE_DIR/files" +cscope -b -k -I include/vendor/libc -I include/ -i "$CSCOPE_DIR/files" CSCOPE_DB="$ROOT/cscope.out" echo "exported CSCOPE_DB to: '$CSCOPE_DB'" diff --git a/compile_commands.json b/compile_commands.json index 68f271e..fac553e 100644 --- a/compile_commands.json +++ b/compile_commands.json @@ -4,12 +4,67 @@ "clang", "-c", "-g", + "-O0", + "-march=native", + "-ffast-math", "-fno-strict-aliasing", "-fwrapv", "-fms-extensions", "-Wno-microsoft-anon-tag", "-Wno-incompatible-function-pointer-types", - "-Iinclude", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libbio/io/newick.o", + "sys/libbio/io/newick.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libbio/io/newick.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-D_GNU_SOURCE", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libmath/test.o", + "sys/libmath/test.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libmath/test.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", "-o", "build/libn/flate.o", "sys/libn/flate.c" @@ -22,30 +77,169 @@ "clang", "-c", "-g", + "-O0", + "-march=native", + "-ffast-math", "-fno-strict-aliasing", "-fwrapv", "-fms-extensions", "-Wno-microsoft-anon-tag", "-Wno-incompatible-function-pointer-types", - "-Iinclude", + "-isystem", + "include/vendor/libc", + "-I", + "include", "-o", - "build/libn/bufio.o", - "sys/libn/bufio.c" + "build/libn/gz.o", + "sys/libn/gz.c" ], "directory": "/home/nolln/root", - "file": "sys/libn/bufio.c" + "file": "sys/libn/gz.c" }, { "arguments": [ "clang", "-c", "-g", + "-O0", + "-march=native", + "-ffast-math", "-fno-strict-aliasing", "-fwrapv", "-fms-extensions", "-Wno-microsoft-anon-tag", "-Wno-incompatible-function-pointer-types", - "-Iinclude", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libn/mmap.o", + "sys/libn/mmap.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libn/mmap.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-D_GNU_SOURCE", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libmath/blas.o", + "sys/libmath/blas.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libmath/blas.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libbio/test.o", + "sys/libbio/test.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libbio/test.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-ffreestanding", + "-fno-builtin", + "-nostdlib", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libc/string.o", + "sys/libc/string.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libc/string.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-ffreestanding", + "-fno-builtin", + "-nostdlib", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libc/stdio.o", + "sys/libc/stdio.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libc/stdio.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", "-o", "build/libn/error.o", "sys/libn/error.c" @@ -58,17 +252,287 @@ "clang", "-c", "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libbio/align.o", + "sys/libbio/align.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libbio/align.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libn/random.o", + "sys/libn/random.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libn/random.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", "-fno-strict-aliasing", "-fwrapv", "-fms-extensions", "-Wno-microsoft-anon-tag", "-Wno-incompatible-function-pointer-types", - "-Iinclude", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libbio/phylo.o", + "sys/libbio/phylo.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libbio/phylo.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libn/memory.o", + "sys/libn/memory.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libn/memory.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libn/io.o", + "sys/libn/io.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libn/io.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libn/string.o", + "sys/libn/string.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libn/string.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", "-o", "build/libn/coro.o", "sys/libn/coro.c" ], "directory": "/home/nolln/root", "file": "sys/libn/coro.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libn/sort.o", + "sys/libn/sort.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libn/sort.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libn/bufio.o", + "sys/libn/bufio.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libn/bufio.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libbio/io/fasta.o", + "sys/libbio/io/fasta.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libbio/io/fasta.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libbio/simulate.o", + "sys/libbio/simulate.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libbio/simulate.c" + }, + { + "arguments": [ + "clang", + "-c", + "-g", + "-O0", + "-march=native", + "-ffast-math", + "-fno-strict-aliasing", + "-fwrapv", + "-fms-extensions", + "-Wno-microsoft-anon-tag", + "-Wno-incompatible-function-pointer-types", + "-isystem", + "include/vendor/libc", + "-I", + "include", + "-o", + "build/libn/test.o", + "sys/libn/test.c" + ], + "directory": "/home/nolln/root", + "file": "sys/libn/test.c" } ] \ No newline at end of file diff --git a/include/libbio.h b/include/libbio.h index 9eebb0b..f802e22 100644 --- a/include/libbio.h +++ b/include/libbio.h @@ -25,14 +25,24 @@ typedef struct bio·Tree void *h; } bio·Tree; -// clade functions +/* clade manipulation */ error phylo·addchild(bio·Node* parent, bio·Node* child); error phylo·rmchild(bio·Node* parent, bio·Node* child); +/* clade statistics */ error phylo·countnodes(bio·Node *node, int *n); error phylo·countleafs(bio·Node *node, int *n); -error phylo·ladderize(bio·Node *root); +/* topological sorting */ +error phylo·ladderize(bio·Node *root); +double phylo·diameter(bio·Tree tree, int *len, bio·Node **path); + +/* generic computation on tree */ +void *phylo·postorder(bio·Node *clade, void *(*op)(bio·Node*, void*), void *ctx); +void *phylo·preorder(bio·Node *clade, void *(*op)(bio·Node*, void*), void *ctx); + +/* simple helpers */ +void phylo·getleafs(bio·Tree tree, bio·Node **leafs); /* newick i/o */ error bio·readnewick(io·Peeker stream, void*, bio·Tree* tree); @@ -66,7 +76,7 @@ enum { aln·K = 20, // kmer size (k <= 32) aln·L = 3, // number of kmers / hash - aln·N = 1000, // number of hashes + aln·N = 10, // number of hashes }; error aln·sketch(byte *seq, int l, uint64 *phi[aln·N], int *locs[aln·N]); diff --git a/include/libn.h b/include/libn.h index 063a408..8e30eae 100644 --- a/include/libn.h +++ b/include/libn.h @@ -311,3 +311,11 @@ void sort·floats(uintptr n, float arr[]); void sort·doubles(uintptr n, double arr[]); void sort·strings(uintptr n, byte* arr[]); + +// ----------------------------------------------------------------------------- +// fast random number generation + +error rng·init(uint64 seed); +double rng·random(); +bool rng·bernoulli(double f); +uint64 rng·randi(int max); diff --git a/include/libn/macro/map.h b/include/libn/macro/map.h index f81dc0f..2e65570 100644 --- a/include/libn/macro/map.h +++ b/include/libn/macro/map.h @@ -1,754 +1,237 @@ -// clang-format off -// -/* The MIT License - Copyright (c) 2008, 2009, 2011 by Attractive Chaos - Permission is hereby granted, free of charge, to any person obtaining - a copy of this software and associated documentation files (the - "Software"), to deal in the Software without restriction, including - without limitation the rights to use, copy, modify, merge, publish, - distribute, sublicense, and/or sell copies of the Software, and to - permit persons to whom the Software is furnished to do so, subject to - the following conditions: - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS - BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN - ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN - CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - SOFTWARE. -*/ - -/* - An example: -#include "khash.h" -HASH_MAP_INIT_INT(32, char) -int main() { - int ret, is·missing; - mapiter_t k; - khash_t(32) *h = map·init(32); - k = map·put(32, h, 5, &ret); - map·value(h, k) = 10; - k = map·get(32, h, 10); - is·missing = (k == map·end(h)); - k = map·get(32, h, 5); - map·del(32, h, k); - for (k = map·begin(h); k != map·end(h); ++k) - if (map·exist(h, k)) map·value(h, k) = 1; - map·destroy(32, h); - return 0; -} -*/ - -/* - 2013-05-02 (0.2.8): - * Use quadratic probing. When the capacity is power of 2, stepping - function i*(i+1)/2 guarantees to traverse each bucket. It is better than - double hashing on cache performance and is more robust than linear probing. In - theory, double hashing should be more robust than quadratic probing. However, - my implementation is probably not for large hash tables, because the second - hash function is closely tied to the first hash function, which reduce the - effectiveness of double hashing. Reference: - http://research.cs.vt.edu/AVresearch/hashing/quadratic.php 2011-12-29 (0.2.7): - * Minor code clean up; no actual effect. - 2011-09-16 (0.2.6): - * The capacity is a power of 2. This seems to dramatically improve the - speed for simple keys. Thank Zilong Tan for the suggestion. Reference: - - http://code.google.com/p/ulib/ - - http://nothings.org/computer/judy/ - * Allow to optionally use linear probing which usually has better - performance for random input. Double hashing is still the default as - it is more robust to certain non-random input. - * Added Wang's integer hash function (not used by default). This hash - function is more robust to certain non-random input. - 2011-02-14 (0.2.5): - * Allow to declare global functions. - 2009-09-26 (0.2.4): - * Improve portability - 2008-09-19 (0.2.3): - * Corrected the example - * Improved interfaces - 2008-09-11 (0.2.2): - * Improved speed a little in map·put() - 2008-09-10 (0.2.1): - * Added map·clear() - * Fixed a compiling error - 2008-09-02 (0.2.0): - * Changed to token concatenation which increases flexibility. - 2008-08-31 (0.1.2): - * Fixed a bug in map·get(), which has not been tested previously. - 2008-08-31 (0.1.1): - * Added destructor -*/ - -// Modified by Nicholas Noll 2019. -// Modified API to be more inline with codebase. - #pragma once -/*! - @header - Generic hash table library. +/* + * Based on khash.h + * Modified to make control more granular and similar to QSORT */ -#define AC_VERSION_HASH_H "0.2.8" - -/* compiler specific configuration */ - -// #if UINT_MAX == 0xffffffffu -// typedef unsigned int mapint32_t; -// #elif ULONG_MAX == 0xffffffffu -// typedef unsigned long mapint32_t; -// #endif -// -// #if ULONG_MAX == ULLONG_MAX -// typedef unsigned long mapint64_t; -// #else -// typedef unsigned long long mapint64_t; -// #endif - -#ifndef map·inline -#ifdef _MSC_VER -#define map·inline __inline -#else -#define map·inline inline -#endif -#endif /* map·inline */ - -#ifndef klib_unused -#if (defined __clang__ && __clang_major__ >= 3) || \ - (defined __GNUC__ && __GNUC__ >= 3) -#define klib_unused __attribute__((__unused__)) -#else -#define klib_unused -#endif -#endif /* klib_unused */ +static const double __ac_HASH_UPPER = 0.77; -typedef int32 mapiter; +#define _roundup32(x) \ + (--(x), (x) |= (x) >> 1, (x) |= (x) >> 2, (x) |= (x) >> 4, (x) |= (x) >> 8, \ + (x) |= (x) >> 16, ++(x)) #define __ac_isempty(flag, i) ((flag[i >> 4] >> ((i & 0xfU) << 1)) & 2) #define __ac_isdel(flag, i) ((flag[i >> 4] >> ((i & 0xfU) << 1)) & 1) #define __ac_iseither(flag, i) ((flag[i >> 4] >> ((i & 0xfU) << 1)) & 3) #define __ac_set_isdel_false(flag, i) \ - (flag[i >> 4] &= ~(1ul << ((i & 0xfU) << 1))) + (flag[i >> 4] &= ~(1ul << ((i & 0xfU) << 1))) #define __ac_set_isempty_false(flag, i) \ - (flag[i >> 4] &= ~(2ul << ((i & 0xfU) << 1))) + (flag[i >> 4] &= ~(2ul << ((i & 0xfU) << 1))) #define __ac_set_isboth_false(flag, i) \ - (flag[i >> 4] &= ~(3ul << ((i & 0xfU) << 1))) + (flag[i >> 4] &= ~(3ul << ((i & 0xfU) << 1))) #define __ac_set_isdel_true(flag, i) (flag[i >> 4] |= 1ul << ((i & 0xfU) << 1)) #define __ac_fsize(m) ((m) < 16 ? 1 : (m) >> 4) -#ifndef _roundup32 -#define _roundup32(x) \ - (--(x), \ - (x) |= (x) >> 1, \ - (x) |= (x) >> 2, \ - (x) |= (x) >> 4, \ - (x) |= (x) >> 8, \ - (x) |= (x) >> 16, \ - ++(x)) -#endif - -static const double __ac_HASH_UPPER = 0.77; - -#define __HASH_TYPE(name, map·key_t, map·val·t) \ - typedef struct map·##name##_s { \ - int32 n_buckets, size, n_occupied, upper_bound; \ - int32* flags; \ - map·key_t* keys; \ - map·val·t* vals; \ - } map·##name##_t; - -#define __HASH_PROTOTYPES(name, map·key_t, map·val·t) \ - extern map·##name##_t* map·init·##name(void); \ - extern void map·destroy_##name(map·##name##_t* h); \ - extern void map·clear_##name(map·##name##_t* h); \ - extern int32 map·get_##name(const map·##name##_t* h, mapkey_t key); \ - extern int map·resize_##name(map·##name##_t* h, int32 new_n_buckets); \ - extern int32 map·put_##name(map·##name##_t* h, mapkey_t key, int* ret); \ - extern void map·del_##name(map·##name##_t* h, int32 x); - -#define __HASH_IMPL( \ - name, SCOPE, mapkey_t, mapval·t, map·is·map, __hash_func, __hash_equal) \ - SCOPE map·##name##_t* map·init·##name(void) \ - { \ - return (map·##name##_t*)calloc(1, sizeof(map·##name##_t)); \ +#define MAP_STRUCT_BODY(key_t, val_t) \ + mem·Allocator heap; \ + void *h; \ + int32 n_buckets, size, n_occupied, upper_bound; \ + int32 *flags; \ + key_t *keys; \ + val_t *vals; + +#define MAP_MAKE(type) \ + type *map; \ + map = heap.alloc((h), 1, sizeof(*map)); \ + map->heap = heap; \ + map->h = h; \ + return map + +#define MAP_FREE(map) \ + map->heap.free(map->h, map->keys); \ + map->heap.free(map->h, map->flags); \ + map->heap.free(map->h, map->vals); \ + map->heap.free(map->h, map); + +#define MAP_RESET(map) \ + if (map && map->flags) { \ + memset(map->flags, 0xaa, __ac_fsize(map->n_buckets) * sizeof(int32)); \ + map->size = map->n_occupied = 0; \ + } + +#define MAP_GET(map, key, hashfunc, equalfunc) \ + if (map->n_buckets) { \ + int32 k = 0; \ + int32 i = 0; \ + int32 last = 0; \ + int32 mask = 0; \ + int32 step = 0; \ + mask = map->n_buckets - 1; \ + k = hashfunc(key); \ + i = k & mask; \ + last = i; \ + while (!__ac_isempty(map->flags, i) && \ + (__ac_isdel(map->flags, i) || !equalfunc(map->keys[i], key))) { \ + i = (i + (++step)) & mask; \ + if (i == last) \ + return map->n_buckets; \ } \ - SCOPE void map·destroy_##name(map·##name##_t* h) \ - { \ - if (h) \ - { \ - free((void*)h->keys); \ - free(h->flags); \ - free((void*)h->vals); \ - free(h); \ + return __ac_iseither(map->flags, i) ? map->n_buckets : i; \ + } else \ + return 0; + +#define MAP_GROW(map, key_t, val_t, new_n_buckets, hashfunc) \ + int32 *new_flags = nil; \ + int32 j = 1; \ + { \ + _roundup32(new_n_buckets); \ + if (new_n_buckets < 4) \ + new_n_buckets = 4; \ + if (map->size >= (int32)(new_n_buckets * __ac_HASH_UPPER + 0.5)) \ + j = 0; \ + else { \ + new_flags = \ + map->heap.alloc(map->h, __ac_fsize(new_n_buckets), sizeof(int32)); \ + if (!new_flags) return -1; \ + memset(new_flags, 0xaa, __ac_fsize(new_n_buckets) * sizeof(int32)); \ + if (map->n_buckets < new_n_buckets) { /* expand */ \ + key_t *new_keys = \ + map->heap.alloc(map->h, new_n_buckets, sizeof(key_t)); \ + if (!new_keys) { \ + map->heap.free(map->h, new_flags); \ + return -1; \ } \ - } \ - SCOPE void map·clear_##name(map·##name##_t* h) \ - { \ - if (h && h->flags) \ - { \ - memset(h->flags, 0xaa, __ac_fsize(h->n_buckets) * sizeof(int32)); \ - h->size = h->n_occupied = 0; \ + map->heap.free(map->h, map->keys); \ + map->keys = new_keys; \ + val_t *new_vals = \ + map->heap.alloc(map->h, new_n_buckets, sizeof(val_t)); \ + if (!new_vals) { \ + map->heap.free(map->h, new_flags); \ + return -1; \ } \ + map->heap.free(map->h, map->vals); \ + map->vals = new_vals; \ + } \ } \ - SCOPE int32 map·get_##name(const map·##name##_t* h, mapkey_t key) \ - { \ - if (h->n_buckets) \ - { \ - int32 k = 0; \ - int32 i = 0; \ - int32 last = 0; \ - int32 mask = 0; \ - int32 step = 0; \ - mask = h->n_buckets - 1; \ - k = __hash_func(key); \ - i = k & mask; \ - last = i; \ - while ( \ - !__ac_isempty(h->flags, i) && \ - (__ac_isdel(h->flags, i) || !__hash_equal(h->keys[i], key))) \ + } \ + if (j) { /* rehashing is needed */ \ + for (j = 0; j != map->n_buckets; ++j) { \ + if (__ac_iseither(map->flags, j) == 0) { \ + key_t key = map->keys[j]; \ + val_t val; \ + int32 new_mask; \ + new_mask = new_n_buckets - 1; \ + val = map->vals[j]; \ + __ac_set_isdel_true(map->flags, j); \ + while (1) { \ + int32 k, i, step = 0; \ + k = hashfunc(key); \ + i = k & new_mask; \ + while (!__ac_isempty(new_flags, i)) \ + i = (i + (++step)) & new_mask; \ + __ac_set_isempty_false(new_flags, i); \ + if (i < map->n_buckets && __ac_iseither(map->flags, i) == 0) { \ { \ - i = (i + (++step)) & mask; \ - if (i == last) return h->n_buckets; \ - } \ - return __ac_iseither(h->flags, i) ? h->n_buckets : i; \ - } else \ - return 0; \ - } \ - SCOPE int map·resize_##name(map·##name##_t* h, int32 new_n_buckets) \ - { /* This function uses 0.25*n_buckets bytes of working space instead of \ - [sizeof(key_t+val·t)+.25]*n_buckets. */ \ - int32* new_flags = 0; \ - int32 j = 1; \ - { \ - _roundup32(new_n_buckets); \ - if (new_n_buckets < 4) new_n_buckets = 4; \ - if (h->size >= (int32)(new_n_buckets * __ac_HASH_UPPER + 0.5)) \ - j = 0; /* requested size is too small */ \ - else \ - { /* hash table size to be changed (shrink or expand); rehash */ \ - new_flags = \ - (int32*)malloc(__ac_fsize(new_n_buckets) * sizeof(int32)); \ - if (!new_flags) return -1; \ - memset( \ - new_flags, 0xaa, __ac_fsize(new_n_buckets) * sizeof(int32)); \ - if (h->n_buckets < new_n_buckets) \ - { /* expand */ \ - mapkey_t* new_keys = (mapkey_t*)realloc( \ - (void*)h->keys, new_n_buckets * sizeof(mapkey_t)); \ - if (!new_keys) \ - { \ - free(new_flags); \ - return -1; \ - } \ - h->keys = new_keys; \ - if (map·is·map) \ - { \ - mapval·t* new_vals = (mapval·t*)realloc( \ - (void*)h->vals, new_n_buckets * sizeof(mapval·t)); \ - if (!new_vals) \ - { \ - free(new_flags); \ - return -1; \ - } \ - h->vals = new_vals; \ - } \ - } /* otherwise shrink */ \ + key_t tmp = map->keys[i]; \ + map->keys[i] = key; \ + key = tmp; \ } \ - } \ - if (j) \ - { /* rehashing is needed */ \ - for (j = 0; j != h->n_buckets; ++j) \ { \ - if (__ac_iseither(h->flags, j) == 0) \ - { \ - mapkey_t key = h->keys[j]; \ - mapval·t val; \ - int32 new_mask; \ - new_mask = new_n_buckets - 1; \ - if (map·is·map) val = h->vals[j]; \ - __ac_set_isdel_true(h->flags, j); \ - while (1) \ - { /* kick-out process; sort of like in Cuckoo hashing */ \ - int32 k, i, step = 0; \ - k = __hash_func(key); \ - i = k & new_mask; \ - while (!__ac_isempty(new_flags, i)) \ - i = (i + (++step)) & new_mask; \ - __ac_set_isempty_false(new_flags, i); \ - if (i < h->n_buckets && \ - __ac_iseither(h->flags, i) == 0) \ - { /* kick out the existing element */ \ - { \ - mapkey_t tmp = h->keys[i]; \ - h->keys[i] = key; \ - key = tmp; \ - } \ - if (map·is·map) \ - { \ - mapval·t tmp = h->vals[i]; \ - h->vals[i] = val; \ - val = tmp; \ - } \ - __ac_set_isdel_true(h->flags, \ - i); /* mark it as deleted in \ - the old hash table */ \ - } else \ - { /* write the element and jump out of the loop */ \ - h->keys[i] = key; \ - if (map·is·map) h->vals[i] = val; \ - break; \ - } \ - } \ - } \ + val_t tmp = map->vals[i]; \ + map->vals[i] = val; \ + val = tmp; \ } \ - if (h->n_buckets > new_n_buckets) \ - { /* shrink the hash table */ \ - h->keys = (mapkey_t*)realloc((void*)h->keys, \ - new_n_buckets * sizeof(mapkey_t)); \ - if (map·is·map) \ - h->vals = (mapval·t*)realloc( \ - (void*)h->vals, new_n_buckets * sizeof(mapval·t)); \ - } \ - free(h->flags); /* free the working space */ \ - h->flags = new_flags; \ - h->n_buckets = new_n_buckets; \ - h->n_occupied = h->size; \ - h->upper_bound = (int32)(h->n_buckets * __ac_HASH_UPPER + 0.5); \ + __ac_set_isdel_true(map->flags, i); \ + } else { \ + map->keys[i] = key; \ + map->vals[i] = val; \ + break; \ + } \ } \ - return 0; \ + } \ } \ - SCOPE int32 map·put_##name(map·##name##_t* h, mapkey_t key, int* ret) \ - { \ - int32 x = 0; \ - if (h->n_occupied >= h->upper_bound) \ - { /* update the hash table */ \ - if (h->n_buckets > (h->size << 1)) \ - { \ - if (map·resize_##name(h, h->n_buckets - 1) < 0) \ - { /* clear "deleted" elements */ \ - *ret = -1; \ - return h->n_buckets; \ - } \ - } else if (map·resize_##name(h, h->n_buckets + 1) < 0) \ - { /* expand the hash table */ \ - *ret = -1; \ - return h->n_buckets; \ - } \ - } /* TODO: to implement automatically shrinking; resize() already \ - support shrinking */ \ - { \ - int32 k, i, site, last, mask = h->n_buckets - 1, step = 0; \ - x = site = h->n_buckets; \ - k = __hash_func(key); \ - i = k & mask; \ - if (__ac_isempty(h->flags, i)) \ - x = i; /* for speed up */ \ - else \ - { \ - last = i; \ - while ( \ - !__ac_isempty(h->flags, i) && \ - (__ac_isdel(h->flags, i) || !__hash_equal(h->keys[i], key))) \ - { \ - if (__ac_isdel(h->flags, i)) site = i; \ - i = (i + (++step)) & mask; \ - if (i == last) \ - { \ - x = site; \ - break; \ - } \ - } \ - if (x == h->n_buckets) \ - { \ - if (__ac_isempty(h->flags, i) && site != h->n_buckets) \ - x = site; \ - else \ - x = i; \ - } \ - } \ - } \ - if (__ac_isempty(h->flags, x)) \ - { /* not present at all */ \ - h->keys[x] = key; \ - __ac_set_isboth_false(h->flags, x); \ - ++h->size; \ - ++h->n_occupied; \ - *ret = 1; \ - } else if (__ac_isdel(h->flags, x)) \ - { /* deleted */ \ - h->keys[x] = key; \ - __ac_set_isboth_false(h->flags, x); \ - ++h->size; \ - *ret = 2; \ - } else \ - *ret = 0; /* Don't touch h->keys[x] if present and not deleted */ \ - return x; \ + if (map->n_buckets > new_n_buckets) { /* shrink the hash table */ \ + key_t *new_keys = map->heap.alloc(map->h, new_n_buckets, sizeof(key_t)); \ + map->heap.free(map->h, map->keys); \ + map->keys = new_keys; \ + \ + val_t *new_vals = map->heap.alloc(map->h, new_n_buckets, sizeof(val_t)); \ + map->heap.free(map->h, map->vals); \ + map->vals = new_vals; \ } \ - SCOPE void map·del_##name(map·##name##_t* h, int32 x) \ - { \ - if (x != h->n_buckets && !__ac_iseither(h->flags, x)) \ - { \ - __ac_set_isdel_true(h->flags, x); \ - --h->size; \ + map->heap.free(map->h, map->flags); /* free the working space */ \ + map->flags = new_flags; \ + map->n_buckets = new_n_buckets; \ + map->n_occupied = map->size; \ + map->upper_bound = (int32)(map->n_buckets * __ac_HASH_UPPER + 0.5); \ + } \ + return 0; + +#define MAP_PUT(map, key, val, hashfunc, equalfunc, resizefunc, err) \ + { \ + int32 x = 0; \ + if (map->n_occupied >= map->upper_bound) { \ + if (map->n_buckets > (map->size << 1)) { \ + if (resizefunc(map, map->n_buckets - 1) < 0) { \ + *err = -1; \ + return map->n_buckets; \ } \ - } - -#define HASH_DECLARE(name, mapkey_t, mapval·t) \ - __HASH_TYPE(name, mapkey_t, mapval·t) \ - __HASH_PROTOTYPES(name, mapkey_t, mapval·t) - -#define MAP·MAKE2( \ - name, SCOPE, mapkey_t, mapval·t, map·is·map, __hash_func, __hash_equal) \ - __HASH_TYPE(name, mapkey_t, mapval·t) \ - __HASH_IMPL( \ - name, SCOPE, mapkey_t, mapval·t, map·is·map, __hash_func, __hash_equal) - -#define MAP·MAKE( \ - name, mapkey_t, mapval·t, map·is·map, __hash_func, __hash_equal) \ - MAP·MAKE2(name, \ - static map·inline klib_unused, \ - mapkey_t, \ - mapval·t, \ - map·is·map, \ - __hash_func, \ - __hash_equal) - -/* --- BEGIN OF HASH FUNCTIONS --- */ - -/*! @function - @abstract Integer hash function - @param key The integer [mapint32_t] - @return The hash value [mapint_t] - */ -#define map·int_hash_func(key) (int32)(key) -/*! @function - @abstract Integer comparison function - */ -#define map·int_hash_equal(a, b) ((a) == (b)) -/*! @function - @abstract 64-bit integer hash function - @param key The integer [mapint64_t] - @return The hash value [mapint_t] - */ -#define map·int64_hash_func(key) (int32)((key) >> 33 ^ (key) ^ (key) << 11) -/*! @function - @abstract 64-bit integer comparison function - */ -#define map·int64_hash_equal(a, b) ((a) == (b)) -/*! @function - @abstract const char* hash function - @param s Pointer to a null terminated string - @return The hash value - */ -static map·inline int32 - __ac_X31_hash_string(const char* s) -{ - int32 h = (int32)(*s); - if (h != 0) { - for (const char* it = s; *it; ++it) - h = (h << 5) - h + (int32)*it; - } - return h; -} -/*! @function - @abstract Another interface to const char* hash function - @param key Pointer to a null terminated string [const char*] - @return The hash value [mapint_t] - */ -#define map·str_hash_func(key) __ac_X31_hash_string(key) -/*! @function - @abstract Const char* comparison function - */ -#define map·str_hash_equal(a, b) (strcmp(a, b) == 0) - -static inline int32 -__ac_Wang_hash(int32 key) -{ - key += ~(key << 15); - key ^= (key >> 10); - key += (key << 3); - key ^= (key >> 6); - key += ~(key << 11); - key ^= (key >> 16); - return key; -} -#define map·int_hash_func2(key) __ac_Wang_hash((int32)key) - -/* --- END OF HASH FUNCTIONS --- */ - -/* Other convenient macros... */ - -/*! - @abstract Type of the hash table. - @param name Name of the hash table [symbol] - */ -#define map·t(name) map·##name##_t - -/*! @function - @abstract Initiate a hash table. - @param name Name of the hash table [symbol] - @return Pointer to the hash table [khash_t(name)*] - */ -#define map·init(name) map·init·##name() - -/*! @function - @abstract Destroy a hash table. - @param name Name of the hash table [symbol] - @param h Pointer to the hash table [khash_t(name)*] - */ -#define map·destroy(name, h) map·destroy_##name(h) - -/*! @function - @abstract Reset a hash table without deallocating memory. - @param name Name of the hash table [symbol] - @param h Pointer to the hash table [khash_t(name)*] - */ -#define map·clear(name, h) map·clear_##name(h) - -/*! @function - @abstract Resize a hash table. - @param name Name of the hash table [symbol] - @param h Pointer to the hash table [khash_t(name)*] - @param s New size [mapint_t] - */ -#define map·resize(name, h, s) map·resize_##name(h, s) - -/*! @function - @abstract Insert a key to the hash table. - @param name Name of the hash table [symbol] - @param h Pointer to the hash table [khash_t(name)*] - @param k Key [type of keys] - @param r Extra return code: -1 if the operation failed; - 0 if the key is present in the hash table; - 1 if the bucket is empty (never used); 2 if the element in - the bucket has been deleted [int*] - @return Iterator to the inserted element [mapint_t] - */ -#define map·put(name, h, k, r) map·put_##name(h, k, r) - -/*! @function - @abstract Retrieve a key from the hash table. - @param name Name of the hash table [symbol] - @param h Pointer to the hash table [khash_t(name)*] - @param k Key [type of keys] - @return Iterator to the found element, or map·end(h) if the element is - absent [mapint_t] - */ -#define map·get(name, h, k) map·get_##name(h, k) - -/*! @function - @abstract Remove a key from the hash table. - @param name Name of the hash table [symbol] - @param h Pointer to the hash table [khash_t(name)*] - @param k Iterator to the element to be deleted [mapint_t] - */ -#define map·del(name, h, k) map·del_##name(h, k) - -/*! @function - @abstract Test whether a bucket contains data. - @param h Pointer to the hash table [khash_t(name)*] - @param x Iterator to the bucket [mapint_t] - @return 1 if containing data; 0 otherwise [int] - */ -#define map·exist(h, x) (!__ac_iseither((h)->flags, (x))) - -/*! @function - @abstract Get key given an iterator - @param h Pointer to the hash table [khash_t(name)*] - @param x Iterator to the bucket [mapint_t] - @return Key [type of keys] - */ -#define map·key(h, x) ((h)->keys[x]) - -/*! @function - @abstract Get value given an iterator - @param h Pointer to the hash table [khash_t(name)*] - @param x Iterator to the bucket [mapint_t] - @return Value [type of values] - @discussion For hash sets, calling this results in segfault. - */ -#define map·val(h, x) ((h)->vals[x]) - -/*! @function - @abstract Alias of map·val() - */ -#define map·value(h, x) ((h)->vals[x]) - -/*! @function - @abstract Get the start iterator - @param h Pointer to the hash table [khash_t(name)*] - @return The start iterator [mapint_t] - */ -#define map·begin(h) (int32)(0) - -/*! @function - @abstract Get the end iterator - @param h Pointer to the hash table [khash_t(name)*] - @return The end iterator [mapint_t] - */ -#define map·end(h) ((h)->n_buckets) - -/*! @function - @abstract Get the number of elements in the hash table - @param h Pointer to the hash table [khash_t(name)*] - @return Number of elements in the hash table [mapint_t] - */ -#define map·len(h) ((h)->size) - -/*! @function - @abstract Get the number of buckets in the hash table - @param h Pointer to the hash table [khash_t(name)*] - @return Number of buckets in the hash table [mapint_t] - */ -#define map·n_buckets(h) ((h)->n_buckets) - -/*! @function - @abstract Iterate over the entries in the hash table - @param h Pointer to the hash table [khash_t(name)*] - @param kvar Variable to which key will be assigned - @param vvar Variable to which value will be assigned - @param code Block of code to execute - */ -#define map·foreach(h, kvar, vvar, code) \ - { \ - int32 __i; \ - for (__i = map·begin(h); __i != map·end(h); ++__i) \ - { \ - if (!map·exist(h, __i)) continue; \ - (kvar) = map·key(h, __i); \ - (vvar) = map·val(h, __i); \ - code; \ - } \ - } - -/*! @function - @abstract Iterate over the values in the hash table - @param h Pointer to the hash table [khash_t(name)*] - @param vvar Variable to which value will be assigned - @param code Block of code to execute - */ -#define map·foreach_val(h, vvar, code) \ - { \ - int32 __i; \ - for (__i = map·begin(h); __i != map·end(h); ++__i) \ - { \ - if (!map·exist(h, __i)) continue; \ - (vvar) = map·val(h, __i); \ - code; \ - } \ - } - -/*! @function - @abstract Iterate over the keys in the hash table - @param h Pointer to the hash table [khash_t(name)*] - @param kvar Variable to which key will be assigned - @param code Block of code to execute - */ -#define map·foreach_key(h, kvar, code) \ + } else if (resizefunc(map, map->n_buckets + 1) < 0) { \ + *err = -1; \ + return map->n_buckets; \ + } \ + } \ + \ { \ - int32 __i; \ - for (__i = map·begin(h); __i != map·end(h); ++__i) \ - { \ - if (!map·exist(h, __i)) continue; \ - (kvar) = map·key(h, __i); \ - code; \ + int32 k, i, site, last, mask = map->n_buckets - 1, step = 0; \ + x = site = map->n_buckets; \ + k = hashfunc(key); \ + i = k & mask; \ + if (__ac_isempty(map->flags, i)) \ + x = i; /* for speed up */ \ + else { \ + last = i; \ + while (!__ac_isempty(map->flags, i) && \ + (__ac_isdel(map->flags, i) || !equalfunc(map->keys[i], key))) { \ + if (__ac_isdel(map->flags, i)) \ + site = i; \ + i = (i + (++step)) & mask; \ + if (i == last) { \ + x = site; \ + break; \ + } \ } \ - } - -/*! @function - @abstract Iterate over the values in the hash table - @param h Pointer to the hash table [khash_t(name)*] - @param code Block of code to execute for free operation. - */ -#define map·free_all(h, code) \ - { \ - int32 __i; \ - for (__i = map·begin(h); __i != map·end(h); ++__i) \ - { \ - if (!map·exist(h, __i)) continue; \ - code(map·val(h, __i)); \ + if (x == map->n_buckets) { \ + if (__ac_isempty(map->flags, i) && site != map->n_buckets) \ + x = site; \ + else \ + x = i; \ } \ - } - -/* More convenient interfaces */ - -/*! @function - @abstract Instantiate a hash set containing integer keys - @param name Name of the hash table [symbol] - */ -#define define_int32_set(name) \ - MAP·MAKE(name, int32, char, 0, map·int_hash_func, map·int_hash_equal) - -/*! @function - @abstract Instantiate a hash map containing integer keys - @param name Name of the hash table [symbol] - @param mapval·t Type of values [type] - */ -#define define_int32_map(name, mapval·t) \ - MAP·MAKE(name, int32, mapval·t, 1, map·int_hash_func, map·int_hash_equal) - -/*! @function - @abstract Instantiate a hash map containing integer keys - @param name Name of the hash table [symbol] - @param mapval·t Type of values [type] - */ -#define define_uint32_set(name) \ - MAP·MAKE(name, uint32, char, 0, map·int_hash_func, map·int_hash_equal) - -/*! @function - @abstract Instantiate a hash map containing integer keys - @param name Name of the hash table [symbol] - @param mapval·t Type of values [type] - */ -#define define_uint32_map(name, mapval·t) \ - MAP·MAKE(name, uint32, mapval·t, 1, map·int_hash_func, map·int_hash_equal) - -/*! @function - @abstract Instantiate a hash set containing 64-bit integer keys - @param name Name of the hash table [symbol] - */ -#define define_int64_set(name) \ - MAP·MAKE(name, int64, char, 0, map·int64_hash_func, map·int64_hash_equal) - -/*! @function - @abstract Instantiate a hash map containing 64-bit integer keys - @param name Name of the hash table [symbol] - @param mapval·t Type of values [type] - */ -#define define_int64_map(name, mapval·t) \ - MAP·MAKE(name, int64, mapval·t, 1, map·int64_hash_func, map·int64_hash_equal) - -/*! @function - @abstract Instantiate a hash set containing 64-bit unsigned integer keys - @param name Name of the hash table [symbol] - */ -#define define_uint64_set(name) \ - MAP·MAKE(name, uint64, char, 0, map·int64_hash_func, map·int64_hash_equal) - -/*! @function - @abstract Instantiate a hash map containing 64-bit unsigned integer keys - @param name Name of the hash table [symbol] - @param mapval·t Type of values [type] - */ -#define define_uint64_map(name, mapval·t) \ - MAP·MAKE(name, uint64, mapval·t, 1, map·int64_hash_func, map·int64_hash_equal) - -/*! @function - @abstract Instantiate a hash map containing const char* keys - @param name Name of the hash table [symbol] - */ -#define define_str_hashset(name) \ - MAP·MAKE(name, const byte*, byte, 0, map·str_hash_func, map·str_hash_equal) - -/*! @function - @abstract Instantiate a hash map containing const char* keys - @param name Name of the hash table [symbol] - @param mapval·t Type of values [type] - */ -#define define_str·map(name, mapval·t) \ - MAP·MAKE(name, const byte*, mapval·t, 1, map·str_hash_func, map·str_hash_equal) + } \ + } \ + if (__ac_isempty(map->flags, x)) { /* not present at all */ \ + map->keys[x] = key; \ + __ac_set_isboth_false(map->flags, x); \ + ++map->size; \ + ++map->n_occupied; \ + *err = 1; \ + } else if (__ac_isdel(map->flags, x)) { /* deleted */ \ + map->keys[x] = key; \ + __ac_set_isboth_false(map->flags, x); \ + ++map->size; \ + } else \ + *err = 0; \ + return x; \ + } + +#define MAP_DEL(map, x) \ + { \ + if (x != map->n_buckets && !__ac_iseither(map->flags, x)) { \ + __ac_set_isdel_true(map->flags, x); \ + --map->size; \ + } \ + } + +static int32 hash_string(byte *s) { + int32 h; + byte *it; + + h = (int32)(*s); + if (h != 0) { + for (it = s; *it; ++it) + h = (h << 5) - h + (int32)*it; + } + return h; +} diff --git a/include/libn/macro/qsort.h b/include/libn/macro/qsort.h index 8daf0fa..cc38bca 100644 --- a/include/libn/macro/qsort.h +++ b/include/libn/macro/qsort.h @@ -1,181 +1,88 @@ -/* - * Copyright (c) 2013, 2017 Alexey Tourbin - * - * Permission is hereby granted, free of charge, to any person obtaining a copy - * of this software and associated documentation files (the "Software"), to deal - * in the Software without restriction, including without limitation the rights - * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - * copies of the Software, and to permit persons to whom the Software is - * furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included in - * all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - * SOFTWARE. - */ +#pragma once -/* - * This is a traditional Quicksort implementation which mostly follows - * [Sedgewick 1978]. Sorting is performed entirely on array indices, - * while actual access to the array elements is abstracted out with the - * user-defined `LESS` and `SWAP` primitives. +/* + * Straight implementation of Sedgewick's median qsort + * + * @LEN: name of parameter length + * @QLESS: name of function that computes array[i] < array[j] + * should return a boolean + * @QSWAP: name of function that swaps array[i] <-> array[j] * - * Synopsis: - * QSORT(N, LESS, SWAP); - * where - * N - the number of elements in A[]; - * LESS(i, j) - compares A[i] to A[j]; - * SWAP(i, j) - exchanges A[i] with A[j]. + * NOTE: This can perform on strided arrays. + * Make sure to use parens liberally to ensure hygeine! */ -#pragma once - -/* Sort 3 elements. */ -#define Q_SORT3(q_a1, q_a2, q_a3, Q_LESS, Q_SWAP) \ -do { \ - if (Q_LESS(q_a2, q_a1)) { \ - if (Q_LESS(q_a3, q_a2)) \ - Q_SWAP(q_a1, q_a3); \ - else { \ - Q_SWAP(q_a1, q_a2); \ - if (Q_LESS(q_a3, q_a2)) \ - Q_SWAP(q_a2, q_a3); \ - } \ - } \ - else if (Q_LESS(q_a3, q_a2)) { \ - Q_SWAP(q_a2, q_a3); \ - if (Q_LESS(q_a2, q_a1)) \ - Q_SWAP(q_a1, q_a2); \ - } \ -} while (0) - -/* Partition [q_l,q_r] around a pivot. After partitioning, - * [q_l,q_j] are the elements that are less than or equal to the pivot, - * while [q_i,q_r] are the elements greater than or equal to the pivot. */ -#define Q_PARTITION(q_l, q_r, q_i, q_j, Q_UINT, Q_LESS, Q_SWAP) \ -do { \ - /* The middle element, not to be confused with the median. */ \ - Q_UINT q_m = q_l + ((q_r - q_l) >> 1); \ - /* Reorder the second, the middle, and the last items. \ - * As [Edelkamp Weiss 2016] explain, using the second element \ - * instead of the first one helps avoid bad behaviour for \ - * decreasingly sorted arrays. This method is used in recent \ - * versions of gcc's std::sort, see gcc bug 58437#c13, although \ - * the details are somewhat different (cf. #c14). */ \ - Q_SORT3(q_l + 1, q_m, q_r, Q_LESS, Q_SWAP); \ - /* Place the median at the beginning. */ \ - Q_SWAP(q_l, q_m); \ - /* Partition [q_l+2, q_r-1] around the median which is in q_l. \ - * q_i and q_j are initially off by one, they get decremented \ - * in the do-while loops. */ \ - q_i = q_l + 1; q_j = q_r; \ - while (1) { \ - do q_i++; while (Q_LESS(q_i, q_l)); \ - do q_j--; while (Q_LESS(q_l, q_j)); \ - if (q_i >= q_j) break; /* Sedgewick says "until j < i" */ \ - Q_SWAP(q_i, q_j); \ - } \ - /* Compensate for the i==j case. */ \ - q_i = q_j + 1; \ - /* Put the median to its final place. */ \ - Q_SWAP(q_l, q_j); \ - /* The median is not part of the left subfile. */ \ - q_j--; \ -} while (0) - -/* Insertion sort is applied to small subfiles - this is contrary to - * Sedgewick's suggestion to run a separate insertion sort pass after - * the partitioning is done. The reason I don't like a separate pass - * is that it triggers extra comparisons, because it can't see that the - * medians are already in their final positions and need not be rechecked. - * Since I do not assume that comparisons are cheap, I also do not try - * to eliminate the (q_j > q_l) boundary check. */ -#define Q_INSERTION_SORT(q_l, q_r, Q_UINT, Q_LESS, Q_SWAP) \ -do { \ - Q_UINT q_i, q_j; \ - /* For each item starting with the second... */ \ - for (q_i = q_l + 1; q_i <= q_r; q_i++) \ - /* move it down the array so that the first part is sorted. */ \ - for (q_j = q_i; q_j > q_l && (Q_LESS(q_j, q_j - 1)); q_j--) \ - Q_SWAP(q_j, q_j - 1); \ -} while (0) - -/* When the size of [q_l,q_r], i.e. q_r-q_l+1, is greater than or equal to - * Q_THRESH, the algorithm performs recursive partitioning. When the size - * drops below Q_THRESH, the algorithm switches to insertion sort. - * The minimum valid value is probably 5 (with 5 items, the second and - * the middle items, the middle itself being rounded down, are distinct). */ -#define Q_THRESH 16 - -/* The main loop. */ -#define Q_LOOP(Q_UINT, Q_N, Q_LESS, Q_SWAP) \ -do { \ - Q_UINT q_l = 0; \ - Q_UINT q_r = (Q_N) - 1; \ - Q_UINT q_sp = 0; /* the number of frames pushed to the stack */ \ - struct { Q_UINT q_l, q_r; } \ - /* On 32-bit platforms, to sort a "char[3GB+]" array, \ - * it may take full 32 stack frames. On 64-bit CPUs, \ - * though, the address space is limited to 48 bits. \ - * The usage is further reduced if Q_N has a 32-bit type. */ \ - q_st[sizeof(Q_UINT) > 4 && sizeof(Q_N) > 4 ? 48 : 32]; \ - while (1) { \ - if (q_r - q_l + 1 >= Q_THRESH) { \ - Q_UINT q_i, q_j; \ - Q_PARTITION(q_l, q_r, q_i, q_j, Q_UINT, Q_LESS, Q_SWAP); \ - /* Now have two subfiles: [q_l,q_j] and [q_i,q_r]. \ - * Dealing with them depends on which one is bigger. */ \ - if (q_j - q_l >= q_r - q_i) \ - Q_SUBFILES(q_l, q_j, q_i, q_r); \ - else \ - Q_SUBFILES(q_i, q_r, q_l, q_j); \ - } \ - else { \ - Q_INSERTION_SORT(q_l, q_r, Q_UINT, Q_LESS, Q_SWAP); \ - /* Pop subfiles from the stack, until it gets empty. */ \ - if (q_sp == 0) break; \ - q_sp--; \ - q_l = q_st[q_sp].q_l; \ - q_r = q_st[q_sp].q_r; \ - } \ - } \ -} while (0) - -/* The missing part: dealing with subfiles. - * Assumes that the first subfile is not smaller than the second. */ -#define Q_SUBFILES(q_l1, q_r1, q_l2, q_r2) \ -do { \ - /* If the second subfile is only a single element, it needs \ - * no further processing. The first subfile will be processed \ - * on the next iteration (both subfiles cannot be only a single \ - * element, due to Q_THRESH). */ \ - if (q_l2 == q_r2) { \ - q_l = q_l1; \ - q_r = q_r1; \ - } \ - else { \ - /* Otherwise, both subfiles need processing. \ - * Push the larger subfile onto the stack. */ \ - q_st[q_sp].q_l = q_l1; \ - q_st[q_sp].q_r = q_r1; \ - q_sp++; \ - /* Process the smaller subfile on the next iteration. */ \ - q_l = q_l2; \ - q_r = q_r2; \ - } \ -} while (0) - -/* And now, ladies and gentlemen, may I proudly present to you... */ -#define QSORT(Q_N, Q_LESS, Q_SWAP) \ -do { \ - if ((Q_N) > 1) \ - /* We could check sizeof(Q_N) and use "unsigned", but at least \ - * on x86_64, this has the performance penalty of up to 5%. */ \ - Q_LOOP(ulong, Q_N, Q_LESS, Q_SWAP); \ -} while (0) +#define QSORT(LEN, QLESS, QSWAP) \ + int i, j, m, r, l; \ + struct { \ + int i, j; \ + } *sp, base[48]; \ + sp = base; \ + \ + l = 1; \ + r = LEN-1; \ + \ + if (LEN <= 14) goto ENDOUTER; \ +OUTERLOOP: \ + SWAP((l+r)/2, l+1); \ + \ + if (QLESS(r, l+1)) QSWAP(r, l+1); \ + if (QLESS(r, l)) QSWAP(r, l); \ + if (QLESS(l, l+1)) QSWAP(l, l+1); \ + \ + i = l+1; \ + j = r; \ + \ +INNERLOOP: \ + do ++i; while(QLESS(i, l)); \ + do --j; while(QLESS(l, j)); \ + if (j < i) goto ENDINNER; \ + QSWAP(i, j); \ + goto INNERLOOP; \ + \ +ENDINNER: \ + QSWAP(l, j); \ + \ + if (MAX(j-l, r-i+1) <= 14) { \ + if (sp == base) { \ + goto ENDOUTER; \ + } else { \ + l = sp[-1].i; \ + r = sp[-1].j; \ + sp--; \ + } \ + } else { \ + if (MIN(j-l, r-i+1) <= 14) { \ + if (j-l > r-i+1) { \ + l = l; \ + r = j-1; \ + } else { \ + l = i; \ + r = r; \ + } \ + } else { \ + if (j-l > r-i+1) { \ + sp->i = l; \ + sp->j = j-1; \ + sp++; \ + \ + l = i; \ + r = r; \ + } else { \ + sp->i = i; \ + sp->j = r; \ + sp++; \ + \ + l = l; \ + r = j-1; \ + } \ + } \ + } \ + goto OUTERLOOP; \ + \ +ENDOUTER: \ + for (i = 1; i < LEN; i++) { \ + for (j = i; j > 0 && QLESS(j, j-1); j--) { \ + QSWAP(j, j-1); \ + } \ + } \ diff --git a/rules.mk b/rules.mk index 5260157..3f2de37 100644 --- a/rules.mk +++ b/rules.mk @@ -39,6 +39,7 @@ targets: $(LIBS) $(BINS) .PHONY: database database: $(LIBS) $(BINS) gentags + ctags -R sys .PHONY: clean clean: diff --git a/sys/libbio/align.c b/sys/libbio/align.c index 8d1c119..0f00fce 100644 --- a/sys/libbio/align.c +++ b/sys/libbio/align.c @@ -87,16 +87,24 @@ sortpos(uintptr len, uint64 vals[], int locs[]) #undef SWAP } +/* + * sketch + * @param seq: '0' terminated string + * @param len: number of sequential sketches to keep + * @param vals: buffer to store hashes of sketch. + * @param locs: buffer to store location of sketch hashes + */ error -aln·sketch(byte *seq, int L, uint64 *vals[aln·N], int *locs[aln·N]) +aln·sketch(byte *seq, int len, uint64 *vals[aln·N], int *locs[aln·N]) { int i, n, l, *loc; - uint64 kmer, tmp[4], h[3], *val; + uint64 kmer, h[3], *val; + int tmpi[2]; + uint64 tmpu[2]; for (n = 0; n < aln·N; n++) { - val = vals[n]; - for (l = 0; l < L; l++) { - val[l] = UINT64_MAX; + for (l = 0; l < len; l++) { + vals[n][l] = UINT64_MAX; } } @@ -112,26 +120,21 @@ aln·sketch(byte *seq, int L, uint64 *vals[aln·N], int *locs[aln·N]) loc = locs[n]; h[2] = (h[0] + n * h[1]) & aln·mask; - for (i = 0; i < L && h[2] < val[i]; i++) { + for (i = 0; i < len && h[2] < val[i]; i++) { ; } - tmp[1] = h[2]; - tmp[3] = l; + tmpu[1] = h[2]; + tmpi[1] = l; for (i -= 1; i >= 0; i--) { - tmp[0] = tmp[1]; - tmp[1] = val[i]; - val[i] = tmp[0]; - - tmp[2] = tmp[3]; - tmp[3] = loc[i]; - loc[i] = tmp[2]; + tmpu[0] = tmpu[1], tmpu[1] = val[i], val[i] = tmpu[0]; + tmpi[0] = tmpi[1], tmpi[1] = loc[i], loc[i] = tmpi[0]; } } } for (n = 0; n < aln·N; n++) { - sortpos(L, vals[n], locs[n]); + sortpos(len, vals[n], locs[n]); } return 0; @@ -139,30 +142,37 @@ aln·sketch(byte *seq, int L, uint64 *vals[aln·N], int *locs[aln·N]) static int -cmparrs(int L, uint64 a[], uint64 b[]) +lessarrs(int len, uint64 a[], uint64 b[]) { int l; - for (l = 0; l < L; l++) { - if (a[l] < b[l]) return -1; + for (l = 0; l < len; l++) { + if (a[l] < b[l]) return 1; + if (a[l] > b[l]) return 0; } - return +1; + return 0; +} + +static +void +swaparrs(int len, uint64 a[], uint64 b[]) +{ + int l; + uint64 tmp; + + for (l = 0; l < len; l++) { + tmp = a[l], a[l] = b[l], b[l] = tmp; + } } -// TODO: replace with inlined version error -aln·sort(uintptr len, int L, uint64 *vals) +aln·sort(uintptr n, int len, uint64 *vals) { - uint64 tmp[128]; -#define LESS(i, j) (cmparrs((L), (vals+(i)*L), (vals+(j)*L))) -#define SWAP(i, j) (memcpy(tmp, (vals+(i)*L), (sizeof(uint64)*L)), \ - memcpy((vals+(i)*L), (vals+(j)*L), (sizeof(uint64)*L)), \ - memcpy((vals+(j)*L), tmp, (sizeof(uint64)*L))) - - QSORT(len, LESS, SWAP); +#define LESS(i, j) (lessarrs(len, vals+((i)*len), vals+((j)*len))) +#define SWAP(i, j) (swaparrs(len, vals+((i)*len), vals+((j)*len))) + QSORT(n, LESS, SWAP); #undef LESS #undef SWAP - // qsort(vals, len, l*sizeof(uint64), &compare); return 0; } diff --git a/sys/libbio/phylo.c b/sys/libbio/phylo.c index 1cc5b5d..1bb5bc7 100644 --- a/sys/libbio/phylo.c +++ b/sys/libbio/phylo.c @@ -1,5 +1,6 @@ #include #include +#include #include // ----------------------------------------------------------------------- @@ -50,6 +51,8 @@ FOUND: } else { prev->sibling = child->sibling; } + + parent->nchild--; return error·nil; } @@ -96,6 +99,55 @@ phylo·countleafs(bio·Node *node, int *n) return 0; } +// ----------------------------------------------------------------------- +// generic operations on tree + +void* +phylo·postorder(bio·Node *clade, void *(*op)(bio·Node*, void*), void *ctx) +{ + bio·Node *it; + + for (it = clade->child; it != nil; it = it->sibling) { + ctx = phylo·postorder(it, op, ctx); + } + + return op(clade, ctx); +} + +void* +phylo·preorder(bio·Node *clade, void *(*op)(bio·Node*, void*), void *ctx) +{ + bio·Node *it; + + ctx = op(clade, ctx); + for (it = clade->child; it != nil; it = it->sibling) { + ctx = phylo·preorder(it, op, ctx); + } + + return ctx; +} + +static +inline +void* +appendleaf(bio·Node *node, void* list) +{ + bio·Node **leafs; + + leafs = list; + if (!node->nchild) { + *leafs++ = node; + } + + return leafs; +} + +void +phylo·getleafs(bio·Tree tree, bio·Node **leafs) +{ + phylo·postorder(tree.root, &appendleaf, leafs); +} + // ----------------------------------------------------------------------- // tree editing @@ -147,64 +199,219 @@ phylo·ladderize(bio·Node *root) return 0; } +/* + * compute all distances from a given node + * must provide a working buffer + */ + +struct Tuple +{ + double *d; + bio·Node **n; +}; + +static +struct Tuple +getdistsfrom(bio·Node *node, bio·Node *prev, double curr, double *dist, bio·Node **list) +{ + bio·Node *it; + struct Tuple ret; + + *dist++ = curr; + *list++ = node; + + ret.d = dist; + ret.n = list; + + if (node->parent && node->parent != prev) { + ret = getdistsfrom(node->parent, node, curr + node->dist, dist, list); + + dist = ret.d; + list = ret.n; + } + + for (it = node->child; it != nil; it = it->sibling) { + if (it != prev) { + ret = getdistsfrom(it, node, curr + it->dist, dist, list); + + dist = ret.d; + list = ret.n; + } + } + + return ret; +} + +int +phylo·getdistsfrom(bio·Node *node, int len, double *dist, bio·Node **list) +{ + struct Tuple ret; + // TODO: Better bounds checking. + + ret = getdistsfrom(node, nil, 0.0, dist, list); + + assert(ret.n - list == len); + assert(ret.d - dist == len); + + return len; +} + +/* +static +void +disttoroot(bio·Node *clade, double anc, double *dists) +{ + double d; + bio·Node *it; + + *dists++ = anc + clade->dist; + d = dists[-1]; + for (it = clade->child; it != nil; it = it->sibling) { + disttoroot(it, d, ++dists); + } +} + +void +phylo·disttoroot(bio·Tree tree, double *dists) +{ + disttoroot(tree.root, 0.0, dists); +} +*/ + +/* + * compute the path constituting the tree diameter + * returns the number of edges in the path + */ + +static +void +sort·nodedists(uintptr len, double fs[], bio·Node* ns[]) +{ + double f; + bio·Node *n; +#define LESS(i, j) (fs[i] < fs[j]) +#define SWAP(i, j) (n = ns[i], f = fs[i], \ + fs[i] = fs[j], ns[i] = ns[j], \ + fs[j] = f, ns[j] = n) + QSORT(len, LESS, SWAP); +#undef LESS +#undef SWAP +} + +#define BUFLEN 4096 +double +phylo·diameter(bio·Tree tree, int *len, bio·Node **path) +{ + // TODO: deal with large tree > BUFLEN gracefully + int n; + double fbuf[BUFLEN]; + bio·Node *nbuf[BUFLEN]; + + n = tree.root->nnode; + + assert(n < BUFLEN); + + n = phylo·getdistsfrom(tree.root, tree.root->nnode, fbuf, nbuf); + sort·nodedists(n, fbuf, nbuf); + + path[0] = nbuf[n-1]; + printf("first end '%s'\n", path[0]->name); + + n = phylo·getdistsfrom(path[0], n, fbuf, nbuf); + sort·nodedists(n, fbuf, nbuf); + printf("second end '%s'\n", nbuf[n-1]->name); + + *len = 0; + + // TODO: Traverse up the tree from each node + // Find MRCA by intersection of nodes hit + + return 0.0; +} +#undef BUFLEN + +/* + * reroot a tree on a new node + */ static error rotateparent(bio·Node *node, bio·Node *to) { error err; - bio·Node *it, *prev; + // NOTE: will this ever be taken? if (node->parent == to) { return 0; } - prev = nil; - for (it = node->child; it != nil; it = it->sibling) { - if (it == to) goto FOUND; - prev = it; + if (!node->parent) { + goto RMCHILD; } - errorf("inconsistent topology: attempting to rotate to non-child"); - return 1; -FOUND: - if (err = rotateparent(node->parent, node), err) { + err = rotateparent(node->parent, node); + if (err) { errorf("failure: broken tree"); return err; } - if (!prev) { - node->child = node->parent; - } else { - prev->sibling = node->parent; + err = phylo·addchild(node, node->parent); + if (err) { + errorf("inconsistent topology: could not add parent '%s' as child of '%s'", node->parent->name, node->name); + return err; } - node->parent->sibling = it->sibling; - node->parent = to; +RMCHILD: + err = phylo·rmchild(node, to); + if (err) { + errorf("inconsistent topology: could not remove child '%s' from '%s'", to->name, node->name); + return err; + } + + node->parent = to; return 0; } -#define FLOAT_PREC .00000001 +#define PREC .00000001 error phylo·reroot(bio·Tree *tree, bio·Node *node, double d) { - bio·Node *old; bio·Node *new; // TODO: should check that node is part of this tree? + // TODO: should we check if node->parent != nil? - old = tree->root; - if (fabs(d) < FLOAT_PREC) new = node; - if (fabs(d-node->dist) < FLOAT_PREC) new = node->parent; - else { + if (fabs(d) < PREC) { + new = node; + rotateparent(node->parent, node); + } else if (fabs(d-node->dist) < PREC) { + new = node->parent; + if (new->parent->parent) { + rotateparent(new->parent->parent, new->parent); + } + } else { new = tree->heap.alloc(tree->h, 1, sizeof(*new)); memset(new, 0, sizeof(*new)); + + phylo·addchild(new, node); + node->parent = new; + + phylo·addchild(new, node->parent); + if (node->parent->parent) { + rotateparent(node->parent->parent, node->parent); + } + node->parent->parent = new; } - // TODO: Use rotateparent() to perform necessary flips + printf("number of children on old root: %d\n", tree->root->nchild); + tree->root = new; + tree->nleaf = 0; + + phylo·countleafs(new, &tree->nleaf); + phylo·countnodes(new, &new->nnode); - tree->root = new; return 0; } +#undef PREC // ----------------------------------------------------------------------- // ancestral inference @@ -216,7 +423,7 @@ struct phylo·InferOpts }; error -phylo·ancestralinfer(bio·Tree *tree, struct phylo·InferOpts opts) +phylo·inferancestral(bio·Tree *tree, struct phylo·InferOpts opts) { return 0; } diff --git a/sys/libbio/rules.mk b/sys/libbio/rules.mk index fd3e066..0e3d482 100644 --- a/sys/libbio/rules.mk +++ b/sys/libbio/rules.mk @@ -21,7 +21,7 @@ LIBS_$(d) := $(d)/libbio.a LIBS_$(d) := $(patsubst $(SRC_DIR)/%, $(OBJ_DIR)/%, $(LIBS_$(d))) LIBS := $(LIBS) $(LIBS_$(d)) -BINS_$(d) := $(d)/test +BINS_$(d) := $(d)/test $(d)/simulate BINS_$(d) := $(patsubst $(SRC_DIR)/%, $(OBJ_DIR)/%, $(BINS_$(d))) BINS := $(BINS) $(BINS_$(d)) diff --git a/sys/libbio/simulate.c b/sys/libbio/simulate.c new file mode 100644 index 0000000..a4567d5 --- /dev/null +++ b/sys/libbio/simulate.c @@ -0,0 +1,117 @@ +#include +#include +#include + +#define SEQLEN 2560 +static byte *SEQ = +"GGCGGCTTCGGTGCGCTGTGTGCATTGCCGCAAAAATATCGTGAACCCGTGCTGGTTTCCGGCACTGACGGCGTAGGTAC" +"CAAGCTGCGTCTGGCAATGGACTTAAAACGTCACGACACCATTGGTATTGATCTGGTCGCCATGTGCGTTAATGACCTGG" +"TGGTGCAAGGTGCGGAACCGCTGTTTTTCCTCGACTATTACGCAACCGGAAAACTGGATGTTGATACCGCTTCAGCGGTG" +"ATCAGCGGCATTGCGGAAGGTTGTCTGCAATCGGGCTGTTCTCTGGTGGGTGGCGAAACGGCAGAAATGCCGGGGATGTA" +"TCACGGTGAAGATTACGATGTCGCGGGTTTCTGCGTGGGCGTGGTAGAAAAATCAGAAATCATCGACGGCTCTAAAGTCA" +"GCGACGGCGATGTGCTGATTGCACTCGGTTCCAGCGGTCCGCACTCGAACGGTTATTCGCTGGTGCGCAAAATTCTTGAA" +"GTCAGCGGTTGTGATCCGCAAACCACCGAACTTGATGGTAAGCCATTAGCCGATCATCTGCTGGCACCGACCCGCATTTA" +"CGTGAAGTCAGTGCTGGAGTTGATTGAAAAGGTCGATGTGCATGCCATTGCGCACCTGACCGGCGGCGGCTTCTGGGAAA" +"ACATTCCGCGCGTATTGCCAGATAATACCCAGGCAGTGATTGATGAATCTTCCTGGCAGTGGCCGGAAGTGTTCAACTGG" +"CTGCAAACGGCAGGTAACGTTGAGCGCCATGAAATGTATCGCACCTTCAACTGCGGCGTCGGGATGATTATCGCCCTGCC" +"TGCTCCGGAAGTGGACAAAGCCCTCGCCCTGCTCAATGCCAACGGTGAAAACGCGTGGAAAATCGGTATCATCAAAGCCT" +"CTGATTCCGAACAACGCGTGGTTATCGAATAATGAATATTGTGGTGCTTATTTCCGGCAACGGAAGTAATTTACAGGCAA" +"TTATTGACGCCTGTAAAACCAACAAAATTAAAGGCACCGTACGGGCAGTTTTCAGCAATAAGGCCGACGCGTTCGGCCTT" +"GAACGCGCCCGCCAGGCGGGTATTGCAACGCATACGCTCATCGCCAGCGCGTTTGACAGTCGTGAAGCCTATGACCGGGA" +"GTTGATTCATGAAATCGACATGTACGCACCCGATGTGGTCGTGCTGGCTGGTTTTATGCGCATTCTCAGCCCGGCGTTTG" +"TCTCCCACTATGCCGGGCGTTTGCTGAACATTCACCCTTCTCTGCTGCCGAAATATCCCGGATTACACACCCATCGTCAA" +"GCGCTGGAAAATGGCGATGAAGAGCACGGTACATCGGTGCATTTCGTCACCGATGAACTGGACGGTGGCCCGGTTATTTT" +"ACAGGCGAAAGTCCCGGTATTTGCTGGTGATACGGAAGATGACGTCACCGCCCGCGTGCAAACCCAGGAACACGCCATTT" +"ATCCACTGGTGATTAGCTGGTTTGCCGATGGTCGTCTGAAAATGCACGAAAACGCCGCGTGGCTGGATGGTCAACGTCTG" +"CCGCCGCAGGGCTACGCTGCCGACGAGTAATGCCCCCGTAGTTAAAGCGCCAGCTCTGCCGCTGGCGTTTTTCAATTCAC" +"CTGTAAATCGCAAGCTCCAGCAGTTTTTTTCCCCCTTTTCTGGCATAGTTGGACATCTGCCAATATTGCTCGCCATAATA" +"TCCAGGCAGTGTCCCGTGAATAAAACGGAGTAAAAGTGGTAATGGGTCAGGAAAAGCTATACATCGAAAAAGAGCTCAGT" +"TGGTTATCGTTCAATGAACGCGTGCTTCAGGAAGCGGCGGACAAATCTAACCCGCTGATTGAAAGGATGCGTTTCCTGGG" +"GATCTATTCCAATAACCTTGATGAGTTCTATAAAGTCCGCTTCGCTGAACTGAAGCGACGCATCATTATTAGCGAAGAAC" +"AAGGCTCCAACTCTCATTCCCGCCATTTACTGGGCAAAATTCAGTCCCGGGTGCTGAAAGCCGATCAGGAATTCGACGGC" +"CTCTACAACGAGCTATTGCTGGAGATGGCGCGCAACCAGATCTTCCTGATTAATGAACGCCAGCTCTCCGTCAATCAACA" +"AAACTGGCTGCGTCATTATTTTAAGCAGTATCTGCGTCAGCACATTACGCCGATTTTAATCAATCCTGACACTGACTTAG" +"TGCAGTTCCTGAAAGATGATTACACCTATCTGGCGGTGGAAATTATCCGTGGCGATACCATCCGTTACGCGCTTCTGGAG" +"ATCCCATCAGATAAAGTGCCGCGCTTTGTGAATTTACCGCCAGAAGCGCCGCGTCGACGCAAGCCGATGATTCTTCTGGA" +"TAACATTCTGCGTTACTGCCTTGATGATATTTTCAAAGGCTTCTTTGATTATGACGCGCTGAATGCCTATTCAATGAAGA" +"TGACCCGCGATGCCGAATACGATTTAGTGCATGAGATGGAAGCCAGCCTGATGGAGTTGATGTCTTCCAGTCTCAAGCAG" +"CGTTTAACTGCTGAGCCGGTGCGTTTTGTTTATCAGCGCGATATGCCCAATGCGCTGGTTGAAGTTTTACGCGAAAAACT"; + +byte* +modify(byte *seq, int *len, double p) +{ + byte *head, *new; + + head = calloc(SEQLEN+1, sizeof(byte)); + new = head; + for (; *seq != '\0'; seq++) { + if (rng·bernoulli(p)) { + switch (rng·randi(5)) { + case 0: *new++ = 'A'; break; + case 1: *new++ = 'C'; break; + case 2: *new++ = 'G'; break; + case 3: *new++ = 'T'; break; + case 4: continue; + } + } else { + *new++ = *seq; + } + } + *new = '\0'; + *len = new - head; + return head; +} + +#define NSEQS 20 +int +main() +{ + int n, i, l, lens[NSEQS]; + byte *seqs[NSEQS]; + + int locs[aln·N][NSEQS][aln·L]; + int *loc[aln·N]; + uint64 vals[aln·N][NSEQS][aln·L]; + uint64 *val[aln·N]; + rng·init(0); + + seqs[0] = SEQ; + lens[0] = SEQLEN; + + for (n = 0; n < aln·N; n++) { + for (i = 0; i < NSEQS; i++) { + for (l = 0; l < aln·L; l++) { + vals[n][i][l] = 0; + } + } + } + + for (i = 1; i < NSEQS; i++) { + seqs[i] = modify(SEQ, lens + i, .01*i); + } + + for (i = 0; i < NSEQS; i++) { + for (n = 0; n < aln·N; n++) { + val[n] = vals[n][i]; + loc[n] = locs[n][i]; + } + aln·sketch(seqs[i], aln·L, val, loc); + } + + // for (n = 0; n < aln·N; n++) { + // printf("iteration %d\n", n); + // printf("[\n"); + // for (i = 0; i < NSEQS; i++) { + // printf(" ["); + // for (l = 0; l < aln·L; l++) { + // printf("%lu,", vals[n][i][l]); + // } + // printf("],\n"); + // } + // printf("]\n"); + // } + + for (n = 0; n < aln·N; n++) { + aln·sort(NSEQS, aln·L, (uint64*)vals[n]); + } +} diff --git a/sys/libbio/test.c b/sys/libbio/test.c index 8054864..da29c84 100644 --- a/sys/libbio/test.c +++ b/sys/libbio/test.c @@ -99,6 +99,8 @@ test·newick() io·Peeker rdr; io·Putter wtr; + bio·Node **end, **it, **list; + heap = mem·makearena(mem·sys, nil); rdr = (io·Peeker){.get = &io·getbyte, .unget = &io·ungetbyte}; wtr = (io·Putter){.put = &io·putbyte, .putstr = &io·putstring}; @@ -109,8 +111,22 @@ test·newick() t.h = heap; t.heap = (mem·Allocator){ .alloc = &mem·arenaalloc, .free = nil, }; - err = bio·readnewick(rdr, fd[0], &t); + if (err = bio·readnewick(rdr, fd[0], &t), err) { + errorf("failed to read newick"); + return 1; + } + printf("number of children: %d\n", t.root->nchild); + phylo·ladderize(t.root); + + list = mem·arenaalloc(heap, t.nleaf, sizeof(**list)); + phylo·getleafs(t, list); + for (it = list, end = list + t.nleaf; it != end; ++it) { + printf("Leaf '%s'\n", (*it)->name); + } + + bio·Node *path[100]; + // phylo·diameter(t, path); printf("Loaded tree with %d leafs and %d nodes\n", t.nleaf, t.root->nnode); err = bio·writenewick(t, wtr, fd[1]); diff --git a/sys/libmath/blas.c b/sys/libmath/blas.c index c2f7e6c..bbbe1c8 100644 --- a/sys/libmath/blas.c +++ b/sys/libmath/blas.c @@ -3,7 +3,6 @@ #include #include - #include // ----------------------------------------------------------------------- @@ -382,7 +381,6 @@ void blas·gemm(int n1, int n2, int n3, double a, double *m1, double *m2, double b, double *m3) { int i, j, k, len; - double prod; // TODO: Is there anyway this computation can be integrated into the one below? for (i = 0; i < n1; i++) { diff --git a/sys/libn/rules.mk b/sys/libn/rules.mk index 4731558..ffb41e4 100644 --- a/sys/libn/rules.mk +++ b/sys/libn/rules.mk @@ -24,7 +24,7 @@ LIBS_$(d) := $(d)/libn.a LIBS_$(d) := $(patsubst $(SRC_DIR)/%, $(OBJ_DIR)/%, $(LIBS_$(d))) LIBS := $(LIBS) $(LIBS_$(d)) -BINS_$(d) := $(d)/test +BINS_$(d) := $(d)/test $(d)/scratch BINS_$(d) := $(patsubst $(SRC_DIR)/%, $(OBJ_DIR)/%, $(BINS_$(d))) BINS := $(BINS) $(BINS_$(d)) @@ -37,7 +37,7 @@ $(LIBS_$(d)): $(OBJS_$(d)) $(ARCHIVE) $(BINS_$(d)): TCLIBS := $(LIBS_$(d)) $(LIB_DIR)/vendor/libz.a -$(BINS_$(d)): $(OBJS_$(d)) $(OBJ_DIR)/libn/test.o +$(BINS_$(d)): $(OBJ_DIR)/libn/scratch.o $(LINK) # ---- Pop off stack ---- diff --git a/sys/libn/test.c b/sys/libn/test.c index 3a21a71..973245f 100644 --- a/sys/libn/test.c +++ b/sys/libn/test.c @@ -1,5 +1,6 @@ #include #include +#include #include @@ -122,7 +123,7 @@ error test·sort() { clock_t t; - int i, test[1000000]; + int i, test[10000]; for (i = 0; i < arrlen(test); i++) { test[i] = rand(); } @@ -154,6 +155,56 @@ test·sort() return 0; } +#define HASH(str) hash_string(str) +#define EQUAL(str, str2) (str == str2) + +typedef struct Map +{ + MAP_STRUCT_BODY(byte*, float) +} Map; + +Map* +makemap(mem·Allocator heap, void* h) +{ + MAP_MAKE(Map); +} + +void +mapfree(Map *map) +{ + MAP_FREE(map); +} + +void +mapreset(Map *map) +{ + MAP_RESET(map); +} + +double +mapget(Map *map, byte *key) +{ + MAP_GET(map, key, HASH, EQUAL); +} + +static +error +mapresize(Map *map, int n) +{ + MAP_GROW(map, byte*, float, n, HASH); +} + +static +int +mapput(Map *map, byte *key, float val, error *err) +{ + MAP_PUT(map, key, val, HASH, EQUAL, mapresize, err); +} + + +#undef HASH +#undef EQUAL + error main() { -- cgit v1.2.1