HeartConfig.hpp

/*

Copyright (C) University of Oxford, 2005-2011

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.

Chaste is free software: you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation, either version 2.1 of the License, or
(at your option) any later version.

Chaste is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details. The offer of Chaste under the terms of the
License is subject to the License being interpreted in accordance with
English Law and subject to any action against the University of Oxford
being under the jurisdiction of the English Courts.

You should have received a copy of the GNU Lesser General Public License
along with Chaste. If not, see <http://www.gnu.org/licenses/>.

*/


#ifndef HEARTCONFIG_HPP_
#define HEARTCONFIG_HPP_

#include <string>
#include <vector>
#include <set>

#include "UblasIncludes.hpp"

#include "ArchiveLocationInfo.hpp"
#include "ChasteParameters_3_0.hpp"

#include "AbstractStimulusFunction.hpp"
// These are needed here for Boost < 1.37
#include "SimpleStimulus.hpp"
#include "RegularStimulus.hpp"

#include "ChastePoint.hpp"
#include "ChasteCuboid.hpp"
#include "ChasteEllipsoid.hpp"
#include "DistributedTetrahedralMeshPartitionType.hpp"

#include <boost/shared_ptr.hpp>

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

namespace cp = chaste::parameters::v3_0;

// 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 );
        }
        PetscTools::Barrier("HeartConfig::save");
    }

    template<class Archive>
    void load(Archive & archive, const unsigned int version)
    {
        /*
         *  This method implements the logic required by HeartConfig to be able to handle resuming a simulation via the executable.
         *
         *  When the control reaches the method mpUserParameters and mpDefaultParameters point to the files specified as resuming parameters.
         *  However SetDefaultsFile() and SetParametersFile() will set those variables to point to the archived parameters.
         *
         *  We make a temporary copy of mpUserParameters so we don't lose its content. At the end of the method we update the new mpUserParameters
         *  with the resuming parameters.
         */
        assert(mpUserParameters.use_count() > 0);
        boost::shared_ptr<cp::chaste_parameters_type> p_new_parameters = mpUserParameters;

        std::string defaults_filename_xml = ArchiveLocationInfo::GetArchiveDirectory() + "ChasteDefaults.xml";
        HeartConfig::Instance()->SetDefaultsFile(defaults_filename_xml);

        /*
         *  When we unarchive a simulation, we load the old parameters file in order to inherit things such
         *  as default cell model, stimuli, heterogeneities, ... This has the side effect of inheriting the
         *  <CheckpointSimulation> element (if defined).
         *
         *  We disable checkpointing definition coming from the unarchived config file. We will enable it again
         *  if defined in the resume config file.
         */
        std::string parameters_filename_xml = ArchiveLocationInfo::GetArchiveDirectory() + "ChasteParameters.xml";
        HeartConfig::Instance()->SetParametersFile(parameters_filename_xml);

        HeartConfig::Instance()->SetCheckpointSimulation(false);

        // If we are resuming a simulation, some parameters can be altered at this point.
        if (p_new_parameters->ResumeSimulation().present())
        {
            UpdateParametersFromResumeSimulation(p_new_parameters);
        }
    }
    BOOST_SERIALIZATION_SPLIT_MEMBER()

    
    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 SetDefaultsFile(const std::string& rFileName);

    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>& intraConductivities) const;
    void GetIntracellularConductivities(c_vector<double, 2>& intraConductivities) const;
    void GetIntracellularConductivities(c_vector<double, 1>& intraConductivities) const;

    void GetExtracellularConductivities(c_vector<double, 3>& extraConductivities) const;
    void GetExtracellularConductivities(c_vector<double, 2>& extraConductivities) const;
    void GetExtracellularConductivities(c_vector<double, 1>& extraConductivities) 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;

    double GetTargetErrorForAdaptivity() const;

    double GetSigmaForAdaptivity() const;

    double GetMaxEdgeLengthForAdaptivity() const;

    double GetMinEdgeLengthForAdaptivity() const;

    double GetGradationForAdaptivity() const;

    unsigned GetMaxNodesForAdaptivity() const;

    unsigned GetNumberOfAdaptiveSweeps() 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> >& conductivityAreas,
            std::vector< c_vector<double,3> >& intraConductivities,
            std::vector< c_vector<double,3> >& extraConductivities);
    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>& tissueIds, const std::set<unsigned>& bathIds);

    //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);

    void SetKSPPreconditioner(const char* kspPreconditioner);

    void SetMeshPartitioning(const char* meshPartioningMethod);

    void SetAdaptivityParameters(double targetError, double sigma, double maxEdgeLength, double minEdgeLength,
                                 double gradation, unsigned maxNodes, unsigned numSweeps );

    void SetTargetErrorForAdaptivity(double targetError);

    void SetSigmaForAdaptivity(double sigma);

    void SetMaxEdgeLengthForAdaptivity(double maxEdgeLength);

    void SetMinEdgeLengthForAdaptivity(double minEdgeLength);

    void SetGradationForAdaptivity(double gradation);

    void SetMaxNodesForAdaptivity(unsigned maxNodes);

    void SetNumberOfAdaptiveSweeps(unsigned numSweeps);

    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();

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> mpUserParameters;

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

    static std::auto_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;

    template<class TYPE>
    TYPE* DecideLocation(TYPE* params_ptr, TYPE* defaults_ptr, const std::string& rNameParameter) const;

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

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


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

#endif /*HEARTCONFIG_HPP_*/