//============================================================================
//  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 2014 Sandia Corporation.
//  Copyright 2014 UT-Battelle, LLC.
//  Copyright 2014 Los Alamos National Security.
//
//  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
//  the U.S. Government retains certain rights in this software.
//
//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
//  Laboratory (LANL), the U.S. Government retains certain rights in
//  this software.
//============================================================================
//  Copyright (c) 2016, Los Alamos National Security, LLC
//  All rights reserved.
//
//  Copyright 2016. Los Alamos National Security, LLC. 
//  This software was produced under U.S. Government contract DE-AC52-06NA25396 
//  for Los Alamos National Laboratory (LANL), which is operated by 
//  Los Alamos National Security, LLC for the U.S. Department of Energy. 
//  The U.S. Government has rights to use, reproduce, and distribute this 
//  software.  NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC 
//  MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE 
//  USE OF THIS SOFTWARE.  If software is modified to produce derivative works, 
//  such modified software should be clearly marked, so as not to confuse it 
//  with the version available from LANL.
//
//  Additionally, 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 Los Alamos National Security, LLC, Los Alamos 
//     National Laboratory, LANL, the U.S. Government, 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 LOS ALAMOS NATIONAL SECURITY, LLC 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 LOS ALAMOS 
//  NATIONAL SECURITY, LLC 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 based on 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.

//=======================================================================================
//
// COMMENTS:
//
//  The old ChainGraph has been abstracted a little further - it still does the same job
//  of carrying most of the intermediate stages.  However, since the chain building is 
//  also needed by the mesh to set up the initial graph input, it has been moved (for now
//  to Types.h)
//
//	There will be no explicit constructor - instead, it's the mesh's job to initialise
//	a valid object of this type
//
//=======================================================================================

#ifndef vtkm_worklet_contourtree_chaingraph_h
#define vtkm_worklet_contourtree_chaingraph_h

#include <vtkm/worklet/contourtree/PrintVectors.h>
#include <vtkm/worklet/contourtree/EdgePeakComparator.h>
#include <vtkm/worklet/contourtree/GoverningSaddleFinder.h>
#include <vtkm/worklet/contourtree/RegularPointTransferrer.h>
#include <vtkm/worklet/contourtree/VertexDegreeUpdater.h>
#include <vtkm/worklet/contourtree/ActiveEdgeTransferrer.h>
#include <vtkm/worklet/contourtree/ChainDoubler.h>
#include <vtkm/worklet/contourtree/SaddleAscentFunctor.h>
#include <vtkm/worklet/contourtree/SaddleAscentTransferrer.h>
#include <vtkm/worklet/contourtree/TrunkBuilder.h>
#include <vtkm/worklet/contourtree/JoinTreeTransferrer.h>

#include <vtkm/cont/ArrayHandlePermutation.h>
#include <vtkm/worklet/DispatcherMapField.h>

#define DEBUG_PRINT 1
//#define DEBUG_FUNCTION_ENTRY 1
//#define DEBUG_TIMING 1

namespace vtkm {
namespace worklet {
namespace contourtree {

#define DEBUG_STRING_TRANSFER_GOVERNING_SADDLES "Extrema should now be assigned"
#define DEBUG_STRING_TRANSFER_SADDLE_STARTS "Transfer Saddle Starts "
#define DEBUG_STRING_TRANSFERRED_SADDLE_STARTS "Saddle Starts Transferred"
#define DEBUG_STRING_TRANSFER_TO_MERGE_TREE "Transfer to Merge Tree"
#define DEBUG_STRING_OUTDEGREE "Outdegree"
#define DEBUG_STRING_CHAINEXT "Chain Ext"
#define DEBUG_STRING_ACTIVE_OUTDEGREE "Active Outdegree"
#define DEBUG_STRING_ACTIVE_CHAINEXT "Active Chain Ext"
#define DEBUG_STRING_FAR_ID "Far"
#define DEBUG_STRING_FAR_INDEX "Far Index"
#define DEBUG_STRING_FAR_VALUE "Far Value"
#define DEBUG_STRING_NEAR_ID "Near"
#define DEBUG_STRING_NEAR_INDEX "Near Index"
#define DEBUG_STRING_NEAR_VALUE "Near Value"
#define DEBUG_STRING_EDGE_FAR_ID "Edge Far"
#define DEBUG_STRING_EDGE_NEAR_ID "Edge Near"
#define DEBUG_STRING_EDGE_NEAR_INDEX "Edge Near Index"
#define DEBUG_STRING_EDGE_NEAR_VALUE "Edge Near Value"
#define DEBUG_STRING_SORTED_NEAR_ID "Sorted Near"
#define DEBUG_STRING_SORTED_NEAR_INDEX "Sorted Near Index"
#define DEBUG_STRING_SORTED_NEAR_VALUE "Sorted Near Value"
#define DEBUG_STRING_SORTED_FAR_ID "Sorted Far"

template <typename T, typename StorageType, typename DeviceAdapter>
class ChainGraph
{
public:
        typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> DeviceAlgorithm;

	// we will want a reference to the original data array
	const vtkm::cont::ArrayHandle<T,StorageType> &values;

	// we will also want a reference to the arc array where we write the output
	vtkm::cont::ArrayHandle<vtkm::Id> &arcArray;
	
	// for each vertex, we need to know where it is in the original data array
	vtkm::cont::ArrayHandle<vtkm::Id> valueIndex;

	// and we also need the orientation of the edges (i.e. is it join or split)
	bool isJoinGraph;

	// and we will store the number of iterations the computation took here
	vtkm::Id nIterations;

	// array recording pruning sequence
	// pseudo-extrema prune to pseudo-saddles
	// all others prune to pseudo-extrema
	vtkm::cont::ArrayHandle<vtkm::Id> prunesTo;
	
	// we also want to keep track of the first edge for each vertex
	vtkm::cont::ArrayHandle<vtkm::Id> firstEdge;
	
	// and the outdegree for each vertex
	vtkm::cont::ArrayHandle<vtkm::Id> outdegree;
	
	// finally, we need to keep track of the chain extremum for each vertex
	vtkm::cont::ArrayHandle<vtkm::Id> chainExtremum;

	// we will also need to keep track of both near and far ends of each edge
	vtkm::cont::ArrayHandle<vtkm::Id> edgeFar;
	vtkm::cont::ArrayHandle<vtkm::Id> edgeNear;

	// we will also keep track of the currently active set of vertices and edges
	vtkm::cont::ArrayHandle<vtkm::Id> activeVertices;
	vtkm::cont::ArrayHandle<vtkm::Id> activeEdges;

	// and an array for sorting edges
	vtkm::cont::ArrayHandle<vtkm::Id> edgeSorter;
	
	// constructor takes necessary references
	ChainGraph(const vtkm::cont::ArrayHandle<T,StorageType> &Values,
                   vtkm::cont::ArrayHandle<vtkm::Id> &ArcArray,
                   bool IsJoinGraph) :
                         values(Values),
                         arcArray(ArcArray),
                         isJoinGraph(IsJoinGraph) {}

	// sets initial size of vertex arrays
	void AllocateVertexArrays(vtkm::Id Size);
   
	// sets initial size of edge arrays
	void AllocateEdgeArrays(vtkm::Id Size);

	// routine that builds the merge graph once the initial vertices & edges are set
	void Compute(vtkm::cont::ArrayHandle<vtkm::Id> &saddles);

	// sorts saddle starts to find governing saddles
	void FindGoverningSaddles();

	// marks now regular points for removal
	void TransferRegularPoints();	
	
	// compacts the active vertex list
	void CompactActiveVertices();

	// compacts the active edge list
	void CompactActiveEdges();
	
	// builds the chains for the new active vertices
	void BuildChains();

	// suppresses non-saddles for the governing saddles pass
	void TransferSaddleStarts();

	// sets all remaining active vertices
	void BuildTrunk();

	// transfers partial results to merge tree array
	void TransferToMergeTree(vtkm::cont::ArrayHandle<vtkm::Id> &saddles);

	// prints the contents of the topology graph in a standard format
	void DebugPrint(const char *message);
} ; // class ChainGraph

// sets initial size of vertex arrays
template<typename T, typename StorageType, typename DeviceAdapter>
void ChainGraph<T,StorageType,DeviceAdapter>::AllocateVertexArrays(vtkm::Id Size)
{
  valueIndex.Allocate(Size);
  prunesTo.Allocate(Size);
  firstEdge.Allocate(Size);
  outdegree.Allocate(Size);
  chainExtremum.Allocate(Size);
  activeVertices.Allocate(Size);
} // AllocateVertexArrays()

// sets initial size of edge arrays
template<typename T, typename StorageType, typename DeviceAdapter>
void ChainGraph<T,StorageType,DeviceAdapter>::AllocateEdgeArrays(vtkm::Id Size)
{
  edgeFar.Allocate(Size);
  edgeNear.Allocate(Size);
  activeEdges.Allocate(Size);
} // AllocateEdgeArrays()

// routine that builds the merge graph once the initial vertices & edges are set
template<typename T, typename StorageType, typename DeviceAdapter>
void ChainGraph<T,StorageType,DeviceAdapter>::Compute(vtkm::cont::ArrayHandle<vtkm::Id> &saddles)
{
#ifdef DEBUG_FUNCTION_ENTRY
		std::cout << std::endl;
		std::cout << "===================" << std::endl;
		std::cout << "Compute Chain Graph" << std::endl;
		std::cout << "===================" << std::endl;
		std::cout << std::endl;
#endif
  DebugPrint("Chain Graph Computation Starting");

  // loop until we run out of active edges
  nIterations = 0;
  while (edgeSorter.GetNumberOfValues() > 0)
  {
    // find & label the extrema with their governing saddles
    FindGoverningSaddles();
		
    // label the regular points 
    TransferRegularPoints();
		
    // compact the active set of vertices & edges
    CompactActiveVertices();
    CompactActiveEdges();
		
    // rebuild the chains
    BuildChains();
		
    // choose the subset of edges for the governing saddles
    TransferSaddleStarts();

    // increment the iteration count
    nIterations++;
  } // main loop
	
  // final pass to label the trunk vertices
  BuildTrunk();

  // we can now release many of the arrays to free up space
  firstEdge.ReleaseResources();
  outdegree.ReleaseResources();
  edgeNear.ReleaseResources();
  edgeFar.ReleaseResources();
  activeEdges.ReleaseResources();
  activeVertices.ReleaseResources();
  edgeSorter.ReleaseResources();

  // and transfer results to mergearcs
  TransferToMergeTree(saddles);
		
  // then release the remaining memory
  chainExtremum.ReleaseResources();
  prunesTo.ReleaseResources();
		
#ifdef DEBUG_PRINT
  DebugPrint("Chain Graph Computed");
#endif
} // Compute()

// sorts saddle ascents to find governing saddles
template<typename T, typename StorageType, typename DeviceAdapter>
void ChainGraph<T,StorageType,DeviceAdapter>::FindGoverningSaddles()
{
#ifdef DEBUG_FUNCTION_ENTRY
		std::cout << std::endl;
		std::cout << "======================" << std::endl;
		std::cout << "Find Governing Saddles" << std::endl;
		std::cout << "======================" << std::endl;
		std::cout << std::endl;
#endif	

  // sort with the comparator
  DeviceAlgorithm::Sort(edgeSorter,
                        EdgePeakComparator<T,StorageType,DeviceAdapter>(
                                           values.PrepareForInput(DeviceAdapter()),
                                           valueIndex.PrepareForInput(DeviceAdapter()),
                                           edgeFar.PrepareForInput(DeviceAdapter()),
                                           edgeNear.PrepareForInput(DeviceAdapter()),
                                           arcArray.PrepareForInput(DeviceAdapter()),
                                           isJoinGraph));

#ifdef DEBUG_PRINT
  DebugPrint("After Sorting");
#endif

  // now loop through the edges
  GoverningSaddleFinder governingSaddleFinder;
  vtkm::worklet::DispatcherMapField<GoverningSaddleFinder>
                 governingSaddleFinderDispatcher(governingSaddleFinder);
  vtkm::Id nEdges = edgeSorter.GetNumberOfValues();
  vtkm::cont::ArrayHandleIndex edgeIndexArray(nEdges);

  governingSaddleFinderDispatcher.Invoke(edgeIndexArray,  // input
                                         edgeSorter,      // input (whole array)
                                         edgeFar,         // input (whole array)
                                         edgeNear,        // input (whole array)
                                         prunesTo,        // output (whole array)
                                         outdegree);      // output (whole array)
#ifdef DEBUG_PRINT
  DebugPrint(DEBUG_STRING_TRANSFER_GOVERNING_SADDLES);
#endif
} // FindGoverningSaddles()

// marks now regular points for removal
template<typename T, typename StorageType, typename DeviceAdapter>
void ChainGraph<T,StorageType,DeviceAdapter>::TransferRegularPoints()
{
#ifdef DEBUG_FUNCTION_ENTRY
		std::cout << std::endl;
		std::cout << "=======================" << std::endl;
		std::cout << "Transfer Regular Points" << std::endl;
		std::cout << "=======================" << std::endl;
		std::cout << std::endl;
#endif
  RegularPointTransferrer<T> regularPointTransferrer(isJoinGraph);
  vtkm::worklet::DispatcherMapField<RegularPointTransferrer<T> >
                 regularPointTransferrerDispatcher(regularPointTransferrer);

  regularPointTransferrerDispatcher.Invoke(activeVertices,          // input
                                           chainExtremum,           // input (whole array)
                                           values,                  // input (whole array)
                                           valueIndex,              // input (whole array)
                                           prunesTo,                // i/o (whole array)
                                           outdegree);              // output (whole array)
#ifdef DEBUG_PRINT
  DebugPrint("Regular Points Should Now Be Labelled");
#endif
} // TransferRegularPoints()


// compacts the active vertex list
template<typename T, typename StorageType, typename DeviceAdapter>
void ChainGraph<T,StorageType,DeviceAdapter>::CompactActiveVertices()
{
#ifdef DEBUG_FUNCTION_ENTRY
		std::cout << std::endl;
		std::cout << "=======================" << std::endl;
		std::cout << "Compact Active Vertices" << std::endl;
		std::cout << "=======================" << std::endl;
		std::cout << std::endl;
#endif
  typedef vtkm::cont::ArrayHandle<vtkm::Id> IdArrayType;
  typedef vtkm::cont::ArrayHandlePermutation<IdArrayType, IdArrayType> PermuteIndexType;

  // create a temporary array the same size
  vtkm::cont::ArrayHandle<vtkm::Id> newActiveVertices;

  // Use only the current activeVertices outdegree to match size on StreamCompact
  vtkm::cont::ArrayHandle<vtkm::Id> outdegreeLookup;
  DeviceAlgorithm::Copy(PermuteIndexType(activeVertices, outdegree), outdegreeLookup);

  // compact the activeVertices array to keep only the ones of interest
  DeviceAlgorithm::StreamCompact(activeVertices, outdegreeLookup, newActiveVertices);

  activeVertices.ReleaseResources();
  DeviceAlgorithm::Copy(newActiveVertices, activeVertices);

#ifdef DEBUG_PRINT
  DebugPrint("Active Vertex List Compacted");
#endif
} // CompactActiveVertices()

// compacts the active edge list
template<typename T, typename StorageType, typename DeviceAdapter>
void ChainGraph<T,StorageType,DeviceAdapter>::CompactActiveEdges()
{
#ifdef DEBUG_FUNCTION_ENTRY
		std::cout << std::endl;
		std::cout << "====================" << std::endl;
		std::cout << "Compact Active Edges" << std::endl;
		std::cout << "====================" << std::endl;
		std::cout << std::endl;
#endif
  // grab the size of the array for easier reference
  vtkm::Id nActiveVertices = activeVertices.GetNumberOfValues();

  // first, we have to work out the first edge for each active vertex
  // we start with a temporary new updegree
  vtkm::cont::ArrayHandle<vtkm::Id> newOutdegree;
  newOutdegree.Allocate(nActiveVertices);

  // do a parallel computation using the vertex degree updater
  // WARNING: Using chainMaximum for I/O in parallel loop
  // See functor description for algorithmic justification of safety
  VertexDegreeUpdater vertexDegreeUpdater;
  vtkm::worklet::DispatcherMapField<VertexDegreeUpdater>
                 vertexDegreeUpdaterDispatcher(vertexDegreeUpdater);

  vertexDegreeUpdaterDispatcher.Invoke(activeVertices,        // input
                                       activeEdges,           // input (whole array)
                                       edgeFar,               // input (whole array)
                                       firstEdge,             // input (whole array)
                                       prunesTo,              // input (whole array)
                                       outdegree,             // input (whole array)
                                       chainExtremum,         // i/o (whole array)
                                       newOutdegree);         // output

  // now we do a reduction to compute the offsets of each vertex
  vtkm::cont::ArrayHandle<vtkm::Id> newPosition;
  DeviceAlgorithm::ScanExclusive(newOutdegree, newPosition);
  vtkm::Id nNewEdges = newPosition.GetPortalControl().Get(nActiveVertices-1) + 
                       newOutdegree.GetPortalControl().Get(nActiveVertices-1);

  // create a temporary vector for copying
  vtkm::cont::ArrayHandle<vtkm::Id> newActiveEdges;
  newActiveEdges.Allocate(nNewEdges);

  // now copy the relevant edges into the active edge array
  // WARNING: Using chainMaximum, edgeHigh, firstEdge, updegree for I/O in parallel loop
  // See functor description for algorithmic justification of safety
  ActiveEdgeTransferrer<DeviceAdapter> activeEdgeTransferrer(
                                       activeEdges.PrepareForInput(DeviceAdapter()),
                                       prunesTo.PrepareForInput(DeviceAdapter()));
  vtkm::worklet::DispatcherMapField<ActiveEdgeTransferrer<DeviceAdapter> > 
                 activeEdgeTransferrerDispatcher(activeEdgeTransferrer);

  activeEdgeTransferrerDispatcher.Invoke(activeVertices,       // input
                                         newPosition,          // input
                                         newOutdegree,         // input
                                         firstEdge,            // i/o (whole array)
                                         outdegree,            // i/o (whole array)
                                         chainExtremum,        // i/o (whole array)
                                         edgeFar,              // i/o (whole array)
                                         newActiveEdges);      // output (whole array)

  // resize the original array and recopy
  DeviceAlgorithm::Copy(newActiveEdges, activeEdges);

#ifdef DEBUG_PRINT
  DebugPrint("Active Edges Now Compacted");
#endif
} // CompactActiveEdges()
	
// builds the chains for the new active vertices
template<typename T, typename StorageType, typename DeviceAdapter>
void ChainGraph<T,StorageType,DeviceAdapter>::BuildChains()
{
#ifdef DEBUG_FUNCTION_ENTRY
		std::cout << std::endl;
		std::cout << "============" << std::endl;
		std::cout << "Build Chains" << std::endl;
		std::cout << "============" << std::endl;
		std::cout << std::endl;
#endif
  // a temporary array the full size of the graph
  vtkm::cont::ArrayHandle<vtkm::Id> tempChainExtremum;
  tempChainExtremum.Allocate(edgeNear.GetNumberOfValues());
	
  // compute the number of log steps required in this pass
  vtkm::Id nActiveVertices = activeVertices.GetNumberOfValues();
  vtkm::Id nLogSteps = 1;
  for (vtkm::Id shifter = nActiveVertices; shifter != 0; shifter >>= 1)
    nLogSteps++;

  ChainDoubler chainDoubler;
  vtkm::worklet::DispatcherMapField<ChainDoubler> chainDoublerDispatcher(chainDoubler);

  // 2.	Use path compression / step doubling to collect vertices along ascending chains
  //		until every vertex has been assigned to *an* extremum
  //		Step two at a time, so that we rock between the original and the temp
  for (vtkm::Id logStep = 0; logStep < nLogSteps; logStep++)
  {
    chainDoublerDispatcher.Invoke(activeVertices,  // input
                                  chainExtremum);  // i/o (whole array)
  }
		
#ifdef DEBUG_PRINT
  DebugPrint("Chains Built");
#endif
} // BuildChains()

// transfers saddle ascent edges into edge sorter
template<typename T, typename StorageType, typename DeviceAdapter>
void ChainGraph<T,StorageType,DeviceAdapter>::TransferSaddleStarts()
{
#ifdef DEBUG_FUNCTION_ENTRY
		std::cout << std::endl;
		std::cout << "=======================" << std::endl;
		std::cout << DEBUG_STRING_TRANSFER_SADDLE_STARTS << std::endl;
		std::cout << "=======================" << std::endl;
		std::cout << std::endl;
#endif

  // grab the size of the array for easier reference
  vtkm::Id nActiveVertices = activeVertices.GetNumberOfValues();

  // reset number of edges to sort
  vtkm::Id nEdgesToSort = 0;

  // in parallel, we need to create a vector to count the first edge for each vertex
  vtkm::cont::ArrayHandle<vtkm::Id> newFirstEdge;
  vtkm::cont::ArrayHandle<vtkm::Id> newOutdegree;
  newFirstEdge.Allocate(nActiveVertices);
  newOutdegree.Allocate(nActiveVertices);

  // 2. now test all active vertices to see if they have only one chain maximum
  SaddleAscentFunctor saddleAscentFunctor;
  vtkm::worklet::DispatcherMapField<SaddleAscentFunctor> 
                 saddleAscentFunctorDispatcher(saddleAscentFunctor);

  saddleAscentFunctorDispatcher.Invoke(activeVertices,           // input
                                       firstEdge,                // input (whole array)
                                       outdegree,                // input (whole array)
                                       activeEdges,              // input (whole array)
                                       chainExtremum,            // input (whole array)
                                       edgeFar,                  // input (whole array)
                                       newOutdegree);            // output
 
  // 3. now compute the new offsets in the newFirstEdge array
  DeviceAlgorithm::ScanExclusive(newOutdegree, newFirstEdge);
  nEdgesToSort = newFirstEdge.GetPortalControl().Get(nActiveVertices-1) + 
                 newOutdegree.GetPortalControl().Get(nActiveVertices-1);


  edgeSorter.ReleaseResources();
  edgeSorter.Allocate(nEdgesToSort);

  SaddleAscentTransferrer saddleAscentTransferrer;
  vtkm::worklet::DispatcherMapField<SaddleAscentTransferrer> 
                 saddleAscentTransferrerDispatcher(saddleAscentTransferrer);

  saddleAscentTransferrerDispatcher.Invoke(activeVertices,       // input
                                           newOutdegree,         // input
                                           newFirstEdge,         // input
                                           activeEdges,          // input (whole array)
                                           firstEdge,            // input (whole array)
                                           edgeSorter);          // output (whole array)

#ifdef DEBUG_PRINT
  DebugPrint(DEBUG_STRING_TRANSFERRED_SADDLE_STARTS);
#endif
} // TransferSaddleStarts()

// sets all remaining active vertices
template<typename T, typename StorageType, typename DeviceAdapter>
void ChainGraph<T,StorageType,DeviceAdapter>::BuildTrunk()
{
#ifdef DEBUG_FUNCTION_ENTRY
		std::cout << std::endl;
		std::cout << "===========" << std::endl;
		std::cout << "Build Trunk" << std::endl;
		std::cout << "============" << std::endl;
		std::cout << std::endl;
#endif

  TrunkBuilder trunkBuilder;
  vtkm::worklet::DispatcherMapField<TrunkBuilder> 
                 trunkBuilderDispatcher(trunkBuilder);

  trunkBuilderDispatcher.Invoke(activeVertices,       // input
                                chainExtremum,        // input (whole array)
                                prunesTo);            // output (whole array)
#ifdef DEBUG_PRINT
  DebugPrint("Trunk Built");
#endif
} // BuildTrunk()

// transfers partial results to merge tree array
template<typename T, typename StorageType, typename DeviceAdapter>
void ChainGraph<T,StorageType,DeviceAdapter>::TransferToMergeTree(vtkm::cont::ArrayHandle<vtkm::Id> &saddles)
{
#ifdef DEBUG_FUNCTION_ENTRY
		std::cout << std::endl;
		std::cout << "=====================" << std::endl;
		std::cout << DEBUG_STRING_TRANSFER_TO_MERGE_TREE << std::endl;
		std::cout << "=====================" << std::endl;
		std::cout << std::endl;
#endif

  // first allocate memory for the target array
  saddles.ReleaseResources();

  // initialise it to the arcArray
  DeviceAlgorithm::Copy(arcArray, saddles);

  JoinTreeTransferrer joinTreeTransferrer;
  vtkm::worklet::DispatcherMapField<JoinTreeTransferrer> 
                 joinTreeTransferrerDispatcher(joinTreeTransferrer);
  vtkm::cont::ArrayHandleIndex valueIndexArray(valueIndex.GetNumberOfValues());

  joinTreeTransferrerDispatcher.Invoke(valueIndexArray,    // input
                                       prunesTo,           // input
                                       valueIndex,         // input (whole array)
                                       chainExtremum,      // input (whole array)
                                       saddles,            // output (whole array)
                                       arcArray);          // output (whole array)
} // TransferToMergeTree()

// prints the contents of the topology graph in standard format
template<typename T, typename StorageType, typename DeviceAdapter>
void ChainGraph<T,StorageType,DeviceAdapter>::DebugPrint(const char *message)
{
  std::cout << "---------------------------" << std::endl;
  std::cout << std::string(message) << std::endl;
  std::cout << "---------------------------" << std::endl;
  std::cout << std::endl;

  typedef vtkm::cont::ArrayHandle<vtkm::Id> IdArrayType;
  typedef vtkm::cont::ArrayHandle<T> ValueArrayType;
  typedef vtkm::cont::ArrayHandlePermutation<IdArrayType, IdArrayType> PermuteIndexType;
  typedef vtkm::cont::ArrayHandlePermutation<IdArrayType, ValueArrayType> PermuteValueType;

  // Full Vertex Arrays
  vtkm::Id nValues = valueIndex.GetNumberOfValues();
  vtkm::cont::ArrayHandle<T,StorageType> vertexValues;

  std::cout << "Full Vertex Arrays - Size:  " << nValues << std::endl;
  printHeader(nValues);
  printIndices("Index", valueIndex);
  DeviceAlgorithm::Copy(PermuteValueType(valueIndex, values), vertexValues);
  printValues("Value", vertexValues);
  printIndices("First Edge", firstEdge);
  printIndices("Outdegree", outdegree);
  printIndices("Chain Ext", chainExtremum);
  printIndices("Prunes To", prunesTo);
  std::cout << std::endl;

  // Active Vertex Arrays
  vtkm::Id nActiveVertices = activeVertices.GetNumberOfValues();
  std::cout << "Active Vertex Arrays - Size: " << nActiveVertices << std::endl;
  if (nActiveVertices > 0) {
    vtkm::cont::ArrayHandle<vtkm::Id> tempIndex;
    vtkm::cont::ArrayHandle<T> tempValue;

    printHeader(nActiveVertices);
    printIndices("Active Vertices", activeVertices);
    DeviceAlgorithm::Copy(PermuteIndexType(activeVertices, valueIndex), tempIndex);
    printIndices("Active Indices", tempIndex);
    DeviceAlgorithm::Copy(PermuteValueType(activeVertices, vertexValues), tempValue);
    printValues("Active Values", tempValue);
    DeviceAlgorithm::Copy(PermuteIndexType(activeVertices, firstEdge), tempIndex);
    printIndices("Active First Edge", tempIndex);
    DeviceAlgorithm::Copy(PermuteIndexType(activeVertices, outdegree), tempIndex);
    printIndices("Active Outdegree", tempIndex);
    DeviceAlgorithm::Copy(PermuteIndexType(activeVertices, chainExtremum), tempIndex);
    printIndices("Active Chain Ext", tempIndex);
    DeviceAlgorithm::Copy(PermuteIndexType(activeVertices, prunesTo), tempIndex);
    printIndices("Active Prunes To", tempIndex);
    std::cout << std::endl;
  }

  // Full Edge Arrays
  vtkm::Id nEdges = edgeNear.GetNumberOfValues();
  std::cout << "Full Edge Arrays - Size:     " << nEdges << std::endl;
  vtkm::cont::ArrayHandle<vtkm::Id> farIndices;
  vtkm::cont::ArrayHandle<vtkm::Id> nearIndices;
  vtkm::cont::ArrayHandle<T,StorageType> farValues;
  vtkm::cont::ArrayHandle<T,StorageType> nearValues;
  if (nEdges > 0) {
    printHeader(nEdges);
    printIndices("Far", edgeFar);
    DeviceAlgorithm::Copy(PermuteIndexType(edgeFar, valueIndex), farIndices);
    printIndices("Far Index", farIndices);
    DeviceAlgorithm::Copy(PermuteValueType(farIndices, values), farValues);
    printValues("Far Value", farValues);

    printHeader(nEdges);
    printIndices("Near", edgeNear);
    DeviceAlgorithm::Copy(PermuteIndexType(edgeNear, valueIndex), nearIndices);
    printIndices("Near Index", nearIndices);
    DeviceAlgorithm::Copy(PermuteValueType(nearIndices, values), nearValues);
    printValues("Near Value", nearValues);
  }

  // Active Edge Arrays
  vtkm::Id nActiveEdges = activeEdges.GetNumberOfValues();
  std::cout << "Active Edge Arrays - Size:   " << nActiveEdges << std::endl;
  if (nActiveEdges > 0) {
    vtkm::cont::ArrayHandle<vtkm::Id> activeFarIndices;
    vtkm::cont::ArrayHandle<vtkm::Id> activeNearIndices;
    vtkm::cont::ArrayHandle<vtkm::Id> activeNearLookup;
    vtkm::cont::ArrayHandle<T,StorageType> activeNearValues;

    printHeader(nActiveEdges);
    printIndices("Active Edges", activeEdges);

    DeviceAlgorithm::Copy(PermuteIndexType(activeEdges, edgeFar), activeFarIndices);
    printIndices("Edge Far", activeFarIndices);
    DeviceAlgorithm::Copy(PermuteIndexType(activeEdges, edgeNear), activeNearIndices);
    printIndices("Edge Near", activeNearIndices);
    DeviceAlgorithm::Copy(PermuteIndexType(activeNearIndices, valueIndex), activeNearLookup);
    printIndices("Edge Near Index", activeNearLookup);
    DeviceAlgorithm::Copy(PermuteValueType(activeNearLookup, values), activeNearValues);
    printValues("Edge Near Value", activeNearValues);
    std::cout << std::endl;
  }

  // Edge Sorter Array
  vtkm::Id nEdgeSorter = edgeSorter.GetNumberOfValues();
  std::cout << "Edge Sorter - Size:          " << nEdgeSorter << std::endl;
  if (nEdgeSorter > 0) {
    vtkm::cont::ArrayHandle<vtkm::Id> tempSortIndex;
    vtkm::cont::ArrayHandle<T> tempSortValue;

    printHeader(nEdgeSorter);
    printIndices("Edge Sorter", edgeSorter);
    DeviceAlgorithm::Copy(PermuteIndexType(edgeSorter, edgeNear), tempSortIndex);
    printIndices("Sorted Near", tempSortIndex);
    DeviceAlgorithm::Copy(PermuteIndexType(edgeSorter, nearIndices), tempSortIndex);
    printIndices("Sorted Near Index", tempSortIndex);
    DeviceAlgorithm::Copy(PermuteIndexType(edgeSorter, edgeFar), tempSortIndex);
    printIndices("Sorted Far", tempSortIndex);
    DeviceAlgorithm::Copy(PermuteValueType(edgeSorter, nearValues), tempSortValue);
    printValues("Sorted Near Value", tempSortValue);
    std::cout << std::endl;
  }

  std::cout << "---------------------------" << std::endl;
  std::cout << std::endl;
} // DebugPrint()

}
}
}

#endif