openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
SF_sort.h
Go to the documentation of this file.
1 // SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
2 // SPDX-License-Identifier: Apache-2.0
3 
12 #ifndef _SF_SORT_H
13 #define _SF_SORT_H
14 
15 #include <type_traits>
16 
17 #include "SF_vector.h"
18 #include "hashmap.hpp"
19 
20 
21 namespace SF {
22 
26 template<class T> inline
27 void _binary_sort(T *_P, T *_Q, T s)
28 {
29  T *P = _P, *Q = _Q;
30  T p=0, q=0;
31  while(P != Q)
32  {
33  if(((p = P[0]) & s) != 0) break;
34  P++;
35  }
36  while(P != Q)
37  {
38  if(((q = Q[-1]) & s) == 0) break;
39  Q--;
40  }
41  while(P != Q)
42  {
43  *P++ = q;
44  *--Q = p;
45  while(((p = P[0]) & s) == 0) P++;
46  while(((q = Q[-1]) & s) != 0) Q--;
47  }
48  s >>= 1;
49  if(s)
50  {
51  if(_Q - P > 1) _binary_sort(P, _Q, s);
52  if(Q - _P > 1) _binary_sort(_P, Q, s);
53  }
54 }
55 
59 template<class T, class S> inline
60 void _binary_sort_copy(T *_P, T *_Q, S *_U, S* _V, T s)
61 {
62  T *P = _P, *Q = _Q;
63  S *U = _U, *V = _V;
64  T p=0, q=0;
65  while(P != Q)
66  {
67  if(((p = P[0]) & s) != 0) break;
68  P++;
69  U++;
70  }
71  while(P != Q)
72  {
73  if(((q = Q[-1]) & s) == 0) break;
74  Q--;
75  V--;
76  }
77  while(P != Q)
78  {
79  *P++ = q;
80  *--Q = p;
81  S u, v;
82  u = *U;
83  v = *--V;
84  *U++ = v;
85  *V = u;
86  while(((p = P[0]) & s) == 0) P++, U++;
87  while(((q = Q[-1]) & s) != 0) Q--, V--;
88  }
89  s >>= 1;
90  if(s)
91  {
92  if(_Q - P > 1) _binary_sort_copy(P, _Q, U, _V, s);
93  if(Q - _P > 1) _binary_sort_copy(_P, Q, _U, V, s);
94  }
95 }
96 
100 template<class T, class S, class R> inline
101 void _binary_sort_copy_copy(T *_P, T *_Q, S *_A, S *_B, R *_U, R *_V, T s)
102 {
103  T *P = _P, *Q = _Q;
104  S *A = _A, *B = _B;
105  R *U = _U, *V = _V;
106  T p=0, q=0;
107  while(P != Q)
108  {
109  if(((p = P[0]) & s) != 0) break;
110  P++;
111  A++;
112  U++;
113  }
114  while(P != Q)
115  {
116  if(((q = Q[-1]) & s) == 0) break;
117  Q--;
118  B--;
119  V--;
120  }
121  while(P != Q)
122  {
123  *P++ = q;
124  *--Q = p;
125  S u, v;
126  u = *U;
127  v = *--V;
128  *U++ = v;
129  *V = u;
130  S a, b;
131  a = *A;
132  b = *--B;
133  *A++ = b;
134  *B = a;
135  while(((p = P[0]) & s) == 0) P++, A++, U++;
136  while(((q = Q[-1]) & s) != 0) Q--, B--, V--;
137  }
138  s >>= 1;
139  if(s)
140  {
141  if(_Q - P > 1) _binary_sort_copy_copy(P, _Q, A, _B, U, _V, s);
142  if(Q - _P > 1) _binary_sort_copy_copy(_P, Q, _A, B, _U, V, s);
143  }
144 }
145 
149 template<class T> inline
150 void _binary_sort_sort(T *_P, T *_Q, T *_A, T *_B, T s, T t)
151 {
152  T *P = _P, *Q = _Q, *A = _A, *B = _B;
153  T p=0, q=0;
154  while(P != Q)
155  {
156  if(((p = P[0]) & s) != 0) break;
157  P++;
158  A++;
159  }
160  while(P != Q)
161  {
162  if(((q = Q[-1]) & s) == 0) break;
163  Q--;
164  B--;
165  }
166  while(P != Q)
167  {
168  *P++ = q;
169  *--Q = p;
170  T a, b;
171  a = *A;
172  b = *--B;
173  *A++ = b;
174  *B = a;
175  while(((p = P[0]) & s) == 0) P++, A++;
176  while(((q = Q[-1]) & s) != 0) Q--, B--;
177  }
178  s >>= 1;
179  if(s)
180  {
181  if(_Q - P > 1) _binary_sort_sort(P, _Q, A, _B, s, t);
182  if(Q - _P > 1) _binary_sort_sort(_P, Q, _A, B, s, t);
183  }
184  else if(t)
185  {
186  if(_Q - P > 1) _binary_sort_sort(A, _B, P, _Q, t, s);
187  if(Q - _P > 1) _binary_sort_sort(_A, B, _P, Q, t, s);
188  }
189 }
190 
194 template<class T, class S> inline
195 void _binary_sort_sort_copy(T *_P, T *_Q, T *_A, T *_B, S *_U, S* _V, T s, T t)
196 {
197  T *P = _P, *Q = _Q, *A = _A, *B = _B;
198  S *U = _U, *V = _V;
199  T p=0, q=0;
200  while(P != Q)
201  {
202  if(((p = P[0]) & s) != 0) break;
203  P++;
204  A++;
205  U++;
206  }
207  while(P != Q)
208  {
209  if(((q = Q[-1]) & s) == 0) break;
210  Q--;
211  B--;
212  V--;
213  }
214  while(P != Q)
215  {
216  *P++ = q;
217  *--Q = p;
218  T a, b;
219  a = *A;
220  b = *--B;
221  *A++ = b;
222  *B = a;
223  S u, v;
224  u = *U;
225  v = *--V;
226  *U++ = v;
227  *V = u;
228  while(((p = P[0]) & s) == 0) P++, A++, U++;
229  while(((q = Q[-1]) & s) != 0) Q--, B--, V--;
230  }
231  s >>= 1;
232  if(s)
233  {
234  if(_Q - P > 1) _binary_sort_sort_copy(P, _Q, A, _B, U, _V, s, t);
235  if(Q - _P > 1) _binary_sort_sort_copy(_P, Q, _A, B, _U, V, s, t);
236  }
237  else if(t)
238  {
239  if(_Q - P > 1) _binary_sort_sort_copy(A, _B, P, _Q, U, _V, t, s);
240  if(Q - _P > 1) _binary_sort_sort_copy(_A, B, _P, Q, _U, V, t, s);
241  }
242 }
243 
244 
248 template<class T> inline
249 T _binary_log(const T *P, const T *Q)
250 {
251  // the mask is walked up one bit at a time, which is only defined on an unsigned type
252  using mask_t = typename std::make_unsigned<T>::type;
253 
254  mask_t s = 0;
255  while(P != Q)
256  {
257  s |= mask_t(*P++);
258  }
259  mask_t t = ~mask_t(0);
260  while(s & t)
261  {
262  s &= t;
263  t <<= 1;
264  }
265  return T(s);
266 }
267 
268 
273 template<class T> inline
275 {
276  if(_V.size() < 2) return;
277  _binary_sort(&_V[0], &_V[0]+_V.size(), _binary_log(&_V[0], &_V[0]+_V.size()));
278 }
279 
285 template<class T, class S> inline
287 {
288  if(_V.size() < 2) return;
289  _binary_sort_copy(&_V[0], &_V[0]+_V.size(), &_W[0], &_W[0]+_W.size(), _binary_log(&_V[0], &_V[0]+_V.size()));
290 }
291 
298 template<class T, class S, class R> inline
300 {
301  if(_V.size() < 2) return;
302  _binary_sort_copy_copy(&_V[0], &_V[0]+_V.size(), &_W[0], &_W[0]+_W.size(), &_A[0], &_A[0]+_A.size(), _binary_log(&_V[0], &_V[0]+_V.size()));
303 }
304 
310 template<class T> inline
312 {
313  if(_V.size() < 2) return;
314  _binary_sort_sort(&_V[0], &_V[0]+_V.size(), &_W[0], &_W[0]+_W.size(),
315  _binary_log(&_V[0], &_V[0]+_V.size()), _binary_log(&_W[0], &_W[0]+_W.size()));
316 }
317 
324 template<class T, class S> inline
326 {
327  if(_V.size() < 2) return;
328  _binary_sort_sort_copy(&_V[0], &_V[0]+_V.size(), &_W[0], &_W[0]+_W.size(), &_A[0], &_A[0]+_A.size(),
329  _binary_log(&_V[0], &_V[0]+_V.size()), _binary_log(&_W[0], &_W[0]+_W.size()));
330 }
331 
332 
337 template<class T> inline
339 {
340  if(_P.size() < 2) return;
341  T *P = &_P[0], *Q = &_P[0] + _P.size();
342 
343  if(P != Q)
344  {
345  T* R = P;
346  ++P;
347  while(P != Q)
348  {
349  if ((*R != *P))
350  {
351  *++R = *P;
352  }
353  ++P;
354  }
355  ++R;
356  _P.resize(int(R - &_P[0]));
357  }
358 }
359 
365 template<class T> inline
367 {
368  if(_P.size() < 2) return;
369  T *P = &_P[0], *Q = &_P[0] + _P.size();
370  T *U = &_U[0];
371 
372  if(P != Q)
373  {
374  T* R = P;
375  T* W = U;
376  ++P; ++U;
377  while(P != Q)
378  {
379  if ((*R != *P) || (*W != *U))
380  {
381  *++R = *P;
382  *++W = *U;
383  }
384  ++P; ++U;
385  }
386  ++R;
387  ++W;
388  _P.resize(int(R - &_P[0]));
389  _U.resize(int(W - &_U[0]));
390  }
391 }
392 
398 template<class T, class S> inline
400 {
401  if(_P.size() < 2) return;
402  T *P = &_P[0], *Q = &_P[0] + _P.size();
403  S *A = &_A[0];
404 
405  if(P != Q)
406  {
407  T* R = P;
408  S* C = A;
409  ++P; ++A;
410  while(P != Q)
411  {
412  if (*R == *P)
413  {
414  *C += *A;
415  }
416  else
417  {
418  *++R = *P;
419  *++C = *A;
420  }
421  ++P; ++A;
422  }
423  ++R;
424  ++C;
425  _P.resize(int(R - &_P[0]));
426  _A.resize(int(C - &_A[0]));
427  }
428 }
429 
436 template<class T, class S> inline
438 {
439  if(_P.size() < 2) return;
440  T *P = &_P[0], *Q = &_P[0] + _P.size();
441  T *U = &_U[0];
442  S *A = &_A[0];
443 
444  if(P != Q)
445  {
446  T* R = P;
447  T* W = U;
448  S* C = A;
449  ++P; ++U; ++A;
450  while(P != Q)
451  {
452  if ((*R == *P) && (*W == *U))
453  {
454  *C += *A;
455  }
456  else
457  {
458  *++R = *P;
459  *++W = *U;
460  *++C = *A;
461  }
462  ++P; ++U; ++A;
463  }
464  ++R;
465  ++W;
466  ++C;
467  _P.resize(int(R - &_P[0]));
468  _U.resize(int(W - &_U[0]));
469  _A.resize(int(C - &_A[0]));
470  }
471 }
472 
483 template<class T> inline
484 void global_to_local_sorted(const vector<T> & glob, vector<T> & data, bool sortedData, bool doWarn)
485 {
486  size_t dsize = data.size(), gsize = glob.size();
487  vector<T> perm;
488 
489  // sort data if necessary
490  if(!sortedData) {
491  perm.resize(dsize);
492  interval(perm, 0, dsize);
493  binary_sort_copy(data, perm);
494  }
495 
496  size_t cidx, ridx, widx=0;
497  for (ridx=0, cidx=0; ridx<dsize; ridx++)
498  {
499  while (glob[cidx] < data[ridx] && cidx < (gsize-1)) cidx++;
500  if(data[ridx] == glob[cidx]) {
501  data[widx] = cidx;
502  widx++;
503  }
504  }
505 
506  // check if all values could be mapped. if its not the case, then permuting back makes no sense
507  if(widx < ridx) {
508  data.resize(widx);
509  if(doWarn)
510  fprintf(stderr, "global_to_local warning: Not all indices could be mapped!\n");
511  }
512  else {
513  if(!sortedData)
514  binary_sort_copy(perm, data);
515  }
516 }
517 
518 
529 template<class T> inline
530 void global_to_local(const vector<T> & glob, vector<T> & data, bool sortedData, bool doWarn)
531 {
532  size_t dsize = data.size(), gsize = glob.size();
533  vector<T> _glob(glob), perm_glob(glob.size()), perm_data;
534 
535  interval(perm_glob, 0, gsize);
536  binary_sort_copy(_glob, perm_glob);
537 
538  // sort data if necessary
539  if(!sortedData) {
540  perm_data.resize(dsize);
541  interval(perm_data, 0, dsize);
542  binary_sort_copy(data, perm_data);
543  }
544 
545  size_t cidx, ridx, widx=0;
546  for (ridx=0, cidx=0; ridx<dsize; ridx++)
547  {
548  while (_glob[cidx] < data[ridx] && cidx < (gsize-1)) cidx++;
549  if(data[ridx] == _glob[cidx]) {
550  data[widx] = perm_glob[cidx];
551  widx++;
552  }
553  }
554 
555  // check if all values could be mapped. if its not the case, then permuting back makes no sense
556  if(widx < ridx) {
557  data.resize(widx);
558  if(doWarn)
559  fprintf(stderr, "global_to_local warning: Not all indices could be mapped!\n");
560  }
561  else {
562  if(!sortedData)
563  binary_sort_copy(perm_data, data);
564  }
565 }
566 
576 template<class T> inline
577 void global_to_local(const hashmap::unordered_map<T,T> & g2l, vector<T> & data, bool doWarn)
578 {
579  size_t widx = 0;
580  for(size_t i=0; i<data.size(); i++) {
581  auto it = g2l.find(data[i]);
582  if(it != g2l.end()) {
583  data[widx++] = it->second;
584  }
585  }
586 
587  // check if all values could be mapped. if its not the case, then permuting back makes no sense
588  if(widx < data.size()) {
589  data.resize(widx);
590  if(doWarn)
591  fprintf(stderr, "%s warning: Not all indices could be mapped!\n", __func__);
592  }
593 }
594 
595 }
596 #endif
597 
The vector class and related algorithms.
A vector storing arbitrary data.
Definition: SF_vector.h:28
size_t size() const
The current size of the vector.
Definition: SF_vector.h:89
void resize(size_t n)
Resize a vector.
Definition: SF_vector.h:194
iterator find(const K &key)
Search for key. Return iterator.
Definition: hashmap.hpp:626
Classes similar to unordered_set and unordered_map, but with better performance.
Definition: dense_mat.hpp:19
void _binary_sort_copy_copy(T *_P, T *_Q, S *_A, S *_B, R *_U, R *_V, T s)
Definition: SF_sort.h:101
void interval(vector< T > &vec, size_t start, size_t end)
Create an integer interval between start and end.
Definition: SF_vector.h:335
void binary_sort_copy_copy(vector< T > &_V, vector< S > &_W, vector< R > &_A)
Definition: SF_sort.h:299
void unique_accumulate(vector< T > &_P, vector< S > &_A)
Definition: SF_sort.h:399
void binary_sort_copy(vector< T > &_V, vector< S > &_W)
Definition: SF_sort.h:286
void unique_resize(vector< T > &_P)
Definition: SF_sort.h:338
void _binary_sort(T *_P, T *_Q, T s)
Definition: SF_sort.h:27
T _binary_log(const T *P, const T *Q)
Definition: SF_sort.h:249
void _binary_sort_sort(T *_P, T *_Q, T *_A, T *_B, T s, T t)
Definition: SF_sort.h:150
void _binary_sort_sort_copy(T *_P, T *_Q, T *_A, T *_B, S *_U, S *_V, T s, T t)
Definition: SF_sort.h:195
void _binary_sort_copy(T *_P, T *_Q, S *_U, S *_V, T s)
Definition: SF_sort.h:60
void global_to_local(const vector< T > &glob, vector< T > &data, bool sortedData, bool doWarn)
Definition: SF_sort.h:530
void binary_sort(vector< T > &_V)
Definition: SF_sort.h:274
void binary_sort_sort_copy(vector< T > &_V, vector< T > &_W, vector< S > &_A)
Definition: SF_sort.h:325
void binary_sort_sort(vector< T > &_V, vector< T > &_W)
Definition: SF_sort.h:311
void global_to_local_sorted(const vector< T > &glob, vector< T > &data, bool sortedData, bool doWarn)
Definition: SF_sort.h:484