NodeMover.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 "NodeMover.h"

#include "bio/Cell.h"
#include "bio/CellAttributes.h"
#include "bio/Mesh.h"
#include "bio/Node.h"
#include "math/DurstenfeldShuffle.h"
#include "math/MeshGeometry.h"
#include "math/RandomEngine.h"
#include "model/ComponentFactoryProxy.h"
#include "sim/CoreData.h"
#include "util/misc/Exception.h"
#include "util/misc/log_debug.h"

#include <trng/uniform_dist.hpp>
#include <trng/uniform01_dist.hpp>
#include <trng/normal_dist.hpp>
#include <algorithm>
#include <array>
#include <set>
#include <unordered_set>

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

namespace SimPT_Sim{

NodeMover::NodeMover()
        : m_mesh(nullptr), m_tissue_energy(0.0)
{
}

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

void NodeMover::CalculateCellEnergies()
{
        m_tissue_energy = 0.0;
        for (auto c : m_mesh->GetCells()) {
                const double e = m_hamiltonian(c);
                m_cell_energies[c] = e;
                m_tissue_energy += e;
        }
        const double e_bp = m_hamiltonian(m_mesh->GetBoundaryPolygon());
        m_cell_energies[m_mesh->GetBoundaryPolygon()] = e_bp;
        m_tissue_energy += e_bp;
}

void NodeMover::CalculateCellEnergiesDelta()
{
        m_tissue_energy = 0.0;
        for (auto c : m_mesh->GetCells()) {
                m_tissue_energy += m_hamiltonian(c);
        }
        m_tissue_energy += m_hamiltonian(m_mesh->GetBoundaryPolygon());
}

void NodeMover::Initialize(const CoreData& cd)
{
        assert((cd.Check() == true) && "CoreData not OK." );

        m_cd                = cd;
        m_cell_energies.clear();

        const auto model_group       = m_cd.m_parameters->get<string>("model.group", "Default");
        const auto factory           = ComponentFactoryProxy::Create(model_group);
        const bool allow_use_delta   = m_cd.m_parameters->get<bool>("cell_mechanics.mc_allow_use_delta", true);

        if (allow_use_delta) {
                m_delta_hamiltonian = factory->CreateDeltaHamiltonian(m_cd);
                assert((m_delta_hamiltonian) && ("DeltaHamiltonian produces nullptr."));
        }
        // Even if allow_use_delta == true, we need the hamiltonian to calculate
        // a baseline energy to determine a convergence rate.
        m_hamiltonian       = factory->CreateHamiltonian(m_cd);
        assert((m_hamiltonian) && ("Hamiltonian produces nullptr."));

        m_mesh              = m_cd.m_mesh.get();
        m_tissue_energy     = 0.0;

        // Initialize the various RNGs.
        for (size_t i = 0; i < 1; ++i) {
                m_mc_move_generator.push_back(factory->CreateMoveGenerator(m_cd));
                m_normal_generator.push_back(m_cd.m_random_engine->GetGenerator(trng::normal_dist<double>(0.0, 1.0), i));
                m_uniform_generator.push_back(m_cd.m_random_engine->GetGenerator(trng::uniform01_dist<double>(), i));
        }

        if (UsesDeltaHamiltonian()) {
                CalculateCellEnergiesDelta();
        } else {
                CalculateCellEnergies();
        }
}

bool NodeMover::IsNonIntersecting(Node* node, const std::array<double, 3>& move) const
{
        const auto& owners = m_mesh->GetNodeOwningNeighbors(node);
        bool non_intersects = true;
        for (const auto& owner : owners) {
                if ((*owner.GetCell()).MoveSelfIntersects(node, *(node)+move)) {
                        non_intersects = false;
                        break;
                }
        }
        // TODO: possibly check all owning cells that have 3 nodes for flipping
        return non_intersects;
}

NodeMover::MetropolisInfo NodeMover::MoveNode(Node* n, double step_size, double temperature)
{
        // Find the node owners (cells) and node owning walls
        const auto& owners          = m_mesh->GetNodeOwningNeighbors(n);
        const auto& owningWalls     = m_mesh->GetNodeOwningWalls(n);

        // Sum the old (current) energy of node owners (cells)
        double old_energy = 0.0;
        double stiff  = 0.0;
        for (const auto& owner : owners) {
                Cell* const cell = owner.GetCell();
                old_energy += m_cell_energies[cell];
                stiff  += cell->GetStiffness();
        }

        // Temporarily move the node.
        const array<double,3> move = m_mc_move_generator[0]() * step_size;
        *n += move;
        for (const auto& owner : owners) {
                owner.GetCell()->SetGeoDirty();
        }
        for (auto& wall : owningWalls) {
                wall->SetLengthDirty();
        }

        // Recompute energies. We do not check validity of the move (i.e. is it
        // non-intersecting) at this point. So, the numbers below might be garbage
        // (in case of intersection the energy caculations produce garbage numbers).
        // But even if, on the basis of the garbage numbers the energy or mc criteria
        // accept the move, it will be rejected because of the intersection.
        double new_energy = 0.0;
        map<Cell*, double> new_energies_cache;
        for (auto& owner : owners) {
                Cell* cell = owner.GetCell();
                const double e = m_hamiltonian(cell);
                new_energies_cache[cell] = e;
                new_energy += e;
        }
        const double dh = new_energy - old_energy;

        // Accept or Reject node movement, provided it is non-intersecting.
        MetropolisInfo info;
        const bool dh_accept = (dh < -stiff);
        const bool mc_accept = m_uniform_generator[0]() < exp((-dh - stiff)/temperature);
        if ((dh_accept || mc_accept) && IsNonIntersecting(n, move)) {
                // ACCEPT: adjust energy vector.
                for (auto& owner : owners) {
                        Cell* cell = owner.GetCell();
                        m_cell_energies[cell] = new_energies_cache[cell];
                }
                m_tissue_energy += dh;
                info = dh_accept ? MetropolisInfo(1U, 0U, dh, 0.0) : MetropolisInfo(0U, 1U, 0.0, dh);
        } else {
                // REJECT: move the node back to original position.
                *n -= move;
                for (auto& owner : owners) {
                        owner.GetCell()->SetGeoDirty();
                }
                for (auto& wall : owningWalls) {
                        wall->SetLengthDirty();
                }
        }

        return info;
}

NodeMover::MetropolisInfo NodeMover::MoveNodeDelta(Node* n, double step_size, double temperature)
{
        // Find the node owners (cells) and node owning walls.
        const auto& owners       = m_mesh->GetNodeOwningNeighbors(n);
        const auto& owningWalls  = m_mesh->GetNodeOwningWalls(n);

        // Sum energy differences of node owners (cells) due to (potential) node move.
        const array<double,3> move = m_mc_move_generator[0]() * step_size;
        double dh    = 0.0;
        double stiff = 0.0;
        for (const auto& owner : owners) {
                dh    += m_delta_hamiltonian(owner, n, move);
                stiff += owner.GetCell()->GetStiffness();
        }

        // Acceptance criteria for Metropolis.
        const bool dh_accept = (dh < -stiff);
        const bool mc_accept = m_uniform_generator[0]() < exp((-dh - stiff)/temperature);

        // Accept (and the execute it) or Reject (and then do nothing).
        MetropolisInfo info;
        if ((dh_accept || mc_accept) && IsNonIntersecting(n, move)) {
                *n += move;
                for (auto& owner : owners) {
                        owner.GetCell()->SetGeoDirty();
                }
                for (auto& wall : owningWalls) {
                        wall->SetLengthDirty();
                }
                m_tissue_energy += dh;
                info = dh_accept ? MetropolisInfo(1U, 0U, dh, 0.0) : MetropolisInfo(0U, 1U, 0.0, dh);
        }

        return info;
}

NodeMover::MetropolisInfo NodeMover::SweepOverNodes(double step_size, double temperature)
{
        // -------------------------------------------------------------------------
        // Select the move method.
        // -------------------------------------------------------------------------
        const auto move_method = UsesDeltaHamiltonian() ? &NodeMover::MoveNodeDelta : &NodeMover::MoveNode;

        // -------------------------------------------------------------------------
        // Shuffle the indices for the sweep.
        // -------------------------------------------------------------------------
        const auto& nodes = m_mesh->GetNodes();
        vector<unsigned int> shuffled_indices(nodes.size());
        DurstenfeldShuffle::Fill(shuffled_indices, *m_cd.m_random_engine);

        // -------------------------------------------------------------------------
        // Loop of all indices, move, accumulate info.
        // -------------------------------------------------------------------------
        MetropolisInfo sweep_info;
        for (const auto& i : shuffled_indices) {
                const auto node = nodes[i];
                if (!node->IsFixed()) {
                        sweep_info += (this->*move_method)(node, step_size, temperature);
                }
        }

        return sweep_info;
}

} // namespace