TreeCompiler.h 17.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
//============================================================================
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt for details.
//
//  This software is distributed WITHOUT ANY WARRANTY; without even
//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
//  PURPOSE.  See the above copyright notice for more information.
//============================================================================
// Copyright (c) 2018, The Regents of the University of California, through
// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
// from the U.S. Dept. of Energy).  All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// (1) Redistributions of source code must retain the above copyright notice, this
//     list of conditions and the following disclaimer.
//
// (2) Redistributions in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
// (3) Neither the name of the University of California, Lawrence Berkeley National
//     Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
//     used to endorse or promote products derived from this software without
//     specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
// OF THE POSSIBILITY OF SUCH DAMAGE.
//
//=============================================================================
//
//  This code is an extension of the algorithm presented in the paper:
//  Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
//  Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
//  Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
//  (LDAV), October 2016, Baltimore, Maryland.
//
//  The PPP2 algorithm and software were jointly developed by
//  Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
//  Oliver Ruebel (LBNL)
//==============================================================================


#ifndef _TREECOMPILER_H_
#define _TREECOMPILER_H_

57
#include <iomanip>
58
59
#include <iostream>
#include <vtkm/Types.h>
60
#include <vtkm/cont/ArrayCopy.h>
61
#include <vtkm/cont/DataSet.h>
62
63
64
65
66
67
68
69
70
71
#include <vtkm/worklet/contourtree_augmented/Types.h>

namespace vtkm
{
namespace worklet
{
namespace contourtree_distributed
{
// Possibly change the following when comapring to PPP prototype
constexpr int PRINT_WIDTH = 12;
72
using dataType = vtkm::Float64;
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
using indexType = vtkm::Id;

// small class for storing the contour arcs
class Edge
{ // Edge
public:
  indexType low, high;

  // constructor - defaults to -1
  Edge(vtkm::Id Low = -1, vtkm::Id High = -1)
    : low(Low)
    , high(High)
  {
  }
}; // Edge

// comparison operator <
inline bool operator<(const Edge& LHS, const Edge& RHS)
{ // operator <
#if 0
  if (LHS.low < RHS.low) return true;
  if (LHS.low > RHS.low) return false;
  if (LHS.high < RHS.high) return true;
  if (LHS.high > RHS.high) return false;
#endif
  if (std::min(LHS.low, LHS.high) < std::min(RHS.low, RHS.high))
    return true;
  else if (std::min(LHS.low, LHS.high) > std::min(RHS.low, RHS.high))
    return false;
  if (std::max(LHS.low, LHS.high) < std::max(RHS.low, RHS.high))
    return true;
  else if (std::max(LHS.low, LHS.high) > std::max(RHS.low, RHS.high))
    return false;
  return false;
} // operator <

// comparison operator ==
inline bool operator==(const Edge& LHS, const Edge& RHS)
{ // operator ==
  return (LHS.low == RHS.low && LHS.high == RHS.high) ||
    (LHS.low == RHS.high && LHS.high == RHS.low);
} // operator ==

// a helper class which stores a single supernode inserted onto a superarc
class SupernodeOnSuperarc
{ // class SupernodeOnSuperarc
public:
  // the global ID of the supernode
  indexType globalID;
  // the data value stored at the supernode
  dataType dataValue;

  // the low and high ends of the superarc it is on (may be itself)
  indexType lowEnd, highEnd;

  // constructor
  SupernodeOnSuperarc(indexType GlobalID = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT,
                      dataType DataValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT,
                      indexType LowEnd = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT,
                      indexType HighEnd = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT)
    : globalID(GlobalID)
    , dataValue(DataValue)
    , lowEnd(LowEnd)
    , highEnd(HighEnd)
  { // constructor
  } // constructor
};  // class SupernodeOnSuperarc

// overloaded comparison operator
// primary sort is by superarc (low, high),
// then secondary sort on datavalue
// tertiary on globalID to implement simulated simplicity
inline bool operator<(const SupernodeOnSuperarc& left, const SupernodeOnSuperarc& right)
{ // < operator
  // simple lexicographic sort
  if (left.lowEnd < right.lowEnd)
    return true;
  if (left.lowEnd > right.lowEnd)
    return false;
  if (left.highEnd < right.highEnd)
    return true;
  if (left.highEnd > right.highEnd)
    return false;
  if (left.dataValue < right.dataValue)
    return true;
  if (left.dataValue > right.dataValue)
    return false;
  if (left.globalID < right.globalID)
    return true;
  if (left.globalID > right.globalID)
    return false;

  // fall-through (shouldn't happen, but)
  // if they're the same, it's false
  return false;
} // < operator

// stream output
std::ostream& operator<<(std::ostream& outStream, SupernodeOnSuperarc& node);

// stream input
std::istream& operator>>(std::istream& inStream, SupernodeOnSuperarc& node);

// the class that compiles the contour tree
class TreeCompiler
{ // class TreeCompiler
public:
  // we want a vector of supernodes on superarcs
  std::vector<SupernodeOnSuperarc> supernodes;

  // and a vector of Edges (the output)
  std::vector<Edge> superarcs;
185

186
187
188
189
190
191
192
193
194
195
196
  // routine to add a known hierarchical tree to it
  // note that this DOES NOT finalise - we don't want too many sorts
  void AddHierarchicalTree(const vtkm::cont::DataSet& addedTree);

  // routine to compute the actual superarcs
  void ComputeSuperarcs();

  // routine to print a superarcs array in our format
  static void PrintSuperarcArray(const std::vector<Edge>& superarc_array);

  // routine to print the superarcs
197
  void PrintSuperarcs(bool) const;
198
199
200
201
202
203
204

  // routine to write out binary file
  void WriteBinary(FILE* outFile) const;

  // routine to read in binary file & append to contents
  void ReadBinary(FILE* inFile);
}; // class TreeCompiler
205
206

// stream output
207
inline std::ostream& operator<<(std::ostream& outStream, SupernodeOnSuperarc& node)
208
209
210
211
212
213
214
{ // stream output
  outStream << node.lowEnd << " " << node.highEnd << " " << node.dataValue << " " << node.globalID
            << std::endl;
  return outStream;
} // stream output

// stream input
215
inline std::istream& operator>>(std::istream& inStream, SupernodeOnSuperarc& node)
216
217
218
219
220
221
222
{ // stream input
  inStream >> node.lowEnd >> node.highEnd >> node.dataValue >> node.globalID;
  return inStream;
} // stream input

// routine to add a known hierarchical tree to it
// note that this DOES NOT finalise - we don't want too many sorts
223
inline void TreeCompiler::AddHierarchicalTree(const vtkm::cont::DataSet& addedTree)
224
225
{ // TreeCompiler::AddHierarchicalTree()
  // Copy relevant tree content to STL arrays
226
  vtkm::cont::UnknownArrayHandle dataValues_array = addedTree.GetField("DataValues").GetData();
227
  std::vector<vtkm::Float64> dataValues(dataValues_array.GetNumberOfValues());
228
  auto dataValues_handle = vtkm::cont::make_ArrayHandle(dataValues, vtkm::CopyFlag::Off);
229
  vtkm::cont::ArrayCopy(dataValues_array, dataValues_handle);
230
231
232
233
234
235
  dataValues_handle.SyncControlArray();

  auto regularNodeGlobalIds_array = addedTree.GetField("RegularNodeGlobalIds").GetData();
  std::vector<vtkm::Id> regularNodeGlobalIds(regularNodeGlobalIds_array.GetNumberOfValues());
  auto regularNodeGlobalIds_handle =
    vtkm::cont::make_ArrayHandle(regularNodeGlobalIds, vtkm::CopyFlag::Off);
236
  vtkm::cont::ArrayCopy(regularNodeGlobalIds_array, regularNodeGlobalIds_handle);
237
238
239
240
241
242
  regularNodeGlobalIds_handle
    .SyncControlArray(); //Forces values to get updated if copy happened on GPU

  auto superarcs_array = addedTree.GetField("Superarcs").GetData();
  std::vector<vtkm::Id> added_tree_superarcs(superarcs_array.GetNumberOfValues());
  auto superarcs_handle = vtkm::cont::make_ArrayHandle(added_tree_superarcs, vtkm::CopyFlag::Off);
243
  vtkm::cont::ArrayCopy(superarcs_array, superarcs_handle);
244
245
246
247
248
  superarcs_handle.SyncControlArray(); //Forces values to get updated if copy happened on GPU

  auto supernodes_array = addedTree.GetField("Supernodes").GetData();
  std::vector<vtkm::Id> added_tree_supernodes(supernodes_array.GetNumberOfValues());
  auto supernodes_handle = vtkm::cont::make_ArrayHandle(added_tree_supernodes, vtkm::CopyFlag::Off);
249
  vtkm::cont::ArrayCopy(supernodes_array, supernodes_handle);
250
251
252
253
254
  supernodes_handle.SyncControlArray(); //Forces values to get updated if copy happened on GPU

  auto superparents_array = addedTree.GetField("Superparents").GetData();
  std::vector<vtkm::Id> superparents(superparents_array.GetNumberOfValues());
  auto superparents_handle = vtkm::cont::make_ArrayHandle(superparents, vtkm::CopyFlag::Off);
255
  vtkm::cont::ArrayCopy(superparents_array, superparents_handle);
256
257
258
  superparents_handle.SyncControlArray(); //Forces values to get updated if copy happened on GPU

  // loop through all of the supernodes in the hierarchical tree
259
260
  for (indexType supernode = 0; supernode < static_cast<indexType>(added_tree_supernodes.size());
       supernode++)
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
  { // per supernode
    // retrieve the regular ID for the supernode
    indexType regularId = added_tree_supernodes[supernode];
    indexType globalId = regularNodeGlobalIds[regularId];
    dataType dataVal = dataValues[regularId];

    // retrieve the supernode at the far end
    indexType superTo = added_tree_superarcs[supernode];

    // now test - if it is NO_SUCH_ELEMENT, there are two possibilities
    if (vtkm::worklet::contourtree_augmented::NoSuchElement(superTo))
    { // no Superto

      // retrieve the superparent
      indexType superparent = superparents[regularId];


      // the root node will have itself as its superparent
      if (superparent == supernode)
        continue;
      else
      { // not own superparent - an attachment point
        // retrieve the superparent's from & to
        indexType regularFrom = added_tree_supernodes[superparent];
        indexType globalFrom = regularNodeGlobalIds[regularFrom];
        indexType superParentTo = added_tree_superarcs[superparent];
        indexType regularTo =
          added_tree_supernodes[vtkm::worklet::contourtree_augmented::MaskedIndex(superParentTo)];
        indexType globalTo = regularNodeGlobalIds[regularTo];

        // test the superTo to see whether we ascend or descend
        // note that we will never have NO_SUCH_ELEMENT here
        if (vtkm::worklet::contourtree_augmented::IsAscending(superParentTo))
        { // ascending
          this->supernodes.push_back(SupernodeOnSuperarc(globalId, dataVal, globalFrom, globalTo));
        } // ascending
        else
        { // descending
          this->supernodes.push_back(SupernodeOnSuperarc(globalId, dataVal, globalTo, globalFrom));
        } // descending
      }   // not own superparent - an attachment point
    }     // no Superto
    else
    { // there is a valid superarc
      // retrieve the "to" and convert to global
      indexType maskedTo = vtkm::worklet::contourtree_augmented::MaskedIndex(superTo);
      indexType regularTo = added_tree_supernodes[maskedTo];
      indexType globalTo = regularNodeGlobalIds[regularTo];
      dataType dataTo = dataValues[regularTo];

      // test the superTo to see whether we ascend or descend
      // note that we will never have NO_SUCH_ELEMENT here
      // we add both ends
      if (vtkm::worklet::contourtree_augmented::IsAscending(superTo))
      { // ascending
        this->supernodes.push_back(SupernodeOnSuperarc(globalId, dataVal, globalId, globalTo));
        this->supernodes.push_back(SupernodeOnSuperarc(globalTo, dataTo, globalId, globalTo));
      } // ascending
      else
      { // descending
        this->supernodes.push_back(SupernodeOnSuperarc(globalId, dataVal, globalTo, globalId));
        this->supernodes.push_back(SupernodeOnSuperarc(globalTo, dataTo, globalTo, globalId));
      } // descending
    }   // there is a valid superarc
  }     // per supernode

} // TreeCompiler::AddHierarchicalTree()

// routine to compute the actual superarcs
330
inline void TreeCompiler::ComputeSuperarcs()
331
332
333
334
335
336
337
338
339
340
{ // TreeCompiler::ComputeSuperarcs()
  // first we sort the vector
  std::sort(supernodes.begin(), supernodes.end());

  // we could do a unique test here, but it's easier just to suppress it inside the loop

  // now we loop through it: note the -1
  // this is because we know a priori that the last one is the last supernode on a superarc
  // and would fail the test inside the loop. By putting it in the loop test, we avoid having
  // to have an explicit if statement inside the loop
341
342
  for (indexType supernode = 0; supernode < static_cast<vtkm::Id>(supernodes.size() - 1);
       supernode++)
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
  { // loop through supernodes
    // this is actually painfully simple: if the (lowEnd, highEnd) don't match the next one,
    // then we're at the end of the group and do nothing.  Otherwise, we link to the next one
    if ((supernodes[supernode].lowEnd != supernodes[supernode + 1].lowEnd) ||
        (supernodes[supernode].highEnd != supernodes[supernode + 1].highEnd))
      continue;

    // if the supernode matches, then we have a repeat, and can suppress
    if (supernodes[supernode].globalID == supernodes[supernode + 1].globalID)
      continue;

    // otherwise, add a superarc to the list
    superarcs.push_back(Edge(supernodes[supernode].globalID, supernodes[supernode + 1].globalID));
  } // loop through supernodes

  // now sort them
  std::sort(superarcs.begin(), superarcs.end());
} // TreeCompiler::ComputeSuperarcs()

// routine to print the superarcs
363
364
365
inline void TreeCompiler::PrintSuperarcArray(const std::vector<Edge>& superarc_array)
{ // TreeCompiler::PrintSuperarcArray()
  for (indexType superarc = 0; superarc < static_cast<indexType>(superarc_array.size()); superarc++)
366
  { // per superarc
367
    if (superarc_array[superarc].low < superarc_array[superarc].high)
368
    { // order by ID not value
369
370
      std::cout << std::setw(PRINT_WIDTH) << superarc_array[superarc].low << " ";
      std::cout << std::setw(PRINT_WIDTH) << superarc_array[superarc].high << std::endl;
371
372
373
    } // order by ID not value
    else
    { // order by ID not value
374
375
      std::cout << std::setw(PRINT_WIDTH) << superarc_array[superarc].high << " ";
      std::cout << std::setw(PRINT_WIDTH) << superarc_array[superarc].low << std::endl;
376
377
378
379
    } // order by ID not value

  } // per superarc

380
381
} // TreeCompiler::PrintSuperarcArray()

382
inline void TreeCompiler::PrintSuperarcs(bool printHeader = false) const
383
{
384
385
386
387
388
389
  if (printHeader)
  {
    std::cout << "============" << std::endl;
    std::cout << "Contour Tree" << std::endl;
  }

390
391
  PrintSuperarcArray(this->superarcs);
}
392
393

// routine to write out binary file
394
inline void TreeCompiler::WriteBinary(FILE* outFile) const
395
396
397
398
399
400
401
{ // WriteBinary()
  // do a bulk write of the entire contents
  // no error checking, no type checking, no nothing
  fwrite(&(supernodes[0]), sizeof(SupernodeOnSuperarc), supernodes.size(), outFile);
} // WriteBinary()

// routine to read in binary file and append
402
inline void TreeCompiler::ReadBinary(FILE* inFile)
403
404
405
406
407
{ // ReadBinary()
  // use fseek to jump to the end
  fseek(inFile, 0, SEEK_END);

  // use fTell to retrieve the size of the file
408
  std::size_t nBytes = ftell(inFile);
409
410
411
412
  // now rewind
  rewind(inFile);

  // compute how many elements are to be read
413
  std::size_t nSupernodes = nBytes / sizeof(SupernodeOnSuperarc);
414
415

  // retrieve the current size
416
  std::size_t currentSize = supernodes.size();
417
418
419
420
421

  // resize to add the right number
  supernodes.resize(currentSize + nSupernodes);

  // now read directly into the right chunk
422
423
  std::size_t nSupernodesRead =
    fread(&(supernodes[currentSize]), sizeof(SupernodeOnSuperarc), nSupernodes, inFile);
424

425
426
427
428
429
430
431
  if (nSupernodesRead != nSupernodes)
  {
    VTKM_LOG_S(vtkm::cont::LogLevel::Error,
               "Error: Expected to read " << nSupernodes << " supernodes but could read only "
                                          << nSupernodesRead << ". Output will be incorrect!"
                                          << std::endl);
  }
432
433
434
} // ReadBinary()

// stream output - just dumps the supernodeonsuperarcs
435
inline std::ostream& operator<<(std::ostream& outStream, TreeCompiler& tree)
436
{ // stream output
437
438
  for (indexType supernode = 0; supernode < static_cast<indexType>(tree.supernodes.size());
       supernode++)
439
440
441
442
443
    outStream << tree.supernodes[supernode];
  return outStream;
} // stream output

// stream input - reads in the supernodeonsuperarcs & appends them
444
inline std::istream& operator>>(std::istream& inStream, TreeCompiler& tree)
445
446
447
448
449
450
451
452
453
454
455
{ // stream input
  while (!inStream.eof())
  {
    SupernodeOnSuperarc tempNode;
    inStream >> tempNode;
    tree.supernodes.push_back(tempNode);
  }
  // we will overshoot, so subtract one
  tree.supernodes.resize(tree.supernodes.size() - 1);
  return inStream;
} // stream input
456
457
458
459
460
461

} // namespace contourtree_distributed
} // namespace worklet
} // namespace vtkm

#endif