text:
// Copyright (C) 2011 by the BEM++ Authors
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
// THE SOFTWARE.

#include "bempp/common/config_ahmed.hpp"

#include "../type_template.hpp"
#include "../check_arrays_are_close.hpp"

#include "assembly/context.hpp"
#include "assembly/discrete_boundary_operator.hpp"
#include "assembly/identity_operator.hpp"
#include "assembly/laplace_3d_adjoint_double_layer_boundary_operator.hpp"
#include "assembly/laplace_3d_double_layer_boundary_operator.hpp"
#include "assembly/laplace_3d_hypersingular_boundary_operator.hpp"
#include "assembly/laplace_3d_single_layer_boundary_operator.hpp"
#include "assembly/numerical_quadrature_strategy.hpp"
#include "grid/grid_factory.hpp"
#include "space/piecewise_constant_scalar_space.hpp"
#include "space/piecewise_linear_continuous_scalar_space.hpp"

#include <boost/test/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/type_traits/is_complex.hpp>
#include "grid/grid.hpp"

using namespace Bempp;

// Tests

namespace
{

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_V_plus_K
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
            context, domain, range, dualToRange) +
                laplace3dDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_V_plus_3_K
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                3. * laplace3dDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_3_V_plus_K
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return 3. * laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                laplace3dDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_3_V_plus_4_K
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return 3. * laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                4. * laplace3dDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_V_plus_K_plus_Kp
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                laplace3dDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                laplace3dAdjointDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_3_V_plus_K_plus_Kp
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return 3. * laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                laplace3dDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                laplace3dAdjointDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_V_plus_3_K_plus_Kp
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                3. * laplace3dDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                laplace3dAdjointDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_V
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_3_V
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return 3. * laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_I
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return identityOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_3_I
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return 3. * identityOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_I_plus_V
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return identityOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_V_plus_I
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                identityOperator<BasisFunctionType, ResultType>(
                                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_3_I_plus_V
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return 3. * identityOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_I_plus_3_V
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return identityOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                3. * laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_2_I_plus_V_plus_I_plus_4_K
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return 2. * identityOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                identityOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange) +
                2. * laplace3dDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange);
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_K_plus_adj_Kp
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return laplace3dDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange, "DLP") +
                adjoint(laplace3dAdjointDoubleLayerBoundaryOperator<
                        BasisFunctionType, ResultType>(
                            context, dualToRange, dualToRange /* hack */, domain, "ADLP"));
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_K_plus_I_plus_4_adj_Kp
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return laplace3dDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange, "DLP") +
                identityOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange, "Id") +
                4. * adjoint(laplace3dAdjointDoubleLayerBoundaryOperator<
                        BasisFunctionType, ResultType>(
                            context, dualToRange, dualToRange /* hack */, domain, "ADLP"));
    }
};

template <typename BasisFunctionType, typename ResultType>
struct CreateOperator_K_plus_I_plus_5_V_plus_4_adj_Kp
{
    BoundaryOperator<BasisFunctionType, ResultType> operator() (
            const shared_ptr<const Context<BasisFunctionType, ResultType> >& context,
            const shared_ptr<const Space<BasisFunctionType> >& domain,
            const shared_ptr<const Space<BasisFunctionType> >& range,
            const shared_ptr<const Space<BasisFunctionType> >& dualToRange) const
    {
        return laplace3dDoubleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange, "DLP") +
                identityOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange, "Id") +
                5. * laplace3dSingleLayerBoundaryOperator<BasisFunctionType, ResultType>(
                    context, domain, range, dualToRange, "SLP") +
                4. * adjoint(laplace3dAdjointDoubleLayerBoundaryOperator<
                        BasisFunctionType, ResultType>(
                            context, dualToRange, dualToRange /* hack */, domain, "ADLP"));
    }
};

template <typename BasisFunctionType, typename ResultType,
          template<typename BFT, typename RT> class CreateOperatorFunctor>
struct TestInDenseMode
{
    void test()
    {
        typedef ResultType RT;
        typedef typename ScalarTraits<ResultType>::RealType RealType;
        typedef BasisFunctionType BFT;

        CreateOperatorFunctor<BFT, RT> createOperator;

        GridParameters params;
        params.topology = GridParameters::TRIANGULAR;
        shared_ptr<Grid> grid = GridFactory::importGmshGrid(
                    params, "meshes/cube-12-reoriented.msh", false /* verbose */);

        shared_ptr<Space<BFT> > pwiseConstants(
                    new PiecewiseConstantScalarSpace<BFT>(grid));
        shared_ptr<Space<BFT> > pwiseLinears(
                    new PiecewiseLinearContinuousScalarSpace<BFT>(grid));

        AccuracyOptions accuracyOptions;
        accuracyOptions.doubleRegular.setRelativeQuadratureOrder(2);
        shared_ptr<NumericalQuadratureStrategy<BFT, RT> > quadStrategy(
                    new NumericalQuadratureStrategy<BFT, RT>(accuracyOptions));

        AssemblyOptions assemblyOptionsRef;
        assemblyOptionsRef.setVerbosityLevel(VerbosityLevel::LOW);
        shared_ptr<Context<BFT, RT> > contextRef(
                    new Context<BFT, RT>(quadStrategy, assemblyOptionsRef));

        BoundaryOperator<BFT, RT> opRef = createOperator(
                    contextRef, pwiseLinears, pwiseLinears, pwiseConstants);
        arma::Mat<RT> weakFormRef = opRef.weakForm()->asMatrix();

        AssemblyOptions assemblyOptionsTest;
        assemblyOptionsTest.setVerbosityLevel(VerbosityLevel::LOW);
        assemblyOptionsTest.enableJointAssembly();
        shared_ptr<Context<BFT, RT> > contextTest(
                    new Context<BFT, RT>(quadStrategy, assemblyOptionsTest));

        BoundaryOperator<BFT, RT> opTest = createOperator(
                    contextTest, pwiseLinears, pwiseLinears, pwiseConstants);
        arma::Mat<RT> weakFormTest = opTest.weakForm()->asMatrix();

        BOOST_CHECK(check_arrays_are_close<RT>(
                        weakFormTest, weakFormRef,
                        10. * std::numeric_limits<RealType>::epsilon()));
    }
};

template <typename BasisFunctionType, typename ResultType,
          template<typename BFT, typename RT> class CreateOperatorFunctor>
struct TestInAcaMode
{
    void test()
    {
        typedef ResultType RT;
        typedef typename ScalarTraits<ResultType>::RealType RealType;
        typedef BasisFunctionType BFT;

        CreateOperatorFunctor<BFT, RT> createOperator;

        GridParameters params;
        params.topology = GridParameters::TRIANGULAR;
        shared_ptr<Grid> grid = GridFactory::importGmshGrid(
                    params, "meshes/cube-12-reoriented.msh", false /* verbose */);

        shared_ptr<Space<BFT> > pwiseConstants(
                    new PiecewiseConstantScalarSpace<BFT>(grid));
        shared_ptr<Space<BFT> > pwiseLinears(
                    new PiecewiseLinearContinuousScalarSpace<BFT>(grid));

        AccuracyOptions accuracyOptions;
        accuracyOptions.doubleRegular.setRelativeQuadratureOrder(2);
        shared_ptr<NumericalQuadratureStrategy<BFT, RT> > quadStrategy(
                    new NumericalQuadratureStrategy<BFT, RT>(accuracyOptions));

        AcaOptions acaOptions;
        acaOptions.minimumBlockSize = 2;

        AssemblyOptions assemblyOptionsRef;
        assemblyOptionsRef.setVerbosityLevel(VerbosityLevel::LOW);
        assemblyOptionsRef.switchToAcaMode(acaOptions);
        assemblyOptionsRef.setMaxThreadCount(1);
        shared_ptr<Context<BFT, RT> > contextRef(
                    new Context<BFT, RT>(quadStrategy, assemblyOptionsRef));

        std::cout << "Assembling ref op.\n";
        BoundaryOperator<BFT, RT> opRef = createOperator(
                    contextRef, pwiseLinears, pwiseLinears, pwiseConstants);
        arma::Mat<RT> weakFormRef = opRef.weakForm()->asMatrix();

        AssemblyOptions assemblyOptionsTest;
        //assemblyOptionsTest.setVerbosityLevel(VerbosityLevel::LOW);
        assemblyOptionsTest.enableJointAssembly();
        assemblyOptionsTest.switchToAcaMode(acaOptions);
        assemblyOptionsTest.setMaxThreadCount(1);
        shared_ptr<Context<BFT, RT> > contextTest(
                    new Context<BFT, RT>(quadStrategy, assemblyOptionsTest));

        std::cout << "Assembling test op.\n";
        BoundaryOperator<BFT, RT> opTest = createOperator(
                    contextTest, pwiseLinears, pwiseLinears, pwiseConstants);
        arma::Mat<RT> weakFormTest = opTest.weakForm()->asMatrix();

        BOOST_CHECK(check_arrays_are_close<RT>(
                        weakFormTest, weakFormRef, 2. * acaOptions.eps));
    }
};

} // namespace

BOOST_AUTO_TEST_SUITE(JointAssembly)

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_V_plus_K,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_V_plus_K> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_3_V_plus_K,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_3_V_plus_K> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_V_plus_3_K,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_V_plus_3_K> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_V,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_V> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_3_V,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_3_V> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_V_plus_K_plus_Kp,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_V_plus_K_plus_Kp> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_3_V_plus_K_plus_Kp,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_3_V_plus_K_plus_Kp> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_V_plus_3_K_plus_Kp,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_V_plus_3_K_plus_Kp> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_I,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_I> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_3_I,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_3_I> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_V_plus_I,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_V_plus_I> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_3_I_plus_V,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_3_I_plus_V> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_I_plus_V,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_I_plus_V> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_I_plus_3_V,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_I_plus_3_V> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_2_I_plus_V_plus_I_plus_4_K,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_2_I_plus_V_plus_I_plus_4_K> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_K_plus_I_plus_4_adj_Kp,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_K_plus_I_plus_4_adj_Kp> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_K_plus_I_plus_5_V_plus_4_adj_Kp,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_K_plus_I_plus_5_V_plus_4_adj_Kp> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_dense_mode_for_K_plus_adj_Kp,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInDenseMode<BFT, RT, CreateOperator_K_plus_adj_Kp> test;
    test.test();
}

#ifdef WITH_AHMED

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_V_plus_K,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_V_plus_K> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_3_V_plus_K,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_3_V_plus_K> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_V_plus_3_K,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_V_plus_3_K> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_V,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_V> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_3_V,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_3_V> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_V_plus_K_plus_Kp,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_V_plus_K_plus_Kp> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_3_V_plus_K_plus_Kp,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_3_V_plus_K_plus_Kp> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_V_plus_3_K_plus_Kp,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_V_plus_3_K_plus_Kp> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_I,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_I> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_3_I,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_3_I> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_V_plus_I,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_V_plus_I> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_3_I_plus_V,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_3_I_plus_V> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_I_plus_V,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_I_plus_V> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_I_plus_3_V,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_I_plus_3_V> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_2_I_plus_V_plus_I_plus_4_K,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_2_I_plus_V_plus_I_plus_4_K> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_K_plus_adj_Kp,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_K_plus_adj_Kp> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_K_plus_I_plus_4_adj_Kp,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_K_plus_I_plus_4_adj_Kp> test;
    test.test();
}

BOOST_AUTO_TEST_CASE_TEMPLATE(joint_assembly_works_in_aca_mode_for_K_plus_I_plus_5_V_plus_4_adj_Kp,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    TestInAcaMode<BFT, RT, CreateOperator_K_plus_I_plus_5_V_plus_4_adj_Kp> test;
    test.test();
}

#endif // WITH_AHMED

BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: nicolas-chaulet/bempp

max_issues_repo_name: nicolas-chaulet/bempp

max_forks_repo_name: nicolas-chaulet/bempp


text:
// Copyright 2021 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. See LICENSE.txt for details.
//
// Author: Denis Bueno

#include "bmc_one_shot_refiner.h"

#include <boost/range/algorithm/transform.hpp>

#include "fp_mus_refiner.h"
#include "supp/reachability_graph.h"
#include "supp/expr_substitution.h"

using namespace std;

namespace {

class CxState {
 public:
  CxState(int64_t i) : index_(i) {}

  void foo() { (void)index_; }

 private:
  int64_t index_;
};

}

//^----------------------------------------------------------------------------^

namespace euforia {

std::ostream& operator<<(std::ostream& os, const BmcStep& s) {
  fmt::print(os, "BmcStep constraints: {}\n", s.constraints());
  return os;
}

class BmcOneShotRefiner::Impl {
 public:
  Impl(BmcOneShotRefiner& parent,
       const ReachabilityGraph& reached)
      : parent_(parent), reached_(reached), vars_(ctx()),
        steps_(reached.cx_length(), BmcStep(ctx(), 0)) {
  }

  //! Returns true if feasible
  bool OneShotQuery();

  z3::expr_vector FindInterpolants(z3::expr_vector core);
  unique_ptr<Counterexample> mk_counterexample();

  z3::context& ctx() const { return parent_.ctx(); }
  const auto& axsys() const { return parent_.axsys(); }
  const auto& cxsys() const { return parent_.cxsys(); }
  Solver& solver() const { return *parent_.solver_; }

  bool lemma_refutes_abstract_cx(const z3::expr&);

 private:
  BmcOneShotRefiner& parent_;
  const ReachabilityGraph& reached_;
  z3::expr_vector vars_;
  // steps of the counterexample
  std::vector<BmcStep> steps_;
  std::unique_ptr<Solver> absolver_;

  ExprSubstitution MakeRenamer(int64_t i);
  ExprSubstitution mk_renamer(int64_t i) const;
  ExprSubstitution mk_input_renamer(int64_t i) const;
  ExprSubstitution mk_next_renamer(int64_t i) const;
  z3::expr get_step_var(int64_t i, const z3::expr& e) const;
};

//^----------------------------------------------------------------------------^

BmcOneShotRefiner::BmcOneShotRefiner(
    AbstractChecker::Impl& achk,
    const std::shared_ptr<Solver>& s, const ReachabilityGraph& acx,
    int64_t nr)
    : Refiner(achk, s, nr),
      pimpl_(std::make_unique<BmcOneShotRefiner::Impl>(*this, acx)) {}

BmcOneShotRefiner::~BmcOneShotRefiner() = default;

unique_ptr<Counterexample>
BmcOneShotRefiner::make_counterexample() {
  return pimpl_->mk_counterexample();
}

//^----------------------------------------------------------------------------^

boost::optional<z3::expr_vector>
BmcOneShotRefiner::operator()() {
  bool feasible = pimpl_->OneShotQuery();
  if (!feasible) {
    auto core = solver().unsat_core();
    return pimpl_->FindInterpolants(core);
  } else {
    return boost::none;
  }
}

//^----------------------------------------------------------------------------^

bool BmcOneShotRefiner::Impl::OneShotQuery() {
  vars_ = z3::expr_vector(ctx()); // clear
  // Creates a formula for one BMC query from the reachability graph.
  int64_t i = 1;
  for (auto it = reached_.cx_begin(), ie = reached_.cx_end(); it != ie;
       ++it, ++i) {
    const auto& reach_step = *it;
    // this loop will fill bmc_step in.
    auto& bmc_step = steps_[i-1] = {ctx(), it->state_index()};
    // Renames current state variables v to v-(i-1). Renames input variables
    // v to v-(i-1).
    const auto& subst_current = bmc_step.subst_ =
        mk_renamer(i-1);
    boost::copy(subst_current.dst(),
                ExprVectorInserter(bmc_step.current_vars_));
    const auto rename_current = [&](const z3::expr& abs_e) {
      return subst_current(axsys().Concretize(abs_e));
    };

    if (i == 1) {
      boost::transform(ExprConjuncts(cxsys().init_state()),
                       bmc_step.constraint_inserter(),
                       rename_current);
    }

    boost::transform(
        reach_step.state(), bmc_step.constraint_inserter(),
        rename_current);

    if (reach_step.has_transition()) {
      // Renames inputs v to v-i
      bmc_step.input_vars_ = z3::expr_vector(ctx());
      ExprSubstitution subst_input = mk_input_renamer(i);
      boost::copy(subst_input.dst(),
                  ExprVectorInserter(*bmc_step.input_vars_));
      // Renames next-state vars v to v-i.
      ExprSubstitution subst_next = mk_next_renamer(i);
      bmc_step.next_vars_ = z3::expr_vector(ctx());
      boost::copy(subst_next.dst(),
                  ExprVectorInserter(*bmc_step.next_vars_));
      bmc_step.subst_ = subst_current & subst_next & subst_input;
      auto rename_both = [&](const z3::expr& abs_e) {
        return bmc_step.subst_(axsys().Concretize(abs_e));
      };
      if (euforia_config.use_osbmc_slicing) {
        auto preds = rename_both(expr_mk_and(
            ctx(), parent_.sliced_trans_predicates(reach_step.step_model())));
        boost::copy(ExprConjuncts(preds), bmc_step.constraint_inserter());
      } else {
        auto new_trans = bmc_step.subst_(cxsys().trans());
        boost::copy(
            ExprConjuncts(new_trans), bmc_step.constraint_inserter());
      }

      auto inputs = rename_both(expr_mk_and(ctx(), reach_step.input()));
      boost::copy(
          ExprConjuncts(inputs), bmc_step.constraint_inserter());
    }
    logger.Log(4, "bmcstep {} constraints size: {}",
               bmc_step.state_index(), bmc_step.constraints().size());
  }

  vector<z3::expr> bmc;
  for (const auto& bmc_step : steps_) {
    logger.Log(5, "bmcstep {} constraints:\n{{{{{{{}}}}}}}",
               bmc_step.state_index(), bmc_step);
    boost::copy(bmc_step.constraints(), back_inserter(bmc));
  }
  if (!LOG_DIR.empty()) {
    auto benchmark_name = fmt::format(
        "{}/ref{}-bmc.{}.smt2", LOG_DIR, parent_.num_refinements_+1,
        parent_.solver_->version());
    ofstream outfile(benchmark_name);
    // fmt::print(outfile, "; {} input constraints\n", i_hat.size());
    // fmt::print(outfile, "; {} target constraints\n", s_next_hat.size());
    fmt::print(outfile, "; {} assumptions\n", bmc.size());
    // fmt::print(outfile, "; {} trans constraints\n", num_trans_predicates);
    solver().DumpBenchmark(outfile, bmc.size(), bmc.data());
    logger.Log(3, "logging benchmark {}", benchmark_name);
  }
  logger.Log(5, "bmc one-shot query:\n{{{{{{{}}}}}}}", bmc);
  TimerDuration d(0);
  ScopedTimeKeeper bmc_time(&d);
  auto feasible = solver().Check(bmc) == CheckResult::kSat;
  auto check_duration = bmc_time.get();
  if (!feasible) {
    logger.Log(3, "one-shot BMC is CUNSAT [{:.3f} s]",
               check_duration.count());
  } else {
    logger.Log(3, "one-shot BMC is CSAT [{:.3f} s]",
               check_duration.count());
  }
  return feasible;
}

z3::expr_vector
BmcOneShotRefiner::Impl::FindInterpolants(z3::expr_vector core) {
    FpMusRefiner get_interpolants(parent_.achk_, parent_.solver_,
                                  steps_, core,
                                  parent_.num_refinements_);
    auto cores = get_interpolants();
    // Checks that assuming the negations of all the cores refutes the abstract
    // BMC query.
    EDEBUG_CODE(
        z3::expr_vector lemmas(ctx());
        boost::transform(cores, ExprVectorInserter(lemmas), expr_not);
        assert(lemma_refutes_abstract_cx(expr_mk_and(lemmas))));
    return cores;
}

unique_ptr<Counterexample> BmcOneShotRefiner::Impl::mk_counterexample() {
  z3::expr_vector vars(ctx());
  for (const auto& bmc_step : steps_) {
    boost::copy(bmc_step.current_vars(), ExprVectorInserter(vars));
  }
  auto cx = std::make_unique<BmcOneShotCounterexample>();
  cx->set_model(solver().get_model(vars));
  // logger.Log(0, "one-shot bmc model:\n{}", *model);
  return cx;
}

ExprSubstitution BmcOneShotRefiner::Impl::mk_renamer(int64_t i) const {
  ExprSubstitution subst(ctx());
  for (const StateVar& var : cxsys().vars()) {
    auto passified_var = get_step_var(i, var.current());
    subst.AddSubstitution(var.current(), passified_var);
  }
  return subst;
}

ExprSubstitution BmcOneShotRefiner::Impl::mk_input_renamer(int64_t i) const {
  ExprSubstitution subst(ctx());
  for (const PrimaryInput& input : cxsys().inputs()) {
    auto passified_var = get_step_var(i, input);
    subst.AddSubstitution(input, passified_var);
  }
  return subst;
}

//! Renames next-state inputs with i suffix
ExprSubstitution BmcOneShotRefiner::Impl::mk_next_renamer(int64_t i) const {
  ExprSubstitution subst(ctx());
  for (const StateVar& var : cxsys().vars()) {
    auto passified_var = get_step_var(i, var.current());
    subst.AddSubstitution(var.next(), passified_var);
  }
  return subst;
}

z3::expr BmcOneShotRefiner::Impl::get_step_var(
    int64_t i, const z3::expr& e) const {
  assert(e.num_args() == 0);
  auto new_name = fmt::format("{}-{}", e.decl().name().str(), i);
  return e.ctx().constant(new_name.c_str(), e.get_sort());
}

bool
BmcOneShotRefiner::Impl::lemma_refutes_abstract_cx(const z3::expr& lemma) {
  // The lemma uses (X, Y, X'). Instantiate it for every step of the cx and
  // assert those instantiations.  Assert the bmc formula. Check for unsat.

  if (!absolver_) {
    // Initializes the solver with the abstract BMC query
    auto s = make_unique<Z3Solver>(ctx(), "QF_UF");
    z3::params p(ctx());
    p.set("model", false);
    p.set("unsat_core", false);
    s->set_params(p);
    absolver_ = move(s);
    for (auto it = steps_.begin(), ie = steps_.end(); it != ie; ++it) {
      absolver_->Add(axsys().Abstract(
              expr_mk_and(ctx(), it->constraints())));
    }
    // Sanity checks that the BMC query is abstractly feasible
    // assert(absolver_->Check() == CheckResult::kSat);
  }

  absolver_->Push();
  logger.Log(4, "lemma_refutes_abstract_cx: lemma:");
  logger.LogFold(4, "{}", axsys().Abstract(lemma));
  for (auto it = steps_.begin(), ie = steps_.end(); it != ie; ++it) {
    auto abs_lemma_bmc = axsys().Abstract(it->subst()(lemma));
    absolver_->Add(abs_lemma_bmc);
  }
  logger.Log(5, "lemma_refutes_abstract_cx: assertions:");
  logger.LogFold(5, "{}", absolver_->assertions());
  logger.Log(3, "checking that fp lemma refutes abstract cx");
  // Checks that the abstracted lemma refutes the BMC query
  auto result = absolver_->Check();
  absolver_->Pop();
  return result == CheckResult::kUnsat;
}

//^----------------------------------------------------------------------------^
} // end namespace euforia

ext: cc

max_stars_repo_name: dbueno/euforia

max_issues_repo_name: dbueno/euforia

max_forks_repo_name: dbueno/euforia