Grow.cpp

Go to the documentation of this file.

/*
 * Copyright 2011-2016 Universiteit Antwerpen
 *
 * Licensed under the EUPL, Version 1.1 or  as soon they will be approved by
 * the European Commission - subsequent versions of the EUPL (the "Licence");
 * You may not use this work except in compliance with the Licence.
 * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the Licence is distributed on an "AS IS" basis,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the Licence for the specific language governing
 * permissions and limitations under the Licence.
 */
#include "Grow.h"

#include "algo/NodeInserter.h"
#include "algo/NodeMover.h"
#include "bio/Mesh.h"
//#include "math/CSRMatrix.h"
//#include "math/MeshTopology.h"
//#include "util/config/ConfigInfo.h"

#include <boost/circular_buffer.hpp>
//#include <omp.h>
#include <algorithm>
#include <cmath>
#include <tuple>

namespace SimPT_Default {
namespace TimeEvolver {

using namespace std;
using namespace boost::property_tree;
using namespace SimPT_Sim::Util;

Grow::Grow(const CoreData& cd)
{
        Initialize(cd);
}

void Grow::Initialize(const CoreData& cd)
{
        assert( cd.Check() && "CoreData not ok.");
        m_cd                      = cd;

        const auto& p             = m_cd.m_parameters;
        m_mc_cell_step_size       = p->get<double>("cell_mechanics.mc_cell_stepsize");

        m_mc_retry_limit          = p->get<unsigned int>("cell_mechanics.mc_retry_limit");
        m_mc_sliding_window       = p->get<unsigned int>("cell_mechanics.mc_sliding_window");
        m_mc_step_size            = p->get<double>("cell_mechanics.mc_stepsize");
        m_mc_sweep_limit          = p->get<double>("cell_mechanics.mc_sweep_limit");
        m_mc_temperature          = p->get<double>("cell_mechanics.mc_temperature");
        m_stationarity_check      = p->get<bool>("termination.stationarity_check");

        const double abs_tol      = abs(p->get<double>("cell_mechanics.mc_abs_tolerance"));
        m_mc_abs_tolerance        = abs_tol;
        const double rel_tol      = abs(p->get<double>("cell_mechanics.mc_rel_tolerance"));
        const bool stringent      = p->get<bool>("cell_mechanics.mc_stringent_tolerance", false);
        if (stringent) {
                m_within_tolerance = [abs_tol, rel_tol](double dh, double e) {
                                        return (dh <= 0.0) && (abs(dh) < min(abs_tol, abs(e * rel_tol)));
                                };
        } else {
                m_within_tolerance = [abs_tol, rel_tol](double dh, double e) {
                                        return (dh <= 0.0) && (abs(dh) < max(abs_tol, abs(e * rel_tol)));
                                };
        }
}

std::tuple<SimTimingTraits::CumulativeTimings, bool> Grow::operator()(double, SimPhase)
{
        // --------------------------------------------------------------------------
        // Intro
        // --------------------------------------------------------------------------
        Stopclock chrono_grow("growth", true);

        // --------------------------------------------------------------------------
        // Node insertion.
        // --------------------------------------------------------------------------
        NodeInserter node_inserter;
        node_inserter.Initialize(m_cd);
        const unsigned int nodes_inserted = node_inserter.InsertNodes();

        // --------------------------------------------------------------------------
        // Node motion (elastic relaxation). No time increase here.
        // --------------------------------------------------------------------------
        NodeMover node_mover;
        node_mover.Initialize(m_cd);
        double temperature = m_mc_temperature;
        const double alpha_temperature = 0.95;  // annealing schedule: T -> alpha * T
        const double e_tissue_before = node_mover.GetTissueEnergy();

        // ----------- Energy sliding window and counters.
        boost::circular_buffer<double> cb(m_mc_sliding_window);
        double mean_dh = 0.0;
        unsigned int sweep_count = 0U;
        unsigned int retry_count = 0U;
        unsigned int total_move_down_count = 0U;
        bool go_on = false;

        // ----------- Calculate cell based node incidence and edge based node incidence if needed.
//      CSRMatrix cell_incidence;
//      CSRMatrix edge_incidence;
//      if (ConfigInfo::HaveOpenMP()) {
//              cell_incidence = MeshTopology::NodeCellNodeIncidence(m_cd.m_mesh);
//              if (node_mover.UsesDeltaHamiltonian()) {
//                      edge_incidence = MeshTopology::NodeEdgeNodeIncidence(m_cd.m_mesh);
//              }
//      }

        // ----------- Sweep over mesh until convergence criteria are satisfied.
        do {
                ++sweep_count;
                ++retry_count;

                // Actual node moves.
                NodeMover::MetropolisInfo info;
//              if (ConfigInfo::HaveOpenMP() && node_mover.UsesDeltaHamiltonian()) {
//                      info = node_mover.SweepOverNodes(cell_incidence, edge_incidence, m_mc_step_size, temperature);
//              } else if (ConfigInfo::HaveOpenMP()) {
//                        info = node_mover.SweepOverNodes(cell_incidence, m_mc_step_size, temperature);
//
//                } else {
//                        info = node_mover.SweepOverNodes(m_mc_step_size, temperature);
//                }

                info = node_mover.SweepOverNodes(m_mc_step_size, temperature);

                // Process the situation.
                total_move_down_count += info.down_count;
                cb.push_back(info.down_dh + info.up_dh);
                mean_dh = std::accumulate(cb.begin(), cb.end(), 0.0) / cb.size();

                // If the energy has dropped a lot we reset the retry counter for
                // to zero because we need to explore further.
                retry_count = mean_dh < -abs(m_mc_abs_tolerance) ? 0U : retry_count;

                // Stopping condition:
                // Criterion 1: We want to decide depending on average energy drop
                //     over a number of sweeps. So be sure you've made that many sweeps.
                // Criterion 2: If the energy has dropped a lot (defined as more than the
                //     tolerance value, either absolute or relative) we want to go on.
                // Criterion 3: If there was an energy increase, we want to retry a number
                //     of times before giving up.
                go_on = (sweep_count < m_mc_sliding_window) || (! m_within_tolerance(mean_dh, e_tissue_before))
                                        || (mean_dh > 0.0 && retry_count < m_mc_retry_limit);

                // Adjust temperature for next sweep.
                temperature *= alpha_temperature;

        } while (go_on && sweep_count < m_mc_sweep_limit);

        const bool is_stationary = (nodes_inserted == 0U) && (total_move_down_count == 0U);

        // -------------------------------------------------------------------------------
        // We are done ...
        // -------------------------------------------------------------------------------
        CumulativeTimings timings;
        timings.Record(chrono_grow.GetName(), chrono_grow.Get());
        return make_tuple(timings, is_stationary);
}

} // namespace
} // namespace