HeartConfig.hpp

/*

Copyright (C) University of Oxford, 2005-2010

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 "UblasIncludes.hpp"

#include "ArchiveLocationInfo.hpp"
//#include "ChasteParameters_1_1.hpp"
#include "ChasteParameters_2_0.hpp"
#include "SimpleStimulus.hpp"
#include "ChasteCuboid.hpp"
#include "AbstractTetrahedralMesh.hpp"

#include <xercesc/dom/DOM.hpp>
#include <xercesc/dom/DOMDocument.hpp>
#include <xercesc/dom/DOMElement.hpp>


#include <boost/shared_ptr.hpp>

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

namespace cp = chaste::parameters::v2_0;


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 coonfiguration 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, update the simulation duration (and any other parameters to come...)
        if (p_new_parameters->ResumeSimulation().present())
        {
            if ( (p_new_parameters->ResumeSimulation().get().SpaceDimension() != HeartConfig::Instance()->GetSpaceDimension())
                 ||(p_new_parameters->ResumeSimulation().get().Domain() != HeartConfig::Instance()->GetDomain()))
            {
                EXCEPTION("Problem type and space dimension should match when restarting a simulation.");
            }

            HeartConfig::Instance()->SetSimulationDuration(p_new_parameters->ResumeSimulation().get().SimulationDuration());

            if (p_new_parameters->ResumeSimulation().get().CheckpointSimulation().present())
            {
                HeartConfig::Instance()->SetCheckpointSimulation(true,
                                                                 p_new_parameters->ResumeSimulation().get().CheckpointSimulation().get().timestep(),
                                                                 p_new_parameters->ResumeSimulation().get().CheckpointSimulation().get().max_checkpoints_on_disk());
            }

            //Visualization parameters are compulsory
            HeartConfig::Instance()->SetVisualizeWithVtk(p_new_parameters->ResumeSimulation().get().OutputVisualizer().vtk() == cp::yesno_type::yes);
            HeartConfig::Instance()->SetVisualizeWithCmgui(p_new_parameters->ResumeSimulation().get().OutputVisualizer().cmgui() == cp::yesno_type::yes);
            HeartConfig::Instance()->SetVisualizeWithMeshalyzer(p_new_parameters->ResumeSimulation().get().OutputVisualizer().meshalyzer() == cp::yesno_type::yes);
        }
    }
    BOOST_SERIALIZATION_SPLIT_MEMBER()

    
    xsd::cxx::xml::dom::auto_ptr<xercesc::DOMDocument> ReadFileToDomDocument(
        const std::string& rFileName,
        ::xml_schema::error_handler& rErrorHandler,
        const ::xml_schema::properties& rProps);

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

private:
    SchemaLocationsMap mSchemaLocations;

    void SetDefaultSchemaLocations();

    std::string EscapeSpaces(const std::string& rPath);

    xercesc::DOMElement* AddNamespace(xercesc::DOMDocument* pDocument,
                                      xercesc::DOMElement* pElement,
                                      const std::string& rNamespace);

    xercesc::DOMElement* AddNamespace(xercesc::DOMDocument* pDocument,
                                      xercesc::DOMElement* pElement,
                                      const XMLCh* pNamespace);


    void TransformIonicModelDefinitions(xercesc::DOMDocument* pDocument,
                                        xercesc::DOMElement* pRootElement);

    void WrapContentInElement(xercesc::DOMDocument* pDocument,
                              xercesc::DOMElement* pElement,
                              const XMLCh* pNewElementLocalName);

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

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

    static void Reset();

    /*
     *  Get methods
     */

    // 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<ChasteCuboid<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<SimpleStimulus> >& rStimuliApplied, std::vector<ChasteCuboid<DIM> >& rStimulatedAreas) const;

    template<unsigned DIM>
    void GetCellHeterogeneities( std::vector<AbstractChasteRegion<DIM>* >& rCellHeterogeneityRegions,
                                 std::vector<double>& rScaleFactorGks,
                                 std::vector<double>& rScaleFactorIto,
                                 std::vector<double>& rScaleFactorGkr);
    bool GetConductivityHeterogeneitiesProvided() const; 
    bool AreCellularTransmuralHeterogeneitiesRequested();

    bool AreCellularlHeterogeneitiesSpecifiedByCuboids();

    double GetEpiLayerFraction();

    double GetEndoLayerFraction();

    double GetMidLayerFraction();

    unsigned GetEpiLayerIndex();

    unsigned GetEndoLayerIndex();

    unsigned GetMidLayerIndex();


    template<unsigned DIM>
    void GetConductivityHeterogeneities(std::vector<ChasteCuboid<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> &outputVariables) const;

    bool GetCheckpointSimulation() const;

    double GetCheckpointTimestep() const;

    unsigned GetMaxCheckpointsOnDisk() const;

    // ResumeSimulation
    std::string 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() const; 
    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; 
    // 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;

    bool IsPostProcessingRequested() const;

    bool IsApdMapsRequested() const;
    void GetApdMaps(std::vector<std::pair<double,double> >& apdMaps) const;

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

    bool IsMaxUpstrokeVelocityMapRequested() const;

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

    bool IsConductionVelocityMapsRequested() const;

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


    // Output visualization

    bool IsOutputVisualizerPresent() const;

    bool GetVisualizeWithMeshalyzer() const;

    bool GetVisualizeWithCmgui() const;

    bool GetVisualizeWithVtk() const;

     /*
     *  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> >& 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 SetCheckpointSimulation(bool checkpointSimulation, double checkpointTimestep=-1.0, unsigned maxCheckpointsOnDisk=UINT_MAX);


    // Physiological
    void SetIntracellularConductivities(const c_vector<double, 3>& intraConductivities);
    void SetIntracellularConductivities(const c_vector<double, 2>& intraConductivities);
    void SetIntracellularConductivities(const c_vector<double, 1>& intraConductivities);

    void SetExtracellularConductivities(const c_vector<double, 3>& extraConductivities);
    void SetExtracellularConductivities(const c_vector<double, 2>& extraConductivities);
    void SetExtracellularConductivities(const c_vector<double, 1>& extraConductivities);

    void SetBathConductivity(double bathConductivity);

    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 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> >& apdMaps);

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

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

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


    // Output visualization

    void EnsureOutputVisualizerExists(void);

    void SetVisualizeWithMeshalyzer(bool useMeshalyzer=true);

    void SetVisualizeWithCmgui(bool useCmgui=true);

    void SetVisualizeWithVtk(bool useVtk=true);


    ~HeartConfig(); 
protected:
    // Only to be accessed by the tests
    friend class TestHeartConfig;

private:
    /*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;

    bool mUseFixedSchemaLocation;

    double mEpiFraction;

    double mEndoFraction;

    double mMidFraction;

    unsigned mIndexMid;

    unsigned mIndexEpi;

    unsigned mIndexEndo;

    bool mUserAskedForCellularTransmuralHeterogeneities;

    bool mUserAskedForCuboidsForCellularHeterogeneities;

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

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

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

};

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

#endif /*HEARTCONFIG_HPP_*/