#include <cassert>
#include <cstdlib>
#include <iostream>
#include <vector>

#include "cctk.h"
#include "cctk_Parameters.h"

#include "defs.hh"
#include "dh.hh"
#include "th.hh"
#include "vect.hh"

#include "gh.hh"

using namespace std;



  // Constructors
gh::
gh (vector<ivect> const & reffacts_, centering const refcent_,
    int const mgfact_,  centering const mgcent_,
    vector<vector<ibbox> > const & baseextents_,
    i2vect const & boundary_width_)
  : reffacts(reffacts_), refcent(refcent_),
    mgfact(mgfact_), mgcent(mgcent_),
    baseextents(baseextents_),
    boundary_width(boundary_width_)
{
  assert (reffacts.size() >= 1);
  assert (all (reffacts.front() == 1));
  for (int rl = 1; rl < (int)reffacts.size(); ++ rl) {
    assert (all (reffacts.AT(rl) >= reffacts.AT(rl-1)));
    assert (all (reffacts.AT(rl) % reffacts.AT(rl-1) == 0));
  }
  assert (refcent == vertex_centered or refcent == cell_centered);
  
  assert (mgfact >= 1);
  assert (mgcent == vertex_centered or mgcent == cell_centered);
  
  assert (baseextents.size() >= 1);
  assert (baseextents.AT(0).size() >= 1);
  assert (baseextents.AT(0).size() == reffacts.size());
  for (int ml = 1; ml < (int)baseextents.size(); ++ ml) {
    assert (baseextents.AT(ml).size() == baseextents.AT(ml-1).size());
  }
  for (int ml = 0; ml < (int)baseextents.size(); ++ ml) {
    for (int rl = 1; rl < (int)baseextents.AT(ml).size(); ++ rl) {
      ibbox const & cbox = baseextents.AT(ml).AT(rl-1);
      ibbox const & fbox = baseextents.AT(ml).AT(rl);
      assert (all (cbox.stride() * reffacts.AT(rl-1) == 
                   fbox.stride() * reffacts.AT(rl)));
    }
  }
  
  assert (all (all (boundary_width >= 0)));
  for (int ml = 0; ml < (int)baseextents.size(); ++ ml) {
    for (int rl = 0; rl < (int)baseextents.AT(ml).size(); ++ rl) {
      ibbox const & box = baseextents.AT(ml).AT(rl);
      // This condition must hold even for zero-sized grid arrays
      assert (all (box.shape() / box.stride() >=
                   boundary_width[0] + boundary_width[1]));
    }
  }
}

// Destructors
gh::
~gh ()
{
}



// Modifiers
void
gh::
regrid (mregs const & regs)
{
  DECLARE_CCTK_PARAMETERS;
  
  // Save the grid hierarchy
  oldregions.clear ();
  swap (oldregions, regions);
  regions = regs;
  
  
  
  // Consistency checks
  
  // Note: there may be zero refinement levels
  
  // Check multigrid consistency
  assert (mglevels()>0);
  for (int ml=1; ml<mglevels(); ++ml) {
    for (int rl=0; rl<reflevels(); ++rl) {
      for (int c=0; c<components(rl); ++c) {
	assert (all(extent(ml,rl,c).stride() ==
                    mgfact * extent(ml-1,rl,c).stride()));
	assert (extent(ml,rl,c)
		.contracted_for(extent(ml-1,rl,c))
		.is_contained_in(extent(ml-1,rl,c)));
      }
    }
  }
  
  // Check component consistency
  for (int ml=0; ml<mglevels(); ++ml) {
    for (int rl=0; rl<reflevels(); ++rl) {
      assert (components(rl)>0);
      for (int c=0; c<components(rl); ++c) {
        ibbox const & b  = extent(ml,rl,c);
        ibbox const & b0 = extent(ml,rl,0);
	assert (all(b.stride() == b0.stride()));
	assert (b.is_aligned_with(b0));
        for (int cc=c+1; cc<components(rl); ++cc) {
          assert ((b & extent(ml,rl,cc)).empty());
        }
      }
    }
  }
  
  // Check base grid extent
  for (int ml=0; ml<mglevels(); ++ml) {
    for (int rl=0; rl<reflevels(); ++rl) {
      for (int c=0; c<components(rl); ++c) {
        assert (extent(ml,rl,c).is_contained_in(baseextents.at(ml).at(rl)));
      }
    }
  }

  // Check proper nesting, i.e., whether the fine grids are contained
  // in the coarser grids
  {
    bool have_error = false;
    for (int ml=0; ml<mglevels(); ++ml) {
      for (int rl=1; rl<reflevels(); ++rl) {
        assert (all (extent(ml,rl,0).stride() * reffacts.AT(rl) ==
                     extent(ml,rl-1,0).stride() * reffacts.AT(rl-1)));
        // Check contained-ness:
        // first take all coarse grids
        ibset coarse;
        for (int c=0; c<components(rl-1); ++c) {
          coarse += extent(ml,rl-1,c);
        }
        coarse.normalize();
        // then check all fine grids
        for (int c=0; c<components(rl); ++c) {
          ibbox const & fine =
            extent(ml,rl,c).contracted_for(extent(ml,rl-1,0));
          if (not (fine <= coarse)) {
            if (not have_error) {
              cout << "The following components are not properly nested, i.e.," << endl
                   << "they are not contained within the next coarser level's components:" << endl;
              have_error = true;
            }
            cout << "   ml " << ml << " rl " << rl << " c " << c << ":   "
                 << fine << endl;
          }
        } // for c
      }   // for rl
    }     // for ml
    if (have_error) {
      cout << "The grid hierarchy is:" << endl;
      for (int ml=0; ml<mglevels(); ++ml) {
        for (int rl=0; rl<reflevels(); ++rl) {
          for (int c=0; c<components(rl); ++c) {
            cout << "   ml " << ml << " rl " << rl << " c " << c << ":   "
                 << extent(ml,rl,c) << endl;
          }
        }
      }
      CCTK_WARN (0, "The refinement hierarchy is not properly nested.");
    }
  }
  
  // Output
  if (output_bboxes) {
    do_output_bboxes (cout);
    do_output_bases (cout);
  }
  
  // Recompose the other hierarchies
  for (list<th*>::iterator t=ths.begin(); t!=ths.end(); ++t) {
    (*t)->regrid();
  }
  
  for (list<dh*>::iterator d=dhs.begin(); d!=dhs.end(); ++d) {
    (*d)->regrid();
  }
}



bool
gh::
recompose (int const rl,
           bool const do_prolongate)
{
  // Handle changes in number of mglevels
  if (oldregions.size() != regions.size()) {
    oldregions.resize (regions.size());
  }
  
  bool const do_recompose = level_did_change(rl);
  
  if (do_recompose) {
    
    // Recompose the other hierarchies
    for (list<dh*>::iterator d=dhs.begin(); d!=dhs.end(); ++d) {
      (*d)->recompose (rl, do_prolongate);
    }
    
    // Overwrite old with new grid hierarchy
    for (int ml=0; ml<mglevels(); ++ml) {
      oldregions.AT(ml).resize (regions.AT(ml).size());
      oldregions.AT(ml).AT(rl) = regions.AT(ml).AT(rl);
    }
    
  }
  
  return do_recompose;
}



bool
gh::
level_did_change (int const rl)
  const
{
  // Find out whether this level changed
  if (regions.size() != oldregions.size()) return true;
  for (int ml=0; ml<mglevels(); ++ml) {
    assert (rl>=0 and rl<reflevels());
    if (rl >= (int)oldregions.AT(ml).size()) return true;
    if (regions.AT(ml).AT(rl).size() != oldregions.AT(ml).AT(rl).size()) {
      return true;
    }
    for (int c=0; c<components(rl); ++c) {
      if (regions.AT(ml).AT(rl).AT(c) != oldregions.AT(ml).AT(rl).AT(c)) {
        return true;
      }
    } // for c
  } // for ml
  return false;
}



// Accessors

int
gh::
local_components (int const rl)
  const
{
  int lc = 0;
  for (int c=0; c<components(rl); ++c) {
    lc +=  is_local(rl,c);
  }
  return lc;
}



// Time hierarchy management

void
gh::
add (th * const t)
{
  ths.push_back (t);
}

void
gh::
remove (th * const t)
{
  ths.remove (t);
}



// Data hierarchy management

void
gh::
add (dh * const d)
{
  dhs.push_back (d);
}

void
gh::
remove (dh * const d)
{
  dhs.remove (d);
}



// Memory usage

size_t
gh::
memory ()
  const
{
  return
    memoryof (reffacts) +
    memoryof (refcent) +
    memoryof (mgfact) +
    memoryof (mgcent) +
    memoryof (baseextents) +
    memoryof (boundary_width) +
    memoryof (regions) +
    memoryof (oldregions) +
    memoryof (ths) +
    memoryof (dhs);
}



// Output

void
gh::
do_output_bboxes (ostream & os)
  const
{
  for (int ml=0; ml<mglevels(); ++ml) {
    for (int rl=0; rl<reflevels(); ++rl) {
      for (int c=0; c<components(rl); ++c) {
        os << eol;
        os << "gh bboxes:" << eol;
        os << "ml=" << ml << " rl=" << rl << " c=" << c << eol;
        os << "extent=" << extent(ml,rl,c) << eol;
        os << "outer_boundaries=" << outer_boundaries(rl,c) << eol;
        os << "processor=" << processor(rl,c) << eol;
      }
    }
  }
}

void
gh::
do_output_bases (ostream & os)
  const
{
  for (int ml=0; ml<mglevels(); ++ml) {
    for (int rl=0; rl<reflevels(); ++rl) {
      if (components(rl)>0) {
        os << eol;
        os << "gh bases:" << eol;
        os << "ml=" << ml << " rl=" << rl << eol;
        os << "base=" << baseextents.AT(ml).AT(rl) << eol;
      }
    }
  }
}

ostream &
gh::
output (ostream & os)
  const
{
  os << "gh:"
     << "reffacts=" << reffacts << ",refcentering=" << refcent << ","
     << "mgfactor=" << mgfact << ",mgcentering=" << mgcent << ","
     << "baseextents=" << baseextents << ","
     << "boundary_width=" << boundary_width << ","
     << "regions=" << regions << ","
     << "dhs={";
  {
    bool isfirst = true;
    for (list<dh*>::const_iterator
           d = dhs.begin(); d != dhs.end(); ++ d, isfirst = false)
    {
      if (not isfirst) os << ",";
      os << *d;
    }
  }
  os << "}";
  return os;
}