SF_sort.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------

#ifndef _SF_SORT_H
#define _SF_SORT_H

#include "SF_vector.h"
#include "hashmap.hpp"


namespace SF {

template<class T> inline
void _binary_sort(T *_P, T *_Q, T s)
{
  T *P = _P, *Q = _Q;
  T p=0, q=0;
  while(P != Q)
  {
    if(((p = P[0]) & s) != 0) break;
    P++;
  }
  while(P != Q)
  {
    if(((q = Q[-1]) & s) == 0) break;
    Q--;
  }
  while(P != Q)
  {
    *P++ = q;
    *--Q = p;
    while(((p = P[0]) & s) == 0) P++;
    while(((q = Q[-1]) & s) != 0) Q--;
  }
  s >>= 1;
  if(s)
  {
    if(_Q - P > 1) _binary_sort(P, _Q, s);
    if(Q - _P > 1) _binary_sort(_P, Q, s);
  }
}

template<class T, class S> inline
void _binary_sort_copy(T *_P, T *_Q, S *_U, S* _V, T s)
{
  T *P = _P, *Q = _Q;
  S *U = _U, *V = _V;
  T p=0, q=0;
  while(P != Q)
  {
    if(((p = P[0]) & s) != 0) break;
    P++;
    U++;
  }
  while(P != Q)
  {
    if(((q = Q[-1]) & s) == 0) break;
    Q--;
    V--;
  }
  while(P != Q)
  {
    *P++ = q;
    *--Q = p;
    S u, v;
    u = *U;
    v = *--V;
    *U++ = v;
    *V = u;
    while(((p = P[0]) & s) == 0) P++, U++;
    while(((q = Q[-1]) & s) != 0) Q--, V--;
  }
  s >>= 1;
  if(s)
  {
    if(_Q - P > 1) _binary_sort_copy(P, _Q, U, _V, s);
    if(Q - _P > 1) _binary_sort_copy(_P, Q, _U, V, s);
  }
}

template<class T, class S, class R> inline
void _binary_sort_copy_copy(T *_P, T *_Q, S *_A, S *_B, R *_U, R *_V, T s)
{
  T *P = _P, *Q = _Q;
  S *A = _A, *B = _B;
  R *U = _U, *V = _V;
  T p=0, q=0;
  while(P != Q)
  {
    if(((p = P[0]) & s) != 0) break;
    P++;
    A++;
    U++;
  }
  while(P != Q)
  {
    if(((q = Q[-1]) & s) == 0) break;
    Q--;
    B--;
    V--;
  }
  while(P != Q)
  {
    *P++ = q;
    *--Q = p;
    S u, v;
    u = *U;
    v = *--V;
    *U++ = v;
    *V = u;
    S a, b;
    a = *A;
    b = *--B;
    *A++ = b;
    *B = a;
    while(((p = P[0]) & s) == 0) P++, A++, U++;
    while(((q = Q[-1]) & s) != 0) Q--, B--, V--;
  }
  s >>= 1;
  if(s)
  {
    if(_Q - P > 1) _binary_sort_copy_copy(P, _Q, A, _B, U, _V, s);
    if(Q - _P > 1) _binary_sort_copy_copy(_P, Q, _A, B, _U, V, s);
  }
}

template<class T> inline
void _binary_sort_sort(T *_P, T *_Q, T *_A, T *_B, T s, T t)
{
  T *P = _P, *Q = _Q, *A = _A, *B = _B;
  T p=0, q=0;
  while(P != Q)
  {
    if(((p = P[0]) & s) != 0) break;
    P++;
    A++;
  }
  while(P != Q)
  {
    if(((q = Q[-1]) & s) == 0) break;
    Q--;
    B--;
  }
  while(P != Q)
  {
    *P++ = q;
    *--Q = p;
    T a, b;
    a = *A;
    b = *--B;
    *A++ = b;
    *B = a;
    while(((p = P[0]) & s) == 0) P++, A++;
    while(((q = Q[-1]) & s) != 0) Q--, B--;
  }
  s >>= 1;
  if(s)
  {
    if(_Q - P > 1) _binary_sort_sort(P, _Q, A, _B, s, t);
    if(Q - _P > 1) _binary_sort_sort(_P, Q, _A, B, s, t);
  }
  else if(t)
  {
    if(_Q - P > 1) _binary_sort_sort(A, _B, P, _Q, t, s);
    if(Q - _P > 1) _binary_sort_sort(_A, B, _P, Q, t, s);
  }
}

template<class T, class S> inline
void _binary_sort_sort_copy(T *_P, T *_Q, T *_A, T *_B, S *_U, S* _V, T s, T t)
{
  T *P = _P, *Q = _Q, *A = _A, *B = _B;
  S *U = _U, *V = _V;
  T p=0, q=0;
  while(P != Q)
  {
    if(((p = P[0]) & s) != 0) break;
    P++;
    A++;
    U++;
  }
  while(P != Q)
  {
    if(((q = Q[-1]) & s) == 0) break;
    Q--;
    B--;
    V--;
  }
  while(P != Q)
  {
    *P++ = q;
    *--Q = p;
    T a, b;
    a = *A;
    b = *--B;
    *A++ = b;
    *B = a;
    S u, v;
    u = *U;
    v = *--V;
    *U++ = v;
    *V = u;
    while(((p = P[0]) & s) == 0) P++, A++, U++;
    while(((q = Q[-1]) & s) != 0) Q--, B--, V--;
  }
  s >>= 1;
  if(s)
  {
    if(_Q - P > 1) _binary_sort_sort_copy(P, _Q, A, _B, U, _V, s, t);
    if(Q - _P > 1) _binary_sort_sort_copy(_P, Q, _A, B, _U, V, s, t);
  }
  else if(t)
  {
    if(_Q - P > 1) _binary_sort_sort_copy(A, _B, P, _Q, U, _V, t, s);
    if(Q - _P > 1) _binary_sort_sort_copy(_A, B, _P, Q, _U, V, t, s);
  }
}


template<class T> inline
T _binary_log(const T *P, const T *Q)
{
  T s = 0;
  while(P != Q)
  {
    s |= *P++;
  }
  T t = ~0;
  while(s & t)
  {
    s &= t;
    t <<= 1;
  }
  return s;
}


template<class T> inline
void binary_sort(vector<T> &_V)
{
  if(_V.size() < 2) return;
  _binary_sort(&_V[0], &_V[0]+_V.size(), _binary_log(&_V[0], &_V[0]+_V.size()));
}

template<class T, class S> inline
void binary_sort_copy(vector<T> &_V, vector<S> &_W)
{
  if(_V.size() < 2) return;
  _binary_sort_copy(&_V[0], &_V[0]+_V.size(), &_W[0], &_W[0]+_W.size(), _binary_log(&_V[0], &_V[0]+_V.size()));
}

template<class T, class S, class R> inline
void binary_sort_copy_copy(vector<T> &_V, vector<S> &_W, vector<R> &_A)
{
  if(_V.size() < 2) return;
  _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()));
}

template<class T> inline
void binary_sort_sort(vector<T> &_V, vector<T> &_W)
{
  if(_V.size() < 2) return;
  _binary_sort_sort(&_V[0], &_V[0]+_V.size(), &_W[0], &_W[0]+_W.size(), 
    _binary_log(&_V[0], &_V[0]+_V.size()), _binary_log(&_W[0], &_W[0]+_W.size()));
}

template<class T, class S> inline
void binary_sort_sort_copy(vector<T> &_V, vector<T> &_W, vector<S> &_A)
{
  if(_V.size() < 2) return;
  _binary_sort_sort_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()), _binary_log(&_W[0], &_W[0]+_W.size()));
}


template<class T> inline
void unique_resize(vector<T> &_P)
{
  if(_P.size() < 2) return;
  T *P = &_P[0], *Q = &_P[0] + _P.size();

  if(P != Q)
  {
    T* R = P;
    ++P;
    while(P != Q)
    {
      if ((*R != *P))
      {
        *++R = *P;
      }
      ++P;
    }
    ++R;
    _P.resize(int(R - &_P[0]));
  }
}

template<class T> inline
void unique_resize(vector<T> &_P, vector<T> &_U)
{
  if(_P.size() < 2) return;
  T *P = &_P[0], *Q = &_P[0] + _P.size();
  T *U = &_U[0];

  if(P != Q)
  {
    T* R = P;
    T* W = U;
    ++P; ++U;
    while(P != Q)
    {
      if ((*R != *P) || (*W != *U))
      {
        *++R = *P;
        *++W = *U;
      }
      ++P; ++U;
    }
    ++R;
    ++W;
    _P.resize(int(R - &_P[0]));
    _U.resize(int(W - &_U[0]));
  }
}

template<class T, class S> inline
void unique_accumulate(vector<T> &_P, vector<S> &_A)
{
  if(_P.size() < 2) return;
  T *P = &_P[0], *Q = &_P[0] + _P.size();
  S *A = &_A[0];

  if(P != Q)
  {
    T* R = P;
    S* C = A;
    ++P; ++A;
    while(P != Q)
    {
      if (*R == *P)
      {
        *C += *A;
      }
      else
      {
        *++R = *P;
        *++C = *A;
      }
      ++P; ++A;
    }
    ++R;
    ++C;
    _P.resize(int(R - &_P[0]));
    _A.resize(int(C - &_A[0]));
  }
}

template<class T, class S> inline
void unique_accumulate(vector<T> &_P, vector<T> &_U, vector<S> &_A)
{
  if(_P.size() < 2) return;
  T *P = &_P[0], *Q = &_P[0] + _P.size();
  T *U = &_U[0];
  S *A = &_A[0];

  if(P != Q)
  {
    T* R = P;
    T* W = U;
    S* C = A;
    ++P; ++U; ++A;
    while(P != Q)
    {
      if ((*R == *P) && (*W == *U))
      {
        *C += *A;
      }
      else
      {
        *++R = *P;
        *++W = *U;
        *++C = *A;
      }
      ++P; ++U; ++A;
    }
    ++R;
    ++W;
    ++C;
    _P.resize(int(R - &_P[0]));
    _U.resize(int(W - &_U[0]));
    _A.resize(int(C - &_A[0]));
  }
}

template<class T> inline
void global_to_local_sorted(const vector<T> & glob, vector<T> & data, bool sortedData, bool doWarn)
{
  size_t dsize = data.size(), gsize = glob.size();
  vector<T> perm;

  // sort data if necessary
  if(!sortedData) {
    perm.resize(dsize);
    interval(perm, 0, dsize);
    binary_sort_copy(data, perm);
  }

  size_t cidx, ridx, widx=0;
  for (ridx=0, cidx=0; ridx<dsize; ridx++)
  {
    while (glob[cidx] < data[ridx] && cidx < (gsize-1)) cidx++;
    if(data[ridx] == glob[cidx]) {
      data[widx] = cidx;
      widx++;
    }
  }

  // check if all values could be mapped. if its not the case, then permuting back makes no sense
  if(widx < ridx) {
    data.resize(widx);
    if(doWarn)
      fprintf(stderr, "global_to_local warning: Not all indices could be mapped!\n");
  }
  else {
    if(!sortedData)
      binary_sort_copy(perm, data);
  }
}


template<class T> inline
void global_to_local(const vector<T> & glob, vector<T> & data, bool sortedData, bool doWarn)
{
  size_t dsize = data.size(), gsize = glob.size();
  vector<T> _glob(glob), perm_glob(glob.size()), perm_data;

  interval(perm_glob, 0, gsize);
  binary_sort_copy(_glob, perm_glob);

  // sort data if necessary
  if(!sortedData) {
    perm_data.resize(dsize);
    interval(perm_data, 0, dsize);
    binary_sort_copy(data, perm_data);
  }

  size_t cidx, ridx, widx=0;
  for (ridx=0, cidx=0; ridx<dsize; ridx++)
  {
    while (_glob[cidx] < data[ridx] && cidx < (gsize-1)) cidx++;
    if(data[ridx] == _glob[cidx]) {
      data[widx] = perm_glob[cidx];
      widx++;
    }
  }

  // check if all values could be mapped. if its not the case, then permuting back makes no sense
  if(widx < ridx) {
    data.resize(widx);
    if(doWarn)
      fprintf(stderr, "global_to_local warning: Not all indices could be mapped!\n");
  }
  else {
    if(!sortedData)
      binary_sort_copy(perm_data, data);
  }
}

template<class T> inline
void global_to_local(const hashmap::unordered_map<T,T> & g2l, vector<T> & data, bool doWarn)
{
  size_t widx = 0;
  for(size_t i=0; i<data.size(); i++) {
    auto it = g2l.find(data[i]);
    if(it != g2l.end()) {
      data[widx++] = it->second;
    }
  }

  // check if all values could be mapped. if its not the case, then permuting back makes no sense
  if(widx < data.size()) {
    data.resize(widx);
    if(doWarn)
      fprintf(stderr, "%s warning: Not all indices could be mapped!\n", __func__);
  }
}

}
#endif