#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/ParallelComm.hpp"
#include "moab/Skinner.hpp"
#include "moab/CN.hpp"
#include "MBTagConventions.hpp"

#include <iostream>
#include <iterator>
#include <list>
#include <algorithm>
#include <sstream>
#include <limits>
#include <string>
#include <time.h>

#include "moab_mpi.h"

using namespace moab;
const unsigned int _dim = 2;
char file[4096] = "./square_10x10_100p_quad4.h5m";

ErrorCode check_ghost_face_ids(Interface& mb);

int main( int argc, char* argv[] )
{
  MPI_Init( &argc, &argv );
  Core moab;
  Interface& instance = moab;

  std::ostringstream opts;
  opts << "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION";
  ErrorCode bval = instance.load_file( file, 0, opts.str().c_str());
  if (MB_SUCCESS != bval) {
    std::cerr << file << " : file read failed (" << instance.get_error_string(bval) << ")" << std::endl;
    return 1;
  }
  
  ErrorCode result;
  result = check_ghost_face_ids(instance);
    assert(MB_SUCCESS==result);
  
  MPI_Finalize();
  return result;
}

ErrorCode check_ghost_face_ids(Interface& instance)
{
  ErrorCode result;
  ParallelComm* pcomm = ParallelComm::get_pcomm(&instance, 0); 
  result = pcomm->resolve_shared_ents( 0, _dim); 
    assert(MB_SUCCESS==result);
  std::vector<int> numFaces(pcomm->size());
  Range shared_faces;
  result = pcomm->get_shared_entities(-1, shared_faces, _dim-1, false, false);
    assert(MB_SUCCESS==result);
  int n_faces = shared_faces.size();
  if (pcomm->size() > 1) 
  {
    int retval = MPI_Allgather(&n_faces, 1, MPI_INT,
                               &numFaces[0], 1, 
                               MPI_INT, pcomm->comm());
    assert (retval == 0);
  }
  int total_elems = 0;
  for (unsigned int proc = 0; proc < pcomm->rank(); proc++) 
    total_elems += numFaces[proc];
  Tag gid_tag;
  int zero = 0;
  result = instance.tag_get_handle("GLOBAL_ID", 1, MB_TYPE_INTEGER, 
                                  gid_tag, MB_TAG_DENSE|MB_TAG_CREAT, &zero);
  assert(MB_SUCCESS==result);
  numFaces.resize(n_faces);
  for (unsigned int i = 0; i < n_faces; i++)
  {
    numFaces[i] = total_elems++;
  }
  unsigned int j = 0;
  for (moab::Range::iterator faceiter = shared_faces.begin(); faceiter != shared_faces.end(); ++faceiter)
  {
    result = instance.tag_set_data(gid_tag, &(*faceiter), 1, &numFaces[j]); 
      assert(MB_SUCCESS==result);
    j++;
  }
  
  std::vector<int> partitionFacesGlobalGhostIndices(n_faces);
  result = pcomm->exchange_tags(gid_tag, shared_faces);
    assert(MB_SUCCESS==result);
  j=0;
  for (moab::Range::iterator faceiter = shared_faces.begin(); faceiter != shared_faces.end(); ++faceiter)
  {
    result = instance.tag_get_data(gid_tag, &(*faceiter), 1, &partitionFacesGlobalGhostIndices[j]); 
      assert(MB_SUCCESS==result);
//    if (pcomm->rank() == 1)
//      std::cout<<" rank "<<pcomm->rank()<<" new id "<<partitionFacesGlobalGhostIndices[j]<<" old id "<<numFaces[j]<<std::endl;
    j++;
  }

  // once again with reduce_tags in place of exchange_tags
  j = 0;
  for (moab::Range::iterator faceiter = shared_faces.begin(); faceiter != shared_faces.end(); ++faceiter)
  {
    result = instance.tag_set_data(gid_tag, &(*faceiter), 1, &numFaces[j]); 
      assert(MB_SUCCESS==result);
    j++;
  }
  std::vector<int> partitionFacesGlobalGhostIndices2(n_faces);
  Tag gid_sum_tag;
  result = instance.tag_get_handle("__global_id_sum", 1, MB_TYPE_INTEGER, 
                                  gid_sum_tag, MB_TAG_DENSE|MB_TAG_CREAT, &zero);
    assert(MB_SUCCESS==result);
  std::vector<Tag> gid_tag_vec;
  gid_tag_vec.push_back(gid_tag);
  std::vector<Tag> gid_sum_tag_vec;
  gid_sum_tag_vec.push_back(gid_sum_tag);
  result = pcomm->reduce_tags(gid_tag_vec, gid_sum_tag_vec, MPI_SUM, shared_faces);
    assert(MB_SUCCESS==result);
  j=0;
  for (moab::Range::iterator faceiter = shared_faces.begin(); faceiter != shared_faces.end(); ++faceiter)
  {
    int global_id;
    int global_id_sum;
    result = instance.tag_get_data(gid_tag, &(*faceiter), 1, &global_id); 
      assert(MB_SUCCESS==result);
    result = instance.tag_get_data(gid_sum_tag, &(*faceiter), 1, &global_id_sum); 
      assert(MB_SUCCESS==result);
    partitionFacesGlobalGhostIndices2[j] = global_id_sum - global_id;
//    if (pcomm->rank() == 1)
//      std::cout<<" rank "<<pcomm->rank()<<" new id "<<partitionFacesGlobalGhostIndices2[j]<<" old id "<<numFaces[j]<<std::endl;
    j++;
  }
  //check
  j=0;
  for (moab::Range::iterator faceiter = shared_faces.begin(); faceiter != shared_faces.end(); ++faceiter)
  {
    unsigned char pstatus;
    pcomm->get_pstatus(*faceiter, pstatus);
    if (pstatus & PSTATUS_NOT_OWNED)
    {
      if (partitionFacesGlobalGhostIndices2[j] != partitionFacesGlobalGhostIndices[j])
	std::cout<<" error, global indices do not match on face "<<j<<" old id "<<numFaces[j]<<" new exchanged id "<<partitionFacesGlobalGhostIndices[j]<<" new reduced id "<<partitionFacesGlobalGhostIndices2[j]<<std::endl;
    }
    j++;
  }
  return result;
}