Browse code

embed samtools sources directly, since Rsamtools has moved to htslib

Michael Lawrence authored on 11/02/2019 20:33:07
Showing1 changed files
1 1
new file mode 100644
... ...
@@ -0,0 +1,295 @@
1
+/* The MIT License
2
+
3
+   Copyright (c) 2008 Genome Research Ltd (GRL).
4
+
5
+   Permission is hereby granted, free of charge, to any person obtaining
6
+   a copy of this software and associated documentation files (the
7
+   "Software"), to deal in the Software without restriction, including
8
+   without limitation the rights to use, copy, modify, merge, publish,
9
+   distribute, sublicense, and/or sell copies of the Software, and to
10
+   permit persons to whom the Software is furnished to do so, subject to
11
+   the following conditions:
12
+
13
+   The above copyright notice and this permission notice shall be
14
+   included in all copies or substantial portions of the Software.
15
+
16
+   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17
+   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
18
+   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
19
+   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
20
+   BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
21
+   ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
22
+   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
23
+   SOFTWARE.
24
+*/
25
+
26
+/* Contact: Heng Li <lh3@sanger.ac.uk> */
27
+
28
+/*
29
+  2012-12-11 (0.1.4):
30
+
31
+    * Defined __ks_insertsort_##name as static to compile with C99.
32
+
33
+  2008-11-16 (0.1.4):
34
+
35
+    * Fixed a bug in introsort() that happens in rare cases.
36
+
37
+  2008-11-05 (0.1.3):
38
+
39
+    * Fixed a bug in introsort() for complex comparisons.
40
+
41
+	* Fixed a bug in mergesort(). The previous version is not stable.
42
+
43
+  2008-09-15 (0.1.2):
44
+
45
+	* Accelerated introsort. On my Mac (not on another Linux machine),
46
+	  my implementation is as fast as std::sort on random input.
47
+
48
+	* Added combsort and in introsort, switch to combsort if the
49
+	  recursion is too deep.
50
+
51
+  2008-09-13 (0.1.1):
52
+
53
+	* Added k-small algorithm
54
+
55
+  2008-09-05 (0.1.0):
56
+
57
+	* Initial version
58
+
59
+*/
60
+
61
+#ifndef AC_KSORT_H
62
+#define AC_KSORT_H
63
+
64
+#include <stdlib.h>
65
+#include <string.h>
66
+
67
+
68
+/* Rsamtools */
69
+#if defined(_WIN32)
70
+#define drand48() ((double)rand() / RAND_MAX)
71
+#else
72
+#ifndef _SVID_SOURCE
73
+#define _SVID_SOURCE            /* Rsamtools: c99  */
74
+#endif
75
+#endif
76
+
77
+typedef struct {
78
+	void *left, *right;
79
+	int depth;
80
+} ks_isort_stack_t;
81
+
82
+#define KSORT_SWAP(type_t, a, b) { register type_t t=(a); (a)=(b); (b)=t; }
83
+
84
+#define KSORT_INIT(name, type_t, __sort_lt)								\
85
+	void ks_mergesort_##name(size_t n, type_t array[], type_t temp[])	\
86
+	{																	\
87
+		type_t *a2[2], *a, *b;											\
88
+		int curr, shift;												\
89
+																		\
90
+		a2[0] = array;													\
91
+		a2[1] = temp? temp : (type_t*)malloc(sizeof(type_t) * n);		\
92
+		for (curr = 0, shift = 0; (1ul<<shift) < n; ++shift) {			\
93
+			a = a2[curr]; b = a2[1-curr];								\
94
+			if (shift == 0) {											\
95
+				type_t *p = b, *i, *eb = a + n;							\
96
+				for (i = a; i < eb; i += 2) {							\
97
+					if (i == eb - 1) *p++ = *i;							\
98
+					else {												\
99
+						if (__sort_lt(*(i+1), *i)) {					\
100
+							*p++ = *(i+1); *p++ = *i;					\
101
+						} else {										\
102
+							*p++ = *i; *p++ = *(i+1);					\
103
+						}												\
104
+					}													\
105
+				}														\
106
+			} else {													\
107
+				size_t i, step = 1ul<<shift;							\
108
+				for (i = 0; i < n; i += step<<1) {						\
109
+					type_t *p, *j, *k, *ea, *eb;						\
110
+					if (n < i + step) {									\
111
+						ea = a + n; eb = a;								\
112
+					} else {											\
113
+						ea = a + i + step;								\
114
+						eb = a + (n < i + (step<<1)? n : i + (step<<1)); \
115
+					}													\
116
+					j = a + i; k = a + i + step; p = b + i;				\
117
+					while (j < ea && k < eb) {							\
118
+						if (__sort_lt(*k, *j)) *p++ = *k++;				\
119
+						else *p++ = *j++;								\
120
+					}													\
121
+					while (j < ea) *p++ = *j++;							\
122
+					while (k < eb) *p++ = *k++;							\
123
+				}														\
124
+			}															\
125
+			curr = 1 - curr;											\
126
+		}																\
127
+		if (curr == 1) {												\
128
+			type_t *p = a2[0], *i = a2[1], *eb = array + n;				\
129
+			for (; p < eb; ++i) *p++ = *i;								\
130
+		}																\
131
+		if (temp == 0) free(a2[1]);										\
132
+	}																	\
133
+	void ks_heapadjust_##name(size_t i, size_t n, type_t l[])			\
134
+	{																	\
135
+		size_t k = i;													\
136
+		type_t tmp = l[i];												\
137
+		while ((k = (k << 1) + 1) < n) {								\
138
+			if (k != n - 1 && __sort_lt(l[k], l[k+1])) ++k;				\
139
+			if (__sort_lt(l[k], tmp)) break;							\
140
+			l[i] = l[k]; i = k;											\
141
+		}																\
142
+		l[i] = tmp;														\
143
+	}																	\
144
+	void ks_heapmake_##name(size_t lsize, type_t l[])					\
145
+	{																	\
146
+		size_t i;														\
147
+		for (i = (lsize >> 1) - 1; i != (size_t)(-1); --i)				\
148
+			ks_heapadjust_##name(i, lsize, l);							\
149
+	}																	\
150
+	void ks_heapsort_##name(size_t lsize, type_t l[])					\
151
+	{																	\
152
+		size_t i;														\
153
+		for (i = lsize - 1; i > 0; --i) {								\
154
+			type_t tmp;													\
155
+			tmp = *l; *l = l[i]; l[i] = tmp; ks_heapadjust_##name(0, i, l); \
156
+		}																\
157
+	}																	\
158
+	static inline void __ks_insertsort_##name(type_t *s, type_t *t)		\
159
+	{																	\
160
+		type_t *i, *j, swap_tmp;										\
161
+		for (i = s + 1; i < t; ++i)										\
162
+			for (j = i; j > s && __sort_lt(*j, *(j-1)); --j) {			\
163
+				swap_tmp = *j; *j = *(j-1); *(j-1) = swap_tmp;			\
164
+			}															\
165
+	}																	\
166
+	void ks_combsort_##name(size_t n, type_t a[])						\
167
+	{																	\
168
+		const double shrink_factor = 1.2473309501039786540366528676643; \
169
+		int do_swap;													\
170
+		size_t gap = n;													\
171
+		type_t tmp, *i, *j;												\
172
+		do {															\
173
+			if (gap > 2) {												\
174
+				gap = (size_t)(gap / shrink_factor);					\
175
+				if (gap == 9 || gap == 10) gap = 11;					\
176
+			}															\
177
+			do_swap = 0;												\
178
+			for (i = a; i < a + n - gap; ++i) {							\
179
+				j = i + gap;											\
180
+				if (__sort_lt(*j, *i)) {								\
181
+					tmp = *i; *i = *j; *j = tmp;						\
182
+					do_swap = 1;										\
183
+				}														\
184
+			}															\
185
+		} while (do_swap || gap > 2);									\
186
+		if (gap != 1) __ks_insertsort_##name(a, a + n);					\
187
+	}																	\
188
+	void ks_introsort_##name(size_t n, type_t a[])						\
189
+	{																	\
190
+		int d;															\
191
+		ks_isort_stack_t *top, *stack;									\
192
+		type_t rp, swap_tmp;											\
193
+		type_t *s, *t, *i, *j, *k;										\
194
+																		\
195
+		if (n < 1) return;												\
196
+		else if (n == 2) {												\
197
+			if (__sort_lt(a[1], a[0])) { swap_tmp = a[0]; a[0] = a[1]; a[1] = swap_tmp; } \
198
+			return;														\
199
+		}																\
200
+		for (d = 2; 1ul<<d < n; ++d);									\
201
+		stack = (ks_isort_stack_t*)malloc(sizeof(ks_isort_stack_t) * ((sizeof(size_t)*d)+2)); \
202
+		top = stack; s = a; t = a + (n-1); d <<= 1;						\
203
+		while (1) {														\
204
+			if (s < t) {												\
205
+				if (--d == 0) {											\
206
+					ks_combsort_##name(t - s + 1, s);					\
207
+					t = s;												\
208
+					continue;											\
209
+				}														\
210
+				i = s; j = t; k = i + ((j-i)>>1) + 1;					\
211
+				if (__sort_lt(*k, *i)) {								\
212
+					if (__sort_lt(*k, *j)) k = j;						\
213
+				} else k = __sort_lt(*j, *i)? i : j;					\
214
+				rp = *k;												\
215
+				if (k != t) { swap_tmp = *k; *k = *t; *t = swap_tmp; }	\
216
+				for (;;) {												\
217
+					do ++i; while (__sort_lt(*i, rp));					\
218
+					do --j; while (i <= j && __sort_lt(rp, *j));		\
219
+					if (j <= i) break;									\
220
+					swap_tmp = *i; *i = *j; *j = swap_tmp;				\
221
+				}														\
222
+				swap_tmp = *i; *i = *t; *t = swap_tmp;					\
223
+				if (i-s > t-i) {										\
224
+					if (i-s > 16) { top->left = s; top->right = i-1; top->depth = d; ++top; } \
225
+					s = t-i > 16? i+1 : t;								\
226
+				} else {												\
227
+					if (t-i > 16) { top->left = i+1; top->right = t; top->depth = d; ++top; } \
228
+					t = i-s > 16? i-1 : s;								\
229
+				}														\
230
+			} else {													\
231
+				if (top == stack) {										\
232
+					free(stack);										\
233
+					__ks_insertsort_##name(a, a+n);						\
234
+					return;												\
235
+				} else { --top; s = (type_t*)top->left; t = (type_t*)top->right; d = top->depth; } \
236
+			}															\
237
+		}																\
238
+	}																	\
239
+	/* This function is adapted from: http://ndevilla.free.fr/median/ */ \
240
+	/* 0 <= kk < n */													\
241
+	type_t ks_ksmall_##name(size_t n, type_t arr[], size_t kk)			\
242
+	{																	\
243
+		type_t *low, *high, *k, *ll, *hh, *mid;							\
244
+		low = arr; high = arr + n - 1; k = arr + kk;					\
245
+		for (;;) {														\
246
+			if (high <= low) return *k;									\
247
+			if (high == low + 1) {										\
248
+				if (__sort_lt(*high, *low)) KSORT_SWAP(type_t, *low, *high); \
249
+				return *k;												\
250
+			}															\
251
+			mid = low + (high - low) / 2;								\
252
+			if (__sort_lt(*high, *mid)) KSORT_SWAP(type_t, *mid, *high); \
253
+			if (__sort_lt(*high, *low)) KSORT_SWAP(type_t, *low, *high); \
254
+			if (__sort_lt(*low, *mid)) KSORT_SWAP(type_t, *mid, *low);	\
255
+			KSORT_SWAP(type_t, *mid, *(low+1));							\
256
+			ll = low + 1; hh = high;									\
257
+			for (;;) {													\
258
+				do ++ll; while (__sort_lt(*ll, *low));					\
259
+				do --hh; while (__sort_lt(*low, *hh));					\
260
+				if (hh < ll) break;										\
261
+				KSORT_SWAP(type_t, *ll, *hh);							\
262
+			}															\
263
+			KSORT_SWAP(type_t, *low, *hh);								\
264
+			if (hh <= k) low = ll;										\
265
+			if (hh >= k) high = hh - 1;									\
266
+		}																\
267
+	}																	\
268
+	void ks_shuffle_##name(size_t n, type_t a[])						\
269
+	{																	\
270
+		int i, j;														\
271
+		for (i = n; i > 1; --i) {										\
272
+			type_t tmp;													\
273
+			j = (int)(drand48() * i);									\
274
+			tmp = a[j]; a[j] = a[i-1]; a[i-1] = tmp;					\
275
+		}																\
276
+	}
277
+
278
+#define ks_mergesort(name, n, a, t) ks_mergesort_##name(n, a, t)
279
+#define ks_introsort(name, n, a) ks_introsort_##name(n, a)
280
+#define ks_combsort(name, n, a) ks_combsort_##name(n, a)
281
+#define ks_heapsort(name, n, a) ks_heapsort_##name(n, a)
282
+#define ks_heapmake(name, n, a) ks_heapmake_##name(n, a)
283
+#define ks_heapadjust(name, i, n, a) ks_heapadjust_##name(i, n, a)
284
+#define ks_ksmall(name, n, a, k) ks_ksmall_##name(n, a, k)
285
+#define ks_shuffle(name, n, a) ks_shuffle_##name(n, a)
286
+
287
+#define ks_lt_generic(a, b) ((a) < (b))
288
+#define ks_lt_str(a, b) (strcmp((a), (b)) < 0)
289
+
290
+typedef const char *ksstr_t;
291
+
292
+#define KSORT_INIT_GENERIC(type_t) KSORT_INIT(type_t, type_t, ks_lt_generic)
293
+#define KSORT_INIT_STR KSORT_INIT(str, ksstr_t, ks_lt_str)
294
+
295
+#endif