bool compare_rank_2_views(const array_type& y,
                          const array_type& y_exp,
                          const scalar_type rel_tol,
                          const scalar_type abs_tol,
                          Teuchos::FancyOStream& out)
{
  typedef typename array_type::size_type size_type;
  typename array_type::HostMirror hy = Kokkos::create_mirror_view(y);
  typename array_type::HostMirror hy_exp = Kokkos::create_mirror_view(y_exp);
  Kokkos::deep_copy(hy, y);
  Kokkos::deep_copy(hy_exp, y_exp);

  size_type num_rows = y.dimension_0();
  size_type num_cols = y.dimension_1();
  bool success = true;
  for (size_type i=0; i<num_rows; ++i) {
    for (size_type j=0; j<num_cols; ++j) {
      scalar_type diff = std::abs( hy(i,j) - hy_exp(i,j) );
      scalar_type tol = rel_tol*std::abs(hy_exp(i,j)) + abs_tol;
      bool s = diff < tol;
      out << "y_expected(" << i << "," << j << ") - "
          << "y(" << i << "," << j << ") = " << hy_exp(i,j)
          << " - " << hy(i,j) << " == "
          << diff << " < " << tol << " : ";
      if (s)
        out << "passed";
      else
        out << "failed";
      out << std::endl;
      success = success && s;
    }
  }

  return success;
}

 static void decode(const array_type& _array, float& _t) {
   float _sum = 0.0;
   for (int i = 0; i < _array.size(); ++i) {
     _sum += (_array[i] / 256.0) * std::pow(256.0,float(i-int(_array.size()/2)));
   }
   _t = _sum;
 }

 static void encode(const float& _t, array_type& _array) {
   for (int i = 0; i < _array.size(); ++i) {
     float _n = _t / std::pow(256.0,float(i-int(_array.size()/2)));
     float _intpart;
     float _fractpart = std::modf (_n , &_intpart);
     _array[i] = uint8_t(256.0*_fractpart);
   }
 }

            virtual void test()
            {
                m_nPushError = 0;

                for ( array_type::const_iterator it = m_arr.begin(); it != m_arr.end(); ++it ) {
                    if ( !m_Queue.push( SimpleValue( *it ) ))
                        ++m_nPushError;
                }
            }

	Integrate(	array_type 			& arg_output , 
				const array_type	& arg_left ,
				const array_type	& arg_right  ) : output(arg_output) , left(arg_left) , right(arg_right)
	{
		numLeft = left.dimension(1);
		numRight = right.dimension(1);
		numPoints = left.dimension(2);
		dim = left.dimension(3);
		if(output.rank() == 2) numLeft = 1;
	}

    key_type key(const array_type & df, const std::size_t index) const
    {
        if (df.shape().second == 0)
        {
            return {0, 0};
        }
        else
        {
            const int prod_id = df[df.row(index)][0]; // TODO, hardcoded
            const char segment = df[df.row(index)][8]; // TODO, hardcoded

            return {prod_id, segment};
        }
    }

    inline void TrBDF2Scheme<TrapezoidalScheme>::step(array_type& fn, Time t) {
        using namespace ext::placeholders;

        QL_REQUIRE(t-dt_ > -1e-8, "a step towards negative time given");

        const Time intermediateTimeStep = dt_*alpha_;

        array_type fStar = fn;
        trapezoidalScheme_->setStep(intermediateTimeStep);
        trapezoidalScheme_->step(fStar, t);

        bcSet_.setTime(std::max(0.0, t-dt_));
        bcSet_.applyBeforeSolving(*map_, fn);

        const array_type f =
            (1/alpha_*fStar - square<Real>()(1-alpha_)/alpha_*fn)/(2-alpha_);

        if (map_->size() == 1) {
            fn = map_->solve_splitting(0, f, -beta_);
        }
        else {
            const ext::function<Disposable<Array>(const Array&)>
                preconditioner(ext::bind(
                    &FdmLinearOpComposite::preconditioner, map_, _1, -beta_));

            const ext::function<Disposable<Array>(const Array&)> applyF(
                ext::bind(&TrBDF2Scheme<TrapezoidalScheme>::apply, this, _1));

            if (solverType_ == BiCGstab) {
                const BiCGStabResult result =
                    QuantLib::BiCGstab(applyF, std::max(Size(10), fn.size()),
                        relTol_, preconditioner).solve(f, f);

                (*iterations_) += result.iterations;
                fn = result.x;
            } else if (solverType_ == GMRES) {
                const GMRESResult result =
                    QuantLib::GMRES(applyF, std::max(Size(10), fn.size()/10u),
                        relTol_, preconditioner).solve(f, f);

                (*iterations_) += result.errors.size();
                fn = result.x;
            }
            else
                QL_FAIL("unknown/illegal solver type");
        }

        bcSet_.applyAfterSolving(fn);
    }

    inline void TRBDF2<Operator>::step(array_type& a, Time t) {
        Size i;
        Array aInit(a.size());
        for (i=0; i<a.size();i++) {
            aInit[i] = a[i];
        }
        aInit_ = aInit;
        for (i=0; i<bcs_.size(); i++)
            bcs_[i]->setTime(t);
        //trapezoidal explicit part
        if (L_.isTimeDependent()) {
            L_.setTime(t);
            explicitTrapezoidalPart_ = I_ - 0.5*alpha_*dt_*L_;
        }
        for (i=0; i<bcs_.size(); i++)
            bcs_[i]->applyBeforeApplying(explicitTrapezoidalPart_);
        a = explicitTrapezoidalPart_.applyTo(a);
        for (i=0; i<bcs_.size(); i++)
            bcs_[i]->applyAfterApplying(a);

        // trapezoidal implicit part
        if (L_.isTimeDependent()) {
            L_.setTime(t-dt_);
            implicitPart_ = I_ + 0.5*alpha_*dt_*L_;
        }
        for (i=0; i<bcs_.size(); i++)
            bcs_[i]->applyBeforeSolving(implicitPart_,a);
        a = implicitPart_.solveFor(a);
        for (i=0; i<bcs_.size(); i++)
            bcs_[i]->applyAfterSolving(a);


        // BDF2 explicit part
        if (L_.isTimeDependent()) {
            L_.setTime(t);
        }
        for (i=0; i<bcs_.size(); i++) {
            bcs_[i]->applyBeforeApplying(explicitBDF2PartFull_);
        }
        array_type b0 = explicitBDF2PartFull_.applyTo(aInit_);
        for (i=0; i<bcs_.size(); i++)
            bcs_[i]->applyAfterApplying(b0);

        for (i=0; i<bcs_.size(); i++) {
            bcs_[i]->applyBeforeApplying(explicitBDF2PartMid_);
        }
        array_type b1 = explicitBDF2PartMid_.applyTo(a);
        for (i=0; i<bcs_.size(); i++)
            bcs_[i]->applyAfterApplying(b1);
        a = b0+b1;

        // reuse implicit part - works only for alpha=2-sqrt(2)
        for (i=0; i<bcs_.size(); i++)
            bcs_[i]->applyBeforeSolving(implicitPart_,a);
        a = implicitPart_.solveFor(a);
        for (i=0; i<bcs_.size(); i++)
            bcs_[i]->applyAfterSolving(a);

    }

std::unique_ptr<void, int (*)(BoosterHandle)>
fit(const array_type & train_data,
    const std::vector<float> & train_y,
    const std::map<const std::string, const std::string> & params,
    _StopCondition stop_condition)
{
    // prepare placeholder for raw matrix later used by xgboost
    std::vector<float> train_vec = train_data.tovector();
    std::cerr << "train_vec size: " << train_vec.size() << std::endl;
//    assert(std::none_of(train_vec.cbegin(), train_vec.cend(), [](float x){return std::isnan(x);}));

    std::unique_ptr<void, int (*)(DMatrixHandle)> tr_dmat(
        XGDMatrixCreateFromMat(
            train_vec.data(),
            train_data.shape().first,
            train_data.shape().second, XGB_MISSING),
        XGDMatrixFree);

    // attach response vector to tr_dmat
    XGDMatrixSetFloatInfo(tr_dmat.get(), "label", train_y.data(), train_y.size());

    const DMatrixHandle cache[] = {tr_dmat.get()};

    // create Booster with attached tr_dmat
    std::unique_ptr<void, int (*)(BoosterHandle)> booster(
            XGBoosterCreate(cache, 1UL),
            XGBoosterFree);

    for (const auto & kv : params)
    {
        std::cerr << kv.first << " => " << kv.second << std::endl;
        XGBoosterSetParam(booster.get(), kv.first.c_str(), kv.second.c_str());
    }


    for (int iter{0}; stop_condition() == false; ++iter)
    {
        XGBoosterUpdateOneIter(booster.get(), iter, tr_dmat.get());
    }

    return booster;
}

//
//
//
/// @param [in] input A previously initialized array with the same rank order.
//
/// @brief The operator @c + applied between two arrays of same rank, add up
/// each element and returns the resulting array by move assignment or move
/// construction. Only if both arrays were previously initialized either by
/// construction or the array::create() member function. Otherwise nothing is
/// done.
//
/// @return An array with the same rank order.
//
array_type operator +(const array_type &input)
{
	array_type output;
	output.init_as(this);
	OMP_STATIC_LOOP_POLICY
	for(INIT_ITER(i, 0); i < output.data_length() && i < input.data_length(); ++i)
	{
		output.data[i] = data[i] + input.data[i];
	}
	return output;
};

    void ImplicitEulerScheme::step(array_type& a, Time t) {
        using namespace ext::placeholders;
        QL_REQUIRE(t-dt_ > -1e-8, "a step towards negative time given");
        map_->setTime(std::max(0.0, t-dt_), t);
        bcSet_.setTime(std::max(0.0, t-dt_));

        bcSet_.applyBeforeSolving(*map_, a);

        if (map_->size() == 1) {
            a = map_->solve_splitting(0, a, -dt_);
        }
        else {
            const ext::function<Disposable<Array>(const Array&)>
                preconditioner(ext::bind(
                    &FdmLinearOpComposite::preconditioner, map_, _1, -dt_));

            const ext::function<Disposable<Array>(const Array&)> applyF(
                ext::bind(&ImplicitEulerScheme::apply, this, _1));

            if (solverType_ == BiCGstab) {
                const BiCGStabResult result =
                    QuantLib::BiCGstab(applyF, std::max(Size(10), a.size()),
                        relTol_, preconditioner).solve(a, a);

                (*iterations_) += result.iterations;
                a = result.x;
            }
            else if (solverType_ == GMRES) {
                const GMRESResult result =
                    QuantLib::GMRES(applyF, std::max(Size(10), a.size()/10u),
                        relTol_, preconditioner).solve(a, a);

                (*iterations_) += result.errors.size();
                a = result.x;
            }
            else
                QL_FAIL("unknown/illegal solver type");
        }
        bcSet_.applyAfterSolving(a);
    }

std::vector<float>
predict(
    BoosterHandle booster,
    const array_type & test_data)
{
    std::vector<float> test_vec = test_data.tovector();
    std::cerr << "test_vec size: " << test_vec.size() << std::endl;

    std::unique_ptr<void, int (*)(DMatrixHandle)> te_dmat(
        XGDMatrixCreateFromMat(
            test_vec.data(),
            test_data.shape().first,
            test_data.shape().second, XGB_MISSING),
        XGDMatrixFree);

    bst_ulong y_hat_len{0};
    const float * y_hat_proba{nullptr};
    XGBoosterPredict(booster, te_dmat.get(), 0, 0, &y_hat_len, &y_hat_proba);
    std::cerr << "Got y_hat_proba of length " << y_hat_len << std::endl;

    std::vector<float> y_hat(y_hat_proba, y_hat_proba + y_hat_len);

    return y_hat;
}

	Transform(	array_type 			& arg_output , 
				const array_type	& arg_input ,
				const array_type	& arg_fields  ) : output(arg_output) , input(arg_input) , fields(arg_fields)
	{
		data_rank = input.rank();
		numDataPts = input.dimension(1);
		in_rank = fields.rank();
		numCells = output.dimension(0);
		numFields = output.dimension(1);
		numPoints = output.dimension(2);
		dim = output.dimension(3);		
	}

    void ImplicitEulerScheme::step(array_type& a, Time t) {
        QL_REQUIRE(t-dt_ > -1e-8, "a step towards negative time given");
        map_->setTime(std::max(0.0, t-dt_), t);
        bcSet_.setTime(std::max(0.0, t-dt_));

        bcSet_.applyBeforeSolving(*map_, a);

        a = BiCGstab(
                boost::function<Disposable<Array>(const Array&)>(
                    boost::bind(&ImplicitEulerScheme::apply, this, _1)), 
                10*a.size(), relTol_,
                boost::function<Disposable<Array>(const Array&)>(
                    boost::bind(&FdmLinearOpComposite::preconditioner, 
                                map_, _1, -dt_))
            ).solve(a).x;
        
        bcSet_.applyAfterSolving(a);
    }

  // Verify:
  KOKKOS_INLINE_FUNCTION
  void operator()( size_t iwork, value_type & errors ) const
  {
    const size_t tile_dim0 = ( m_array.dimension_0() + TileLayout::N0 - 1 ) / TileLayout::N0;
    const size_t tile_dim1 = ( m_array.dimension_1() + TileLayout::N1 - 1 ) / TileLayout::N1;

    const size_t itile = iwork % tile_dim0;
    const size_t jtile = iwork / tile_dim0;

    if ( jtile < tile_dim1 ) {
      tile_type tile = Kokkos::tile_subview( m_array, itile, jtile );

      if ( tile( 0, 0 ) != ptrdiff_t( ( itile + jtile * tile_dim0 ) * TileLayout::N0 * TileLayout::N1 ) ) {
        ++errors;
      }
      else {
        for ( size_t j = 0; j < size_t( TileLayout::N1 ); ++j ) {
          for ( size_t i = 0; i < size_t( TileLayout::N0 ); ++i ) {
            const size_t iglobal = i + itile * TileLayout::N0;
            const size_t jglobal = j + jtile * TileLayout::N1;

            if ( iglobal < m_array.dimension_0() && jglobal < m_array.dimension_1() ) {
              if ( tile( i, j ) != ptrdiff_t( tile( 0, 0 ) + i + j * TileLayout::N0 ) ) ++errors;

              //printf( "tile(%d, %d)(%d, %d) = %d\n", int( itile ), int( jtile ), int( i ), int( j ), int( tile( i, j ) ) );
            }
          }
        }
      }
    }
  }

 void prepare( size_t nStart, size_t nEnd )
 {
     m_arr.reserve( nEnd - nStart );
     for ( size_t i = nStart; i < nEnd; ++i )
         m_arr.push_back( i );
     std::random_shuffle( m_arr.begin(), m_arr.end() );
 }

 // Kernel
 void run()
 {
     dist = sqrt(
           (pow(points1, 2).rowwise() * alpha).rowwise().sum().replicate(1, nsample2).rowwise()
         + (pow(points2, 2).rowwise() * alpha).rowwise().sum().transpose()
         - 2 * ((points1.rowwise() * alpha).matrix() * points2.matrix().transpose()).array()
     );
 }

		//-----------------------------------------------------------------//
		handle_set(bool zhe = true, uint32_t fas = 0) : array_(), erase_set_(),
														zero_handle_enable_(zhe) {
			if(fas) {
				array_.reserve(fas);
				array_.clear();
			}
			if(zero_handle_enable_) array_.push_back(T());
		}

		//-----------------------------------------------------------------//
		handle_type insert(const T& st) {
			handle_type h;
			if(erase_set_.empty()) {
				h = static_cast<handle_type>(array_.size());
				array_.push_back(st);
			} else {
				set_it it = erase_set_.begin();
				h = *it;
				array_[h] = st;
				erase_set_.erase(it);
			}
			return h;
		}

	void insert_right(const char c) {
		auto insert_point = s.end() - 1;

		if( justify == Justify::Left ) {
			insert_point = std::find_if(s.begin(), s.end(), [](const char& a) {
				return a == ' ';
			});
		}

		if( *insert_point != ' ' ) {
			insert_point = shift_left();
		}

		*insert_point = c;
	}