Thyra_DefaultMultiPeriodModelEvaluator.hpp

// @HEADER
// *****************************************************************************
//    Thyra: Interfaces and Support for Abstract Numerical Algorithms
//
// Copyright 2004 NTESS and the Thyra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef THYRA_DEFAULT_MULTI_PERIOD_MODEL_EVALUATOR_HPP
#define THYRA_DEFAULT_MULTI_PERIOD_MODEL_EVALUATOR_HPP


#include "Thyra_ModelEvaluatorDefaultBase.hpp"
#include "Thyra_DefaultProductVectorSpace.hpp"
#include "Thyra_DefaultBlockedTriangularLinearOpWithSolveFactory.hpp" // Default implementation
#include "Thyra_DefaultBlockedLinearOp.hpp"
#include "Thyra_ModelEvaluatorDelegatorBase.hpp"
#include "Thyra_ProductVectorBase.hpp"
#include "Teuchos_implicit_cast.hpp"
#include "Teuchos_AbstractFactory.hpp" // Interface
#include "Teuchos_AbstractFactoryStd.hpp" // Implementation


namespace Thyra {


template<class Scalar>
class DefaultMultiPeriodModelEvaluator
  : virtual public ModelEvaluatorDefaultBase<Scalar>
{
public:


  DefaultMultiPeriodModelEvaluator();

  DefaultMultiPeriodModelEvaluator(
    const int N,
    const Array<RCP<ModelEvaluator<Scalar> > > &periodModels,
    const Array<int> &z_indexes,
    const Array<Array<RCP<const VectorBase<Scalar> > > > &z,
    const int g_index,
    const Array<Scalar> g_weights,
    const RCP<const ProductVectorSpaceBase<Scalar> > &x_bar_space = Teuchos::null,
    const RCP<const ProductVectorSpaceBase<Scalar> > &f_bar_space = Teuchos::null,
    const RCP<LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory = Teuchos::null
    );

  void initialize(
    const int N,
    const Array<RCP<ModelEvaluator<Scalar> > > &periodModels,
    const Array<int> &z_indexes,
    const Array<Array<RCP<const VectorBase<Scalar> > > > &z,
    const int g_index,
    const Array<Scalar> g_weights,
    const RCP<const ProductVectorSpaceBase<Scalar> > &x_bar_space = Teuchos::null,
    const RCP<const ProductVectorSpaceBase<Scalar> > &f_bar_space = Teuchos::null,
    const RCP<LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory = Teuchos::null
    );

  void reset_z(
    const Array<Array<RCP<const VectorBase<Scalar> > > > &z
    );
  


  int Np() const;
  int Ng() const;
  RCP<const VectorSpaceBase<Scalar> > get_x_space() const;
  RCP<const VectorSpaceBase<Scalar> > get_f_space() const;
  RCP<const VectorSpaceBase<Scalar> > get_p_space(int l) const;
  RCP<const Array<std::string> > get_p_names(int l) const;
  RCP<const VectorSpaceBase<Scalar> > get_g_space(int j) const;
  ArrayView<const std::string> get_g_names(int j) const;
  ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
  ModelEvaluatorBase::InArgs<Scalar> getLowerBounds() const;
  ModelEvaluatorBase::InArgs<Scalar> getUpperBounds() const;
  RCP<LinearOpBase<Scalar> > create_W_op() const;
  RCP<PreconditionerBase<Scalar> > create_W_prec() const;
  RCP<const LinearOpWithSolveFactoryBase<Scalar> > get_W_factory() const;
  ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
  void reportFinalPoint(
    const ModelEvaluatorBase::InArgs<Scalar> &finalPoint,
    const bool wasSolved
    );


private:



  RCP<LinearOpBase<Scalar> > create_DfDp_op_impl(int l) const;
  RCP<LinearOpBase<Scalar> > create_DgDx_dot_op_impl(int j) const;
  RCP<LinearOpBase<Scalar> > create_DgDx_op_impl(int j) const;
  RCP<LinearOpBase<Scalar> > create_DgDp_op_impl(int j, int l) const;
  ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
  void evalModelImpl(
    const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
    const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
    ) const;


private:

  // //////////////////////////
  // Private types

  typedef Array<Scalar> g_weights_t;
  typedef Array<Array<RCP<const VectorBase<Scalar> > > > z_t;

  // /////////////////////////
  // Private data members

  RCP<ModelEvaluator<Scalar> > periodModel_;
  Array<RCP<ModelEvaluator<Scalar> > > periodModels_;
  Array<int> z_indexes_;
  Array<int> period_l_map_;
  z_t z_; // size == N
  int g_index_;
  g_weights_t g_weights_; // size == N
  RCP<const ProductVectorSpaceBase<Scalar> > x_bar_space_;
  RCP<const ProductVectorSpaceBase<Scalar> > f_bar_space_;
  RCP<LinearOpWithSolveFactoryBase<Scalar> > W_bar_factory_;
  int Np_;
  int Ng_;
  ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
  ModelEvaluatorBase::InArgs<Scalar> lowerBounds_;
  ModelEvaluatorBase::InArgs<Scalar> upperBounds_;

  // /////////////////////////
  // Private member functions

  void set_z_indexes_and_create_period_l_map( const Array<int> &z_indexes );

  void wrapNominalValuesAndBounds();

  static
  RCP<ProductVectorBase<Scalar> >
  createProductVector(
    const RCP<const ProductVectorSpaceBase<Scalar> > &prodVecSpc
    );

  // Return the index of a "free" parameter in the period model given its
  // index in this mulit-period model.
  int period_l(const int l) const
    {
      TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= l && l < Np_) );
      return period_l_map_[l];
    }

  int numPeriodZs() const { return z_indexes_.size(); }

  int N() const { return x_bar_space_->numBlocks(); }
  
};


// /////////////////////////////////
// Implementations


// Constructors/Intializers/Accessors


template<class Scalar>
DefaultMultiPeriodModelEvaluator<Scalar>::DefaultMultiPeriodModelEvaluator()
  :g_index_(-1), Np_(-1), Ng_(-1)
{}


template<class Scalar>
DefaultMultiPeriodModelEvaluator<Scalar>::DefaultMultiPeriodModelEvaluator(
  const int N,
  const Array<RCP<ModelEvaluator<Scalar> > > &periodModels,
  const Array<int> &z_indexes,
  const Array<Array<RCP<const VectorBase<Scalar> > > > &z,
  const int g_index,
  const Array<Scalar> g_weights,
  const RCP<const ProductVectorSpaceBase<Scalar> > &x_bar_space,
  const RCP<const ProductVectorSpaceBase<Scalar> > &f_bar_space,
  const RCP<LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory
  )
  :g_index_(-1), Np_(-1), Ng_(-1)
{
  initialize(
    N, periodModels, z_indexes, z, g_index, g_weights,
    x_bar_space, f_bar_space, W_bar_factory
    );
}


template<class Scalar>
void DefaultMultiPeriodModelEvaluator<Scalar>::initialize(
  const int N,
  const Array<RCP<ModelEvaluator<Scalar> > > &periodModels,
  const Array<int> &z_indexes,
  const Array<Array<RCP<const VectorBase<Scalar> > > > &z,
  const int g_index,
  const Array<Scalar> g_weights,
  const RCP<const ProductVectorSpaceBase<Scalar> > &x_bar_space,
  const RCP<const ProductVectorSpaceBase<Scalar> > &f_bar_space,
  const RCP<LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory
  )
{

  using Teuchos::implicit_cast;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef ModelEvaluatorBase MEB;

  TEUCHOS_TEST_FOR_EXCEPT( N <= 0 );
  TEUCHOS_TEST_FOR_EXCEPT( implicit_cast<int>(periodModels.size()) != N );
  MEB::InArgs<Scalar> periodInArgs = periodModels[0]->createInArgs();
  MEB::OutArgs<Scalar> periodOutArgs = periodModels[0]->createOutArgs();
  for ( int k = 0; k < implicit_cast<int>(z_indexes.size()); ++k ) {
    TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= z_indexes[k] && z_indexes[k] < periodInArgs.Np() ) );
  }
  TEUCHOS_TEST_FOR_EXCEPT( implicit_cast<int>(z.size())!=N );
  TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= g_index && g_index < periodOutArgs.Ng() ) );
  TEUCHOS_TEST_FOR_EXCEPT( implicit_cast<int>(g_weights.size())!=N );
  // ToDo: Assert that the different period models are compatible!

  Np_ = periodInArgs.Np() - z_indexes.size();

  Ng_ = 1;

  set_z_indexes_and_create_period_l_map(z_indexes);
  
  periodModel_ = periodModels[0]; // Assume basic structure!
  
  periodModels_ = periodModels;
  
  z_.resize(N);
  reset_z(z);
  
  g_index_ = g_index;
  g_weights_ = g_weights;
  
  if ( periodInArgs.supports(MEB::IN_ARG_x) ) {
    if( !is_null(x_bar_space) ) {
      TEUCHOS_TEST_FOR_EXCEPT(!(x_bar_space->numBlocks()==N));
      // ToDo: Check the constituent spaces more carefully against models[]->get_x_space().
      x_bar_space_ = x_bar_space;
    }
    else {
      x_bar_space_ = productVectorSpace(periodModel_->get_x_space(),N);
      // ToDo: Update to build from different models for different x spaces!
    }
  }

  if ( periodOutArgs.supports(MEB::OUT_ARG_f) ) {
    if( !is_null(f_bar_space) ) {
      TEUCHOS_TEST_FOR_EXCEPT(!(f_bar_space->numBlocks()==N));
      // ToDo: Check the constituent spaces more carefully against models[]->get_f_space().
      f_bar_space_ = f_bar_space;
    }
    else {
      f_bar_space_ = productVectorSpace(periodModel_->get_f_space(),N);
      // ToDo: Update to build from different models for different f spaces!
    }
  }

  if ( periodOutArgs.supports(MEB::OUT_ARG_W) ) {
    if ( !is_null(W_bar_factory) ) {
      // ToDo: Check the compatability of the W_bar objects created using this object!
      W_bar_factory_ = W_bar_factory;
    }
    else {
      W_bar_factory_ =
        defaultBlockedTriangularLinearOpWithSolveFactory<Scalar>(
          periodModel_->get_W_factory() );
    }
  }

  wrapNominalValuesAndBounds();
  
}


template<class Scalar>
void DefaultMultiPeriodModelEvaluator<Scalar>::reset_z(
  const Array<Array<RCP<const VectorBase<Scalar> > > > &z
  )
{

  using Teuchos::implicit_cast;
  
  const int N = z_.size();

#ifdef TEUCHOS_DEBUG
  TEUCHOS_TEST_FOR_EXCEPT( N == 0 && "Error, must have called initialize() first!" );
  TEUCHOS_TEST_FOR_EXCEPT( implicit_cast<int>(z.size()) != N ); 
#endif 

  for( int i = 0; i < N; ++i ) {
    const Array<RCP<const VectorBase<Scalar> > >  &z_i = z[i];
#ifdef TEUCHOS_DEBUG
    TEUCHOS_TEST_FOR_EXCEPT( z_i.size() != z_indexes_.size() ); 
#endif 
    z_[i] = z_i;
  }
  
}


// Public functions overridden from ModelEvaulator


template<class Scalar>
int DefaultMultiPeriodModelEvaluator<Scalar>::Np() const
{
  return Np_;
}


template<class Scalar>
int DefaultMultiPeriodModelEvaluator<Scalar>::Ng() const
{
  return Ng_;
}


template<class Scalar>
RCP<const VectorSpaceBase<Scalar> >
DefaultMultiPeriodModelEvaluator<Scalar>::get_x_space() const
{
  return x_bar_space_;
}


template<class Scalar>
RCP<const VectorSpaceBase<Scalar> >
DefaultMultiPeriodModelEvaluator<Scalar>::get_f_space() const
{
  return f_bar_space_;
}


template<class Scalar>
RCP<const VectorSpaceBase<Scalar> >
DefaultMultiPeriodModelEvaluator<Scalar>::get_p_space(int l) const
{
  return  periodModel_->get_p_space(period_l(l));
}


template<class Scalar>
RCP<const Array<std::string> >
DefaultMultiPeriodModelEvaluator<Scalar>::get_p_names(int l) const
{
  return  periodModel_->get_p_names(period_l(l));
}


template<class Scalar>
RCP<const VectorSpaceBase<Scalar> >
DefaultMultiPeriodModelEvaluator<Scalar>::get_g_space(int j) const
{
  TEUCHOS_TEST_FOR_EXCEPT(j!=0);
  return periodModel_->get_g_space(g_index_);
}


template<class Scalar>
ArrayView<const std::string>
DefaultMultiPeriodModelEvaluator<Scalar>::get_g_names(int j) const
{
  return  periodModel_->get_g_names(j);
}


template<class Scalar>
ModelEvaluatorBase::InArgs<Scalar>
DefaultMultiPeriodModelEvaluator<Scalar>::getNominalValues() const
{
  return nominalValues_;
}


template<class Scalar>
ModelEvaluatorBase::InArgs<Scalar>
DefaultMultiPeriodModelEvaluator<Scalar>::getLowerBounds() const
{
  return lowerBounds_;
}


template<class Scalar>
ModelEvaluatorBase::InArgs<Scalar>
DefaultMultiPeriodModelEvaluator<Scalar>::getUpperBounds() const
{
  return upperBounds_;
}


template<class Scalar>
RCP<LinearOpBase<Scalar> >
DefaultMultiPeriodModelEvaluator<Scalar>::create_W_op() const
{
  // Set up the block structure ready to have the blocks filled!
  const RCP<PhysicallyBlockedLinearOpBase<Scalar> >
    W_op_bar = defaultBlockedLinearOp<Scalar>();
  W_op_bar->beginBlockFill(f_bar_space_,x_bar_space_);
  const int N = x_bar_space_->numBlocks();
  for ( int i = 0; i < N; ++i ) {
    W_op_bar->setNonconstBlock( i, i, periodModel_->create_W_op() );
  }
  W_op_bar->endBlockFill();
  return W_op_bar;
}


template<class Scalar>
RCP<PreconditionerBase<Scalar> >
DefaultMultiPeriodModelEvaluator<Scalar>::create_W_prec() const
{
  return Teuchos::null;
}


template<class Scalar>
RCP<const LinearOpWithSolveFactoryBase<Scalar> >
DefaultMultiPeriodModelEvaluator<Scalar>::get_W_factory() const
{
  return W_bar_factory_;
}


template<class Scalar>
ModelEvaluatorBase::InArgs<Scalar>
DefaultMultiPeriodModelEvaluator<Scalar>::createInArgs() const
{
  typedef ModelEvaluatorBase MEB;
  MEB::InArgs<Scalar> periodInArgs = periodModel_->createInArgs();
  MEB::InArgsSetup<Scalar> inArgs;
  inArgs.setModelEvalDescription(this->description());
  inArgs.set_Np(Np_);
  inArgs.setSupports( MEB::IN_ARG_x, periodInArgs.supports(MEB::IN_ARG_x) );
  return inArgs;
}


template<class Scalar>
void DefaultMultiPeriodModelEvaluator<Scalar>::reportFinalPoint(
  const ModelEvaluatorBase::InArgs<Scalar>      &finalPoint
  ,const bool                                   wasSolved
  )
{
  // We are just going to ignore the final point here.  It is not clear how to
  // report a "final" point back to the underlying *periodModel_ object since
  // we have so many different "points" that we could return (i.e. one for
  // each period).  I guess we could report back the final parameter values
  // (other than the z parameter) but there are multiple states x[i] and
  // period parameters z[i] that we can report back.
}


// Public functions overridden from ModelEvaulatorDefaultBase


template<class Scalar>
RCP<LinearOpBase<Scalar> >
DefaultMultiPeriodModelEvaluator<Scalar>::create_DfDp_op_impl(int l) const
{
  TEUCHOS_TEST_FOR_EXCEPT("This class does not support DfDp(l) as a linear operator yet.");
  return Teuchos::null;
}


template<class Scalar>
RCP<LinearOpBase<Scalar> >
DefaultMultiPeriodModelEvaluator<Scalar>::create_DgDx_dot_op_impl(int j) const
{
  TEUCHOS_TEST_FOR_EXCEPT("This class does not support DgDx_dot(j) as a linear operator yet.");
  return Teuchos::null;
}


template<class Scalar>
RCP<LinearOpBase<Scalar> >
DefaultMultiPeriodModelEvaluator<Scalar>::create_DgDx_op_impl(int j) const
{
  TEUCHOS_TEST_FOR_EXCEPT("This class does not support DgDx(j) as a linear operator yet.");
  return Teuchos::null;
}


template<class Scalar>
RCP<LinearOpBase<Scalar> >
DefaultMultiPeriodModelEvaluator<Scalar>::create_DgDp_op_impl(int j, int l) const
{
  TEUCHOS_TEST_FOR_EXCEPT("This class does not support DgDp(j,l) as a linear operator yet.");
  return Teuchos::null;
}


template<class Scalar>
ModelEvaluatorBase::OutArgs<Scalar>
DefaultMultiPeriodModelEvaluator<Scalar>::createOutArgsImpl() const
{

  typedef ModelEvaluatorBase MEB;

  MEB::OutArgs<Scalar> periodOutArgs = periodModel_->createOutArgs();
  MEB::OutArgsSetup<Scalar> outArgs;

  outArgs.setModelEvalDescription(this->description());

  outArgs.set_Np_Ng(Np_,Ng_);

  // f
  if (periodOutArgs.supports(MEB::OUT_ARG_f) ) {
    outArgs.setSupports(MEB::OUT_ARG_f);
  }

  // W_op
  if (periodOutArgs.supports(MEB::OUT_ARG_W_op) ) {
    outArgs.setSupports(MEB::OUT_ARG_W_op);
    outArgs.set_W_properties(periodOutArgs.get_W_properties());
  }
  // Note: We will not directly support the LOWSB form W as we will let the
  // default base class handle this given our W_factory!

  // DfDp(l)
  for ( int l = 0; l < Np_; ++l ) {
    const int period_l = this->period_l(l);
    const MEB::DerivativeSupport period_DfDp_l_support
      = periodOutArgs.supports(MEB::OUT_ARG_DfDp,period_l);
    if (!period_DfDp_l_support.none()) {
      outArgs.setSupports( MEB::OUT_ARG_DfDp, l, period_DfDp_l_support );
      outArgs.set_DfDp_properties(
        l, periodOutArgs.get_DfDp_properties(period_l) );
    }
  }

  // DgDx_dot
  const MEB::DerivativeSupport
    period_DgDx_dot_support = periodOutArgs.supports(MEB::OUT_ARG_DgDx_dot,g_index_);
  if (!period_DgDx_dot_support.none()) {
    outArgs.setSupports( MEB::OUT_ARG_DgDx_dot, 0, period_DgDx_dot_support );
    outArgs.set_DgDx_dot_properties(
      0, periodOutArgs.get_DgDx_dot_properties(g_index_) );
  }

  // DgDx
  const MEB::DerivativeSupport
    period_DgDx_support = periodOutArgs.supports(MEB::OUT_ARG_DgDx,g_index_);
  if (!period_DgDx_support.none()) {
    outArgs.setSupports( MEB::OUT_ARG_DgDx, 0, period_DgDx_support );
    outArgs.set_DgDx_properties(
      0, periodOutArgs.get_DgDx_properties(g_index_) );
  }

  // DgDp(l)
  for ( int l = 0; l < Np_; ++l ) {
    const int period_l = this->period_l(l);
    const MEB::DerivativeSupport period_DgDp_l_support
      = periodOutArgs.supports(MEB::OUT_ARG_DgDp, g_index_, period_l);
    if (!period_DgDp_l_support.none()) {
      outArgs.setSupports( MEB::OUT_ARG_DgDp, 0, l, period_DgDp_l_support );
      outArgs.set_DgDp_properties(
        0, l, periodOutArgs.get_DgDp_properties(g_index_,period_l) );
    }
  }
  
  return outArgs;

}


template<class Scalar>
void DefaultMultiPeriodModelEvaluator<Scalar>::evalModelImpl(
  const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
  const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
  ) const
{

  using Teuchos::rcp_dynamic_cast;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef ModelEvaluatorBase MEB;
  typedef Teuchos::VerboseObjectTempState<ModelEvaluatorBase> VOTSME;

  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
    "DefaultMultiPeriodModelEvaluator",inArgs,outArgs,periodModel_ );
  // ToDo: You will have to set the verbosity level for each of the
  // periodModels_[i] individually below!
  
  const int N = x_bar_space_->numBlocks();
  const int Np = this->Np_;
  //const int Ng = this->Ng_;

  //
  // A) Setup InArgs
  //

  RCP<const ProductVectorBase<Scalar> > x_bar;
  if (inArgs.supports(MEB::IN_ARG_x)) {
    x_bar = rcp_dynamic_cast<const ProductVectorBase<Scalar> >(
      inArgs.get_x(), true );
    TEUCHOS_TEST_FOR_EXCEPTION(
      is_null(x_bar), std::logic_error,
      "Error, if x is supported, it must be set!"
      );
  }

  //
  // B) Setup OutArgs
  //
  
  RCP<ProductVectorBase<Scalar> > f_bar;
  if (outArgs.supports(MEB::OUT_ARG_f)) {
    f_bar = rcp_dynamic_cast<ProductVectorBase<Scalar> >(
      outArgs.get_f(), true );
  }

  Array<MEB::Derivative<Scalar> > DfDp_bar(Np);
  Array<RCP<ProductMultiVectorBase<Scalar> > > DfDp_bar_mv(Np);
  for ( int l = 0; l < Np; ++l ) {
    if (!outArgs.supports(MEB::OUT_ARG_DfDp,l).none()) {
      MEB::Derivative<Scalar>
        DfDp_bar_l = outArgs.get_DfDp(l);
      DfDp_bar[l] = DfDp_bar_l;
      DfDp_bar_mv[l] = rcp_dynamic_cast<ProductMultiVectorBase<Scalar> >(
        DfDp_bar_l.getMultiVector(), true );
      TEUCHOS_TEST_FOR_EXCEPTION(
        (
          !DfDp_bar_l.isEmpty()
          &&
          (
            is_null(DfDp_bar_mv[l])
            ||
            DfDp_bar_l.getMultiVectorOrientation() != MEB::DERIV_MV_BY_COL
            )
          ),
        std::logic_error,
        "Error, we currently can only handle DfDp as an column-based multi-vector!"
        );
    }
  }

  RCP<BlockedLinearOpBase<Scalar> > W_op_bar;
  if (outArgs.supports(MEB::OUT_ARG_W_op)) {
    W_op_bar = rcp_dynamic_cast<BlockedLinearOpBase<Scalar> >(
      outArgs.get_W_op(), true
      );
  }

  RCP<VectorBase<Scalar> >
    g_bar = outArgs.get_g(0);

  MEB::Derivative<Scalar> DgDx_dot_bar;
  RCP<ProductMultiVectorBase<Scalar> > DgDx_dot_bar_mv;
  if (!outArgs.supports(MEB::OUT_ARG_DgDx_dot,0).none()) {
    DgDx_dot_bar = outArgs.get_DgDx_dot(0);
    DgDx_dot_bar_mv = rcp_dynamic_cast<ProductMultiVectorBase<Scalar> >(
      DgDx_dot_bar.getMultiVector(), true );
    TEUCHOS_TEST_FOR_EXCEPTION(
      (
        !DgDx_dot_bar.isEmpty()
        &&
        (
          is_null(DgDx_dot_bar_mv)
          ||
          DgDx_dot_bar.getMultiVectorOrientation() != MEB::DERIV_TRANS_MV_BY_ROW
          )
        ),
      std::logic_error,
      "Error, we currently can only handle DgDx_dot as an row-based multi-vector!"
      );
  }

  MEB::Derivative<Scalar> DgDx_bar;
  RCP<ProductMultiVectorBase<Scalar> > DgDx_bar_mv;
  if (!outArgs.supports(MEB::OUT_ARG_DgDx,0).none()) {
    DgDx_bar = outArgs.get_DgDx(0);
    DgDx_bar_mv = rcp_dynamic_cast<ProductMultiVectorBase<Scalar> >(
      DgDx_bar.getMultiVector(), true );
    TEUCHOS_TEST_FOR_EXCEPTION(
      (
        !DgDx_bar.isEmpty()
        &&
        (
          is_null(DgDx_bar_mv)
          ||
          DgDx_bar.getMultiVectorOrientation() != MEB::DERIV_TRANS_MV_BY_ROW
          )
        ),
      std::logic_error,
      "Error, we currently can only handle DgDx as an row-based multi-vector!"
      );
  }

  Array<MEB::Derivative<Scalar> > DgDp_bar(Np);
  Array<RCP<MultiVectorBase<Scalar> > > DgDp_bar_mv(Np);
  for ( int l = 0; l < Np; ++l ) {
    if (!outArgs.supports(MEB::OUT_ARG_DgDp,0,l).none()) {
      MEB::Derivative<Scalar>
        DgDp_bar_l = outArgs.get_DgDp(0,l);
      DgDp_bar[l] = DgDp_bar_l;
      DgDp_bar_mv[l] = DgDp_bar_l.getMultiVector();
      TEUCHOS_TEST_FOR_EXCEPTION(
        !DgDp_bar_l.isEmpty() && is_null(DgDp_bar_mv[l]),
        std::logic_error,
        "Error, we currently can only handle DgDp as some type of multi-vector!"
        );
    }
  }

  //
  // C) Evaluate the model
  //

  // C.1) Set up storage and do some initializations

  MEB::InArgs<Scalar>
    periodInArgs = periodModel_->createInArgs();
  // ToDo: The above will have to change if you allow different structures for
  // each period model!
  
  // Set all of the parameters that will just be passed through 
  for ( int l = 0; l < Np; ++l )
    periodInArgs.set_p( period_l(l), inArgs.get_p(l) ); // Can be null just fine
  
  MEB::OutArgs<Scalar>
    periodOutArgs = periodModel_->createOutArgs();
  // ToDo: The above will have to change if you allow different structures for
  // each period model!
  
  // Create storage for period g's that will be summed into global g_bar
  periodOutArgs.set_g(
    g_index_, createMember<Scalar>( periodModel_->get_g_space(g_index_) ) );

  // Zero out global g_bar that will be summed into below
  if (!is_null(g_bar) )
    assign( g_bar.ptr(), ST::zero() );
  
  // Set up storage for peroid DgDp[l] objects that will be summed into global
  // DgDp_bar[l] and zero out DgDp_bar[l] that will be summed into.
  for ( int l = 0; l < Np; ++l ) {
    if ( !is_null(DgDp_bar_mv[l]) ) {
      assign(DgDp_bar_mv[l].ptr(), ST::zero());
      periodOutArgs.set_DgDp(
        g_index_, period_l(l),
        create_DgDp_mv(
          *periodModel_, g_index_, period_l(l),
          DgDp_bar[l].getMultiVectorOrientation()
          )
        );
    }
  }
  
  // C.2) Loop over periods and assemble the model

  for ( int i = 0; i < N; ++i ) {

    VOTSME thyraModel_outputTempState(periodModels_[i],out,verbLevel);

    // C.2.a) Set period-speicific InArgs and OutArgs

    for ( int k = 0; k < numPeriodZs(); ++k )
      periodInArgs.set_p( z_indexes_[k], z_[i][k] );

    if (!is_null(x_bar))
      periodInArgs.set_x(x_bar->getVectorBlock(i));

    if (!is_null(f_bar))
      periodOutArgs.set_f(f_bar->getNonconstVectorBlock(i)); // Updated in place!

    if ( !is_null(W_op_bar) )
      periodOutArgs.set_W_op(W_op_bar->getNonconstBlock(i,i));

    for ( int l = 0; l < Np; ++l ) {
      if ( !is_null(DfDp_bar_mv[l]) ) {
        periodOutArgs.set_DfDp(
          period_l(l),
          MEB::Derivative<Scalar>(
            DfDp_bar_mv[l]->getNonconstMultiVectorBlock(i),
            MEB::DERIV_MV_BY_COL
            )
          );
      }
    }
    
    if ( !is_null(DgDx_dot_bar_mv) ) {
      periodOutArgs.set_DgDx_dot(
        g_index_,
        MEB::Derivative<Scalar>(
          DgDx_dot_bar_mv->getNonconstMultiVectorBlock(i),
          MEB::DERIV_TRANS_MV_BY_ROW
          )
        );
    }
    
    if ( !is_null(DgDx_bar_mv) ) {
      periodOutArgs.set_DgDx(
        g_index_,
        MEB::Derivative<Scalar>(
          DgDx_bar_mv->getNonconstMultiVectorBlock(i),
          MEB::DERIV_TRANS_MV_BY_ROW
          )
        );
    }

    // C.2.b) Evaluate the period model

    periodModels_[i]->evalModel( periodInArgs, periodOutArgs );

    // C.2.c) Process output arguments that need processed

    // Sum into global g_bar
    if (!is_null(g_bar)) {
      Vp_StV( g_bar.ptr(), g_weights_[i], *periodOutArgs.get_g(g_index_) );
    }

    // Sum into global DgDp_bar
    for ( int l = 0; l < Np; ++l ) {
      if ( !is_null(DgDp_bar_mv[l]) ) {
        update(
          g_weights_[i],
          *periodOutArgs.get_DgDp(g_index_,period_l(l)).getMultiVector(),
          DgDp_bar_mv[l].ptr()
          );
      }
    }
    
    // Scale DgDx_dot_bar_mv[i]
    if ( !is_null(DgDx_dot_bar_mv) ) {
      scale( g_weights_[i],
        DgDx_dot_bar_mv->getNonconstMultiVectorBlock(i).ptr() );
    }
    
    // Scale DgDx_bar_mv[i]
    if ( !is_null(DgDx_bar_mv) ) {
      scale( g_weights_[i],
        DgDx_bar_mv->getNonconstMultiVectorBlock(i).ptr() );
    }

  }

  // ToDo: We need to do some type of global sum of g_bar and DgDp_bar to
  // account for other clusters of processes.  I might do this with a separate
  // non-ANA class.

  // Once we get here, all of the quantities should be updated and we should
  // be all done!

  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();

}


// private


template<class Scalar>
void DefaultMultiPeriodModelEvaluator<Scalar>::set_z_indexes_and_create_period_l_map(
  const Array<int> &z_indexes
  )
{
#ifdef TEUCHOS_DEBUG
  TEUCHOS_TEST_FOR_EXCEPT( Np_ <= 0 && "Error, Np must be set!" );
#endif
  z_indexes_ = z_indexes;
  period_l_map_.resize(0);
  const int numTotalParams = Np_ + z_indexes_.size();
  Array<int>::const_iterator
    z_indexes_itr = z_indexes_.begin(),
    z_indexes_end = z_indexes_.end();
  int last_z_index = -1;
  for ( int k = 0; k < numTotalParams; ++k ) {
    if ( z_indexes_itr == z_indexes_end || k < *z_indexes_itr ) {
      // This is a "free" parameter subvector
      period_l_map_.push_back(k);
    }
    else {
      // This is a "fixed" period parameter subvector so increment
      // z_indexes iterator.
#ifdef TEUCHOS_DEBUG
      TEUCHOS_TEST_FOR_EXCEPT( k != *z_indexes_itr && "This should never happen!" );
#endif
      const int tmp_last_z_index = *z_indexes_itr;
      ++z_indexes_itr;
      if ( z_indexes_itr != z_indexes_end ) {
#ifdef TEUCHOS_DEBUG
        if ( last_z_index >= 0 ) {
          TEUCHOS_TEST_FOR_EXCEPTION(
            *z_indexes_itr <= last_z_index, std::logic_error,
            "Error, the z_indexes array = " << toString(z_indexes_)
            << " is not sorted or contains duplicate entries!"
            );
        }
#endif
        last_z_index = tmp_last_z_index;
      }
    }
  }
}


template<class Scalar>
void DefaultMultiPeriodModelEvaluator<Scalar>::wrapNominalValuesAndBounds()
{

  using Teuchos::rcp_dynamic_cast;
  typedef ModelEvaluatorBase MEB;

  nominalValues_ = this->createInArgs();
  lowerBounds_ = this->createInArgs();
  upperBounds_ = this->createInArgs();

  const MEB::InArgs<Scalar>
    periodNominalValues = periodModel_->getNominalValues(),
    periodLowerBounds = periodModel_->getLowerBounds(),
    periodUpperBounds = periodModel_->getUpperBounds();
  
  if (periodNominalValues.supports(MEB::IN_ARG_x)) {

    if( !is_null(periodNominalValues.get_x()) ) {
      // If the first peroid model has nominal values for x, then all of them
      // must also!
      RCP<Thyra::ProductVectorBase<Scalar> >
        x_bar_init = createProductVector(x_bar_space_);
      const int N = this->N();
      for ( int i = 0; i < N; ++i ) {
        assign(
          x_bar_init->getNonconstVectorBlock(i).ptr(),
          *periodModels_[i]->getNominalValues().get_x()
          );
      }
      nominalValues_.set_x(x_bar_init);
    }
      
    if( !is_null(periodLowerBounds.get_x()) ) {
      // If the first peroid model has lower bounds for for x, then all of
      // them must also!
      RCP<Thyra::ProductVectorBase<Scalar> >
        x_bar_l = createProductVector(x_bar_space_);
      const int N = this->N();
      for ( int i = 0; i < N; ++i ) {
        assign(
          x_bar_l->getNonconstVectorBlock(i).ptr(),
          *periodModels_[i]->getLowerBounds().get_x()
          );
      }
      lowerBounds_.set_x(x_bar_l);
    }
      
    if( !is_null(periodUpperBounds.get_x()) ) {
      // If the first peroid model has upper bounds for for x, then all of
      // them must also!
      RCP<Thyra::ProductVectorBase<Scalar> >
        x_bar_u = createProductVector(x_bar_space_);
      const int N = this->N();
      for ( int i = 0; i < N; ++i ) {
        assign(
          x_bar_u->getNonconstVectorBlock(i).ptr(),
          *periodModels_[i]->getUpperBounds().get_x()
          );
      }
      upperBounds_.set_x(x_bar_u);
    }

  }

  // There can only be one set of nominal values for the non-period parameters
  // so just take them from the first period!
  for ( int l = 0; l < Np_; ++l ) {
    const int period_l = this->period_l(l);
    nominalValues_.set_p(l,periodNominalValues.get_p(period_l));
    lowerBounds_.set_p(l,periodLowerBounds.get_p(period_l));
    upperBounds_.set_p(l,periodUpperBounds.get_p(period_l));
  }
  
}


template<class Scalar>
RCP<ProductVectorBase<Scalar> >
DefaultMultiPeriodModelEvaluator<Scalar>::createProductVector(
  const RCP<const ProductVectorSpaceBase<Scalar> > &prodVecSpc
  )
{
  return Teuchos::rcp_dynamic_cast<ProductVectorBase<Scalar> >(
    Thyra::createMember<Scalar>(prodVecSpc), true );
}


} // namespace Thyra


#endif // THYRA_DEFAULT_MULTI_PERIOD_MODEL_EVALUATOR_HPP