midpath_check.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//include/bertini2/nag_algorithms/midpath_check.hpp is free software: you can redistribute it and/or modify
//it under the terms of the GNU General Public License as published by
//the Free Software Foundation, either version 3 of the License, or
//(at your option) any later version.
//
//include/bertini2/nag_algorithms/midpath_check.hpp 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 General Public License for more details.
//
//You should have received a copy of the GNU General Public License
//along with include/bertini2/nag_algorithms/midpath_check.hpp.  If not, see <http://www.gnu.org/licenses/>.
//
// Copyright(C) Bertini2 Development Team
//
// See <http://www.gnu.org/licenses/> for a copy of the license,
// as well as COPYING.  Bertini2 is provided with permitted
// additional terms in the b2/licenses/ directory.

 
 
#pragma once
 
#include "bertini2/nag_algorithms/common/config.hpp"
#include "bertini2/detail/configured.hpp"
 
namespace bertini{
    namespace algorithm{
 

        template<typename RealT, typename ComplexT, typename MetaDataType>
        struct MidpathChecker : public detail::Configured<MidPathConfig>
        {
            using MidPathConfT = MidPathConfig;                 
            using AlgConf = detail::Configured<MidPathConfig>;  
            using BoundaryData = SolnCont< MetaDataType >;      
            using PathIndT = unsigned long long;                
 
 
            MidpathChecker() = default;

            template<typename ...T>
            MidpathChecker( T const& ...config) : AlgConf(config...)
            {
                passed_ = true;
            }
 

            struct CrossedPath{

                CrossedPath(PathIndT index, PathIndT crossed_with, bool same_start)
                {
                    index_ = index;
                    crossed_with_.push_back( std::make_pair(crossed_with, same_start) );
                    rerun_ = !same_start;
                }
 

                PathIndT index() const
                {
                    return index_;
                }

                void crossed_with(std::pair<PathIndT, bool> crossed)
                {
                    crossed_with_.push_back(crossed);
                }

                void rerun(bool same_start)
                {
                    rerun_ = rerun_ || !same_start;
                }

                bool rerun() const
                {
                    return rerun_;
                }
 
 
            private:
                PathIndT index_;
                std::vector< std::pair<PathIndT, bool> > crossed_with_; // first is index of path crossed with, second is bool = true if had same starting point as crossed path
                bool rerun_;
 
 
            }; //re: struct CrossedPath
 
 

            bool Passed() const
            {
                return passed_;
            }
 

            const auto& SamePointTol() const
            {
                return this->template Get<MidPathConfT>().same_point_tolerance;
            }

            template <typename StartSystemT>
            bool Check(BoundaryData const& boundary_data, StartSystemT const& start_system)
            {
                // Each Check is independent: the algorithm calls this repeatedly (once before any
                // resolve, then again after each MidpathResolve).  Reset the accumulated state so a
                // later call reports *this* call's crossings, not the union with prior calls.
                // Without this, passed_ (only ever set false below) would stay false forever once
                // any crossing was seen, and crossed_paths_ would grow without bound -- the resolve
                // loop in EGBoundaryAction would then never see a clean pass and would re-track
                // stale indices.
                passed_ = true;
                crossed_paths_.clear();
 
                for (PathIndT ii = 0; ii < boundary_data.size(); ++ii)
                {
                    if ( boundary_data[ii].success_code != SuccessCode::Success)
                            continue;
 
                    const Vec<ComplexT>& solution_ii = boundary_data[ii].path_point;
                    const auto start_ii = start_system.template StartPoint<ComplexT>(ii);
 
                    for (PathIndT jj = ii+1; jj < boundary_data.size(); ++jj)
                    {
                        if ( boundary_data[jj].success_code != SuccessCode::Success)
                            continue;
 
                        const Vec<ComplexT>& solution_jj = boundary_data[jj].path_point;
                        const Vec<ComplexT> diff_sol = solution_ii - solution_jj;
 
                        if ((diff_sol.template lpNorm<Eigen::Infinity>()/solution_ii.template lpNorm<Eigen::Infinity>()) < SamePointTol())
                        {
                            bool i_already_stored = false;
                            bool j_already_stored = false;
                            // Check if start points are the same
 
                            const auto start_jj = start_system.template StartPoint<ComplexT>(jj);
                            auto diff_start = start_ii - start_jj;
                            // same_start is true when the two paths *began at the same start point*
                            // (their start points differ by less than the tolerance).  That governs
                            // whether this coincidence is a genuine path crossing worth re-tracking:
                            // CrossedPath sets rerun_ = !same_start, so two *distinct* starts that
                            // have collided (same_start == false) get re-tracked, while two paths
                            // that legitimately share a start point do not.  NOTE: this comparison
                            // was historically inverted (>), which flipped rerun and meant genuine
                            // crossings were never re-tracked.
                            bool same_start = (diff_start.template lpNorm<Eigen::Infinity>() < SamePointTol());
 
 
                            // Check if path has already been stored in crossed_paths_
                            for(auto& v : crossed_paths_)
                            {
                                if(v.index() == ii)
                                {
                                    v.crossed_with(std::make_pair(jj,same_start));
                                    i_already_stored = true;
                                    v.rerun(same_start);
                                }
                                if(v.index() == jj)
                                {
                                    v.crossed_with(std::make_pair(ii,same_start));
                                    j_already_stored = true;
                                    v.rerun(same_start);
                                }
                            }
 
                            // If not already stored, create CrossedPath object and add to crossed_paths_
                            if(!i_already_stored)
                            {
                                CrossedPath tempPath(ii, jj, same_start);
                                crossed_paths_.push_back(tempPath);
                            }
                            if(!j_already_stored)
                            {
                                CrossedPath tempPath(jj, ii, same_start);
                                crossed_paths_.push_back(tempPath);
                            }
 
                            passed_ = false;
 
                        }
                    }
                }
                return passed_;
            };
 

            const std::vector<CrossedPath> GetCrossedPaths() const
            {
                return crossed_paths_;
            }
 
 
 
        private:
 
            std::vector<CrossedPath> crossed_paths_; // Data for all paths that crossed on last check
            bool passed_;  // Did the check pass?
 
 
        }; //re: struct MidpathChecker
 
    } // re: namespace algorithm
}// re: namespace bertini