     int comp_size = min(bsize, n_rows - offset);
     if (tid < comp_size) {
       for (unsigned int j=0; j < n_cols; ++j) {
-        smem[tid*n_cols+j] = sorted_coords[(tid+offset)*n_cols+j];
-      }
+        smem_coords[tid*n_cols+j] = sorted_coords[(tid+offset)*n_cols+j];
+      }
+      smem_ref_states[tid] = ref_states[(tid+offset)];
     }
     __syncthreads();
     if (gid < i_to) {
-      unsigned int self_id = gid + 1;
-      unsigned int cluster_id = clustering[gid];
-      if (cluster_id == 0) {
-        cluster_id = self_id;
-      }
+      unsigned int min_id = clustering[gid];
+      // load reference coordinates for re-use into shared memory
       unsigned int ref_id = tid+bsize;
-      // load reference coordinates for re-use into shared memory
       for (unsigned int j=0; j < n_cols; ++j) {
-        smem[ref_id*n_cols+j] = sorted_coords[gid*n_cols+j];
-      }
+        smem_coords[ref_id*n_cols+j] = sorted_coords[gid*n_cols+j];
+      }
+      // check against reference structures: if close enough, they'll be
+      // considered as the same state (choosing min. state id).
+      // state lumping will be done on host.
       for (unsigned int i=0; i < comp_size; ++i) {
-        float dist2 = 0.0f;
-        for (unsigned int j=0; j < n_cols; ++j) {
-          float c = smem[ref_id*n_cols+j] - smem[i*n_cols+j];
-          dist2 = fma(c, c, dist2);
-        }
-        if (dist2 < max_dist2) {
-          unsigned int min_id = min(cluster_id, i+offset+1);
-          if (min_id == self_id) {
-            cluster_id = max(cluster_id, i+offset+1);
-          } else {
-            cluster_id = min_id;
+        if (min_id != smem_ref_states[i]
+         && smem_ref_states[i] != 0) {
+          float dist2 = 0.0f;
+          for (unsigned int j=0; j < n_cols; ++j) {
+            float c = smem_coords[ref_id*n_cols+j] - smem_coords[i*n_cols+j];
+            dist2 = fma(c, c, dist2);
           }
-        }
-      }
-      clustering[gid] = cluster_id;
-    }
-  }
-
-//  unsigned int
-//  follow(const std::vector<unsigned int>& clustering
-//       , unsigned int id
-//       , unsigned int min_id = 0) {
-//    if (id == 0
-//     || id == min_id) {
-//      return id;
-//    } else {
-//      return follow(clustering
-//                  , clustering[id-1]
-//                  , (min_id == 0)
-//                      ? id
-//                      : std::min(id, min_id));
-//    }
-//  }
+          if (dist2 < max_dist2
+           && smem_ref_states[i] < min_id) {
+            min_id = smem_ref_states[i];
+          }
+        }
+      }
+      clustering[gid] = min_id;
+    }
+  }
 
   std::vector<unsigned int>
-  sanitize_state_names(std::vector<unsigned int> clustering) {
-
-    //TODO: lump frames into microstates.
-    //      be careful: cluster ids do not necessarily
-    //                  denote frame ids (especially not with previous
-    //                  clustering)
-
-    // map names to ascending index
-    unsigned int id=0;
-    std::map<unsigned int, unsigned int> name_dict;
-    for (unsigned int name: unique_names) {
-      ++id;
-      name_dict[name] = id;
-    }
-    // relabel states
-    for (unsigned int& i: clustering) {
-      i = name_dict[i];
+  lumped_clusters(std::vector<unsigned int> clustering
+                , std::vector<unsigned int> ref_states
+                , unsigned int first_frame_above_threshold) {
+    std::unordered_set<unsigned int> sstates;
+    sstates.insert(clustering.begin()
+                 , clustering.begin() + first_frame_above_threshold);
+    sstates.insert(ref_states.begin()
+                 , ref_states.begin() + first_frame_above_threshold);
+    sstates.erase(0);
+    std::vector<unsigned int> states(sstates.begin(), sstates.end());
+    unsigned int n_states = states.size();
+    // identify lumps
+    unsigned int i, j, s, s1, s2;
+    std::map<unsigned int, unsigned int> state_mapping;
+    for (i=0; i < n_states; ++i) {
+      state_mapping[states[i]] = states[i];
+    }
+    #pragma omp parallel for\
+      default(none)\
+      private(i,j,s1,s2,s)\
+      firstprivate(n_states,first_frame_above_threshold)\
+      shared(state_mapping,states,clustering,ref_states)
+    for (i=0; i < n_states; ++i) {
+      s = states[i];
+      for (j=0; j < first_frame_above_threshold; ++j) {
+        s1 = clustering[j];
+        s2 = ref_states[j];
+        //TODO: no loop over states, just over frames,
+        //      then select 's' from {s1,s2}, compute min and update
+        //      state_mapping in critical section
+        if (s1 == s || s2 == s) {
+          state_mapping[s] = std::min(state_mapping[s]
+                                    , std::min(s1, s2));
+        }
+      }
+    }
+    // assign new (ordered) names to lumps
+    std::set<unsigned int> distinct_states;
+    for (auto sm: state_mapping) {
+      distinct_states.insert(sm.second);
+    }
+    std::map<unsigned int, unsigned int> remapped_distinct;
+    unsigned int i_state = 0;
+    for (unsigned int s: distinct_states) {
+      remapped_distinct[s] = ++i_state;
+    }
+    for (auto& sm: state_mapping) {
+      sm.second = remapped_distinct[sm.second];
+    }
+    // remap states to lumps
+    #pragma omp parallel for\
+      default(none)\
+      private(i)\
+      firstprivate(first_frame_above_threshold)\
+      shared(clustering,state_mapping)
+    for (i=0; i < first_frame_above_threshold; ++i) {
+      clustering[i] = state_mapping[clustering[i]];
     }
     return clustering;
   }
 
-//TODO remove old code
-//    for (unsigned int i=0; i < clustering.size(); ++i) {
-//      clustering[i] = follow(clustering
-//                           , clustering[i]);
-//    }
-//    return clustering;
-//  }
-
+  std::vector<unsigned int>
+  lumped_clusters__old(std::vector<unsigned int> clustering
+                , std::vector<unsigned int> ref_states
+                , unsigned int first_frame_above_threshold) {
+    unsigned int i_state = 0;
+    std::map<unsigned int, std::unordered_set<unsigned int>> lumps;
+    std::set<std::pair<unsigned int, unsigned int>> pairings;
+    for (unsigned int i=0; i < first_frame_above_threshold; ++i) {
+      unsigned int s1 = clustering[i];
+      unsigned int s2 = ref_states[i];
+      pairings.emplace(s1, s2);
+    }
+    std::unordered_set<unsigned int> seen_states;
+    for (auto p: pairings) {
+      unsigned int s1 = p.first;
+      unsigned int s2 = p.second;
+      if (s1 == s2 && seen_states.count(s1) > 0) {
+        continue;
+      } else {
+        seen_states.insert(s1);
+        seen_states.insert(s2);
+      }
+      unsigned int l1=0, l2=0;
+      // find lumps of states
+      for (auto l: lumps) {
+        if (l.second.count(s1) > 0) {
+          l1 = l.first;
+        }
+        if (l.second.count(s2) > 0) {
+          l2 = l.first;
+        }
+        if (l1 != 0 && l2 != 0) {
+          break;
+        }
+      }
+      if (l1 == 0 && l2 == 0) {
+        // no lump found: make new one
+        lumps[++i_state] = {s1, s2};
+      } else if (l1 == 0 && l2 != 0) {
+        // s2 has lump: add s1
+        lumps[l2].insert(s1);
+      } else if (l1 != 0 && l2 == 0) {
+        // s1 has lump: add s2
+        lumps[l1].insert(s2);
+      } else if (l1 != l2) {
+        // s1 and s2 have different lumps: join them
+        auto mm = std::minmax(l1,l2);
+        lumps[mm.first].insert(lumps[mm.second].begin()
+                             , lumps[mm.second].end());
+        lumps.erase(mm.second);
+      }
+    }
+    // recluster trajectory
+    i_state = 0;
+    std::unordered_map<unsigned int, unsigned int> dict;
+    for (auto l: lumps) {
+      ++i_state;
+      for (auto s: l.second) {
+        dict[s] = i_state;
+      }
+    }
+    for (unsigned int i=0; i < clustering.size(); ++i) {
+      if (clustering[i] > 0) {
+        clustering[i] = dict[clustering[i]];
+      }
+    }
+    return clustering;
+  }
 
   std::vector<std::size_t>
   initial_density_clustering(const std::vector<float>& free_energy

                                , unsigned int n_cols
                                , float max_dist2
                                , unsigned int* clustering
                                , unsigned int i_from
                                , unsigned int i_to) {
    extern __shared__ float smem[];
    unsigned int bid = blockIdx.x;
    unsigned int tid = threadIdx.x;
    unsigned int bsize = blockDim.x;

#include "logger.hpp"

#include <algorithm>
#include <unordered_set>
#include <unordered_map>

#include <cuda.h>
#include <omp.h>
