CMlib
Cell mapping algorithms in C++
cscm.h
Go to the documentation of this file.
1 #ifndef CELL_MAPPING_CPP_CSCM_H
2 #define CELL_MAPPING_CPP_CSCM_H
3 
4 #include "scm.h"
5 
6 namespace cm {
7 
8  template <class IDType>
9  class CellTree {
10  private:
11  IDType cmid;
12  std::vector<IDType> cells;
13  CellState state;
14  IDType imageCell;
15  IDType imageCmid; // TODO: Rename CMID to ClasterID
16  public:
18  state = CellState::Untouched;
19  cmid = 0;
20  imageCell = 0;
21  }
22  CellTree(IDType cmid,
23  const std::vector<IDType> &cells,
24  const CellState &state,
25  IDType imageCell) : cmid(cmid), cells(cells), state(state), imageCell(imageCell) {
26  /* Nothing to do yet */
27  }
28  CellTree(IDType cmid,
29  const std::vector<IDType> &cells,
30  const CellState &state) : cmid(cmid), cells(cells), state(state) {
31  /* Nothing to do yet */
32  }
33  void setClusterID(IDType cmid) {
34  CellTree::cmid = cmid;
35  }
36  void setCells(const std::vector<IDType> &cells) {
37  CellTree::cells = cells;
38  }
39  void setState(const CellState &state) {
40  CellTree::state = state;
41  }
42  IDType getClusterID() const {
43  return cmid;
44  }
45  const std::vector<IDType> &getCells() const {
46  return cells;
47  }
48  const CellState &getState() const {
49  return state;
50  }
51  IDType getImageCell() const {
52  return imageCell;
53  }
54  void setImageCell(IDType imageCell) {
55  CellTree::imageCell = imageCell;
56  }
57  IDType getImageCmid() const {
58  return imageCmid;
59  }
60  void setImageCmid(IDType imageCmid) {
61  CellTree::imageCmid = imageCmid;
62  }
63  };
64 
65  template <class CellType, class IDType, class StateVectorType>
66  class ClusteredSCM {
67  private:
70  std::vector<CellTree<IDType>> cellTrees;
71  void join_stage1(DynamicalSystemBase<StateVectorType>* systemp,
74  std::vector<IDType>& sinkDoA1) {
83  IDType cmid1 = scm1p->getCss().getCell(0).getClusterID();
84  IDType cmid2 = scm2p->getCss().getCell(0).getClusterID();
85  IDType z; IDType imz; IDType cmidz;
86  std::vector<IDType> seq;
87  CellTree<IDType> cellTree;
88  for(size_t i=0; i<sinkDoA1.size(); i++) {
89  seq.resize(0);
90  z = sinkDoA1[i];
91  // Skip cell if already processed
92  CellState stz = scm1p->getCss().getCell(z).getState();
93  if (stz == CellState::Processed || stz == CellState::UnderProcessing) {
94  // skip cell
95  } else {
96  // State is UNTOUCHED
97  seq.push_back(z);
98  // Generate sequence from this cell until it leaves its own SCM's state space
99  bool left = false;
100  while(!left) {
101  imz = scm1p->getCss().getCell(z).getImage();
102  if (imz != 0) {
103  // Image is not the sink cell, we are still in the original CSS
104  cmidz = scm1p->getCss().getCell(imz).getClusterID();
105  if (cmidz == cmid1) {
106  // Check for cell state before accumulating cellTrees
107  CellState stimz = scm1p->getCss().getCell(imz).getState();
108  if(stimz == CellState::UnderProcessing) {
109  left = true;
110  // Get the index of the celltree which contains the touched cell
111  size_t ct = scm1p->getCss().getCell(imz).getCellTreeID();
112  // Tag cells
113  for (size_t k = 0; k < seq.size(); k++) {
114  if (scm1p->getCss().getCell(seq[k]).getState() == CellState::UnderProcessing) {
115  //cout << "ERROR while appending!\n";
116  }
117  scm1p->getCss().getCell(seq[k]).setState(CellState::UnderProcessing);
118  scm1p->getCss().getCell(seq[k]).setCellTreeID(ct);
119  }
120  // Append current seq to that cellpath
121  std::vector<IDType> cpcells = cellTrees[ct].getCells();
122  for (size_t w=0; w<cpcells.size(); w++) {
123  seq.push_back(cpcells[w]);
124  }
125  // Now the original seq contains the current seq plus the cells of the touched seq
126  // Overwrite cells of that seq
127  cellTrees[ct].setCells(seq);
128  } else if (stimz == CellState::Processed) {
129  // Set properties of the current cycle
130  left = true;
131  size_t otherG = scm1p->getCss().getCell(imz).getGroup();
132  size_t otherCmid = scm1p->getCss().getCell(imz).getClusterID();
133  for (size_t k = 0; k < seq.size(); k++) {
134  scm1p->getCss().getCell(seq[k]).setClusterID(otherCmid);
135  scm1p->getCss().getCell(seq[k]).setGroup(otherG);
136  scm1p->getCss().getCell(seq[k]).setState(CellState::Processed);
137  }
138  } else {
139  // Add cell to seq
140  seq.push_back(imz);
141  z = imz;
142  }
143  } else {
144  // We are still in the original CSS, but touching a transiting seq
145  left = true;
146  // Set properties of the current seq as the target
147  size_t otherG = scm1p->getCss().getCell(imz).getGroup();
148  size_t otherCmid = scm1p->getCss().getCell(imz).getClusterID();
149  for (size_t k = 0; k < seq.size(); k++) {
150  scm1p->getCss().getCell(seq[k]).setClusterID(otherCmid);
151  scm1p->getCss().getCell(seq[k]).setGroup(otherG);
152  scm1p->getCss().getCell(seq[k]).setState(CellState::Processed);
153  }
154  // Terminating current processing cycle
155  } // end if imz != 0
156  } else {
157  // We have left the original CSS
158  left = true;
159  // Check whether this state is in the other SCM's state space
160  StateVectorType center = systemp->step(scm1p->getCss().getCenter(z));
161  imz = scm2p->getCss().getID(center);
162  if (imz != 0) {
163  // The target is a regular cell of the other CSS, get its group number
164  size_t otherG = scm2p->getCss().getCell(imz).getGroup();
165  if (otherG != 0) {
166  // Tag cells and terminate current processing cycle
167  for (size_t k = 0; k < seq.size(); k++) {
168  scm1p->getCss().getCell(seq[k]).setClusterID(cmid2);
169  scm1p->getCss().getCell(seq[k]).setGroup(otherG);
170  scm1p->getCss().getCell(seq[k]).setState(CellState::Processed);
171  }
172  } else {
173  // Tag cells and Save this tree for later analysis, terminate proc. cycle
174  // Also set cellTreeIndex of cells in this sequence
175  size_t ctIndex = cellTrees.size();
176  for (size_t k = 0; k < seq.size(); k++) {
177  if (scm1p->getCss().getCell(seq[k]).getState() == CellState::UnderProcessing) {
178  //cout << "ERROR while saving cycle!\n";
179  }
180  scm1p->getCss().getCell(seq[k]).setState(CellState::UnderProcessing);
181  scm1p->getCss().getCell(seq[k]).setCellTreeID(ctIndex);
182  }
183  // Set cell tree properties
184  cellTree.setCells(seq);
185  cellTree.setClusterID(cmid1);
187  cellTree.setImageCell(imz);
188  cellTree.setImageCmid(cmid2);
189  cellTrees.push_back(cellTree);
190  }
191  } else {
192  // This seq leads to the (real) sink cell, tag cells, terminate proc. cycle
193  for (size_t k = 0; k < seq.size(); k++) {
194  // Since the group number was 0, it is not necessary to overwrite
195  scm1p->getCss().getCell(seq[k]).setState(CellState::Processed);
196  }
197  }
198  } /* end if */
199  } /* end while */
200  } /* end if already processed or under_processing */
201  } /* end for */
202  }
203  void join_stage2(SCM<CellType, IDType, StateVectorType>* scm1p,
205  IDType cmid1 = scm1p->getCss().getCell(0).getClusterID();
206  IDType cmid2 = scm2p->getCss().getCell(0).getClusterID();
207  /*
208  * STAGE 2: Do an SCM on the cellTrees, their state should be currently UNTOUCHED
209  */
210  for (size_t i=0; i<cellTrees.size(); i++) {
211  CellState cs = cellTrees[i].getState();
212  // Select i-th cellTree, check its state
213  if (cs == CellState::Processed) {
214  // skip
215  } else if (cs == CellState::Untouched) {
216  // Start processing cycle
217  //cout << "Processing cycle at cellTree: " << i << endl;
218  std::vector<size_t> cycle;
219  cycle.push_back(i);
220  cellTrees[i].setState(CellState::UnderProcessing);
221  size_t cpImage = i;
222  CellState cpImageState;
223  bool processing = true;
224  while (processing) {
225  // Before continuing the processing cycle, check if this cellpath leads to an already processed cell
226  size_t cpImageCell = cellTrees[cpImage].getImageCell();
227  size_t cpImageCmid = cellTrees[cpImage].getImageCmid();
228  CellState imageCellState;
229  size_t imageCellGroup;
230  // TODO: Rework this part
231  if (cpImageCmid == cmid1) {
232  imageCellState = scm1->getCss().getCell(cpImageCell).getState();
233  imageCellGroup = scm1->getCss().getCell(cpImageCell).getGroup();
234  } else {
235  imageCellState = scm2->getCss().getCell(cpImageCell).getState();
236  imageCellGroup = scm2->getCss().getCell(cpImageCell).getGroup();
237  }
238  if (imageCellState == CellState::Processed) {
239  // Tag current cell paths in cycle as the last processed image
240  // Tag internal cells (in paths) as well
241  for (size_t k=0; k<cycle.size(); k++) {
242  cellTrees[cycle[k]].setState(CellState::Processed);
243  std::vector<IDType> cellsk = cellTrees[cycle[k]].getCells();
244  if(cellTrees[cycle[k]].getClusterID() == cmid1) {
245  for (size_t ck=0; ck<cellsk.size(); ck++) {
246  scm1->getCss().getCell(cellsk[ck]).setGroup(imageCellGroup);
247  scm1->getCss().getCell(cellsk[ck]).setState(CellState::Processed);
248  scm1->getCss().getCell(cellsk[ck]).setClusterID(cmid1);
249  }
250  } else {
251  for (size_t ck=0; ck<cellsk.size(); ck++) {
252  scm2->getCss().getCell(cellsk[ck]).setGroup(imageCellGroup);
253  scm2->getCss().getCell(cellsk[ck]).setState(CellState::Processed);
254  scm2->getCss().getCell(cellsk[ck]).setClusterID(cmid1);
255  }
256  }
257  }
258  processing = false;
259  //cout << "Known destination (cell) found, terminating cycle at cellTree: " << i << endl;
260  } else {
261  // get the cellpath image index
262  cpImage = cellPathImage(cellTrees, cpImage);
263  // get cellpath images state
264  cpImageState = cellTrees[cpImage].getState();
265  if (cpImageState == CellState::Untouched) {
266  // append current cell seq to the cycle, and also tag it as UNDER_PROCESSING
267  cellTrees[cpImage].setState(CellState::UnderProcessing);
268  cycle.push_back(cpImage);
269  //cout << "Cycle size: " << cycle.size() << endl;
270  } else if (cpImageState == CellState::Processed) {
271  // TODO: Find the group number of image
272  size_t imageGroup = 1;
273  // Tag current cell paths in cycle as the last processed image
274  // Tag internal cells (in paths) as well
275  for (size_t k=0; k<cycle.size(); k++) {
276  cellTrees[cycle[k]].setState(CellState::Processed);
277  std::vector<IDType> cellsk = cellTrees[cycle[k]].getCells();
278  if(cellTrees[cycle[k]].getClusterID() == cmid1) {
279  for (size_t ck=0; ck<cellsk.size(); ck++) {
280  scm1->getCss().getCell(cellsk[ck]).setGroup(imageGroup);
281  scm1->getCss().getCell(cellsk[ck]).setState(CellState::Processed);
282  }
283  } else {
284  for (size_t ck=0; ck<cellsk.size(); ck++) {
285  scm2->getCss().getCell(cellsk[ck]).setGroup(imageGroup);
286  scm2->getCss().getCell(cellsk[ck]).setState(CellState::Processed);
287  }
288  }
289  }
290  processing = false;
291  //cout << "Known destination (cellpath) found, terminating cycle at cellTree: " << i << endl;
292  } else if (cpImageState == CellState::UnderProcessing) {
293  // new periodic group found, tag all cells in the cycle
294  // TODO: Create new group number for clusterID1
295  size_t newGroup = scm1->getPeriodicGroups() + 1;
296  scm1->setPeriodicGroups(newGroup);
297  // Tag internal cells (in paths) as well, also set cmid to cmid1
298  for (size_t k=0; k<cycle.size(); k++) {
299  cellTrees[cycle[k]].setState(CellState::Processed);
300  std::vector<IDType> cellsk = cellTrees[cycle[k]].getCells();
301  if(cellTrees[cycle[k]].getClusterID() == cmid1) {
302  for (size_t ck=0; ck<cellsk.size(); ck++) {
303  scm1->getCss().getCell(cellsk[ck]).setGroup(newGroup);
304  scm1->getCss().getCell(cellsk[ck]).setState(CellState::Processed);
305  scm1->getCss().getCell(cellsk[ck]).setClusterID(cmid1);
306  }
307  } else {
308  for (size_t ck=0; ck<cellsk.size(); ck++) {
309  scm2->getCss().getCell(cellsk[ck]).setGroup(newGroup);
310  scm2->getCss().getCell(cellsk[ck]).setState(CellState::Processed);
311  scm2->getCss().getCell(cellsk[ck]).setClusterID(cmid1); // To have the same color
312  }
313  }
314  }
315  processing = false;
316  //std::cout << "New " << cycle.size() <<" PG found, terminating cycle at cellTree: " << i << std::endl;
317  }
318  }
319  }
320  } else {
321  std::cout << "ERROR: CellTree with UNDER_PROCESSING state found!\n";
322  // This should not happen
323  }
324  }
325  }
326  public:
329  ClusteredSCM::scm1 = scm1p;
330  ClusteredSCM::scm2 = scm2p;
331  }
335  size_t cellPathImage(std::vector<CellTree<IDType>>& cellPaths, size_t cellPathId) {
336  IDType cellPathImageCell = cellPaths[cellPathId].getImageCell();
337  IDType cellPathImageCmid = cellPaths[cellPathId].getImageCmid();
338  // Loop on cellPaths and try to find the image
339  IDType targetCellPath = cellPaths.size();
340  std::vector<IDType> cells;
341  for (IDType cp = 0; cp < cellPaths.size(); cp++) {
342  if (cellPaths[cp].getClusterID() == cellPathImageCmid) {
343  cells = cellPaths[cp].getCells();
344  for (IDType cpc = 0; cpc < cells.size(); cpc++) {
345  if (cells[cpc] == cellPathImageCell) {
346  targetCellPath = cp;
347  break;
348  }
349  }
350  }
351  if (targetCellPath != cellPaths.size()) { break; }
352  }
353  if(targetCellPath == cellPaths.size()) {
354  //cout << "ERROR: Cell path image not found!\n";
355  }
356  return targetCellPath;
357  }
361  void join(bool verbose = false) {
362  /* TODO: Check overlap between the two SCM state spaces, raise error
363  * TODO: Check whether the two systems are the same!
364  */
365  // Get the DoA of sink cell for both SCMs
366  std::vector<IDType> sinkDoA1;
367  std::vector<IDType> sinkDoA2;
368  if (verbose) std::cout << "\nInitialization: " << std::endl;
369  size_t scm1count = scm1->getCss().getCellSum();
370  size_t scm2count = scm2->getCss().getCellSum();
371  IDType cmid1 = scm1->getCss().getCell(0).getClusterID();
372  IDType cmid2 = scm2->getCss().getCell(0).getClusterID();
374  if (cmid2 == cmid1) {
375  cmid2++;
376  #pragma omp parallel for
377  for(size_t i=0; i<scm2count; i++) {
378  scm2->getCss().getCell(i).setClusterID(cmid2);
379  }
380  }
381  // TODO: Check if two SCMs use the same dynamical system!
383  for(size_t i=1; i<scm1count; i++) {
384  if (scm1->getCss().getCell(i).getGroup() == 0) {
385  sinkDoA1.push_back(i);
386  scm1->getCss().getCell(i).setState(CellState::Untouched);
387  }
388  }
389  for(size_t i=1; i<scm2count; i++) {
390  if (scm2->getCss().getCell(i).getGroup() == 0) {
391  sinkDoA2.push_back(i);
392  scm2->getCss().getCell(i).setState(CellState::Untouched);
393  }
394  }
395 
396  size_t sinkCount = sinkDoA1.size() + sinkDoA2.size();
397  if (verbose) {
398  std::cout << " SCM1's sink cell DoA contains: " << sinkDoA1.size() << " cells." << std::endl;
399  std::cout << " SCM2's sink cell DoA contains: " << sinkDoA2.size() << " cells." << std::endl;
400  std::cout << " Total number of cells in the joining procedure: " << sinkCount << std::endl;
401  //scm1->generateImage("scm1_st0.jpg", &coloringMethod);
402  //scm2->generateImage("scm2_st0.jpg", &coloringMethod);
403  std::cout << "\nStage 1:" << std::endl;
404  }
405 
406  cellTrees.resize(0);
407  join_stage1(systemp, scm1, scm2, sinkDoA1);
408  join_stage1(systemp, scm2, scm1, sinkDoA2);
409  IDType pathsum = 0;
410  IDType cellStates[3];
411  if(verbose) {
412  cellStates[0] = 0;
413  cellStates[1] = 0;
414  cellStates[2] = 0;
415  for (IDType i = 0; i < sinkDoA1.size(); i++) {
416  cellStates[(int) scm1->getCss().getCell(sinkDoA1[i]).getState()]++;
417  }
418  for (IDType i = 0; i < sinkDoA2.size(); i++) {
419  cellStates[(int) scm2->getCss().getCell(sinkDoA2[i]).getState()]++;
420  }
421  //std::cout << "Cell States:\n";
422  //std::cout << " Untouched:\t\t" << cellStates[(int) CellState::Untouched] << std::endl;
423  //std::cout << " UnderProcessing:\t" << cellStates[(int) CellState::UnderProcessing] << std::endl;
424  //std::cout << " Processed:\t\t" << cellStates[(int) CellState::Processed] << std::endl;
425  //std::cout << (cellStates[0] + cellStates[1] + cellStates[2]) << "/" << sinkCount << "\n";
426  std::cout << " Processed " << cellStates[(int) CellState::Processed] << " cells (";
427  std::cout << double(cellStates[(int) CellState::Processed])/sinkCount*100.0 << "%)" << std::endl;
428  }
429  // List seq statistics, also reset their state to UNTOUCHED
430  pathsum = 0;
431  for (IDType p=0; p < cellTrees.size(); p++) {
432  //std::cout << "Path " << p << ": clusterID = " << cellTrees[p].getClusterID() << ", cells: " << cellTrees[p].getCells().size() << std::endl;
433  pathsum += cellTrees[p].getCells().size();
434  cellTrees[p].setState(CellState::Untouched);
435  }
436  if (verbose) {
437  std::cout << " Cell trees constructed: " << cellTrees.size() << std::endl;
438  std::cout << " Number of cells in cell trees: " << pathsum << std::endl;
439  //scm1->generateImage("scm1_st1.jpg", &coloringMethod);
440  //scm2->generateImage("scm2_st1.jpg", &coloringMethod);
441  std::cout << "\nStage 2:" << std::endl;
442  }
443 
444  join_stage2(scm1, scm2);
445 
446  if (verbose) {
447  cellStates[0] = 0;
448  cellStates[1] = 0;
449  cellStates[2] = 0;
450  for (size_t i = 0; i < sinkDoA1.size(); i++) {
451  cellStates[(int) scm1->getCss().getCell(sinkDoA1[i]).getState()]++;
452  }
453  for (size_t i = 0; i < sinkDoA2.size(); i++) {
454  cellStates[(int) scm2->getCss().getCell(sinkDoA2[i]).getState()]++;
455  }
456  size_t processedPaths = 0;
457  for (size_t p=0; p < cellTrees.size(); p++) {
458  if(cellTrees[p].getState() == CellState::Processed) {processedPaths++;}
459  }
460  //std::cout << "Cell States:\n";
461  //std::cout << " Untouched:\t\t" << cellStates[(int) CellState::Untouched] << std::endl;
462  //std::cout << " UnderProcessing:\t" << cellStates[(int) CellState::UnderProcessing] << std::endl;
463  //std::cout << " Processed:\t\t" << cellStates[(int) CellState::Processed] << std::endl;
464  //std::cout << (cellStates[0] + cellStates[1] + cellStates[2]) << "/" << sinkCount << "\n";
465  std::cout << " Processed " << cellStates[(int) CellState::Processed] << " cells (";
466  std::cout << double(cellStates[(int) CellState::Processed])/sinkCount*100.0 << "%)" << std::endl;
467  std::cout << " Processed cell trees: " << processedPaths << std::endl;
468  }
469 
470  // Shift SCM2's group numbers
471  IDType groupshift = scm1->getPeriodicGroups();
472  IDType totalgroups = scm1->getPeriodicGroups() + scm2->getPeriodicGroups();
473  scm1->setPeriodicGroups(totalgroups);
474  scm2->setPeriodicGroups(totalgroups);
475  #pragma omp parallel for
476  for (IDType z=1; z<scm1->getCss().getCellSum(); z++) {
477  if (scm1->getCss().getCell(z).getClusterID() == cmid2 && scm1->getCss().getCell(z).getGroup() != 0) {
478  scm1->getCss().getCell(z).setGroup(scm1->getCss().getCell(z).getGroup()+groupshift);
479  }
480  }
481  #pragma omp parallel for
482  for (IDType z=1; z<scm2->getCss().getCellSum(); z++) {
483  if (scm2->getCss().getCell(z).getClusterID() == cmid2 && scm2->getCss().getCell(z).getGroup() != 0) {
484  scm2->getCss().getCell(z).setGroup(scm2->getCss().getCell(z).getGroup()+groupshift);
485  }
486  }
487  if (verbose) {
488  scm1->generateImage("scm1_joined.jpg", &coloringMethod);
489  scm2->generateImage("scm2_joined.jpg", &coloringMethod);
490  }
491  }
492 
493  };
494 
495 }
496 
497 #endif //CELL_MAPPING_CPP_CSCM_H
CellTree()
Definition: cscm.h:17
Definition: coloring.h:50
void generateImage(std::string filepath, SCMColoringMethod< CellType, IDType > *coloringMethod=nullptr, IDType x0=0, IDType y0=0, IDType xw=0, IDType yw=0)
Definition: scm.h:123
CellState
Used during cell mapping algorithms.
Definition: cell.h:11
const CellState & getState() const
Definition: cscm.h:48
size_t cellPathImage(std::vector< CellTree< IDType >> &cellPaths, size_t cellPathId)
Finds the index of the image of a CellTree within the container cellPaths.
Definition: cscm.h:335
void setImageCmid(IDType imageCmid)
Definition: cscm.h:60
DynamicalSystemBase< StateVectorType > * getSystemPointer() const
Definition: scm.h:179
virtual StateVectorType step(const StateVectorType &state) const =0
const std::vector< IDType > & getCells() const
Definition: cscm.h:45
void setImageCell(IDType imageCell)
Definition: cscm.h:54
void setPeriodicGroups(IDType periodicGroups)
Definition: scm.h:191
void setClusterID(IDType cmid)
Definition: cscm.h:33
Definition: system.h:7
Definition: scm.h:15
IDType getImageCmid() const
Definition: cscm.h:57
CellTree(IDType cmid, const std::vector< IDType > &cells, const CellState &state)
Definition: cscm.h:28
void join(bool verbose=false)
Joins the two SCM solutions.
Definition: cscm.h:361
CellTree(IDType cmid, const std::vector< IDType > &cells, const CellState &state, IDType imageCell)
Definition: cscm.h:22
IDType getImageCell() const
Definition: cscm.h:51
Definition: cscm.h:9
Definition: cell.h:6
void setState(const CellState &state)
Definition: cscm.h:39
void setCells(const std::vector< IDType > &cells)
Definition: cscm.h:36
Definition: cscm.h:66
ClusteredSCM(SCM< CellType, IDType, StateVectorType > *scm1p, SCM< CellType, IDType, StateVectorType > *scm2p)
Definition: cscm.h:327
const SCMUniformCellStateSpace< CellType, IDType, StateVectorType > & getCss() const
Definition: scm.h:182
IDType getPeriodicGroups() const
Definition: scm.h:188
IDType getClusterID() const
Definition: cscm.h:42