hawp_commons.hpp

Go to the documentation of this file.

#pragma once

#include <stdexcept>
#include <memory>

#include "hawp_paramset.hpp"
#include "hawp_evaluator.hpp"
#include "shapes/shape_extension_cache.hpp"


namespace waveblocks {
    namespace wavepackets {
        template<dim_t D, class MultiIndex>
        class AbstractScalarHaWpBasis
        {
        public:
            virtual real_t eps() const = 0;

            virtual HaWpParamSet<D> const& parameters() const = 0;

            virtual shapes::ShapeEnumSharedPtr<D, MultiIndex> shape() const = 0;

            template<int N> HaWpEvaluator<D,MultiIndex,N>
            create_evaluator(CMatrix<D,N> const& grid) const
            {
                return {eps(), &parameters(), shape().get(), grid};
            }

            template<int N> HaWpBasisVector<N>
            evaluate_basis(CMatrix<D,N> const& grid) const
            {
                return create_evaluator(grid).all();
            }

            template<int N> HaWpBasisVector<N>
            evaluate_basis(RMatrix<D,N> const& rgrid) const
            {
                CMatrix<D,N> cgrid = rgrid.template cast <complex_t>();
                return evaluate_basis(cgrid);
            }

            //     virtual HaWpBasisVector<Eigen::Dynamic> evaluate_basis(ComplexGrid<D,Eigen::Dynamic> const& grid) const
            //     {
            //         return create_evaluator<Eigen::Dynamic>(grid).all();
            //     }

            shapes::ShapeEnumSharedPtr<D,MultiIndex>
            extended_shape() const
            {
                return shape_extension_cache_.get_extended_shape( this->shape() );
            }

        private:
            shapes::ShapeExtensionCache<D,MultiIndex> shape_extension_cache_;
        };


        template<dim_t D, class MultiIndex>
        class AbstractScalarHaWp : public AbstractScalarHaWpBasis<D,MultiIndex>
        {
        public:
            virtual real_t eps() const = 0;
            virtual HaWpParamSet<D> const& parameters() const = 0;
            virtual shapes::ShapeEnumSharedPtr<D, MultiIndex> shape() const = 0;

            virtual Coefficients const& coefficients() const = 0;

            template<int N> CArray<1,N>
            evaluate(CMatrix<D,N> const& grid) const
            {
                if (this->shape()->n_entries() != (std::size_t)coefficients().size())
                    throw std::runtime_error("shape.size() != coefficients.size()");

                return this->template create_evaluator<N>(grid).reduce(coefficients());
            }

            template<int N> CArray<1,N>
            evaluate(RMatrix<D,N> const& rgrid) const
            {
                CMatrix<D,N> cgrid = rgrid.template cast <complex_t>();
                return evaluate(cgrid);
            }

            //     virtual HaWpBasisVector<Eigen::Dynamic> evaluate(ComplexGrid<D,Eigen::Dynamic> const& grid) const
            //     {
            //         return create_evaluator<Eigen::Dynamic>().reduce(coefficients());
            //     }

            complex_t prefactor() const
            {
                return real_t(1) / this->parameters().sdQ();
            }

            complex_t phasefactor() const
            {
                return std::exp(complex_t(0,1) * this->parameters().S() / eps() / eps());
            }
        };


        template<dim_t D, class MultiIndex>
        class ScalarHaWp : public AbstractScalarHaWp<D, MultiIndex>
        {
        public:
            real_t & eps()
            {
                return eps_;
            }

            real_t eps() const override
            {
                return eps_;
            }

            HaWpParamSet<D> & parameters()
            {
                return parameters_;
            }

            HaWpParamSet<D> const& parameters() const override
            {
                return parameters_;
            }

            shapes::ShapeEnumSharedPtr<D, MultiIndex> & shape()
            {
                return shape_;
            }

            shapes::ShapeEnumSharedPtr<D, MultiIndex> shape() const override
            {
                return shape_;
            }

            Coefficients & coefficients()
            {
                return coefficients_;
            }

            Coefficients const& coefficients() const override
            {
                return coefficients_;
            }

        private:
            real_t eps_;
            HaWpParamSet<D> parameters_;
            shapes::ShapeEnumSharedPtr<D, MultiIndex> shape_;
            Coefficients coefficients_;
        };


        template<dim_t D, class MultiIndex>
        class HomogeneousHaWp
        {
        public:
            class Component : public AbstractScalarHaWp<D,MultiIndex>
            {
            public:
                Component(HomogeneousHaWp const* const owner)
                    : owner_(owner)
                    , coefficients_()
                { }

                Component(Component&& that)
                    : owner_(that.owner_)
                    , shape_(that.shape_)
                    , coefficients_(std::move(that.coefficients_))
                { }

                Component(Component const& that)
                    : owner_(that.owner_)
                    , shape_(that.shape_)
                    , coefficients_(that.coefficients_)
                { }

                Component & operator=(Component&& that)
                {
                    shape_ = that.shape_;
                    coefficients_ = std::move(that.coefficients_);
                    return *this;
                }

                Component & operator=(Component const& that)
                {
                    shape_ = that.shape_;
                    coefficients_ = that.coefficients_;
                    return *this;
                }

                real_t eps() const override
                {
                    return owner_->eps();
                }

                HaWpParamSet<D> const& parameters() const override
                {
                    return owner_->parameters();
                }

                shapes::ShapeEnumSharedPtr<D, MultiIndex> & shape()
                {
                    return shape_;
                }

                shapes::ShapeEnumSharedPtr<D, MultiIndex> shape() const override
                {
                    return shape_;
                }

                Coefficients & coefficients()
                {
                    return coefficients_;
                }

                Coefficients const& coefficients() const override
                {
                    return coefficients_;
                }

            private:
                HomogeneousHaWp const* const owner_;

                shapes::ShapeEnumSharedPtr<D, MultiIndex> shape_;
                Coefficients coefficients_;
            };

            HomogeneousHaWp(std::size_t n)
                : eps_(0.0)
                , parameters_()
                , components_(n, Component(this))
            { }

            real_t & eps()
            {
                return eps_;
            }

            real_t eps() const
            {
                return eps_;
            }

            HaWpParamSet<D> & parameters()
            {
                return parameters_;
            }

            HaWpParamSet<D> const& parameters() const
            {
                return parameters_;
            }

            std::vector<Component> & components()
            {
                return components_;
            }

            std::vector<Component> const& components() const
            {
                return components_;
            }

            Component & component(std::size_t n)
            {
                return components_[n];
            }

            Component const& component(std::size_t n) const
            {
                return components_[n];
            }

            Component & operator[](std::size_t n)
            {
                return component(n);
            }

            Component const& operator[](std::size_t n) const
            {
                return component(n);
            }

            std::size_t n_components() const
            {
                return components_.size();
            }

            shapes::ShapeEnumSharedPtr<D,MultiIndex> union_shape() const
            {
                // check cache status
                bool rebuild_cache = false;

                if (union_cache_snapshot_.size() != n_components()) {
                    rebuild_cache = true;
                    union_cache_snapshot_.resize(n_components());
                }

                for (std::size_t n = 0; n < n_components() && !rebuild_cache; n++) {
                    if (union_cache_snapshot_[n] != component(n).shape().get())
                        rebuild_cache = true;
                }

                // rebuild cache if necessary
                if (rebuild_cache) {
                    std::vector< ShapeEnum<D,MultiIndex> const* > list(n_components());
                    for (std::size_t c = 0; c < n_components(); c++) {
                        list[c] = union_cache_snapshot_[c] = component(c).shape().get();
                    }
                    cached_shape_union_ = std::make_shared< ShapeEnum<D,MultiIndex> >(shapes::shape_enum::strict_union(list));
                }

                return cached_shape_union_;
            }

            template<int N>
            CArray<Eigen::Dynamic,N> evaluate(CMatrix<D,N> const& grid) const
            {
                ScalarHaWp<D,MultiIndex> unionwp;

                unionwp.eps() = eps();
                unionwp.parameters() = parameters();
                unionwp.shape() = union_shape();

                std::vector< ShapeEnum<D,MultiIndex>* > shape_list(n_components());
                std::vector< complex_t const* > coeffs_list(n_components());

                for (std::size_t n = 0; n < n_components(); n++) {
                    shape_list[n] = component(n).shape().get();
                    coeffs_list[n] = component(n).coefficients().data();
                }

                return unionwp.template create_evaluator<N>(grid).vector_reduce(shape_list.data(), coeffs_list.data(), n_components());
            }

            template<int N>
            CArray<Eigen::Dynamic,N> evaluate(RMatrix<D,N> const& rgrid) const
            {
                CMatrix<D,N> cgrid = rgrid.template cast<complex_t>();
                return evaluate(cgrid);
            }

        private:
            real_t eps_;
            HaWpParamSet<D> parameters_;
            std::vector<Component> components_;

            mutable std::vector< ShapeEnum<D,MultiIndex>* > union_cache_snapshot_;
            mutable shapes::ShapeEnumSharedPtr<D,MultiIndex> cached_shape_union_;
        }; // class HomogeneousHaWp


        template<dim_t D, class MultiIndex>
        class InhomogeneousHaWp
        {
        public:
            class Component : public AbstractScalarHaWp<D,MultiIndex>
            {
            public:
                Component(InhomogeneousHaWp const* const owner)
                    : owner_(owner)
                    , shape_()
                    , coefficients_()
                { }

                Component(Component&& that)
                    : owner_(that.owner_)
                    , parameters_(std::move(that.parameters_))
                    , shape_(that.shape_)
                    , coefficients_(std::move(that.coefficients_))
                { }

                Component(Component const& that)
                    : owner_(that.owner_)
                    , parameters_(that.parameters_)
                    , shape_(that.shape_)
                    , coefficients_(that.coefficients_)
                { }

                Component & operator=(Component&& that)
                {
                    parameters_ = std::move(that.parameters_);
                    shape_ = that.shape_;
                    coefficients_ = std::move(that.coefficients_);
                    return *this;
                }

                Component & operator=(Component const& that)
                {
                    parameters_ = that.parameters_;
                    shape_ = that.shape_;
                    coefficients_ = that.coefficients_;
                    return *this;
                }

                real_t eps() const override
                {
                    return owner_->eps();
                }

                HaWpParamSet<D> & parameters()
                {
                    return parameters_;
                }


                HaWpParamSet<D> const& parameters() const override
                {
                    return parameters_;
                }

                shapes::ShapeEnumSharedPtr<D, MultiIndex> & shape()
                {
                    return shape_;
                }

                shapes::ShapeEnumSharedPtr<D, MultiIndex> shape() const override
                {
                    return shape_;
                }

                Coefficients & coefficients()
                {
                    return coefficients_;
                }

                Coefficients const& coefficients() const override
                {
                    return coefficients_;
                }

            private:
                InhomogeneousHaWp const* const owner_;

                HaWpParamSet<D> parameters_;
                shapes::ShapeEnumSharedPtr<D, MultiIndex> shape_;
                Coefficients coefficients_;
            };

            InhomogeneousHaWp(std::size_t n)
                : eps_()
                , components_(n, Component(this))
            { }

            real_t & eps()
            {
                return eps_;
            }

            real_t eps() const
            {
                return eps_;
            }

            std::vector<Component> & components()
            {
                return components_;
            }

            std::vector<Component> const& components() const
            {
                return components_;
            }

            Component & component(std::size_t n)
            {
                return components_[n];
            }

            Component const& component(std::size_t n) const
            {
                return components_[n];
            }

            Component & operator[](std::size_t n)
            {
                return component(n);
            }

            Component const& operator[](std::size_t n) const
            {
                return component(n);
            }

            std::size_t n_components() const
            {
                return components_.size();
            }

            template<int N>
            CArray<Eigen::Dynamic,N> evaluate(CMatrix<D,N> const& grid) const
            {
                CArray<Eigen::Dynamic,N> result(n_components(),grid.cols());

                for (std::size_t c = 0; c < n_components(); c++) {
                    result.row(c) = component(c).evaluate(grid);
                }

                return result;
            }

            template<int N>
            CArray<Eigen::Dynamic,N> evaluate(RMatrix<D,N> const& rgrid) const
            {
                CMatrix<D,N> cgrid = rgrid.template cast<complex_t>();
                return evaluate(cgrid);
            }

        private:
            real_t eps_;
            std::vector<Component> components_;
        };
    }
}