/**
 * \file
 * \date 2014-09-19
 * \brief Implementation of heuristic search length strategy.
 *
 * \copyright
 * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 */

#include "HeuristicSearchLength.h"

#include "BaseLib/Logging.h"
#include "MeshLib/Elements/Element.h"

namespace MeshGeoToolsLib
{
HeuristicSearchLength::HeuristicSearchLength(MeshLib::Mesh const& mesh,
                                             LengthType length_type)
    : _mesh(mesh)
{
    double sum(0.0);
    double sum_of_sqr(0.0);  // total length of edges

    std::size_t n_sampling(0);  // total length of edges squared
    std::vector<MeshLib::Element*> const& elements(_mesh.getElements());

    if (length_type == LengthType::Edge)
    {
        for (auto element : elements)
        {
            std::size_t const n_edges(element->getNumberOfEdges());
            for (std::size_t k(0); k < n_edges; k++)
            {
                auto edge =
                    element->getEdge(k);  // allocation inside getEdge().
                double const len = edge->getContent();
                delete edge;
                sum += len;
                sum_of_sqr += len * len;
            }
            n_sampling += n_edges;
        }
    }
    else
    {
        for (const MeshLib::Element* e : elements)
        {
            auto const [min, max] =
                MeshLib::computeSqrNodeDistanceRange(*e, true);
            sum += std::sqrt(min);
            sum_of_sqr += min;
        }
        n_sampling = _mesh.getNumberOfElements();
    }

    const double mean(sum / n_sampling);
    const double variance((sum_of_sqr - (sum * sum) / n_sampling) /
                          (n_sampling - 1));

    // Set the search length for the case of non-positive variance (which can
    // happen due to numerics).
    _search_length = mean / 2;

    if (variance > 0)
    {
        if (variance < mean * mean / 4)
        {
            _search_length -= std::sqrt(variance);
        }
        else
        {
            _search_length = std::numeric_limits<double>::epsilon();
        }
    }

    DBUG(
        "[MeshNodeSearcher::MeshNodeSearcher] Calculated search length for "
        "mesh '{:s}' is {:f}.",
        _mesh.getName(), _search_length);
}

}  // end namespace MeshGeoToolsLib