void eval_integrand(UserVector & output, std::vector<Scalar> & input) {
   int    dimension = (int)alpha.size();
   Scalar total     = 0.0;
   Scalar point     = 0.0;
   for (int i=0; i<dimension; i++) {
     point     = 0.5*input[i]+0.5;
     total    += powl(alpha[i]*(point-beta[i]),(long double)2.0);
   }
   output.clear(); output.resize(1,std::exp(-total));
 }  

 Scalar error_indicator(UserVector & input) {
   int dimension = (int)input.size();
   Scalar norm2  = 0.0;
   for (int i=0; i<dimension; i++)
     norm2 += input[i]*input[i];
   
   Scalar ID = AdaptiveSparseGridInterface<Scalar,UserVector>::
     getInitialDiff();
   norm2 = std::sqrt(norm2)/ID;
   return norm2;
 }

void AdaptiveSparseGridInterface<Scalar,UserVector>::eval_cubature(
	       UserVector & output, 
	       CubatureTensorSorted<Scalar> & cubRule) {

  //int dimf      = 0;                      // Dimension of the integrand
  Scalar weight = 0.0;
  std::vector<Scalar> point(dimension_,(Scalar)0.0);
  //std::vector<Scalar> f(1,0.0);
  Teuchos::RCP<UserVector> f = output.Create(); output.Update(-1.0,output);

  typename std::map<std::vector<Scalar>,int>::iterator it;
  for (it=cubRule.begin(); it!=cubRule.end(); it++) {
    // Evaluate Function
    point.assign((it->first).begin(),(it->first).end()); // Extract point
    f->Update(-1.0,*f);
    eval_integrand(*f,point);     // Evaluate Integrand at point
 
    // Update integral
    weight = cubRule.getWeight(it->second);
    output.Update(weight,*f);
  }
}

void AdaptiveSparseGridInterface<Scalar,UserVector>::init(UserVector & output) {
  std::vector<int> index(dimension_,1);
  CubatureTensorSorted<Scalar> cubRule(
            dimension_,index,rule1D_,growth1D_,isNormalized_);
  
  // Evaluate the initial contribution to the integral
  initialDiff_ = 1.0;
  output.Update(-1.0,output);
  eval_cubature(output,cubRule);

  // Compute the initial error indicator
  initialDiff_ = error_indicator(output);
  if (fabs(initialDiff_)<INTREPID_TOL) 
    initialDiff_ = 1.0;
}

Scalar AdaptiveSparseGrid<Scalar,UserVector>::refine_grid(
	typename std::multimap<Scalar,std::vector<int> > & activeIndex, 
	std::set<std::vector<int> > & oldIndex, 
	UserVector & integralValue,
	CubatureTensorSorted<Scalar> & cubRule,
	Scalar globalErrorIndicator,
	AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {

  TEUCHOS_TEST_FOR_EXCEPTION((activeIndex.empty()),std::out_of_range,
              ">>> ERROR (AdaptiveSparseGrid): Active Index set is empty.");  

  int dimension = problem_data.getDimension();
  std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
  std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);

  // Initialize Flags
  bool maxLevelFlag     = true;
  bool isAdmissibleFlag = true;

  // Initialize Cubature Rule
  Teuchos::RCP<UserVector> s = integralValue.Create();

  // Initialize iterator at end of inOldIndex
  std::set<std::vector<int> >::iterator it1(oldIndex.end());  

  // Initialize iterator at end of inActiveIndex
  typename std::multimap<Scalar,std::vector<int> >::iterator it;

  // Obtain Global Error Indicator as sum of key values of inActiveIndex
  Scalar eta = globalErrorIndicator;

  // Select Index to refine
  it = activeIndex.end(); it--;        // Decrement to position of final value 
  Scalar G               = it->first;  // Largest Error Indicator is at End
  eta                   -= G;          // Update global error indicator
  std::vector<int> index = it->second; // Get Corresponding index
  activeIndex.erase(it);               // Erase Index from active index set
  // Insert Index into old index set
  oldIndex.insert(it1,index); it1 = oldIndex.end(); 
  
  // Refinement process
  for (int k=0; k<dimension; k++) {
    index[k]++; // index + ek
    // Check Max Level
    maxLevelFlag = problem_data.max_level(index);
    if (maxLevelFlag) {
      // Check Admissibility
      isAdmissibleFlag = isAdmissible(index,k,oldIndex,problem_data);
      if (isAdmissibleFlag) { // If admissible
	// Build Differential Quarature Rule
	CubatureTensorSorted<Scalar> diffRule(0,dimension);
	build_diffRule(diffRule,index,problem_data);
	
	// Apply Rule to function
	problem_data.eval_cubature(*s,diffRule);
	
	// Update integral value
	integralValue.Update(*s);
	
	// Update local error indicator and index set
	G  = problem_data.error_indicator(*s); 	
	if (activeIndex.end()!=activeIndex.begin()) 
	  activeIndex.insert(activeIndex.end()--,
			   std::pair<Scalar,std::vector<int> >(G,index));
	else
	  activeIndex.insert(std::pair<Scalar,std::vector<int> >(G,index));
		
	// Update global error indicators
	eta += G;

	// Update adapted quadrature rule nodes and weights
	cubRule.update(1.0,diffRule,1.0);
      }
    }
    else { // Max Level Exceeded 
      //std::cout << "Maximum Level Exceeded" << std::endl;
    }
    index[k]--;
  }
  return eta;
}

Scalar AdaptiveSparseGrid<Scalar,UserVector>::refine_grid(
   	         typename std::multimap<Scalar,std::vector<int> >  & indexSet, 
		 UserVector & integralValue,
		 AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
  
  int dimension = problem_data.getDimension();
  std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
  std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);
  
  // Copy Multimap into a Set for ease of use
  typename std::multimap<Scalar,std::vector<int> >::iterator it;
  std::set<std::vector<int> > oldSet;
  std::set<std::vector<int> >::iterator it1(oldSet.begin());
  for (it=indexSet.begin(); it!=indexSet.end(); it++) {
    oldSet.insert(it1,it->second);
    it1++;
  }
  indexSet.clear();
  
  // Find Possible Active Points
  int flag = 1;
  std::vector<int> index(dimension,0);
  typename std::multimap<Scalar,std::vector<int> > activeSet;
  for (it1=oldSet.begin(); it1!=oldSet.end(); it1++) {
    index = *it1;
    for (int i=0; i<dimension; i++) {
      index[i]++;
      flag = (int)(!oldSet.count(index));
      index[i]--;
      if (flag) {
	activeSet.insert(std::pair<Scalar,std::vector<int> >(1.0,index));
	oldSet.erase(it1);
	break;
      }
    }
  }

  // Compute local and global error indicators for active set
  typename std::multimap<Scalar,std::vector<int> >::iterator it2;
  Scalar eta = 0.0;
  Scalar G   = 0.0;
  Teuchos::RCP<UserVector> s = integralValue.Create();
  for (it2=activeSet.begin(); it2!=activeSet.end(); it2++) {
    // Build Differential Quarature Rule
    index = it2->second;
    CubatureTensorSorted<Scalar> diffRule(0,dimension);
    build_diffRule(diffRule,index,problem_data);
    
    // Apply Rule to function
    problem_data.eval_cubature(*s,diffRule);
    
    // Update local error indicator and index set
    G  = problem_data.error_indicator(*s);
    activeSet.erase(it2);
    activeSet.insert(it2,std::pair<Scalar,std::vector<int> >(G,index));
    eta += G;
  }

  // Refine Sparse Grid
  eta = refine_grid(activeSet,oldSet,integralValue,eta,
		    dimension,rule1D,growth1D);

  // Insert New Active and Old Index sets into indexSet
  indexSet.insert(activeSet.begin(),activeSet.end());
  for (it1=oldSet.begin(); it1!=oldSet.end(); it1++) {
    index = *it1;
    indexSet.insert(std::pair<Scalar,std::vector<int> >(-1.0,index));
  }
    
  return eta;
}

  void eval_integrand(UserVector & output, std::vector<Scalar> & input) {
    output.clear(); output.resize(1,std::exp(-input[0]*input[0])
				  +10.0*std::exp(-input[1]*input[1]));
  }  

 void eval_integrand(UserVector & output, std::vector<Scalar> & input) {
   output.clear(); output.resize(1,powl(input[0]+input[1],(long double)6.0));
 }