HeartConfig.hpp

/*

Copyright (c) 2005-2018, University of Oxford.
All rights reserved.

University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.

This file is part of Chaste.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
 * Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.
 * Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.
 * Neither the name of the University of Oxford nor the names of its
   contributors may be used to endorse or promote products derived from this
   software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

*/


#ifndef HEARTCONFIG_HPP_
#define HEARTCONFIG_HPP_

#include <string>
#include <vector>
#include <set>
#include <map>
#include <boost/shared_ptr.hpp>

#include "UblasVectorInclude.hpp"

#include "ChasteParameters_2017_1.hpp"

#include "AbstractStimulusFunction.hpp"
#include "AbstractChasteRegion.hpp"
#include "ChastePoint.hpp"
#include "ChasteCuboid.hpp"
#include "ChasteEllipsoid.hpp"
#include "DistributedTetrahedralMeshPartitionType.hpp"
#include "PetscTools.hpp"
#include "FileFinder.hpp"

#include "ChasteSerialization.hpp"
#include "ChasteSerializationVersion.hpp"
#include <boost/serialization/split_member.hpp>
#include <boost/serialization/map.hpp>
#include <boost/serialization/set.hpp>

namespace cp = chaste::parameters::v2017_1;

// Forward declaration to avoid circular includes
class HeartFileFinder;


class HeartConfig
{
private:
    void CheckTimeSteps() const;

    friend class boost::serialization::access;
    template<class Archive>
    void save(Archive & archive, const unsigned int version) const
    {
        //Only the Master should be writing the configuration file
        if (PetscTools::AmMaster())
        {
            mpInstance->Write( true );
        }

        // Archive other member variables that don't appear in the XML
        if (version > 1)
        {
            archive & mEpiFraction;
            archive & mEndoFraction;
            archive & mMidFraction;
            archive & mIndexMid;
            archive & mIndexEpi;
            archive & mIndexEndo;
            archive & mUserAskedForCellularTransmuralHeterogeneities;
            archive & mUseMassLumping;
            archive & mUseMassLumpingForPrecond;
            archive & mUseReactionDiffusionOperatorSplitting;
            archive & mBathConductivities;
            archive & mTissueIdentifiers;
            archive & mBathIdentifiers;
            archive & mUseFixedNumberIterations;
            archive & mEvaluateNumItsEveryNSolves;
        }

        PetscTools::Barrier("HeartConfig::save");
    }

    template<class Archive>
    void load(Archive & archive, const unsigned int version)
    {
        LoadFromCheckpoint();

        // Load other member variables
        if (version > 1)
        {
            archive & mEpiFraction;
            archive & mEndoFraction;
            archive & mMidFraction;
            archive & mIndexMid;
            archive & mIndexEpi;
            archive & mIndexEndo;
            archive & mUserAskedForCellularTransmuralHeterogeneities;
            archive & mUseMassLumping;
            archive & mUseMassLumpingForPrecond;
            archive & mUseReactionDiffusionOperatorSplitting;
            archive & mBathConductivities;
            archive & mTissueIdentifiers;
            archive & mBathIdentifiers;
            archive & mUseFixedNumberIterations;
            archive & mEvaluateNumItsEveryNSolves;
        }
    }
    BOOST_SERIALIZATION_SPLIT_MEMBER()

    
    void LoadFromCheckpoint();

    void UpdateParametersFromResumeSimulation(boost::shared_ptr<cp::chaste_parameters_type> pResumeParameters);

public:
    typedef std::map<std::string, std::string> SchemaLocationsMap;

private:
    SchemaLocationsMap mSchemaLocations;

    void SetDefaultSchemaLocations();

public:
    static HeartConfig* Instance();

    void SetUseFixedSchemaLocation(bool useFixedSchemaLocation);

    void SetFixedSchemaLocations(const SchemaLocationsMap& rSchemaLocations);

    void SetParametersFile(const std::string& rFileName);

    void Write(bool useArchiveLocationInfo=false, std::string subfolderName="output");

    void CopySchema(const std::string& rToDirectory);

    boost::shared_ptr<cp::chaste_parameters_type> ReadFile(const std::string& rFileName);

    static void Reset();

    ~HeartConfig(); 
    unsigned GetVersionFromNamespace(const std::string& rNamespaceUri);

    //
    //  Get methods
    //

    FileFinder GetParametersFilePath();

    // Methods for asking the configuration file about which sections are defined.

    bool IsSimulationDefined() const;

    bool IsSimulationResumed() const;

    // Simulation
    unsigned GetSpaceDimension() const; 
    double GetSimulationDuration() const; 
    cp::domain_type GetDomain() const;

    cp::ionic_model_selection_type GetDefaultIonicModel() const;

    template<unsigned DIM>
    void GetIonicModelRegions(std::vector<boost::shared_ptr<AbstractChasteRegion<DIM> > >& rDefinedRegions,
                              std::vector<cp::ionic_model_selection_type>& rIonicModels) const;

    void SetIonicModelRegions(std::vector<ChasteCuboid<3> >& rDefinedRegions,
                              std::vector<cp::ionic_model_selection_type>& rIonicModels) const;

    bool IsMeshProvided() const; 
    bool GetCreateMesh() const; 
    bool GetCreateSlab() const; 
    bool GetCreateSheet() const; 
    bool GetCreateFibre() const; 
    bool GetLoadMesh() const; 
    void GetSlabDimensions(c_vector<double, 3>& slabDimensions) const;
    void GetSheetDimensions(c_vector<double, 2>& sheetDimensions) const;
    void GetFibreLength(c_vector<double, 1>& fibreLength) const;
    double GetInterNodeSpace() const; 
    std::string GetMeshName() const;
    cp::media_type GetConductivityMedia() const;
    template<unsigned DIM>
    void GetStimuli(std::vector<boost::shared_ptr<AbstractStimulusFunction> >& rStimuliApplied,
                    std::vector<boost::shared_ptr<AbstractChasteRegion<DIM> > >& rStimulatedAreas) const;

    template<unsigned DIM>
    void GetCellHeterogeneities(std::vector<boost::shared_ptr<AbstractChasteRegion<DIM> > >& rCellHeterogeneityRegions,
                                std::vector<double>& rScaleFactorGks,
                                std::vector<double>& rScaleFactorIto,
                                std::vector<double>& rScaleFactorGkr,
                                std::vector<std::map<std::string, double> >* pParameterSettings);

    bool GetConductivityHeterogeneitiesProvided() const; 
    bool AreCellularTransmuralHeterogeneitiesRequested();

    double GetEpiLayerFraction();

    double GetEndoLayerFraction();

    double GetMidLayerFraction();

    unsigned GetEpiLayerIndex();

    unsigned GetEndoLayerIndex();

    unsigned GetMidLayerIndex();


    template<unsigned DIM>
    void GetConductivityHeterogeneities(std::vector<boost::shared_ptr<AbstractChasteRegion<DIM> > >& conductivitiesHeterogeneityAreas,
                                        std::vector< c_vector<double,3> >& intraConductivities,
                                        std::vector< c_vector<double,3> >& extraConductivities) const;
    std::string GetOutputDirectory() const; 
    std::string GetOutputFilenamePrefix() const;

    bool GetOutputVariablesProvided() const;

    void GetOutputVariables(std::vector<std::string>& rOutputVariables) const;

    bool GetOutputUsingOriginalNodeOrdering();

    bool GetCheckpointSimulation() const;

    double GetCheckpointTimestep() const;

    unsigned GetMaxCheckpointsOnDisk() const;

    // ResumeSimulation
    HeartFileFinder GetArchivedSimulationDir() const;


    // Physiological
    void GetIntracellularConductivities(c_vector<double, 3>& rIntraConductivities) const;

    void GetIntracellularConductivities(c_vector<double, 2>& rIntraConductivities) const;

    void GetIntracellularConductivities(c_vector<double, 1>& rIntraConductivities) const;

    void GetExtracellularConductivities(c_vector<double, 3>& rExtraConductivities) const;

    void GetExtracellularConductivities(c_vector<double, 2>& rExtraConductivities) const;

    void GetExtracellularConductivities(c_vector<double, 1>& rExtraConductivities) const;

    double GetBathConductivity(unsigned bathRegion=UINT_MAX) const;

    const std::set<unsigned>& rGetTissueIdentifiers();

    const std::set<unsigned>& rGetBathIdentifiers();


    double GetSurfaceAreaToVolumeRatio() const; 
    double GetCapacitance() const; 
    // Numerical
    double GetOdeTimeStep() const; 
    double GetPdeTimeStep() const; 
    double GetPrintingTimeStep() const; 
    bool GetUseAbsoluteTolerance() const; 
    double GetAbsoluteTolerance() const; 
    bool GetUseRelativeTolerance() const; 
    double GetRelativeTolerance() const;  
    const char* GetKSPSolver() const; 
    const char* GetKSPPreconditioner() const; 
    DistributedTetrahedralMeshPartitionType::type GetMeshPartitioning() const; 
    // Adaptivity
    bool IsAdaptivityParametersPresent() const;

    // Post processing
    bool IsPostProcessingSectionPresent() const;

    void EnsurePostProcessingSectionPresent();

    bool IsPostProcessingRequested() const;

    bool IsApdMapsRequested() const;

    void GetApdMaps(std::vector<std::pair<double,double> >& rApdMaps) const;

    bool IsUpstrokeTimeMapsRequested() const;
    void GetUpstrokeTimeMaps (std::vector<double>& rUpstrokeTimeMaps) const;

    bool IsMaxUpstrokeVelocityMapRequested() const;

    void GetMaxUpstrokeVelocityMaps(std::vector<double>& rUpstrokeVelocityMaps) const;

    bool IsConductionVelocityMapsRequested() const;

    void GetConductionVelocityMaps(std::vector<unsigned>& rConductionVelocityMaps) const;

    bool IsAnyNodalTimeTraceRequested() const;

    void GetNodalTimeTraceRequested(std::vector<unsigned>& rRequestedNodes) const;

    bool IsPseudoEcgCalculationRequested() const;

    template<unsigned SPACE_DIM>
    void GetPseudoEcgElectrodePositions(std::vector<ChastePoint<SPACE_DIM> >& rPseudoEcgElectrodePositions) const;

    bool GetUseStateVariableInterpolation() const;


    // Output visualization

    bool IsOutputVisualizerPresent() const;

    bool GetVisualizeWithMeshalyzer() const;

    bool GetVisualizeWithCmgui() const;

    bool GetVisualizeWithVtk() const;

    bool GetVisualizeWithParallelVtk() const;

    unsigned GetVisualizerOutputPrecision();

    bool IsElectrodesPresent() const;


    void GetElectrodeParameters(bool& rGroundSecondElectrode,
                                unsigned& rIndex, double& rMagnitude,
                                double& rStartTime, double& rDuration );

    bool GetUseMassLumping();

    bool GetUseMassLumpingForPrecond();

    bool GetUseReactionDiffusionOperatorSplitting();

    bool GetUseFixedNumberIterationsLinearSolver();

    unsigned GetEvaluateNumItsEveryNSolves();


    //
    //  Set methods
    //

    // Simulation
    void SetSpaceDimension(unsigned spaceDimension);

    void SetSimulationDuration(double simulationDuration);

    void SetDomain(const cp::domain_type& rDomain);

    void SetDefaultIonicModel(const cp::ionic_models_available_type& rIonicModel);

    void SetSlabDimensions(double x, double y, double z, double inter_node_space);

    void SetSheetDimensions(double x, double y, double inter_node_space);

    void SetFibreLength(double x, double inter_node_space);

    void SetMeshFileName(std::string meshPrefix, cp::media_type fibreDefinition=cp::media_type::NoFibreOrientation);

    void SetConductivityHeterogeneities(std::vector<ChasteCuboid<3> >& rConductivityAreas,
                                        std::vector< c_vector<double,3> >& rIntraConductivities,
                                        std::vector< c_vector<double,3> >& rExtraConductivities);
    void SetConductivityHeterogeneitiesEllipsoid(std::vector<ChasteEllipsoid<3> >& rConductivityAreas,
                                                 std::vector< c_vector<double,3> >& rIntraConductivities,
                                                 std::vector< c_vector<double,3> >& rExtraConductivities);
    void SetOutputDirectory(const std::string& rOutputDirectory);

    void SetOutputFilenamePrefix(const std::string& rOutputFilenamePrefix);

    void SetOutputVariables(const std::vector<std::string>& rOutputVariables);

    void SetOutputUsingOriginalNodeOrdering(bool useOriginal);

     void SetCheckpointSimulation(bool checkpointSimulation, double checkpointTimestep=-1.0, unsigned maxCheckpointsOnDisk=UINT_MAX);

    // Physiological
    void SetIntracellularConductivities(const c_vector<double, 3>& rIntraConductivities);

    void SetIntracellularConductivities(const c_vector<double, 2>& rIntraConductivities);

    void SetIntracellularConductivities(const c_vector<double, 1>& rIntraConductivities);

    void SetExtracellularConductivities(const c_vector<double, 3>& rExtraConductivities);

    void SetExtracellularConductivities(const c_vector<double, 2>& rExtraConductivities);

    void SetExtracellularConductivities(const c_vector<double, 1>& rExtraConductivities);

    void SetBathConductivity(double bathConductivity);

    void SetBathMultipleConductivities(std::map<unsigned, double> bathConductivities);

    void SetTissueAndBathIdentifiers(const std::set<unsigned>& rTissueIds, const std::set<unsigned>& rBathIds);

    //void SetTissueIdentifiers(const std::set<unsigned>& tissueIds);

    void SetSurfaceAreaToVolumeRatio(double ratio);

    void SetCapacitance(double capacitance);

    // Numerical
    void SetOdePdeAndPrintingTimeSteps(double odeTimeStep, double pdeTimeStep, double printingTimeStep);

    void SetOdeTimeStep(double odeTimeStep);

    void SetPdeTimeStep(double pdeTimeStep);

     void SetPrintingTimeStep(double printingTimeStep);

    void SetUseRelativeTolerance(double relativeTolerance);

    void SetUseAbsoluteTolerance(double absoluteTolerance);

    void SetKSPSolver(const char* kspSolver, bool warnOfChange=false);

    void SetKSPPreconditioner(const char* kspPreconditioner);

    void SetMeshPartitioning(const char* meshPartioningMethod);

    void SetApdMaps(const std::vector<std::pair<double,double> >& rApdMaps);

    void SetUpstrokeTimeMaps(std::vector<double>& rUpstrokeTimeMaps);

    void SetMaxUpstrokeVelocityMaps(std::vector<double>& rMaxUpstrokeVelocityMaps);

    void SetConductionVelocityMaps(std::vector<unsigned>& rConductionVelocityMaps);

    void SetRequestedNodalTimeTraces(std::vector<unsigned>& requestedNodes);

    template<unsigned SPACE_DIM>
    void SetPseudoEcgElectrodePositions(const std::vector<ChastePoint<SPACE_DIM> >& rPseudoEcgElectrodePositions);


    // Output visualization

    void EnsureOutputVisualizerExists(void);

    void SetVisualizeWithMeshalyzer(bool useMeshalyzer=true);

    void SetVisualizeWithCmgui(bool useCmgui=true);

    void SetVisualizeWithVtk(bool useVtk=true);

    void SetVisualizeWithParallelVtk(bool useParallelVtk=true);

    void SetVisualizerOutputPrecision(unsigned numberOfDigits);

    void SetElectrodeParameters(bool groundSecondElectrode,
                                unsigned index, double magnitude,
                                double startTime, double duration);

    void SetUseStateVariableInterpolation(bool useStateVariableInterpolation = true);

    void SetUseMassLumping(bool useMassLumping = true);

    void SetUseMassLumpingForPrecond(bool useMassLumping = true);

    void SetUseReactionDiffusionOperatorSplitting(bool useOperatorSplitting = true);

    void SetUseFixedNumberIterationsLinearSolver(bool useFixedNumberIterations = true, unsigned evaluateNumItsEveryNSolves=UINT_MAX);

    bool HasDrugDose() const;

    double GetDrugDose() const;

    void SetDrugDose(double drugDose);

    void SetIc50Value(const std::string& rCurrentName, double ic50, double hill=1.0);

    std::map<std::string, std::pair<double, double> > GetIc50Values();

    //
    // Purkinje-related methods
    //

    bool HasPurkinje();

    double GetPurkinjeCapacitance();

    void SetPurkinjeCapacitance(double capacitance);

    double GetPurkinjeSurfaceAreaToVolumeRatio();

    void SetPurkinjeSurfaceAreaToVolumeRatio(double ratio);

    double GetPurkinjeConductivity();

    void SetPurkinjeConductivity(double conductivity);

private:
    // Only to be accessed by the tests
    friend class TestHeartConfig;

    /*Constructor is private, since the class is only accessed by the singleton instance() method*/
    HeartConfig();

    boost::shared_ptr<cp::chaste_parameters_type> mpParameters;

    static boost::shared_ptr<HeartConfig> mpInstance;

    FileFinder mParametersFilePath;

    bool mUseFixedSchemaLocation;

    double mEpiFraction;

    double mEndoFraction;

    double mMidFraction;

    unsigned mIndexMid;

    unsigned mIndexEpi;

    unsigned mIndexEndo;

    bool mUserAskedForCellularTransmuralHeterogeneities;

    bool mUseMassLumping;

    bool mUseMassLumpingForPrecond;

    bool mUseReactionDiffusionOperatorSplitting;

    std::map<unsigned, double> mBathConductivities;

    std::set<unsigned> mTissueIdentifiers;

    std::set<unsigned> mBathIdentifiers;

    bool mUseFixedNumberIterations;

    unsigned mEvaluateNumItsEveryNSolves;

    void CheckSimulationIsDefined(std::string callingMethod="") const;

    void CheckResumeSimulationIsDefined(std::string callingMethod="") const;
};


BOOST_CLASS_VERSION(HeartConfig, 2)
#include "SerializationExportWrapper.hpp"
// Declare identifier for the serializer
CHASTE_CLASS_EXPORT(HeartConfig)

#endif /*HEARTCONFIG_HPP_*/